Hodgkin-Huxley Neuron — Action Potential from Ionic Currents

1. The Neuron as an Electrical Circuit

The plasma membrane separates charge across its ~7 nm thickness, acting as a capacitor with capacitance C_m ≈ 1 µF/cm². Ion channels — protein pores that open and close in response to voltage — act as voltage-dependent resistors. The driving force on each ion species is the difference between the membrane potential V and the Nernst (equilibrium) potential E.

In Hodgkin and Huxley's voltage-clamp experiments on squid giant axons, they measured the conductances g_Na and g_K directly and extracted the numerical parameters that fit the data. The parameter set below reproduces the original 1952 results.

2. The Four Hodgkin-Huxley Equations

The full system is four coupled ODEs: one for voltage V and three for the gating variables m, h, n — each representing the probability that a particular gate in the ion channel protein is in the open state.

Each gate variable obeys first-order kinetics: it relaxes toward its steady-state value m_∞(V) = α_m/(α_m + β_m) with a time constant τ_m(V) = 1/(α_m + β_m). The rates α and β are empirical functions of V fit to voltage-clamp data.

The m³ term in I_Na reflects the fact that the sodium channel requires three independent activation gates to be open simultaneously. The h factor is an inactivation gate: it starts open (~0.6) but closes as the membrane depolarises, cutting off the Na⁺ current and terminating the action potential. The n⁴ term in I_K reflects four identical subunits.

3. Gate Variables m, h, n — Alpha/Beta Kinetics

The rate functions below are valid with voltage V measured in mV. They contain removable singularities at specific voltages (when the denominator → 0), which must be handled numerically with L'Hôpital's rule.

// All voltages in mV; rates in ms⁻¹
// V is the deviation from rest: use V_shifted = V_membrane − (−65)

function alpha_m(V) {
  const dv = 10 - V;
  return Math.abs(dv) < 1e-7
    ? 1                                    // L'Hôpital limit
    : 0.1 * dv / (Math.exp(dv / 10) - 1);
}
function beta_m(V)  { return 4    * Math.exp(-V / 18); }

function alpha_h(V) { return 0.07 * Math.exp(-V / 20); }
function beta_h(V)  { return 1    / (Math.exp((30 - V) / 10) + 1); }

function alpha_n(V) {
  const dv = 10 - V;
  return Math.abs(dv) < 1e-7
    ? 0.1
    : 0.01 * dv / (Math.exp(dv / 10) - 1);
}
function beta_n(V)  { return 0.125 * Math.exp(-V / 80); }

4. Numerical Integration with RK4

The action potential lasts ~1 ms and the fastest gate variable (m) has τ_m ≈ 0.1 ms at threshold. A safe time step is dt = 0.01 ms. The Euler method is unstable for this system at that step size; RK4 is both stable and fourth-order accurate.

function derivatives(V, m, h, n, Iext) {
  const INa = 120 * m**3 * h * (V - 50);
  const IK  =  36 * n**4     * (V + 77);
  const IL  = 0.3             * (V + 54.4);
  return {
    dV: Iext - INa - IK - IL,             // Cm = 1
    dm: alpha_m(V) * (1 - m) - beta_m(V) * m,
    dh: alpha_h(V) * (1 - h) - beta_h(V) * h,
    dn: alpha_n(V) * (1 - n) - beta_n(V) * n,
  };
}

function stepRK4(V, m, h, n, Iext, dt) {
  const k1 = derivatives(V, m, h, n, Iext);
  const k2 = derivatives(V + dt/2*k1.dV, m + dt/2*k1.dm,
                          h + dt/2*k1.dh, n + dt/2*k1.dn, Iext);
  const k3 = derivatives(V + dt/2*k2.dV, m + dt/2*k2.dm,
                          h + dt/2*k2.dh, n + dt/2*k2.dn, Iext);
  const k4 = derivatives(V + dt*k3.dV, m + dt*k3.dm,
                          h + dt*k3.dh, n + dt*k3.dn, Iext);
  return {
    V: V + dt/6 * (k1.dV + 2*k2.dV + 2*k3.dV + k4.dV),
    m: m + dt/6 * (k1.dm + 2*k2.dm + 2*k3.dm + k4.dm),
    h: h + dt/6 * (k1.dh + 2*k2.dh + 2*k3.dh + k4.dh),
    n: n + dt/6 * (k1.dn + 2*k2.dn + 2*k3.dn + k4.dn),
  };
}

5. Action Potential Phases

Apply a brief current pulse (I_ext = 10 µA/cm² for 0.5 ms) and the membrane potential traces through four distinct phases:

The all-or-nothing property arises from the threshold: small stimuli cause a sub-threshold response and decay; stimuli above ~−55 mV trigger the full spike. This is an emergent property of the non-linear coupling between m, h, and V.

6. Complete JavaScript Implementation

class HHNeuron {
  constructor() {
    this.V  = -10;   // deviation from rest: −65 mV absolute → 0 here
    this.m  = this._minf(this.V);
    this.h  = this._hinf(this.V);
    this.n  = this._ninf(this.V);
    this.dt = 0.01;   // ms
  }

  _alpha_m(V) { const d=10-V; return Math.abs(d)<1e-7 ? 1 : 0.1*d/(Math.exp(d/10)-1); }
  _beta_m(V)  { return 4    * Math.exp(-V/18); }
  _alpha_h(V) { return 0.07 * Math.exp(-V/20); }
  _beta_h(V)  { return 1    / (Math.exp((30-V)/10)+1); }
  _alpha_n(V) { const d=10-V; return Math.abs(d)<1e-7 ? 0.1 : 0.01*d/(Math.exp(d/10)-1); }
  _beta_n(V)  { return 0.125 * Math.exp(-V/80); }

  _minf(V) { const a=this._alpha_m(V); return a/(a+this._beta_m(V)); }
  _hinf(V) { const a=this._alpha_h(V); return a/(a+this._beta_h(V)); }
  _ninf(V) { const a=this._alpha_n(V); return a/(a+this._beta_n(V)); }

  _deriv(V, m, h, n, I) {
    return {
      dV: I - 120*m**3*h*(V-50) - 36*n**4*(V+77) - 0.3*(V+54.4),
      dm: this._alpha_m(V)*(1-m) - this._beta_m(V)*m,
      dh: this._alpha_h(V)*(1-h) - this._beta_h(V)*h,
      dn: this._alpha_n(V)*(1-n) - this._beta_n(V)*n,
    };
  }

  step(Iext = 0) {
    const dt = this.dt;
    let {V, m, h, n} = this;
    const k1 = this._deriv(V, m, h, n, Iext);
    const k2 = this._deriv(V+dt/2*k1.dV, m+dt/2*k1.dm, h+dt/2*k1.dh, n+dt/2*k1.dn, Iext);
    const k3 = this._deriv(V+dt/2*k2.dV, m+dt/2*k2.dm, h+dt/2*k2.dh, n+dt/2*k2.dn, Iext);
    const k4 = this._deriv(V+dt*k3.dV, m+dt*k3.dm, h+dt*k3.dh, n+dt*k3.dn, Iext);
    this.V = V + dt/6*(k1.dV+2*k2.dV+2*k3.dV+k4.dV);
    this.m = m + dt/6*(k1.dm+2*k2.dm+2*k3.dm+k4.dm);
    this.h = h + dt/6*(k1.dh+2*k2.dh+2*k3.dh+k4.dh);
    this.n = n + dt/6*(k1.dn+2*k2.dn+2*k3.dn+k4.dn);
    return this.V + (-65);  // convert to absolute mV
  }
}

// Run a 30 ms simulation with a 10 µA/cm² pulse from t=1 to t=2 ms
const neuron = new HHNeuron();
const trace  = [];
for (let t = 0; t <= 30; t += 0.01) {
  const I = (t >= 1 && t <= 2) ? 10 : 0;
  trace.push({ t, V: neuron.step(I) });
}
// trace now contains [{ t: 0, V: -65 }, ..., { t: 5.2, V: +38 }, ...]

7. 2D Excitable Tissue and Spiral Waves

When many HH neurons are coupled on a grid via gap junctions (diffusive coupling), the system becomes an excitable medium. Voltage spreads from excited cells to resting neighbours: this is the cardiac or cortical wave. The cable equation adds a diffusion term to the voltage equation:

In 2D, a broken planar wave curls into a re-entrant spiral — the computational analogue of dangerous cardiac arrhythmias like ventricular fibrillation. The spiral tip rotates with a period determined by the refractory time of the tissue. Cross-field stimulation during the vulnerable window can induce or terminate spirals.

For a simple implementation, each cell is a separate HH neuron. Coupling is added as an extra current: I_diff = D × Σ(V_neighbour − V_self). A 100×100 grid with D ≈ 0.5 cm²/s and dt = 0.05 ms runs comfortably in real-time on modern hardware.

8. Extensions and Applications

FitzHugh-Nagumo: A Reduced Model

For large-scale tissue simulations, the full 4-variable HH model is expensive. The FitzHugh-Nagumo model reduces it to 2 variables (fast activator u, slow recovery variable w) while preserving excitability and the bifurcation structure. It runs 4× faster and captures spiral waves, travelling pulses, and oscillations.

Multi-Compartment Neurons

Real neurons are not point objects. Dendritic arbours consist of many cylindrical compartments coupled by axial resistance. The Cable Equation (Rall's compartmental model) treats each segment as an HH-like patch connected to neighbours by a conductance g_c = A/(R_a·Δx), giving a linear system of ODEs that can be solved with implicit methods (Crank-Nicolson) for stability.

Stochastic Ion Channels

Individual channels are discrete — they either open or close — and their behaviour in small neurons is probabilistic. Stochastic HH replaces ODEs for m, h, n with binomial Markov chain updates, generating channel noise and spontaneous spiking at subthreshold currents. This is implemented with Binomial(N_channel, α·dt) transitions.