Pseudo-Lyapunov exponents and predictability of Hodgkin-Huxley neuronal network dynamics

We present a numerical analysis of the dynamics of all-to-all coupled Hodgkin-Huxley (HH) neuronal networks with Poisson spike inputs. It is important to point out that, since the dynamical vector of the system contains discontinuous variables, we propose a so-called pseudo-Lyapunov exponent adapted from the classical definition using only continuous dynamical variables, and apply it in our numerical investigation. The numerical results of the largest Lyapunov exponent using this new definition are consistent with the dynamical regimes of the network. Three typical dynamical regimes—asynchronous, chaotic and synchronous, are found as the synaptic coupling strength increases from weak to strong. We use the pseudo-Lyapunov exponent and the power spectrum analysis of voltage traces to characterize the types of the network behavior. In the nonchaotic (asynchronous or synchronous) dynamical regimes, i.e., the weak or strong coupling limits, the pseudo-Lyapunov exponent is negative and there is a good numerical convergence of the solution in the trajectory-wise sense by using our numerical methods. Consequently, in these regimes the evolution of neuronal networks is reliable. For the chaotic dynamical regime with an intermediate strong coupling, the pseudo-Lyapunov exponent is positive, and there is no numerical convergence of the solution and only statistical quantifications of the numerical results are reliable. Finally, we present numerical evidence that the value of pseudo-Lyapunov exponent coincides with that of the standard Lyapunov exponent for systems we have been able to examine.

The work was supported by NSF grants DMS-0506396, DMS-0507901 and a grant from the Swartz foundation.

Authors and Affiliations

1. Courant Institute of Mathematical Sciences, New York University, New York, NY, 10012, USA

   Yi Sun, Douglas Zhou, Aaditya V. Rangan & David Cai

2. Department of Mathematics, Shanghai Jiao Tong University, Shanghai, 200240, People’s Republic of China

   David Cai

Corresponding author

Correspondence to Yi Sun.

Action Editor: David Terman

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):

Numerical method for a single neuron

Here we provide details of the numerical method for evolving the dynamics of a single neuron. For simplicity, we use vector X _i to represent all the variables in the solution of the ith neuron:

Given an initial time t ₀ and time step Δt, initial values X _i (t ₀), and spike times \(T^{\rm F}_{i,k}\) and \(T^{\rm S}_{j\neq i,k}\) from the rest of the network, our method computes a numerical solution of all variables X _i (t ₀ + Δt) as well as the intervening spike times \(T^{\rm S}_{i,k}\) (if any occurred) for the ith neuron as follows:

Algorithm 1. (Single neuron scheme)

Step 1: Input: an initial time t₀, a time step Δt, a set of spike times \(T^{\rm F}_{i,k}\) and \(T^{\rm S}_{j\neq i,k}\) and associated strengths \(F_i^{\rm Q}\) and \(S_{i,j}^{\rm Q}\).

Step 2: Consider the time interval [t₀, t₀ + Δt]. Let M denote the total number of feedforward and presynaptic spikes within this interval. Sort these spikes into an increasing list of M spike times \(T^{\rm sorted}_m\) with corresponding spike strengths \(S^{\rm sorted, Q}_m\). In addition, we extend this notation such that \(T^{\rm sorted}_0\) : = t₀, \(T^{\rm sorted}_{M+1}:= t_0 + \Delta t\) and \(S^{\rm sorted, Q}_0 =S^{\rm sorted, Q}_{M+1}:= 0\).

Step 3: For m = 1, ..., M + 1, advance the equations for the HH neuron model and its conductances (Eqs. (1)–(7)) from \(T^{\rm sorted}_{m-1}\) to \(T^{\rm sorted}_m\) using the standard RK4 scheme to obtain \(X_i(T^{\rm sorted}_m)\); Then, update the conductance \(\widetilde{G}_i^{\rm Q}(T^{\rm sorted}_m)\) by adding the appropriate strengths \(S^{\rm sorted, Q}_m\).

Step 4: If the calculated values for \(V_i(T^{\rm sorted}_m)\) are each less than V^th, then we can accept \(X_i(T^{\rm sorted}_{M+1})\) as the solution X _i (t₀ + Δt). We update t₀ ←t₀ + Δt and return to step 2 and continue.

Step 5: Otherwise, let \(V_i(T^{\rm sorted}_m)\) be the first calculated voltage greater than V^th. We know that the ith neuron fired somewhere during the interval \([T^{\rm sorted}_{m-1}, T^{\rm sorted}_m]\).

Step 6: In this case we use a high-order polynomial interpolation to find an approximation of the spike time t^fire in the interval \([T^{\rm sorted}_{m-1}, T^{\rm sorted}_m]\). For example, we can use the numerical values of \(V_i(T^{\rm sorted}_{m-1})\), \(V_i(T^{\rm sorted}_m)\), \(\frac{d}{dt}V_i(T^{\rm sorted}_{m-1})\), \(\frac{d}{dt}V_i(T^{\rm sorted}_m)\) to form a cubic polynomial. We record t^fire as the (k + 1)th postsynaptic spike time \(T^{\rm S}_{i,k+1}\) of the ith neuron. We update t₀ ←t₀ + Δt and return to step 2 and continue.

Numerical method for the HH network

Here we briefly outline an algorithm which accounts for spike-spike interactions, and can accurately evolve recurrent HH networks governed by Eqs. (1)–(7).

Algorithm 2. (Network model scheme)

Step 1: Given initial values initial values X _i (t₀) and feedforward input spike times \(T^{\rm F}_{i,k}\) for all i, we can apply Algorithm 1 and time step Δt to obtain an estimate of neuronal trajectories on the interval [t₀, t₀ + Δt].

Step 2: We use this rough estimate to classify the neurons into two groups , —those that are estimated to fire within [t₀, t₀ + Δt] and those that are not.

Step 3: We sort the approximate spike times of neurons within into a list with corresponding coupling strengths .

Step 4: We use the feedforward spikes \(T^{\rm F}_{i,k}\) as well as the approximate spike times and Algorithm 1 and the same time step Δt to correct the neuronal trajectories of the subnetwork over the interval [t₀, t₀ + Δt], and obtain a more accurate approximation to the spike times . This spike-spike correction procedure is equivalent to stepping through the list and computing the effect of each spike on all future spikes. We step through this correction process until the neuronal trajectories and spike times of neurons in converge.

Step 5: Finally we use the corrected estimates of the spike times as well as the feedforward spikes \(T^{\rm F}_{i,k}\) and Algorithm 1 and time step Δt to evolve the remainder of the neurons from t₀ to t₀ + Δt.

Reprints and permissions