Next: 6.4 Asynchronous firing Up: 6. Population Equations Previous: 6.2 Density Equations

Subsections

# 6.3 Integral Equations for the Population Activity

In this section, we derive, starting from a small set of assumptions, an integral equation for the population activity. The essential idea of the mathematical formulation is that we work as much as possible on the macroscopic level without reference to a specific model of neuronal dynamics. We will see that the interval distribution PI(t|) that has already been introduced in Chapter 5 plays a central role in the formulation of the population equation. Both the activity variable A and the interval distribution PI(t|) are macroscopic' spike-based quantities that could, in principle, be measured by extracellular electrodes. If we have access to the interval distribution PI(t|) for arbitrary input I(t), then this knowledge is enough to formulate the population equations. In particular, there is no need to know anything about the internal state of the neuron, e.g., the current values of the membrane potential of the neurons.

Integral formulations of the population dynamics have been developed by (Gerstner, 1995,2000b; Wilson and Cowan, 1972; Gerstner and van Hemmen, 1992; Knight, 1972a). The integral equation (6.44) that we have derived in Section 6.2 via integration of the density equations turns out to be a specific instance of the general framework presented in this section.

## 6.3.1 Assumptions

We consider a homogeneous and fully connected network of spiking neurons in the limit of N. We aim for a dynamic equation that describes the evolution of the population activity A(t) over time.

We have seen in Eq. (6.8) that, given the population activity A(t') and the external input Iext(t') in the past, we can calculate the current input potential hPSP(t|) of a neuron that has fired its last spike at , but we have no means yet to transform the potential hPSP back into a population activity. What we need is therefore another equation which allows us to determine the present activity A(t) given the past. The equation for the activity dynamics will be derived from three observations:

(i)
The total number of neurons in the population remains constant. We exploit this fact to derive a conservation law.

(ii)
The model neurons are supposed to show no adaptation. According to Eq. (6.6), the state of neuron i depends explicitly only on the most recent firing time (and of course on the input hPSP), but not on the firing times of earlier spikes of neuron i. This allows us to work in the framework of an input-dependent renewal theory; cf. Chapter 5. In particular, the probability density that neuron i fires again at time t given that its last spike was at time and its input for t't was I(t') is given by the input-dependent interval distribution PI(t|).

(iii)
On a time scale t that is shorter than the axonal transmission delay, all N neurons in the population can be considered as independent. The number of spikes nact(t;t + t) that the network emits in a short time window t is therefore the sum of independent random variables that, according to the law of large numbers, converges (in probability) to its expectation value. For a large network it is thus sufficient to calculate expectation values.

## 6.3.2 Integral equation for the dynamics

Because of observation (ii) we know that the input-dependent interval distribution PI(t | ) contains all relevant information. We recall that PI(t | ) gives the probability density that a neuron fires at time t given its last spike at and an input I(t') for t' < t. Integration of the probability density over time PI(s | ) ds gives the probability that a neuron which has fired at fires its next spike at some arbitrary time between and t. Just as in Chapter 5, we can define a survival probability,

 SI(t | ) = 1 - PI(s | ) ds , (6.71)

i.e., the probability that a neuron which has fired its last spike at survives' without firing up to time t.

We now return to the homogeneous population of neurons in the limit of N and use observation (iii). We consider the network state at time t and label all neurons by their last firing time . The proportion of neurons at time t which have fired their last spike between t0 and t1 < t (and have not fired since) is expected to be

 = SI(t | ) A() d . (6.72)

For an interpretation of the integral on the right-hand side of Eq. (6.72), we recall that A() is the fraction of neurons that have fired in the interval [, +  ]. Of these a fraction SI(t|) are expected to survive from to t without firing. Thus among the neurons that we observe at time t the proportion of neurons that have fired their last spike between t0 and t1 is expected to be SI(t | A() d ; cf. Fig. 6.7.

Finally, we make use of observation (i). All neurons have fired at some point in the past6.2. Thus, if we extend the lower bound t0 of the integral on the right-hand side of Eq. (6.72) to - and the upper bound to t, the left-hand side becomes equal to one,

 1 = SI(t | ) A() d , (6.73)

because all N neurons have fired their last spike in the interval [- , t]. Since the number of neurons remains constant, the normalization (6.73) must hold at arbitrary times t. Eq. (6.73) is an implicit equation for the population activity A. It is the starting point for the discussions in this and the following chapters.

Since Eq. (6.73) is rather abstract, we try to put it into a form that is easier to grasp intuitively. To do so, we take the derivative of Eq. (6.73) with respect to t. We find

 0 = SI(t| t) A(t) +  A() d . (6.74)

We now use PI(t|) = - SI(t|) and SI(t| t) = 1 which is a direct consequence of Eq. (6.71). This yields the activity dynamics

 A(t) = PI(t |) A() d (6.75)

A different derivation of Eq. (6.75) has been given in Section 6.2; cf. Eq. (6.44).

Eq. (6.75) is easy to understand. The kernel PI(t | ) is the probability density that the next spike of a neuron which is under the influence of an input I occurs at time t given that its last spike was at . The number of neurons which have fired at is proportional to A() and the integral runs over all times in the past. The interval distribution PI(t|) depends upon the total input (both external input and synaptic input from other neurons in the population) and hence upon the postsynaptic potential (6.8). Eqs. (6.8) and (6.75) together with an appropriate noise model yield a closed system of equations for the population dynamics. Eq. (6.75) is exact in the limit of N. Corrections for finite N, have been discussed by Meyer and van Vreeswijk (2001) and Spiridon and Gerstner (1999).

An important remark concerns the proper normalization of the activity. Since Eq. (6.75) is defined as the derivative of Eq. (6.73), the integration constant on the left-hand side of (6.73) is lost. This is most easily seen for constant activity A(t) A0. In this case the variable A0 can be eliminated on both sides of Eq. (6.75) so that Eq. (6.75) yields the trivial statement that the interval distribution is normalized to unity. Equation (6.75) is therefore invariant under a rescaling of the activity A0 c A0. with any constant c. To get the correct normalization of the activity A0 we have to go back to Eq. (6.73).

We conclude this section with a final remark on the form of Eq. (6.75). Even though (6.75) looks linear, it is in fact a highly non-linear equation because the kernel PI(t | ) depends non-linearly on hPSP, and hPSP in turn depends on the activity via Eq. (6.8).

### 6.3.2.1 Absolute Refractoriness and the Wilson-Cowan integral equation

Let us consider a population of Poisson neurons with an absolute refractory period . A neuron that is not refractory, fires stochastically with a rate f[h(t)] where h(t) is the total input potential, viz., the sum of the postsynaptic potentials caused by presynaptic spike arrival. After firing, a neuron is inactive during the time . The population activity of a homogeneous group of Poisson neurons with absolute refractoriness is (Wilson and Cowan, 1972)

 A(t) = f[h(t)] 1 - A(t') dt' (6.76)

We will show below that Eq. (6.76) is a special case of the population activity equation (6.75).

The Wilson-Cowan integral equation (6.76) has a simple interpretation. Neurons stimulated by a total postsynaptic potential h(t) fire with an instantaneous rate f[h(t)]. If there were no refractoriness, we would expect a population activity A(t) = f[h(t)]. Not all neurons may, however, fire since some of the neurons are in the absolute refractory period. The fraction of neurons that participate in firing is 1 - A(t') dt' which explains the factor in curly brackets.

For constant input potential, h(t) = h0, the population activity has a stationary solution,

 A0 = = g(h0) . (6.77)

For the last equality sign we have used the definition of the gain function of Poisson neurons with absolute refractoriness in Eq. (5.40). Equation (6.77) tells us that in a homogeneous population of neurons the population activity in a stationary state is equal to the firing rate of individual neurons. This is a rather important result since it allows us to calculate the stationary population activity from single-neuron properties. A generalization of Eq. (6.77) to neurons with relative refractoriness will be derived in Section 6.4.

The function f in Eq. (6.76) was motivated by an instantaneous escape rate' due to a noisy threshold in a homogeneous population. In this interpretation, Eq. (6.76) is the exact equation for the population activity of neurons with absolute refractoriness. In their original paper, Wilson and Cowan motivated the function f by a distribution of threshold values in an inhomogeneous population. In this case, the population equation (6.76) is an approximation, since correlations are neglected (Wilson and Cowan, 1972).

### 6.3.2.2 Derivation of the Wilson-Cowan integral equation (*)

We apply the population equation (6.75) to SRM0 neurons with escape noise; cf. Chapter 5. The escape rate f (u - ) is a function of the distance between the membrane potential and the threshold. For the sake of notational convenience, we set = 0. The neuron model is specified by a refractory function as follows. During an absolute refractory period 0 < s, we set (s) to - . For s, we set (s) = 0. Thus the neuron exhibits absolute refractoriness only; cf. Eq. (5.57). The total membrane potential is u(t) = (t - ) + h(t) with

 h(t) = J0(s) A(t - s) ds + (s) Iext(t - s) ds (6.78)

Given it seems natural to split the integral in the activity equation (6.75) into two parts

 A(t) = PI(t | ) A() d + PI(t | ) A() d . (6.79)

The interval distribution PI(t | ) for the noisy threshold model has been derived in Chapter 5 and is repeated here for convenience

 PI(t | ) = f[h(t) + (t - )] exp - f[h(t') + (t' - )]dt' . (6.80)

Let us evaluate the two terms on the right-hand side of Eq. (6.79). Since spiking is impossible during the absolute refractory time, i.e., f[- ] = 0, the second term in Eq. (6.79) vanishes. In the first term we can move a factor f[h(t) + (t - )] = f[h(t)] in front of the integral since vanishes for t - > . The exponential factor is the survivor function of neurons with escape noise as defined in Eq. (5.7); cf. Chapter 5. Therefore Eq. (6.79) reduces to

 A(t) = f[h(t)]SI(t | ) A() d (6.81)

Let us now recall the normalization of the population activity defined in Eq. (6.73), i.e.,

 SI(t | ) A() d = 1 . (6.82)

The integral in (6.81) can therefore be rewritten as

 A(t) = f[h(t)] 1 - SI(t | ) A() d . (6.83)

During the absolute refractory period we have a survival probability SI(t | ) = 1 since the neurons can not fire. This yields

 A(t) = f[h(t)] 1 - A(t') dt' , (6.84)

which is the Wilson-Cowan integral equation (6.76) for neurons with absolute refractoriness (Wilson and Cowan, 1972).

### 6.3.2.3 Quasi-stationary dynamics (*)

Integral equations are often difficult to handle. Wilson and Cowan aimed therefore at a transformation of their equation into a differential equation (Wilson and Cowan, 1972). To do so, they had to assume that the population activity changes slowly during the time and adopted a procedure of time coarse-graining. Here we present a slightly modified version of their argument (Pinto et al., 1996; Gerstner, 1995).

 h(t) = J0(s) A(t - s) ds + (s) Iext(t - s) ds , (6.85)

contains the low-pass filters' and . If the postsynaptic potential is broad, h changes only slowly. In particular the dynamics of h is limited by the membrane time constant . The activity A(t), however, could still change rapidly. As a first, and essential, step in the procedure of `time coarse-graining', we simply assume that A changes only slowly over the time . In this case the integral A(t') dt' on the right-hand side of Eq. (6.76) can be approximated by A(t. With this approximation, we can solve (6.76) for A and find

 A(t) = = : g[h(t)] (6.86)

In other words, time coarse-graining implies that we replace the instantaneous activity A(t) by its stationary value (6.77). As an aside we note that the activity A is given as a function of the input potential h rather than the input current.

As a second step, we transform Eq. (6.85) into a differential equation. If the response kernels are exponentials, i.e., (s) = (s) =  exp(- s/), the derivative of Eq. (6.85) is

 = - h(t) + J0 A(t) + Iext(t) . (6.87)

The evolution of the activity A is given by Eq. (6.86) and follows the evolution of h as defined by Eq. (6.87). In particular, the coarse-grained activity A(t) cannot be faster than the membrane potential h and is therefore limited by the membrane time constant . If the membrane time constant is much larger than the refractory time , the approximation (6.86) is good. For , the initial transient of the population activity after a change in the input Iext can be faster than . The difference between the numerical solution of the Wilson-Cowan integral equation (6.76) and that of the coarse-grained equation (6.86) is shown in Fig. 6.8B for a step current input. Whereas the approximative solution approaches the new stationary solution asymptotically from below, the solution of the integral equation exhibits an oscillatory component. In particular, the initial response during the first 2ms is faster than that of the approximative solution. In Chapter 7 we will study transient responses in more detail and show that in the limit of low noise, the population activity can respond much faster than the membrane potential.

Equation (6.87) is a differential equation for the membrane potential. Alternatively, the integral equation (6.76) can also be approximated by a differential equation for the activity A. We start from the observation that in a stationary state the activity A can be written as a function of the input current, i.e, A0 = g(I0) where I0 = Iext + J0 A0 is the sum of the external driving current and the postsynaptic current caused by lateral connections within the population. What happens if the input current changes? The population activity of neurons with a large amount of escape noise will not react to rapid changes in the input instantaneously, but follow slowly with a certain delay similar to the characteristics of a low-pass filter. An equation that qualitatively reproduces the low-pass behavior is

 = - A(t) + gIext + J0 A(t)  . (6.88)

This is the Wilson-Cowan differential equation. Note that, in contrast to Eq. (6.87) the sum Iext + J0 A(t) appears inside the argument of the gain function g. If synaptic dynamics is slow, the term J0 A(t) should be replaced by J0(sA(t - s) ds where (s) is the time course of the postsynaptic current; cf. Eq. (5.138). The time constant is arbitrary but is often identified with the membrane time constant .

### 6.3.2.4 Example: Finite refractoriness and escape noise (*)

The Wilson-Cowan integral equation that has been discussed above is valid for neurons with absolute refractoriness only. We now generalize some of the arguments to a Spike Response Model with relative refractoriness. We suppose that refractoriness is over after a time so that (s) = 0 for s. For 0 < s < , the refractory kernel may have any arbitrary shape. Furthermore we assume that for t > the kernels (t - , s) and (t - , s) do not depend on t - . For 0 < t - < we allow for an arbitrary time-dependence. Thus, the postsynaptic potential is

 (t - , s) =  , (6.89)

and similarly for . As a consequence, the input potential hPSP(t|) defined in Eq. (6.8) depends on for t - , but is independent of for t - > .

We start from the population equation (6.75) and split the integral into two parts

 A(t) = PI(t |) A() d  + PI(t |) A() d (6.90)

Since (s) vanishes for s, refractoriness plays no role in the second term on the right-hand side of (6.90). We assume escape noise with an escape rate f and use the same arguments as with Eqs. (6.81) - (6.83). The result is

 A(t) = PI(t |) A() d  + f[h(t)] 1 - SI(t|) A() d ; (6.91)

cf. Wilson and Cowan (1972), Gerstner (2000b). Here we have used hPSP(t|) = h(t) for t - > . The term in square brackets is the fraction of neurons that are not refractory. These neurons form a homogeneous group and fire with instantaneous rate f[h(t)]. For the term in square brackets vanishes and we retrieve the standard population equation (6.75).

The integrals on the right-hand side of (6.91) have finite support which makes Eq. (6.91) more convenient for numerical implementation than the standard formulation (6.75). For a rapid implementation scheme, it is convenient to introduce discretized refractory densities as discussed in Section 6.2; cf. Eq. (6.57).

Next: 6.4 Asynchronous firing Up: 6. Population Equations Previous: 6.2 Density Equations
Gerstner and Kistler
Spiking Neuron Models. Single Neurons, Populations, Plasticity
Cambridge University Press, 2002