Coarse-grained event tree analysis for quantifying Hodgkin-Huxley neuronal network dynamics

Abstract

We present an event tree analysis of studying the dynamics of the Hodgkin-Huxley (HH) neuronal networks. Our study relies on a coarse-grained projection to event trees and to the event chains that comprise these trees by using a statistical collection of spatial-temporal sequences of relevant physiological observables (such as sequences of spiking multiple neurons). This projection can retain information about network dynamics that covers multiple features, swiftly and robustly. We demonstrate that for even small differences in inputs, some dynamical regimes of HH networks contain sufficiently higher order statistics as reflected in event chains within the event tree analysis. Therefore, this analysis is effective in discriminating small differences in inputs. Moreover, we use event trees to analyze the results computed from an efficient library-based numerical method proposed in our previous work, where a pre-computed high resolution data library of typical neuronal trajectories during the interval of an action potential (spike) allows us to avoid resolving the spikes in detail. In this way, 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 statistical accuracy in terms of average firing rate and power spectra of voltage traces. Our numerical simulation results show that the library method is efficient in the sense that the results generated by using this numerical method with much larger time steps contain sufficiently high order statistical structure of firing events that are similar to the ones obtained using a regular HH solver. We use our event tree analysis to demonstrate these statistical similarities.

Appendices

Appendix A: 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}{lll} && 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}, \ \ \widetilde{S}^{\rm E}=0.05\sim1.0~{\rm mS/cm^2}, \\ && F^{\rm I}=0.01\sim0.05~{\rm mS/cm^2}, \ \ \widetilde{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}$$

Appendix B: 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:

$$ X_i(t)=\big(V_i(t), m_i(t), h_i(t), n_i(t), G_i^{\rm Q}(t), \widetilde{G}_i^{\rm Q}(t)\big). $$

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 other neurons in the network, a preset threshold value V ^th at which an action potential starts, 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 $\widetilde{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.

Appendix C: The library-based numerical method

Here, we briefly outline Algorithm 2, which uses the library to recover the spike. 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 other neurons in the network, our method calculates 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.

In order to build the data library, we choose N _I different values of I ^input equally distributed in its range from 7.5 to 50.0 μA/cm², which can essentially cover all typical values of I ^input of the spiking neurons at the moments when they fire in our network simulations. We also equally distribute N _lib points of m ^input, h ^input, n ^input in their parameter ranges, respectively. These ranges should cover all typical values of m, h, n of a neuron when it fires in the HH network.

Algorithm 2.  (Library algorithm)

Step 0: Pre-compute the data library of V, m, h, n for N _I different values of the constant input current I ^input using very fine time step δ t. For each I ^input , we isolate an action potential form that starts exactly from V ^th to a later time (after T ^ref  ms) where the membrane potential drops down around its minimum. Then we use this action potential form (i.e., the intermediate replica of the membrane potential) with different values of m ^input, h ^input, n ^input as initial values to compute the intermediate replica for m , h , n individually. In particular, we can obtain the reset values m ^re , h ^re , n ^re at the end point where the form of the action potential terminates. For each case of I ^input , there are N _lib data sets of m, h, n, respectively.

Step 1: Input: the library, an initial time t ₀ , a large 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 $\widetilde{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 to 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.

Step 7: We compute the values of $I^{\rm th} = - \sum_{\rm Q} G_i^{\rm Q}(t^{\rm f\/ire})\times$ $(V^{\rm th} - V_{\rm G}^{\rm Q}),$ as well as the gating variables m ^th , h ^th and n ^th at this time. Then, we perform a linear interpolation to find the corresponding reset values of V ^re, m ^re, h ^re, n ^re in the library. Meanwhile, we stop evolving Eqs. ( 1 )–( 5 ) for the next T ^ref  ms, but evolve Eqs. ( 6 ) and ( 7 ) for the conductance terms $G_i^{\rm Q}(t)$ and $\widetilde{G}_i^{\rm Q}(t)$ as usual. We update $t_0\leftarrow \min(t^{\rm f\/ire}+T^{\rm ref}$ , t + Δt) and return to step 2 and continue with the reset values V ^re, m ^re, h ^re, n ^re as the initial values V _i (t ₀), m _i (t ₀), h _i (t ₀) and n _i (t ₀).

Appendix D: More details on discriminability

Here we provide details of several definitions we used in Section 3.2. Although the occurrence count rate n takes only discrete values as shown in the x-axis of Fig. 6, it is more convenient to treat n as a continuous variable for the following discussion. Given a single sample T _obs observation of the network dynamics, our discrimination task is to determine a possible candidate stimulus by using the two distributions P ₁(n) and P ₂(n) of an event chain. A simple procedure is to obtain the occurrence rate n of the event chain in the sample observation and compare it to a threshold number z. If n > z, we report I ₁; otherwise we report I ₂. Figure 6 shows that if we choose z to lie somewhere between the two distributions, say, the location where P ₁(z) = P ₂(z), this procedure can give the correct answer when the two distributions are fairly well separated, but will have difficulty discriminating the stimuli if their distributions merge together. This difficulty is clearly related to the degree to which the two distributions overlap.

First, we show how to define and calculate the hit rate A. The probability that the procedure will make the correct choice (called a hit) under the stimulus I ₁ is the conditional probability that n > z given the stimulus I ₁ was presented, $\alpha:=P(n>z|I_1)=\int_z^\infty P_1(n)dn$. The probability that it will give the answer I ₁ when the stimulus was actually I ₂ (called a false alarm) is similarly $\beta:=P(n>z|I_2)=\int_z^\infty P_2(n)dn$. These two probabilities completely determine the performance of the discrimination procedure because the probabilities for the other two cases (choosing I ₂ when the correct answer is I ₁, choosing I ₂ when the correct answer is I ₂) are $1-\alpha=\int_0^z P_1(n)dn$ and $1-\beta=\int_0^z P_2(n)dn$, respectively. Therefore, by applying this procedure to all different independent single samples in which the two stimuli I ₁ and I ₂ occur with equal probability, the overall probability of making a correct choice with the even chain $\{\sigma^{j_1} \rightarrow \ldots \rightarrow \sigma^{j_m}\}$, i.e., the hit rate is given by

$$\begin{array}{rll} A &=& \frac{1}{2}(\alpha+1-\beta) = \frac{1}{2}\left(\int_z^\infty P_1(n)dn + \int_0^z P_2(n)dn\right) \\ &=& \frac{1}{2}\int_0^\infty \max\big(P_1(n),P_2(n)\big)dn, \end{array}$$

and the false alarm rate is

$$\begin{array}{rll} B &=& 1-A = \frac{1}{2}(1-\alpha+\beta)\\ &=& \frac{1}{2}\int_0^\infty \min\big(P_1(n),P_2(n)\big)dn. \end{array}$$

With A and B in hand, we measure the separation between the two distributions by using the information ratio, which is defined as $I_{I_1,I_2}^{\sigma^{j_1} \rightarrow \ldots \rightarrow \sigma^{j_m}} \equiv A/B$. The information ratio is equal to one if the two distributions totally overlap (i.e., $A=B=\frac{1}{2}$), and it approaches to + ∞ when the two distributions are totally separated (i.e., A = 1, B = 0). Therefore, when we classify the stimulus by using the occurrence count of every event chain within an event tree, the “vote” of each contributing event chain for either stimulus I ₁ or I ₂ is weighted with the factor $\log\left(I_{I_1,I_2}^{\sigma^{j_1} \rightarrow \ldots \rightarrow \sigma^{j_m}}\right)$. Then the weighted votes across the entire event tree are summed up to determine the candidate stimulus underlying the sample observation. This is similar to the concept of entropy in the information theory. For instance, if there are total N _c event chains within an event tree and the vote of each chain $\mathcal C_i$ is denoted by v _i , which takes either value + 1 for stimulus I ₁ or − 1 for I ₂, the summed vote v _s is given by

$$ v_{\rm s} = \sum\limits_{i=1}^{N_{\rm c}} v_i\log\left(I_{I_1,I_2}^{\mathcal C_i}\right). $$

We choose I ₁ if v _s > 0, otherwise we choose I ₂. If the two distributions of an event chain $\mathcal C_i$ totally merge, its vote is void because its weight $\log\left(I_{I_1,I_2}^{\mathcal C_i}\right) =0$. Thus, the final choice is made by the votes with relatively significant weights.

Finally, we apply the above procedure with the information of the entire tree to all different independent single samples. In a similar way as defining the hit rate above, the discriminability of the event tree is defined as the percentage of sample observations that were correctly classified under our voting procedure (choosing I ₁ when the correct answer is I ₁, choosing I ₂ when the correct answer is I ₂). The formula is given by D = N _correct/N _sample where N _sample is the total number of sample observations and N _correct is the number of times of making correct choices. In Figs. 10–13 in Section 4.2, we present the results by showing the discriminability as a function of m _max.