The FitzHugh-Nagumo Model

1. Origin and reduction from Hodgkin-Huxley

The Hodgkin-Huxley (HH) model (1952) accurately reproduces the squid giant axon action potential using four coupled ODEs: membrane voltage V, and gating variables m, h, and n for sodium and potassium channel dynamics. The model's complexity, however, makes geometric intuition difficult — a 4D phase space has no simple picture.

Richard FitzHugh (1961) noticed that in HH dynamics, m (the fast sodium activation) tracks V almost instantaneously, while h and n are slower and correlated. He proposed a 2D caricature using a fast membrane voltage variable v with cubic nullcline (replacing the fast V and m dynamics) and a slow recovery variable w (replacing the combined slow dynamics of h and n).

Independently, J. Nagumo, S. Arimoto, and S. Yoshizawa (1962) built an electronic analogue circuit implementing the same equations, demonstrating action potential firing with a tunnel diode. The combined model bears both names.

2. Equations and parameters

The standard FitzHugh-Nagumo equations are:

The cubic term v − v³/3 creates the N-shaped v-nullcline (dv/dt = 0), which is the geometric heart of the model. The slow linear recovery equation keeps trajectories bounded and generates the return path of the action potential.

The τ formulation

An alternative parameterisation uses time constants explicitly:

The ratio τ_w / τ_v ≈ 12.5 is the source of the sharp spike: the voltage variable v moves ~12 times faster than recovery w, producing the characteristic rapid upstroke followed by a slow plateau and return.

3. Phase plane analysis

The phase plane is the (v, w) plane. Each point in the plane represents a state of the system; time evolution traces a trajectory through that plane. Phase plane analysis asks: what are the qualitative trajectories for all initial conditions?

Key tools for phase plane analysis:

• Nullclines: curves where dv/dt = 0 (v-nullcline) and dw/dt = 0 (w-nullcline); trajectories cross the v-nullcline vertically and cross the w-nullcline horizontally.
• Fixed points: intersections of the two nullclines; where the system is at rest.
• Vector field: the direction of (dv/dt, dw/dt) at every point; arrows show where trajectories go.
• Limit cycles: closed trajectories that are stable (attracting); correspond to periodic oscillation (repetitive firing).
• Separatrix: boundary in phase space separating regions that converge to the fixed point from those that first execute a large excursion (the action potential).

The power of the FHN phase plane is that a single action potential appears as a trajectory that: starts near the resting fixed point, is pushed past the separatrix by I_ext, races up the fast cubic v-nullcline, slowly traverses the top plateau, falls down the right branch, slowly traverses the bottom plateau, and finally spirals back to rest.

4. Nullclines and fixed points

Setting each equation to zero defines the nullclines:

As I_ext increases, the v-nullcline shifts upward. The fixed point slides up the left branch of the cubic. When I_ext exceeds a critical value I_c ≈ 0.34, the fixed point moves onto the unstable middle branch of the cubic — triggering a supercritical Hopf bifurcation and the onset of repetitive firing.

Stability of the fixed point

The Jacobian at the fixed point (v*, w*) is:

5. Limit cycle and Hopf bifurcation

When I_ext exceeds the Hopf bifurcation threshold I_c, the stable fixed point becomes unstable and a stable limit cycle appears — the neuron fires action potentials periodically without external input. The Hopf bifurcation in FHN is supercritical: the cycle amplitude grows continuously from zero as I_ext crosses I_c.

The limit cycle grows with I_ext: higher current → shorter period → higher firing frequency. The f-I curve (firing frequency vs. injected current) starts at a finite frequency at the Hopf bifurcation (Type II excitability) or grows continuously from 0 (Type I excitability), depending on the bifurcation type.

Phase portrait

For I_ext < I_c, all trajectories spiral to the stable fixed point. For I_ext > I_c, all trajectories (except those starting on an unstable manifold) settle onto the limit cycle. The transition from excitable to oscillatory is sharp — it's a bifurcation, not a smooth change.

6. Classes I and II excitability

Neurons are classified by how their firing frequency depends on the applied current. This classification was proposed by Hodgkin (1948) and is directly readable from the FHN phase plane:

Class I — SNIC bifurcation

A Class I neuron transitions from rest to repetitive firing via a Saddle-Node on Invariant Circle (SNIC) bifurcation. The key signature: the firing frequency starts at zero and increases continuously with current — f(I_c) = 0. Class I neurons can integrate inputs over time (low-pass filter behaviour). FHN does not naturally produce SNIC; it can be mimicked by tuning parameters.

Class II — Hopf bifurcation

A Class II neuron transitions via a Hopf bifurcation. The firing frequency at the threshold is non-zero — the neuron jumps discontinuously from silence to firing at a minimum frequency. Standard FHN is Class II. Class II neurons are better resonators (bandpass filter behaviour) and respond most strongly to inputs near the natural Hopf frequency.

7. Coupled FHN oscillators

Two FHN oscillators coupled by a linear diffusive (gap junction) term exhibit rich synchronisation behaviour depending on coupling strength and the intrinsic frequency difference:

Spatially extended FHN — reaction-diffusion waves

Adding spatial diffusion of v across a 1D axon gives the FHN PDE:

The spatially-extended FHN is one of the primary models for cardiac arrhythmia. Spiral wave breakup in 2D FHN represents the transition from reentrant tachycardia to fibrillation. The parameter regions supporting spiral waves vs. planar waves vs. irregular patterns map directly to clinical arrhythmia classification.

8. JavaScript implementation

// FitzHugh-Nagumo ODE solver using RK4
// Plots v(t) and w(t); also draws phase plane trajectory

const params = {
  a:    0.7,
  b:    0.8,
  eps:  0.08,
  Iext: 0.5,    // try 0 (silent), 0.34 (threshold), 0.5–1.5 (repetitive firing)
};

const dt    = 0.1;    // ms time step
const T_end = 300;    // ms total simulation time
const tSteps = Math.floor(T_end / dt);

// Initial conditions (near rest)
let v = -1.2, w = -0.6;

// RK4 step for FHN
function fhnRHS(v, w) {
  const dv = v - (v * v * v) / 3 - w + params.Iext;
  const dw = params.eps * (v + params.a - params.b * w);
  return [dv, dw];
}

function rk4Step(v, w, dt) {
  const [k1v, k1w] = fhnRHS(v, w);
  const [k2v, k2w] = fhnRHS(v + dt/2*k1v, w + dt/2*k1w);
  const [k3v, k3w] = fhnRHS(v + dt/2*k2v, w + dt/2*k2w);
  const [k4v, k4w] = fhnRHS(v + dt*k3v, w + dt*k3w);
  return [
    v + dt/6 * (k1v + 2*k2v + 2*k3v + k4v),
    w + dt/6 * (k1w + 2*k2w + 2*k3w + k4w),
  ];
}

// Generate trajectory
const traj = [];
for (let i = 0; i < tSteps; i++) {
  traj.push([v, w]);
  [v, w] = rk4Step(v, w, dt);
}

// Draw v-nullcline and w-nullcline on phase plane canvas
function drawNullclines(ctx, W, H, scale) {
  const { a, b, Iext } = params;
  const vRange = [-2.5, 2.5];

  // v-nullcline: w = v - v³/3 + Iext
  ctx.strokeStyle = '#6366f1';
  ctx.lineWidth   = 2;
  ctx.beginPath();
  for (let px = 0; px <= W; px++) {
    const vv = vRange[0] + (px / W) * (vRange[1] - vRange[0]);
    const ww = vv - (vv*vv*vv)/3 + Iext;
    const x = px;
    const y = H/2 - ww * scale;
    px === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  }
  ctx.stroke();

  // w-nullcline: w = (v + a) / b
  ctx.strokeStyle = '#f59e0b';
  ctx.beginPath();
  for (let px = 0; px <= W; px++) {
    const vv = vRange[0] + (px / W) * (vRange[1] - vRange[0]);
    const ww = (vv + a) / b;
    const x = px;
    const y = H/2 - ww * scale;
    px === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  }
  ctx.stroke();
}

// Draw trajectory on phase plane
function drawTrajectory(ctx, traj, W, H, scale) {
  const vRange = [-2.5, 2.5];
  ctx.strokeStyle = '#34d399';
  ctx.lineWidth   = 1.5;
  ctx.beginPath();
  traj.forEach(([vv, ww], i) => {
    const x = (vv - vRange[0]) / (vRange[1] - vRange[0]) * W;
    const y = H/2 - ww * scale;
    i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  });
  ctx.stroke();
}

// Draw v(t) time series on a separate canvas
function drawTimeSeries(ctx, traj, dt, W, H) {
  ctx.strokeStyle = '#6366f1';
  ctx.lineWidth   = 1.5;
  ctx.beginPath();
  traj.forEach(([vv], i) => {
    const x = (i / traj.length) * W;
    const y = H/2 - vv * (H/6);
    i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  });
  ctx.stroke();
}
        

Detecting action potential spikes

// Count spikes: threshold crossing detection
function countSpikes(traj, threshold = 0) {
  let spikes = 0;
  let above  = false;
  for (const [vv] of traj) {
    if (!above && vv > threshold)  { spikes++; above = true; }
    if ( above && vv < threshold)  { above = false; }
  }
  return spikes;
}

// Build f-I curve: firing frequency vs. applied current
function buildFICurve(IValues, T = 500) {
  return IValues.map(Iext => {
    params.Iext = Iext;
    v = -1.2; w = -0.6;   // reset initial conditions
    const tSteps = Math.floor(T / dt);
    const traj = [];
    for (let i = 0; i < tSteps; i++) {
      traj.push([v, w]);
      [v, w] = rk4Step(v, w, dt);
    }
    const spikes = countSpikes(traj.slice(tSteps/2));  // ignore transient
    return (spikes / (T / 2)) * 1000;   // Hz (spikes per 1000 ms)
  });
}

// Usage: const Ivals = Array.from({length:50}, (_,i) => i*0.05);
//        const freqs = buildFICurve(Ivals);