Hodgkin-Huxley Model: Neuronal Excitability

• Hodgkin-Huxley Model is a foundational mathematical framework for describing action potential generation using nonlinear ODEs and voltage-dependent ionic conductances.
• It models fast and slow gating kinetics of sodium, potassium, and leakage channels to simulate the dynamics of excitable membranes.
• Extensions of the model include reduced dimensional formulations, parameter inference methods, and stochastic approaches to address biological variability.

The Hodgkin–Huxley (HH) model is a canonical mathematical framework for quantitatively describing the genesis and propagation of action potentials in excitable biological membranes, notably in neurons such as the squid giant axon. It consists of a coupled system of nonlinear ordinary differential equations accounting for voltage-dependent ionic conductances through sodium, potassium, and leakage channels, with explicit representation of the fast and slow gating kinetics underlying excitability. The model introduced a paradigm that remains foundational in computational neuroscience, electrophysiology, and the theory of excitable media.

1. Mathematical Formulation and Biophysical Principles

The original HH model is formulated as a four-dimensional ODE system capturing the membrane voltage $V(t)$ and three gating variables $m(t)$, $h(t)$ (sodium activation and inactivation), and $n(t)$ (potassium activation), each governed by first-order voltage-dependent kinetics. The system is:

$$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$

with canonical rate functions: $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$ where $C$ is the membrane capacitance, the $\bar g$ are maximal conductances, $E_{i}$ are reversal potentials, and $I_{\mathrm{ext}}$ is the externally applied current (Ori et al., 2018, Chialva et al., 2022).

Each gating variable $m(t)$0 reflects the probability of channel subunits being in the "open" state and modulates channel conductance multiplicatively.

2. Dimensional Reductions and Reduced Models

Despite the HH model's four-dimensionality, its key dynamical behaviors can be preserved in reduced descriptions. Several systematic reductions leverage the observation that $m(t)$1 is much faster than $m(t)$2 or $m(t)$3, allowing $m(t)$4. Moreover, $m(t)$5 and $m(t)$6 often evolve on similar, slower time scales and can be coupled via an approximate linear relation $m(t)$7 (Branco et al., 2024). This enables reduction to:

• Three-dimensional $m(t)$8 or $m(t)$9 models (steady-state for $h(t)$0 or $h(t)$1)
• Two-dimensional $h(t)$2 or $h(t)$3 models (collapsing $h(t)$4 onto $h(t)$5)
• Heuristic models like FitzHugh–Nagumo and Morris–Lecar that abstract the cubic-like $h(t)$6-nullcline geometry (Chialva et al., 2022)

Reduced models retain the ability to reproduce action potentials, threshold phenomena, and fundamental bifurcation structure. For propagation, reduced cable formulations recover the canonical scaling relation for conduction speed, $h(t)$7, where $h(t)$8 is effective diffusivity, $h(t)$9 is axial resistivity, and $n(t)$0 (Branco et al., 2024).

3. Parameter Robustness, Sensitivity, and $n(t)$1 Reduction

Simulation and analytical studies demonstrate that, although the HH model has $n(t)$210 independent parameters (capacitance, conductances, and kinetic rates), its qualitative excitability is largely determined by two composite variables:

• Structural coordinate $n(t)$3, quantifying the balance between inward and outward maximal conductance.
• Kinetic coordinate $n(t)$4, measuring the relative speed of restorative (e.g., $n(t)$5) versus excitatory (e.g., $n(t)$6) gates (Ori et al., 2018).

In the $n(t)$7 plane, parameter variations cluster into three sharply distinguished excitability regimes: nonexcitable, excitable, and oscillatory. The primary transition boundary is given approximately by $n(t)$8, with $n(t)$9 thresholds at $$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$0 (oscillatory/nonoscillatory) and $$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$1 (nonexcitable). This reduction captures 89% of model outcome variance, enabling analytic and machine-learning-based regime classification.

Sensitivity analysis in this context indicates sodium and potassium rate constants ($$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$2) and maximal conductances strongly determine excitability, whereas leak, sodium inactivation dynamics, and capacitance are much less influential (Ori et al., 2018).

4. Homeostatic Stabilization via Slow Inactivation

Long-timescale adaptation in excitability arises via slow inactivation processes, particularly affecting the sodium current. Modeling divides the total sodium-channel population into a fast-gating "available" pool $$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$3 (HH gating applies) and a slowly inactivated pool $$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$4 (timescale much longer than action-potential duration). The effective sodium conductance thus tracks $$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$5, and hence $$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$6. Activity-dependent feedback equations of the form

$$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$7

($$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$8 measures integrated depolarization)

and its projection into $$\begin{aligned} C\,\frac{dV}{dt} &= I_{\mathrm{ext}} - \Big[ \bar g_{\mathrm{Na}}\,m^3h\,(V-E_{\mathrm{Na}}) + \bar g_{K}\,n^4\,(V-E_{K}) + \bar g_L(V-E_L) \Big] \ \frac{dm}{dt} &= \alpha_m(V)\,[1-m] - \beta_m(V)\,m \ \frac{dh}{dt} &= \alpha_h(V)\,[1-h] - \beta_h(V)\,h \ \frac{dn}{dt} &= \alpha_n(V)\,[1-n] - \beta_n(V)\,n \ \end{aligned}$$9 yield a robust homeostatic mechanism. This slow inactivation drives $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$0 adaptively toward the nerve's excitable regime on the $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$1 diagram, restoring spike competency after structural or kinetic drift (Ori et al., 2018).

Simulations show that, in the excitable band, moderate fluctuations in $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$2 and $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$3 are re-centered by this feedback within seconds to minutes; in extreme regimes, excitability is robust to even order-unity fluctuations in $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$4.

5. Parameter Inference and Identification

Inference of HH model parameters from measured voltage traces is a major practical challenge, conventionally addressed by labor-intensive voltage-clamp protocols. Modern strategies employ derivative-based optimization such as Landweber regularized iteration. This approach computes gradient steps based on the Fréchet derivative and its adjoint, optimizing either conductance values $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$5 or gating exponents $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$6 to minimize the $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$7 discrepancy between observed and model-predicted voltage (Valle et al., 2019).

Empirically, this method achieves $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$8 parameter recovery accuracy in the presence of $$\begin{aligned} \alpha_m(V) &= 0.1\,(V+40)/(1 - e^{-(V+40)/10}), & \beta_m(V) &= 4\,e^{-(V+65)/18} \ \alpha_h(V) &= 0.07\,e^{-(V+65)/20}, & \beta_h(V) &= 1/(1 + e^{-(V+35)/10}) \ \alpha_n(V) &= 0.01\,(V+55)/(1 - e^{-(V+55)/10}), & \beta_n(V) &= 0.125\,e^{-(V+65)/80} \end{aligned}$$9–$C$0 noise, requiring tens to hundreds of thousands of forward/adjoint ODE integrations for convergence. The technique is extensible to models with more complex kinetic structures or temperature dependencies, and to higher-dimensional (multi-compartment) HH systems.

6. Extensions: Mechanistic, Reduced, and Quantum Models

Numerous directions extend or modify the HH formalism:

• Mechanistic alternatives replace empirical $C$1 fits with explicit energetic or exponential laws for channel activation, voltage-gating, and ion pumps (Deng, 2015), yielding models with fewer parameters and emergent hard thresholds (all-or-none dynamics).
• Electrodiffusive models build currents fundamentally from time-dependent Goldman–Hodgkin–Katz theory, yielding a resonator-type system with richer temperature and ionic concentration dependence, and correcting limitations of the original HH approach under constant-current drive or in rebound spiking (Stiles et al., 2019).
• Data-driven and hybrid methodologies embed learnable function approximators (e.g., shallow artificial neural networks) into the HH equations to directly estimate voltage-dependent rates from current-clamp recordings, enabling parameterization without extensive dedicated voltage-clamp data (Estienne, 2023).
• Quantum extensions interpret channel kinetics through the lens of quantum coherence and entanglement (e.g., selectivity-filter models), introduce quantum contributions to gating, and motivate quantum neuromorphic hardware analogs (Moradi et al., 2014, Gonzalez-Raya et al., 2018, Gonzalez-Raya et al., 2019).

7. Stochasticity, Channel Noise, and Robustness

The deterministic HH equations represent the mean-field limit of a large number $C$2 of stochastically gating ion channels. At finite $C$3, channel noise produces conductance fluctuations and stochasticity in spike timing, which can be captured by either direct Markov-chain kinetics or accurate stochastic differential equation (SDE) approximations (Goldwyn et al., 2011). The most accurate among these, such as the Fox–Lu system-size SDEs, reproduce both the correct stationary and time-dependent conductance variance and dynamics.

Stochastic HH models with noise-driven inputs can realize any mixture of deterministic spiking and quiescence over finite times, a property arising from local Hörmander hypoellipticity; this underscores how small channel-number fluctuations or synaptic variability can shift the neuron between firing and silent regimes in vivo (Höpfner et al., 2012). This noise-induced flexibility underlies both the observed variability of real neurons and the mathematical robustness of the excitable regime in the face of biophysical fluctuation.

---

References: (Ori et al., 2018, Chialva et al., 2022, Branco et al., 2024, Valle et al., 2019, Goldwyn et al., 2011, Höpfner et al., 2012, Moradi et al., 2014, Deng, 2015, Stiles et al., 2019, Estienne, 2023, Gonzalez-Raya et al., 2019).