13 Continuity Equation and the Fokker-Planck Approach

13.1 Continuity equation

In a population of neurons, each neuron may be in a different internal state. In this section we derive partial differential equations that describe how the distribution of internal states evolves as a function of time. We start in Section 13.1.1 with a population of integrate-and-fire neurons. Since the state of an integrate-and-fire neuron is characterized by its membrane potential, we describe the dynamics of the population as the evolution of membrane potential densities.

Population activity, A(t)A(t), was introduced in Chapter 7 as the fraction of neurons that fire at time tt in a finite population. In this chapter and the next, the population activity is the expected population activity A(t)A(t)A(t)\equiv\langle A(t)\rangle. Since A(t) is self-averaging, it can also be said that we consider the limit of large populations. The finite-size effects will be discussed in Section 14.6.

The formulation of the dynamics of a population of integrate-and-fire neurons on the level of membrane potential densities has been developed by Abbott and van Vreeswijk (4), Brunel and Hakim (78), Fusi and Mattia (163), Nykamp and Tranchina (367), Brunel (79) and Omurtag et al. (371). The closely related formulation in terms of refractory densities has been studied by Wilson and Cowan (552), Gerstner and van Hemmen (180), Bauer and Pawelzik (41), and Gerstner (183). Generalized density equations have been discussed by Knight (265) and Fourcaud and Brunel (154).

13.1.1 Distribution of membrane potentials

We study a homogeneous population of integrate-and-fire neurons. The internal state of a neuron ii is determined by its membrane potential which changes according to

τmddtui=f(ui)+RIi(t)forui<θreset.\tau_{m}{{\text{d}}\over{\text{d}}t}u_{i}=f(u_{i})+R\,I_{i}(t)\quad{\rm for}% \quad u_{i}<\theta_{\rm reset}\,. (13.1)

Here RR is the input resistance, τm=RC\tau_{m}=R\,C the membrane time constant, and Ii(t)I_{i}(t) the total input (external driving current and synaptic input). If uiθresetu_{i}\geq\theta_{\rm reset} the membrane potential is reset to ui=ur<θresetu_{i}=u_{r}<\theta_{\rm reset}. Here f(ui)f(u_{i}) is an arbitrary function of uu; cf. Eq. (5.2) in Chapter 5. For f(ui)=-(ui-urest)f(u_{i})=-(u_{i}-u_{\rm rest}) the equation reduces to the standard leaky integrate-and-fire model. For the moment we keep the treatment general and restrict it to the leaky integrate-and-fire model only later.

In a population of NN integrate-and-fire neurons, we may ask how many of the neurons have at time tt a given membrane potential. For NN\to\infty the fraction of neurons ii with membrane potential u0<ui(t)u0+Δuu_{0}<u_{i}(t)\leq u_{0}+\Delta u is

limN{neuronswithu0<ui(t)u0+ΔuN}=u0u0+Δup(u,t)du,\lim_{N\to\infty}\left\{{{\rm neurons~{}with~{}}u_{0}<u_{i}(t)\leq u_{0}+% \Delta u\over N}\right\}=\int_{u_{0}}^{u_{0}+\Delta u}p(u,t)\,{\text{d}}u\,, (13.2)

where p(u,t)p(u,t) is the membrane potential density; cf. Chapter 8. The integral over this density remains constant over time, i.e.,

-θresetp(u,t)du=1 .\int_{-\infty}^{\theta_{\rm reset}}p(u,t)\,{\text{d}}u=1\,. (13.3)

The normalization to unity expresses the fact that all neurons have a membrane potential below or equal to threshold.

The aim of this section is to describe the evolution of the density p(u,t)p(u,t) as a function of time. As we will see, the equation that describes the dynamics of p(u,t)p(u,t) is nearly identical to that of a single integrate-and-fire neuron with diffusive noise; cf. Eqs. (8.40) and (8.41) in Chapter 8.

13.1.2 Flux and continuity equation

Let us consider the portion of neurons with a membrane potential between u0u_{0} and u1u_{1},

n(u0;u1)N=u0u1p(u,t)du.{n(u_{0};u_{1})\over N}=\int_{u_{0}}^{u_{1}}p(u^{\prime},t)\,{\text{d}}u^{% \prime}\,. (13.4)

The fraction of neurons with u0<u<u1u_{0}<u<u_{1} increases if neurons enter from below through the boundary u0u_{0} or from above through the boundary u1u_{1}; cf. Fig. 13.1. Since there are many neurons in the population, we expect that in each short time interval Δt\Delta t, many trajectories cross one of the boundaries. The flux J(u,t)J(u,t) is the net fraction of trajectories per unit time that crosses the value uu. A positive flux J(u,t)>0J(u,t)>0 is defined as a flux toward increasing values of uu. In other words, in a finite population of NN neurons, the quantity NJ(u0,t)ΔtN\,J(u_{0},t)\,\Delta t describes the number of trajectories that cross in the interval Δt\Delta t the boundary u0u_{0} from below, minus the number of trajectories crossing u0u_{0} from above. Note that u0u_{0} here is an arbitrary value that we chose as the lower bound of the integral in Eq. (13.4) and as such has no physical meaning for the neuron model.

Fig. 13.1: The number of trajectories in the interval [u0,u1][u_{0},u_{1}] changes if one of the trajectories crosses the boundary u0u_{0} or u1u_{1}. For a large number of neurons this fact is described by the continuity equation; cf. Eq. (13.6). Schematic figure where only three trajectories are shown.

Since trajectories cannot simply end, a change in the number n(u0;u1)n(u_{0};u_{1}) of trajectories in the interval u0<u<u1u_{0}<u<u_{1} can be traced back to the flux of trajectories in and out of that interval. We therefore have the conservation law

tu0u1p(u,t)du=J(u0,t)-J(u1,t)    .{\partial\over\partial t}\int_{u_{0}}^{u_{1}}p(u^{\prime},t)\,{\text{d}}u^{% \prime}=J(u_{0},t)-J(u_{1},t)\qquad\,. (13.5)

Taking the derivative with respect to the upper boundary u1u_{1} and changing the name of the variable from u1u_{1} to uu yields the continuity equation,

tp(u,t)=-uJ(u,t)  foruuranduθreset,\qquad{\partial\over\partial t}p(u,t)=-{\partial\over\partial u}J(u,t)\quad{% \rm~{}for~{}}u\neq u_{r}{\rm~{}and~{}}u\neq\theta_{\rm reset}\,, (13.6)

which expresses the conservation of the number of trajectories. In integrate-and-fire models, however, there are two special voltage values, uru_{r} and θreset\theta_{\rm reset}, where the number of trajectories is not conserved, because of the fire-and-reset mechanism.

Since neurons that have fired start a new trajectory at uru_{r}, we have a ‘source of new trajectories’ at u=uru=u_{r}, i.e., new trajectories appear in the interval [ur-ϵ,ur+ϵ][u_{r}-\epsilon,u_{r}+\epsilon] that have not entered the interval through one of the borders. Adding a term A(t)δ(u-ur)A(t)\,\delta(u-u_{r}) on the right-hand side of (13.6) accounts for this source of trajectories. The trajectories that appear at uru_{r} disappear at θreset\theta_{\rm reset}, so that we have

tp(u,t)=-uJ(u,t)+A(t)δ(u-ur)-A(t)δ(u-θreset).\qquad{\partial\over\partial t}p(u,t)=-{\partial\over\partial u}J(u,t)+A(t)\,% \delta(u-u_{r})-A(t)\,\delta(u-\theta_{\rm reset})\,. (13.7)

The density p(u,t)p(u,t) vanishes for all values u>θresetu>\theta_{\rm reset}.

The population activity A(t)A(t) is, by definition, the fraction of neurons that fire, i.e., those that pass through the threshold. Therefore we find

A(t)=J(θreset,t).A(t)=J(\theta_{\rm reset},t)\,. (13.8)

Equations (13.7) and (13.8) describe the evolution of the the membrane potential densities and the resulting population activity as a function of time. We now specify the neuron model so as to have an explicit expression for the flux.