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

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.

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 [1,[14][15][16]. 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 [1,2,[17][18][19][20][21][22][23][24][25]. Synchronization of collective spatiotemporal rhythms has been investigated experimentally in chemical systems [26,27] and may be of significant practical importance, e.g., in biomedical engineering [14][15][16]. 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 [2,24,25,[28][29][30][31][32][33] and for rotating spirals [34,35]. 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 approxi-mate 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.

A. Phase reduction of limit-cycle solutions in reaction-diffusion systems
We consider weakly perturbed RD systems exhibiting stable rhythmic dynamics, described by ∂ ∂t X(r, t) = F(X, r) + D∇ 2 X + p(r, t).
Here, the vector field X(r, t) represents the state (e.g., concentrations of chemical species) of the RD medium at point r at time t, F(X, r) represents the local reaction dynamics at r, D∇ 2 X represents the diffusion of X over the medium with a matrix D of diffusion constants, and p(r, t) represents weak spatiotemporal perturbations. Explicit dependence of F on r, such as medium heterogeneity, may exist. We assume that the RD system (1) without perturbation (p = 0) exhibits a stable rhythmic dynamics, i.e., it possesses a stable limit-cycle solution χ : X 0 (r, t) = X 0 (r, t + T ) of period T = 2π/ω, where ω denotes frequency, and that this solution persists and deforms only slightly even if the system is weakly perturbed (p = 0). 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 (p = 0). As in the ODE case (see Appendix A), we first introduce a phase θ = ωt (mod 2π) to a system state X 0 (r, t) on the limit cycle χ so thatθ(t) = ω constantly holds, and denote the system state as X 0 (r; θ) using the phase θ. To perform phase reduction, we also need to assign a phase to a system state X(r, t) 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 θ = Θ{X(r, t)} that maps X(r, t) in the basin of χ to a scalar phase θ such thatθ(t) = ω constantly holds. This leads to the notion of isochrons [1-5, 36, 37], 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 Θ{X(r)} of a system state X(r) near X 0 (r; θ) can be linearly approximated, using a certain function Q(r; θ), as around χ, where [A(r), B(r)] = A(r) · B(r)dr is the inner product between two functions. For a system state X(r) = X 0 (r; θ) on χ with phase θ, an identity Θ{X 0 (r; θ)} = θ should hold by the above definition of the phase. Moreover, for any system state X(r, t) near χ evolving under Eq. (1) with p = 0, we require that θ(t) = Θ{X(r, t)} satisfiesθ(t) = ω constantly within linear approximation. As we will derive in Appendix B, if Q(r; θ) is a periodic solution to a generalized adjoint equation with a normalization condition the functional Θ{X(r)} assumed in Eq. (2) satisfies the above requirements for the phase, and that such Q(r; θ) plays the role of the phase sensitivity function [1][2][3][4][5] for the RD system. Here, J(θ) = J(X 0 (r; θ)) is a Jacobi matrix of F estimated at X = X 0 (r; θ) on χ. Namely, we can show that the phase θ(t) = Θ{X(r, t)} of the infinite-dimensional RD system (1) approximately obeys a simple one-dimensional phase equatioṅ which is correct up to the lowest order of the perturbation p(r, t). Thus, once Q(r; θ) 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 X 0 (r; θ) with phase θ on χ is instantaneously perturbed by a weak spatial stimulus s(r), the response of the system phase after relaxation, namely, the phase response curve (PRC) [1,5], is given by within linear approximation. This PRC R(θ) can be directly measured in numerical simulations by applying impulsive perturbations to the RD system as we will illustrate in Section III.

B. 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, where X 1,2 represent the system states. We assume local and linear mutual coupling, G{X, Y} = K(Y(r, t)−X(r, t)), with a diagonal matrix K representing the intensity of the weak mutual coupling. Experimental systems like Eq. (7) have been realized by coupling a pair of photosensitive Belousov-Zhabotinsky chemical reactions via video cameras and projectors [26], and by coupling a pair of electrochemical oscillators via electrodes [27].
Denoting the phase variables of the two systems as θ 1,2 and considering the coupling term G as weak perturbations, we can approximate Eq. (7) by a pair of coupled phase equations, θ 1 (t) = ω + [Q(r; θ 1 ), G{X 0 (r; θ 1 ), X 0 (r; θ 2 )}], where X 1,2 in G are approximated by X 0 (r; θ 1,2 ) as the lowest-order approximation [1,2]. Note that the two infinitedimensional RD systems are reduced to just two one-dimensional phase equations. The coupled phase equations (8) can be analyzed in the same way as those for ordinary limit cycles [1][2][3][4][5]. Since the coupling term G is small, we can apply the averaging method to Eqs. (8), which yieldṡ where the phase-coupling function Γ is given by By subtraction, the phase difference φ = θ 1 − θ 2 obeyṡ where Γ a (φ) is a 2π-periodic, anti-symmetric function. Thus, the phase difference φ = θ 1 − θ 2 between the two RD systems approximately obeys a quite simple onedimensional equation. By examining the zeros of Γ a (φ) 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, Γ a (φ) always vanishes at φ = 0 and at φ = ±π, so that the existence of in-phase (φ = 0) and anti-phase (φ = ±π) synchronized states is assured. Their stability is determined by the slope of Γ a (φ).
FIG. 1: Circulating pulses. The system is a 1D ring of length L = 300 with periodic boundary conditions. The system parameters are α = 0, τ −1 = 0.018, γ = 1, κ = 1, and δ = 0.02. With these values, the pulse exhibits a wavy tail [19]. The oscillation period (time needed for the pulse to go around the ring) is T ≈ 1125. (a) Snapshots of the stable circulating pulse with a wavy tail, X0(x; θ = 0) = (u(x), v(x)), and the corresponding phase sensitivity function, Q(x; θ = 0) = (Qu(x), Qv(x)). (b) Phase response curves R(θ) of the circulating pulse normalized by the stimulus intensity ε. Either bell-shaped [s(x) = ε exp{−(x − 150) 2 /90}] or cosine [s(x) = ε cos(4πx/300)] perturbation is given to the activator (u) component. (c) Evolution of phase differences between two systems coupled through the u component with the coupling intensity matrix K = diag(0.001, 0) and the anti-symmetric part of the phase-coupling function Γa(φ) (rescaled by the coupling intensity 0.001). 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 [2,25,[28][29][30][31][32][33][34][35]. 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.
In numerical simulations, the size of the system is L = 80 − 600 for 1D cases and discretized using ∆x = 0.5 − 1.0 spatial grids. For 2D cases, the system size is L x × L y = 80 × 80 − 120 × 120, and discretized with ∆x = ∆y = 0.5 − 1.0 spatial grids. The explicit Euler method with a time step ∆t = 0.01 − 0.05 is used for numerical simulations of the RD system.
To numerically obtain the phase sensitivity function Q(r; θ), the adjoint equation (3) is integrated backward in time [5]. 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.

B. Circulating pulses
Our first example is a circulating-pulse solution of the FHN model with a wavy tail on a 1D ring of length L [43]. 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 X 0 (x; θ) and the corresponding phase sensitivity function Q(x; θ) for θ = 0, both propagating to the right. Results for other values of θ can simply be obtained by translating Fig. 1(a) in the x direction. It is observed that Q(x; θ) 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 Q(x; θ) 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 Q(x; θ). Figure 1(b) compares the PRCs R(θ) of the system to the weak spatial stimulus s(x) obtained by direct numerical simulations (DNS) with the theoretical results, R(θ) = [Q(x; θ), s(x)], where Q(x; θ) 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 φ = θ 1 − θ 2 from various initial conditions obtained by the DNS, and compares them with the theoretical function Γ a (φ). Reflecting the wavy shapes of X 0 and Q, Γ a (φ) 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 Γ a (φ) with negative Γ ′ a (φ). 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 [41].

C. Oscillating spots
Our second example is an oscillating spot solution of the 1D FHN system of length L with no-flux boundaries [20]. 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 X 0 (x; θ) and the phase sensitivity function Q(x; θ) for θ = 0. The activator component of Q(x; θ) 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 X 0 (x; θ) and corresponding Q(x; θ) for one oscillation period (0 ≤ θ < 2π). 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 Q(x; θ) also reflects the oscillation of the spot. Figure 2(c) shows the PRCs R(θ) to the weak stimulus s(x), 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 Γ a (φ) and evolution of the phase difference φ are plotted. For comparison, two different system sizes, L = 80 and L = 120, are used. Since the parameter α is spatially heterogeneous, the shape and oscillation period of the spot vary with L. When L = 80, in-phase synchronization (φ = 0) is linearly stable because Γ ′ a (φ = 0) < 0. In contrast, when L = 120, 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).

D. Target waves
As the third example, we consider a target wave solution [1,2,18] of the 2D FHN model on a square of side L 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 X 0 (x, y; θ) and the corresponding phase sensitivity function Q(x, y; θ) for θ = 0. As θ increases, X 0 (x, y; θ) undergoes oscillations corresponding to the emission of concentric target waves from the pacemaker, and Q(x, y; θ) oscillates accordingly. Reflecting that the pacemaker dominates overall rhythms of the system, Q(x, y; θ) is localized at the pacemaker. Figure 3(b) shows the PRCs R(θ) obtained by applying weak cosine spatial stimulus s(x, y) 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 x direction as shown in Fig. 3(e). The function Γ a (φ) has five zeros, with the in-phase (φ = 0) 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.

E. Rotating spirals
Our final example is a rotating-spiral solution [1,2,17,18,26] of the FHN model on a 2D square of side L with no-flux boundaries. Synchronization between a pair of rotating spirals was experimentally studied in [26]. 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.  3: Target waves. The system is a 2D square of side L = 100 with no-flux boundary conditions. To create a pacemaker region, the parameter α is assumed to possess localized circular heterogeneity, i.e., α(x, y) = α0 + (α1 − α0) exp(−r 4 /r 4 0 ), where r = [(x − x0) 2 + (y − y0) 2 ] 1/2 is the distance from the pacemaker center at (x0, y0) and r0 is the radius of the pacemaker region, so that α(x, y) → α1 as r → 0 and α(x, y) → α0 as r → ∞. The parameters are α0 = 0.1, α1 = −0.1, r0 = 10, and (x0, y0) = (80, 50). With these values, the system is self-oscillatory near the pacemaker center, and is excitable otherwise.  4(a) shows snapshots of the spiral solution X 0 (x, y; θ) and the corresponding phase sensitivity function Q(x, y; θ) at θ = 0. 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, Q(x, y; θ) 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, R(θ) = [Q(x, y; θ), s(x, y)], to a checkerboard-like stimulus s(x, y) 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 Γ a (φ) with five zeros, the two spirals can exhibit either in-phase (φ = 0) or anti-phase (φ = ±π) synchronization as determined by the initial conditions. Typical time sequences of the activator component measured at x = L/4, y = L/2 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 [1-3, 5, 6, 11-13, 38-40] of spatiotemporal rhythms, and will be a basis for controlling and designing spatiotemporal rhythmics in various systems.  4: (a) Rotating spirals. The system is a 2D square of side L = 120 with no-flux boundary conditions. To pin the spiral at the center, a localized circular heterogeneity of radius r0 = 4 is introduced to the parameter α(x, y) as α(x, y) = α0 + (α1 − α0) exp(−r 4 /r 4 0 ), where r = [(x − L/2) 2 + (y − L/2) 2 ] 1/2 is a distance from center of system. We assume α0 = 0.05 and α1 = 0.5, so that excitability is the highest at center. Other parameters are fixed at τ −1 = 0.005, γ = 2.5, κ = 0.15, and δ = 0. With these parameters, the oscillation period of the spiral is T = 217.37. (a) Spiral solution X0(x, y; θ) = (u(x, y), v(x, y)) and the corresponding phase sensitivity functions Q(x, y; θ) = (Qu(x, y), Qv(x, y)) at θ = 0. (b) Phase response curves R(θ) of the spiral normalized by the stimulus intensity ε. Checkerboard-like spatial perturbation is given either to the activator (u) or inhibitor (v) component, where s(x, y) = ε for x, y > L/2 or x, y < L/2, and s(x, y) = 0 otherwise. (c) Evolution of phase differences between two systems coupled through the u component with the coupling intensity matrix K = diag(2 × 10 −4 , 0), compared with the anti-symmetric part of the phase-coupling function Γa(φ) (rescaled by the coupling intensity 2 × 10 −4 ). 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.
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 finitedimensional ODEs. See Refs. [1][2][3][4][5] for details of the theory and applications to various synchronization phenomena in systems of coupled oscillators.

Geometric formulation of the phase reduction theory
We consider a limit-cycle oscillator described by the ODĖ where X(t) is a d ≥ 2 dimensional vector representing the oscillator state at time t and F determines its dynamics. Suppose that Eq. (A1) has a stable limit-cycle solution of period T , which is denoted as χ. We introduce a phase θ(t) ∈ [0, 2π) to the state X 0 (t) on χ in such a way that θ(t) increases with a constant frequency ω = 2π/T as X 0 (t) evolves along χ under Eq. (A1). This can be performed by choosing a certain state X 0 (t = 0) on χ as the origin of the phase, i.e., θ = 0, and assigning a phase value to the oscillator state X 0 (t) on χ (t ≥ 0) evolving under Eq. (A1) from the phase origin X 0 (t = 0). 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 X 0 (θ) henceforth. The above definition of the phase on χ can be extended to the whole basin of χ by assigning the same phase value θ(t) to the set of oscillator states {X(t)} that asymptotically approach the oscillator state X 0 (θ(t)) on χ under Eq. (A1), i.e., lim t→+∞ |X(t) − X 0 (θ(t))| = 0, where | · · · | represents the ordinary vector norm. This defines a phase function that maps a given oscillator state X in the basin of χ to a scalar phase θ. It is clear that the phase θ(t) = Θ(X(t)) of the state X(t) evolving under Eq. (A1) obeys a simple phase equation, 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θ where ∂Θ(X)/∂X| X=X(t) is the gradient of the phase function Θ(X) at X = X(t). Thus, the phase function Θ(X) should satisfy in the basin of χ. The set of oscillator states sharing the same phase value is called the isochron [36,37] and is the fundamental concept in the analysis of limit-cycle oscillators [1][2][3][4][5]. The whole basin of χ is foliated by such isochrons. See Fig. 5(a) for a schematic illustration of the isochrons.
Now we consider the case that the limit-cycle oscillator is weakly perturbed aṡ where the perturbation p(X, t), generally a function of the oscillator state X and time t, is assumed to be sufficiently weak so that the original limit-cycle solution χ is only slightly deformed. From Eq. (A7), the phase θ(t) = Θ(X(t)) of the perturbed oscillator obeyṡ However, this is not a closed equation for θ(t) because the gradient ∂Θ(X)/∂X| X=X(t) and the perturbation p(X(t), t) still depend on X(t). To obtain a closed equation for θ(t), X(t) in these terms are replaced by X 0 (θ(t)) at the lowest order approximation, assuming that the perturbation p(X(t), t) is sufficiently weak so that X(t) does not significantly deviate from X 0 (θ(t)) on χ. This yields an approximate closed phase equation for θ(t), which is correct up to O(|p|). Thus, by denoting p(θ, t) = p(X 0 (θ), t) and introducing a function the d-dimensional ODE (A9) describing a perturbed limit cycle can be reduced to a simple one-dimensional phase equationθ at the lowest order in the perturbation.
The key quantity for this approximation is the phase sensitivity function Z(θ) defined in Eq. (A12), which is the gradient of the isochron estimated at X = X 0 (θ) on the limit-cycle solution χ. The function Z(θ) quantifies linear phase response property of the oscillator state X 0 (θ) with phase θ on χ to infinitesimal perturbations. If X 0 (θ) is instantaneously perturbed by a weak stimulus s, the resulting phase response is given by under linear approximation. This R(θ) is called the phase response curve (PRC) of the limit-cycle oscillator described by Eq. (A1). The PRC can be obtained by applying impulsive perturbations to a limit-cycle oscillator and has been measured in various experimental systems [1].

Linear theory around the limit-cycle solution
Though we have developed a formal geometric theory by assuming the existence of the phase function Θ(X), it is generally impossible to obtain Θ(X) explicitly, except for a few simple models of limit-cycle oscillators. However, to obtain the lowest-order phase equation (A12) for weak perturbations, only the phase sensitivity function Z(θ) is actually necessary. As we show below, the function Z(θ) can be obtained as the 2π-periodic solution to the following adjoint equation [3,5]: with the constraint Z(θ) · F(X 0 (θ)) = ω, or equivalently, for 0 ≤ θ < 2π, where J(θ) = J(X 0 (θ)) is the Jacobi matrix of F(X) at X = X 0 (θ) on χ and † indicates matrix transpose.
The adjoint equation (A15) and the normalization condition (A16) can be derived in several different ways. We here use a simple argument as in [4,5] with an emphasis on the linear approximation of the isochrons near the limit cycle (see Eq. (A18) 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.
Therefore, Z(θ(t)) should satisfy the following adjoint equation: which is equivalent to Eq. (A15) by the relation d/dt = ωd/dθ (note that θ = ωt). To obtain the normalization condition Eq. (A16), we differentiate the identity θ(t) = Θ(X 0 (θ(t))) = ωt by t, which yields This gives the normalization condition Eq. (A16), again by the relation d/dt = ωd/dθ. Thus, the function Z(θ) can be obtained by solving the adjoint equation (A15) under the normalization condition Eq. (A16), and the phase function near χ is given by Eq. (A18) within linear approximation. It can also be shown that Z(θ) is the unique solution to Eq. (A15) by using the Floquet theorem characterizing the linear stability of the limit cycle, since Z(θ) is closely related to the Floquet eigenvector with the zero Floquet exponent [2][3][4][5]. In actual numerical calculations, it is useful to integrate Eq. (A15) backward in time to avoid numerical overflow, with occasional normalization using Eq. (A16) [5]. Then, by virtue of the Floquet theorem, only the functional component corresponding to Z(θ) remains numerically.
Once we obtain the frequency ω and the phase sensitivity function Z(θ), we can write down the approximate phase equation (A13) for a weakly perturbed limit-cycle oscillator described by Eq. (A9). 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 [1][2][3][4][5].
Eq. (B1). This is performed by identifying the phase with the time multiplied by the frequency. Namely, we choose a certain system state X 0 (r, t = 0) on χ as the origin of the phase, θ = 0, and assign a phase value to a state X 0 (r, t) on χ (t ≥ 0) evolving under Eq. (B1) from the phase origin X 0 (r, t = 0). We will denote the system state on χ with the phase value θ as X 0 (r; θ) henceforth. 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 θ(t) to the set of system states {X(r, t)} that eventually converge to the system state X 0 (r; θ(t)) on χ under Eq. (B1), namely, lim t→+∞ ||X(r, t) − X 0 (r; θ(t))|| = 0. (B4) Here, || · · · || denotes the L 2 norm of a spatial pattern defined as ||A(r)|| 2 = [A(r), A(r)], and the inner product between two spatial patterns A(r) and B(r) is defined as The integral is taken over the considered spatial domain with appropriate boundary conditions. This introduces a phase functional that maps a given system state X(r) in the basin of χ to a scalar phase θ. Then, the phase θ(t) = Θ{X(r, t)} of the state X(r, t) evolving under Eq. (B1) will constantly obeẏ 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 aṡ = δΘ{X(r)} δX(r) X(r)=X(r,t) , F(X(r, t), r) + D∇ 2 X(r, t) = ω, where δΘ{X(r)}/δX(r) is the functional derivative of Θ{X(r)} with respect to X(r) = X(r, t). Thus, the phase functional Θ{X(r)} should satisfy δΘ{X(r)} δX(r) , F(X(r), r) + D∇ 2 X(r) = ω (B9) 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 ∂ ∂t X(r, t) = F(X(r, t), r) + D∇ 2 X(r, t) + p{X(r, t), r, t}, where the perturbation p{X(r, t), r, t} is generally a functional of the state X(r, t), location r, and time t. We assume that the original limit-cycle solution χ is only slightly deformed by the perturbation p. If the phase functional Θ{X(r)} is given, then from Eq. (B8), the phase θ(t) = Θ{X(r, t)} of the perturbed system obeyṡ , F(X(r, t), r) + D∇ 2 X(r, t) + p{X(r, t), r, t} However, this is not a closed equation for θ(t) because the functional derivative of Θ{X(r, t)} and the perturbation p{X(r, t), r, t} still depend on X(r, t). Therefore, as in the ODE case, we approximate X(r, t) in these terms by X 0 (r; θ) on χ, assuming the perturbation p{X(r, t), r, t} to be weak enough so that the system state X(r, t) deviates from X 0 (r; θ) on χ only slightly. Then, at the lowest order approximation, a closed phase equation for θ(t) can be obtained asθ (t) ≃ ω + δΘ{X(r)} δX(r) X(r)=X0(r;θ(t)) , p{X 0 (r; θ(t)), r, t} , which is correct up to O(||p||). By denoting p(θ, r, t) = p{X 0 (r; θ), r, t} and introducing a phase sensitivity function Q(r; θ) = δΘ{X(r)} δX(r) X(r)=X0(r;θ) , the reduced phase equation (B12) can be concisely written aṡ θ(t) = ω + [Q(r; θ), p(θ, r, t)] at the lowest order in the perturbation. The function Q(r; θ) is the (functional) gradient of the isochron estimated at X(r) = X 0 (r; θ) on the limit-cycle solution χ and plays the key role in the present theory.

Linear theory around the limit-cycle solution
Though we have formally developed a geometric theory of phase reduction for the RD system, it is impossible to obtain Θ{X(r)} explicitly for general RD systems. However, only the phase sensitivity function Q(r; θ) is actually necessary to write down the lowest-order phase equation (B14) for the weakly perturbed RD systems Eq. (B10). We thus try to derive the equation for Q(r; θ) as in the ODE case, focusing only on the vicinity of the limit cycle χ.
Note that this condition can also be expressed, using again the relation ∂/∂t = ω∂/∂θ, as Q(r; θ), ∂ ∂θ X 0 (r; θ) = 1, which yields the normalization condition Eq. (4) given in Section II. Thus, if the function Q(r; θ) is the 2π-periodic solution to Eq. (B23) with the constraint (B26), the approximate phase function Eq. (B15) will satisfy the desired condition Eq. (B19) in the vicinity of χ, and Q(r; θ) will play the role of the phase sensitivity function of the limit-cycle solution χ of the RD system. Note that Eqs. (B23) and (B26) are straightforward generalizations of the conventional adjoint method for the ODE, Eq. (A15) and Eq. (A16). Generally, the function Q(r; θ) should be calculated numerically by solving the adjoint Eq. (B23) with Eq. (B26). In numerical calculations, it is useful to integrate the adjoint equation (B23) backward in time to avoid numerical overflow with occasional normalization by Eq. (B26), as we explained in Section III A.
Once we obtain the frequency ω and the phase sensitivity function Q(r; θ), we can write down the approximate phase equation (B14) for a slightly perturbed RD system described by Eq. (B10). Note that the infinite-dimensional RD system subjected to weak perturbations is reduced to a single one-dimensional phase equation, which drastically simplifies the analysis of weakly perturbed rhythmic spatiotemporal patterns. As a simple example of this phase reduction theory for RD systems, we analyzed synchronization dynamics of a pair of coupled RD systems exhibiting rhythmic patterns in Section III, i.e., the circulating pulses, oscillating spots, target waves, and rotating spirals of the FitzHugh-Nagumo model.