Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Rev E Stat Nonlin Soft Matter Phys. Author manuscript; available in PMC 2010 December 28.
Published in final edited form as:
Phys Rev E Stat Nonlin Soft Matter Phys. 2007 November; 76(5 Pt 2): 056110.
Published online 2007 November 15. doi:  10.1103/PhysRevE.76.056110
PMCID: PMC3010866

Solving the Fokker-Planck equation with the finite-element method

An example studying stochastic synchronization of neuronal oscillators


We apply an efficient approach from computational engineering, the finite-element method, to numerically solve the Fokker-Planck equation in two dimensions. This approach permits us to find the solution to stochastic problems that cannot be solved analytically. We illustrate our strategy with an example from neuroscience that recently has attracted considerable attention - synchronization of neural oscillators. In particular, we show that resonators (type II neural oscillators) respond and synchronize more reliably when provided correlated stochastic inputs than do integrators (type I neural oscillators). This result is consistent with recent experimental and computational work. We briefly discuss its relevance for neuroscience.

Keywords: stochastic synchronization, neuronal reliability, integrators, resonators, computational engineering


The interest in stochastic processes has increased remarkably in the last few years, in part motivated by the need to account for the effects of noise in biological systems [1]. A quantitative description of these phenomena often requires the solution of the Fokker-Planck equation (FPE), which is derived from the equations of motion modeling stochastic processes [2]. The FPE, however, rarely has an analytical solution in more than one dimensions and when it has an analytical solution, it usually involves integrals of rapidly increasing functions which make even numerical evaluation difficult. Furthermore, in many situations, especially those of current interest in biological sciences, boundary conditions must be added to the FPE, as well as spatial or temporal correlations of the stochastic driving forces. It is usually in these cases, however, when an unexpected constructive role for noise may emerge (e.g., from temporal correlations: [3-5]; from spatial correlations: [6-9]).

For these reasons, an efficient numerical approach to solving the FPE and related equations must be easy to generalize to more than one dimensions; it must be adaptive, so that integration can be accurate in those parts of the spatiotemporal domain where the actual solution changes quickly; and its computer implementation should be as efficient as possible. The finite-element method (FEM), an approach for solving partial differential equations (PDE) that is widely used in engineering [10], satisfies these requirements. Surprisingly, this approach is rarely used in physics and only recently has been applied to solve the FPE in simple problems [11]. Here we use the FEM to calculate the numerical solution of a two-dimensional stochastic problem with boundary conditions that cannot be solved analytically and where the application of perturbation theory in some limit cases is not sufficient to gain insight into the general behavior of the underlying physical phenomenon, namely, noise-induced synchronization of phase oscillators.

In the FEM, the spatial domain [should you say what you mean by spatial domain?] of interest is covered with a mesh of N knots, which need not be evenly spaced. Then, a set of N basis functions is defined over the mesh. The basis typically consists of “tent functions”: the n-th tent function, νn takes value one on the n-th knot and zero anywhere else on the mesh. Then, the PDE is projected onto the basis set by multiplying the equation by νn and integrating over the entire domain. Thus, one ends up with a matrix equation that in the case of the FPE has the general form:


where ui is the numerical solution to the PDE on the i-th knot, with i=1,..,N. The value of the solution on an arbitrary point of the spatial domain that is not a knot can be found by simple interpolation. The matrix Aij (which is very sparse) and the vector fi (which is derived from the inhomogeneous terms of the PDE and of the Neumann’s boundary conditions) are time independent if so are the coefficients of the PDE and the boundary conditions. If one is interested only in the stationary solution, the left hand side of equation (1) can be directly set to zero.

At present, there are several commercial and freely available software packages of the FEM in two and three dimension. Here we have used the package freeFEM++ developed by F. Hecht, O. Pironneau, A. Le Hyaric and K. Ohtsuka, which can be freely downloaded from the Internet [12].<<I would avoid this text and just list this as a reference when you say that we solved this system...


Populations of periodically firing neurons will fire synchronized action potentials when they receive correlated fluctuations as input [8].<<This first sentence is unclear. I am not sure that my rephrasing is what you mean. This phenomenon, called stochastic synchrony, can be modeled with two identical phase oscillators that are driven by correlated stochastic inputs starting with different initial conditions [13]:


where [var phi]i(0)= [var phi]i(2π), i=1,2, is the instantaneous phase of the oscillating neuron i; ω is the average angular firing frequency; Z([var phi]) is the phase-dependent sensitivity [14] (a.k.a. phase response or phase-resetting curve) of the neuronal oscillator, which can be estimated experimentally [15]; ηi(t) are zero-mean, white-noise stochastic inputs, <ηi(t)ηi(tτ)>=σi2δ(τ), that are spatially correlated <η1(t)η2(t)>=c, thus, the correlation coefficient of the inputs is r=c/(σ1σ2). As inputs become more correlated (i.e. in the limit of r→1) stochastic synchrony is equivalent to spike-time reliability [16]: two repetitions of the same rapidly fluctuating stimulus, η1(t)=η2(t), yield highly reproducible responses [var phi]1(t), [var phi]2(t), in a single neuron. This limit has been recently studied in detail by J. Teramae and D. Tanaka [7].

An interesting question for neuroscience is how intrinsic cell properties embodied in Z([var phi]) influence the amount of stochastic synchronization for each level of input correlation, r. In a first, but useful approximation, one can classify single neural dynamics in two major types: integrators, or type I, whose phase response Z([var phi]) is nonnegative, and resonators or type II, whose phase response is partially positive and partially negative. To study the effect of the phase response on stochastic synchronization, one can numerically integrate system (2) and, e.g., investigate the histograms of the phase difference, [var phi]1-[var phi]2 as a function of r, for different phase responses. Alternatively, for a more thorough study, one can use stochastic theory. Consider the FPE for stochastic system (2) that determines the probability P[equivalent]P([var phi]1,[var phi]2,t) of finding the system at point ([var phi]1,[var phi]2) at time t:


with initial condition at point ([var phi]10,[var phi]20): P([var phi]1,[var phi]2,0) = δ([var phi]1 - [var phi]10)δ([var phi]2-[var phi]20); periodic boundary conditions: P(0,[var phi]2,t) = P(2π,[var phi]2,t), P([var phi]1,0,t) = P([var phi]1,2π,t) ; and the normalization condition over the square domain Ω [subset or is implied by] [0,2π) × [0,2π), at any time t: ΩP(φ1,φ2,t)dφ1dφ2=1. This problem can readily be solved with the FEM at each point in time. For our purposes, however, we will focus on the stationary probability distribution, which satisfies eq. (3) setting [partial differential]P/[partial differential]t=0.

Note that eq. (3) involves mixed derivatives and non-constant coefficients. In principle, a change of variables can be applied to remove the cross derivatives. However, the boundary conditions would then mix the new spatial variables. Furthermore, the coefficients of the derivatives are complicated functions of the spatial variables. As a result, problem (3) does not have an obvious analytical solution. Nevertheless, it can be efficiently solved with the FEM. To do so, we first project eq. (3) onto the finite elements of a square mesh with N=50×50 regularly spaced knots covering the domain, Ω. We thus obtain a homogeneous algebraic equation satisfied by the finite-element representation of the stationary distribution. In matrix notation, we have Au=0. There are two solutions of this problem. One is the trivial solution, u=0, and the other is any vector belonging to the null space of A, i.e. an eigenvector with zero eigenvalue, which in this case is a constant vector, A1=0. This property is a consequence of lacking Dirichlet boundary conditions, i.e., it is due to the fact that the exact solution is not given along any part of the boundary a priori. This in turn implies that the solution to the linear problem with any right hand side, b is determined only up to an additive, arbitrary constant vector, a:


Since the projection of the solution to the FPE, P onto the FEM space, u is defined up to an additive constant, it seems reasonable to go back to the original problem and try the change of variables: Q=P+1. This leads to an inhomogeneous FPE in Q which in turn leads to an inhomogeneous algebraic problem, Au=b, whose solution is actually not determined because A is not invertible due to the zero eigenvalue. However, we can perturb the diagonal of A with a small constant, ε. This shifts the zero eigenvalue to ε, turning A into a full-rank matrix, which allows us to solve the linear problem. The solution obtained is neither zero nor constant and has a bias term proportional to 1/ε, but since u is defined only up to an arbitrary constant, we can ignore this term. In summary, with the FEM we obtain the stationary probability distribution up to an additive constant (which also precludes the application of the normalization condition). Nevertheless, this solution, which for simplicity we will still call P, is sufficient to study stochastic synchronization of system (2) for different phase-response curves, Z([var phi]) and increasing input correlation, r.

Figure 1 shows P([var phi]1,[var phi]2) for r=0.6, for a prototypical phase response of type I neurons, Z([var phi]) = N(1 + cos([var phi] +π)) and for a prototypical phase response of type II neurons, Z ([var phi]) =-M sin[var phi], being N and M normalization constants, such that the integral of the phase response squared is one (equation parameters: ω=1, σ12=1). At this level of intermediate input correlation, it is clear that both, type I and type II pairs tend to synchronize on average as indicated by the white band along the diagonal. However, stochastic synchronization is larger for type II. Although not necessary, we have made an educated guess on the undetermined additive constant of P([var phi]1, [var phi]2). We know that in the limit of r=1, synchrony must be maximal whereas anti-synchrony must be close to zero. Thus, we took the minimal value of P([var phi]1,[var phi]1) for r=1 as the undetermined offset and assuming to be the same for all r, we subtracted it from P.

Figure 1
Probability density (up to an additive, arbitrary constant) P([var phi]1,[var phi]2). The probability of the synchronous states, [var phi]1=[var phi]2 is larger for type II than for type I. The correlation coefficient of the stochastic inputs was r=0.6. ...

From P([var phi]1,[var phi]2) we can easily obtain the “probability density” P(Δ) of the phase difference Δ=[var phi]2-[var phi]1. As shown in Fig. 2, the probability of synchrony, P(Δ=0) increases with increasing input correlation for both neural types, but it is larger for type II. Figure 2 is in agreement with the phase-difference histograms obtained from numerical integration of system (2) presented in [13]. To quantify the increase of stochastic synchronization as a function of the input correlation, we calculate the area of P(Δ) between -π/4 and π/4 relative to the total area and plot this result in Fig. 3. Clearly, synchronization monotonically improves with increasing correlation of the stochastic inputs in both cases, but type II neurons are more efficient at synchronizing than type I are for each r, except for the limits r=0 and r=1 where both are equal.

Figure 2
Probability density of the phase difference, P(Δ). Stochastic synchronization increases with increasing input correlation, but it is larger for type II than for type I. This result is in agreement with computer simulations of system (2) presented ...
Figure 3
Stochastic synchronization improves with increasing input correlation, r but it is systematically larger for type II than for type I.


We have presented an application of the FEM to the solution of stochastic problems that frequently emerge in the study of physical and biological processes. In particular, the FEM has permitted us to investigate an important phenomenon of neuroscience, namely, stochastic synchronization of neurons that cannot be treated analytically in the general case. We have shown that neuronal resonators (type II excitable neurons) are more reliable and more susceptible to synchronize through stochastic correlated inputs than integrators (type I excitable neurons) are. These correlated inputs closely resemble the synaptic currents observed in real neurons [8], suggesting that neurons that have to reliably encode sensory information, as well as those that are likely to process or route information through synchrony and accurate spike timing, should have evolved to resonators. On the other hand, neurons that relay information through firing rates rather than spike timing might have evolved to integrators. Interestingly, we have previously shown that principal neurons in the olfactory system are indeed resonators [15,17]. Moreover, recent work by T. Tateno and H.P. Robinson [18] shows that fast-spiking cells in the neocortex, which are thought to synchronize generating fast network oscillations in cortical areas, are type II, whereas regularly-spiking cells that communicate across cortical columns are type I. Furthermore, in agreement with the predictions presented here, the cited authors report a remarkably higher reliability of the responses to fluctuating inputs in fast-spiking cells than in regularly-spiking cells [18].

In conclusion, the FEM has permitted us to investigate in detail a fundamental problem in current neuroscience. We believe that the physicists’ community will considerably profit in the future from the application of this approach of computational engineering to physical and biological problems.


[1] Springer M, Paulsson J. Nature. 2006;439:27. [PubMed]
[2] Gardiner CW. Handbook of stochastic methods for physics, chemistry, and the natural sciences. Springer; Berlin: 2004.
[3] Doering CR, Dontcheva LA, Klosek MM. Chaos. 1998;8:643. [PubMed]
[4] Middleton JW, Chacron MJ, Lindner B, Longtin A. Physical Review e. 2003;68 [PubMed]
[5] Salinas E, Sejnowski TJ. Neural Comput. 2002;14:2111. [PMC free article] [PubMed]
[6] Zhou T, Chen L, Aihara K. Phys. Rev. Lett. 2005;95:178103. [PubMed]
[7] Teramae J, Tanaka D. Phys. Rev. Lett. 2004;93:204103. [PubMed]
[8] Galán RF, Fourcaud-Trocmé N, Ermentrout GB, Urban NN. J. Neurosci. 2006;26:3646. [PubMed]
[9] Galán RF, Ermentrout GB, Urban NN. Sensors and Actuators B: Chemical. 2006;116:168.
[10] Zienkiewicz OC, Taylor RL, Zhu JZ. The Finite Element Method: Its Basis and Fundamentals. Elsevier Butterworth-Heinemann; 2005.
[11] Kumar P, Narayanan S. Sadhana. 2006;31:445.
[13] Galán RF, Ermentrout GB, Urban NN. Neurocomp. 2007 to be published.
[14] Kuramoto Y. Chemical Oscillations, Waves, and Turbulence. Dover Publications, Inc.; Mineola, New York: 2003.
[15] Galán RF, Ermentrout GB, Urban NN. Phys. Rev. Lett. 2005;94:158101. [PMC free article] [PubMed]
[16] Mainen ZF, Sejnowski TJ. Science. 1995;268:1503. [PubMed]
[17] Galán RF, Ermentrout GB, Urban NN. Neurocomp. 2006;69:1112.
[18] Tateno T, Robinson HP. J. Neurophysiol. 2006;95:2650. [PubMed]