Dynamics of Networks of Leaky-Integrate-and-Fire Neurons

Abstract

The dynamics of pulse-coupled leaky-integrate-and-fire neurons is discussed in networks with arbitrary structure and in the presence of delayed interactions. The evolution equations are formally recasted as an event-driven map in a general context where the pulses are assumed to have a finite width. The final structure of the mathematical model is simple enough to allow for an easy implementation of standard nonlinear dynamics tools. We also discuss the properties of the transient dynamics in the presence of quenched disorder (and δ-like pulses). We find that the length of the transient depends strongly on the number N of neurons. It can be as long as 10⁶–10⁷ inter-spike intervals for relatively small networks, but it decreases upon increasing N because of the presence of stable clustered states. Finally, we discuss the same problem in the presence of randomly fluctuating synaptic connections (annealed disorder). The stationary state turns out to be strongly affected by finite-size corrections, to the extent that the number of clusters depends on the network size even for N≈20,000.

Appendix: Rescaling the Equations of Motion

Models of LIF neurons are often introduced by referring to different normalizations. In order to facilitate the comparison of the results obtained by different groups, in the following we illustrate the rescaling needed to express the equations in the adimensional setting adopted in (11.1).

11.1.1 A.1 Setup A

Zillmer et al. [44] analyze a network of N LIF neurons by referring to the equations

$$\beta\dot{V}_{i}=I_{\mathrm{ext}}-V_{i}+\beta\frac{J}{K}\sum_{l=1}^{N}\mu_{i,l}\sum_{m}F\bigl(t-t_{l}^{(m)}-D\bigr),\quad i=1,\dots,N$$

(11.26)

(some variable names have been changed to avoid the confusion arising from the overlap between symbols which denote different quantities), where K is the number of incoming links to each neuron (f=K/N is the dilution of the network); β is the membrane time constant of the neuron; I _ext is an external current and D is the synaptic delay. Moreover, the model definition includes the reset potential V _r and the firing threshold V _t . The parameter J with J<0 (resp., J>0) for inhibitory (resp., excitatory) coupling represents the coupling strength, while the topology of the network is defined by the connectivity matrix μ _i,l . The distribution P(μ) is chosen to be dichotomic, i.e. P(μ)=(1−f)δ(μ)+f δ(μ−1), which implies that the average is 〈μ〉=f, while the variance is σ ²(μ)=f(1−f). Finally, the function F(t) (that becomes a Dirac’s δ-function in the case of zero-width pulses) describes the shape of the single pulse, while $t_{l}^{(m)}$ represents the mth firing time of the lth neuron.

The above equations transform into (11.1), once the following changes of variables are introduced,

• t→t/β,

• D→τ=D/β,

• V→v=(V−V _r )/(V _t −V _r ),

• I _ext→a=(I _ext−V _r )/(V _t −V _r ),

• J→g=J/(V _t −V _r ),

• μ _i,j →S _i,j =μ _i,j /〈μ〉=μ _i,j /f,

where the new distribution P(S) of connections strengths writes $P(S)= (1-\nobreak f)\delta(S)+f\delta(S-1/f)$, so that its average is equal to one as required, while the variance is σ ²(S)=(1−f)/f. The parameter values corresponding to this setup are summarized in Table 11.1.

11.1.2 A.2 Setup B

In Jahnke et al. [20, 21] the model is defined as

$$\dot{V}_{i}=I_{\mathrm{ext}}-\gamma V_{i}+\sum_{j=1}^{N}\epsilon_{i,j}\sum_{m}F\bigl(t-t_{j}^{(m)}-D\bigr),\quad i=1,\dots,N,$$

(11.27)

where V _r =0, V _t =1 and $\sum_{j=1}^{N}\epsilon_{i,j}=J_{T}$. The coupling strengths ε _i,j are randomly chosen according to the distribution

$$P(\epsilon)=(1-f)\delta(\epsilon)+\frac{f^3N}{2J_T}H(\epsilon)H\biggl(\frac{2J_T}{Nf^2}-\epsilon\biggr)$$

(11.28)

where H(x) is the Heaviside step function.

By performing the transformations,

• t→t γ,

• D→τ=D γ,

• I _ext→a=I _ext/γ,

• J _T →g=J _T ,

• ε _i,j →S _i,j =ε _i,j /〈ε〉=ε _i,j Nf/J _T ,

the model (11.27) can be rewritten in the form (11.1). The probability of connections strengths becomes

$$P(S)=(1-f)\delta(S)+\frac{f^2}{2}H(S)H(2/f-S)$$

so that 〈S〉=1, and σ ²(S)=4/(3f)−1. The parameter values corresponding to this setup are summarized in Table 11.1.