Exponential Time Differencing Algorithm For Pulse-coupled Hodgkin-Huxley Neuronal Networks

The exponential time differencing (ETD) method allows using a large time step to efficiently evolve the stiff system such as Hodgkin-Huxley (HH) neural networks. For pulse-coupled HH networks, the synaptic spike times cannot be predetermined and are convoluted with neuron's trajectory itself. This presents a challenging issue for the design of an efficient numerical simulation algorithm. The stiffness in the HH equations are quite different between the spike and non-spike regions. Here, we design a second-order adaptive exponential time differencing algorithm (AETD2) for the numerical evolution of HH neural networks. Compared with the regular second-order Runge-Kutta method (RK2), our AETD2 method can use time steps one order of magnitude larger and improve computational efficiency more than ten times while excellently capturing accurate traces of membrane potentials of HH neurons. This high accuracy and efficiency can be robustly obtained and do not depend on the dynamical regimes, connectivity structure or the network size.


Introduction
The Hodgkin-Huxley (HH) model [7,13,16] is a classical neuron model, originally proposed to describe the behaviors of action potentials of the squid's giant axon. It provides a useful mechanism that accounts for the detailed generation of action potentials and the existence of the absolute refractory periods. It also serves as the foundation for other neuron models such as the one that can describe the behaviors of bursting and adaption [21]. However, the HH equations are so complicated that it is difficult to study its properties analytically such as the Hopf bifurcation and chaotic dynamics [1,10,12,19]. Therefore, its investigation often relies on numerical simulations, for example, by the Runge-Kutta (RK) methods.
There are several difficulties to design an efficient and accurate numerical algorithm for the HH neural network, especially when the network size is large. First, when an HH neuron driven by external input generates an action potential (the interval of action potential is called spike period in this work), the HH neuron equations become stiff. Regular RK methods have to use very small time step to satisfy the requirement of numerical stability [2,3,10,17]. This small time step will significantly increase the computational cost when studying long time behavior of large-scale HH networks such as chaotic attractor dynamics or collecting reliable statistical information of HH neurons such as the distribution of inter-spike-intervals.
The exponential time differencing (ETD) method [5,8,14,15,17,20] is proposed for efficient simulation of stiff ordinary differential equations (ODEs). The basic idea is to decompose the ODEs into a linear stiff part and a nonlinear non-stiff part. Then, the linear stiff part can be solved by using the integrating factor method, while the nonlinear non-stiff part can be approximated by numerical quadrature. A second-order ETD method for HH neural networks has been proposed in a recent work [3], which allows using a large time step to raise computational efficiency. However, the HH neural network considered in previous works is a special case where each neuron is driven by constant input and the synaptic conductance is described by a smooth function [9]. For realistic situations, the neurons are generally driven by stochastic spike input and the interaction term is usually modeled by a Dirac delta function (pulsecoupled), while the spike-induced conductance dynamics are modeled by an α function [11,23,24]. These make the system become non-smooth and event-driven, while providing challenges for the design of efficient numerical simulation algorithms. For instance, it is impossible to predetermine the synaptic spike times since they are convoluted with neurons' trajectories themselves. As a result, one has to evolve the HH network by ignoring the spike interactions among neurons and then use spike-spike interaction to amend the neurons' trajectories at the end of the time step [4,11]. Without a careful recalibration for the neuronal spikes, the numerical algorithm often suffers from the issue of instability or relatively low numerical accuracy.
In this work, we first provide a second-order ETD method (ETD2) to evolve a pulse-coupled HH neural network driven by stochastic spike input. Note that the stiffness of HH equations are quite different between the spike and non-spike periods, and we find that the ETD method may introduce a relatively large error in the membrane potentials in the non-spike period if using the same time step as that in the spike period. We then design an adaptive ETD2 method (AETD2) that using different decompositions of the linear and nonlinear parts in spike and non-spike periods. In addition, for the situation where neurons generate spikes in the time step, the effects of the spikes are carefully recalibrated in our AETD2 method to achieve a second-order numerical accuracy. Our AETD2 method is capable of using a large time step, while achieving the same high accurate traces of membrane potentials of each neuron as the RK2 method using a very small time step. It can improve computational efficiency more than one order of magnitude compared with the RK2 method. This high numerical accuracy and computational efficiency can be achieved over a wide range of dynamical regimes and does not depend on the network connectivity and size.

The model
The dynamics of the ith neuron of an HH neural network is governed by where C is the cell membrane capacitance, V i is the membrane potential, m i , h i , and n i are gating variables for sodium and potassium currents, respectively [6]. The parameters V N a , V K , and V L are the reversal potentials for the sodium, potassium, and leak currents, respectively, G N a , G K and G L are the corresponding maximum conductances. The form of α z and β z are set as [6]: is given by where G E i and G I i are excitatory and inhibitory conductances, respectively, V E G and V I G are the corresponding reversal potentials. The dynamics of conductances can be explicitly expressed as where s il is the spike time of the feedforward Poisson input with strength f and rate ν, and τ jl is the lth spike time of the jth neuron. The spike-induced conductance change is described as where σ d and σ r are slow decay and fast rise time scale, respectively, and Θ(·) is the Heaviside function. Each neuron is either excitatory or inhibitory and its coupling strength is labeled by its type E or I, respectively. For example, S E ij (S I ij ) is the coupling strength from the jth excitatory (inhibitory) neuron to its postsynaptic ith neuron. The model parameters are C = 1µF·cm −2 , V N a = 50 mV, .0 ms, σ I r = 0.5 ms, and σ I d = 7.0 ms [6]. The voltage V i evolves continuously according to Equations (1) and (2). When it reaches the threshold V th = −50 mV, we say the ith neuron generates a spike at this time, say τ il . Then it will trigger its postsynaptic jth neuron's conductance change in the form of S Q ji H(σ Q d , σ Q r , t − τ il ), Q = E, I. For the ease of discussion about our algorithm design, we consider an all-to-all connected network with S Q ij = S/N , where Q = E, I, S is the coupling strength and N is the total number of neurons in the network. Note that our algorithm can be easily extended to networks with more complicated connectivity structure.

Runge-Kutta method
Without loss of generality, we consider the second-order Runge-Kutta method (RK2) as the benchmark and compare it with the ETD methods. We first introduce the RK2 method to evolve the HH neural network with a fixed time step ∆t, for example, to evolve the system from time t = t k = k∆t to t = t k+1 = (k + 1)∆t. Since the synaptic spike times in [t k , t k+1 ] can not be predetermined, one has to evolve the network without considering synaptic spike interactions and reconsider their effects by using spike-spike interactions at the end of time step [4,11].
Due to the pulse-coupled dynamics in Equation (6), the numerical accuracy may be very low if the spike timing is not well estimated. For example, suppose that a presynaptic spike fired att between t k and t k+1 . If one simply assigns it to be the end of time step t k+1 , then the error of the spike-induced conductance change is Therefore, the error with the magnitude of ∆t will be introduced when the system evolves to t = t k+1 . We now solve the above issue arising from the pulse-coupled dynamics to achieve a second-order numerical accuracy. First, we evolve the HH neural network without considering the feedforward and synaptic spikes during the time interval [t k , t k+1 ]. Then, at time t = t k+1 , some neuron's voltage may be above the threshold, i.e., generating a spike, say neuron i, if V i (t k ) < V th and V i (t k+1 ) ≥ V th . The spike time, say τ il , can be estimated following the idea proposed in Ref. [11,22]. The neuron's membrane potential during the time interval can be approximated by a linear interpolation: and the spike time τ il can be estimated by solving the equation: Since there may be some neurons firing and some feedforward spikes emitting during the time interval and they will induce the conductance change, one should update the conductance by considering the effects of both the feedforward and synaptic spikes. When they are not considered, the conductance, denoted byG Q , is and it should be recalibrated as for j = 1, 2, ..., N . A detailed algorithm of the RK2 method is given in Algorithm 1.

Algorithm 1: RK2 algorithm
Input: an initial time t k and feedforward input times {s il } Output: Solutions at time t k+1 for i = 1 to N do Solve the HH eqautions for the ith neuron without considering spike input using RK2 scheme.
Estimate the spike time τ il by Equation (9).
We show that the above algorithm can indeed achieve a second-order numerical accuracy as follows. If there are no feedforward or synaptic spikes, then all the dependent variables are infinitely differentiable and the RK2 method can achieve an error of order O(∆t 2 ). For the time step that contains feedforward or synaptic spikes, an error of order O(∆t) is introduced in the conductance with the form of G Q −G Q , Q = E, I. Nevertheless, the dependent variables of V, m, h, and n can have an error of order O(∆t 2 ). The synaptic spike times are estimated by a linear interpolation and also have an error of order O(∆t 2 ). After recalibration shown in Equations (12) and (13), the conductances can achieve numerical accuracy of second-order at the end of the time step. Therefore, all the dependent variables V, m, h, n, G E , and G I have an error of order O(∆t 2 ) (see below for verification of numerical results).

Exponential time differencing method
Exponential time differencing method is proposed to solve the stiff problem in differential equations by decomposing the system into a linear stiff term and a nonlinear non-stiff term [5,14,17,20]. Following this idea, we propose the ETD schemes for HH Equations (1) and (2) below. As illustrated in Algorithm 1, each neuron in the HH network is evolved independently and their conductances are recalibrated at the end of time step. Thus, one can first derive an ETD scheme for a single HH neuron and then consider the spike interactions among neurons, and obtain an ETD scheme for the numerical evolution of an HH neural network.
Consider the evolution of a single HH neuron from t k to t k+1 . We use z k to represent z(t k ) for z = V, m, h, n of this neuron and rewrite Equations (1) and (2) as where and Here, F z (t, V, m, h, n) is actually a function of t, V , and z for z = m, h, n, but we write in this way for ease of illustration. Note that the linear coefficient c z in Equation (14) is a constant value in the kth time step [t k , t k+1 ] and is updated with respect to k. Multiplying Equation (14) by an integrating factor e −czt and taking integral from t k to t k+1 , we obtain (19) for z = V, m, h and n.
The essence of the ETD method is to derive proper approximations to the above integration. We take a second-order ETD formulae with RK time stepping as follows [5]. Let and approximate F z during the time interval [t k , t k+1 ] by for z = V, m, h, and n, where F z,k represents F z (t k , V k , m k , h k , n k ). Substituting the above approximation into Equation (19) yields the ETD2 scheme which is given by for z = V, m, h, and n. The procedure of the ETD2 algorithm for an HH neural network is similar to that of the RK2 algorithm given in Algorithm 1, but the RK2 scheme is replaced by the ETD2 scheme in Equation (22).

Adaptive Exponential time differencing method
The ETD2 method can indeed use a large time step to improve computational efficiency, but we find that it will introduce relatively large error in the trajectories of neurons' membrane potentials and even lead to the missing of action potentials (see below for numerical results). In addition, the number of the missed action potentials in the ETD2 method can grow with the increase of time steps compared with the RK2 method using a small time step. Thus, it is important to design an efficient but also reliable ETD method to solve this issue.  (14) for the non-spike period.
As shown in Figure 1A, the slope of voltage has a very large value when the neuron generates an action potential (spike period) and quickly reduces to a value around zero in the non-spike period until the next spike time. Therefore, the stiffness of HH equations is quite different between spike and non-spike periods. In the non-spike period, the slope of voltage is almost zero, while the linear and nonlinear part in Equation (14) have a much larger absolute value and nearly cancel each other out as shown in Figure 1B. Therefore, the decomposition in Equation (14) may not be appropriate in the non-spike period since both the linear and nonlinear parts become stiff while the summation of them is indeed non-stiff. Based on this, we propose a different decomposition in the non-spike period from that in the spike period: taking c z = 0 and the whole right hand side of Equation (14) as the nonlinear part. For such a decomposition, the ETD2 scheme reduces to the RK2 scheme in the non-spike period.  Figure 2: Illustration of the AETD2 method. After neuron 1 fires a spike, we use the ETD2 scheme to evolve the HH equations for neuron 1 during the spike period indicated by the red vertical lines, while the HH equations for neuron 2 is evolved using RK2 scheme since neuron 2 is in the non-spike period. The starting point of the spike period can be determined as the spike time and it lasts for the following about 3.5 ms. The circles and dots indicate the time nodes where we use the ETD2 and RK2 schemes, respectively.
Based on the above observation, we give our AETD2 method for HH neural network as following: each neuron is evolved using ETD2 scheme if it is in the spike period and use the reduced ETD2 scheme, the RK2 scheme, otherwise, as shown in Figure 2. The starting point of the spike period can be determined as the spike time and the interval of spike period can be determined as 3.5 ms which is sufficient large to cover the highly stiff region of the spike. Detailed AETD2 algorithm is given in Algorithm 2.

Algorithm 2: AETD2 algorithm
Input: an initial time t k , feedforward input times {s il } Output: Solutions at time t k+1 for i = 1 to N do Solve the HH equaions for the ith neuron without considering spike input: if The ith neuron is inside spike period then use ETD2 scheme The ith neuron spiked in [t k , t k+1 ]. Estimate the spike time τ il by Equation (9).

RESULTS
We consider an all-to-all connected network of 80 excitatory and 20 inhibitory neurons driven by Poisson feedforward input. For the ease of illustration, we choose the Poisson input strength f = 0.06 mS·cm −2 and input rate ν = 300 Hz, and the coupling strength between neurons are chosen as S = 0.2 mS·cm −2 throughout this work, unless indicated otherwise. However, our algorithm can be applied to HH neural networks under a variety of dynamical regimes. First, we verify the second-order numerical accuracy by performing convergence tests. A high precision solution is obtained by using RK2 method with a sufficiently small time step ∆t = 1 × 10 −6 ms and is denoted by a superscript "high". It is compared with the solutions computed by the RK2, ETD2, and AETD2 methods with various values of larger time steps ∆t = 2 −4 , 2 −5 , ..., 2 −12 ms which is denoted by a superscript "∆t". Errors of membrane potentials at final run time T = 2000 ms and the last spike time of each neuron are computed: where τ il * indicates the last spike time of the ith neuron during the run time interval. As shown in Figure  3, if one naively assigns the end of time step as the spike times in the RK2 method, the numerical accuracy of the membrane potentials and spike times can only be of the first-order. In contrast, if one determines the spike times by linear interpolation and recalibrate the conductances accordingly, all the RK2, ETD2, and AETD2 methods can achieve a second-order numerical accuracy. In addition, we find that the ETD2 method has much larger error compared with the RK2 and AETD2 methods using the same time step as shown in Figure 3. When using a time step larger than ∆t = 2 −6 = 0.0156 ms, the ETD2 method performs even worse than the naive RK2 method. The underlying reason is that the HH equations are almost non-stiff in the non-spike period, but the decomposition in Equation (14) induces a relatively large stiffness for the nonlinear term as discussed previously. We next discuss the numerical performance of our AETD2 method and compare it with other different numerical methods. As shown in the top panel of Figure 4, the AETD2 method with large time steps (maximum time step ∆t = 0.277 ms) can obtain the same high accuracy in membrane potentials as the RK2 method using a very small time step ∆t = 0.01 ms. The bottom panel of Figure 4 shows the raster plots (neuron index versus its spike time) of the spike events in the network. It can be seen that the spike times are well captured by the AETD2 method with large time steps. The ETD2 method is proved to be unconditionally stable for the HH system in Ref. [3]. Hence, it can use a much larger time step compared with AETD2 method in principle. However, as shown in Figure  5, the ETD2 method is highly inaccurate when the time step ∆t = 0.277 ms is used (the maximum time step in AETD2 method). There are some spikes missing as shown in Figure 5A and Figure 5B, in contrast, the AETD2 method does not encounter such issues when using the same time step. Figure  5C shows the relative error in the mean firing rate (the average number of synaptic spikes per unit time) between the RK2 and ETD2 (AETD2) methods over different values of coupling strength. It can be seen that the ETD2 method can achieve only one digit of numerical accuracy while the AETD2 method can robustly achieve more than two digits of numerical accuracy when the time step ∆t = 0.277 ms is used in both methods. Therefore, the ETD2 method has worse numerical accuracy in terms of voltage traces and firing rates compared with the AETD2 method. To demonstrate the efficiency of our AETD2 method, we compare the simulation time that RK2, ETD2, and AETD2 methods take for a common total run time. We simulate the all-to-all connected network by RK2, ETD2, and AETD2 methods on Visual Studio using an Intel i7 2.6 GHz processor, and the simulation time and numerical accuracy of mean firing rate are given in Table 1. The AETD2 method can achieve over an order of magnitude of speedup compared with the RK2 method while achieving the same high accuracy in terms of the mean firing rate.  Table 1: Simulation of the all-to-all connected network with a total run time T = 10 seconds. The simulation time is measured in seconds. The relative error in the mean firing rate between each method using different time steps and the RK2 method using a very small time step ∆t = 1×10 −6 ms is measured in percentage and the mean firing rate is measured in Hz given inside the parentheses. Asterisks indicate overflow errors.
In addition, we define the efficiency ratio of the AETD2 method over the RK2 method as where T RK2 and T AETD2 indicate the simulation times of the RK2 and AETD2 methods, respectively, for the HH neural network to evolve the run time T . Note that the RK2 and ETD2 methods take almost the same simulation time when using the same small time step as shown in Table 1. Thus, the above efficiency ratio can be approximated by the ratio of the total number of time steps each method requires as where ∆t RK2 and ∆t AETD2 indicate the time steps used in the RK2 and AETD2 methods, respectively. To demonstrate that the above efficiency ratio is independent of the network connectivity and size, we evolve the all-to-all connected network of 80 excitatory and 20 inhibitory neurons and a randomly connected network of 800 excitatory and 200 inhibitory neurons. As shown in Figure 6A, the efficiency ratio approximated by Equation (26) agrees well with the one measured by the ratio of simulation time between the RK2 and AETD2 methods in both two networks. These two networks are further evolved with a variety choice of coupling strength to cover a wide range of different dynamical regimes as shown in Figure 6B. Nevertheless, the efficiency ratio measured by the ratio of simulation time is still consistent with the one approximated by Equation (26) in both networks. Hence, the efficiency ratio of the AETD2 method relies on only the number of evolved time steps and appears to be independent of parameter choice of coupling strength, connectivity structure and network size.

Discussion
We have presented an adaptive second-order ETD method to evolve the pulse-coupled HH neural network. Our AETD2 method can solve the stiff problem in the HH equations when an HH neuron generates an action potential (spike period). It can use a large time step to raise computational efficiency while accurately capturing dynamical properties of HH neurons such as the trace of membrane potentials, spike times of each neuron and the mean firing rate. We point out that our AETD2 method can robustly enlarge time steps and raise computational efficiency over one order of magnitude compared with the RK2 method. This high efficiency seems to be independent of parameter choice of connectivity structure, dynamical regimes or network size. We should point out that our decomposition of HH equations in Equation (14) is an example to apply the ETD method. Other forms of decomposition can be similarly derived under the framework of the ETD method in future. Here, we use the ETD2 scheme derived by approximating the integration in Equation (19) with RK time stepping. Other forms of numerical schemes can also be used to approximate the integration. For example, one can use a liner interpolation to approximate the nonlinear part in Equation (14) to obtain another form of ETD2 scheme. The high efficiency and numerical accuracy can be obtained in our additional numerical experiments.
In this work, the numerical accuracy of our AETD2 method is second-order. In some situations, high accurate traces of membrane potentials may be required, especially the accurate shape of action potentials [18,25]. Therefore, one future work may be the design of the fourth-order ETD method. As illustrated above, due to the discontinuity arising from the pulse-coupled dynamics, an even more careful recalibration needs to be designed to achieve a fourth-order numerical accuracy.
Finally, we point out that our AETD2 method can be easily extended to networks of other HH type neurons [21]. And our AETD2 method can also robustly achieve high numerical accuracy and efficiency. In addition, our method is naturally a parallel algorithm which can be applied to simulations of large-scale neural network dynamics.