A guided climb from an everyday picture up to the full phase-plane equations. No calculus is needed for the first steps.
START — the mousetrap idea (no math)
A neuron behaves like a mousetrap or a toilet flush: a gentle nudge does nothing — it just settles back. But push past a tipping point and it snaps all the way, releasing one full "spike," then slowly resets before it can fire again. The FitzHugh-Nagumo (FHN) model is the simplest cartoon that captures this rule — ignore small pushes, fire all-out for a big one, then recover.
BUILD — two numbers, one fast, one slow
FHN tracks just two quantities: a fast one, $V$ (the voltage — how excited the cell is), and a slow one, $w$ (the recovery — a brake that builds up after firing). The fast variable obeys $\frac{dV}{dt}=V-\frac{V^3}{3}-w+I$. The cubic term $-V^3/3$ is the whole trick: for small $V$ it pushes the cell back toward rest, but in a middle range it pushes the other way — that runaway push is the all-or-nothing spike. The current $I$ is simply how hard you are poking the cell.
WORKED NUMBER — crossing threshold
Start at rest, $V\approx-1.2$. A tiny stimulus lifts $V$ a little and the cubic drags it straight back — no spike. But cross the tipping point near $V\approx 0$ and the cubic flips sign: $V$ shoots up toward $+2$, and only then does the slow brake $w$ catch up and haul it back down. One poke, one clean spike.
DEEPEN — the phase plane (AP / intro-college)
Pair the fast equation with the slow one, $\frac{dw}{dt}=\frac{1}{\tau}(V+a-bw)$, where $\tau\gg1$ makes $w$ lag far behind $V$ — that timescale gap is what carves each trajectory into a sharp spike plus a long recovery. In the $(V,w)$ plane, the cubic V-nullcline $w=V-V^3/3+I$ and the straight w-nullcline $w=(V+a)/b$ cross at the fixed point. If they meet on the cubic's down-sloping middle branch, that point is unstable and the state circles forever — a limit cycle, i.e. repetitive firing. The sliders map straight onto these symbols: $I$ slides the cubic up and down, $a$ and $b$ shift and tilt the line, and $\tau$ stretches the spike.
TRY THIS — in the sim above
① Load the Excitable preset and set $V_0$ just below, then just above, the threshold — one nudge dies out, the other launches a full spike. ② Raise $I$ from about $0.4$ upward: the resting dot turns into an endless loop in the Phase Plane — a Hopf bifurcation, visible as the firing rate lifting off zero in the Bifurcation tab. ③ Drag $\tau$ from $3$ up to $40$ and watch the action potential go from a quick wobble to a long, lazy spike.
§ 03
Equation Derivation — FitzHugh-Nagumo (1961)
▸ The FitzHugh-Nagumo System
Proposed by Richard FitzHugh (1961) and independently by Nagumo, Arimoto & Yoshizawa (1962), the FHN model is a 2D reduction of the Hodgkin-Huxley model that preserves all qualitative dynamical features (excitability, oscillations, threshold) while being analytically tractable.
$$\boxed{\frac{dV}{dt} = V - \frac{V^3}{3} - w + I_{\text{ext}}}$$
$$\boxed{\frac{dw}{dt} = \frac{1}{\tau}(V + a - b\,w)}$$
The 4D HH system (V,m,h,n) is reduced to 2D via two observations: (1) m relaxes ~10× faster than V — approximate m ≈ m∞(V) (instantaneous activation); (2) h and n evolve together with roughly opposite changes — define w = 1 - h ≈ n (single recovery variable). Result: a 2D fast-slow system (V fast, w slow).
STEP 2 — Van der Pol Oscillator Connection
The FHN V-equation \(\dot V = V - V^3/3 - w + I\) is related to the Van der Pol oscillator. The cubic term V³/3 provides the nonlinear "negative resistance" that drives the excursion, analogous to the Na⁺ channel's positive feedback. The w-equation provides the linear restoring force that returns the system to rest.
STEP 3 — Three Operating Regimes
Depending on I: (1) Excitable (I < 0): stable fixed point, sub-threshold perturbations decay; large perturbations trigger a single action-potential-like excursion; (2) Oscillatory (0 < I < ~1.2): Hopf bifurcation, limit cycle, repetitive firing; (3) High-amplitude oscillation (large I): limit cycle with large amplitude. The threshold is the unstable branch of the V-nullcline.
STEP 4 — Phase Plane Analysis
The V-nullcline (cubic: dV/dt=0) and w-nullcline (linear: dw/dt=0) intersect at the fixed point(s). If the fixed point lies on the middle (negative-slope) branch of the cubic, it is unstable → limit cycle (Hopf bifurcation). If on the outer branches, it is stable → excitable or depolarised rest. The τ parameter controls how far the trajectory overshoots before being pulled back (larger τ → larger, longer action potential).
STEP 5 — RK4 Integration
Both equations are integrated with RK4 at dt = 0.05 (dimensionless time units). The FHN system is not stiff (no fast singularities like αm in HH), so dt can be larger than in HH. The singularity-handling and gating variable clamping of HH are not needed. Noise is added as Gaussian white noise σξ(t)dt to the V equation.
🔬 SimulationWhat does the Phase Plane tab show, and what are the nullclines?▼
The Phase Plane tab shows the trajectory of the system state (V,w) over time as an animated parametric curve. The cyan V-nullcline (cubic S-curve) is the set of (V,w) where dV/dt = 0 — on this curve, V is momentarily unchanging. The magenta w-nullcline (straight line) is where dw/dt = 0. Their intersection(s) are the fixed points. The trajectory moves: (1) horizontally rightward where it is above the V-nullcline (dV/dt > 0), leftward below; (2) vertically upward above the w-nullcline, downward below. During an action potential excursion, the trajectory makes a large counter-clockwise loop — the "spike" is the right half of this loop, the "repolarisation" is the return. A limit cycle is a closed loop that the trajectory circles indefinitely.
Key takeaway: Every feature of the action potential is visible in the phase plane — threshold is the unstable fixed point, the AP is the large loop, and the refractory period is the slow return along the lower branch of the V-nullcline.
🧠 ConceptualWhy is FHN considered a "reduced" model — what exactly was reduced from HH?▼
The Hodgkin-Huxley model has 4 variables (V, m, h, n) and 12 nonlinear functions (α, β rates for each gate). FHN reduces this to 2 variables and 2 smooth functions. The reduction exploits timescale separation: m activates ~10× faster than V, so it can be slaved to V (quasi-static approximation m ≈ m∞(V)). The two slow variables h (Na⁺ inactivation) and n (K⁺ activation) both contribute to "recovery" and move in opposite directions during an AP — they can be combined into a single recovery variable w. The result preserves: threshold/excitability, all-or-nothing firing, refractory period, repetitive firing, phase-locking to stimuli, and anode break excitation. What is lost: the exact AP shape, temperature dependence, absolute timing precision, and subthreshold oscillations near rest.
Key takeaway: FHN is "right for the right reasons" — it preserves the topological structure of excitable and oscillatory dynamics while eliminating all ion-channel-specific detail.
🌍 AppliedWhere is the FHN model used in medicine, engineering, and neurotechnology?▼
FHN is used wherever the qualitative behaviour of excitable/oscillatory systems matters more than quantitative biophysical accuracy. In cardiac modelling, FHN-type equations describe excitable wave propagation in heart tissue — spiral wave breakup (fibrillation) and reentry circuits are studied with FHN reaction-diffusion models. In neural prosthetics, FHN networks model how populations of neurons synchronise (relevant for DBS optimisation). In chemical pattern formation, the Belousov-Zhabotinsky reaction obeys FHN-like equations — it produces spirals and target patterns identical to cardiac arrhythmias. In ecological modelling, predator-prey cycles (Lotka-Volterra) have similar fast-slow structure. The universality of FHN means any qualitative analysis carries over across all these domains.
Key takeaway: The FHN model is a canonical model of excitability that appears across biology, chemistry, and physics — any phenomenon that shows threshold, all-or-nothing response, and recovery follows FHN-like dynamics.
💡 Non-ObviousWhat is a Hopf bifurcation, and where does it occur in the FHN model?▼
A Hopf bifurcation is a qualitative change in dynamics where a stable fixed point loses stability and a limit cycle (oscillation) is born. In FHN, as I_ext increases, the fixed point moves from the outer (stable) branch of the V-nullcline toward the middle (negative-slope) branch. When the fixed point crosses onto the middle branch, the Jacobian trace changes sign from negative to positive — the fixed point becomes an unstable spiral, and a small limit cycle appears around it. This is a supercritical Hopf bifurcation if the limit cycle appears continuously with small amplitude; it is subcritical if a large-amplitude cycle appears abruptly. In FHN with standard parameters, the transition to oscillation is supercritical — you can see this in the Bifurcation tab: firing rate increases continuously from 0 Hz as I crosses the bifurcation point.
Key takeaway: The Hopf bifurcation in FHN is the mathematical equivalent of the neuron's firing threshold — it marks the transition from silence to repetitive firing and determines the Type II excitability character (non-zero onset frequency).
📐 ComputationalWhy is FHN easier to simulate than HH, and what time units does it use?▼
FHN uses dimensionless time and voltage units — V and w are normalised (typical range −2 to +2), not millivolts. Time is in dimensionless units where one unit ≈ 10–20 ms in physical time. This makes the system non-stiff: the fastest time constant is ~1 (τ_V ≈ 1/(1-V*²) near rest), and the slowest is τ (typically 12.5). The stiffness ratio is only ~12.5, compared to ~100 in HH. Forward Euler is stable for dt < 0.1 in FHN; RK4 with dt = 0.05 gives excellent accuracy. There are no singularities in the rate functions (no α_m at V=−40 mV issue), no temperature-dependent factors, and no need to clamp variables to [0,1]. This makes FHN ideal for educational purposes and for large-scale network simulations where speed matters.
Key takeaway: FHN is 10–100× faster to simulate than HH and requires no special numerical handling — but remember that "1 time unit" in FHN ≈ 10–20 ms in real time, so axis labels should clarify the unit system being used.
🎓 DeepWhat is the connection between FHN, Van der Pol, and the Bonhoeffer model?▼
The FHN model is a modified Van der Pol oscillator. The Van der Pol equation \(\ddot x - \mu(1-x^2)\dot x + x = 0\) describes a limit cycle oscillator with nonlinear damping. Writing it as a 2D system gives \(\dot v = v - v^3/3 - w\) and \(\dot w = v\) — almost identical to FHN but without the separate timescale τ and without the recovery parameters a, b. FitzHugh added the recovery variable parameters to model the excitable (non-oscillating) regime, making it applicable to resting neurons that only fire in response to stimuli. The Bonhoeffer-Van der Pol (BVP) model is an earlier related model. Nagumo's electronic circuit implementation of FHN (using tunnel diodes and inductors) was the first electronic neural model and directly inspired neuromorphic computing. The modern LIF neuron can be seen as a linearised, threshold-rule version of FHN.
Key takeaway: FHN is the biological specialisation of the Van der Pol oscillator — adding recovery parameters transforms a pure oscillator into an excitable system that can model both resting and periodically firing neurons.
🔬 SimulationWhat does the τ parameter control, and why must it be > 1?▼
The τ parameter (default 12.5) sets the timescale separation between V (fast) and w (slow): the w equation is \(\dot w = (1/τ)(V+a-bw)\), so larger τ → w changes more slowly. τ must be greater than 1 (typically 5–20) to maintain the fast-slow structure that produces action potential-like excursions. If τ ≈ 1, V and w evolve on the same timescale — the system loses the slow recovery phase and action potential-like waveforms disappear, replaced by small, fast oscillations or dampened spirals. If τ is too large (> 100), the recovery becomes negligible and the system fires indefinitely without returning to rest (like HH without K⁺ channels). In the simulation, observe: small τ (1–3) → small, fast spirals; τ = 12.5 → clear AP-like spikes; τ > 50 → very long, slow action potentials.
Key takeaway: τ ≫ 1 is the fundamental assumption of FHN — without this timescale separation, the reduction from HH is invalid and the system loses its excitable/oscillatory character.
§ 04 Best Resources
Izhikevich — Dynamical Systems in Neuroscience — dynamicalsystems.org/ds/ds.pdf (free, Ch. 4–6)
❌"The V variable in FHN is in millivolts and directly comparable to measured membrane potentials."
✅FHN uses dimensionless normalised variables. V ≈ 0 in FHN does not correspond to 0 mV — it corresponds roughly to −50 mV (near threshold). V = −1.2 (rest) ≈ −65 mV, V = +1.5 (AP peak) ≈ +40 mV. To compare with experiments, you must apply a linear rescaling: V_mV ≈ V_FHN × 50 + (−30) mV (approximate). Similarly, one FHN time unit ≈ 10–20 ms. Forgetting this makes FHN simulations appear to produce instantaneous APs with zero-millivolt resting potentials.
📖 FitzHugh (1961), §2; Izhikevich — DSN Ch. 4
❌"A stable fixed point in FHN always means the neuron is resting and will never fire."
✅A stable fixed point means only that small perturbations decay back to rest — the system is excitable, not permanently silent. A suprathreshold perturbation (or current step) can push the state over the threshold (unstable middle branch of V-nullcline), triggering a full action potential excursion before returning to the stable rest. This is the all-or-nothing property in phase space: below threshold, trajectory returns directly; above threshold, trajectory makes the full large-amplitude loop first. The FHN preset "Excitable" demonstrates this — click "AP state" to push V over threshold.
❌"The FHN model is a simplification of HH, so it is less accurate and less useful than HH."
✅FHN is a different tool, not a worse one. For understanding the qualitative dynamical mechanisms (bifurcations, synchrony, excitability types, network oscillations), FHN is often superior to HH because it is analytically tractable — fixed points, nullclines, bifurcation points can all be computed exactly. Many fundamental results in computational neuroscience (Type I vs Type II excitability, Rinzel-Ermentrout classification, synchrony theory) were first derived using FHN/2D models. HH is superior for quantitative matching to experimental recordings and pharmacological predictions.
❌Using dt = 0.5 or larger (thinking FHN is non-stiff so Euler is fine with any step size), causing visible phase plane trajectory distortion and incorrect limit cycle size.
✅While FHN is far less stiff than HH, RK4 with dt = 0.05 gives trajectories indistinguishable from exact. With forward Euler at dt = 0.5, the trajectory spirals outward (Euler adds artificial energy) and the limit cycle amplitude and period are wrong by ~15–30%. Use RK4 with dt ≤ 0.1 dimensionless time units for accurate phase portraits.
🔍 Why: "Non-stiff" is misinterpreted as "any step size is fine"; Euler's energy-adding error is small per step but accumulates over an oscillation period.
❌Computing the w-nullcline as w = V + a - bw (not solving for w): w_null = V + a - b*w — a circular definition that is undefined unless solved explicitly.
✅The w-nullcline is dw/dt = 0, so \((1/τ)(V+a-bw) = 0\) gives \(w = (V+a)/b\). Code: w_null = (V + a) / b. The V-nullcline requires more care: dV/dt = 0 gives \(w = V - V³/3 + I\), which is straightforward. Both are functions of V, parametrically defining a curve in (V,w) space. Never use the current value of w to compute the nullcline — it is a geometric object in phase space, independent of the trajectory.
🔍 Why: Students confuse the nullcline equation (a geometric curve) with the time derivative equation (a differential equation evaluated along the trajectory).
❌Forgetting that τ appears in the w equation only: writing dV = (V - V³/3 - w + I)/tau instead of only applying τ to dw/dt.
✅The standard FHN form is: dV/dt = V - V³/3 - w + I (no τ) and dw/dt = (1/τ)(V + a - b*w) (τ only in w equation). Some formulations write τ dV/dt = ... and dw/dt = ... — in this case τ appears in the V equation. Check which convention is used. This simulation uses the first (τ in w equation only), consistent with FitzHugh's original paper and most textbooks.
🔍 Why: Different textbooks use different conventions for where τ appears; copying code from mixed sources introduces this error.
§ 05 References
FitzHugh (1961) — Biophys. J. 1:445–466
Izhikevich — Dynamical Systems in Neuroscience, MIT Press, 2007. Ch. 4–6. Free: dynamicalsystems.org