13 Continuity Equation and the Fokker-Planck Approach


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), Omurtag et al. (371), and Knight (265), but the Fokker-Planck equation has been used much earlier for the probabilistic description of a single neuron driven by stochastic spike arrivals (243; 88; 430). The classic application of the Fokker-Planck approach to a network of excitatory and inhibitory leaky integrate-and-fire neurons is Brunel (79). For the general theory of Fokker-Planck equations see Risken (440).

Efficient numerical solutions of the Fokker-Planck equation, both for stationary input and periodic input have been developed by M.J.E. Richardson (2007,2009). These methods can be used to find activity in networks of nonlinear integrate-and-fire models, but also generalized neuron models with slow spike-triggered currents or conductances (435). For the treatment of colored noise, see Fourcaud and Brunel (154).


  1. 1.

    Diffusion limit in the Quadratic integrate-and-fire model . Consider a population of quadratic integrate-and-fire models. Assume that spikes from external sources arrive at excitatory and inhibitory synapses stochastically, and independently for different neurons, at a rate νE(t)\nu_{E}(t) and νI(t)\nu_{I}(t), respectively.

    (i) Write down the membrane potential density equations assuming that each spike causes a voltage jump by an amount ±Δu\pm\Delta u.

    (ii) Take the diffusion limit so as to arrive at a Fokker-Planck equation.

  2. 2.

    Voltage distribution of the Quadratic integrate-and-fire model . Find the stationary solution of the membrane potential distribution for the quadratic integrate-and-fire model with white diffusive noise.

  3. 3.

    Non-leaky integrate-and-fire model Consider a non-leaky integrate-and fire model subject to stochastic spike arrival

       𝑑   u   𝑑   t=1CI(t)=qCfδ(t-t(f)){{\text{d}}u\over{\text{d}}t}={1\over C}I(t)={q\over C}\sum_{f}\delta(t-t^{(f)}) (13.76)

    where qq is the charge that each spike puts on the membrane and the spike arrival is constant and equal to ν\nu. At u=ϑ=1u=\vartheta=1 the membrane potential is reset to u=ur=0u=u_{r}=0.

    (i) Formulate the continuity equation for the membrane potential density equation.

    (ii) Make the diffusion approximation

    (iii) Solve for the stationary membrane potential density distribution under the assumption that the flux through uru_{r} vanishes.

  4. 4.

    Linear response. A population is driven by a current I0+I1(t)I_{0}+I_{1}(t). The response of the population is described by

    A(t)=A0+A1(t)=A0+0G(s)I1(t-s)   𝑑   sA(t)=A_{0}+A_{1}(t)=A_{0}+\int_{0}^{\infty}G(s)\,I_{1}(t-s)\,{\text{d}}s\, (13.77)

    where GG is called the linear response filter.

    (i) Take the Fourier transformation and show that the convolution with the filter GG turns into a simple multiplication:

    A^1(ω)=G^(ω)I^1(ω)\hat{A}_{1}(\omega)=\hat{G}(\omega)\,\hat{I}_{1}(\omega) (13.78)

    where the hats denote the Fourier transformed variable.

    Hint: Replace G(s)G(s) in Eq. (13.77) by a causal filter Gc(s)=0G_{c}(s)=0 for s0s\leq 0 and Gc(s)=G(s)G_{c}(s)=G(s) for s>0s>0, extend the lower integral bound to --\infty, and apply standard rules for Fourier transforms.

    (ii) The squared quantity |I^1(ω)|2|\hat{I}_{1}(\omega)|^{2} is the power of the input at frequency ω\omega. What is the power |A^1(ω)|2|\hat{A}_{1}(\omega)|^{2} of the population activity at frequency ω\omega?

    (iii) Assume that the filter is given by G(s)=exp[-(s-Δ)/τg]G(s)=\exp[-(s-\Delta)/\tau_{g}] for s>Δs>\Delta and zero otherwise. Calculate G^(ω)\hat{G}(\omega).

    (iv) Assume that the input current has a power spectrum |I^1(ω)|2=(c/ω)|\hat{I}_{1}(\omega)|^{2}=(c/\omega) with c>0c>0 for ω>ω0\omega>\omega_{0} and zero otherwise.

    What is the power spectrum of the population activity with a linear filter as in (iii)?

  5. 5.

    Stability of stationary state .

    The response of the population is described by

    A(t)=A0+A1(t)=A0+0G(s)I1(t-s)   𝑑   sA(t)=A_{0}+A_{1}(t)=A_{0}+\int_{0}^{\infty}G(s)\,I_{1}(t-s)\,{\text{d}}s\, (13.79)

    where GG is called the linear response filter. Set G(s)=exp[-(s-Δg)/τg]G(s)=\exp[-(s-\Delta_{g})/\tau_{g}] for s>Δgs>\Delta_{g} and zero otherwise.

    Assume that the input arises due to self-coupling with the population:

    I1(t)=0α(s)A1(t-s)   𝑑   sI_{1}(t)=\int_{0}^{\infty}\alpha(s)\,A_{1}(t-s)\,{\text{d}}s\, (13.80)

    Set α(s)=(α0/τα)exp[-(s-Δα)/τα]\alpha(s)=(\alpha_{0}/\tau_{\alpha})\,\exp[-(s-\Delta_{\alpha})/\tau_{\alpha}] for s>Δαs>\Delta_{\alpha} and zero otherwise.

    (i) Search for solutions A1(t)exp[λ(ω)t]cos(ωt)A_{1}(t)\propto\exp[\lambda(\omega)t]\cos(\omega t). The stationary state A0A_{0} is stable if λ<0\lambda<0 for all frequencies ω\omega.

    (ii) Analyze the critical solutions λ=0\lambda=0 as a function of the delays ΔG\Delta_{G} and Δα\Delta_{\alpha} and the feedback strength α0\alpha_{0}.

  6. 6.

    Conductance input .

    Consider NEN_{E} excitatory and NIN_{I} inhibitory leaky integrate-and-fire neurons in the subthreshold regime

    Cdudt=-gL(u-EL)-gE(t)(u-EE)-gI(t)(u-EI)C{du\over dt}=-g_{L}\,(u-E_{L})-g_{E}(t)\,(u-E_{E})-g_{I}(t)\,(u-E_{I}) (13.81)

    where CC is the membrane capacity, gLg_{L} the leak conductance and EL,EE,EIE_{L},E_{E},E_{I} are the reversal potentials for leak, excitation, and inhibition, respectively. Assume that input spikes at excitatory arrive at a rate νE\nu_{E} and lead to a conductance change

    gE(t)=ΔgEjfexp[-(t-tj(f))/τE]fort>tj(f)g_{E}(t)=\Delta g_{E}\,\sum_{j}\sum_{f}\exp[-(t-t_{j}^{(f)})/\tau_{E}]\quad{% \rm for~{}}t>t_{j}^{(f)} (13.82)

    (and zero otherwise) with amplitude ΔgE\Delta g_{E} and decay time constant τE\tau_{E}. A similar expression holds for inhibition with ΔgI=2ΔgE\Delta g_{I}=2\Delta g_{E} and τI=τE/2\tau_{I}=\tau_{E}/2. Spike arrival rates are identical νI=νE\nu_{I}=\nu_{E}.

    (i) Determine the mean potential μ\mu.

    (ii) Introduce

    αE(t-tj(f))=aEexp[-(t-tj(f))/τE](t-tj(f))(μ-EE)\alpha_{E}(t-t_{j}^{(f)})=a_{E}\,\exp[-(t-t_{j}^{(f)})/\tau_{E}]{\mathcal{H}}(% t-t_{j}^{(f)})\,(\mu-E_{E}) (13.83)

    and an analogous expression for inhibitory input currents αI\alpha_{I}.

    Show that the membrane with conductance-based synapses Eq. (13.81) can be approximated by a model with current-based synapses

    τeffdudt=-(u-EL)+jNEfαE(t-tj(f))+jNIfαI(t-tj(f))\tau_{\rm eff}{du\over dt}=-(u-E_{L})+\sum_{j\in N_{E}}\sum_{f}\alpha_{E}(t-t_% {j}^{(f)})+\sum_{j\in N_{I}}\sum_{f}\alpha_{I}(t-t_{j}^{(f)}) (13.84)

    where ELE_{L} is the leak potential defined earlier in Eq. (13.81). Determine aE,aIa_{E},a_{I} and τeff\tau_{\rm eff}. What are the terms that are neglected in this approximation? Why are they small?

    (iii) Assume that the reversal potential for inhibition and leak are the same EI=ELE_{I}=E_{L}. What is the mean potential μ\mu in this case? How does inhibitory input manifest itself? What would change if you replaced inhibition by a constant current that sets the mean membrane potential (in the presence of the same amount of excitation as before) to μ\mu?

  7. 7.

    Firing rate of leaky integrate-and-fire neurons in the Brunel-Hakim formulation .

    Show that the Siegert formula of Eq. ( 13.30 ) can be also written in the form ( 78 )

    1A0τm=20   𝑑   ue-u2[e2y2u-e2y1uu]{1\over A_{0}\,\tau_{m}}=2\int_{0}^{\infty}{\text{d}}u\,e^{-u^{2}}\left[{e^{2y% _{2}u}-e^{2y_{1}u}\over u}\right] (13.85)

    with y2=(ϑ-h0)/σy_{2}=(\vartheta-h_{0})/\sigma and y1=(ur-h0)/σy_{1}=(u_{r}-h_{0})/\sigma .

    Hint: Use the definition of the error function erf{\rm erf} given above Eq. ( 13.30 )