Simulate a Hodgkin-Huxley Neuron in the Browser

1. The Four Hodgkin-Huxley Equations

The membrane behaves like a capacitor in parallel with three variable resistors (ion channels). Current balance gives the master equation for voltage V, plus one first-order kinetic equation per gating variable:

Standard squid-axon parameters (Hodgkin & Huxley, 1952):

• C_m = 1 μF/cm² — membrane capacitance
• g_Na = 120 mS/cm², E_Na = 50 mV — sodium conductance and reversal potential
• g_K = 36 mS/cm², E_K = −77 mV — potassium conductance and reversal potential
• g_L = 0.3 mS/cm², E_L = −54.4 mV — leak conductance and reversal potential

m³h models fast Na⁺ activation (m, three independent gates) combined with slower Na⁺ inactivation (h, one gate); n⁴ models delayed-rectifier K⁺ activation (four gates). These specific exponents were fit empirically by Hodgkin and Huxley, and they still reproduce real axon data remarkably well.

2. Rate Functions α/β

Each gate's opening rate α and closing rate β are voltage-dependent functions, empirically fitted from voltage-clamp experiments:

function alphaM(V) { return 0.1 * (V + 40) / (1 - Math.exp(-(V + 40) / 10)); }
function betaM(V)  { return 4.0 * Math.exp(-(V + 65) / 18); }

function alphaH(V) { return 0.07 * Math.exp(-(V + 65) / 20); }
function betaH(V)  { return 1.0 / (1 + Math.exp(-(V + 35) / 10)); }

function alphaN(V) { return 0.01 * (V + 55) / (1 - Math.exp(-(V + 55) / 10)); }
function betaN(V)  { return 0.125 * Math.exp(-(V + 65) / 80); }

3. State Vector and Derivative Function

Pack [V, m, h, n] into one array and write a single function returning all four derivatives — this is what the RK4 integrator will call four times per step:

const params = {
  Cm: 1.0, gNa: 120, gK: 36, gL: 0.3,
  ENa: 50, EK: -77, EL: -54.4
};

function derivatives(state, Iext) {
  const [V, m, h, n] = state;
  const { Cm, gNa, gK, gL, ENa, EK, EL } = params;

  const iNa = gNa * m**3 * h * (V - ENa);
  const iK  = gK  * n**4 * (V - EK);
  const iL  = gL  * (V - EL);

  const dV = (Iext - iNa - iK - iL) / Cm;
  const dm = alphaM(V)*(1-m) - betaM(V)*m;
  const dh = alphaH(V)*(1-h) - betaH(V)*h;
  const dn = alphaN(V)*(1-n) - betaN(V)*n;

  return [dV, dm, dh, dn];
}

Initial resting state at V = −65 mV: the gates sit at their steady states, m₀ ≈ 0.05, h₀ ≈ 0.6, n₀ ≈ 0.32 — computed once from x_∞(−65) = α(−65)/(α(−65)+β(−65)) for each gate.

4. RK4 Integration

The Na⁺ activation gate m opens in under a millisecond — explicit Euler needs an impractically tiny timestep to stay stable there. Fourth-order Runge-Kutta (RK4) gets fourth-order accuracy per step, letting you use dt ≈ 0.01 ms comfortably:

function rk4Step(state, Iext, dt) {
  const add = (a, b, scale) => a.map((v, i) => v + b[i] * scale);

  const k1 = derivatives(state, Iext);
  const k2 = derivatives(add(state, k1, dt/2), Iext);
  const k3 = derivatives(add(state, k2, dt/2), Iext);
  const k4 = derivatives(add(state, k3, dt), Iext);

  return state.map((v, i) =>
    v + (dt/6) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])
  );
}

5. Driving With External Current

The interesting behaviour appears once you inject current. Below a threshold (roughly 6–7 μA/cm² for a brief pulse with these parameters), the perturbation just decays. Above it, the positive feedback loop between depolarisation and Na⁺ activation (m rises → more inward current → V rises faster → m rises further) fires a full, stereotyped spike — the all-or-nothing law:

function stimulusCurrent(t) {
  // 1 ms pulse of 10 µA/cm² starting at t = 5 ms
  return (t >= 5 && t < 6) ? 10 : 0;
}

let state = [-65, 0.05, 0.6, 0.32]; // [V, m, h, n]
let t = 0;
const dt = 0.01; // ms — small enough for stable Na+ kinetics

function tick() {
  const Iext = stimulusCurrent(t);
  state = rk4Step(state, Iext, dt);
  t += dt;
}

6. Rendering on Canvas

Push V(t) into a rolling buffer and draw it as a scrolling line plot, exactly like an oscilloscope trace. Overlay m, h, n scaled to [0,1] on a second axis so you can see the gates open and close in sync with the spike:

const ctx = canvas.getContext('2d');
const history = []; // { V, m, h, n } samples

function render() {
  ctx.clearRect(0, 0, canvas.width, canvas.height);
  const w = canvas.width, h = canvas.height;

  // Voltage trace: map V in [-90, 40] mV to canvas Y
  ctx.beginPath();
  history.forEach((s, i) => {
    const x = (i / history.length) * w;
    const y = h - ((s.V + 90) / 130) * h;
    i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  });
  ctx.strokeStyle = '#60a5fa';
  ctx.stroke();
}

7. Refractory Period and Repetitive Firing

Right after a spike, h (Na⁺ inactivation) is near zero and n (K⁺ activation) is elevated — the neuron is temporarily unable to fire again even with a strong stimulus. This is the absolute refractory period. As h recovers and n relaxes back down over a few milliseconds, the neuron enters a relative refractory period where only a stronger-than-normal stimulus can trigger another spike.

Drive the model with a sustained (rather than pulsed) current above a second, higher threshold and you get periodic repetitive firing — the frequency increases with Iext up to a saturation point set by how quickly the gating variables can recover, a simplified version of the firing-rate coding real neurons use to represent stimulus intensity.