A coarse-graining framework for spiking neuronal networks: from strongly-coupled conductance-based integrate-and-fire neurons to augmented systems of ODEs

Abstract

Homogeneously structured, fluctuation-driven networks of spiking neurons can exhibit a wide variety of dynamical behaviors, ranging from homogeneity to synchrony. We extend our partitioned-ensemble average (PEA) formalism proposed in Zhang et al. (Journal of Computational Neuroscience, 37(1), 81–104, 2014a) to systematically coarse grain the heterogeneous dynamics of strongly coupled, conductance-based integrate-and-fire neuronal networks. The population dynamics models derived here successfully capture the so-called multiple-firing events (MFEs), which emerge naturally in fluctuation-driven networks of strongly coupled neurons. Although these MFEs likely play a crucial role in the generation of the neuronal avalanches observed in vitro and in vivo, the mechanisms underlying these MFEs cannot easily be understood using standard population dynamic models. Using our PEA formalism, we systematically generate a sequence of model reductions, going from Master equations, to Fokker-Planck equations, and finally, to an augmented system of ordinary differential equations. Furthermore, we show that these reductions can faithfully describe the heterogeneous dynamic regimes underlying the generation of MFEs in strongly coupled conductance-based integrate-and-fire neuronal networks.

Appendices

Appendix A: spike resolution

An important feature of the neuronal network we considered in this paper is that the synapses are instantaneous. In other words, the postsynaptic impulse response functions are taken to be δ-functions with no synaptic delay. This δ-function impulse response allows us to determine whether or not one neuron ‘caused’ another neuron to fire, i.e., the instantaneous synapses ensure that cascade-induced firing-events transpire instantaneously. Thus the feedforward inputs will not be confused with the network spikes. To clarify the dynamics that occurs whenever multiple neurons cross threshold simultaneously at time t_spk, we consider our system to be a specific limit of the following fast conductance system

$$\begin{array}{@{}rcl@{}} \frac{d{V}^{{Q}}_{i}(t)}{dt}&=&- g^{E} \left( V^{{Q}}_{i}-{V}_{E}\right) -g^{I} \left( {V}^{{Q}}_{i}-{V}_{{I}}\right),\\ \tau_{{E}} \frac{d g^{{E}}(t)}{dt} &=& - g^{{E}} + \sum\limits_{k}{S}^{{QE}} \delta \left( t-t_{jk}^{{E}}\right), \\ \tau_{I}\frac{dg^{{I}}(t)}{dt} &=& - g^{{I}} +\sum\limits_{k}{S}^{{QI}} \delta \left( t-t_{jk}^{{I}}\right), \end{array} $$

(34)

where the g^Q represent synaptic conductances, and the decay times τ_E ≈ τ_I → 0. Now at the time t_spk, we freeze the macroscopic time t of the original system and create an “infinitesimal-time” system which evolves according to an infinitesimal time τ_Q. Since the MFE is initiated by one excitatory neuron firing before others, at t_spk the state-variables of the original system are given by

$$ {V}_{j}^{{Q}}(t^{-}_{\text{{spk}}}) \text{ for } j = 1,{\cdot}s, {N}_{{Q}} \text{ with } {V}_{1}^{{E}}(t_{\text{{spk}}}^{-}) \geq {V}_{{T}}. $$

(35)

This infinitesimal-time system will use the infinitesimal time τ to describe the effects of the excitatory conductances in Eq. (36) by taking τ_E ≈ τ_I → 0. This is, at this macro-time t_spk, we solve the dynamics of the infinitesimal system

$$\begin{array}{@{}rcl@{}} \frac{d{V}^{{Q}}_{i}(\tau)}{d\tau}&=&- g^{{E}} \left( V^{{Q}}_{i}-{V}_{{E}}\right) -g^{{I}} \left( V^{{Q}}_{i}-{V}_{{I}}\right),\\ \frac{d g^{{E}}(\tau)}{d\tau} &=& - g^{{E}} + \sum\limits_{j\neq i}{S}^{{QE}} \delta \left( \tau-\tau_{j}^{{E}}\right) \\ \frac{dg^{{I}}(\tau)}{d\tau} &=& - g^{{I}} +\sum\limits_{j\neq i}{S}^{{QI}} \delta \left( \tau-\tau_{jk}^{{I}}\right), \end{array} $$

(36)

by fixing τ_E = τ_I = 1, letting the micro-time τ →∞, and using the same conditions as Eq. (35). The spike-times $\tau _{j}^{{Q}}$ record the infinitesimal time at which neurons j fires. By denoting r = g^I/g^E and g(τ) = g^E⋅e^−τ, without ambiguity, we delete the subscript for a while, and rewrite the first equation in (36) as

$$\frac{dv}{d\tau} = -g(\tau) (1+r)(v-{V}_{s}), $$

where V_s = (V_E + rV_I)/(1 + r). Noteing that $\frac {dv}{d \tau } = \frac {dv}{dg}\frac {dg}{d\tau }$, we arrive at

$$\frac{dv}{dg} = (1+r)(v -{V}_{s}) $$

with the solution v(τ) = V_s + g^E(V₀ − V_s)e^1+r(1 − e^−τ). The voltage will evolve to a steady state of

$$\max(v) = \lim\limits_{\tau\rightarrow \infty}v(\tau) = {V}_{s} + g^{{E}}({V}_{0}-{V}_{s}){e}^{1+r}. $$

If max(v) ≥ V_T, the neuron will fire at time

$$\tau_{j}^{{Q}} = \tau^{\prime} - \log \left[1 + \frac{1}{(1+r)g^{{E}}} \log \frac {{V}_{{T}}-{V}_{s}}{{V}_{0}-{V}_{s}}\right], $$

where τ^′ is the time that the last previous neuron fires. By using this analytical formula, we can resolve the trajectory to machine precision during this infinitesimal time until the time $\tau = \min \limits _{j,{ Q}} \tau _{j}^{{Q}}$. Once a single neuron fires, it is clamped at V_R for a refractory period, and never fires again during this infinitesimal time. On the other hand, we can update the voltages of the other neurons that do not fire at this macro-time t_spk. After this infinitesimal time, we go back to the macro-time t to evolve the system. Thus we obtain a well-posed system of ODEs. One can see other similar implementation in Appendix B of Ref. Rangan and Young (2013a) and Appendix 12 of Ref. Zhang and Rangan (2015).

Appendix B: derivation of the Master equation f _{B, single}

In this section, the Master equation f_{B, single} is used to evolve the single-neuron distributions $\boldsymbol {{\rho }}_{{{ single}}}\left (v,t\right ) = ({\rho }^{{E}}_{{{single}}}, {\rho }^{{I}}_{{{single}}})$. Over the small voltage interval [V₁, V₂], the total probability changes as a function of time. The difference from time t₁ to t₁ + Δt is determined by the probability flux at the boundaries V₁ and V₂

$${\int}_{v={V}_{1}}^{v={V}_{2}}\left[{\rho}^{{Q}}\left( v,t\right)\right]_{t=t_{1}}^{t=t_{1}+{\Delta} t}dv=f_{{{B, single}}}^{{Q}}\left( {\rho}^{{Q}},v,t\right) {,} $$

where f_{B, single} = (fB, singleE, fB, singleI) is given by

$$\begin{array}{@{}rcl@{}} && f^{{Q}}_{{B},{{single}}}\left( {\rho}^{{Q}},v,t\right) ={\int}_{t_{1}}^{t_{1}+{\Delta} t}\left[{g}_{{L}}\left( v-{V}_{{L}}\right) {\rho}^{{Q}}\left( v,t\right) \right]_{v={V}_{1}}^{v={V}_{2}}dt \\ &&\quad\quad+\sum\limits_{{R}\in \left\{ E,I,Y\right\}}{\int}_{t_{1}}^{t_{1}+{\Delta} t}{\int}_{T({V}_{1},{S}^{{QR}})}^{T({V}_{2},{S}^{{QR}})}m_{{{ net}}}^{{QR}}\left( t\right) {\cdot} {\rho}^{{Q}}\left( v,t\right) dvdt \\ &&\quad\quad-\sum\limits_{{R}\in \left\{ { E,I,Y}\right\}}{\int}_{t_{1}}^{t_{1}+{\Delta} t}{\int}_{v={V}_{1}}^{v={V}_{2}}m_{{{net}}}^{{QR}}\left( t\right) {\cdot} \mathbf{1}_{\left[-\infty ,T({V}_{T},{S}^{{QR}})\right] }\left( v\right) {\cdot} {\rho}^{{Q}}\left( v,t\right) dvdt \\ &&\quad\quad-{\int}_{t_{1}}^{t_{1}+{\Delta} t}{\int}_{v={V}_{1}}^{v={V}_{2}}m_{{{net}}}^{{QY}}\left( t\right) {\cdot} \mathbf{{\gamma}}^{{Q}}\left( t\right) {\cdot} \mathbf{1}_{\left[T({V}_{{T}},{S}^{{QY}}),{V}_{{T}}\right] }\left( v\right) {\cdot} {\rho}^{{Q}}\left( v,t\right) dvdt \\ &&\quad\quad+\mathbf{1}_{\left[{V}_{1},{V}_{2}\right] }\left( {V}_{{L}}\right) {\int}_{t_{1}}^{t_{1}+{\Delta} t}{\int}_{T({V}_{{T}},{S}^{{EY}})}^{{V}_{{T}}}m_{{{net}}}^{{QY}}\left( t\right) {\cdot} \mathbf{{\gamma}}^{{Q}}\left( t\right) {\cdot} {\rho}^{{Q}}\left( v,t\right) dvdt{,} \end{array} $$

(37)

where the fraction γ^I, and γ^E are used to condition f_{B, single} not to produce MFEs. The motivation to introduce γ^I and γ^E can be found in Zhang et al. (2014a). Generally, γ^I ≡ 1 and γ^E is given by the fraction of time that:

$$\begin{array}{@{}rcl@{}} {\gamma}^{E}\left( t^{+}\right) &= & \left[{\int}_{- \infty}^{T({V}_{{T}},{S}^{{EE}})}{\rho}^{{E}}_{{{single}}}\left( v,t^{-}\right) dv\right]^{{N}_{{ E}}-1}\\&&\times\left[{\int}_{-\infty}^{T({V}_{{T}},{S}^{{IE}})}{\rho}^{{I}}_{{{single}}}\left( v,t^{-}\right) dv\right]^{{N}_{{I}}}{.} \end{array} $$

The network firing rates are given by:

$$\begin{array}{@{}rcl@{}} m_{{net}}^{{QE}}\left( t^{+}\right) &=&\left( {N}_{{Q}}-\delta_{{QE}}\right) m_{{{single}}}^{{E}}\left( t^{-}\right),\\ \text{and } m_{{{ net}}}^{{QI}}\left( t^{+}\right) &=&\left( {N}_{{Q}}-\delta_{{QI}}\right) m_{{{single}}}^{{I}}\left( t^{-}\right) {,} \end{array} $$

(38)

where the single-neuron firing rates msingleQ (t) represents the instantaneous firing rate of a single neuron in an ensemble driven by the various mnetQR (t), and is calculated by the flux of ρ^Q (v, t) over the threshold point V_T at time t in form of

$$\begin{array}{@{}rcl@{}} m_{{{single}}}^{{Q}}\left( t^{-}\right) = {\int}_{t_{1}}^{t_{1}+{\Delta} t}{\int}_{T({V}_{{T}},{S}^{{EY}})}^{{V}_{{T}}}m_{{{net}}}^{{QY}}\left( t\right) {\cdot} \mathbf{{\gamma}}^{{Q}}\left( t\right) {\cdot} {\rho}^{{Q}}\left( v,t\right). \end{array} $$

Appendix C: derivation of the standard Master equation

Here we derive the Master equation from system of conductance-based IF neuronal networks (1). We start by considering what can happen during a single time-interval of length Δt? The coarsest approximation involves choosing from amongst 5 possibilities: (i) the neuron starts out held in refractory and remains at V_L, (ii) the neuron starts out in (V_I, V_T) and decays naturally, (iii) the neuron starts out in (V_I, V_T) and gets a excitatory kick by the external input, or (iv,v) the neuron starts out in (V_I, V_T) and gets a kick of type-Q with synaptic coupling. The simplest possibility is (i), which can be treated by a straightforward computation. In possibility (ii), the pure decay case, v gets mapped to

$$\begin{array}{@{}rcl@{}} { R}_{{L}}^{-1}\left( v,{g}_{{L}}{\Delta} t\right) ={V}_{{L}}+\left( v-{V}_{{L}}\right){\cdot} {e}^{{g}_{{L}}{\Delta} t}={T}_{{L}}\left( v,{g}_{{ L}}{\Delta} t\right). \end{array} $$

Thus, for a small dv, the distribution ρ (v, t)dv within the interval [v, v + dv] is replaced by

$$\begin{array}{@{}rcl@{}} {\rho} \left( {T}_{{L}}\left( v,{g}_{{L}}{\Delta} t\right),t\right)\!{\cdot}\! {T}_{{L}}^{\prime }\left( v,{g}_{{L}}{\Delta} t\right) dv = {\rho}\left( {V}_{{L}} + \left( v - {V}_{{L}}\right){\cdot} {e}^{{g}_{{L}}{\Delta} t},t\right)\!{\cdot}\! {e}^{{g}_{{L}}{\Delta} t}dv. \end{array} $$

(39)

Now let us discuss possibilities (iii,iv,v) which correspond to getting a kick from the external input and other excitatory or inhibitory neurons in the network. For example, given a single kick of the external Poisson input, the distribution concentrated in the interval [V₁, V₂]gets mapped to the interval [R (V₁, S^QY) , R (V₂, S^QY)], and is mapped into by the interval [T (V₁, S^QY) , T (V₂, S^QY)]. Thus, if V₁ = v, and V₂ = v + dv for dv small, then ρ (v)dv is replaced by ρ (T (v, S^QY) , t)⋅T^′(v, S^QY) dv, implying that, for each voltage v, the distribution ρ (v) is replaced by the distribution ${\rho } \left ({T}\left (v,{S}^{{QY}}\right ),t\right ) {\cdot } {T}^{\prime }\left (v,{S}^{{QY}}\right ) ={\rho } \left ({V}_{{E}}+\left ({V}-{V}_{{E}}\right ){\cdot } {e}^{{S}^{{QY}}},t\right ) {\cdot } {e}^{{S}^{{QY}}}$. Similarly, we have corresponding maps for excitatory and inhibitory network spikes.

Thus, combining these 5 possibilities, ρ (v, t) evolves according to

$$\begin{array}{@{}rcl@{}} {\rho}^{{Q}} \left( v,t+{\Delta} t\right) &=&\left[1-({m}^{{QY}}_{{net}} + {N}_{{E}}{m}^{{E}} + {N}_{{I}}{m}^{I}){\Delta} t\right] {\cdot} {\rho}^{{Q}} \left( {{V}}_{{L}}+\left( {V}-{V}_{{L}}\right) {\cdot} {e}^{{g}_{{L}}{\Delta} t},t\right) {\cdot} {e}^{{g}_{{L}}{\Delta} t} \\ &&+{m}^{{QY}}_{{net}}(t){\Delta} t{\cdot} {\rho}^{{Q}} \left( {V}_{{E}}+\left( v-{V}_{{E}}\right){\cdot} {e}^{{S}^{{QY}}},t\right) {\cdot} {e}^{{S}^{{QY}}} \\ &&+{N}_{{E}}{m}^{{E}}(t){\Delta} t{\cdot} {\rho}^{{Q}} \left( {V}_{{E}}+\left( v-{V}_{{E}}\right){\cdot} {e}^{{S}^{{QE}}},t\right) {\cdot} {e}^{{S}^{{QE}}} \\ &&+{N}_{{I}}{m}^{{I}}(t){\Delta} t{\cdot} {\rho}^{{Q}} \left( {V}_{{I}}+\left( v-{V}_{{I}}\right) {\cdot} {e}^{{S}^{{QI}}},t\right) {\cdot} {e}^{{S}^{{QI}}}+\mathcal{{S}}^{{Q}}\left( t\right) {\cdot} \delta\left( {V}-{V}_{{R}}\right) \end{array} $$

(40)

with Q ∈{E, I}. Taylor expanding,

$$\begin{array}{@{}rcl@{}} {\rho}^{{Q}} \left( {V}_{{L}}+\left( {V}-{V}_{{L}}\right) {\cdot} {e}^{{g}_{{L}}{\Delta} t},t\right){\cdot} {e}^{{g}_{{L}}{\Delta} t,t} = {\rho}^{{Q}}(v,t) + {g}_{{L}}{\Delta} t {\partial}_{v}[(v-{V}_{{L}}){\rho}^{{Q}}(v,t)] + o({\Delta} t). \end{array} $$

Substituting the above into Eq. (40) and ignoring the terms of order o(Δt) and higher, we arrive at

$$\begin{array}{@{}rcl@{}} {\rho}^{{Q}} \left( v,t+{\Delta} t\right) - {\rho}^{{Q}}(v,t) &=&{g}_{L}{\Delta} t {\partial}_{v} \left[(v-{V}_{{L}}){\rho}^{{Q}}(v,t)\right] +{m}^{{QY}}_{{net}}(t){\Delta} t {\cdot} \left[{\rho}^{{Q }}\left( T(v,{S}^{{QY}}),t\right) {\cdot} {T}^{\prime}(v,{S}^{{QY}}) - {\rho}^{{ Q}}(v,t)\right] \\ &&+{N}_{{E}}{m}^{{E}}(t){\Delta} t{\cdot} \left[{\rho}^{{Q }}\left( T(v,{S}^{{QE}}),t\right){\cdot} T^{\prime}(v,{S}^{{QE}}) - {\rho}^{{Q}}(v,t)\right] \\ &&+{N}_{{I}}{m}^{{I}}(t){\Delta} t {\cdot} \left[{\rho}^{{Q}} \left( T(v,{S}^{{QI}}),t\right) {\cdot} {T}^{\prime}(v,{S}^{{QI}}) - {\rho}^{{Q}}(v,t)\right] +\mathcal{{S}}^{{Q}}\left( t\right) {\cdot} \delta\left( {V}-{V}_{{R}}\right), \end{array} $$

(41)

where $\mathcal {S}^{Q}(t)$ is the flux across the threshold and is given by

$$\begin{array}{@{}rcl@{}} \mathcal{S}^{Q}\left( t\right) &=&{m}^{QY}_{net}(t) {\int}_{T({V}_{T},{S}^{QY})}^{{V}_{T}}{\rho}^{Q}(v,t)dv +{N}_{E}{m}^{E}(t){\int}_{T({V}_{T},{S}^{QE})}^{{V}_{T}} {\rho}^{Q}(v,t)dv. \end{array} $$

(42)

Let Δt → 0, we have the Master equation in the following form

$$\begin{array}{@{}rcl@{}} \frac{d{\rho}^{Q} \left( v,t\right)}{dt} &=&{g}_{L}{\partial}_{v} \left[(v-{V}_{L}){\rho}^{Q}(v,t)\right] +{m}^{QY}_{net}(t) {\cdot} \left[{\rho}^{Q} \left( T(v,{S}^{QY}),t\right) {\cdot} T^{\prime}(v,{S}^{QY}) - {\rho}^{Q}(v,t)\right] \\ &&+{N}_{E}{m}^{E}(t){\cdot} \left[{\rho}^{Q} \left( T(v,{S}^{QE}),t\right){\cdot} T^{\prime}(v,{S}^{QE}) - {\rho}^{Q}(v,t)\right] \\ &&+{N}_{I}{m}^{I}(t) {\cdot} \left[{\rho}^{Q} \left( T(v,{S}^{QI}),t\right) {\cdot} T^{\prime}(v,{S}^{QI}) - {\rho}^{Q}(v,t)\right]+\mathcal{S}^{Q}\left( t\right) {\cdot} \delta\left( V-{V}_{R}\right). \end{array} $$

(43)

For a similar derivation of Eq. (43), see Nykamp (2000) and Cai et al. (2006).

Appendix D: derivation of Fokker-Planck equation

Similar to the reductions performed in Rangan and Cai (2006) and Helias et al. (2010), we take a single neuron driven only by the external Poisson input for example. More complicated case can be computed similarly. Let $s ={e}^{{S}^{{QY}}} - 1$, we can rewrite

$$T\left( v,{S}^{{QY}}\right) = v+s(v-{V}_{{E}}), \text{and} T^{\prime}(v,{S}^{{QY}})= 1+s. $$

Ignoring the network synaptic connection for now, and considering only the external Poisson input, thus the Master Eq. (8) is given by

$$\begin{array}{@{}rcl@{}} {\partial}_{t} {\rho}^{{Q}}(v,t) & = & {g}_{{L}}{\partial}_{v} \left[(v - {{V}}_{{L}}){\rho}^{{Q}}(v,t)\right] + {m}^{{QY}}_{{net}}(t) \\&&{\cdot} \left[{\rho}^{{Q}} \left( v + (v - {V}_{{E}})s,t\right)(1 + s) - {\rho}^{{Q}}(v,t)\right]. \end{array} $$

(44)

Taylor expanding and keeping the first order in s, we get

$$ {\rho}^{{Q}} (v+(v-{{V}}_{{E}})s)\approx {\rho}^{{Q}}(v) + (v-{V}_{{E}})s{\partial}_{v}{\rho}^{{Q}} + \frac{(v-{V}_{{ E}})^{2}{S}^{2}}{2}{\partial}_{vv}{\rho}^{{Q}}, $$

(45)

thus we can approximate (44) by

$$\begin{array}{@{}rcl@{}} {\partial}_{t} {\rho}^{{Q}}(v,t) &=& {g}_{{L}}{\partial}_{v} \left[(v-{{V}}_{{L}}){\rho}^{{Q}}(v,t)\right] \\&&+{m}^{{QY}}_{{net}}(t){\cdot} \left[(v-{V}_{{E}})s(1+s){\partial}_{v}{\rho}^{{Q}}+ \frac{(v-{V}_{{E}})^{2}{S}^{2}(1+s)}{2}{\partial}_{vv}{\rho}^{{Q}} + s{\rho}^{{Q}}\right]. \end{array} $$

(46)

Noting the relationships

$$\begin{array}{@{}rcl@{}} (v-{V}_{E})s{\partial}_{v}{\rho}^{{Q}} &= &{\partial}_{v} \left[s(v-{V}_{{E}}){\rho}^{{Q}}\right] -s{\rho}^{{Q}}, \\ \frac{(v-{V}_{{E}})^{2}{S}^{2}}{2}{\partial}_{vv}{\rho}^{{Q}} &=&{\partial}_{v}\left[\frac{(v-{V}_{{E}})^{2}{S}^{2}}{2}{\partial}_{v}{\rho}^{{Q}}\right] \\&&- {\partial}_{v}\left[{S}^{2}(v-{{V}}_{{E}}){\rho}^{{Q}}\right] + {S}^{2}{\rho}^{{Q}}, \end{array} $$

we have

$$\begin{array}{@{}rcl@{}} &&(v-{V}_{{E}})s{\partial}_{v}{\rho}^{{Q}} + \frac{(v-{V}_{{E}})^{2}{S}^{2}}{2}{\partial}_{vv}{\rho}^{{Q}} \\&=&{\partial}_{v}\left[(s-{S}^{2})(v-{V}_{{E}}){\rho}^{{Q}} + \frac{(v-{V}_{{E}})^{2}{S}^{2}}{2}{\partial}_{v}{\rho}^{{Q}}\right] - (s-{S}^{2}){\rho}^{{Q}}. \end{array} $$

(47)

Substituting (47) into (46), we arrive at the Fokker-Planck equation

$$\begin{array}{@{}rcl@{}} {\partial}_{t} {\rho}^{{Q}}(v,t) &= &{\partial}_{v} \{\left[{g}_{{L}}(v-{V}_{{L}}) + {m}^{{QY}}_{{net}}(t) (s-{S}^{3})(v-{V}_{{E}})\right]{\rho}^{{Q}}(v,t) \} \\ && +{\partial}_{v}\{\frac{{m}^{{QY}}_{{net}}(t)(1+s){S}^{2}(v-{V}_{{E}})^{2}}{2}{\partial}_{v}{\rho}^{{Q}} \} + {m}^{{QY}}_{{net}}(t){S}^{3}{\rho}^{{Q}}. \end{array} $$

(48)

We remark that the last term ${m}^{{QY}}_{net}(t){S}^{3}{\rho }^{{Q}}$ in Eq. (48) has a small influence on the distribution when the parameter s is small. Hence, if we only keep the second order of s in the Taylor’s expansion of Eq. (45), we can ignore this term in most practical analysis. Similarly, for the full IF network, we let $s = {e}^{{S}^{{QY}}} - 1, {s}_{1} = {e}^{{S}^{{QE}}} - 1, {s}_{2} = {e}^{{S}^{{QI}}} - 1$, and can approximate the Master Eq. (8) by a FP-type Eq. (9).