8 Noisy Input Models: Barrage of Spike Arrivals

8.4 Diffusion limit and Fokker-Planck equation(*)

In this section we analyze the model of stochastic spike arrival defined in Eq. (8.20) and show how to map it to the diffusion model defined in Eq. (8.7) (190; 243; 88). Suppose that the neuron has fired its last spike at time t^\hat{t}. Immediately after the firing the membrane potential was reset to uru_{r}. Because of the stochastic spike arrival, we cannot predict the membrane potential for t>t^t>\hat{t}, but we can calculate its probability density, p(u,t)p(u,t). The evolution of the probability density as a function of time is described, in the diffusion limit, by the Fokker-Planck equation which we derive here in the context of a single neuron subject to noisy input. In part III (Chapter 13.1.1) the Fokker-Planck equation will be introduced systematically in the context of populations of neurons.

For the sake of simplicity, we set for the time being Iext=0I^{\rm ext}=0 in Eq. (8.20). The input spikes at synapse kk are generated by a Poisson process and arrive stochastically with rate νk(t)\nu_{k}(t). The probability that no spike arrives in a short time interval Δt\Delta t is therefore

Prob{nospikein[t,t+Δt]}=1-kνk(t)Δt.{\rm Prob}\big\{{\rm no~{}spike~{}in~{}}[t,t+\Delta t]\big\}=1-\sum_{k}\nu_{k}% (t)\,\Delta t\,. (8.36)

If no spike arrives in [t,t+Δt][t,t+\Delta t], the membrane potential changes from u(t)=uu(t)=u^{\prime} to u(t+Δt)=uexp(-Δt/τm)u(t+\Delta t)=u^{\prime}\,\exp(-\Delta t/\tau_{m}). On the other hand, if a spike arrives at synapse kk, the membrane potential changes from uu^{\prime} to uexp(-Δt/τm)+wku^{\prime}\,\exp(-\Delta t/\tau_{m})+w_{k}. Given a value of uu^{\prime} at time tt, the probability density of finding a membrane potential uu at time t+Δtt+\Delta t is therefore given by

Ptrans(u,t+Δt|u,t)\displaystyle P^{\rm trans}(u,t+\Delta t|u^{\prime},t) =\displaystyle= [1-Δtkνk(t)]δ(u-ue-Δt/τm)\displaystyle\big[1-\Delta t\sum_{k}\nu_{k}(t)\big]\,\delta\left(u-u^{\prime}% \,e^{-\Delta t/\tau_{m}}\right) (8.37)
+Δtkνk(t)δ(u-ue-Δt/τm-wk).\displaystyle+\Delta t\sum_{k}\nu_{k}(t)\delta\left(u-u^{\prime}\,e^{-\Delta t% /\tau_{m}}-w_{k}\right)\,.

We will refer to PtransP^{\rm trans} as the transition law. Since the membrane potential is given by the differential equation (8.20) with input spikes generated by a Poisson distribution, the evolution of the membrane potential is a Markov Process (i.e., a process without memory) and can be described by (529)

p(u,t+Δt)=Ptrans(u,t+Δt|u,t)p(u,t)du.p(u,t+\Delta t)=\int P^{\rm trans}(u,t+\Delta t|u^{\prime},t)\,p(u^{\prime},t)% \,{\text{d}}u^{\prime}\,. (8.38)

We insert Eq. (8.37) in Eq. (8.38). To perform the integration, we have to recall the rules for δ\delta functions, viz., δ(au)=a-1δ(u)\delta(a\,u)=a^{-1}\,\delta(u). The result of the integration is

p(u,t+Δt)\displaystyle p(u,t+\Delta t) =\displaystyle= [1-Δtkνk(t)]eΔt/τmp(eΔt/τmu,t)\displaystyle\big[1-\Delta t\sum_{k}\nu_{k}(t)\big]\,e^{\Delta t/\tau_{m}}\,p% \left(e^{\Delta t/\tau_{m}}\,u,t\right) (8.39)
+Δtkνk(t)eΔt/τmp(eΔt/τmu-wk,t).\displaystyle+\Delta t\sum_{k}\nu_{k}(t)\,e^{\Delta t/\tau_{m}}\,p\left(e^{% \Delta t/\tau_{m}}\,u-w_{k},t\right)\,.

Since Δt\Delta t is assumed to be small, we expand Eq. (8.39) about Δt=0\Delta t=0 and find to first order in Δt\Delta t

p(u,t+Δt)-p(u,t)Δt\displaystyle{p(u,t+\Delta t)-p(u,t)\over\Delta t} =\displaystyle= 1τmp(u,t)+1τmuup(u,t)\displaystyle{1\over\tau_{m}}\,p(u,t)+{1\over\tau_{m}}\,u\,{\partial\over% \partial u}p(u,t) (8.40)

For Δt0\Delta t\to 0, the left-hand side of Eq. (8.40) turns into a partial derivative p(u,t)/t\partial p(u,t)/\partial t. Furthermore, if the jump amplitudes wkw_{k} are small, we can expand the right-hand side of Eq. (8.40) with respect to uu about p(u,t)p(u,t):

τmtp(u,t)\displaystyle\tau_{m}{\partial\over\partial t}p(u,t) =\displaystyle= -u[-u+τmkνk(t)wk]p(u,t)\displaystyle-{\partial\over\partial u}\left[-u+\tau_{m}\sum_{k}\nu_{k}(t)\,w_% {k}\right]\,p(u,t) (8.41)
+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),

where we have neglected terms of order wk3w_{k}^{3} and higher. The expansion in wkw_{k} is called the Kramers-Moyal expansion. Eq. (8.41) is an example of a Fokker-Planck equation (529), i.e., a partial differential equation that describes the temporal evolution of a probability distribution. The right-hand side of Eq. (8.41) has a clear interpretation: The first term in rectangular brackets describes the systematic drift of the membrane potential due to leakage (-u\propto-u) and mean background input (kνk(t)wk\propto\sum_{k}\nu_{k}(t)\,w_{k}). The second term in rectangular brackets corresponds to a ‘diffusion constant’ and accounts for the fluctuations of the membrane potential. The Fokker-Planck Eq. (8.41) is equivalent to the Langevin equation (8.7) with RI(t)=τmkνk(t)wkR\,I(t)=\tau_{m}\sum_{k}\nu_{k}(t)\,w_{k} and time-dependent noise amplitude

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

The specific process generated by the Langevin-equation (8.7) with constant noise amplitude σ\sigma is called the Ornstein-Uhlenbeck process (527), but Eq. (8.42) indicates that, in the context of neuroscience, the effective noise amplitude generated by stochastic spike arrival is in general time-dependent. We will return to the Fokker-Planck equation in Ch. 13.

For the transition from Eq. (8.40) to (8.41) we have suppressed higher-order terms in the expansion. The missing terms are

n=3(-1)nn!An(t)nunp(u,t)\sum_{n=3}^{\infty}{(-1)^{n}\over n!}\,A_{n}(t){\partial^{n}\over\partial u^{n% }}p(u,t) (8.43)

with An=τmkνk(t)wknA_{n}=\tau_{m}\sum_{k}\nu_{k}(t)\,w_{k}^{n}. What are the conditions that these terms vanish? As in the example of Figs. 8.6A and B, we consider a sequence of models where the size of the weights wkw_{k} decreases so that An0A_{n}\to 0 for n3n\geq 3 while the mean kνk(t)wk\sum_{k}\nu_{k}(t)\,w_{k} and the second moment kνk(t)wk2\sum_{k}\nu_{k}(t)\,w_{k}^{2} remain constant. It turns out, that, given both excitatory and inhibitory input, it is always possible to find an appropriate sequence of models (286; 287). For wk0w_{k}\to 0, the diffusion limit is attained and Eq. (8.41) is exact. For excitatory input alone, however, such a sequence of models does not exist (406).

8.4.1 Threshold and firing

The Fokker-Planck equation (8.41) and the Langevin equation (8.7) are equivalent descriptions of drift and diffusion of the membrane potential. Neither of these describe spike firing. To turn the Langevin equation (8.7) into a sensible neuron model, we have to incorporate a threshold condition. In the Fokker-Planck equation (8.41), the firing threshold is incorporated as a boundary condition

p(ϑ,t)0forallt.p(\vartheta,t)\equiv 0\,\quad{\rm for~{}all~{}}t\,. (8.44)

The boundary condition reflects the fact that, under a white noise model, in each short interval Δt\Delta t many excitatory and inhibitory input spikes arrive which each cause a tiny jump of size ±δ\pm\delta of the membrane potential. Any finite density p(u,t)p(u,t) at a value ϑ-δ<uϑ\vartheta-\delta<u\leq\vartheta would be rapidly removed, because one of the many excitatory spikes which arrive at each moment would push the membrane potential above threshold. The white noise limit corresponds to infinite spike arrival rate and jump size δ0\delta\to 0, as discussed above. As a consequence there are, in each short interval Δt\Delta t infinitely many ’attempts’ to push the membrane potential above threshold. Hence the density at u=ϑu=\vartheta must vanish. The above argument also shows that for colored noise the density at threshold is finite, because the the effective frequency of ’attempts’ is limited by the cut-off frequency of the noise.

Before we continue the discussion of the diffusion model in the presence of a threshold, let us study the solution of Eq. (8.41) without threshold.

Example: Free distribution

The solution of the Fokker-Planck equation (8.41) with initial condition p(u,t^)=δ(u-ur)p(u,\hat{t})=\delta(u-u_{r}) is a Gaussian with mean u0(t)u_{0}(t) and variance Δu2(t)\langle\Delta u^{2}(t)\rangle, i.e.,

p(u,t)=12πΔu2(t)exp{-[u(t|t^)-u0(t)]22Δu2(t)}p(u,t)={1\over\sqrt{2\pi\,\langle\Delta u^{2}(t)\rangle}}\,\exp\left\{-{[u(t|% \hat{t})-u_{0}(t)]^{2}\over 2\,\langle\Delta u^{2}(t)\rangle}\right\} (8.45)

as can be verified by inserting Eq. (8.45) into (8.41). In particular, the stationary distribution that is approached in the limit of tt\to\infty for constant input I0I_{0} is

p(u,)=1π1σexp{[u-RI0]2σ2},p(u,\infty)={1\over\sqrt{\pi}}{1\over\sigma}\exp\left\{{[u-R\,I_{0}]^{2}\over% \sigma^{2}}\right\}\,, (8.46)

which describes a Gaussian distribution with mean u=RI0u_{\infty}=R\,I_{0} and variance σ/2\sigma/\sqrt{2}.

8.4.2 Interval distribution for the diffusive noise model

Let us consider a leaky integrate-and-fire neuron that starts at time t^\hat{t} with a membrane potential uru_{r} and is driven for t>t^t>\hat{t} by a known input I(t)I(t). Because of the diffusive noise generated by stochastic spike arrival, we cannot predict the exact value of the neuronal membrane potential u(t)u(t) at a later time t>t^t>\hat{t}, but only the probability that the membrane potential is in a certain interval [u0,u1][u_{0},u_{1}]. Specifically, we have

Prob{u0<u(t)<u1| u(t^)=ur}=u0u1p(u,t)du{\rm Prob}\left\{u_{0}\!<\!u(t)\!<\!u_{1}\,|\,u(\hat{t})\!=\!u_{r}\right\}\,=% \,\int_{u_{0}}^{u_{1}}p(u,t)\,{\text{d}}u (8.47)

where p(u,t)p(u,t) is the probability density of the membrane potential at time tt. In the diffusion limit, p(u,t)p(u,t) can be found by solution of the Fokker-Planck equation Eq. (8.41) with initial condition p(u,t^)=δ(u-ur)p(u,\hat{t})=\delta(u-u_{r}) and boundary condition p(ϑ,t)=0p(\vartheta,t)=0.

The boundary is absorbing. In other words, in a simulation of a single realization of the stochastic process, the simulation is stopped when the trajectory passes the threshold for the first time. To be concrete, imagine that we run 100 simulation trials, i.e., 100 realization of a leaky integrate-and-fire model with diffusive noise. Each trial starts at the same value u(t^)=uru(\hat{t})=u_{r} and uses the same input current I(t)I(t^{\prime}) for t>t^t^{\prime}>\hat{t}. Some of the trials will exhibit trajectories that reach the threshold at some point t<tt^{\prime}<t. Others stay below threshold for the whole period t^<t<t\hat{t}<t^{\prime}<t. The expected fraction of simulations that have not yet reached the threshold and therefore still ’survive’ up to time tt is given by the survivor function,

SI(t|t^)=-ϑp(u,t)du.S_{I}(t|\hat{t})=\int_{-\infty}^{\vartheta}p(u,t)\,{\text{d}}u\,. (8.48)

In other words, the survivor function in the diffusive noise model is equal to the probability that the membrane potential has not yet reached the threshold between t^\hat{t} and tt.

In view of Eq. (7.24), the input-dependent interval distribution is therefore

PI(t|t^)=-ddt-ϑp(u,t)du.P_{I}(t|\hat{t})=-\frac{{\text{d}}}{{\text{d}}t}\int_{-\infty}^{\vartheta}p(u,% t)\,{\text{d}}u\,. (8.49)

We recall that PI(t|t^)ΔtP_{I}(t|\hat{t})\,\Delta t for Δt0\Delta t\to 0 is the probability that a neuron fires its next spike between tt and t+Δtt+\Delta t given a spike at t^\hat{t} and input II. In the context of noisy integrate-and-fire neurons PI(t|t^)P_{I}(t|\hat{t}) is called the distribution of ‘first passage times’. The name is motivated by the fact, that firing occurs when the membrane potential crosses ϑ\vartheta for the first time. Unfortunately, no general solution is known for the first passage time problem of the Ornstein-Uhlenbeck process. For constant input I(t)=I0I(t)=I_{0}, however, it is at least possible to give a moment expansion of the first passage time distribution. In particular, the mean of the first passage time can be calculated in closed form.

Fig. 8.9: Without a threshold, several trajectories can reach at time tt the same value u=ϑu=\vartheta from above or below.

Example: Numerical evaluation of PI(t|t^)P_{I}(t|\hat{t})

We have seen that, in the absence of a threshold, the Fokker-Planck Equation (8.41) can be solved; cf. Eq. (8.45). The transition probability from an arbitrary starting value uu^{\prime} at time tt^{\prime} to a new value uu at time tt is

Ptrans(u,t|u,t)=12πΔu2(t)exp{-[u-u0(t)]22Δu2(t)}P^{\rm trans}(u,t|u^{\prime},t^{\prime})={1\over\sqrt{2\pi\,\langle\Delta u^{2% }(t)\rangle}}\,\exp\left\{-{[u-u_{0}(t)]^{2}\over 2\,\langle\Delta u^{2}(t)% \rangle}\right\} (8.50)


u0(t)\displaystyle u_{0}(t) =\displaystyle= ue-(t-t)/τm+0t-te-s/τmI(t-s)   d   s\displaystyle u^{\prime}\,e^{-(t-t^{\prime})/\tau_{m}}+\int_{0}^{t-t^{\prime}}% e^{-s^{\prime}/\tau_{m}}\,I(t-s^{\prime})\,{\text{d}}s (8.51)
Δu2(t)\displaystyle\langle\Delta u^{2}(t)\rangle =\displaystyle= σ22[1-e-2(t-s)/τm].\displaystyle{\sigma^{2}\over 2}\left[1-e^{-2\,(t-s)/\tau_{m}}\right]\,. (8.52)

A method due to Schrödinger uses the solution of the unbounded problem in order to calculate the input-dependent interval distribution PI(t|t^)P_{I}(t|\hat{t}) of the diffusion model with threshold (459; 404; 83). The idea of the solution method is illustrated in Fig. 8.9. Because of the Markov property, the probability density to cross the threshold (not necessarily for the first time) at a time tt, is equal to the probability to cross it for the first time at t<tt^{\prime}<t and to return back to ϑ\vartheta at time tt, that is,

Ptrans(ϑ,t|ur,t^)=t^tPI(t|t^)Ptrans(ϑ,t|ϑ,t)   d   t.P^{\rm trans}(\vartheta,t|u_{r},\hat{t})=\int_{\hat{t}}^{t}P_{I}(t^{\prime}|% \hat{t})\,P^{\rm trans}(\vartheta,t|\vartheta,t^{\prime})\,{\text{d}}t^{\prime% }\,. (8.53)

This integral equation can be solved numerically for the distribution PI(t|t^)P_{I}(t^{\prime}|\hat{t}) for arbitrary input current I(t)I(t) (405). An example is shown in Fig. 8.10. The probability to emit a spike is high whenever the noise-free trajectory is close to the firing threshold. Very long intervals are unlikely, if the noise-free membrane potential was already several times close to the threshold before, so that the neuron has had ample opportunity to fire earlier.

Fig. 8.10: A time-dependent input current I(t)I(t) generates a noise-free membrane potential u0(t)u_{0}(t) shown in the lower part of the figure. In the presence of diffusive noise, spikes can be triggered although the reference trajectory stays below the threshold (dashed line). This gives rise to an input-dependent interval distribution PI(t|0)P_{I}(t|0) shown in the upper panel. Taken from (403).

8.4.3 Mean interval and mean firing rate (diffusive noise)

For constant input I0I_{0} the mean interspike interval is s=0sPI0(s|0)ds=0sP0(s)ds\langle s\rangle=\int_{0}^{\infty}s\,P_{I_{0}}(s|0){\text{d}}s=\int_{0}^{% \infty}s\,P_{0}(s){\text{d}}s; cf. Eq. (7.13). For the diffusion model Eq. (8.7) with threshold ϑ\vartheta, reset potential uru_{r}, and membrane time constant τm\tau_{m}, the mean interval is

s=τmπur-h0σϑ-h0σduexp(u2)[1+erf(u)],\langle s\rangle=\tau_{m}\,\sqrt{\pi}\int_{{u_{r}-h_{0}\over\sigma}}^{{% \vartheta-h_{0}\over\sigma}}{\text{d}}u\,\exp\left(u^{2}\right)\,\left[1+{\rm erf% }(u)\right]\,, (8.54)

where h0=RI0h_{0}=R\,I_{0} is the input potential caused by the constant current I0I_{0} (476; 243). Here ’erf’ denotes the error function erf(x)=22pi0xexp(-u2)du{\rm erf}(x)=\frac{2}{\sqrt{2pi}}\int_{0}^{x}\exp(-u^{2})du. This expression, sometimes called the Siegert-formula, can be derived by several methods; for reviews see, e.g., (529). The inverse of the mean interval is the mean firing rate. Hence, Eq. (8.54) enables us the express the mean firing rate of a leaky integrate-and-fire model with diffusive noise as a function of a (constant) input I0I_{0} (Fig.  8.11). We will derive Eq. (8.54) in Chapter 13 in the context of populations of spiking neurons.

Fig. 8.11: Mean firing rate of a leaky integrate-and-fire model as a function of constant input evaluated for different levels of diffusive noise, using the Siegert formula, Eq. (8.54). From top to bottom: σ=1.0\sigma=1.0, σ=0.5\sigma=0.5, σ=0.2\sigma=0.2 (solid line), σ=0.1\sigma=0.1, σ=0.0\sigma=0.0.