Library-based numerical reduction of the Hodgkin–Huxley neuron for network simulation

Abstract

We present an efficient library-based numerical method for simulating the Hodgkin–Huxley (HH) neuronal networks. The key components in our numerical method involve (i) a pre-computed high resolution data library which contains typical neuronal trajectories (i.e., the time-courses of membrane potential and gating variables) during the interval of an action potential (spike), thus allowing us to avoid resolving the spikes in detail and to use large numerical time steps for evolving the HH neuron equations; (ii) an algorithm of spike-spike corrections within the groups of strongly coupled neurons to account for spike-spike interactions in a single large time step. By using the library method, we can evolve the HH networks using time steps one order of magnitude larger than the typical time steps used for resolving the trajectories without the library, while achieving comparable resolution in statistical quantifications of the network activity, such as average firing rate, interspike interval distribution, power spectra of voltage traces. Moreover, our large time steps using the library method can break the stability requirement of standard methods (such as Runge–Kutta (RK) methods) for the original dynamics. We compare our library-based method with RK methods, and find that our method can capture very well phase-locked, synchronous, and chaotic dynamics of HH neuronal networks. It is important to point out that, in essence, our library-based HH neuron solver can be viewed as a numerical reduction of the HH neuron to an integrate-and-fire (I&F) neuronal representation that does not sacrifice the gating dynamics (as normally done in the analytical reduction to an I&F neuron).

Appendix: Parameter values for the Hodgkin–Huxley equations

Parameter values or ranges and function definitions of the Hodgkin–Huxley model are as follows (Dayan and Abbott 2001):

$$\begin{array}{l} G_{\rm Na}=120{\rm mS/cm^2}, \ \ V_{\rm Na}=50{\rm mV}, \\ G_{\rm K}=36{\rm mS/cm^2}, \ \ V_{\rm K}=-77{\rm mV}, \\ G_{\rm L}=0.3{\rm mS/cm^2}, \ \ V_{\rm L}=-54.387{\rm mV}, \\ C =1{\rm \mu F/cm^2}, \ \ V_{\rm G}^{\rm E} =0{\rm mV}, \ \ V_{\rm G}^{\rm I} =-80{\rm mV}, \\ F^{\rm E}=0.05\sim0.1{\rm mS/cm^2}, \ \ S^{\rm E}=0.05\sim1.0{\rm mS/cm^2}, \\ F^{\rm I}=0.01\sim0.05{\rm mS/cm^2}, \ \ S^{\rm I}=0.05\sim1.0{\rm mS/cm^2}, \\ \sigma_{\rm r}^{\rm E}=0.5{\rm ms}, \ \ \sigma_{\rm d}^{\rm E}=3.0{\rm ms}, \\ \sigma_{\rm r}^{\rm I}=0.5{\rm ms}, \ \ \sigma_{\rm d}^{\rm I}=7.0{\rm ms}, \\ \alpha_m(V) = 0.1(V+40)/(1-\exp{(-(V+40)/10)}), \\ \beta_m(V) = 4 \exp{(-(V+65)/18)}, \\ \alpha_h(V) = 0.07 \exp{(-(V+65)/20)}, \\ \beta_h(V) = 1/(1+\exp{(-(35+V)/10})), \\ \alpha_n(V) = 0.01(V+55)/(1-\exp{(-(V+55)/10)}), \\ \beta_n(V) = 0.125 \exp{(-(V+65)/80)}. \end{array}$$