Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3010866

Formats

Article sections

Authors

Related links

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.056110PMCID: PMC3010866

NIHMSID: NIHMS58906

An example studying stochastic synchronization of neuronal oscillators

See other articles in PMC that cite the published article.

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 [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

$$\frac{d{u}_{i}}{dt}=\sum _{j=1}^{N}{A}_{ij}{u}_{j}+{f}_{i},$$

(1)

where *u _{i}* is the numerical solution to the PDE on the

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]:

$$\{\begin{array}{c}\frac{d{\phi}_{1}}{dt}=\omega +Z\left({\phi}_{1}\right){\eta}_{1}\left(t\right)\hfill \\ \frac{d{\phi}_{2}}{dt}=\omega +Z\left({\phi}_{2}\right){\eta}_{2}\left(t\right)\hfill \end{array}\phantom{\}}$$

(2)

where * _{i}(0)= _{i}(2π)*,

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

$$\begin{array}{cc}\hfill \frac{\partial P}{\partial t}=& -\frac{\partial}{\partial {\phi}_{1}}\left(\left(\omega +{\sigma}_{1}^{2}\frac{\partial Z\left({\phi}_{1}\right)}{\partial {\phi}_{1}}Z\left({\phi}_{1}\right)\right)P\right)-\frac{\partial}{\partial {\phi}_{2}}\left(\left(\omega +{\sigma}_{2}^{2}\frac{\partial Z\left({\phi}_{2}\right)}{\partial {\phi}_{2}}Z\left({\phi}_{2}\right)\right)P\right)+\hfill \\ \hfill & +\frac{{\sigma}_{1}^{2}}{2}\frac{{\partial}^{2}}{\partial {\phi}_{1}^{2}}\left({Z}^{2}\left({\phi}_{1}\right)P\right)+\frac{{\sigma}_{2}^{2}}{2}\frac{{\partial}^{2}}{\partial {\phi}_{2}^{2}}\left({Z}^{2}\left({\phi}_{2}\right)P\right)+c\frac{{\partial}^{2}}{\partial {\phi}_{1}\partial {\phi}_{2}}\left(Z\left({\phi}_{1}\right)Z\left({\phi}_{2}\right)P\right)\hfill \end{array}$$

(3)

with initial condition at point (* _{10},_{20}*):

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 $A\overrightarrow{u}=\overrightarrow{0}$. There are two solutions of this problem. One is the trivial solution, $\overrightarrow{u}=\overrightarrow{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, $A\overrightarrow{1}=\overrightarrow{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, $\overrightarrow{b}$ is determined only up to an additive, arbitrary constant vector, $\overrightarrow{a}$:

$$A(\overrightarrow{u}+\overrightarrow{a})=\overrightarrow{b}\to A\overrightarrow{u}+A\overrightarrow{a}=A\overrightarrow{u}+\overrightarrow{0}=A\overrightarrow{u}=\overrightarrow{b}.$$

Since the projection of the solution to the FPE, *P* onto the FEM space, $\overrightarrow{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, $A\overrightarrow{u}=\overrightarrow{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 $\overrightarrow{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()* 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

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

From *P*(* _{1},_{2}*) we can easily obtain the “probability density”

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.

[12] www.freefem.org.

[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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |