14 The Integral-equation Approach

14.4 Density equations vs. integral equations

In this section we relate the integral equation (14.5) to the membrane potential density approach for integrate-and-fire neurons that we have discussed in Chapter 13. The two approaches are closely related. For noise-free neurons driven by a constant supra-threshold stimulus, the two mathematical formulations are, in fact, equivalent and related by a simple change of variables. Even for noisy neurons with subthreshold stimulation, the two approaches are comparable. Both methods are linear in the densities and amenable to efficient numerical implementations. The formal mathematical relation between the two approaches is shown in Section 14.4.3.

We use Section 14.4.1 to introduce a ’refractory density’ which is analogous to the membrane potential density that we have seen in Chapter 13. In particular, the dynamics of refractory density variables follow a continuity equation. The solution of the continuity equation (14.4.2) transforms the partial differential equation into the integral equations that we have seen in Section 14.1.

14.4.1 Refractory densities

In this section we develop a density formalism for spike response neurons, similar to the membrane potential density approach for integrate-and-fire neurons that we have discussed in Chapter 13. The main difference is that we replace the membrane potential density p(u,t)p(u,t) by a refractory density q(r,t)q(r,t), to be introduced below.

We study a homogeneous population of neurons with escape noise. The neuron model should be consistent with time-dependent renewal theory. In this case, we can write the membrane potential of a neuron ii in general form as

ui(t,t^)=η(t-t^i)+h(t,t^i).u_{i}(t,\hat{t})=\eta(t-\hat{t}_{i})+h(t,\hat{t}_{i})\,. (14.61)

For example, a nonlinear integrate-and-fire model which has fired that last time at t^i\hat{t}_{i} and follows a differential equation u˙(t)=[f(u)+RI(t)]\dot{u}(t^{\prime})=[f(u)+R\,I(t^{\prime})] has for t>t^it>\hat{t}_{i} a potential ui(t,t^i)=h(t,t^i)=ur+t^itu˙(t)dtu_{i}(t,\hat{t}_{i})=h(t,\hat{t}_{i})=u_{r}+\int_{\hat{t}_{i}}^{t}\dot{u}(t^{% \prime}){\text{d}}t^{\prime} and η=0\eta=0. A SRM0{}_{0} neuron has a potential ui(t,t^)=η(t-t^i)+h(t)u_{i}(t,\hat{t})=\eta(t-\hat{t}_{i})+h(t) with an arbitrary spike-afterpotential η(t-t^i)\eta(t-\hat{t}_{i}) and an input potential hh.

The notation in Eq. (14.61) emphasizes the importance of the last firing time t^\hat{t}. We denote the refractory state of a neuron by the variable

r=t-t^0,r=t-\hat{t}\geq 0\,, (14.62)

i.e., by the time that has passed since the last spike. If we know rr and the total input current in the past, we can calculate the membrane potential, u(t)=η(r)+h(t,t-r)u(t)=\eta(r)+h(t,t-r). Given the importance of the refractory variable rr, we may wonder how many neurons in the population have a refractory state between r0r_{0} and r0+Δrr_{0}+\Delta r. For a large population (NN\to\infty) the fraction of neurons with a momentary value of rr in the interval [r0,r0+Δr][r_{0},r_{0}+\Delta r] is given by

limN{neurons with r0<r(t)r0+ΔrN}=r0r0+Δrq(r,t)dr,\lim_{N\to\infty}\left\{{\text{neurons with }r_{0}<r(t)\leq r_{0}+\Delta r% \over N}\right\}=\int_{r_{0}}^{r_{0}+\Delta r}q(r,t)\,{\text{d}}r\,, (14.63)

where q(r,t)q(r,t) is the refractory density. The aim of this section is to describe the dynamics of a population of SRM neurons by the evolution of q(r,t)q(r,t).

We start from the continuity equation (cf. Ch. 13),

tq(r,t)=-rJrefr(r,t),{\partial\over\partial t}q(r,t)=-{\partial\over\partial r}J_{\rm refr}(r,t)\,, (14.64)

where we have introduced the flux JrefrJ_{\rm refr} along the axis of the refractory variable rr. As long as the neuron does not fire, the variable r=t-t^r=t-\hat{t} increases at a speed of dr/dt=1{\text{d}}r/{\text{d}}t=1. The flux is the density qq times the velocity, hence

Jrefr(r,t)=q(r,t)drdt=q(r,t).J_{\rm refr}(r,t)=q(r,t)\,{{\text{d}}r\over{\text{d}}t}=q(r,t)\,. (14.65)

The continuity equation (14.64) expresses the fact that, as long as a neuron does not fire, its trajectories r(t)r(t) can neither start nor end. On the other hand, if a neuron fires, the trajectory stops at the current value of rr and ‘reappears’ at r=0r=0. In the escape rate formalism, the instantaneous firing rate of a neuron with refractory variable rr is given by the hazard

ρ(t|t-r)=f[η(r)+h(t|t-r)].\rho(t|t-r)=f[\eta(r)+h(t|t-r)]\,. (14.66)

If we multiply the hazard (14.66) with the density q(r,t)q(r,t), we get the loss per unit of time,

Jloss=-ρ(t|t-r)q(r,t).J_{\rm loss}=-\rho(t|t-r)\,q(r,t)\,. (14.67)

The total number of trajectories that disappear at time tt due to firing is equal to the population activity, i.e.,

A(t)=0ρ(t|t-r)q(r,t)dr.A(t)=\int_{0}^{\infty}\rho(t|t-r)\,q(r,t){\text{d}}r\,. (14.68)

The loss (14.67) has to be added as a ‘sink’ term on the right-hand side of the continuity equation, while the activity A(t)A(t) acts as a source at r=0r=0. The full dynamics is

tq(r,t)=-[rq(r,t)]-ρ(t|t-r)q(r,t)+δ(r)A(t).{\partial\over\partial t}q(r,t)=-\left[{\partial\over\partial r}q(r,t)\right]-% \rho(t|t-r)\,q(r,t)+\delta(r)\,A(t)\,. (14.69)

This partial differential equation is the analog of the Fokker-Planck equation (13.16) for the membrane potential density of integrate-and-fire neurons. The relation between the two equations will be discussed in Section 14.4.3.

Equation (14.69) can be integrated analytically and rewritten in the form of an integral equation for the population activity. The mathematical details of the integration will be presented below. The final result is

A(t)=-tPI(t|t^)A(t^)dt^,A(t)=\int_{-\infty}^{t}P_{I}(t|\hat{t})\,A(\hat{t})\,{\text{d}}\hat{t}\,, (14.70)

where

PI(t|t^)=ρ(t|t^)exp[-t^tρ(t|t^)dt]P_{I}(t|\hat{t})=\rho(t|\hat{t})\,\exp\left[-\int_{\hat{t}}^{t}\rho(t^{\prime}% |\hat{t})\,{\text{d}}t^{\prime}\right] (14.71)

is the interval distribution of neurons with escape noise; cf. Eq. (7.28). Thus, neurons that have fired their last spike at time t^\hat{t} contribute with weight PI(t|t^)P_{I}(t|\hat{t}) to the activity at time tt. Integral equations of the form (14.70) are the starting point for the formal theory of population activity presented in Section 14.1.

14.4.2 From refractory densities to the integral equation (*)

All neurons that have fired together at time t^\hat{t} form a group that moves along the rr-axis at constant speed. To solve Eq. (14.69) we turn to a frame of reference that moves along with the group. We replace the variable rr by t-rt^t-r\equiv\hat{t} and define a new density

Q(t^,t)=q(t-t^,t),Q(\hat{t},t)=q(t-\hat{t},t)\,, (14.72)

with t^t\hat{t}\leq t. The total derivative of QQ with respect to tt is

ddtQ(t^,t)=rq(r,t)|r=t-t^drdt+tq(r,t)|r=t-t^{{\text{d}}\over{\text{d}}t}Q(\hat{t},t)={\partial\over\partial r}\left.q(r,t)% \right|_{r=t-\hat{t}}\,{{\text{d}}r\over{\text{d}}t}+{\partial\over\partial t}% \left.q(r,t)\right|_{r=t-\hat{t}}\, (14.73)

with dr/dt=1{\text{d}}r/{\text{d}}t=1. We insert Eq. (14.69) on the right-hand side of (14.73) and obtain for t>t^t>\hat{t}

ddtQ(t^,t)=-ρ(t|t^)Q(t^,t).{{\text{d}}\over{\text{d}}t}Q(\hat{t},t)=-\rho(t|\hat{t})\,Q(\hat{t},t)\,. (14.74)

The partial differential equation (14.69) has thus been transformed into an ordinary differential equation that is solved by

Q(t^,t)=Q(t^,t0)exp[-t0tρ(t|t^)dt],Q(\hat{t},t)=Q(\hat{t},t_{0})\,\exp\left[-\int_{t_{0}}^{t}\rho(t^{\prime}|\hat% {t})\,{\text{d}}t^{\prime}\right]\,, (14.75)

where Q(t^,t0)Q(\hat{t},t_{0}) is the initial condition, which is still to be fixed.

From the definition of the refractory density q(r,t)q(r,t) we conclude that q(0,t)q(0,t) is the proportion of neurons at time tt that have just fired, i.e., q(0,t)=A(t)q(0,t)=A(t) or, in terms of the new refractory density, Q(t,t)=A(t)Q(t,t)=A(t). We can thus fix the initial condition in Eq. (14.75) at t0=t^t_{0}=\hat{t} and find

Q(t^,t)=A(t^)exp[-t^tρ(t|t^)dt].Q(\hat{t},t)=A(\hat{t})\,\exp\left[-\int_{\hat{t}}^{t}\rho(t^{\prime}|\hat{t})% \,{\text{d}}t^{\prime}\right]\,. (14.76)

On the other hand, from (14.68) we have

A(t)=-tρ(t|t^)Q(t^,t)dt^.A(t)=\int_{-\infty}^{t}\rho(t|\hat{t})\,Q(\hat{t},t)\,{\text{d}}\hat{t}\,. (14.77)

If we insert (14.76) into (14.77), we find

A(t)=-tρ(t|t^)exp[-t^tρ(t|t^)dt]A(t^)dt^A(t)=\int_{-\infty}^{t}\rho(t|\hat{t})\,\exp\left[-\int_{\hat{t}}^{t}\rho(t^{% \prime}|\hat{t})\,{\text{d}}t^{\prime}\right]\,A(\hat{t})\,{\text{d}}\hat{t} (14.78)

which is the population equation (14.70) mentioned above. It has been derived from refractory densities in Gerstner and van Hemmen (180) and similarly in the appendix of the papers by Wilson and Cowan (552). If we insert Eq. (14.76) into the normalization condition 1=-tQ(t^,t)dt^1=\int_{-\infty}^{t}Q(\hat{t},t)\,{\text{d}}\hat{t} we arrive at

1=-texp[-t^tρ(t|t^)dt]A(t^)dt^.1=\int_{-\infty}^{t}\exp\left[-\int_{\hat{t}}^{t}\rho(t^{\prime}|\hat{t})\,{% \text{d}}t^{\prime}\right]\,\,A(\hat{t})\,{\text{d}}\hat{t}\,. (14.79)

Both the population equation (14.78) and the normalization condition (14.79) play an important role in the general theory outlined in Section 14.1.

For a numerical implementation of population dynamics, it is more convenient to take a step back from the integral equations to the iterative updates that describe the survival Q(t^,t)Q(\hat{t},t) of the group of neurons that has fired the last time at time t^\hat{t}. In other words, efficient numerical implementation schemes work directly on the level of the density equations (14.74). A simple discretization scheme for numerical implementations has been discussed above in Section 14.1.5.

14.4.3 From refractory densities to membrane potential densities (*)

A B
Fig. 14.11: A. In a noise-free (nonlinear) integrate-and-fire neuron we define a refractory variable r=t-t^r=t-\hat{t} and write the trajectory as u(t,t^)u(t,\hat{t}). B. The membrane potential density p(u,t)p(u,t) at time tt is plotted to the left as a function of the membrane potential uu (vertical axis). The fraction of neurons p(u,t)Δup(u,t)\,\Delta u with membrane potentials around uu (gray shaded area) is proportional to the refractory density Q(t,t^)Δt^Q(t,\hat{t})\Delta\hat{t} (shaded area at the bottom).

In this section we want to show the formal relation between the dynamics of p(u,t)p(u,t) and the evolution of the refractory densities q(r,t)q(r,t). We focus on a population of nonlinear or leaky integrate-and-fire neurons with escape noise. For a known time-dependent input we can calculate the membrane potential u(t,t^)u(t,\hat{t}) where the notation with t^\hat{t} highlights that, in addition to external input, we also need to know the last firing time (and the reset value uru_{r}) in order to predict the voltage at time tt. We require that the input is constant or only weakly modulated so that the trajectory is increasing monotonously (u/t>0\partial u/\partial t>0).

To stay concrete, we focus on leaky integrate-and-fire neuron for which the membrane potential is

u(t,t^)=urexp(-t-t^τm)+Rτmt^texp(-t-tτm)I(t)dtu(t,\hat{t})=u_{r}\,\exp\left(-{t-\hat{t}\over\tau_{m}}\right)+{R\over\tau_{m}% }\int_{\hat{t}}^{t}\exp\left(-{t-t^{\prime}\over\tau_{m}}\right)\,I(t^{\prime}% )\,{\text{d}}t^{\prime} (14.80)

Knowledge of u(t,t^)u(t,\hat{t}) can be used to define a transformation from voltage to refractory variables: ur=t-t^u\longleftrightarrow r=t-\hat{t}; cf. Fig. 14.11A. It turns out that the final equations are even simpler if we take t^\hat{t} instead of rr as our new variable. We therefore consider the transformation ut^u\longrightarrow\hat{t}.

Before we start, we calculate the derivatives of Eq. (14.80). The derivative with respect to tt yields u/t=[-u+RI(t)]/τm{\partial u/\partial t}=[-u+R\,I(t)]/\tau_{m} as expected for integrate-and-fire neurons. According to our assumption, u/t>0{\partial u/\partial t}>0 or RI(t)>uRI(t)>u (for all neurons in the population). The derivative with respect to t^\hat{t} is

-ut^=RI(t)-urτmexp(-t-t^τm)=F(t,t^)>0,-{\partial u\over\partial\hat{t}}={R\,I(t)-u_{r}\over\tau_{m}}\,\exp\left(-{t-% \hat{t}\over\tau_{m}}\right)=F(t,\hat{t})>0\,, (14.81)

where the function FF is defined by Eq. (14.81).

The densities in the variable t^\hat{t} are denoted as Q(t^,t)Q(\hat{t},t). From Q(t^,t)dt^=p(u,t)duQ(\hat{t},t)\,{\text{d}}\hat{t}=p(u,t)\,{\text{d}}u we have

Q(t^,t)=p[u(t,t^),t]F(t,t^).Q(\hat{t},t)=p[u(t,\hat{t}),t]\,F(t,\hat{t})\,. (14.82)

We now want to show that the differential equation for the density Q(t^,t)Q(\hat{t},t) that we derived in (14.74),

tQ(t^,t)=-ρ(t|t^)Q(t^,t),fort^<t,{\partial\over\partial t}Q(\hat{t},t)=-\rho(t|\hat{t})\,Q(\hat{t},t)\,\,,\quad% {\rm for~{}}\hat{t}<t\,, (14.83)

is equivalent to the partial differential equation for the membrane potential densities. If we insert Eq. (14.82) into Eq. (14.83) we find

puutF+ptF+pFt=-ρpF.{\partial p\over\partial u}\,{\partial u\over\partial t}\,F+{\partial p\over% \partial t}\,F+p\,{\partial F\over\partial t}=-\rho\,p\,F\,. (14.84)

For the linear integrate-and-fire neuron we have F/t=-F/τm\partial F/\partial t=-F/\tau_{m}. Furthermore, according to our assumption of monotonously increasing trajectories, we have F0F\neq 0. Thus we can divide (14.84) by FF and rewrite Eq. (14.84) in the form

p(u,t)t=-u[-u+RI(t)τmp(u,t)]-f(u-ϑ)p(u,t), forur<u<ϑ,{\partial p(u,t)\over\partial t}=-{\partial\over\partial u}\left[{-u+R\,I(t)% \over\tau_{m}}\,p(u,t)\right]-f(u-\vartheta)\,p(u,t)\,,\quad{\rm for~{}}u_{r}<% u<\vartheta\,, (14.85)

where we have used the definition of the hazard via the escape function ρ(t|t^)=f[u(t,t^)-ϑ]\rho(t|\hat{t})=f[u(t,\hat{t})-\vartheta] and the definition of the reset potential ur=η0u_{r}=\eta_{0}. If we compare Eq. (14.85) with the Fokker-Planck equation (13.16), we see that the main difference is the treatment of the noise. For noise-free integrate-and-fire neurons (i.e., ρ(t|t^)=0\rho(t|\hat{t})=0 for uϑu\neq\vartheta) the equation (13.16) for the membrane potential densities is therefore equivalent to the density equation Q(t^,t)/t=0\partial Q(\hat{t},t)/\partial t=0; cf. Eq. (14.83).