Stochastic dynamics of uncoupled neural oscillators: Fokker-Planck studies with the ﬁnite element method

We have investigated the effect of the phase response curve on the dynamics of oscillators driven by noise in two limit cases that are especially relevant for neuroscience. Using the ﬁnite element method to solve the Fokker-Planck equation we have studied (cid:1) i (cid:2) the impact of noise on the regularity of the oscillations quantiﬁed as the coefﬁcient of variation, (cid:1) ii (cid:2) stochastic synchronization of two uncoupled phase oscillators driven by correlated noise, and (cid:1) iii (cid:2) their cross-correlation function. We show that, in general, the limit of type II oscillators is more robust to noise and more efﬁcient at synchronizing by correlated noise than type I.


INTRODUCTION
The interest in stochastic processes has greatly increased in the last few years, in part motivated by attempts to understand 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,3͔. In many situations, especially those of current interest in biological sciences, boundary conditions as well as spatial or temporal correlations of the stochastic driving forces must be taken into account. It is usually in these cases when an unexpected constructive role for noise may emerge ͑e.g., from temporal correlations ͓4-6͔ and from spatial correlations ͓7-11͔͒. Motivated by recent advances in experimental and computational neuroscience we have studied several properties of the stochastic dynamics of uncoupled phase oscillators. In particular, we have investigated the effect of two general cases of the phase response curve ͑PRC͒ on the dynamics of oscillatory neurons driven by noise: how it contributes to the robustness of the oscillations and how it helps synchronize uncoupled neurons. Whereas important properties of noise-driven phase oscillators can be studied analytically, like relative-phase densities ͓7͔ and Liapunov exponents ͓9͔, we show here that the calculation of quantities like the joint probability of the phases or their crosscorrelation function benefits from the application of an efficient numerical approach: the finite element method ͑FEM͒. Although the FEM has been applied before to solve the FPE in other contexts ͓12-15͔, it is still scarcely used in physics.
Here we show its potential by studying the effect of correlated noise in stochastic systems whose dynamics cannot be derived from a generalized energy function or where the application of perturbation theory in some limit cases is not sufficient to gain insight into the general case.
In the FEM, the spatial region 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 nth tent function v n takes a value of 1 on the nth knot and 0 anywhere else on the mesh. Then, the partial differential equation ͑PDE͒ is projected onto this non-orthogonal basis set by multiplying the equation by v 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 u i is the numerical solution to the PDE on the ith 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 interpolation. The matrix A ij and the vector f i ͑which is derived from the inhomogeneous terms of the PDE and Neumann's boundary conditions͒ are sparse. In general, for the FPE of stochastic systems with either additive or multiplicative noise, f i = 0, which leads to a homogeneous, algebraic equation for the steady state of Eq. ͑1͒. Several software packages are available that implement the FEM in two and three dimensions. Here we have used FREEFEM + + developed by Hecht et al. ͓16͔.

RESULTS
Ensembles of uncoupled oscillating neurons will fire synchronously when they receive correlated fast fluctuations as input ͓10͔. 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 ͓17͔: where i ͑0͒ = i ͑2͒, i =1,2, is the instantaneous phase of the oscillating neuron i; i is the average angular firing fre-*galan@cnbc.cmu.edu quency; and Z͑͒ is the phase-dependent sensitivity ͓18͔ ͑also known as phase response or phase-resetting curve͒ of the neuronal oscillator, which can be estimated experimentally ͓19͔. Z͑͒ describes the change in phase as a function of the phase at which an input arrives. The i ͑t͒ are zeromean, white-noise stochastic inputs, which 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 in the case of identical neural oscillators is equivalent to spike-time reliability ͓20͔: two repetitions of the same rapidly fluctuating stimulus, 1 ͑t͒ = 2 ͑t͒, yield highly reproducible responses 1 ͑t͒ and 2 ͑t͒ in a single neuron.
We have analyzed how the intrinsic membrane properties embodied in Z͑͒ influence the degree of stochastic synchronization generated by a given level of input correlation, r. To a first approximation, one can classify the dynamics of a single neuron into two categories based on their phase response: integrators or type-I neurons have non-negative phase response curves Z͑͒ ͓Fig. 1͑a͒, left͔ whereas resonators or type-II neurons have phase response curves that are partially positive and partially negative ͓Fig. 1͑a͒, right͔. To study the effect of the phase response on stochastic synchronization, one can numerically integrate the stochastic differential equations ͑SDEs͒ in Eqs. ͑2͒, as shown in Fig. 1͑b͒, and, e.g., investigate the histograms of the phase difference 1 − 2 as a function of r for different phase response curves ͓17͔. Alternatively, for a more efficient calculation of probability distributions and other statistical quantities, one can use stochastic theory and the FEM. In effect, the calculation of the probability distributions below required only a few seconds with the FEM whereas their estimation from the integration of the SDE required several minutes to obtain comparable accuracy on the same computer. Furthermore, the FEM does not require any a priori knowledge of the solution along the boundaries if an appropriate procedure is applied, as shown below. This is a major advantage with respect to other numerical schemes to solve partial differential equations, like the finite-difference method, where the solution over the inner domain is obtained by propagating the values from the boundaries.
The forward Fokker-Planck differential operator for the stochastic system ͑2͒ is given by and the backward ͑adjoint͒ Fokker-Planck operator is given by where the differential operators act to the right, as usual. Before studying the two-dimensional problem, we will focus on the effect of noise on a single phase oscillator ͓therefore, ignoring the terms in ͑3͒ and ͑4͒ with 2 , 2 , and c͔. In particular, we quantify the impact of noise on the regularity of the oscillations as the coefficient of variation ͑CV͒ of the oscillator dynamics. The first, T 1 , and second, T 2 , temporal moments of the backward FPE, L FP  1. ͑a͒ Phase response curves for two limit cases from neuroscience. Left, type-I excitability or neural integrator; right, type-II excitability or neural resonator. ͑b͒ Phase evolution of two oscillators ͑black and gray lines͒ and threshold crossings ͑i.e., spikes as black and gray dots͒ for different levels of input correlation r. The number of synchronous spikes clearly increases with r and appears to be larger for type II ͑right͒ than for type I ͑left͒, as confirmed in Fig. 3. Equation pa- Here, T 1 ͑͒ is the mean time required to go from to the end of the cycle ͑ =2͒ and T 2 ͑͒ − T 1 2 ͑͒ is the variance. Thus, the CV is calculated as the standard deviation of the cycle duration divided by the mean of the cycle duration: . Figure 2͑a͒ displays the results of this calculation as a function of the noise amplitude and compares it with the CV obtained from integration of the SDE ͑2͒, showing good agreement between both approaches. Clearly, for each the CV is lower for a type-II oscillator, indicating more robustness to noise than a type-I oscillator.
We now study the two-dimensional case of two phase oscillators driven by correlated noise. The probability P ϵ P͑ 1 , 2 ; t͒ of finding the system at point ͑ 1 , 2 ͒ at time t is determined by the forward FPE L FP P͑ 1 , 2 ;t͒ = ‫ץ‬ P͑ 1 , 2 ;t͒/‫ץ‬t, ͑5͒ with the 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. ͑5͒, setting ‫ץ‬P / ‫ץ‬t =0. Note that Eq. ͑5͒ involves mixed derivatives and nonconstant 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 second derivatives are complicated functions of the spatial variables. As a result, problem ͑5͒ does not have an obvious analytical solution. Nevertheless, it can be efficiently solved with the FEM. To do so, we first project Eq. ͑5͒ onto the finite elements of a square mesh with N =50ϫ 50 ͑also with N = 100ϫ 100, obtaining identical results͒ 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 ជ ͑note that the FPE leads to this homogeneous, algebraic equation for stochastic systems with either additive or multiplicative noise͒. 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. To avoid ending with an inhomogeneous algebraic system that only leads to these trivial solutions, we go back to the original stationary problem and perform the change of variables Q = P + k, where k is a constant that is determined below. This leads to an inhomogeneous FPE in Q, which in turn leads to an inhomogeneous algebraic problem Au ជ = b ជ . The solution of this system is not determined because A is not invertible due to the zero eigenvalue. However, if we remove an arbitrary equation of the system and set the corresponding component of u ជ to be zero, the solution will be determined. In addition,  2. ͑a͒ Coefficient of variation of the cycle duration as a function of the noise amplitude ͑lines, FEM; circles, SDE͒. ͑b͒ Stochastic synchronization quantified as the correlation coefficient of the phases, R, increases with increasing input correlation r, but it is consistently larger for type II than for type I ͑dashed lines, R = r for reference͒. ͑c͒ Effect of the difference of the intrinsic frequencies on stochastic synchronization. As the difference increases, synchrony deteriorates quickly, but more rapidly for type I. ͑d͒ The noise amplitude does not affect synchrony as long as the phase oscillator approximation is valid ͑ = 1, larger than the average of the term with ͒.
we require the integral of Q over ⍀ to vanish. Then, applying the normalization condition ͑6͒ for P, we observe that k must be chosen as k =−1/͑2͒ 2 . Figure 3͑a͒ 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 , N and M being normalization constants, such that the integral of the absolute value of the phase response curves is 1 ͑equation parameters 1 = 2 =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 light gray band along the diagonal.
However, stochastic synchronization is more pronounced for type-II oscillators.
From P͑ 1 , 2 ͒ we can easily obtain the probability density P͑⌬͒ of the phase difference ⌬ = 2 − 1 . As shown in Fig. 3͑b͒ ͑thicker lines͒, the probability of synchrony, P͑⌬ =0͒, increases with increasing input correlation for both neural types, but it is larger for type II. Note the good agreement with the results of numerical integration of system ͑2͒, shown as thin lines with markers. To quantify the increase of stochastic synchronization as a function of the input correlation, we compute the cross-correlation coefficient R of the phases from P͑ 1 , 2 ͒ as Clearly, the synchronization R monotonically improves with increasing correlation of the stochastic inputs in both cases ͓Fig. 2͑b͔͒, 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. When r = 1 the diffusion matrix of the FPE is degenerate and the applicability of the equation is no longer valid. However, this limit has been recently studied with another approach by Teramae and Tanaka ͓9͔, showing that identical oscillators driven by identical random forces will perfectly synchronize ͑R =1͒, provided that Z͑͒ is continuous.
The former analysis can be easily generalized to other interesting conditions. In particular, we have studied the effect of the difference of the intrinsic frequencies on stochastic synchronization. Figure 2͑c͒ shows that as this difference increases, synchrony deteriorates more quickly in type-I than type-II pairs. In addition, we have investigated the effect of the amplitude of the correlated noise on stochastic synchronization ͓Fig. 2͑d͔͒. Interestingly, as long as the phase oscillator approximation is valid-i.e., if ͗Z͑͒͑t͒͘-the amplitude of the input noise, , plays no role on stochastic synchronization, but only its mutual correlation r.
The efficacy of the FEM also allows one to calculate the cross-correlation function C͑͒ easily. The cross-correlation function quantifies the ability to forecast the phase of one oscillator by knowing the current phase of the other. Stochastic theory demonstrates that C͑͒ can be computed from the eigenfunctions of the Fokker-Planck operator P ͑ 1 , 2 ͒ and its adjoint P + ͑ 1 , 2 ͒ ͓3͔: L FP P ͑ 1 , 2 ͒ = P ͑ 1 , 2 ͒, Using the normalization condition we define the autocorrelation function as where P 0 ͑ 1 , 2 ͒ is the stationary solution of the forward FPE calculated above. Note that since the Fokker-Plack operator ͑3͒ is not Hermitian, the eigenvalues can be complex. The complex eigenvalues are responsible for the oscillatory shape of C͑͒. In addition, the real part of the eigenvalues cannot be positive ͓3͔ so that the envelope of C͑͒ will decay with . Figure 3͑c͒ displays C͑͒ for type-I and type-II pairs of uncoupled oscillators ͑2͒ for several values of the input correlation r. In both cases, the peak increases with increasing r, whereas the decay rate decreases. However, for the same r, type II has a higher peak than type I. This result shows again that type-II oscillators are not only more efficient at synchronizing by correlated inputs than type I, but also tend to keep synchronized for longer times.

DISCUSSION
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. To illustrate this approach we have used the FEM to investigate an important phenomenon of neuroscience: namely, stochastic synchronization of neurons. In particular, we have shown that neural resonators ͑type-II excitable neurons͒ are more reliable and more readily to synchronize through correlated stochastic inputs than integrators ͑type-I excitable neurons͒. Our results can be interpreted intuitively as follows: consider two uncoupled neural oscillators at opposite extremes of the intrinsic period. According to the phase response ͓Fig. 1͑a͔͒, if they receive a correlated fluctuation, the phase difference of the integrators and of the resonators will evolve differently: whereas both integrators will move in the same direction ͓same sign of Z͔͑͒, and therefore without remarkably changing their phase difference, both resonators will move in opposite directions. However, because the phase is periodic ͑with period 2͒, moving in opposite directions actually means coming closer to each other. Thus, correlated fluctuations will tend to diminish the phase difference between resonators. This suggests that resonating neurons will more reliably encode the timing of sensory information and may process or route information through synchrony and accurate spike timing. On the other hand, neurons that relay information through firing rates rather than spike timing might have adapted to integrators. In fact, we have previously shown that principal neurons ͑mitral cells͒ in the olfactory system are indeed resonators ͓19͔. Although realistic synaptic inputs do not always resemble the noisy inputs considered here, we have previously demonstrated in a detailed experimental and computational study that realistic synaptic inputs and continuous noisy inputs lead to similar statistics of stochastic synchronization ͓10͔. Interestingly, synchronization of neural oscillations have been observed in vivo in the olfactory system of rodents during odor processing ͓21͔. The robustness of the oscillations around 40 Hz in mitral cells suggests that the stochastic synaptic inputs they receive are weak. Moreover, because mitral cells are only sparsely connected, the model presented here is a simple but realistic description of that system. Recent work by Tateno and Robinson ͓22,23͔ shows that fast-spiking cells in the neocortex, which are thought to synchronize generating fast network oscillations in cortical areas, are type II, whereas regular-spiking cells that communicate across cortical columns are type I. Furthermore, in agreement with the predictions presented here, the cited authors report a higher reliability of the responses to fluctuating inputs in fast-spiking cells than in regular-spiking cells ͓22͔. . ͑a͒ Probability density P͑ 1 , 2 ͒ on a gray color scale. The probability of the synchronous states, 1 = 2 , is larger for type II ͑right panels͒ than for type I ͑left panels͒. The correlation coefficient of the stochastic inputs here is r = 0.6 ͑parameters 1 = 2 =1, 1 = 2 =1͒, where 1 means small omega with subindex 1. ͑b͒ Probability density of the phase difference, P͑⌬͒ ͑thick lines͒.The probability of being around ⌬ = 0 is larger for type II than for type I for any intermediate value of r ͑only three representative values shown͒. This result is in agreement with the results obtained from numerical integration of system ͑2͒, shown as thin lines with markers. ͑c͒ Cross-correlation function C͑͒ normalized in such a way that C͑0͒ gives the correlation coefficient R ͑see text͒.