13 Continuity Equation and the Fokker-Planck Approach

13.3 Fokker-Planck equation

The equation for the membrane potential density in a population of integrate-and-fire neurons (see Eq. (13.14) in the previous section) can, in the limit of small jump amplitudes wkw_{k}, be approximated by a diffusion equation. To show this, we expand the right-hand side of Eq. (13.14) into a Taylor series up to second order in wkw_{k}. The result is the Fokker-Planck equation,

τmtp(u,t)\displaystyle\tau_{m}{\partial\over\partial t}p(u,t) =\displaystyle= -u{[f(u)+RIext(t)+τmkνk(t)wk]p(u,t)}\displaystyle-{\partial\over\partial u}\left\{\left[f(u)+R\,I^{\rm ext}(t)+% \tau_{m}\sum_{k}\nu_{k}(t)\,w_{k}\right]\,p(u,t)\right\} (13.16)
+12[τmkνk(t)wk2]2u2p(u,t)\displaystyle+{1\over 2}\,\left[\tau_{m}\sum_{k}\nu_{k}(t)\,w_{k}^{2}\right]\,% {\partial^{2}\over\partial u^{2}}p(u,t)
+τmA(t)δ(u-ur)-τmA(t)δ(u-θreset)+𝒪(ωk3).\displaystyle+\tau_{m}\,A(t)\,\delta(u-u_{r})-\tau_{m}\,A(t)\,\delta(u-\theta_% {\rm reset})+{\mathcal{O}}(\omega_{k}^{3})\,.

The term with the second derivative describes a ‘diffusion’ in terms of the membrane potential.

It is convenient to define the total ‘drive’ in voltage units as

μ(t)=RIext(t)+τmkνk(t)wk\mu(t)=R\,I^{\rm ext}(t)+\tau_{m}\sum_{k}\nu_{k}(t)\,w_{k} (13.17)

and the amount of diffusive noise (again in voltage units) as

σ2(t)=τmkνk(t)wk2.\sigma^{2}(t)=\tau_{m}\sum_{k}\nu_{k}(t)\,w_{k}^{2}\,. (13.18)

The firing threshold acts as an absorbing boundary so that the density at threshold vanishes:

p(θreset,t)=0 .p(\theta_{\rm reset},t)=0\,. (13.19)

The boundary condition p(θreset,t)=0p(\theta_{\rm reset},t)=0 arises from the mathematical limiting process: We consider that all weights go to zero, wk0w_{k}\to 0, while the total drive μ\mu and the amount of diffusive noise remains constant. This is only possible if we have at least two different types of synapse (excitatory and inhibitory) and we increase the rate of both while decreasing the weights of synapses. In other words, the spike arrival rates at the synapses go to infinity so that the ‘shot noise’ generated by spike arrivals turns into a white Gaussian distributed noise. Any neuron with a membrane potential just below threshold would then immediately fire because of the noise. Therefore the density at threshold has to be zero. This argument is no longer valid, if spike arrival rates are finite, or if synapses react slowly so that the noise is colored (cf. Section 13.6.4 below).

In order to calculate the flux through the reset threshold we expand Eq. (13.13) in wkw_{k} about u=θresetu=\theta_{\rm reset} and obtain

A(t)=-σ2(t)2τmp(u,t)u|u=θreset,A(t)=-\left.{\sigma^{2}(t)\over 2\tau_{m}}\,{\partial p(u,t)\over\partial u}% \right|_{u=\theta_{\rm reset}}\,, (13.20)

Eqs. (13.16) – (13.20) together with the normalization (13.3) define the dynamics of a homogeneous population of integrate-and-fire neurons with ‘diffusive’ noise. For a more detailed discussion of the diffusion limit see Chapter 8, in particular Eqs. (8.41) and (8.43).

Example: Flux in the diffusion limit

If we compare the continuity equation (13.7) with the explicit form of the Fokker-Planck equation (13.16), we can identify the flux caused by stochastic spike arrival and external current in the diffusion limit

Jdiff(u,t)=f(u)+μ(t)τmp(u,t)-12σ2(t)τmup(u,t).\,J^{\rm diff}(u,t)={f(u)+\mu(t)\over\tau_{m}}\,p(u,t)-{1\over 2}\,{\sigma^{2}% (t)\over\tau_{m}}{\partial\over\partial u}p(u,t). (13.21)

We emphasize that stochastic spike arrival contributes to the mean drive μ(t)=RIext(t)+τmkνk(t)wk\mu(t)=R\,I^{\rm ext}(t)+\tau_{m}\sum_{k}\nu_{k}(t)\,w_{k} as well as to the diffusive noise σ2(t)=τmkνk(t)wk2\sigma^{2}(t)=\tau_{m}\sum_{k}\nu_{k}(t)\,w_{k}^{2}.

13.3.1 Stationary solution for leaky integrate-and-fire neurons(*)

The stationary distribution p(u)p(u) of the membrane potential is of particular interest, since it is experimentally accessible (86; 126; 108). We now derive the stationary solution p(u,t)p(u)p(u,t)\equiv p(u) of the Fokker-Planck equation (13.16) for a population of leaky integrate-and-fire neurons. Neurons are described by Eq. (13.1) with f(u)=-uf(u)=-u, i.e., the voltage scale is shifted so that the equilibrium potential is at zero. For leaky integrate-and-fire neurons, the reset threshold is the same as the rheobase firing threshold and will be denoted by ϑ=θreset=ϑrh\vartheta=\theta_{\rm reset}=\vartheta_{rh}.

We assume that the total input h0=RIext+τmkνkwkh_{0}=R\,I^{\rm ext}+\tau_{m}\sum_{k}\nu_{k}\,w_{k} is constant. In the stationary state, the temporal derivative on the left-hand-side of Eq. (13.16) vanishes. The terms on the right-hand side can be transformed such that the stationary Fokker-Planck equation reads (for u<ϑu<\vartheta)

0=-uJ(u)+A0δ(u-ur),0=-{\partial\over\partial u}J(u)+A_{0}\,\delta(u-u_{r})\,, (13.22)

where A0A_{0} is the population activity (or mean firing rate) in the stationary state and

J(u)=-u+h0τmp(u)-12σ2τmup(u)J(u)={-u+h_{0}\over\tau_{m}}\,p(u)-{1\over 2}{\sigma^{2}\over\tau_{m}}{% \partial\over\partial u}p(u) (13.23)

is the total flux; cf. Eq. (13.6). The meaning of Eq. (13.22) is that the flux is constant except at u=uru=u_{r} where it jumps by an amount A0A_{0}. Similarly, the boundary condition p(ϑ,t)=0p(\vartheta,t)=0 implies a second discontinuity of the flux at u=ϑu=\vartheta.

With the results from Chapter 8 in mind, we expect that the stationary solution approaches a Gaussian distribution for u-u\to-\infty. In fact, we can check easily that for any constant c1c_{1}

p(u)=c1σexp[-(u-h0)2σ2]  foruurp(u)={c_{1}\over\sigma}\exp\left[-{(u-h_{0})^{2}\over\sigma^{2}}\right]\quad{% \rm~{}for~{}}u\leq u_{r} (13.24)

is a solution of Eq. (13.22) with flux J(u)=0J(u)=0. However, for u>uru>u_{r} a simple Gaussian distribution cannot be a solution since it does not respect the boundary condition p(ϑ)=0p(\vartheta)=0. Nevertheless, we can make an educated guess and try a modified Gaussian, (188; 78)

p(u)=c2σ2exp[-(u-h0)2σ2]uϑexp[(x-h0)2σ2]dx  forur<uϑ,p(u)={c_{2}\over\sigma^{2}}\exp\left[-{(u-h_{0})^{2}\over\sigma^{2}}\right]\,% \cdot\,\int_{u}^{\vartheta}\exp\left[{(x-h_{0})^{2}\over\sigma^{2}}\right]\,{% \text{d}}x\quad{\rm~{}for~{}}u_{r}<u\leq\vartheta\,, (13.25)

with some constant c2c_{2}. We have written the above expression as a product of two terms. The first factor on the right-hand side is a standard Gaussian while the second factor guarantees that p(u)0p(u)\to 0 for uϑu\to\vartheta. If we insert Eq. (13.25) in (13.22) we can check that it is indeed a solution. The constant c2c_{2} is proportional to the flux,

c2=2τmJ(u)  forur<uϑ.c_{2}=2\,\tau_{m}\,J(u)\quad{\rm~{}for~{}}u_{r}<u\leq\vartheta\,. (13.26)

The solution defined by Eqs. (13.24) and (13.25) must be continuous at u=uru=u_{r}. Hence

c1=c2σurϑexp[(x-h0)2σ2]dx.c_{1}={c_{2}\over\sigma}\int_{u_{r}}^{\vartheta}\exp\left[{(x-h_{0})^{2}\over% \sigma^{2}}\right]\,{\text{d}}x\,. (13.27)

Finally, the constant c2c_{2} is determined by the normalization condition (13.3). We use Eqs. (13.24), (13.25), and (13.27) in (13.3) and find

1c2=-ururϑf(x,u)dxdu+urϑuϑf(x,u)dxdu=urϑ-xf(x,u)dudx,{1\over c_{2}}=\int_{-\infty}^{u_{r}}\int_{u_{r}}^{\vartheta}f(x,u)\,{\text{d}% }x\,{\text{d}}u+\int_{u_{r}}^{\vartheta}\int_{u}^{\vartheta}f(x,u)\,{\text{d}}% x\,{\text{d}}u=\int_{u_{r}}^{\vartheta}\int_{-\infty}^{x}f(x,u)\,{\text{d}}u\,% {\text{d}}x\,, (13.28)


f(x,u)=1σ2exp[-(u-h0)2σ2]exp[(x-h0)2σ2].f(x,u)={1\over\sigma^{2}}\exp\left[-{(u-h_{0})^{2}\over\sigma^{2}}\right]\,% \exp\left[{(x-h_{0})^{2}\over\sigma^{2}}\right]\,. (13.29)

Figure 13.5B shows the resulting stationary density p(u)p(u) for different noise amplitudes.

The activity A0A_{0} is identical to the flux J(u)J(u) between uru_{r} and ϑ\vartheta and therefore proportional to the constant c2c_{2}; cf. Eq. (13.26). If we express the integral over uu in Eq. (13.28) in terms of the error function, erf(x)=2π0xexp(-u2)du{\rm erf}(x)={2\over\sqrt{\pi}}\int_{0}^{x}\exp(-u^{2})\,{\text{d}}u, we obtain

A0-1=τmπur-h0σϑ-h0σexp(x2)[1+erf(x)]dx,A_{0}^{-1}=\tau_{m}\,\sqrt{\pi}\,\int_{{u_{r}-h_{0}\over\sigma}}^{{\vartheta-h% _{0}\over\sigma}}\,\exp\left({x}^{2}\right)\,\left[1+{\rm erf}(x)\right]\,{% \text{d}}x\,, (13.30)

which is identical to the Siegert formula (476) of the single neuron firing rate [cf. Eqs. (8.54) or (12.26)] that we used previously in Chapters 8 and 12. An alternative formulation (see Exercises) of the remaining integral in Eq. (13.30) is also possible (78).

Fig. 13.5: A. Membrane potential trajectories of 5 neurons (R=1R=1 and τm=10\tau_{m}=10\,ms) driven by a constant background current I0=0.8I_{0}=0.8 and stochastic background input with ν+=ν-=0.8\nu_{+}=\nu_{-}=0.8\,kHz and w±=±0.05w_{\pm}=\pm 0.05. These parameters correspond to h0=0.8h_{0}=0.8 and σ=0.2\sigma=0.2 in the diffusive noise model. B. Stationary membrane potential distribution in the diffusion limit for σ=0.2\sigma=0.2 (solid line), σ=0.1\sigma=0.1 (short-dashed line), and σ=0.5\sigma=0.5 (long-dashed line). (Threshold ϑ=1\vartheta=1). C. Mean activity of a population of integrate-and-fire neurons with diffusive noise as a function of h0h_{0} for four different noise levels, viz. (from top to bottom) σ=1.0,σ=0.5,σ=0.2\sigma=1.0,\sigma=0.5,\sigma=0.2 (solid line), σ=0.1,σ=0.0\sigma=0.1,\sigma=0.0.