FitzHugh-Nagumo: Spike Nucleation & Propagation of Excitation

1. From four equations to two

The Hodgkin-Huxley model tracks membrane voltage V plus three gating variables m, h, n — each with its own voltage-dependent rate functions. It reproduces real axon data almost perfectly, but its five-dimensional-ish structure hides the qualitative reason why a nerve spikes at all.

FitzHugh noticed that m (Na⁺ activation) is much faster than h and n, and that h and n move roughly in lockstep. Collapsing the fast variables into one "activator" and the slow variables into one "recovery" variable gives a 2-variable relaxation oscillator — the same mathematical family as the Van der Pol oscillator used to describe vacuum-tube circuits.

The payoff: the entire state of a FitzHugh-Nagumo (FHN) neuron is a single point (v, w) in a 2D plane, and its dynamics are just two nullcline curves crossing at an equilibrium. You can see excitability instead of only computing it.

2. The FitzHugh-Nagumo equations

The classical dimensionless form couples a fast activator v (playing the role of membrane voltage) with a slow recovery variable w:

where the parameters have this role:

• v — fast variable: excitation / membrane voltage analogue
• w — slow variable: recovery / net (Na⁺ inactivation + K⁺ activation)
• I_ext — external stimulus current
• ε (epsilon, typically 0.08) — time-scale separation; small ε ⇒ w moves much slower than v
• a, b (typically a = 0.7, b = 0.8) — shape and position of the recovery nullcline

The cubic term v − v³/3 is what makes the system excitable rather than simply damped: it creates an S-shaped (N-shaped) nullcline with three intersections with any nearly-flat line, which is the geometric root of the threshold behaviour below.

3. Phase plane and nullclines

A nullcline is the set of points where one derivative is zero. Setting dv/dt = 0 and dw/dt = 0 gives two curves in the (v, w) plane:

Where the two curves cross is a fixed point — the resting state, if it is stable. Because the v-nullcline is fast and the w-nullcline is slow, trajectories snap almost horizontally onto the cubic curve, then crawl slowly along it, jumping to the other branch whenever they reach a fold (the local max/min of the cubic). This fast-jump / slow-crawl pattern is the shape of an action potential: depolarisation (fast jump up), plateau (slow crawl), repolarisation (fast jump down), refractory recovery (slow crawl back).

4. Threshold and excitability

Because the v-nullcline is cubic, a small perturbation from rest that stays on the "near" branch simply decays back — the system is sub-threshold. But a perturbation that pushes v past the local fold point (the threshold) gets picked up by the fast dynamics and launched onto the far branch: a full spike fires, independent of exactly how big the extra push was, as long as it crossed the fold. This is the mathematical version of the all-or-nothing law of nerve firing.

Depending on parameters, the FHN system exhibits three regimes:

• Excitable — single stable resting point; a supra-threshold stimulus fires exactly one spike, then relaxes back.
• Oscillatory — the resting point becomes unstable (crosses a Hopf bifurcation) and the system fires spontaneously and periodically — a simple model of a pacemaker cell.
• Bistable — for some parameter ranges, two stable fixed points coexist, giving switch-like behaviour instead of spiking.

5. Spatial coupling: travelling waves

A single FHN unit is a point. Real tissue — nerve axon, cardiac muscle — is spatially extended, and neighbouring cells are electrically coupled (gap junctions or the axoplasm itself). Add a diffusion term to the fast variable and you get a reaction-diffusion PDE that supports self-sustaining travelling pulses:

The diffusion coefficient D sets how far the fast jump spreads to neighbours before they, too, cross threshold and fire — this is exactly how an action potential propagates along an axon, or how a depolarisation wave sweeps across cardiac tissue during a heartbeat. Because the recovery variable w lags behind, the tissue just behind the wavefront is refractory (temporarily unexcitable) — this is what prevents the wave from immediately reversing on itself and what allows stable spiral waves to form in 2D tissue, a mechanism believed to underlie cardiac arrhythmias such as atrial fibrillation.

6. Pseudocode

One explicit-Euler step of the spatially coupled system on a 2D grid:

function stepFHN(v, w, dt, dx):
  // v, w are 2D arrays (grid), same shape
  const a = 0.7, b = 0.8, eps = 0.08, D = 1.0

  for each cell (i, j):
    // 5-point Laplacian (diffusion of the fast variable)
    lap = (v[i+1][j] + v[i-1][j] + v[i][j+1] + v[i][j-1]
           - 4*v[i][j]) / (dx*dx)

    dv = D*lap + v[i][j] - v[i][j]**3/3 - w[i][j] + Iext(i, j)
    dw = eps * (v[i][j] + a - b*w[i][j])

    v_new[i][j] = v[i][j] + dv*dt
    w_new[i][j] = w[i][j] + dw*dt

  return v_new, w_new

Stability requires D·dt/dx² < 0.25 in 2D (standard explicit-diffusion CFL condition). For real-time browser simulation, dt ≈ 0.05–0.1 with a coarse grid (100–200 cells per side) is usually enough to see clean spiral waves and pulse annihilation on collision.