|Home | About | Journals | Submit | Contact Us | Français|
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.
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 . 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 . 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 , satisfies these requirements. Surprisingly, this approach is rarely used in physics and only recently has been applied to solve the FPE in simple problems . 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 .<<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 .<<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 :
where i(0)= i(2π), i=1,2, is the instantaneous phase of the oscillating neuron i; ω is the average angular firing frequency; Z() is the phase-dependent sensitivity  (a.k.a. phase response or phase-resetting curve) of the neuronal oscillator, which can be estimated experimentally ; ηi(t) are zero-mean, white-noise stochastic inputs, , 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 : two repetitions of the same rapidly fluctuating stimulus, η1(t)=η2(t), yield highly reproducible responses 1(t), 2(t), in a single neuron. This limit has been recently studied in detail by J. Teramae and D. Tanaka .
An interesting question for neuroscience is how intrinsic cell properties embodied in Z() 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() 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, 1-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 PP(1,2,t) of finding the system at point (1,2) at time t:
with initial condition at point (10,20): P(1,2,0) = δ(1 - 10)δ(2-20); periodic boundary conditions: P(0,2,t) = P(2π,2,t), P(1,0,t) = P(1,2π,t) ; and the normalization condition over the square domain Ω [0,2π) × [0,2π), at any time t: . 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 P/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 . There are two solutions of this problem. One is the trivial solution, , 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, . 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, is determined only up to an additive, arbitrary constant vector, :
Since the projection of the solution to the FPE, P onto the FEM space, 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, , 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 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() and increasing input correlation, r.
Figure 1 shows P(1,2) for r=0.6, for a prototypical phase response of type I neurons, Z() = N(1 + cos( +π)) and for a prototypical phase response of type II neurons, Z () =-M sin, being N and M normalization constants, such that the integral of the phase response squared is one (equation parameters: ω=1, σ1=σ2=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(1, 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(1,1+π) for r=1 as the undetermined offset and assuming to be the same for all r, we subtracted it from P.
From P(1,2) we can easily obtain the “probability density” P(Δ) of the phase difference Δ=2-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 . 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.
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 , 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  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 .
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.