Channel-noise-induced critical slowing in the subthreshold Hodgkin-Huxley neuron

The dynamics of a spiking neuron approaching threshold is investigated in the framework of Markov-chain models describing the random state-transitions of the underlying ion-channel proteins. We characterize subthreshold channel-noise-induced transmembrane potential fluctuations in both type-I (integrator) and type-II (resonator) parametrizations of the classic conductance-based Hodgkin-Huxley equations. As each neuron approaches spiking threshold from below, numerical simulations of stochastic trajectories demonstrate pronounced growth in amplitude simultaneous with decay in frequency of membrane voltage fluctuations induced by ion-channel state transitions. To explore this progression of fluctuation statistics, we approximate the exact Markov treatment with a 12-variable channel-based stochastic differential equation (SDE) and its Ornstein-Uhlenbeck (OU) linearization and show excellent agreement between Markov and SDE numerical simulations. Predictions of the OU theory with respect to membrane potential fluctuation variance, autocorrelation, correlation time, and spectral density are also in agreement and illustrate the close connection between the eigenvalue structure of the associated deterministic bifurcations and the observed behavior of the noisy Markov traces on close approach to threshold for both integrator and resonator point-neuron varieties.

Figure 1

Markov-chain model (kinetic scheme) of the Hodgkin-Huxley potassium channel with four activation $n$ gates. There are five states in total, each characterized by a different number of open gates, and eight possible state transitions. The only fully open and conducting state is $n_4$. State subscripts indicate the number of open $n$ gates and state transition rates are given next to the associated arrow.

Figure 2

Markov-chain model (kinetic scheme) of the Hodgkin-Huxley sodium channel with three activation $m$ gates and one inactivation $h$ gate. There are a total of 8 states, each characterized by a different number of open gates, and 20 possible state transitions. The only fully open and conducting state is $m_3{h}_1$. State subscripts indicate the number of open $m$ and $h$ gates and state transition rates are given next to the associated arrow.

Figure 3

Markov-chain stochastic simulations of a Hodgkin-Huxley neuron configured as (a) a type-I integrator and (b) a type-II resonator for a sequence of four distinct subthreshold current-density values $I_{\text{DC}}$ (in $\mu \text{A/cm}{}^2$) for each neuron type. See Table 1 for parameter values and Figs. 1 and 2 for Markov-chain kinetics. Blue (gray) traces show membrane potential, $V_m$, fluctuations induced by stochastic channel-state transitions; fluctuations become larger and more spectrally colored as spiking threshold is approached from below. In both cases, the top trace captures a single spiking event, shown in thick red (thick gray), that has been cropped to highlight baseline fluctuations. For each run, all 14 state variables (five ${\mathrm{K}}^{+}$ channel states, eight ${\text{Na}}^{+}$ channel states, and membrane potential) were initialized at the deterministic steady state and then iterated using the Chow and White [26] implementation of the Gillespie stochastic simulation algorithm (SSA) [25] and sampled at deterministic time intervals $\Delta t=0.005$ ms for 10 s of simulated time; 1 s extracts are shown here. Traces in (b) have been offset vertically to reduce overlap (offset in mV = [0, 4, 8, 12] counting from bottom to top). These Markov simulations are computationally demanding (matlab run times $\sim 64$ h per 10 s trace on an Intel Xeon 2.80-GHz CPU).

Figure 4

Distribution densities for subthreshold voltage fluctuations in (a) type-I integrator and (b) type-II resonator Hodgkin-Huxley neurons, comparing Markov [blue (gray)] and SDE [orange (light gray)] numerical simulations with linear Ornstein-Uhlenbeck (OU) theory (thick black). The Markov and SDE probability densities are estimated from their respective $V_m$ time series by computing transmembrane potential histograms and then dividing by $N\Delta V$ with $N=2001$ being the number of bins and $\Delta V$ the bin width. For each $I_{\text{DC}}$ simulation, the thick black reference curve is the OU prediction that the subthreshold fluctuations will be Gaussian distributed with zero mean. Vertical scale markers carry units of ${\text{mV}}^{-1}$ and the total area under each curve is a dimensionless unity. Distributions from numerical simulations accurately follow theoretical prediction, even at the highest currents that result in occasional action potential generation. The broadening of the distributions corresponds to increased fluctuation amplitude as instability threshold is approached from below. (The “thorns” riding on the top of the maximum-current Markov and SDE histograms are a signature artifact generated by the zeroing of action-potential events detected in the time series.)

Figure 5

Autocorrelation functions for subthreshold transmembrane potential fluctuations in the (a) type-I integrator and (b) type-II resonator Hodgkin-Huxley neurons, comparing Markov [blue (gray)] and SDE [orange (light gray)] numerical simulations with linear Ornstein-Uhlenbeck (OU) theory (thick black). For both neuron types, there is a dramatic surge in variance (zero-lag autocorrelation) and temporal persistence (time for the decay-envelope to drop to $1/e$ of the zero-lag maximum) as spiking threshold is approached from below. We note the excellent agreement between the linear OU prediction and the Markov and SDE simulations; although the SDE numerics run about $200\times$ faster than the Markov SSA codes, both algorithms produce output that is statistically indistinguishable. In (a), the bottom three traces are calibrated to the ${\text{0.001-mV}}^2$ vertical scale marker, while the top trace has been rescaled to a ${\text{0.01-mV}}^2$ interval. The Markov correlations were computed from the full 10-s Markov recordings (see Fig. 3) by averaging autocorrelations from five 2000-ms nonoverlapping segments; the SDE traces were computed from a corresponding set of SDE simulation runs (not shown); the OU curves are predictions from OU theory [see Eqs. (7) and (8)]. Spike events occurred occasionally (about 12 spikes per 10 s) at the highest level of $I_{\text{DC}}$ drive for the integrator neuron; these events were detected and suppressed by setting the transmembrane potential fluctuation to zero for the spike duration.

Figure 6

Power spectral densities for subthreshold voltage fluctuations in the (a) type-I integrator and (b) type-II resonator Hodgkin-Huxley neurons, comparing Markov [blue (gray)] and SDE [orange (light gray)] numerical simulations with linear Ornstein-Uhlenbeck (OU) theory (thick black). These curves are the Fourier-domain equivalent of the autocorrelation functions plotted in Fig. 5, but the Markov and SDE spectra were computed directly from Markov (Fig. 3) and SDE numerical simulations (not shown) via the matlab spectrogram function ($10000$ ms recordings at 200-kHz sample rate, 100 ms epochs with 50% overlap, and top-hat window); the OU theoretical spectra are linear predictions from Eqs. (7) and (9). Vertical scale markers carry units of ${\text{mV}}^2/\text{Hz}$. Total area under each curve estimates the $V_m$ time-series fluctuation variance (i.e., zero-lag correlation peak in Fig. 5). Agreement between numerically simulated and theoretical spectra is generally excellent, except for the topmost trace of (a) where the 100-ms epoch lengths are too short to adequately capture the energy at zero frequency. (Much longer epochs can retrieve a more accurate low-frequency peak but at the cost of degraded signal-to-noise from the diminished number of spectral samples available to be averaged.)

Figure 7

Divergence of correlation decay times as a function of proximity to instability threshold in the (a) type-I integrator and (b) type-II resonator Hodgkin-Huxley neurons. Curves show characteristic times ${\tau }_i=1/{\alpha }_i$ (in ms), with $i=1,2,3$ and $\alpha \equiv -\text{Re}\left(\text{eig}\right)$ being the decay rate predicted by the first 3 (of 12) Ornstein-Uhlenbeck (OU) eigenvalues. The Markov and SDE correlation times are measured as the time required for the autocorrelation functions of their $V_m$ fluctuations to decay to $1/e$ of the zero-lag peak; for (b), we revise this to measure the time for the oscillations envelope to decay to $1/e$. Evidently these correlation times follow power-law growth trends with (a) $\tau \sim {\epsilon }^{-0.5}$ for integrator and (b) $\tau \sim {\epsilon }^{-1}$ for resonator, where $\epsilon \equiv (I_{\text{DC}}^{\text{crit}}-I_{\text{DC}})/I_{\text{DC}}^{\text{crit}}$ (see Table 1 for $I_{\text{DC}}^{\text{crit}}$ values) measures proximity to the instability point as determined by a sign-change in the real part of the dominant eigenvalue in the deterministic equations.

Figure 8

Testing for power-law growth of membrane-voltage fluctuation variance as a function of proximity to instability threshold in (a) type-I integrator and (b) type-II resonator Hodgkin-Huxley neurons. Ornstein-Uhlenbeck (OU) theory predicts that voltage fluctuation power will scale asymptotically as (a) ${\epsilon }^{-0.5}$ for integrator and (b) ${\epsilon }^{-1}$ for resonator, where $\epsilon \equiv (I_{\text{DC}}^{\text{crit}}-I_{\text{DC}})/I_{\text{DC}}^{\text{crit}}$ is a dimensionless measure of the distance from the critical point. Squares (Markov) correspond to the four sets of currents shown in previous figures, and asterisks (channel-SDE) show spike-suppressed fluctuation variances for stochastic simulations over a much wider range of $I_{\text{DC}}$ values; down-triangles (Markov) and up-triangles (SDE) show fluctuation power when spikes are not suppressed in the variance calculations. The collapse in spike-suppressed fluctuation power for $\epsilon \lesssim {10}^{-0.6}$ in (b) arises because resonator spikes tend to occur in bursts, leading to prolonged intervals of zeroed activity, whereas in the integrator (a), intermittent spikes remain brief, solitary events. Prior to emergence of spikes, both Markov and SDE simulations follow OU predictions closely. (Note that the region of concordance would extend further to the left if the membrane area—and hence the number of Na and K channels—were increased.)