Picture a single brain cell as a tiny bucket that fills up with electric charge. Current trickles in, the level rises, and the instant it tips past the brim the cell fires a spike — a sharp pulse — then empties and starts filling again. What makes one neuron tick along steadily while another fires in bursts comes down to two things: how fast it refills, and how hard it snaps shut after each spike.
Two numbers tell the whole story. The membrane voltage $v$ is how full the bucket is (in millivolts), and a recovery variable $u$ is a slowly-building brake that drags the voltage back down. Eugene Izhikevich (2003) showed that just these two, with one compact rule, reproduce more than 20 real cortical firing patterns:
$$\frac{dv}{dt} = 0.04v^2 + 5v + 140 - u + I$$
Read it plainly: the cell charges itself (the $v^2$ term is the runaway that builds the spike), the brake $u$ holds it back, and $I$ is the current you inject. Once $v$ shoots past $+30$ mV the spike is finished, so we snap $v$ back down to $c$ and bump the brake up by $d$.
ONE WORKED NUMBER
Start a resting cell at $v=-65$ mV, $u=-13$, and inject $I=10$. Plug in: $0.04(-65)^2 + 5(-65) + 140 - (-13) + 10 = 169 - 325 + 140 + 13 + 10 = +7$ mV/ms. The level is climbing, so the cell will rise toward threshold and fire — predicted from a single line of arithmetic.
▸ What the Four Knobs Do
The brake has its own rule, $\frac{du}{dt} = a(bv - u)$. Here $a$ sets how quickly the brake reacts (small $a$ = sluggish = bursting), $b$ sets how tightly the brake tracks the voltage, $c$ is the voltage the cell drops to right after a spike, and $d$ is the extra brake added on every spike — which is why a steadily-driven cell gradually slows down, an effect called spike-frequency adaptation. Mathematically the $0.04v^2$ term is the normal form of a saddle-node bifurcation: the universal shape any neuron takes right at its firing threshold, which is why so few numbers buy so many behaviours. The sliders in the panel above are exactly $a$, $b$, $c$, $d$ and $I_\text{ext}$.
TRY THIS IN THE SIM ABOVE
1. Load the RS preset, then drag $d$ from 8 up toward 20 and watch each spike push the next one further apart — adaptation happening in real time. 2. Switch to FS ($a=0.1$) and notice the firing barely slows: a fast brake fully resets between spikes. 3. Pick IB and step $I_\text{ext}$ up slowly to watch the cell move from silent, to a burst-then-tonic train, to rapid tonic firing.
§ 03
Equation Derivation — Izhikevich Model (2003)
▸ The Izhikevich System
Published by Eugene Izhikevich in 2003, this 2D model reproduces 20+ cortical firing patterns with a simple quadratic integrate-and-fire equation plus a linear recovery variable.
$$\boxed{\frac{dv}{dt} = 0.04v^2 + 5v + 140 - u + I}$$
$$\boxed{\frac{du}{dt} = a(bv - u)}$$
$$\text{if } v \geq 30 \text{ mV}: \quad v \leftarrow c, \quad u \leftarrow u + d$$
The four parameters a, b, c, d are sufficient to reproduce all known cortical firing patterns:
Recovery timescale (small a → slow recovery → bursting)
ms⁻¹
Burst vs. tonic
b
Sensitivity of u to subthreshold fluctuations
—
Resonance vs. integrator
c
After-spike reset voltage
mV
AHP depth
d
After-spike reset of u (step increase)
—
Adaptation strength
0.04v²+5v+140
Quadratic approximation to Na⁺ activation + K⁺ rectification
mV/ms
Spike initiation
STEP 1 — Quadratic Integrate-and-Fire Origin
The term 0.04v² + 5v + 140 approximates the fast Na⁺ current near threshold. Near threshold (~−55 mV), the Na⁺ activation curve is approximately quadratic in V. The Izhikevich model fits the coefficients (0.04, 5, 140) to match the subthreshold V-nullcline of the Hodgkin-Huxley model when reduced to 2D. The quadratic term provides the positive feedback (threshold crossing) that the linear LIF model lacks.
STEP 2 — Recovery Variable u
u combines Na⁺ inactivation and K⁺ activation into a single variable (as in FHN). The equation du/dt = a(bv − u) drives u toward bv with time constant 1/a. When v rises during a spike, u rises, providing negative feedback (recovery). The parameter b controls how tightly u tracks v subthreshold — large b creates resonance (u follows v oscillations). After each spike, u jumps by d (additional recovery force), creating spike-frequency adaptation.
STEP 3 — Why This Model Is Computationally Superior
Izhikevich (2004) compared models on two axes: biological plausibility (number of firing patterns reproduced) vs. computational efficiency (spikes per second simulated). HH: biologically accurate but slow (1000× more computation than Izhikevich). LIF: fast but only 1 pattern. Izhikevich: 20+ patterns at LIF-like speed (~1000× faster than HH). This makes the Izhikevich model the de facto standard for large-scale brain simulations (100,000+ neurons in real time).
STEP 4 — Euler Integration (Midpoint Method)
The Izhikevich model is typically integrated with a midpoint (modified Euler) method at dt = 0.5 ms: dv = 0.5*(0.04v²+5v+140−u+I); du = 0.5*a*(b*v−u); v+=dv; v+=dv; u+=du; u+=du. This doubles each step but maintains stability. The quadratic term can cause v to diverge rapidly at threshold — the reset rule (v ≥ 30 → v=c, u+=d) prevents this and represents the spike peak + repolarisation.
The neuron is not at rest (dv/dt = 7 mV/ms > 0) — it will fire. V rises toward 30 mV (the spike cutoff), triggers reset: v←−65, u←−13+8=−5. After reset, dv/dt is recalculated with the same v=−65 but the new u=−5: 0.04(−65)² + 5(−65) + 140 − (−5) + 10 = 169 − 325 + 140 + 5 + 10 = −1 mV/ms. Because u sits higher, the −u term now pulls dv/dt down (from +7 to −1), so the cell must charge up from rest again → slower next spike (adaptation). This is spike-frequency adaptation: each spike increases u by d=8, progressively slowing subsequent spikes.
▸ References
[IZH03]
Izhikevich — Simple model of spiking neurons, IEEE Trans. Neural Netw. 14(6):1569–1572, 2003
[IZH04]
Izhikevich — Which model to use for cortical spiking neurons? IEEE Trans. Neural Netw. 15(5):1063, 2004
[IZH07]
Izhikevich — Dynamical Systems in Neuroscience, MIT Press, 2007. Free: dynamicalsystems.org
🔬 SimulationWhat do the four parameters a, b, c, d actually control — what changes when I move each slider?▼
a (recovery timescale, 0.001–0.3): small a → u recovers slowly → long bursts or no adaptation. Large a → fast recovery → fast spike-frequency adaptation. Try a=0.02 (RS, slow) vs a=0.1 (FS, fast). b (sensitivity, −0.5 to 0.5): positive b → u increases with v → more adaptation; negative b → u decreases with v → rebound potential (LTS). b near 0.25 creates resonance. c (reset mV, −80 to −40): determines where V resets after a spike. c=−65 (RS) → normal AHP; c=−50 (CH) → quick reset → chattering. d (reset increment, −20 to 20): each spike adds d to u. Large d → strong after-spike adaptation (u increases, slows firing). d=8 (RS) produces gradual frequency adaptation; d=4 (IB) produces bursting.
Key takeaway: a sets temporal dynamics, b sets subthreshold coupling, c sets spike reset location, d sets cumulative adaptation strength. Four numbers, 20+ biological patterns.
🧠 ConceptualWhy is the Izhikevich model considered the "best of both worlds" compared to HH and LIF?▼
Izhikevich (2004) explicitly quantified this: on a plot of computational cost (spikes/second simulated on a PC) vs. biological plausibility (number of firing patterns reproduced), all models fall on a Pareto frontier. HH: highest plausibility (~5 patterns well-reproduced, biophysically accurate), highest cost (~1000× slower than Izhikevich). LIF: lowest cost, only 1 firing pattern (regular spiking). Izhikevich: 20+ firing patterns, 1000× faster than HH, only 2 ODEs + 4 parameters. It also reproduces phenomena LIF cannot: bursting, chattering, spike-frequency adaptation, rebound spikes, bistability. The tradeoff: no biophysical parameter interpretation (a,b,c,d have no direct channel correlates), no AP shape.
Key takeaway: Izhikevich is the efficiency Pareto optimum — maximum biological richness per unit of computational cost. Use it for large-scale network simulations where pattern diversity matters but not AP shape.
🌍 AppliedHow is the Izhikevich model used in large-scale brain simulations and neuromorphic hardware?▼
The Izhikevich model is used in the largest cortical simulations ever published: Izhikevich & Edelman (2008) simulated 100 billion synaptic connections of a thalamocortical model on a 27-processor cluster. Modern large-scale simulations (Blue Brain Project, Human Brain Project, Allen Brain) use Izhikevich-type models for most excitatory neurons where computational efficiency is critical. In neuromorphic hardware, the Intel Loihi chip implements a generalised Izhikevich-type neuron (quadratic integrate-and-fire with adaptation) rather than LIF, because the additional pattern diversity is worth the minor extra cost. The SpiNNaker 2 chip (University of Manchester) runs Izhikevich neurons natively on its ARM cores at biological real-time speed.
Key takeaway: The Izhikevich model is the current practical standard for large-scale brain simulation. Its 4-parameter flexibility enables network models to reproduce the diverse firing pattern statistics of real cortex.
💡 Non-ObviousWhat produces the burst-then-tonic transition in Intrinsic Bursting (IB) neurons?▼
In IB neurons (a=0.02, b=0.2, c=−55, d=4): the first stimulus produces a burst of 3–5 spikes in rapid succession, then the neuron transitions to tonic (regular) spiking. The mechanism: at spike onset, u starts low (after reset, u accumulates). During the first burst, each spike adds only d=4 to u — small enough that u doesn't rise much. But the reset voltage c=−55 mV (higher than RS c=−65) means the neuron immediately re-reaches threshold quickly, firing again. After several spikes, u has accumulated enough to slow threshold approach, and the bursting gives way to tonic. In the phase portrait: the trajectory orbits the limit cycle close to the v-nullcline initially (burst), then settles into a regular limit cycle orbit (tonic). This burst-then-tonic pattern is a hallmark of layer 5 pyramidal neurons and is important for the initiation of cortical UP states.
Key takeaway: IB burst→tonic transition is controlled by the interplay of c (reset position, determines burst intra-spike interval) and cumulative u accumulation. It reproduces the signature of L5 pyramidal cells that initiate cortical network states.
📐 ComputationalWhy does the Izhikevich model use Euler integration instead of RK4 — isn't RK4 always more accurate?▼
The Izhikevich model is designed to use a modified Euler (midpoint) method at dt = 0.5 ms, which Izhikevich specifically validated against the full HH model. The rationale: (1) The reset rule at v=30 mV is a discontinuous event — any smooth ODE solver (RK4, CVODE) must handle this event detection separately, eliminating their advantage over Euler for the near-threshold regime; (2) The quadratic term 0.04v² produces rapid divergence near threshold, but since we reset at v=30 mV before any unbounded growth, stability is not an issue; (3) At dt=0.5 ms, modified Euler gives <2% error in ISI timing for all 20 firing patterns. Using RK4 at dt=0.5 ms gives negligible improvement for 4× more computation. The correct approach: Euler + dt=0.5 ms + exact event detection at v=30 mV.
Key takeaway: For the Izhikevich model, Euler at dt=0.5 ms is optimal — RK4 adds computation without accuracy benefit because the dominant numerical error comes from the discontinuous reset rule, not the ODE integration method.
🎓 DeepHow does the Izhikevich model relate to normal forms of neuron bifurcations?▼
The quadratic term in the Izhikevich v-equation is precisely the normal form of the saddle-node bifurcation: near threshold, the membrane voltage dynamics behave as \(\dot v \approx \mu + v^2\), where μ encodes the distance from threshold. The quadratic integrate-and-fire (QIF) model \(\dot v = v^2 + I\) is the exact normal form, and the Izhikevich model adds the coefficients and offsets to match HH units (mV, ms). This means Izhikevich neurons are inherently Type I excitable (saddle-node/SNIC onset bifurcation), with f–I curve starting from 0 Hz. The linear recovery variable u extends this to a 2D system capable of Type II behaviour by creating negative eigenvalue feedback. The theta model \(\dot\theta = 1 - \cos\theta + I(1+\cos\theta)\) is the exact normal form of SNIC; Izhikevich is its Cartesian-coordinates approximation.
Key takeaway: The quadratic term 0.04v² is the normal form of the saddle-node bifurcation — the Izhikevich model is not ad hoc but is theoretically grounded as the universal QIF near-threshold approximation.
🧠 ConceptualWhat is the difference between Fast Spiking (FS) and Regular Spiking (RS) neurons in cortical circuits?▼
RS neurons (a=0.02, d=8) are excitatory pyramidal neurons: they show spike-frequency adaptation (firing rate decreases over a sustained stimulus), relatively depolarised reset (c=−65 mV), and a gradual f–I curve. FS neurons (a=0.1, d=2) are inhibitory parvalbumin (PV) interneurons: they fire at very high rates (100–200 Hz) with minimal adaptation (d=2), narrow spikes, and constant firing frequency under sustained input. The fast recovery (a=0.1) means u returns to baseline quickly between spikes, preventing adaptation. In cortical circuits, RS-FS interactions generate gamma oscillations (PING): RS cells drive FS cells; FS cells inhibit RS cells; the E-I loop creates gamma. The ratio of RS-to-FS input determines gamma frequency and amplitude. Try both presets and observe the dramatic difference in adaptation.
❌"The Izhikevich parameters a, b, c, d have direct biophysical interpretations — a is the membrane time constant, b is the conductance ratio, etc."
✅The Izhikevich parameters are phenomenological fitting parameters, not biophysical quantities. They were chosen to reproduce firing pattern statistics, not to match ion channel parameters. a is a timescale ratio (related to the slow K⁺ current time constant), b controls the slope of the u-nullcline, c is the post-spike reset (related to AHP depth), and d is the after-spike u increment (related to adaptation current amplitude). You cannot directly set a=1/τ_K or b=g_K/g_Na. To add biophysical meaning, use the AdEx model (which has τ_m, a, τ_w, b parameters with clear biophysical interpretations).
❌"The Izhikevich model can reproduce any firing pattern by tuning a, b, c, d — it is a universal neuron model."
✅The Izhikevich model is very flexible but NOT universal. It cannot reproduce: (1) Subthreshold oscillations that depend on h-current (I_h) — the 2D structure lacks the right eigenvalue structure; (2) Graded responses (non-spike outputs from dendrites); (3) Calcium-dependent bursting with precise burst envelope shapes; (4) Class III excitability (neurons that respond only to onset of inputs, never tonically). For models requiring these, use the adaptive exponential (AdEx) or multi-compartmental Hodgkin-Huxley models. The 20 patterns it reproduces are all spiking patterns of cortical neurons.
📖 Izhikevich (2007) — DSN, Ch. 8; Naud et al. (2008) — AdEx model
❌"The threshold in the Izhikevich model is 30 mV — the neuron fires when v reaches 30 mV."
✅30 mV is the numerical spike peak cutoff, not the biophysical threshold. The actual firing threshold (where positive feedback becomes unstable) is at v ≈ −55 to −50 mV, where dv/dt begins to grow faster than u can compensate. The Izhikevich equations are structurally equivalent to a QIF neuron — there is an unstable fixed point (threshold) and the reset at v=30 mV merely terminates the numerical spike before v diverges to infinity. Setting the cutoff to +30 mV (rather than, say, +100 mV) is a convention that produces realistic-looking spike heights.
❌Using dt = 1 ms instead of 0.5 ms for Izhikevich integration — produces spike timing errors of 5–15% and incorrect firing patterns (RS neurons look like FS).
✅Izhikevich (2003) validates the model with dt = 0.5 ms. The modified Euler method runs the half-step update twice: for (let i=0;i<2;i++){v+=0.5*(0.04*v*v+5*v+140-u+I); u+=0.5*a*(b*v-u);} then check if v≥30 → reset. At dt=1 ms, the quadratic spike initiation is poorly resolved — ISIs are underestimated, adaptation appears faster, and IB neurons may appear as RS. Always use dt ≤ 0.5 ms.
🔍 Why: The quadratic term 0.04v² changes rapidly near threshold — coarse time steps miss the inflection point and shift the effective threshold.
❌Not checking for v ≥ 30 BEFORE the ODE update — allowing v to overshoot to very large values (1000+ mV) before reset, causing NaN or Inf.
✅Check the reset condition AFTER each integration step: if (v >= 30) { v = c; u += d; recordSpike(t); }. Do NOT check before updating — you need the step to push v over 30. If you use large dt and v overshoots 30 mV significantly (e.g., v=150 mV), apply a retrospective reset: record the spike at the interpolated crossing time t_cross = t − dt × (v_prev − 30)/(v − v_prev). This improves spike timing accuracy by sub-ms interpolation.
🔍 Why: The quadratic term causes rapid acceleration near threshold — large dt means v can jump from 25 mV to 300 mV in one step without hitting the 30 mV reset check.