Distribution of correlated spiking events in a population-based approach for Integrate-and-Fire networks

Randomly connected populations of spiking neurons display a rich variety of dynamics. However, much of the current modeling and theoretical work has focused on two dynamical extremes: on one hand homogeneous dynamics characterized by weak correlations between neurons, and on the other hand total synchrony characterized by large populations firing in unison. In this paper we address the conceptual issue of how to mathematically characterize the partially synchronous “multiple firing events” (MFEs) which manifest in between these two dynamical extremes. We further develop a geometric method for obtaining the distribution of magnitudes of these MFEs by recasting the cascading firing event process as a first-passage time problem, and deriving an analytical approximation of the first passage time density valid for large neuron populations. Thus, we establish a direct link between the voltage distributions of excitatory and inhibitory neurons and the number of neurons firing in an MFE that can be easily integrated into population–based computational methods, thereby bridging the gap between homogeneous firing regimes and total synchrony.

1. Donsker’s theorem states that the fluctuations of an empirical CDF about its theoretical CDF converge to Gaussian random variables with zero mean and certain variance. The sequence of independent Gaussian random variables can be formulated in terms of a standard Brownian bridge, a continuous-time stochastic process on the unit interval, conditioned to begin and end at zero.

The authors would like to thank David Cai for useful discussions. J. Z. is partially supported by NSF grant DMS-1009575, K. N. is supported by the Courant Institute. D. Z. is supported by Shanghai Pujiang Program (Grant No. 10PJ1406300), NSFC (Grant No. 11101275 and No. 91230202), as well as New York University Abu Dhabi Research Grant G1301. A. R. is supported by NSF Grant DMS-0914827.

Authors and Affiliations

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

   Jiwei Zhang, Katherine Newhall & Aaditya Rangan

2. Department of Mathematics, MOE-LSC and Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai, 200240, China

   Douglas Zhou

Corresponding author

Correspondence to Katherine Newhall.

Action Editor: Gaute T. Einevoll

Appendixes

1.1 Appendix A: Firing events in I&F networks

In order to accurately simulate the dynamics of Eq. (1a), we need to resolve the neurons’ firing sequence correctly despite the fact that all neurons fire during the same instance of time in a single MFE. The algorithm used must produce an MFE magnitude consistent with the cascade condition (6) in Sec. 3.1. We present one such algorithm here in which the voltages of the remaining neurons are updated based on the neuron with the highest voltage first. We also emphasize that neurons having fired previously in an MFE do not receive input from neurons firing after it within the same MFE, because the neurons have a short refractory period, thus neurons can never fire more than once in an MFE.

To resolve the firing sequence of an MFE, we start from the sets \(\left \{v_{j}\right \}_{j=1}^{N_{E}}\) and \(\left \{w_{j}\right \}_{j=1}^{N_{I}}\) of the excitatory and inhibitory neuronal voltages at the time one excitatory voltage is above threshold, indicating that this neuron is about to fire. Then, the following algorithm is used:

1. 1.

   Find the k _E excitatory and k _I inhibitory voltages within the sets v _j and w _j that are above the threshold, V _T .

2. 2.

   If both k _E and k _I are zero, stop; the MFE has ended. Otherwise, find the neuron with largest voltage, V _max, out of the k _E excitatory and k _I inhibitory voltages above the threshold.

3. 3.

   The neuron with voltage V _max fires; it is reset to V _R . If it is type E (I), it is removed from the list v _j {w _j } and the remaining voltages of neurons that have not fired are updated by adding S ^EEto (subtracting S ^EIfrom) the voltages in {v _j } and adding S ^IEto (subtracting S ^IIfrom) the voltages in w _j .

4. 4.

   If either set v _j or w _j is non-empty, return to Step 1. Otherwise, stop.

The MFE magnitudes for the excitatory and inhibitory populations are obtained by subtracting the number of neurons that did not fire in the MFE from the total numbers of neurons in the two populations,

where # denotes the number of elements in the set. Note that the sets in Eq. (36) only include the voltages of non-firing neurons, as detailed in Step 3 of the algorithm.

1.2 Appendix B: Derivation of the geometric method

Starting from condition (6) in the main text, we describe here how to obtain a single condition, and therefore the MFE magnitude described by the single intersection point of G(v) defined in Eq. (14) and the line l(v) in Eq. (5). From the last points v < V _T and w < V _T where condition (6) is satisfied,

we have

Subtracting Eq. (40) from Eq. (39) we have

which can be written as

where δ = S ^II S ^EE / S ^IE − S ^EI. The first two terms on the RHS suggest the transformation of the inhibitory voltages defined in Eq. (7) in the main text,

This, together with the fact that one can show the empirical CDF

allows Eq. (41) to be written as

Next, we must consider the two cases of δ > 0 and δ < 0 separately.

Case 1

If δ > 0, we further transform \(\hat {w}\) by

appearing as Eq. (9) in the main text. Then, Eq. (44) is simply

and we substitute this into Eq. (39), together with Eq. (43) to obtain

What we want to do next is replace the function of \(\hat {w}\) by a function of \(\bar {w}\). Suppose we can invert the transformation \(\bar {w}= \hat {w} - \delta N_{I} \left (1-\hat {F}_{I}(\hat {w})\right )\) to obtain

(Note that \(h^{-1}(x) = x - \delta N_{I}(1-\hat {F}_{I}(x))\), we have \(\left [h^{-1}(x)\right ]' = 1+\hat {F}'_{I}(x)>0\). Therefore, the transformation h is a monotonic increasing map and has only one root.) Then we have

Changing the integration to the variable y defined by y = h ^{− 1}(z), we have

Taylor expanding h(y) and only keeping the linear term, we have \(h(y)-h(\bar {w}) = C_{0}(y-\bar {w})\) and h ^″(y) = C ₀, resulting in

Substituting Eq. (49) into Eq. (47), we have

corresponding to Eq. (12) in the main text, where we defined \(\bar {v}=v\) so that \(\bar {F}_{E}(v) = F_{E}(v)\). Solving this for \(\bar {w}\) is equivalent to finding the intersection point between the line \(l(v) = 1+\frac {1}{N_ES^{EE}}(v-V_T)\) and the function

Case 2

If δ < 0, then \(\bar {w}\) defined in Eq. (45) may be larger than V _T , which we would like to avoid. Rather, we solve Eq. (44) for \(\hat {w}\),

and transform v to (Eq. (10) in the main text)

so that Eq. (51) becomes

Using Eq. (53) to rewrite Eq. (52) in terms of \(\bar {v}\) and solving for v gives

Using Eq. (54) and that \(F_{I}(w) = \hat {F}_{I}(\hat {w})\), Eq. (40) becomes

corresponding to Eq. (12) in the main text, where we defined \(\bar {w}=\hat {w}\), so that \(\bar {F}_{I}(v) = \hat {F}_{I}(v)\).

Now, what we want to do is replace the function of v by a function of \(\bar {v}\). If we denote Eq. (54) by

then change the integration variable in

to y using z = g(y) and \(v_{j} = g(\bar {v}_j)\), we have that

following the same argument used to derive Eq. (49). Now, Eq. (55) becomes

Solving this for \(\bar {v}\) is equivalent to finding the intersection point between the line \(l(v) = 1 + \frac {1}{N_ES^{EE}} (v-V_T)\) and the function

Combining these two cases, we arrive at the MFE magnitude defined via the intersection of the function G(v) defined in Eq. (14) in the main text, and the line l(v) defined in Eq. (5) in the main text.

1.3 Appendix C: Analytical formula of first passage time

In this appendix, we explain how to reformulate the passage of a generic 2-dimensional anisotropic Brownian motion across a moving boundary to the viewpoint of one-dimensional Brownian motion as used in Durbin’s papers (Durbin 1985; Durbin and Williams 1992).

We define x = (x, y) to be two-dimensional anisotropic brownian motion with linear drift, given by the solution to

with initial starting point x(0) = 0 for simplicity. If β _x ≠ β _y , then we define the isotropic version of this process via \(\mathbf {z}=( \hat {x},\hat {y})\), with \(\hat {x}=x/\beta _{x}\), \(\hat {y}=y/\beta _{y}\), or z = β ^{− 1} x, with

We will also define

We consider the first crossing density of x to a linear boundary \( \mathcal {A}\) described by the normal constraint

where \(\bar {\eta }=\left [\bar {\eta }_{x},\bar {\eta }_{y},\bar {\eta }_{t}\right ]^{\intercal }\) is the unit normal to \(\mathcal {A}\) in (x, y, t) space and k is the constant ‘offset’ of the linear plane (i.e., the distance between \(\mathcal {A}\) and the origin (0, 0, 0)). We also define \(\eta _{\perp }=\left [ \bar {\eta }_{x},\bar {\eta }_{y} \right ]^{\intercal }\) (not a unit vector) to be the projection of \(\bar {\eta }\) onto the x, y plane. So we know the process (xt), y(t) crosses \(\mathcal {A}\) when

or where

or when

That is, x crosses \(\mathcal {A}\) when z crosses \(\bar {\beta } ^{-1}\mathcal {A}\) (i.e., the boundary perpendicular to \(\bar {\beta }\bar {\eta } \)). Now z is isotropic brownian motion (plus a linear drift), so z can be described in terms of z _⊥ and z _∥, the directions parallel and perpendicular to η _⊥. Each of the processes z _⊥ and z _∥ are given by independent brownian motion (plus a linear drift). So the probability of z hitting the rescaled boundary \(\bar {\beta }^{-1}\mathcal {A}\) at time t is just the probability of z _⊥ reaching the moving (scalar) boundary \( \left [ \bar {\beta }^{-1}\mathcal {A}\right ]_{\perp }=\left ( k-\bar {\eta } _{t}t\right ) /\left \vert \beta \eta _{\perp }\right \vert \) at time t. That is to say

According to Durbin (1985), this probability is simply

Now if z crosses at time t at a particular point \(\bar {\beta } ^{-1}a( t) \), then two things must happen. First, z _⊥ must equal \(\left [ \bar {\beta }^{-1}a( t) \right ]_{\perp }\) which is also \(\left ( k-\bar {\eta }_{t}t\right ) /\left \vert \beta \eta _{\perp }\right \vert \) at time t. Second, z _∥ must equal \(\left [ \bar {\beta } ^{-1}a( t) \right ]_{\parallel }\) at time t, where \(\left [ \bar {\beta }^{-1}a( t) \right ]\) is the component of a(t) perpendicular to η _⊥. Thus, the probability that the process z crosses at time t at a particular point \(\bar {\beta }^{-1}a( t) \) on \(\bar {\beta }^{-1} \mathcal {A}\) is given by

Now, as we discussed above,

For the second term, we have

and together, since z _⊥ and z _∥ are independent,

As we discussed above, if z crosses at time t at a point \(\bar {\beta } ^{-1}a( t) \), then x crosses at time t at at. However, the density \( p\left ( t\text {, }\mathbf {z}\text { crosses at }\bar {\beta } ^{-1}a( t) \right )\) is not equivalent to the density p (t, x crosses at a(t)). Rather, we have

where \(\bar {\eta }_{\parallel }=\left [ -\bar {\eta }_{y},\bar {\eta }_{x},0 \right ]^{\intercal }\), and \(\eta _{\parallel }=\left [ -\bar {\eta }_{y},\bar {\eta }_{x}\right ]^{\intercal }\) (not a unit vector) lies along \(\mathcal {A}\) and lies in the x-y plane, and is perpendicular to η _⊥. Comparing the lengths of these line segments, we see that length \(\left | \bar {\beta }^{-1}a( t) ,\bar {\beta } ^{-1}( a\left ( t) +\delta \cdot \bar {\eta }_{\parallel }\right ) \right | =\delta \cdot \left \vert \beta ^{-1}\eta _{\parallel }\right \vert \) and length \(\left | a( t) ,a( t) +\delta \cdot \bar {\eta }_{\parallel }\right |=\delta \cdot \left \vert \eta _{\parallel }\right \vert \text {.} \) Using this relationship we can conclude the density

This in turn implies that the probability that the process x crosses \(\mathcal {A}\) at time t at a particular point a(t) on \( \mathcal {A}\) is given by the density

Finally, since

where β is the determinant of β. The density

Using the language described in Section 3.3, the above can be written as

leading to

where f(t, a | 0, x0) is a (unconditioned) density for the process x at t, a, see details in Appendix B.

We can now fill in the geometry for the specific problem at hand. The points (x, y) on the surface \(\mathcal {A}\) at time t was given in Eq. (23) and comes from the constraint of \(\bar {G}(t) = \bar {l}(t)\). The vector

is normal to the surface in 3D space at the point (x, y, t). The 2-component normal to \(\mathcal {A}\), restricted to the x-y plane, is

and the 2-component vector perpendicular to n _⊥ is

where | η | is the length of the vector η. One can also think of \(\mathcal {A}\) as a line parallel to − γx + αy in the x-y plane, with a time dependent offset of \(k/N_{E}+\gamma \bar {f}_{E}( t) -\alpha \bar {f}_{I}( t) -\frac {1}{ N_{E}S^{EE}}t\). So, around any point (a(t), t) = a _x , a _y , t on \(\mathcal {A}\) , the point

is on the linearized boundary tangent to \(\mathcal {A}\) at the point a(t). The factor

can be thought of as the speed at which the boundary propagates in the normal direction, as s (and hence t − s) changes and the unit vector in the x-y plane \( \mathbf {u} = \frac {[-\gamma ,\alpha ] }{\sqrt {\gamma ^{2}+\alpha ^{2}}} \) is perpendicular to \(\widehat {\mathcal {A}}\) and therefore parallel to n _⊥. Therefore,

We now change variables to the n _⊥ and n _∥directions by first changing to variables in which the stochastic processes (ϕ _E (t), ϕ _I (t)) obeys isotropic Brownian motion in 2D. This is done with the matrix

and requires the change of variables term given by

where \(\zeta =\sqrt {1+1/\nu _{n}^{2}}\).

Using the definition \( l( s,\mathbf {x}|t,\mathbf {a}) =\left ( \widehat {\mathcal {A}}_{a\left ( t\right ) }|_{s}-\mathbf {x}\right ) \cdot \mathbf {u}\; \zeta \; , \) we can calculate the density of crossing at time t and position a on \(\mathcal {A}\) via

In the above expression the term f(t, a) refers to the unconditioned density of the process (ϕ _E , ϕ _I , t) at time t and position a (presumed to lie on \(\mathcal {A}\)). The first term in the sum is straightforward, the second term is a little trickier. Note that the conditional first passage density at time r to point b given that the process also crosses the point a at time t can be written as q(r, b)f(r, b | t, a) = qr, bf(r, b, t, a) / f(t, a), where q(r, b) is defined via p(r, b) = q(r, b)f(r, b). Following Durbin’s method, we achieve the analytical formula (32) by taking the first two terms of Eq. (59) in the main text.

1.4 Appendix D: Probability density of stochastic process

In this appendix, we derive the distributions f(t, a) in Eq. (33) and f(r, b | t, a) in Eq. (34) in the main text. Recall that \(\bar {f}_{E}(t)\) and \(\bar {f}_{I}(t)\) are the smooth CDFs for the transformed voltage variables, defined in Eq. (18). These define the stochastic processes ϕ _E (t) and ϕ _I (t) in Eq. (20).

The distribution f(t, a) refers to the probability of observing the unconstrained process (ϕ _E , ϕ _I ) at at at time t. In other words, f(t, a) is the probability that

at time t, without assuming anything about how ϕ crosses \(\mathcal {A}.\) This probability is given by:

where each component of ϕ evolves independently. Now of course the probability that ϕ _Q (t) = a _Q (t) at the time t is the same as the probability that \(\bar {\phi }_{Q}\left ( \tau \right ) =a_{Q}( t) \) at time τ = f _Q (t). This in turn is the same as \(\sqrt {N_{Q}}\) times the probability \(P\left ( B_{\tau }=\sqrt {N_{Q}}a_{Q}\right ) \), the probability that a standard brownian bridge (connecting (0, 0) to (1, 0)) reaches the value \(\sqrt {N_{Q}}a_{Q}( t) \) at time τ = f _Q t. According to Durbin and Williams (1992), this probability is given by

therefore f(t, a) is precisely Eq. (33) in the main text.

In a similar vein the distribution f(r, b | t, a) is the conditional distribution of the process ϕ at r, b, given that the point t, a is reached. Using logic similar to the above, this is given by

Note that B _τ conditioned to hit \(\sqrt {N_{Q}}a_{Q}\) at time f _Q (t) is just a brownian bridge connecting 0, 0 to \(( f_{Q}\left ( t) ,\sqrt {N_{Q}}a_{Q}\right ) \). One can think of this new brownian bridge as a rescaled version of a standard brownian bridge:

Then

and formula (34) for f(r, b | t, a) in the main text follows. In the computation of MFE magnitude, we must simply account for not selecting the k neurons that initiate the MFE, so we replace N _E by N _E − k.

Reprints and permissions