In order to formulate a first version of an integral equation for a population of neurons, we start in the framework of time-dependent renewal theory (Chapter 7). To do so, we have to assume the state of a neuron $i$ at time $t$ to be completely described by (i) its last firing time $\hat{t}_{i}$; (ii) the input $I(t^{\prime})$ it received for times $t^{\prime}<t$; (iii) the characteristics of potential noise sources, be it noise in the input, noise in the neuronal parameters, or noise in the output.

Are these assumptions overly restrictive or can we still place an interesting and rich set of neuron models in the class of models consistent with assumptions (i) - (iii) of time-dependent renewal theory?

The state of a single-variable integrate-and-fire neuron $i$, for example one of the nonlinear integrate-and-fire models of Chapter 5, is completely characterized by its momentary membrane potential $u_{i}(t)$. After firing at a time $\hat{t}_{i}$, the membrane potential is reset, so that all information that arrived before the firing is ‘forgotten’. Integration continues with $u=u_{r}$. Therefore knowledge of last firing time $\hat{t}_{i}$ and the input current $I(t^{\prime})$ for $\hat{t}_{i}<t^{\prime}\leq t$ is sufficient to predict the momentary state of neuron $i$ at time $t$ and therefore (for a deterministic model) its next firing time. If the integrate-and-fire model is furthermore subject to noise in the form of stochastic spike arrivals of known rate (but unknown spike arrival times), we will not be able to predict the neuron’s exact firing time. Instead we are interested in the probability density

that the next spike occurs around time $t$ given that the last spike was at time $\hat{t}_{i}$ and the neuron was subject to an input $I(t^{\prime})$ for $t^{\prime}<t$. In practice, it might be difficult to write down an exact analytical formula for $P_{I}(t|\hat{t}_{i})$ for noise models corresponding to stochastic spike arrival or diffusive noise (cf. Chapter 8), but the quantity $P_{I}(t|\hat{t}_{i})$ is still well defined. We call $P_{I}(t|\hat{t}_{i})$ the generalized interval distribution (cf. Chapter 9).

For several other noise models, there exist explicit mathematical expressions for $P_{I}(t|\hat{t}_{i})$. Consider, for example, integrate-and-fire models with slow noise in the parameters, defined as follows. After each firing time, the value $u_{r}$ for the reset or the value $\tau_{m}$ of the membrane time constant is drawn from a predefined distribution. Between two resets, the membrane potential evolves deterministically. As a consequence, the distribution $P_{I}(t|\hat{t}_{i})$ of the next firing time can be predicted from the distribution of parameters and the deterministic solution of the threshold crossing equations (183). Such a model of slow noise in the parameters can also be considered as an approximation to a heterogeneous population of neurons where different neurons in the population have slightly different, but fixed parameters.

Another tractable noise model is the ‘escape noise’ model which is also the basis for the Generalized Linear Models of spiking neurons already discussed in Chapter 9. The short-term memory approximation SRM${}_{0}$ of the Spike Response Model with escape noise (cf. Eq. (9.19) in Ch. 9) is a prime example of a model that fits into the framework of renewal theory specified by assumptions (i) - (iii).

The escape noise model can be used in combination with any linear or nonlinear integrate-and-fire model. For a suitable choice of the escape function $f$ it provides also an excellent approximation to diffusive noise in the input; cf. Section 9.4 in Ch. 9. An example of how to formulate a nonlinear integrate-and-fire model with escape noise is given below.

The major limitation of time-dependent renewal theory is that, for the moment, we need to exclude adaptation effects, because the momentary state of adaptive neurons not only on the last firing time (and the input in the past) but also on the firing times of earlier spikes. Section 14.1.4 shows that adaptation can be included by extending the renewal theory to a quasi-renewal theory. The treatment of adaptive neurons is important because most cortical neurons exhibit adaptation (Ch. 6).

The integral equation (182; 183; 552) for activity dynamics with time-dependent renewal theory states that the activity at time $t$ depends on the fraction of active neurons at earlier times $\hat{t}$ times the probability of observing a spike at $t$ given a spike at $\hat{t}$:

Equation (14.5) is easy to understand. The kernel $P_{I}(t|\hat{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 $\hat{t}$. The number of neurons which have fired at $\hat{t}$ is proportional to $A(\hat{t})$ and the integral runs over all times in the past. The interval distribution $P_{I}(t|\hat{t})$ depends on the total input (both the external input and the synaptic input from other neurons in the population). For an unconnected population, $I(t)$ corresponds to the external drive. For connected neurons, $I(t)$ is the sum of the external input and the recurrent input from other neurons in the population; the case of connected neurons will be further analyzed in Section 14.2.

We emphasize that $A(t)$ on the left-hand side of Eq. (14.5) is the expected activity at time $t$ while $A(\hat{t})$ on the right-hand side is the observed activity in the past. In the limit of the number $N$ of neurons in the population going to infinity, the fraction of neurons that actually fire in a short time $\Delta t$ is the same as its expected value $A(t)\,\Delta t$. Therefore, Eq. (14.5) becomes exact in the limit of $N\to\infty$, so that we can use the same symbol for $A(t)$ on both sides of the equation. Finite-size effects are discussed in Section 14.6.

We conclude this section with some final remarks on the form of Eq. (14.5). First, we observe that we can multiply the activity value $A$ on both sides of the equation with a constant $c$ and introduce a new variable $A^{\prime}=c\,A$. If $A(t)$ solves the equation, then $A^{\prime}(t)$ will solve it as well. Thus, Eq. (14.5) cannot predict the correct normalization of the activity $A$. This reflects the fact that, instead of defining the activity by a spike count divided by the number $N$ of neurons, we could have chosen to work directly with the spike count per unit of time or any other normalization. The proper normalization consistent with our definition of the population activity $A(t)$ is derived below in Section 14.1.3.

Second, even though Eq. (14.5) is linear in the variable $A$, it is in fact a highly nonlinear equation in the drive, because the kernel $P_{I}(t|\hat{t})$ depends nonlinearly on the input $I(t)$. However, numerical implementations of the integral-equations lead to rapid schemes that predict the population activity in response to changing input, even in the highly nonlinear regime when strong input transiently synchronizes neurons in the population (Fig. 14.1).

To derive Eq. (14.5), we recall that $P_{I}(t|\hat{t})$ is the probability density that a neuron fires at time $t$ given its last spike at $\hat{t}$ and an input $I(t^{\prime})$ for $t^{\prime}<t$. Integration of the probability density over time $\int_{\hat{t}}^{t}P_{I}(s|\hat{t})ds$ gives the probability that a neuron which has fired at $\hat{t}$ fires its next spike at some arbitrary time between $\hat{t}$ and $t$. Just as in Chapter 7, we can define a survival probability,

i.e., the probability that a neuron which has fired its last spike at $\hat{t}$ “survives” without firing up to time $t$.

We now return to the homogeneous population of neurons in the limit of $N\rightarrow\infty$ and assume that the firing of different neurons at time $t$ is independent, given that we know the history (the input and last spike) of each neuron. The technical term for such a situation is ’conditional independence’. We consider the network state at time $t$ and label all neurons by their last firing time $\hat{t}$. The proportion of neurons at time $t$ which have fired their last spike between $t_{0}$ and $t_{1}<t$ (and have not fired since) is expected to be

For an interpretation of the integral on the right-hand side of Eq. (14.7), we recall that $A(\hat{t})d\hat{t}$ is the fraction of neurons that have fired in the interval $[\hat{t},\hat{t}+\Delta\hat{t}]$. Of these a fraction $S_{I}(t|\hat{t})$ are expected to survive from $\hat{t}$ 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 $t_{0}$ and $t_{1}$ is expected to be $\int_{t_{0}}^{t_{1}}S_{I}(t|\hat{t})A(\hat{t})d\hat{t}$; cf. Fig. 14.2.

Finally, we use the fact that the total number of neurons remains constant. All neurons have fired at some point in the past.⁸⁸Neurons which have never fired before are assigned a formal firing time $\hat{t}=-\infty$. Thus, if we extend the lower bound $t_{0}$ of the integral on the right-hand side of Eq. 14.7 to $-\infty$ and the upper bound to $t$, the left-hand side becomes equal to one,

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

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

We now use $P_{I}(t|\hat{t})=-\frac{\partial}{\partial t}S_{I}(t|\hat{t})$ and $S_{I}(t|t)=1$ which is a direct consequence of Eq. (14.6). This yields the activity dynamics of Eq. (14.5).

We repeat an important remark concerning the normalization of the activity. Since Eq. (14.5) is defined as the derivative of Eq. (14.8), the integration constant on the left-hand side of Eq. (14.8) is lost. This is most easily seen for constant activity $A(t)=A_{0}$. In this case the variable $A_{0}$ can be eliminated on both sides of Eq. (14.5) so that it yields the trivial statement that the interval distribution is normalized to unity. Equation (14.5) is therefore invariant under a rescaling of the activity $A_{0}\rightarrow cA_{0}$ with any constant $c$, as mentioned earlier. To get the correct normalization of the activity we have to return to Eq. (14.8).

The integral equations presented so far are restricted to non-adapting neurons, because we assumed that the knowledge of the input and of the last firing time was sufficient to characterize the state of a neuron. For adapting neurons, however, the whole spiking history can shape the interspike interval distribution. In this section we extend the integral equation (14.5) to the case of adapting neurons.

For an isolated adaptive neuron in the presence of noise, the probability density of firing around time $t$ will, in general, depend on its past firing times $\hat{t}_{n}<\hat{t}_{n-1}\dots<\hat{t}_{2}<\hat{t}_{1}=\hat{t}<t$ where $\hat{t}=\hat{t}_{1}$ denotes the most recent spike time. The central insight is that, in a population of neurons, we can approximate the past firing times by the population activity $A(\hat{t})$ in the past. Indeed, the probability of one of the neurons having a past spike time $\hat{t}_{k}=t^{\prime}$ in the interval $\hat{t}<t^{\prime}<\hat{t}+\Delta t$ is $A(\hat{t})\Delta t$.

Just as in the time-dependent renewal theory, we will treat the most recent spike $\hat{t}_{1}=\hat{t}$ of each neuron explicitly. For all spikes $\hat{t}_{k}$ with $k\geq 2$ we approximate the actual spike train history of an individual by the average history as summarized in the population activity $A(t)$. Let $P_{I,A}(t|\hat{t})$ be the probability of observing a spike at time $t$ given the last spike at time $\hat{t}$, the input-current and the activity history $A(t)$ until time $t$, then we can rewrite equation (14.5) in a format suitable for adaptive neurons

In Section 14.5 we explain the methods to describe populations of adaptive neurons in more detail. Appropriate numerical schemes (Section 14.1.5) make it possible to describe the response of a population of adaptive neurons to arbitrary time-dependent input; see Fig. 14.4.

Integral equations can be solved numerically. However, the type of equation we face in Eq. (14.5) or Eq. (14.12) can not be cast in the typical Volterra or Fredholm forms for which efficient numerical methods have been extensively studied (301; 31). In the following, we describe a method that can be used to solve Eq. (14.5), Eq. (14.12) or Eq. (14.8), but we take as an example the quasi-renewal equivalent of Eq. (14.8)

To derive a numerical algorithm, we proceed in three steps. First, truncation of the integral in Eq. (14.13) at some lower bound. Second, the discretization of this integral and, third, the discretization of the integral defining the survivor function. Here we will use the rectangle rule for the discretization of every integral. Note that adaptive quadratures or Monte-Carlo methods could present more efficient alternatives.

The first step is the truncation. The probability to survive a very long time is essentially zero. Let $\tau_{c}$ be a period of time such that the survivor $S_{I,A}(t|t-\tau_{c})$ is very small. Then, we can truncate the infinite integral in 14.13

Next, we proceed to discretize the integral on small bins of size $\Delta t$. Let $\mathbf{m}^{(t)}$ be the vector made of the fraction of neurons at $t$ with last spike within $\hat{t}$ and $\hat{t}+\Delta t$. This definition means that the $k$th element is $m_{k}^{(t)}=S(t|t-k\,\Delta t)A(t-k\,\Delta t)\Delta t$. The element with $k=0$ is then the momentary population activity in the time step starting at time $t$, $m^{(t)}_{0}=A_{t}\Delta t$, since $S(t|t)=1$. Therefore we arrive at

Finally, we discretize the integral defining the survival function to obtain $m_{k}^{(t)}$ as a function of $m_{k}^{(t-\Delta t)}$. Because of $S(t|\hat{t})=\exp[-\int_{\hat{t}}^{t}\rho(t^{\prime}|\hat{t}){\text{d}}t^{\prime}]=\exp[-\int_{t-\Delta t}^{t}\rho(t^{\prime}|\hat{t}){\text{d}}t^{\prime}]S(t-\Delta t|\hat{t})$, we find for sufficiently small time steps

Note that we work in a moving coordinate because the index $k$ always counts time steps backward starting from the present time $t$. Therefore $m_{k}^{(t)}$ and $m_{k-1}^{(t-\Delta t)}$ refer to the same group of neurons, i.e., those that have fired their last spike around time $\hat{t}=t-k\Delta t=(t-\Delta t)-(k-1)(\Delta t))$. Eq. (14.16) indicates that the number of neurons who survive from $\hat{t}$ up to $t$ decrease in each time step by an exponential factor. In Eq. (14.16), $\rho(t|t^{\prime})$ can be either the renewal conditional intensity of a non-adaptive neuron model (e.g., Eq. (14.3)) or the effective ’quasi-renewal’ intensity of adaptive neurons (see Section 14.5). Together, Eq. (14.15) and (14.16) can be used to solve $A_{t}$ iteratively.