Phase reduction approach to synchronization of spatiotemporal rhythms in reaction-diffusion systems

Phase reduction approach to synchronization of spatiotemporal rhythms in reaction-diffusion systems

Hiroya Nakao nakao@mei.titech.ac.jp Graduate School of Information Science and Engineering, Tokyo Institute of Technology, Tokyo 152-8552, Japan    Tatsuo Yanagita Osaka Electro-Communication University, Neyagawa 572-8530, Japan    Yoji Kawamura Department of Mathematical Science and Advanced Technology, Japan Agency for Marine-Earth Science and Technology, Yokohama 236-0001, Japan
July 3, 2019
Abstract

Reaction-diffusion systems can describe a wide class of rhythmic spatiotemporal patterns observed in chemical and biological systems, such as circulating pulses on a ring, oscillating spots, target waves, and rotating spirals. These rhythmic dynamics can be considered limit cycles of reaction-diffusion systems. However, the conventional phase-reduction theory, which provides a simple unified framework for analyzing synchronization properties of limit-cycle oscillators subjected to weak forcing, has mostly been restricted to low-dimensional dynamical systems. Here, we develop a phase-reduction theory for stable limit-cycle solutions of infinite-dimensional reaction-diffusion systems. By generalizing the notion of isochrons to functional space, the phase sensitivity function — a fundamental quantity for phase reduction — is derived. For illustration, several rhythmic dynamics of the FitzHugh-Nagumo model of excitable media are considered. Nontrivial phase response properties and synchronization dynamics are revealed, reflecting their complex spatiotemporal organization. Our theory will provide a general basis for the analysis and control of spatiotemporal rhythms in various reaction-diffusion systems.

I Introduction

The phase-reduction theory provides a general framework to simplify multidimensional ordinary differential equations (ODEs) describing weakly perturbed limit-cycle oscillators to one-dimensional approximate phase equations winfree80 (); kuramoto (); hoppensteadt97 (); brown04 (); ermentrout10 (). It has drastically facilitated theoretical and experimental analysis of the synchronization properties of weakly interacting nonlinear oscillators such as chemical oscillators and spiking neurons winfree80 (); kuramoto (); hoppensteadt97 (); brown04 (); ermentrout10 (); pikovsky01 (); hudson (); taylor (); tinsley (). Methods for controlling limit-cycle oscillators have also been developed on the basis of the phase reduction theory kiss (); moehlis (); harada (); zlotnik ().

In real-world systems, rhythmic dynamics often arise collectively from a number of spatially distributed interacting elements, rather than from a single isolated oscillator, e.g., heartbeats generated by an ensemble of pulsating cardiac cells winfree80 (); glass0 (); glass1 (); bub (). Such systems are often modeled by reaction-diffusion (RD) systems, and the collective spatiotemporal rhythms are described by stable limit-cycle solutions of the RD systems winfree80 (); kuramoto (); rabinovich (); mikhailov (); hastings (); hagberg (); nomura (); yanagita (); ertl (); tyson (); maneville (). Synchronization of collective spatiotemporal rhythms has been investigated experimentally in chemical systems hildebrand01 (); fukushima () and may be of significant practical importance, e.g., in biomedical engineering glass0 (); glass1 (); bub (). In order to analyze and control the dynamics of collective spatiotemporal rhythms, it is desirable to develop a phase-reduction theory for the RD systems.

Various types of low-dimensional phase equations have been derived for RD systems, in particular for traveling pulses kuramoto (); tyson (); ei (); maneville (); cross (); mori (); ohta (); ermentrout (); lober () and for rotating spirals sandstede (); biktashev (). In most cases, however, it is assumed that the system is symmetric with respect to continuous spatial translation or rotation and the spatial structure is rigidly translating or rotating, so that their location or angle is simply identified as the phase. However, such assumptions exclude various intriguing rhythmic dynamics of RD systems that lack continuous spatial symmetry. Since limit cycles are essentially associated with temporal translational symmetry, we should be able to derive phase equations from general RD systems without recourse to spatial symmetry.

Our goal in the present study is to develop, without assuming any spatial symmetry or rigidity, a phase-reduction theory for weakly perturbed RD systems exhibiting stable rhythmic dynamics. We solve this problem by generalizing the conventional phase-reduction theory for ODEs to RD systems. Our theory gives a systematic method to approximate rhythmic dynamics of infinite-dimensional RD systems by one-dimensional phase equations, thereby facilitating detailed analysis of the synchronization dynamics of rhythmic spatiotemporal patterns. As a simple example, we analyze mutual synchronization between two interacting layers of RD systems exhibiting rhythmic dynamics. The proposed theory provides a simple, unified description of rhythmic spatiotemporal patterns and will be the basis for developing methods to control and design rhythmic spatiotemporal patterns in RD systems.

Ii Phase description of spatiotemporal rhythms

In this section, we summarize the essential results of the proposed phase reduction theory for RD systems and apply it to mutual synchronization of a pair of coupled RD systems. Full derivation of the theory will be given in Appendix B. See also Appendix A for a review of the phase reduction theory for ordinary limit-cycle oscillators described by ODEs.

ii.1 Phase reduction of limit-cycle solutions in reaction-diffusion systems

We consider weakly perturbed RD systems exhibiting stable rhythmic dynamics, described by

(1)

Here, the vector field represents the state (e.g., concentrations of chemical species) of the RD medium at point at time , represents the local reaction dynamics at , represents the diffusion of over the medium with a matrix of diffusion constants, and represents weak spatiotemporal perturbations. Explicit dependence of on , such as medium heterogeneity, may exist. We assume that the RD system (1) without perturbation () exhibits a stable rhythmic dynamics, i.e., it possesses a stable limit-cycle solution of period , where denotes frequency, and that this solution persists and deforms only slightly even if the system is weakly perturbed (). Such a limit cycle includes the circulating pulses on a ring, oscillating spots, target waves, and rotating spirals that we will analyze in Section III (see Figs. 1, 2, 3, and 4).

The purpose of the phase reduction theory is to derive a simple closed equation for the phase approximately describing limit-cycle oscillations of Eq. (1) under weak perturbation (). As in the ODE case (see Appendix A), we first introduce a phase to a system state on the limit cycle so that constantly holds, and denote the system state as using the phase . To perform phase reduction, we also need to assign a phase to a system state that is not on the limit cycle but eventually converges to , because the system state can deviate from due to perturbations. Specifically, we need a functional that maps in the basin of to a scalar phase such that constantly holds. This leads to the notion of isochrons winfree80 (); kuramoto (); hoppensteadt97 (); brown04 (); ermentrout10 (); winfree67 (); guckenheimer (), i.e., equal-phase contours of the system state around . The notion of the isochrons is at the core of the conventional phase reduction theory for ODEs and should be generalized to RD systems. It is, however, generally impossible to obtain such a functional explicitly.

To proceed, we use the assumption that the perturbation is weak and focus on the vicinity of . We make an ansatz that the phase of a system state near can be linearly approximated, using a certain function , as

(2)

around , where is the inner product between two functions. For a system state on with phase , an identity should hold by the above definition of the phase. Moreover, for any system state near evolving under Eq. (1) with , we require that satisfies constantly within linear approximation.

As we will derive in Appendix B, if is a periodic solution to a generalized adjoint equation

(3)

with a normalization condition

(4)

the functional assumed in Eq. (2) satisfies the above requirements for the phase, and that such plays the role of the phase sensitivity function winfree80 (); kuramoto (); hoppensteadt97 (); ermentrout10 (); brown04 () for the RD system. Here, is a Jacobi matrix of estimated at on .

Namely, we can show that the phase of the infinite-dimensional RD system (1) approximately obeys a simple one-dimensional phase equation

(5)

which is correct up to the lowest order of the perturbation . Thus, once is obtained from Eqs. (3) and (4), rhythmic dynamics of the RD system subjected to weak spatiotemporal perturbations, Eq. (1), can easily be analyzed using Eq. (5). This is the main result of the present study.

The phase equation (5) also shows that, if a system state with phase on is instantaneously perturbed by a weak spatial stimulus , the response of the system phase after relaxation, namely, the phase response curve (PRC) winfree80 (); ermentrout10 (), is given by

(6)

within linear approximation. This PRC can be directly measured in numerical simulations by applying impulsive perturbations to the RD system as we will illustrate in Section III.

ii.2 Mutual synchronization of spatiotemporal rhythms

As a simple example of the phase reduction approach, let us consider synchronization of a pair of weakly coupled RD systems exhibiting rhythmic dynamics,

(7)
(8)

where represent the system states. We assume local and linear mutual coupling, , with a diagonal matrix representing the intensity of the weak mutual coupling. Experimental systems like Eq. (8) have been realized by coupling a pair of photosensitive Belousov-Zhabotinsky chemical reactions via video cameras and projectors hildebrand01 (), and by coupling a pair of electrochemical oscillators via electrodes fukushima ().

Denoting the phase variables of the two systems as and considering the coupling term as weak perturbations, we can approximate Eq. (8) by a pair of coupled phase equations,

(9)
(10)

where in are approximated by as the lowest-order approximation winfree80 (); kuramoto (). Note that the two infinite-dimensional RD systems are reduced to just two one-dimensional phase equations.

The coupled phase equations (10) can be analyzed in the same way as those for ordinary limit cycles winfree80 (); kuramoto (); hoppensteadt97 (); brown04 (); ermentrout10 (). Since the coupling term is small, we can apply the averaging method to Eqs. (10), which yields

(11)

where the phase-coupling function is given by

(12)

By subtraction, the phase difference obeys

(13)

where is a -periodic, anti-symmetric function.

Thus, the phase difference between the two RD systems approximately obeys a quite simple one-dimensional equation. By examining the zeros of and their stability, we can predict the stable phase differences at which phase synchronization occurs between the two limit-cycle solutions of the coupled RD systems under the phase-reduction approximation. Since we are considering symmetrically coupled identical RD systems, always vanishes at and at , so that the existence of in-phase () and anti-phase () synchronized states is assured. Their stability is determined by the slope of .

Figure 1: Circulating pulses. The system is a 1D ring of length with periodic boundary conditions. The system parameters are , , , , and . With these values, the pulse exhibits a wavy tail hastings (). The oscillation period (time needed for the pulse to go around the ring) is . (a) Snapshots of the stable circulating pulse with a wavy tail, , and the corresponding phase sensitivity function, . (b) Phase response curves of the circulating pulse normalized by the stimulus intensity . Either bell-shaped [] or cosine [] perturbation is given to the activator () component. (c) Evolution of phase differences between two systems coupled through the component with the coupling intensity matrix and the anti-symmetric part of the phase-coupling function (rescaled by the coupling intensity ). The two pulses show multimodal phase synchronization. (d) Snapshots of phase-locked pulses with four different stable phase phase differences. Graphs (A–D) correspond to the stable phase differences shown in (c).

Iii Examples

In this section, we illustrate the phase reduction theory for RD systems by numerical simulations. We analyze phase response properties and synchronization dynamics of circulating pulses on a ring, oscillating spots, target waves, and rotating spirals of the FitzHugh-Nagumo model of excitable media (Figs. 1-4). Among these rhythmic patterns, the circulating pulses and rotating spirals are rigid and spatially symmetric, so that they may in principle be analyzed using the conventional methods kuramoto (); ei (); maneville (); cross (); mori (); ohta (); ermentrout (); lober (); sandstede (); biktashev (). In contrast, the oscillating spots and target waves are not rigid and lack translational or rotational symmetry; therefore, they cannot be treated by the conventional methods that rely on such assumptions. In any case, the phase reduction can provide a simple, unified approach to the synchronization properties of spatiotemporal rhythms. As we will see, complex spatiotemporal profiles of the rhythmic patterns can lead to interesting synchronization dynamics.

iii.1 The FitzHugh-Nagumo model

The FitzHugh-Nagumo (FHN) reaction-diffusion model is a classical model of neural spike transmission, whose dynamics is described by

(14)

where and are activator and inhibitor variables, respectively. By appropriately choosing the parameters , , , and the diffusion constants and , the FHN model can exhibit various types of rhythmic spatiotemporal dynamics yanagita (); nomura (); hastings (); hagberg (), such as the circulating pulses on a ring (Fig. 1), oscillating spots (Fig. 2), target waves (Fig. 3), and rotating spirals (Fig. 4).

In numerical simulations, the size of the system is for 1D cases and discretized using spatial grids. For 2D cases, the system size is , and discretized with spatial grids. The explicit Euler method with a time step is used for numerical simulations of the RD system.

To numerically obtain the phase sensitivity function , the adjoint equation (3) is integrated backward in time ermentrout10 (). Namely, one period of the limit-cycle oscillation is recorded by integrating the original RD system forward with sufficiently small time grids; then, the adjoint equation is spatially discretized and numerically integrated backward using the recorded time sequence of the limit cycle, with occasional normalization of the solution so that Eq. (4) is satisfied. Owing to the assumed stability of the limit-cycle solution, all modes other than the zero mode corresponding to temporal translational invariance eventually decay (the Floquet theorem), and the resultant solution gives the phase sensitivity function.

Figure 2: Oscillating spots. The system is a 1D interval of length or with no-flux boundary conditions. The parameter is space-dependent, i.e., with and , so that is the largest at the center () and the smallest at the boundaries (). Other parameter values are , , , . With these conditions, an oscillating spot constrained at the center can be generated. The oscillation period is () or (). (a) Snapshots of the oscillating spot solution and the corresponding phase sensitivity function for . (b) Evolution of and during . (c) Phase response curves of the oscillating spot normalized by the stimulus intensity . Perturbation is either bell-shaped [] or sinusoidal [] and is given to the activator () component. Results obtained by direct numerical simulations are compared with the theory, , where the inner product is taken over the 1D interval (). (d) Evolution of phase differences between two systems coupled through the component with the intensity matrix , showing in-phase synchronization for (A) and anti-phase synchronization for (B). The anti-symmetric part of the phase-coupling function (rescaled by the coupling intensity ) is shown for comparison. (e) Snapshots of activator patterns of both systems in the in-phase (A) and anti-phase (B) synchronized states.

iii.2 Circulating pulses

Our first example is a circulating-pulse solution of the FHN model with a wavy tail on a 1D ring of length  111Preliminary result on this system was partially presented in our conference proceedings (not refereed) nakao () without detailed derivation of the phase-reduction theory.. Since the pattern is rigid and the system is translationally symmetric, the phase can simply be identified as the pulse location in this case.

Figure 1(a) shows snapshots of the limit-cycle solution and the corresponding phase sensitivity function for , both propagating to the right. Results for other values of can simply be obtained by translating Fig. 1(a) in the direction. It is observed that is localized near the pulse, indicating that perturbations given only in this region can affect the phase of the pulse. It is also seen that has a wavy front, reflecting the wavy tail of the pulse. This counterintuitive result can be explained as follows. The system exhibits localized damped oscillations when it is perturbed at some spatial point. If the pulse propagates into such a region, the pulse location (i.e., the system phase) is either advanced or retarded depending on the timing of the collision, yielding the wavy front of . Figure 1(b) compares the PRCs of the system to the weak spatial stimulus obtained by direct numerical simulations (DNS) with the theoretical results, , where is obtained from the adjoint equation. The stimulus is either a bell shape localized at the center or a cosine curve, and it is given only to the activator component for the sake of simplicity. When the intensity of the stimulus is sufficiently small, good agreement is obtained.

Figure 1(c) shows the synchronization dynamics of the two RD systems, i.e., the evolution of the phase difference from various initial conditions obtained by the DNS, and compares them with the theoretical function . Reflecting the wavy shapes of and , is also wavy with many zeros, which implies the coexistence of multiple stable phase-locking points for the two-coupled circulating pulses. This is confirmed by DNS, which shows that the final phase differences are in good agreement with the zero-crossing points of with negative . Figure 1(d) shows several pairs of stably phase-locked pulses obtained by evolving the system from four different initial conditions. The two pulses synchronize where their wavy tails match, yielding multiple stable phase differences as predicted by the phase-reduction analysis. Similar multi-modal phase locking is also observed in complex oscillations of delay-differential systems kotani ().

iii.3 Oscillating spots

Our second example is an oscillating spot solution of the 1D FHN system of length with no-flux boundaries hagberg (). To pin the spot at the center, the parameter of the model is assumed to be spatially heterogeneous, namely, the excitability of the system is the largest at the center and the smallest at the boundaries. Note that the pattern is not rigid and the system lacks spatial symmetry.

Figure 2(a) shows snapshots of the limit-cycle solution and the phase sensitivity function for . The activator component of is sharply localized at both fronts of the spot; namely, the phase of the system is sensitive only to perturbations near the fronts. Figure 2(b) shows and corresponding for one oscillation period (). Perturbations given to the pulse fronts result in an advance or a delay in phase, depending on the timing, i.e., whether the spot is expanding or shrinking. The inhibitor component of also reflects the oscillation of the spot. Figure 2(c) shows the PRCs to the weak stimulus , which is either a bell shape or a cosine curve and is only applied to the activator. There is good agreement between the results of numerical simulations and theory.

The synchronization properties of a pair of oscillating spots coupled through the activator component are shown in Fig. 2(d), where the function and evolution of the phase difference are plotted. For comparison, two different system sizes, and , are used. Since the parameter is spatially heterogeneous, the shape and oscillation period of the spot vary with . When , in-phase synchronization () is linearly stable because . In contrast, when , in-phase synchronization is unstable and anti-phase synchronization () becomes stable. This prediction is confirmed by numerical simulations with various initial phase differences. Typical snapshots of the synchronized patterns are shown in Fig. 2(e).

Figure 3: Target waves. The system is a 2D square of side with no-flux boundary conditions. To create a pacemaker region, the parameter is assumed to possess localized circular heterogeneity, i.e., , where is the distance from the pacemaker center at and is the radius of the pacemaker region, so that as and as . The parameters are , , , and . With these values, the system is self-oscillatory near the pacemaker center, and is excitable otherwise. Other parameters are , , , and . The temporal oscillation period is . (a) Target wave solution and the corresponding phase sensitivity function at . (b) Phase response curves of the target wave normalized by the stimulus intensity . Sinusoidal perturbation is given either to activator () or inhibitor () component. Results obtained by direct numerical simulations are compared with the theory, , where the inner product is now taken over the 2D square (). (c) Evolution of phase differences between two systems coupled through the component with the intensity matrix , compared with the anti-symmetric part of the phase-coupling function (rescaled by the coupling intensity ). Both in-phase (A) and anti-phase (B) synchronization can occur depending on initial conditions. (d) Evolution of the activator measured at the center of systems 1 and 2 in the in-phase (A) and anti-phase (B) synchronized states. Solid line corresponds to system 1, and the dashed line corresponds to system 2. (e) Snapshots of the activator patterns of both systems in the in-phase (A) and anti-phase (B) synchronized states.

Figure 4: (a) Rotating spirals. The system is a 2D square of side with no-flux boundary conditions. To pin the spiral at the center, a localized circular heterogeneity of radius is introduced to the parameter as , where is a distance from center of system. We assume and , so that excitability is the highest at center. Other parameters are fixed at , , , and . With these parameters, the oscillation period of the spiral is . (a) Spiral solution and the corresponding phase sensitivity functions at . (b) Phase response curves of the spiral normalized by the stimulus intensity . Checkerboard-like spatial perturbation is given either to the activator () or inhibitor () component, where for or , and otherwise. (c) Evolution of phase differences between two systems coupled through the component with the coupling intensity matrix , compared with the anti-symmetric part of the phase-coupling function (rescaled by the coupling intensity ). Both in-phase synchronization (A) and anti-phase synchronization (B) can occur depending on initial conditions. Solid line corresponds to system 1, and the dashed line corresponds to system 2. (d) Snapshots of the in-phase and anti-phase synchronized states.

iii.4 Target waves

As the third example, we consider a target wave solution winfree80 (); kuramoto (); mikhailov () of the 2D FHN model on a square of side with no-flux boundaries. A circular pacemaker region is created by assuming the parameter to be heterogeneous. The system is rotationally symmetric around the pacemaker region in this case, but the target pattern is not rigid; the phase should be associated with the temporal dynamics of the pattern.

Figure 3(a) shows the limit-cycle solution and the corresponding phase sensitivity function for . As increases, undergoes oscillations corresponding to the emission of concentric target waves from the pacemaker, and oscillates accordingly. Reflecting that the pacemaker dominates overall rhythms of the system, is localized at the pacemaker. Figure 3(b) shows the PRCs obtained by applying weak cosine spatial stimulus to either the activator or inhibitor component. The numerical results are in good agreement with the theory.

Figure 3(c) shows synchronization between two target waves. Here, we consider counter-propagating target waves, i.e., one of the RD systems is inverted in the direction as shown in Fig. 3(e). The function has five zeros, with the in-phase () and anti-phase () synchronized states both being stable. Therefore, depending on initial conditions, the two target waves can exhibit both types of synchronization, as confirmed by numerical simulations. Figure 3(d) shows time sequences of the activator at the center of the two systems corresponding to the in-phase and anti-phase synchronized states, and Fig. 3(e) shows corresponding snapshots.

iii.5 Rotating spirals

Our final example is a rotating-spiral solution winfree80 (); kuramoto (); rabinovich (); mikhailov (); hildebrand01 () of the FHN model on a 2D square of side with no-flux boundaries. Synchronization between a pair of rotating spirals was experimentally studied in hildebrand01 (). Here, to pin the core of the spiral at the center of the system, circular heterogeneity in the parameter is introduced. The spiral rigidly rotates around this pinning region without changing its shape.

Figure 4(a) shows snapshots of the spiral solution and the corresponding phase sensitivity function at . Both rotate in the clockwise direction as increases. Since the system is symmetric with respect to spatial rotation around the center and the pattern is rigid, the phase simply corresponds to the rotation angle and the results for other values of can be obtained by rotating Fig. 4(a). As in the other cases, is strongly localized near the core of the spiral, indicating that the spiral tip dominates the overall phase of the system; perturbations given only to this region can affect the overall system phase. Figure 4(b) compares the PRCs obtained by DNS with the theory, , to a checkerboard-like stimulus applied either to the activator or inhibitor, showing good agreement.

Figure 4(c) shows the synchronization process between two spirals. As expected from the function with five zeros, the two spirals can exhibit either in-phase () or anti-phase () synchronization as determined by the initial conditions. Typical time sequences of the activator component measured at in the in-phase and anti-phase synchronized states are shown in Fig. 4(d), and typical snapshots of the synchronized spirals are shown in Fig. 4(e).

Iv Conclusions

We developed a phase-reduction theory for limit-cycle solutions of infinite-dimensional RD systems and illustrated its validity by analyzing mutual synchronization of a pair of RD systems exhibiting rhythmic dynamics. Our theory does not assume rigidity and spatial symmetry; therefore, it is generally applicable to a wide class of rhythmic spatiotemporal dynamics in RD systems. The theory can readily be applied, for example, to the analysis of phase locking to periodic external stimulus and noise-induced synchronization winfree80 (); kuramoto (); hoppensteadt97 (); ermentrout10 (); pikovsky01 (); moehlis (); harada (); zlotnik (); noisesync0 (); noisesync1 (); noisesync2 () of spatiotemporal rhythms, and will be a basis for controlling and designing spatiotemporal rhythmics in various systems.

References

  • (1) A. T. Winfree, The Geometry of Biological Time (Springer, 1980).
  • (2) Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer, 1984).
  • (3) F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks (Springer, 1997).
  • (4) E. Brown, J. Moehlis, and P. Holmes, On the phase reduction and response dynamics of neural oscillator populations. Neural Computation 16, 673-715 (2004).
  • (5) G. B. Ermentrout and D. H. Terman, Mathematical Foundations of Neuroscience (Springer, 2010).
  • (6) A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge University Press, 2001).
  • (7) I. Z. Kiss, Y. Zhai, and J. L. Hudson, Emerging coherence in a population of chemical oscillators, Science 296, 1676-1678 (2002).
  • (8) A. F. Taylor, M. R. Tinsley, F. Wang, Z. Huang, and K. Showalter, Dynamical quorum sensing and synchronization in large populations of chemical oscillators, Science 323, 614-617 (2009).
  • (9) M. R. Tinsley, S. Nkomo, and K. Showalter, Chimera and phase-cluster states in populations of coupled chemical oscillators, Nature Physics 8, 662-665 (2012).
  • (10) I. Z. Kiss, C. G. Rusin, H. Kori, and J. L. Hudson, Engineering complex dynamical structures: Sequential patterns and desynchronization, Science 316, 1886-1889 (2007).
  • (11) J. Moehlis, E. Shea-Brown, and H. Rabitz, Optimal inputs for phase models of spiking neurons, Journal of Computational and Nonlinear Dynamics 1, 358-367 (2006).
  • (12) T. Harada, H. A. Tanaka, M. J. Hankins, and I. Z. Kiss, Optimal waveform for the entrainment of a weakly forced oscillator, Physical Review Letters 105, 088301 (2010).
  • (13) A. Zlotnik, Y. Chen, I. Z. Kiss, H. A. Tanaka, and J. S. Li, Optimal waveform for fast entrainment of weakly forced nonlinear oscillators, Physical Review Letters 111, 024102 (2013).
  • (14) L. Glass, Synchronization and rhythmic processes in physiology, Nature 410, 277-284 (2001).
  • (15) L. Glass, L, Multistable spatiotemporal patterns of cardiac activity, PNAS 102, 10409-10410 (2005).
  • (16) G. Bub, A. Shrier, and L. Glass, Spiral wave generation in heterogeneous excitable media, Physical Review Letters 88, 058101 (2002).
  • (17) M. I. Rabinovich, A. B. Ezersky, and P. D. Weidman, The Dynamics of Patterns (World Scientific, Singapore, 2000).
  • (18) A. S. Mikhailov and K. Showalter, Control of waves, patterns and turbulence in chemical systems, Physics Reports 425, 79-194 (2006).
  • (19) S. P. Hastings, On the existence of homoclinic and periodic orbits for the Fitzhugh-Nagumo equations, The Quarterly Journal of Mathematics 27, 123-134 (1976).
  • (20) A. Hagberg and E. Meron, Pattern formation in non-gradient reaction-diffusion systems: the effects of front bifurcations, Nonlinearity 7, 805 (1994).
  • (21) T. Nomura and L. Glass, Entrainment and termination of reentrant wave propagation in a periodically stimulated ring of excitable media, Physical Review E 53, 6353 (1996).
  • (22) T. Yanagita, H. Suetani, and K. Aihara, Bifurcation analysis of solitary and synchronized pulses and formation of reentrant waves in laterally coupled excitable fibers, Physical Review E 78, 056208 (2008).
  • (23) A. S. Mikhailov and G. Ertl (Eds.), Engineering of Chemical Complexity (World Scientific Lecture Notes in Complex Systems Vol. 11). (World Scientific Publishing, 2013).
  • (24) J. J. Tyson and J. P. Keener, Singular perturbation theory of traveling waves in excitable media (a review), Physica D 32, 327-361 (1988).
  • (25) P. Manneville, Dissipative structures and weak turbulence, (Academic Press, 1990).
  • (26) M. Hildebrand, J. Cui, E. Mihaliuk, J. Wang, and K. Showalter, Synchronization of spatiotemporal patterns in locally coupled excitable media, Physical Review E 68, 026205 (2003).
  • (27) S. Fukushima, S. Nakanishi, K. Fukami, S. I. Sakai, T. Nagai, T. Tada, and Y. Nakato, Observation of synchronized spatiotemporal reaction waves in coupled electrochemical oscillations of an NDR type, Electrochemistry Communications 7, 411-415 (2005).
  • (28) S. I. Ei, The motion of weakly interacting pulses in reaction-diffusion systems, Journal of Dynamics and Differential Equations 14, 85-137 (2002).
  • (29) M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Reviews of Modern Physics 65, 851-1112 (1993).
  • (30) H. Mori and Y. Kuramoto, Dissipative structures and chaos, (Springer, 1997).
  • (31) T. Ohta, Pulse dynamics in a reaction-diffusion system, Physica D 151, 61-72 (2001).
  • (32) G. B. Ermentrout, J. Z. Jalics, and J. E. Rubin, Stimulus-driven traveling solutions in continuum neuronal models with a general smooth firing rate function, SIAM Journal on Applied Mathematics 70, 3039-3064 (2010).
  • (33) J. Löber, M. Bär, and H. Engel, Front propagation in one-dimensional spatially periodic bistable media, Physical Review E 86, 066210 (2012).
  • (34) B. Sandstede, A. Scheel, and C. Wulff, Dynamics of spiral waves on unbounded domains using center-manifold reductions, Journal of Differential Equations 141, 122-149 (1997).
  • (35) V. N. Biktashev, D. Barkley, and I. V. Biktasheva, Orbital motion of spiral waves in excitable media, Physical Review Letters 104, 058302 (2010).
  • (36) A. T. Winfree, Biological rhythms and the behavior of populations of coupled oscillators, Journal of Theoretical Biology 16, 15-42 (1967).
  • (37) J. Guckenheimer, Isochrons and phaseless sets, Journal of Mathematical Biology 1, 259-273 (1975).
  • (38) J. N. Teramae and D. Tanaka, Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators, Physical Review Letters 93, 204103 (2004).
  • (39) H. Nakao, K. Arai, and Y. Kawamura, Noise-induced synchronization and clustering in ensembles of uncoupled limit-cycle oscillators, Physical Review Letters 98, 184101 (2007).
  • (40) Y. Kawamura, H. Nakao, K. Arai, H. Kori, and Y. Kuramoto, Collective phase sensitivity, Physical Review Letters 101, 024101 (2008).
  • (41) K. Kotani, I. Yamaguchi, Y. Ogawa, Y. Jimbo, H. Nakao, and G. B. Ermentrout, Adjoint method provides phase response functions for delay-induced oscillations, Physical Review Letters 109, 044101 (2012).
  • (42) H. Nakao, T. Yanagita, and K. Kawamura, Phase description of stable limit-cycle solutions in reaction-diffusion systems, Procedia IUTAM 5, 227-233 (2012) (not refereed).

Appendix A Phase reduction of ordinary limit-cycle oscillators

In this section, we review the classical phase reduction theory for ordinary limit-cycle oscillators described by finite-dimensional ODEs. See Refs. winfree80 (); kuramoto (); hoppensteadt97 (); ermentrout10 (); brown04 () for details of the theory and applications to various synchronization phenomena in systems of coupled oscillators.

a.1 Geometric formulation of the phase reduction theory

We consider a limit-cycle oscillator described by the ODE

(15)

where is a dimensional vector representing the oscillator state at time and determines its dynamics. Suppose that Eq. (15) has a stable limit-cycle solution of period ,

(16)

which is denoted as . We introduce a phase to the state on in such a way that increases with a constant frequency as evolves along under Eq. (15). This can be performed by choosing a certain state on as the origin of the phase, i.e., , and assigning a phase value

(17)

to the oscillator state on () evolving under Eq. (15) from the phase origin . Namely, we identify the oscillator phase with the time multiplied by the frequency. We will denote the oscillator state on with the phase value as henceforth.

The above definition of the phase on can be extended to the whole basin of by assigning the same phase value to the set of oscillator states that asymptotically approach the oscillator state on under Eq. (15), i.e.,

(18)

where represents the ordinary vector norm. This defines a phase function

(19)

that maps a given oscillator state in the basin of to a scalar phase . It is clear that the phase of the state evolving under Eq. (15) obeys a simple phase equation,

(20)

not only on the limit-cycle solution but also in the whole basin of . Using the chain rule for the derivatives, it can be shown that

(21)

where is the gradient of the phase function at . Thus, the phase function should satisfy

(22)

in the basin of . The set of oscillator states sharing the same phase value is called the isochron winfree67 (); guckenheimer () and is the fundamental concept in the analysis of limit-cycle oscillators winfree80 (); kuramoto (); hoppensteadt97 (); ermentrout10 (); brown04 (). The whole basin of is foliated by such isochrons. See Fig. 5(a) for a schematic illustration of the isochrons.

Figure 5: (a) Isochrons of a limit cycle. The same phase value is assigned to the oscillator states that asymptotically converge to the same state on the limit cycle. (b) Linear approximation of the phase function near the limit cycle.

Now we consider the case that the limit-cycle oscillator is weakly perturbed as

(23)

where the perturbation , generally a function of the oscillator state and time , is assumed to be sufficiently weak so that the original limit-cycle solution is only slightly deformed. From Eq. (21), the phase of the perturbed oscillator obeys

(24)

However, this is not a closed equation for because the gradient and the perturbation still depend on . To obtain a closed equation for , in these terms are replaced by at the lowest order approximation, assuming that the perturbation is sufficiently weak so that does not significantly deviate from on . This yields an approximate closed phase equation for ,

(25)

which is correct up to . Thus, by denoting and introducing a function

(26)

the -dimensional ODE (23) describing a perturbed limit cycle can be reduced to a simple one-dimensional phase equation

(27)

at the lowest order in the perturbation.

The key quantity for this approximation is the phase sensitivity function defined in Eq. (26), which is the gradient of the isochron estimated at on the limit-cycle solution . The function quantifies linear phase response property of the oscillator state with phase on to infinitesimal perturbations. If is instantaneously perturbed by a weak stimulus , the resulting phase response is given by

(28)

under linear approximation. This is called the phase response curve (PRC) of the limit-cycle oscillator described by Eq. (15). The PRC can be obtained by applying impulsive perturbations to a limit-cycle oscillator and has been measured in various experimental systems winfree80 ().

a.2 Linear theory around the limit-cycle solution

Though we have developed a formal geometric theory by assuming the existence of the phase function , it is generally impossible to obtain explicitly, except for a few simple models of limit-cycle oscillators. However, to obtain the lowest-order phase equation (26) for weak perturbations, only the phase sensitivity function is actually necessary. As we show below, the function can be obtained as the -periodic solution to the following adjoint equation hoppensteadt97 (); ermentrout10 ():

(29)

with the constraint , or equivalently,

(30)

for , where is the Jacobi matrix of at on and indicates matrix transpose.

The adjoint equation (29) and the normalization condition (30) can be derived in several different ways. We here use a simple argument as in ermentrout10 (); brown04 () with an emphasis on the linear approximation of the isochrons near the limit cycle (see Eq. (34) below). We use the same idea to develop a phase reduction theory for the limit-cycle solutions of reaction-diffusion systems in Appendix B. See Fig. 5(b) for a schematic illustration.

We first note that, when is sufficiently small, the phase function for the ODE can be expanded in a Taylor series around as

(31)
(32)
(33)

using the phase sensitivity function in Eq. (26). Therefore, if the oscillator state is close to the oscillator state with phase on the limit-cycle solution , can be linearly approximated as

(34)

Suppose a initial state on , and a slightly perturbed initial state

(35)

near , where is a small perturbation given to . We evolve these two states without applying further perturbations. Then, from Eq. (15), the linearized equation for the small perturbation is given by

(36)

where is the Jacobi matrix of at on . From Eq. (34), the phase of the unperturbed state is given by , and that of the perturbed state can be expressed as

(37)

under the linear approximation. Note here that should also increase with a constant frequency , i.e., within the linear approximation, because no perturbation is given after . Thus, the following equation should hold for evolving from arbitrary :

(38)

Therefore, should satisfy the following adjoint equation:

(39)

which is equivalent to Eq. (29) by the relation (note that ). To obtain the normalization condition Eq. (30), we differentiate the identity by , which yields

(40)

This gives the normalization condition Eq. (30), again by the relation .

Thus, the function can be obtained by solving the adjoint equation (29) under the normalization condition Eq. (30), and the phase function near is given by Eq. (34) within linear approximation. It can also be shown that is the unique solution to Eq. (29) by using the Floquet theorem characterizing the linear stability of the limit cycle, since is closely related to the Floquet eigenvector with the zero Floquet exponent kuramoto (); ermentrout10 (); brown04 (); hoppensteadt97 (). In actual numerical calculations, it is useful to integrate Eq. (29) backward in time to avoid numerical overflow, with occasional normalization using Eq. (30ermentrout10 (). Then, by virtue of the Floquet theorem, only the functional component corresponding to remains numerically.

Once we obtain the frequency and the phase sensitivity function , we can write down the approximate phase equation (27) for a weakly perturbed limit-cycle oscillator described by Eq. (23). This approximation, called the phase reduction, greatly simplifies theoretical analysis of weakly perturbed limit cycles and has been extensively used for analyzing synchronization dynamics of weakly interacting nonlinear oscillators winfree80 (); kuramoto (); hoppensteadt97 (); ermentrout10 (); brown04 ().

Appendix B Phase reduction theory for reaction-diffusion systems

In this section, we give a full derivation of the phase reduction theory for RD systems that we briefly summarized in Section II A. Our aim is to derive a simple one-dimensional phase equation for rhythmic spatiotemporal patterns described as limit-cycle solutions of RD systems without recourse to spatial symmetry of the patterns. We do not require the patterns to be rigidly translating or rotating in the RD medium without changing their spatial profiles, as typically assumed in the conventional derivation of the phase equations for RD systems. Such rhythmic patterns can thus include oscillating spots and target waves, which vary their spatial profiles periodically. Rigidly circulating waves or rotating spirals with spatial translational or rotational symmetry are also limit-cycle solutions of RD systems, and thus they can also be treated in the same framework as we showed in Section III.

Our strategy is to generalize the conventional phase reduction theory for ordinary limit cycles described by ODEs (see Appendix A for comparison), which assumes only temporal translational symmetry of the oscillator dynamics, to limit-cycle solutions of infinite-dimensional RD systems, thereby avoiding the assumptions on spatial symmetry. We can develop the theory almost in parallel with the ODE case by noticing that the finite-dimensional vector is replaced by a vector field , and correspondingly the ordinary dot product of two vectors is replaced by the inner product of two vector fields.

b.1 Geometric formulation of the phase reduction theory

We consider a RD equation of the form

(41)

where the -dimensional vector represents the state of the RD medium at point in the -dimensional space at time , specifies local reaction dynamics at point , and represents diffusion of over the medium with a constant diffusion matrix . Explicit dependence of on , such as heterogeneity of the medium, may exist. Appropriate boundary conditions (e.g., periodic or no-flux) for the problem under consideration are introduced. We assume that Eq. (41) has a stable limit-cycle solution of period ,

(42)

which is denoted by . As in the ODE case, we first define a phase of the system state on the limit-cycle solution so that increases with a constant frequency as evolves on under Eq. (41). This is performed by identifying the phase with the time multiplied by the frequency. Namely, we choose a certain system state on as the origin of the phase, , and assign a phase value

(43)

to a state on () evolving under Eq. (41) from the phase origin . We will denote the system state on with the phase value as henceforth.

Figure 6: (a) Isochrons of a limit-cycle solution of a reaction-diffusion system. The same phase value is assigned to the system states (represented by vector fields) that asymptotically converge to the same state on the limit-cycle solution. (b) Linear approximation of the phase near the limit-cycle solution.

Next, we need to extend the definition of the phase to the basin of . As in the ODE case, we assign the same phase value to the set of system states that eventually converge to the system state on under Eq. (41), namely,

(44)

Here, denotes the norm of a spatial pattern defined as , and the inner product between two spatial patterns and is defined as

(45)

The integral is taken over the considered spatial domain with appropriate boundary conditions. This introduces a phase functional

(46)

that maps a given system state in the basin of to a scalar phase . Then, the phase of the state evolving under Eq. (41) will constantly obey

(47)

not only on the limit-cycle solution but also in the whole basin of . Using the chain rule for the functional derivatives, the above equation can be written as

(48)
(49)

where is the functional derivative of with respect to . Thus, the phase functional should satisfy

(50)

in the basin of . We call a set of system states sharing the same phase value the isochron of the RD system, generalizing the same notion for ODEs (Appendix A). The whole basin of is foliated by such isochrons. See Fig. 6(a) for a schematic illustration.

Now we consider the case that the RD system is weakly perturbed as

(51)

where the perturbation is generally a functional of the state , location , and time . We assume that the original limit-cycle solution is only slightly deformed by the perturbation . If the phase functional is given, then from Eq. (49), the phase of the perturbed system obeys

(52)
(53)

However, this is not a closed equation for because the functional derivative of and the perturbation still depend on . Therefore, as in the ODE case, we approximate in these terms by on , assuming the perturbation to be weak enough so that the system state deviates from on only slightly. Then, at the lowest order approximation, a closed phase equation for can be obtained as

(54)

which is correct up to . By denoting