# Stochastic synchronization of neural activity waves

###### Abstract

We demonstrate that waves in distinct layers of a neuronal network can become phase-locked by common spatiotemporal noise. This phenomenon is studied for stationary bumps, traveling waves, and breathers. A weak noise expansion is used to derive an effective equation for the position of the wave in each layer, yielding a stochastic differential equation with multiplicative noise. Stability of the synchronous state is characterized by a Lyapunov exponent, which we can compute analytically from the reduced system. Our results extend previous work on limit-cycle oscillators, showing common noise can synchronize waves in a broad class of models.

Introduction. Nonlinear waves arise in many physical, biological, and chemical systems including non-equilibrium reactions Jakubith et al. (1990), shallow water Huang et al. (1999), bacterial populations Danino et al. (2010), epidemics Cummings et al. (2004), and cortical tissue Wang (2010). Phase synchronization of multiple waves can occur when their dynamics are coupled. For example, spiral waves in chemical systems become entrained when coupled diffusively through a membrane Winston et al. (1991). Experiments on amoeba populations have also demonstrated entrainment via interactions between wave-emitting centers and spiral waves of cell density Lee et al. (1996). Common forcing can also synchronize waves at distinct spatial locations. Spatiotemporal analysis of epidemics reveals that both seasonality and vaccination schedule can entrain the nucleation of outbreak waves across geographical space Grenfell et al. (2001). Furthermore, activity recordings from primary visual cortex show that triggering switches during binocular rivalry leads to synchronized wave initiation Lee et al. (2005). In total, experimental studies demonstrate a wide array of mechanisms for synchronizing the onset and propagation of waves.

Our goal in this Rapid Communication is to show that stochastic forcing can also entrain the phases of distinct waves. We focus on neural activity waves that arise due to distance-dependent synaptic interactions Ermentrout and Kleinfeld (2001); Wang (2010). Waves of neural activity underlie sensory processing Xu et al. (2007), motor action Rubino et al. (2006), and sleep states Massimini et al. (2004). Proposed computational roles of neural activity waves include heightening the responsiveness of specific portions of a network and labeling incoming signals with a distinct phase Ermentrout and Kleinfeld (2001). Thus, it may be advantageous for waves to be coordinated across multiple brain areas, and we propose that correlated fluctuations may underlie such coordination. Many theoretical and experimental studies have identified ways noise correlations can degrade neural information encoding Cohen and Kohn (2011), but recent work has shown correlated noise can reliably synchronize activity across populations of neurons Ermentrout et al. (2008).

Our analysis extends previous work which showed common noise can synchronize the phases of limit cycle oscillators Pikovsky et al. (2001); Teramae and Tanaka (2004). A key observation of these studies is that the synchronous state, where all oscillators have the same phase, is absorbing when each oscillator receives identical noise. Stability of the phase-locked state can then be determined by computing an associated Lyapunov exponent, which is negative for nontrivial phase response curves Teramae and Tanaka (2004). As we show in this Rapid Communication, these principles can be applied to waves forced by common spatiotemporal noise. Waves are driven to an attracting synchronous state, where the waves’ phases are identical.

Synchronization of neural activity waves across distinct network locations can be important in several behavioral and sensory contexts. For instance, waves of activity in the visual cortex may function like a “bar code scanner,” ensuring a portion of the network is always maximally sensitive to external inputs Ermentrout and Kleinfeld (2001). Since locations in visual space are represented by multiple layers of a network, coordinating background waves across layers could ensure the network is always sensitive at the same visual location in each layer. Thus, our findings implicate an important potential role for large-scale correlations in nervous system fluctuations, which are often deemed a nuisance to cognitive performance Cohen and Kohn (2011).

In this Rapid Communication, we analyze the stochastic dynamics of waves in a pair of uncoupled neural field models driven by common noise. Neural fields are nonlinear integrodifferential equations whose integral term describes the connectivity of a neuronal network Amari (1977); Bressloff (2012). Recent studies have considered stochastic versions of neural field equations, formulating the dynamics as Langevin equations with spatiotemporal noise Bressloff (2012); Hutt et al. (2007):

(1) |

where is neural activity of population at at time , synaptic connectivity is described by the convolution , and is a nonlinearity describing the fraction of active neurons.
Small amplitude () spatiotemporal noise is white in time and correlated in space, so and . As noise correlations must thus be even and -periodic, we write .

Stationary bumps. We begin by demonstrating noise-induced wave synchronization for a single realization of the system Eq. (1) in Fig. 1A num (). Here is an even symmetric lateral inhibitory weight function, known to lead to stable stationary bump solutions in the unperturbed system () Amari (1977). Additive spatiotemporal noise causes each bump to wander diffusively about the spatial domain. Once both waves’ positions and meet, they are phase-locked for the remainder of the simulation, since both layers receive identical noise and obey the identical governing Eq. (1).

To analyze this behavior, we first note that the unperturbed system, in Eq. (1), has stationary bump solutions () satisfying Ermentrout (1998); Kilpatrick and Ermentrout (2013). There is a degeneracy in the position of each bump’s center of mass , since the waves are neutrally stable to translating perturbations Kilpatrick and Ermentrout (2013), and noise will cause waves to stochastically wander about their mean position Panja (2004); Bressloff and Webber (2012).

We now derive effective equations for the bumps’ response to noise () by utilizing a perturbation expansion that tracks the stochastically varying position of each bump and fluctuations in the wave profiles , so . Bumps in each layer may begin at distinct locations , but we will show that when receiving common noise, their locations become synchronized in the long time limit .

Plugging this ansatz into Eq. (1) and expanding to , we have the linear stochastic equations

(2) |

for , where is a linear operator defined for integrable functions . Both layers receiving identical noise, so we must keep track of each wave’s relative phase . A solution to Eq. (2) exists if we require the right hand side be orthogonal to the nullspace of the adjoint operator . Define to be a one-dimensional basis of . Taking the inner product of Eq. (2) with , we find

(3) |

so obeys a Langevin equation with multiplicative noise

(4) |

We can represent using the Fourier expansion

(5) |

where and are white noise processes and ; is the Kronecker delta function. Using trigonometric identities, we can express

(6) |

where , and the vanishes since and

Thus, we have reduced Eq. (1) to a system describing two phase oscillators perturbed by weak noise, Eq. (6). The coefficients describe the relative contributions of each term of the Fourier series, Eq. (5), to the phase-dependent sensitivity of the oscillators to noise.

The synchronized solution to Eq. (6) is absorbing since the the right-hand sides of both equations () will subsequently be identical. To assess stability of the absorbing state, we compute the associated Lyapunov exponent . Proper calculation requires translating Eq. (6) into its equivalent Ito formulation Gardiner (2004)

(7) |

where the Ito Eq. (7) introduces the drift term

(8) |

accounting for the fact that the correlation between state variables and noise terms, present in the Stratonovich Eq. (6), subsequently vanishes. We proceed by formulating the variational equation for the perturbative phase difference (), which can be derived from Eq. (7), so

(9) |

where and obeys Eq. (7). Defining and appealing to Ito’s formula, we can rewrite Eq. (9) as

(10) |

where

(11) |

Subsequently, we can integrate Eq. (10) to determine the mean drift of

(12) |

which is also the mean rate of growth of . The phase difference will tend to decay (grow) if () and synchrony will be stable (unstable). Utilizing ergodicity of Eq. (10), we can equivalently compute with the ensemble average across realizations of , so Teramae and Tanaka (2004)

(13) |

where is the steady state distribution of . Since noise is weak (), we can approximate the distribution as constant, . Applying this to Eq. (13), we find the first term of the integrand vanishes since according to Eq. (8), so the Lyapunov exponent is approximated by the formula

(14) |

Note that as long as or for some , we expect , so the phase-locked state will be linearly stable.

We now compare the analytical result Eq. (14) to results from numerical simulations of Eq. (1). Explicit calculations are straightforward in the case of a Heaviside nonlinearity ; cosine weight function ; and cosine spatial noise correlations . Stable stationary bump solutions are given by the formulas and , and the null vector Kilpatrick and Ermentrout (2013). Coefficients of the noise Fourier components in Eq. (5) are , so and , . Finally, utilizing Eq. (11) along with Eq. (14), we have

(15) |

Results from numerical simulations in Fig. 1B,C corroborate with Eq. (15), showing the Lyapunov exponent’s magnitude increases with noise intensity and threshold . We note that our theoretical approximation breaks down as increases and the system nears a saddle-node bifurcation at which the stable/unstable branch of bump solutions annihilate in the noise-free system Amari (1977); Ermentrout (1998); Kilpatrick and Ermentrout (2013). Furthermore, the amplitude increases as the parameter is increased towards this bifurcation.

We show the robustness of these results by studying the impact of independent noise (Fig. 2A). To do so, we consider a modified version of Eq. (1), , where has an independent component in each layer (see Appendix). Extending analysis of limit cycle oscillators Nakao et al. (2007), we derive an expression for the stationary density

(16) |

of the phase difference . Here, is the degree of noise correlation, (), and is a normalization constant. As is decreased from unity, the stationary density widens, representing the effects of independent noise in each layer. However, the density will still tend to be peaked at (see Appendix). Indeed, our theory, Eq. (16), is corroborated by numerical simulations (Fig. 2B).

Traveling waves. Our results for stationary bumps can be extended to address stochastic synchronization of traveling waves in networks with asymmetric weights (), as in Fig. 3A. Thus, the unperturbed system, in Eq. (1), will have traveling wave solutions , , with wave speed (), so Ermentrout (1998); Kilpatrick and Ermentrout (2012). Furthermore, waves will be neutrally stable to translating perturbations, so spatiotemporal noise will cause an effective diffusion of their phases. For , we apply the ansatz with and will show . At , we find Eq. (2) with corresponding linear operator . Solvability is enforced by ensuring the right hand side of Eq. (2) is orthogonal to the nullspace of the adjoint , yielding the Langevin equation, Eq. (4). The Lyapunov exponent associated with the stability of the absorbing state is then approximated by Eq. (14).

To compare our analytical results for traveling waves with numerical simulations, we compute from Eq. (14) when , , and . Stable traveling waves have profile , width defined by thresholds where and , and speed Kilpatrick and Ermentrout (2012). The null vector can also be computed explicitly:

Fourier coefficients of in Eq. (6) are thus given and , , so we compute Eq. (14), finding

(17) |

comparing with numerical results in Fig. 3B.

Breathers. Lastly, we show that uncoupled oscillatory waves can also become phase-locked due to a common noise source. We extend Eq. (1) by incorporating linear adaptation as an auxiliary variable () and a spatially varying external input Pinto and Ermentrout (2001); Folias and Bressloff (2005); Ermentrout et al. (2014)

(18a) | ||||

(18b) |

where and are the rate and strength of adaptation. A detailed analysis of the onset of breathers in Eq. (18), via a Hopf bifurcation, can be found in Folias and Bressloff (2005); Ermentrout et al. (2014).

Our main interest is the rate at which breathers in a pair of uncoupled neural fields synchronize their phases when subject to common noise as in Eq. (18). An example of this phenomenon is shown in Fig. 4A. Note, we must take care in interpreting how the centers of mass of each oscillating bump relate to the phase of the underlying oscillation. In the case of bumps and waves, there was a one-to-one mapping between the wave positions and the phase of the stochastically-driven oscillator. Here, we must track the activity and adaptation variables to resolve the phases of the underlying oscillations. Assuming breathers have period , and so , and there is a mapping for all values along the trajectory of a breather. We save a more detailed analytical determination of this mapping for future work. Here, we numerically determine , average, and compute the rate of decay using the approximation using least squares (Fig. 4B). Note, as in the case of bumps, the Lyapunov exponent increases in amplitude as it nears the pattern-generating (Hopf) bifurcation.

Discussion. Our results demonstrate common spatiotemporal fluctuations in neuronal networks can synchronize the phases of waves. We have shown this for stationary bumps, traveling waves, and breathers. Since our derivations mainly rely on our ability to derive an effective equation for the relative position of a noise-driven wave, we suspect we could extend them to the case of traveling fronts Bressloff and Webber (2012) or Turing patterns Hutt et al. (2007) in stochastic neural fields. Patterns on two-dimensional (2D) domains could also be addressed by deriving multidimensional effective equations for each pattern’s position. For instance, a bump’s position would be represented with a 2D vector Poll and Kilpatrick (2014), so there would be two Lyapunov exponents associated with bumps’ phase-locked state. On the other hand, spiral waves would be characterized by a 2D position vector and a scalar phase Laing (2005), and phase/position-locked states would then have three Lyapunov exponents. We could also consider interlaminar coupling Kilpatrick (2014), exploring competing impacts of common noise and coupling on phase synchronization as in García-Álvarez
et al. (2009). Furthermore, these results should be applicable to nonlinear PDE models of reaction-diffusion systems Panja (2004); Sagués et al. (2007). Overall, our results suggest a novel mechanism for generating coherent waves in laminar media, presenting a testable hypothesis that could be probed experimentally.

Appendix. Here, we analyze the stochastic dynamics of bumps in a pair of uncoupled neural field models driven by both common and independent noise sources, extending previous work on limit cycle oscillators Nakao et al. (2007). We incorporate an independent noise term into each layer of the stochastic neural field model Kilpatrick (2014):

(19) |

Small amplitude () spatiotemporal noise terms () are white in time and correlated in space, so ; () with . The degree of correlation between layers is controlled by the parameter .

Our analysis proceeds by considering stationary bumps in a network with even symmetric connectivity (). As in the main text, we characterize stochastic bump motion by applying the ansatz , and . Plugging this ansatz into Eq. (19), expanding to , and applying a solvability condition, we find that each () obeys the Langevin equation

, where the first and second term correspond to correlated and independent noise. Here, is a one-dimensional basis of , where . Since is even symmetric, all components of the nullspace of are necessarily odd symmetric Kilpatrick and Ermentrout (2013). Note, we can represent (), where and are normalized white noise processes. We can thus use trigonometric expansions to express

(20) |

, where are multiplicative noise terms defined

and since is odd symmetric, vanishes, and

Eq. (20) can be reformulated as an Ito equation , where () has correlations defined (), and the drift (). Components of the correlation matrix are given

where and , so it is straightforward to compute ().

The corresponding Fokker-Planck equation, describing the coevolution of the position variables is thus

(21) | ||||

where . Note, since , then for . We can write Eq. (21) as a separable equation by employing a change of variables that tracks the average and phase difference of the phase variables and

(22) | ||||

where . Eq. (22) can be decoupled by plugging in the ansatz and noting the equation will be satisfied by the system

(23) | ||||

Thus, we can solve for the stationary solution of the system, Eq. (23), by setting and requiring periodic boundary conditions. The stationary distribution for the position average is . Furthermore, we can integrate the stationary equation for to find that the stationary density of the phase difference is

(24) |

where is a normalization factor. When noise to each layer is independent (, uncorrelated), then is constant in space. Since no common noise source entrains the phase of each bump, the bumps diffuse independently of one another. However, when noise is totally correlated between layers (), then . Thus, all initial conditions eventually result in the phase-locked state . The stationary distribution broadens as is decreased, with a peak still remaining at .

To compare our results to numerical simulations, we compute the stationary density explicitly by using ; ; and (). Stable stationary bumps satisfy the threshold condition , with half-width . We can thus compute the null vector and find and , . Therefore

(25) |

ZPK was funded by NSF-DMS-1311755. We thank Oliver Langhorne for helpful conversations.

## References

- Jakubith et al. (1990) S. Jakubith, H. Rotermund, W. Engel, A. Von Oertzen, and G. Ertl, Phys Rev Lett 65, 3013 (1990).
- Huang et al. (1999) N. E. Huang, Z. Shen, and S. R. Long, Annu Rev Fluid Mech 31, 417 (1999).
- Danino et al. (2010) T. Danino, O. Mondragón-Palomino, L. Tsimring, and J. Hasty, Nature 463, 326 (2010).
- Cummings et al. (2004) D. A. Cummings, R. A. Irizarry, N. E. Huang, T. P. Endy, A. Nisalak, K. Ungchusak, and D. S. Burke, Nature 427, 344 (2004).
- Wang (2010) X.-J. Wang, Physiol Rev 90, 1195 (2010).
- Winston et al. (1991) D. Winston, M. Arora, J. Maselko, V. Gáspár, and K. Showalter, Nature 351, 132 (1991).
- Lee et al. (1996) K. J. Lee, E. C. Cox, and R. E. Goldstein, Phys Rev Lett 76, 1174 (1996).
- Grenfell et al. (2001) B. Grenfell, O. Bjørnstad, and J. Kappey, Nature 414, 716 (2001).
- Lee et al. (2005) S.-H. Lee, R. Blake, and D. J. Heeger, Nat Neurosci 8, 22 (2005).
- Ermentrout and Kleinfeld (2001) G. B. Ermentrout and D. Kleinfeld, Neuron 29, 33 (2001).
- Xu et al. (2007) W. Xu, X. Huang, K. Takagaki, and J.-y. Wu, Neuron 55, 119 (2007).
- Rubino et al. (2006) D. Rubino, K. A. Robbins, and N. G. Hatsopoulos, Nat Neurosci 9, 1549 (2006).
- Massimini et al. (2004) M. Massimini, R. Huber, F. Ferrarelli, S. Hill, and G. Tononi, J Neurosci 24, 6862 (2004).
- Cohen and Kohn (2011) M. R. Cohen and A. Kohn, Nat Neurosci 14, 811 (2011).
- Ermentrout et al. (2008) G. B. Ermentrout, R. F. Galán, and N. N. Urban, Trends Neurosci 31, 428 (2008).
- Pikovsky et al. (2001) A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge Nonlinear Series, 2001).
- Teramae and Tanaka (2004) J.-n. Teramae and D. Tanaka, Phys Rev Lett 93, 204103 (2004).
- Amari (1977) S. Amari, Biol Cybern 27, 77 (1977).
- Bressloff (2012) P. C. Bressloff, J Phys. A: Math. Theor. 45, 033001 (2012).
- Hutt et al. (2007) A. Hutt, A. Longtin, and L. Schimansky-Geier, Phys Rev Lett 98, 230601 (2007).
- (21) Numerical simulations used Euler-Maruyama method along with a direct spectral method for evaluating integrals. Timesteps and spatial steps . Numerically evaluated Lyapunov exponents were determined by averaging across realizations and fitting to using MATLAB’s polyfit function.
- Ermentrout (1998) B. Ermentrout, Rep Prog Phys 61, 353 (1998).
- Kilpatrick and Ermentrout (2013) Z. P. Kilpatrick and B. Ermentrout, SIAM J Appl Dyn Syst 12, 61 (2013).
- Panja (2004) D. Panja, Phys Rep 393, 87 (2004).
- Bressloff and Webber (2012) P. C. Bressloff and M. A. Webber, SIAM J Appl Dyn Syst 11, 708 (2012).
- Gardiner (2004) C. W. Gardiner, Handbook of stochastic methods for physics, chemistry, and the natural sciences (Springer-Verlag, Berlin, 2004).
- Nakao et al. (2007) H. Nakao, K. Arai, and Y. Kawamura, Phys Rev Lett 98, 184101 (2007).
- Kilpatrick and Ermentrout (2012) Z. P. Kilpatrick and B. Ermentrout, Phys Rev E 85, 021910 (2012).
- Pinto and Ermentrout (2001) D. J. Pinto and G. B. Ermentrout, SIAM J Appl Math 62, 226 (2001).
- Folias and Bressloff (2005) S. E. Folias and P. C. Bressloff, Phys Rev Lett 95, 208107 (2005).
- Ermentrout et al. (2014) G. B. Ermentrout, S. E. Folias, and Z. P. Kilpatrick, in Neural Fields: Theory and Applications (Springer, 2014).
- Poll and Kilpatrick (2014) D. B. Poll and Z. P. Kilpatrick, arXiv:1412.3410 (2014).
- Laing (2005) C. R. Laing, SIAM Journal on Applied Dynamical Systems 4, 588 (2005).
- Kilpatrick (2014) Z. P. Kilpatrick, Physical Review E 89, 022706 (2014).
- García-Álvarez et al. (2009) D. García-Álvarez, A. Bahraminasab, A. Stefanovska, and P. McClintock, EPL (Europhysics Letters) 88, 30005 (2009).
- Sagués et al. (2007) F. Sagués, J. M. Sancho, and J. García-Ojalvo, Rev Mod Phys 79, 829 (2007).