Hodgkin–Huxley Model

The simulation numerically integrates the four coupled ODEs of the Hodgkin-Huxley model in real time using the 4th-order Runge-Kutta method with dt = 0.025 ms by default. The main canvas shows membrane voltage V(t) as a waveform — each upswing is a true action potential generated by the interplay of Na⁺ activation (m gate), Na⁺ inactivation (h gate), and K⁺ activation (n gate). The I_ext slider is your stimulus: below ~6.5 μA/cm², the neuron is subthreshold; above it, you see repetitive firing whose rate increases with current. The ḡ_Na and ḡ_K sliders change the peak conductances — reducing ḡ_Na below ~60 mS/cm² eliminates action potentials entirely. The Noise σ slider adds stochastic current that can trigger spontaneous spikes or irregular firing.

The HH model captures the macroscopic ionic currents responsible for the action potential — specifically how voltage-dependent conductances create a self-reinforcing depolarization followed by active repolarization. It successfully explains threshold, the all-or-nothing law, the absolute and relative refractory periods, repetitive firing, and anode break excitation. What it leaves out is equally important: it does not model individual channel stochasticity (it treats conductance as a continuous, averaged variable), it does not include calcium channels (critical for bursting and synaptic release), it does not model spatial propagation along the axon (only a point neuron or voltage-clamped patch), and it uses squid axon parameters that differ from mammalian neurons. Modern extensions include the Traub model (CA3 pyramidal cell) and the Blue Brain / Allen Brain cell models.

The HH model is the foundation of virtually all quantitative neurology and neural engineering. In epilepsy research, HH-type models explain how hyperexcitability arises from altered Na⁺ channel inactivation (SCN1A mutations in Dravet syndrome). In cardiac medicine, the Luo-Rudy and O'Hara-Rudy heart cell models are HH-type models — the same formalism predicts long-QT arrhythmias and guides drug design targeting K⁺ channels. In deep brain stimulation (DBS) for Parkinson's disease, HH-based network models help optimize stimulation parameters. Cochlear implant pulse sequence design uses HH-type auditory nerve models to maximize speech intelligibility. In AI and neuromorphic computing, Intel's Loihi chip implements LIF neurons inspired by HH simplifications, and spiking neural networks (SNNs) aim to match the energy efficiency of biological spike-based computation.

Anode break excitation (also called "post-inhibitory rebound spike") is one of the most counterintuitive phenomena in neuroscience: apply a strong hyperpolarizing (negative) current for 50–100 ms, then abruptly stop — and the neuron fires a spike, even though no depolarizing current was applied. The mechanism: prolonged hyperpolarization drives the h gate (Na⁺ inactivation) toward 1 (fully de-inactivated) while simultaneously driving n (K⁺ activation) toward 0 (K⁺ channels close). When the current turns off, the membrane springs back toward rest, but now the Na⁺ channel is unusually "primed" (high h, low n) — the inward Na⁺ current overwhelms the reduced K⁺ current and triggers a spike. To see this in the simulation: set I_ext to −5 μA/cm², run for ~60 ms, then snap I_ext back to 0.

The HH equations are stiff ODEs: the fastest time constant (\(\tau_m\) ≈ 0.1 ms at V = −40 mV) is ~100× smaller than the slowest (\(\tau_n\) ≈ 10 ms near rest). Forward Euler requires \(dt \ll \tau_{\min}\) for stability, meaning \(dt \lesssim 0.01\) ms. With \(dt = 0.1\) ms, Euler produces catastrophic instability — the gating variables blow up to values outside [0,1] and V diverges to ±∞ within milliseconds. RK4 uses four derivative evaluations per step, achieving 4th-order accuracy: the local truncation error is O(dt⁵) versus O(dt²) for Euler. In practice, RK4 with \(dt = 0.025\) ms gives accurate, stable results for the HH system. For production simulations, CVODE (implicit BDF method) is preferred because it adaptively controls step size on stiff ODEs. In this simulation, we enforce dt ≤ 0.05 ms and warn if a larger value is set.

The transition from rest to repetitive firing as I_ext increases is a bifurcation — a qualitative change in the system's dynamics. For the standard HH parameters, the onset of repetitive firing occurs via a subcritical Hopf bifurcation: the resting state loses stability and a large-amplitude limit cycle (the action potential train) appears abruptly, with a finite minimum firing frequency (~50 Hz). This makes HH a "Type II" excitable system (Rinzel-Ermentrout classification): the f-I curve has a discontinuous onset. In contrast, a SNIC (saddle-node on invariant circle) bifurcation produces Type I excitability with a continuous f-I curve starting from 0 Hz (as in the Morris-Lecar model). The distinction has profound implications: Type II neurons tend to synchronize in networks, while Type I neurons encode a wider dynamic range of input intensities. These bifurcation structures are invisible to the time-domain simulation but are revealed by continuation methods (AUTO, XPPAUT) that trace equilibria and limit cycles as parameters vary.

Immediately after an action potential, the Na⁺ inactivation gate h has been driven to nearly 0 (Na⁺ channels are fully inactivated), and the K⁺ activation gate n is near its maximum (K⁺ channels are wide open, hyperpolarizing the membrane). During this absolute refractory period (~1–2 ms), no amount of current can generate another action potential because the inward Na⁺ current is completely suppressed. The duration is set by the recovery time of h, governed by \(\tau_h(V)\) — roughly how long it takes h to recover from ~0 back toward 0.6. The subsequent relative refractory period (2–10 ms) occurs while h is partially recovered and n is still elevated; a spike can be triggered but requires a larger-than-normal I_ext and produces a smaller, slower action potential.