A finite volume method for stochastic integrate-and-fire models

The stochastic integrate and fire neuron is one of the most commonly used stochastic models in neuroscience. Although some cases are analytically tractable, a full analysis typically calls for numerical simulations. We present a fast and accurate finite volume method to approximate the solution of the associated Fokker-Planck equation. The discretization of the boundary conditions offers a particular challenge, as standard operator splitting approaches cannot be applied without modification. We demonstrate the method using stationary and time dependent inputs, and compare them with Monte Carlo simulations. Such simulations are relatively easy to implement, but can suffer from convergence difficulties and long run times. In comparison, our method offers improved accuracy, and decreases computation times by several orders of magnitude. The method can easily be extended to two and three dimensional Fokker-Planck equations.

1. I modified the proof.

We thank Cheng Ly, Roberto Fernández Galán and Brent Doiron for helpful discussions. This research was supported by NSF grants DMS-0604429 and DMS-0817649, and a Texas ARP/ATP award.

Authors and Affiliations

1. Department of Mathematics, University of Houston, Houston, TX, 77204-3008, USA

   Fabien Marpeau, Aditya Barua & Krešimir Josić

Corresponding author

Correspondence to Fabien Marpeau.

Action Editor: G. Bard Ermentrout

Appendix A: Construction of numerical scheme (11–12) and Proof of Proposition 1

^{Footnote 1}First, consider the second order Taylor expansion

(26)

Then, integrating over one cell using \(\partial_{t}P=-\partial_{V}(Pf)\), \(\partial_{tt}P=\partial_{V }\big( f\partial_{V }(Pf)\big) =f\partial_{V V }(Pf)\!+\!\partial_{V }f\partial_{V }(Pf)\) and the integration by parts formula, one has

(27)

For the time being, let us neglect the terms containing Δt ², so that (26) reduces to the first order Taylor expansion. Then, one of the simplest stable schemes can be derived by using the first order upwind approximation \(P(t_{n},V_{i+{\frac{1}{2}} })f\big(V_{i+{\frac{1}{2}} }\big)\thickapprox P_{i}^{n}f_{i+{\frac{1}{2}} }^{+}+P_{i+1}^{n}f_{i+{\frac{1}{2}} }^{-}\) in the first term of the right-hand side of (27). We obtain the numerical scheme (11–12) with φ ≡ 0, which is first order accurate in space and time. Such accuracy does not properly capture some possible drift effects, such as sharp fronts (see for instance Bruneau et al. 2005; Godlewski and Raviart 1990). It has been observed that using a second order approximation of the space derivatives improves drastically the accuracy of numerical solutions to hyperbolic equations, even if the accuracy in time is still of first order. For this reason, we use instead the centered discretization \(P\big(t_{n},V_{i+{\frac{1}{2}} }\big)f\big(V_{i+{\frac{1}{2}} }\big)\thickapprox f_{i+{\frac{1}{2}} }\big( \gamma_{i+{\frac{1}{2}} }P_{i}^{n}+\big(1-\gamma_{i+{\frac{1}{2}} }\big)P_{i+1}^{n}\big) \) in the first term of the right-hand side of (27), where the coefficient \(\gamma _{i+{\frac{1}{2}} }=\frac{1}{2}\frac{\Delta V_{i+1}}{\Delta V_{i+{\frac{1}{2}} }}\) interpolates the numerical solution at the interface \(V_{i+{\frac{1}{2}} }\), thus allowing us to handle variable size mesh elements. Finally, we further improve the accuracy of our scheme by including the term \(\Big[ f^{2}\partial_{V } P(t_{n},.) \Big]_{V_{i-{\frac{1}{2}} }}^{V_{i+{\frac{1}{2}}}}\) and using a centered approximation, \(f^{2}\big(V_{i+{\frac{1}{2}} }\big)\partial_{V } P\big(t_{n},V_{i+{\frac{1}{2}} }\big) \thickapprox f_{i+{\frac{1}{2}} }^{2}\frac{P_{i+1}^{n}-P_{i}^{n}}{\Delta V }\) for all i. This results in the numerical scheme (11–12) with φ ≡ 1, which is more accurate in space and time than the upwind scheme, but unstable. Also notice that it becomes second order accurate in time if the drift term is constant, hence cancelling out the terms containing \(\partial_{V}f\) in Eq. (27).

Indeed, the flux (12) is a perturbation of the upwind scheme, so a compromise between accuracy and stability is offered by selecting φ satisfying (13).

In order to prove Proposition 1, the resulting scheme (11–12) is written as \(P_{i}^{n+{\frac{1}{2}} }=P_{i}^{n}\Big( 1- \Delta t\frac{f_{i+{\frac{1}{2}} }-f_{i-{\frac{1}{2}} }}{\Delta V_{i}}\Big) + A_{i+{\frac{1}{2}} }^{n}\Delta P_{i+{\frac{1}{2}} }^{n} - B_{i-{\frac{1}{2}} }^{n}\Delta P_{i-{\frac{1}{2}} }^{n}\), with

Using property (13) and re-phrasing

one has \(A_{i+{\frac{1}{2}} }^{n}\!\ge\! 0\), \(B_{i-{\frac{1}{2}} }^{n}\!\ge\! 0\) and \(1 \!-\! \frac{\Delta t}{\Delta V_{i}} ( f_{i+{\frac{1}{2}} }-f_{i-{\frac{1}{2}} })^{+} - A_{i+{\frac{1}{2}} }^{n} - B_{i-{\frac{1}{2}} }^{n} \ge 0\) under conditions (15) and (16), implying that the scheme (11–12) is nonnegativity preserving and l ^∞–stable.

Appendix B: Proof of Proposition 2

Recursively assuming that \(P_{i}^{n}\ge 0\) for a given index n and all i, which is true for n = 0, Proposition 1 and boundary condition imply that \(P_{i}^{n+{\frac{1}{2}} }\ge 0\) for all i under conditions (15–16). Then, multiplying (11) by ΔV _i and summing from i = 1 to i = N, one has \(\| P^{n+{\frac{1}{2}} }\|_{l^{1}}=\| P^{n}\|_{l^{1}}\).

Next, P ^n + 1 is obtained by inverting the linear system MP ^n + 1 = S whose matrix M is given in Fig. 3 and the right-hand side S is defined by \(S_{i}=\Delta V_{i}P_{i}^{n+{\frac{1}{2}} }\) if i ≠ i _R , and \(S_{i_{R}}=\Delta V_{i_{R}}P_{i_{R}}^{n+{\frac{1}{2}} }-\int_{t_{n}-\tau_{R}}^{ \min (t_{n+1}-\tau_{R},t_{n}) }F(t)\, dt\).

Indeed, the matrix M is strongly column diagonally dominant with positive diagonal coefficients and nonpositive off-diagonal coefficients. So it is invertible and M ^− 1 is a positive matrix. Since the right-hand side S has nonnegative entries, the solution of the system, P ^n + 1, is component-wise nonnegative.

Finally, sum (24) and (19) from i = 1 to i = N to obtain

thus recursively implying (25).

Reprints and permissions