13 Continuity Equation and the Fokker-Planck Approach

13.4 Networks of leaky integrate-and-fire neurons

In the previous section, the formalism of membrane potential densities has been applied to a single population driven by spikes arriving at some rate νk\nu_{k} at the synapse of type kk, where k=1,,Kk=1,\dots,K. In a population that is coupled to itself, the spikes that drive a given neuron are, at least partially, generated within the same population. For a homogeneous population with self-coupling the feedback is then proportional to the population activity A(t)A(t); cf. Chapter 12.

We now formulate the interaction between several coupled populations using the Fokker-Planck equation for the membrane potential density (Section 13.4.1) and apply it subsequently to the special cases of a population of excitatory integrate-and-fire neurons interacting with a population of inhibitory ones (Section 13.4.2).

13.4.1 Multiple populations

We consider multiple populations k=1,,Kk=1,\dots,K. The population with index kk contains NkN_{k} neurons and its activity is denoted by AkA_{k}. We recall that NkAk(t)ΔtN_{k}\,A_{k}(t)\,\Delta t is the total number of spikes emitted by population kk in a short interval Δt\Delta t.

Let us suppose that population kk sends its spikes to another population nn. If each neuron in population nn receives input from all neurons in kk the total spike arrival rate from population kk is therefore νk(t)=NkAk(t)\nu_{k}(t)=N_{k}\,A_{k}(t) (‘full connectivity’). If each neuron in population nn receives connections only from a subset of CnkC_{nk} randomly chosen neurons of population kk, then the total spike arrival rate is νk(t)=CnkAk(t)\nu_{k}(t)=C_{nk}\,A_{k}(t) (’random connectivity’ with a connection probability pnk=Cnk/Nkp_{nk}=C_{nk}/N_{k} from population kk to population nn). We assume that all connections from kk to nn have the same weight wnkw_{nk}.

For each population n=1,,Kn=1,\dots,K we write down a Fokker-Planck equation analogous to Eq. (13.16). Neurons are leaky integrate-and-fire neurons. Within a population nn, all neurons have the same parameters τn\tau_{n}, RnR_{n}, urnu_{r}^{n}, and in particular the same firing threshold ϑn\vartheta_{n}. For population nn the Fokker-Planck equation for the evolution of membrane potential densities is then

τntpn(u,t)\displaystyle\tau_{n}{\partial\over\partial t}p_{n}(u,t) =\displaystyle= -u{[-u+RnInext(t)+τnkCnkAk(t)wnk]pn(u,t)}\displaystyle-{\partial\over\partial u}\left\{\left[-u+R_{n}\,I_{n}^{\rm ext}(% t)+\tau_{n}\sum_{k}C_{nk}A_{k}(t)\,w_{nk}\right]\,p_{n}(u,t)\right\} (13.31)
+12[τnkCnkAk(t)wnk2]2u2pn(u,t)\displaystyle+{1\over 2}\,\left[\tau_{n}\sum_{k}C_{nk}A_{k}(t)\,w_{nk}^{2}% \right]\,{\partial^{2}\over\partial u^{2}}p_{n}(u,t)
+τnAn(t)δ(u-urn)-τnAn(t)δ(u-ϑn).\displaystyle+\tau_{n}\,A_{n}(t)\,\delta(u-u_{r}^{n})-\tau_{n}\,A_{n}(t)\,% \delta(u-\vartheta_{n})\,.

The population activity AnA_{n} is the flux through the threshold (cf. Eq. (13.20)), which gives in our case

An(t)=-12[kCnkAk(t)wnk2](pn(u,t)u)u=ϑn.A_{n}(t)=-{1\over 2}\,\left[\sum_{k}C_{nk}A_{k}(t)\,w_{nk}^{2}\right]\,\left({% \partial p_{n}(u,t)\over\partial u}\right)_{u=\vartheta_{n}}\,. (13.32)

Thus populations interact with each other via the variable Ak(t)A_{k}(t); cf. Fig. 13.6.

Fig. 13.6: Interacting Populations. A. Five populations interact via their population activities AkA_{k}. The parameter wnkw_{nk} gives the synaptic weight of a connection from a presynaptic neuron in population kk to a postsynaptic neuron in population nn. Not all populations are coupled with each other. B. Brunel network: An excitatory population of leaky integrate-and-fire neurons is coupled to itself and to an inhibitory population. Neurons in both populations receive also external input from a third population with population activity Aext=νextA_{\rm ext}=\nu^{\rm ext} described as a homogeneous Poisson process.c

Example: Background input

Sometimes it is useful to focus on a single population, coupled to itself, and replace the input arising from other populations by background input. For example, we may focus on population n=1n=1 which is modeled by Eq. (13.31), but use as the inputs AkA_{k} spikes generated by a homogeneous or inhomogeneous Poisson process with rate Ak=νkA_{k}=\nu_{k}.

Example: Full coupling and random coupling

Full coupling can be retrieved by setting Cnk=NkC_{nk}=N_{k} and wnk=Jnk/Nkw_{nk}=J_{nk}/N_{k}. In this case the amount of diffusive noise in population nn

σn2(t)=τn[kCnkAk(t)wnk2]\sigma_{n}^{2}(t)=\tau_{n}\,\left[\sum_{k}C_{nk}A_{k}(t)\,w_{nk}^{2}\right]\, (13.33)

decreases with increasing population size σn21/Nk\sigma^{2}_{n}\propto 1/N_{k}.

13.4.2 Synchrony, oscillations, and irregularity

In Chapter 12 we have already discussed the stationary state of the population activity in a population coupled to itself. With the methods discussed in Chapter 12 we were, however, unable to study the stability of the stationary state; nor were we able to make any prediction about potential time-dependent solutions.

We now give a more complete characterization of a network of leaky integrate-and-fire neurons consisting of two populations: a population of NEN_{E} excitatory neurons coupled to a population with NI=NE/4N_{I}=N_{E}/4 inhibitory neurons (79). The structure of the network is shown in Fig. 13.6B. Each neuron (be it excitatory or inhibitory) receives CEC_{E} connections from the excitatory population each one with weight wEE=wIE=w0w_{EE}=w_{IE}=w_{0}; it also receives CI=CE/4C_{I}=C_{E}/4 connections from the inhibitory population (wEI=wII=-gw0w_{EI}=w_{II}=-g\,w_{0}) and furthermore CEC_{E} connections from an external population (weight wEEw_{EE}) with neurons that fire at a fixed rate νext\nu^{\rm ext}. Each spike causes, after a delay of Δ=1.5ms\Delta=1.5ms a voltage jump of w0=0.1w_{0}=0.1mV and the threshold is 20mV above resting potential. Note that each neuron receives four times as many excitatory than inhibitory inputs so that the total amount of inhibition balances excitation if inhibition is four times stronger (g=4g=4), but the relative strength gg is kept as a free parameter. Also note that, in contrast to Chapter 12, the weight here has units of voltage and directly gives the amplitude of the voltage jump: w0=ΔuEw_{0}=\Delta u_{E}.

Fig. 13.7: Population of pulse-coupled leaky integrate-and-fire neurons (Brunel network). A. Phase diagram of the Brunel network. The population activity can be in a state of asynchronous irregular (AI), synchronous regular (SR) or synchronous irregular (SI) activity. The horizontal axis indicates the relative importance of inhibition (g=4g=4 corresponds to balance between excitation and inhibition). The vertical axis is the amount of external input each neuron receives. Input I=1I=1 corresponds to a mean input just sufficient to reach the neuronal firing threshold (in the absence of recurrent input from the network). Stability of the asynchronous irregular firing state (AI) breaks down at the dashed or solid lines and can lead to synchronous regular (SR), or synchronous irregular (SI, either fast or slow) activity, or to a near-quiescent state (Q); redrawn after Brunel (79). B. Typical time course of the population activity A(t)A(t) (bottom part of each subgraph, units in kHz) and associated spike patterns (50 randomly chosen neurons are plotted along the vertical axis in the top part of each subgraph, each dot indicates a spike) in different firing regimes. ’SI fast’ refers to fast and regular oscillations, but individual neurons fire irregularly with intervals at random multiples of the oscillation period. Horizontal black bar: 50ms. Simulation of 10 000 excitatory and 2 500 inhibitory neurons with connection probability p=0.1p=0.1. Each spike causes, after a delay of Δ=1.5\Delta=1.5ms a voltage jump of 0.10.1mV; distance from equilibrium potential to threshold is 20mV; absolute refractory period 2ms; membrane time constant 2020ms. Parameters are top left: g=3g=3, input =2 (SR); top right: g=6g=6, input =4 (SI fast); bottom left: g=5g=5, input =2 (AI) bottom right g=4.5g=4.5, input =0.9 (SI slow). Adapted from Brunel (79).

The population can be in a state of asynchronous irregular activity (AI), where neurons in the population fire at different times (’asynchronous’ firing) and the distribution of interspike intervals is fairly broad (’irregular’ firing of individual neurons); cf. Fig. 13.7B. This is the state that corresponds to the considerations of stationary activity in Chapter 12. However, with a slight change of parameters, the same network can also be in a state of fast synchronous regular (SR) oscillations. It is characterized by periodic oscillations of the population activity and a sharply peaked interval distribution of individual neurons. The network can also be in the state of synchronous irregular firing (SI) either with fast (SI fast) or with slow (SI slow) oscillations of the population activity. The oscillatory temporal structure emerges despite the fact that the input has a constant spike arrival rate. It can be traced back to an instability of the asynchronous firing regime toward oscillatory activity.

With the mathematical approach of the Fokker-Planck equations, it is possible to determine the instabilities analytically. The mathematical approach will be presented in Section 13.5.2, but Fig. 13.7A already shows the result. Stability of a stationary state of asynchronous firing is most easily achieved in the regime where inhibition dominates excitation. Since each neuron receives four times as many excitatory as inhibitory inputs, inhibition must be at least four times as strong (g>4g>4) as excitation. In order to get non-zero activity despite the strong inhibition, the external input alone must be sufficient to make the neurons fire. ’Input=1’ corresponds to an average external input just sufficient to reach the firing threshold, without additional input from the network.

Consider now the regime of strong inhibition (g>4g>4) and strong input, say spikes from the external source arrive at a rate leading to a mean input of amplitude 4. To understand how the network can run into an instability, let us consider the following intuitive argument. Suppose a momentary fluctuation leads to an increase in the total amount of activity in the excitatory population. This causes, after a transmission delay Δ\Delta an increase in inhibition and, after a further delay Δ\Delta a suppression of excitation. If this feedback loop is strong enough, an oscillation with period 4Δ4\Delta may appear, leading to fast-frequency (’SI fast’) oscillations (79).

If the external input is, on its own, not sufficient to keep the network going (’input <1<1’), then a similar small fluctuation may eventually turn the network from a self-activated state into a quiescent state where the only activity is that caused by external input. It then needs another fluctuation (caused by variations of spike arrivals from the external source) to kick it back into a short burst of activity (79). This leads to slow irregular, but synchronous bursts of firing (’SI slow’).

Finally, if inhibition is globally weak compared to excitation (g<4g<4), then the population is in a state of high activity where each neuron fires close to its maximal rate (set by the inverse of an absolute refractory period). This high-activity state is unstable compared to fast, but regular oscillations (’SR’). Typically, the population can split up into two or more subgroups, the number of which depends on the transmission delay (79; 181).

Example: Analysis of the Brunel network

All neurons have the same parameters, so that we can assume that they all fire at the same time-dependent firing rate ν(t)=A(t)\nu(t)=A(t). We assume that the network is large NCEN\gg C_{E}, so that it is unlikely that neurons share a large fraction of presynaptic neurons; therefore inputs can be considered as uncorrelated, except for the trivial correlations induced by their common modulation of the firing rate A(t)A(t).

The mean input at time tt to a neuron in either the excitatory or the inhibitory population is (in voltage units)

μ(t)=τmw0CE[1-g/4]A(t-Δ)+τmw0CEνext,\mu(t)=\tau_{m}\,w_{0}\,C_{E}\,[1-g/4]\,A(t-\Delta)\,+\tau_{m}\,w_{0}\,C_{E}\,% \nu^{\rm ext}, (13.34)

and the input also generates diffusive noise (in voltage units) of strength

σ2(t)=τmw02CE[(1+g2/4)A(t-Δ)+νext].\sigma^{2}(t)=\tau_{m}\,{w_{0}^{2}\,C_{E}}[(1+g^{2}/4)A(t-\Delta)+\nu_{\rm ext% }]\,. (13.35)

The mean and σ2\sigma^{2} are inserted in the Fokker-Planck equation (cf. Eq. 13.31)

τmtp(u,t)\displaystyle\tau_{m}{\partial\over\partial t}p(u,t) =\displaystyle= -u{[-u+μ(t)]p(u,t)}\displaystyle-{\partial\over\partial u}\left\{\left[-u+\mu(t)\right]\,p(u,t)\right\} (13.36)
+12σ2(t)2u2p(u,t)\displaystyle+{1\over 2}\,\sigma^{2}(t)\,{\partial^{2}\over\partial u^{2}}p(u,t)
+τmA(t)δ(u-ur)-τmA(t)δ(u-ϑ)+𝒪(ωk3).\displaystyle+\tau_{m}\,A(t)\,\delta(u-u_{r})-\tau_{m}\,A(t)\,\delta(u-% \vartheta)+{\mathcal{O}}(\omega_{k}^{3})\,.

The main difference to the stationary solution considered in Section 13.3.1 and Chapter 12 is that mean input and noise amplitude are kept time-dependent.

We first solve to find the stationary solution of the Fokker-Planck equation, denoted as p0(u)p_{0}(u) and activity A0A_{0}. This step is completely analogous to the calculation in Section 13.3.1.

In order to analyze the stability of the stationary solution, we search for solutions of the form A(t)=A0+A1eλtcos(ωt)A(t)=A_{0}+A_{1}\,e^{\lambda t}\,\cos(\omega\,t) (for details, see Section 13.5.2). Parameter combinations that lead to a value λ>0\lambda>0 indicate an instability of the stationary state of asynchronous firing for a network with this set of parameters.

Figure 13.7A indicates that there is a broad regime of stability of the asynchronous irregular firing state (AI). This holds true even if the network is completely deterministic (79) and driven by a constant input I0I_{0} (i.e., we set νext=0\nu^{\rm ext}=0 in the noise term, Eq. (13.35), and replace spike arrivals by a constant current RI0ext=w0CEνextRI_{0}^{\rm ext}=w_{0}C_{E}\,\nu^{\rm ext} in Eq. (13.34).

Stability of the stationary state of asynchronous irregular activity is most easily achieved if the network is in the regime where inhibition slightly dominates excitation and external input is sufficient to drive neurons to threshold. This is called the ‘inhibition dominating’ regime. The range of parameters where asynchronous firing is stable can, however, be increased into the range of dominant excitation, if delays are not fixed at D=1.5D=1.5ms, but drawn from range 0<D<30<D<3ms. Finite size effects can also be taken into account (79).