Optimization of linear and nonlinear interaction schemes for stable synchronization of weakly coupled limit-cycle oscillators
OOptimization of linear and nonlinear interaction schemes for stable synchronization ofweakly coupled limit-cycle oscillators
Nobuhiro Watanabe, Yuzuru Kato, ∗ Sho Shirasaka, and Hiroya Nakao Department of Systems and Control Engineering,Tokyo Institute of Technology, Tokyo 152-8552, Japan Department of Information and Physical Sciences,Graduate School of Information Science and Technology,Osaka University 1-5 Yamadaoka, Suita, Osaka 565-0871, Japan
Optimization of mutual synchronization between a pair of limit-cycle oscillators with weak sym-metric coupling is considered in the framework of the phase reduction theory. By generalizing aprevious study [34] on the optimization of cross-diffusion coupling matrices between the oscillators,we consider optimization of mutual coupling signals to maximize the linear stability of the synchro-nized state, which are functionals of the past time sequences of the oscillator states. For the case oflinear coupling, optimization of the delay time and linear filtering of coupling signals are considered.For the case of nonlinear coupling, general drive-response coupling is considered, and the optimalresponse and driving functions are derived. The theoretical results are exemplified by numericalsimulations.
PACS numbers: 89.75.Hc, 89.75.Fb, 05.10.-a
I. INTRODUCTION
Synchronization of rhythmic dynamical elements ex-hibiting periodic oscillations is widely observed in thereal world [1–3]. Collective oscillations arising from themutual synchronization of dynamical elements often playimportant functional roles, such as the synchronized se-cretion of insulin from pancreatic beta cells [1, 4] and syn-chronized oscillation of power generators [3, 5, 6]. Clar-ifying the mechanisms of collective synchronization anddevising efficient methods of mutual synchronization arethus both fundamentally and practically important.The stable periodic dynamics of rhythmic elements areoften modeled as limit-cycle oscillators [1–3, 7]. Whenmutual interactions between limit-cycle oscillators areweak, synchronization dynamics of the oscillators canbe analyzed using the phase reduction theory [1, 8–12].In this approach, nonlinear multi-dimensional dynam-ics of an oscillator is reduced to a simple approximatephase equation, characterized by the natural frequencyand phase sensitivity of the oscillator. The phase reduc-tion theory facilitates systematic and detailed analysis ofsynchronization dynamics. It has been used to explainnontrivial synchronization dynamics of coupled oscilla-tors, such as the collective synchronization transition ofan ensemble of coupled oscillators [1, 8–12]. General-ization of the method for non-conventional limit-cyclingsystems, such as time-delayed oscillators [13, 14], hybridoscillators [15, 16], collectively oscillating networks [17],and rhythmic spatiotemporal patterns [18, 19], has alsobeen discussed.Recently, the phase reduction theory has been appliedfor the control and optimization of synchronization dy-namics in oscillatory systems. For example, Moehlis et ∗ Corresponding author: [email protected] al. [20], Harada et al. [21], Dasanayake and Li [22], Zlot-nik et al. [23–25], Pikovsky [26], Tanaka et al. [27–29],Wilson et al. [30], Pyragas et al. [31], and Monga etal. [32, 33] have used the phase reduction theory (as wellas the phase-amplitude reduction theory) to derive opti-mal driving signals for the stable entrainment of nonlin-ear oscillators in various physical situations.In a similar spirit, in our previous study [34], we con-sidered a problem of improving the linear stability of syn-chronized state between a pair of limit-cycle oscillatorsby optimizing a cross-diffusion coupling matrix betweenthe oscillators, where different components of the oscilla-tors are allowed to interact. We also considered a pair ofmutually interacting reaction-diffusion systems exhibit-ing rhythmic spatiotemporal patterns and derived opti-mal spatial filters for stable mutual synchronization [35].In this study, we consider this problem in a more gen-eral setting, whereby the oscillators can interact not onlyby their present states but also through the time se-quences of their past states. We first consider linearcoupling with time delay or temporal filtering and de-rive the optimal delay time or linear filter. We thenconsider general nonlinear coupling with a mutual drive-response configuration and derive the optimal responsefunction and driving function for stable synchronization.We argue that, although we consider general couplingthat can depend on the past time sequences of the oscil-lators, the optimal mutual coupling can be obtained asa function of the present phase values of the oscillatorsin the framework of the phase-reduction approximation.The results are illustrated by numerical simulations us-ing Stuart-Landau and FitzHugh-Nagumo oscillators asexamples.This paper is organized as follows. In Sec. II, we in-troduce a general model of coupled limit-cycle oscillatorsand reduce it to coupled phase equations. In Sec. III,we consider the case with linear coupling and derive the a r X i v : . [ n li n . AO ] J a n optimal time delay and optimal linear filter for couplingsignals. In Sec. IV, we consider nonlinear coupling ofthe drive-response type and derive the optimal responsefunction and driving function. In Sec. V, a summary isprovided. II. MODELA. Pair of weakly coupled oscillators
In this study, we consider a pair of weakly and sym-metrically coupled limit-cycle oscillators with identicalproperties, where the oscillators can mutually interactnot only through their present states but also throughtheir past time sequences. We assume that the oscilla-tors are generally described by the following equations:˙ X ( t ) = F ( X ( t )) + (cid:15) ˆ H { X ( t )1 ( · ) , X ( t )2 ( · ) } , ˙ X ( t ) = F ( X ( t )) + (cid:15) ˆ H { X ( t )2 ( · ) , X ( t )1 ( · ) } , (1)where X , ∈ R N are N -dimensional state vectors of theoscillators 1 and 2 at time t , F : R N → R N is a suf-ficiently smooth vector field representing the dynamicsof individual oscillators, and (cid:15) ˆ H represents weak mutualcoupling between the oscillators. Here, ˆ H : C × C → R N ( C is a function space of the time sequences of length T )is a functional of two vector functions; i.e., the past timesequences of X , ( t ), and 0 < (cid:15) (cid:28) X = F ( X ), hasan exponentially stable limit cycle ˜ X ( t ) = ˜ X ( t + T ) ofperiod T and frequency ω = 2 π/T , and the system statedeviates only slightly from this limit cycle if weak pertur-bations are applied. The time sequences X ( t ) i ( · ) ( i = 1 , H represent functions X i ( τ )on the interval [ t − T, t ] and are defined as X ( t ) i ( τ ) = X i ( t + τ ) ( − T ≤ τ ≤ . (2)In this study, we fix the length of the time sequences usedfor the coupling as the natural period T of the limit-cycleoscillator.Under the assumptions stated above, we can employthe standard method of phase reduction [1, 8–12] andintroduce a phase function Θ( X ) : R N → [0 , π ), whichsatisfies ˙Θ( X ) = F ( X ) · ∇ Θ( X ) = ω in the basin ofthe limit cycle. Using Θ( X ), the phase variable of thisoscillator can be defined as θ = Θ( X ), which constantlyincreases with time as ˙ θ = ω , both on and in the basinof the limit cycle (2 π is identified with 0 as usual). Theoscillator state on the limit cycle can be represented asa function of θ as X ( θ ) = ˜ X ( t = θ/ω ), which is a 2 π -periodic function of θ , X ( θ ) = X ( θ + 2 π ). Similarly,to represent a time sequence on the limit cycle, we alsointroduce the notation X ( θ )0 ( τ ) = X ( θ + ωτ ) ( − T ≤ τ ≤
0) (3) and abbreviate this as X ( θ )0 ( · ). We hereafter use thesenotations to represent the system states and their timesequences on the limit cycle. B. Phase reduction and averaging
Assuming that perturbations applied to the oscillatorare sufficiently weak, the state vector of the oscillator canbe represented approximately by X ( t ) ≈ X ( θ ( t )) as afunction of the phase θ ( t ). Then, the phase sensitivityfunction (PSF), defined by Z ( θ ) = ∇ Θ( X ) | X = X ( θ ) :[0 , π ) → R N , characterizes the linear response prop-erty of the oscillator phase to weak perturbations. Whenthe oscillator is weakly driven by a perturbation, p , as˙ X = F ( X ) + (cid:15) p , the reduced phase equation is givenby ˙ θ = ω + (cid:15) Z ( θ ) · p up to O ( (cid:15) ). The PSF can becalculated as a 2 π -periodic solution to an adjoint equa-tion ωd Z ( θ ) /dθ = − J † ( X ( θ )) Z ( θ ), with a normaliza-tion condition Z ( θ ) · d X ( θ ) /dθ = 1, where J( X ) : R N → R N × N is a Jacobian matrix of F at X = X ( θ ) and † denotes transpose. In the numerical analysis, Z ( θ ) canbe calculated easily using the backward time-evolution ofthe adjoint equation as proposed by Ermentrout [10].Defining the phase values of the oscillators 1 , θ , = Θ( X , ), the following pair of equations can bederived from Eq. (1):˙ θ ( t ) = ω + (cid:15) Z ( θ ( t )) · ˆ H { X ( t )1 ( · ) , X ( t )2 ( · ) } , ˙ θ ( t ) = ω + (cid:15) Z ( θ ( t )) · ˆ H { X ( t )2 ( · ) , X ( t )1 ( · ) } , (4)which are correct up to O ( (cid:15) ). Next, we use the fact thatthe deviation of each oscillator state from the limit cycleis small and of O ( (cid:15) ): X ( t )1 , ( τ ) = X ( θ , ( t ))0 ( τ ) + O ( (cid:15) ) ( − T ≤ τ ≤ . (5)Substituting into Eq. (4) and ignoring errors of O ( (cid:15) ), weobtain a pair of reduced phase equations,˙ θ = ω + (cid:15) Z ( θ ) · ˆ H { X ( θ )0 ( · ) , X ( θ )0 ( · ) } , ˙ θ = ω + (cid:15) Z ( θ ) · ˆ H { X ( θ )0 ( · ) , X ( θ )0 ( · ) } , (6)which are also correct up to O ( (cid:15) ).The coupling term, ˆ H { X ( θ )0 ( · ) , X ( θ )0 ( · ) } , is a func-tional of the two time sequences X ( θ )0 ( · ) and X ( θ )0 ( · ).However, because we can neglect the deviations of theoscillator states from the limit cycle at the lowest orderapproximation, the functionals of X ( θ )0 ( · ) and X ( θ )0 ( · )are actually determined solely by the two phase values θ and θ . Therefore, we can regard the coupling termas an ordinary function of θ and θ , defined by H ( θ , θ ) = ˆ H { X ( θ )0 ( · ) , X ( θ )0 ( · ) } , (7)and rewrite the phase equations as˙ θ = ω + (cid:15) Z ( θ ) · H ( θ , θ ) , ˙ θ = ω + (cid:15) Z ( θ ) · H ( θ , θ ) . (8)Thus, though we started from Eq. (1) with general cou-pling functionals that depend on the past time sequencesof the oscillators, the coupled system reduces to a pair ofsimple ordinary differential equations that depend onlyon the present phase values θ and θ of the oscillatorswithin the phase-reduction approximation.Following the standard averaging procedure [8, 9], weintroduce slow phase variables φ , = θ , − ωt , rewritethe equations as˙ φ = (cid:15) Z ( φ + ωt ) · H ( φ + ωt, φ + ωt ) , ˙ φ = (cid:15) Z ( φ + ωt ) · H ( φ + ωt, φ + ωt ) , (9)and average the small right-hand side of these equationsover one-period of oscillation. This yields the followingaveraged phase equations, which are correct up to O ( (cid:15) ):˙ θ = ω + (cid:15) Γ( θ − θ ) , ˙ θ = ω + (cid:15) Γ( θ − θ ) , (10)where Γ( φ ) is the phase coupling function defined asΓ( φ ) = 12 π (cid:90) π Z ( φ + ψ ) · H ( φ + ψ, ψ ) dψ = (cid:68) Z ( φ + ψ ) · H ( φ + ψ, ψ ) (cid:69) ψ = (cid:68) Z ( ψ ) · H ( ψ, ψ − φ ) (cid:69) ψ . (11)Here, we have defined an average of a function f ( ψ ) overone period of oscillation as (cid:68) f ( ψ ) (cid:69) ψ = 12 π (cid:90) π f ( ψ ) dψ. (12) C. Linear stability of the in-phase synchronizedstate
From the coupled phase equations (10), the dynamicsof the phase difference φ = θ − θ obeys˙ φ = (cid:15) [Γ( φ ) − Γ( − φ )] , (13)where the right-hand side is the antisymmetric part ofthe phase coupling function Γ( φ ). This equation alwayshas a fixed point at φ = 0 corresponding to the in-phasesynchronized state. In a small vicinity of φ = 0, theabove equation can be linearized as˙ φ ≈ (cid:15) Γ (cid:48) (0) φ. (14)The derivative of the phase coupling function is given byΓ (cid:48) ( φ ) = − (cid:68) Z ( ψ ) · H (cid:48) ( ψ, ψ − φ ) (cid:69) ψ , (15)where the partial derivative of H with respect to thesecond argument is denoted as H (cid:48) ( ψ , ψ ) = ∂ H ( ψ , ψ ) ∂ψ . (16) Thus, the linear stability of this state is characterized bythe exponent 2 (cid:15) Γ (cid:48) (0), and a larger − Γ (cid:48) (0) yields a higherlinear stability of the in-phase synchronized state.In this study, we consider optimization of the mutualcoupling term, H , or either the parameters or functionsincluded in it, so that the linear stability − Γ (cid:48) (0) = (cid:68) Z ( ψ ) · H (cid:48) ( ψ, ψ ) (cid:69) ψ (17)of the in-phase synchronized state, φ = 0, is maximizedunder appropriate constraints on the intensity of the mu-tual coupling. (a)(b) FIG. 1. Limit-cycle orbit and phase sensitivity function of theFitzHugh-Nagumo oscillator. Time sequences of the x and y components are plotted for one period of oscillation, 0 ≤ θ < π . (a) Limit cycle ( x ( θ ) , y ( θ )). (b) Phase sensitivityfunction ( Z x ( θ ) , Z y ( θ )). D. Examples of limit-cycle oscillators
The Stuart-Landau (SL) oscillator is used in the fol-lowing numerical examples, which is a normal form of thesupercritical Hopf bifurcation and is described by X = (cid:18) xy (cid:19) ∈ R , (18) F = (cid:18) F x F y (cid:19) = (cid:18) x − ay − (cid:0) x + y (cid:1) ( x − by ) ax + y − (cid:0) x + y (cid:1) ( bx + y ) (cid:19) , (19)where the parameters are fixed at a = 2 and b = 1. Thisoscillator has a stable limit cycle with a natural frequency ω = a − b = 1 and period T = 2 π , represented by X ( θ ) = (cid:18) x ( θ ) y ( θ ) (cid:19) = (cid:18) cos θ sin θ (cid:19) , (0 ≤ θ < π ) . (20)The phase function can be explicitly represented byΘ( x, y ) = arctan yx − b x + y ) (21)on the whole xy -plane except (0 , Z ( θ ) = (cid:18) Z x Z y (cid:19) = (cid:18) − sin θ − b cos θ cos θ − b sin θ (cid:19) . (22)As another example, we use the FitzHugh-Nagumo(FHN) oscillator, described by X = (cid:18) xy (cid:19) ∈ R , (23) F = (cid:18) F x F y (cid:19) = (cid:18) x ( x − c )(1 − x ) − yµ − ( x − dy ) (cid:19) , (24)where the parameters are fixed at c = − . d = 0 . µ = 100. As µ is large, this oscillator is a slow-fast system whose x variable evolves much faster thanthe y variable, leading to relaxation oscillations. Withthese parameters, the natural period of the oscillation is T ≈ . ω ≈ . X ( θ ) = ( x ( θ ) , y ( θ )) † and phase functionΘ( X ) cannot be obtained analytically for this model, butthe PSF Z ( θ ) can be obtained by numerically solving theadjoint equation. Figure 1 shows the time sequences ofthe limit-cycle orbit X ( θ ) and PSF Z ( θ ) for one periodof oscillation, 0 ≤ θ < π . III. LINEAR COUPLINGA. Time-delayed coupling
First, we consider a simple time-delayed coupling,where each oscillator is driven by the past state of theother oscillator. The model is given by˙ X = F ( X ) + (cid:15) √ P K X ( t − τ ) , ˙ X = F ( X ) + (cid:15) √ P K X ( t − τ ) , (25)where 0 < (cid:15) (cid:28) P > ∈ R N × N isa constant matrix specifying which components of theoscillator states X , ( t ) are coupled, and τ (0 ≤ τ ≤ T )is a time delay. In our previous study [34], we consideredoptimization of the matrix K for the case where the twooscillators are nearly identical and coupled without timedelay. Here, we consider two oscillators with identicalproperties, fix the matrix K constant, and vary the timedelay, τ , to improve the linear stability of the in-phasesynchronized state.In this case, the coupling functionals are given byˆ H { X ( t )1 , X ( t )2 } = √ P K X ( t − τ ) , ˆ H { X ( t )2 , X ( t )1 } = √ P K X ( t − τ ) , (26) which can be expressed as functions of θ and θ as H ( θ , θ ) = √ P K X ( θ − ωτ ) , H ( θ , θ ) = √ P K X ( θ − ωτ ) , (27)after phase reduction. The phase coupling function isΓ( φ ) = (cid:68) Z ( ψ ) · H ( ψ, ψ − φ ) (cid:69) ψ = (cid:68) Z ( ψ ) · √ P K X ( ψ − φ − ωτ ) (cid:69) ψ , (28)and the linear stability is characterized by − Γ (cid:48) (0) = (cid:68) Z ( ψ ) · H (cid:48) ( ψ, ψ ) (cid:69) ψ = (cid:68) Z ( ψ ) · √ P K X (cid:48) ( ψ − ωτ ) (cid:69) ψ , (29)where X (cid:48) ( θ ) = d X ( θ ) /dθ .The maximum stability is attained only when τ satis-fies ∂∂τ {− Γ (cid:48) (0) } = − ω (cid:68) Z ( ψ ) · √ P K X (cid:48)(cid:48) ( ψ − ωτ ) (cid:69) ψ = 0 , (30)where X (cid:48)(cid:48) ( θ ) = d X ( θ ) /dθ . We denote the value of τ satisfying the above equation as τ ∗ , i.e., (cid:68) Z ( ψ ) · √ P K X (cid:48)(cid:48) ( ψ − ωτ ∗ ) (cid:69) ψ = 0 . (31)By partial integration using the 2 π -periodicity of Z ( θ )and X ( θ ), this can also be expressed as (cid:68) Z (cid:48) ( ψ ) · √ P K X (cid:48) ( ψ − ωτ ∗ ) (cid:69) ψ = 0 , (32)which has the form of a cross-correlation function be-tween Z (cid:48) ( θ ) and √ P K X (cid:48) ( θ ). Because both of these func-tions are 2 π -periodic with zero-mean, the left-hand sideof Eq. (32) is also 2 π -periodic with zero mean. Thus,there are at least two values of τ satisfying the aboveequation, as long as Z (cid:48) ( θ ) and √ P K X (cid:48) ( θ ) are non-constant functions (which holds generally for ordinarylimit cycles). By choosing an appropriate value of τ ∗ ,the maximum stability is given by − Γ (cid:48) (0) = (cid:114) P (cid:68) Z ( ψ ) · K X (cid:48) ( ψ − ωτ ∗ ) (cid:69) ψ . (33) B. Coupling via linear filtering
Generalizing the time-delayed coupling, we consider acase in which the past time sequences of both oscillatorstates are linearly filtered and used as driving signals forthe other oscillators. The model is given by˙ X = F ( X ) + (cid:15) (cid:90) T h ( τ )K X ( t − τ ) dτ, ˙ X = F ( X ) + (cid:15) (cid:90) T h ( τ )K X ( t − τ ) dτ, (34)where h ( τ ) : [0 , T ] → R is a T -periodic real scalar func-tion representing a linear filter, and K ∈ R N × N is aconstant matrix specifying which components of X arecoupled. We optimize the linear filter h ( τ ) for a givencoupling matrix K under a constraint specified below.The coupling functionals are given byˆ H { X ( t )1 , X ( t )2 } = (cid:90) T h ( τ )K X ( t − τ ) dτ, ˆ H { X ( t )2 , X ( t )1 } = (cid:90) T h ( τ )K X ( t − τ ) dτ, (35)which simplify to ordinary functions H ( θ , θ ) = (cid:90) T h ( τ )K X ( θ − ωτ ) dτ, H ( θ , θ ) = (cid:90) T h ( τ )K X ( θ − ωτ ) dτ, (36)after phase reduction. The phase coupling function isgiven byΓ( φ ) = (cid:68) Z ( φ + ψ ) · (cid:90) T h ( τ )K X ( ψ − ωτ ) dτ (cid:69) ψ = (cid:68) (cid:90) T Z ( ψ ) · h ( τ )K X ( ψ − ωτ − φ ) dτ (cid:69) ψ (37)and the linear stability of the in-phase synchronized stateis characterized by − Γ (cid:48) (0) = (cid:68) (cid:90) T Z ( ψ ) · h ( τ )K X (cid:48) ( ψ − ωτ ) dτ (cid:69) ψ . (38)We constrain the L -norm (cid:107) h ( τ ) (cid:107) = (cid:113)(cid:82) T h ( τ ) dτ ofthe linear filter, h ( τ ), as (cid:107) h ( τ ) (cid:107) = Q , where Q > h ( τ ) that maximizes the linear stability, − Γ (cid:48) (0). That is,we consider an optimization problem:maximize − Γ (cid:48) (0) subject to (cid:107) h ( τ ) (cid:107) = Q. (39)To this end, we define an objective functional as S { h, λ } = − Γ (cid:48) (0) + λ ( (cid:107) h ( ψ ) (cid:107) − Q )= (cid:68) (cid:90) T Z ( ψ ) · h ( τ )K X (cid:48) ( − ωτ + ψ ) dτ (cid:69) ψ + λ (cid:32)(cid:90) T h ( τ ) dτ − Q (cid:33) , (40)where λ is a Lagrange multiplier. From the extremumcondition of S , the functional derivative of S with respectto h ( τ ) should satisfy δSδh ( τ ) = (cid:68) Z ( ψ ) · K X (cid:48) ( − ωτ + ψ ) (cid:69) ψ + 2 λh ( τ ) = 0 (41) and the partial derivative of S by λ should satisfy ∂S∂λ = (cid:90) T h ( τ ) dτ − Q = 0 . (42)Thus, the optimal linear filter h ( τ ) is given by h ( τ ) = − λ (cid:68) Z ( ψ ) · K X (cid:48) ( ψ − ωτ ) (cid:69) ψ . (43)The Lagrange multiplier λ is determined from the con-straint (cid:107) h ( τ ) (cid:107) = Q , i.e.,14 λ (cid:90) T (cid:68) Z ( ψ ) · K X (cid:48) ( ψ − ωτ ) (cid:69) ψ dτ = Q, (44)as λ = − (cid:115) Q (cid:90) T (cid:68) Z ( ψ ) · K X (cid:48) ( ψ − ωτ ) (cid:69) ψ dτ , (45)where the negative sign should be chosen for the in-phasesynchronized state to be linearly stable, − Γ (cid:48) (0) >
0. Themaximum linear stability with the optimized h ( τ ) is − Γ (cid:48) (0) = (cid:115) Q (cid:90) T (cid:68) Z ( ψ ) · K X (cid:48) ( ψ − ωτ ) (cid:69) ψ dτ . (46) C. Numerical examples
1. Time-delayed coupling
We use the SL and FHN oscillators in the followingnumerical illustrations. For both models, the couplingmatrix is assumed to beK = (cid:18) (cid:19) . (47)We compare the optimized case with the non-optimizedcase, i.e., ˙ X = F ( X ) + (cid:15) √ P K X ( t ) , ˙ X = F ( X ) + (cid:15) √ P K X ( t ) , (48)where (cid:15) is a small parameter that determines the couplingstrength and P control the norm of the coupling signal.The mean square of the coupling term over one-period ofoscillation is the same in both cases, that is, (cid:68) |√ P K X ( ψ ) | (cid:69) ψ = (cid:68) |√ P K X ( ψ − ωτ ∗ ) | (cid:69) ψ . (49)First, for the SL oscillator, we can analytically calcu-late the optimal time delay. The linear stability of thein-phase synchronized state, Eq. (29), is given by − Γ (cid:48) (0) = √ P (cid:68) Z x ( ψ ) x (cid:48) ( ψ − ωτ ) (cid:69) ψ = √ P ωτ ) − b sin( ωτ )] . (50) no delayoptimal no delay (phase)no delay (DNS)optimal (phase)optimal (DNS) (b)(a) FIG. 2. Synchronization of two Stuart-Landau oscillators cou-pled with time delay. The results with the optimal time delayare compared with those without time delay. (a) Evolution ofthe difference ∆ x between x variables of the two oscillators.(b) Evolution of the phase difference φ between the oscillators. (a) (b) no delayoptimal (c) no delayoptimal (d) no delay(phase)no delay(DNS)optimal(phase)optimal(DNS) FIG. 3. Synchronization of two FitzHugh-Nagumo oscilla-tors coupled with time delay. In (b–d), the results with theoptimal time delay are compared with those without time de-lay. (a) Linear stability − Γ (cid:48) (0) and its derivative − ∂ Γ (cid:48) (0) /∂τ vs. time delay τ . The crosses indicate the values of τ where − ∂ Γ (cid:48) (0) /∂τ = 0. (b) Antisymmetric part of the phase cou-pling function, Γ( φ ) − Γ( − φ ). (c) Evolution of the difference∆ x between x variables of the two oscillators. (d) Evolutionof the phase difference φ between the oscillators. The optimal time delay τ = τ ∗ is determined fromEq. (30), or equivalently from (cid:68) Z x ( ψ ) x (cid:48)(cid:48) ( ψ − ωτ ) (cid:69) ψ = b cos( ωτ ) + sin( ωτ )2 = 0 . (51)For the parameter values b = 1 and ω = 1, this equationis satisfied when τ = 3 π/ τ = 7 π/
4. Substitutingthis into Eq. (50), we find that τ ∗ = 7 π/ − Γ (cid:48) (0) = √ P / √
2. For the case with no time delay,the linear stability is − Γ (cid:48) (0) = √ P /
2. Thus, by ap-propriately choosing the time delay, the linear stabilityimproves by a factor of √ (cid:15) = 0 . P = 1, and the initial phase difference is φ (0) = π/
4. In Fig. 2(a), the difference ∆ x between the x variables of the two oscillators, obtainedby direct numerical simulations of the coupled SL oscil-lators, is plotted as a function of t . It can be seen thatthe in-phase synchronized state is established faster inthe optimized case because of the higher linear stability.Figure 2(b) shows the convergence of the phase difference φ to 0. It can be seen that the results of the reducedphase equation agree well with those of direct numericalsimulations.Figure 3 shows the results for two FHN oscillators,where (cid:15) = 0 . P = 1, and the initial phase differ-ence is φ = π/
4. Figure 3(a) plots the linear stability − Γ (cid:48) (0) and its derivative − ∂ Γ (cid:48) (0) /∂τ as functions of thetime delay τ , where there are two extrema of − Γ (cid:48) (0). Wechoose the larger extremum, which is attained at the op-timal time delay τ ∗ ≈ .
6. The antisymmetric part ofthe phase coupling function, Γ( φ ) − Γ( − φ ), is shown inFig 3(b) for the cases with the optimal delay and withoutdelay. It can be seen that the stability of the in-phasesynchronized state φ = 0 is improved, as indicated by thestraight lines in Fig. 3(b), where − Γ (cid:48) (0) ≈ .
654 with theoptimized time delay and − Γ (cid:48) (0) ≈ .
221 without thetime delay. The evolution of the difference ∆ x betweenthe x variables of the oscillators is plotted as a function of t in Fig. 3(c). The phase differences φ converging toward0, obtained from the phase equation and direct numericalsimulations of the original model, are shown in Fig. 3(d).It can be seen that the convergence to the in-phase syn-chronization is faster with the optimized time delay, andthe results of the reduced phase equation agree well withdirect numerical simulations.
2. Coupling via linear filtering
We also assume that the coupling matrix K is givenby Eq. (47), and compare the results for the optimizedcase with linear filtering with those for the non-filteredcase given by Eq. (48). We choose the parameter Q thatconstrains the norm of the linear filter so that the squaredaverage of the coupling term over one-period of oscillationbecomes equal to that in the non-filtered case given byEq. (48), i.e., (cid:68) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) T h ( τ )K X ( ψ − ωτ ) dτ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:69) ψ = (cid:68) |√ P K X ( ψ ) | (cid:69) ψ . (52)For the SL oscillators, the optimal filter h ( τ ), Eq. (43),is explicitly calculated as h ( τ ) = (cid:115) Qωπ (1 + b ) [cos( ωτ ) − b sin( ωτ )] . (53)The optimal phase coupling function, Eq. (37), and op-timized linear stability, Eq. (46), are expressed asΓ( φ ) = − (cid:114) π (1 + b ) Qω sin φ, (54)and − Γ (cid:48) (0) = 12 (cid:114) π (1 + b ) Qω , (55)respectively. We take Q = ωP/π so that Eq. (52) is satis-fied. The linear stability is then − Γ (cid:48) (0) = (cid:112) (1 + b ) P / − Γ (cid:48) (0) = √ P / √ b = 1.It is important to note that, in the SL oscillator case, X ( ψ ), Z ( ψ ), and hence the linear filter h ( τ ) containonly the fundamental frequency, i.e., they are purely si-nusoidal. Thus, the linear filtering can only shift thephase of the coupling signal and gives the same result asthe previous case with the simple time delay. It is alsointeresting to note that the stability cannot be improved(it is already optimal without filtering) when the param-eter b , which characterizes non-isochronicity of the limitcycle, is zero.Figure 4 shows the synchronization of two SL oscilla-tors, with and without linear filtering, where (cid:15) = 0 . P = 1, and the initial phase difference is φ = π/
4. Fig-ure 4(a) shows evolution of the difference ∆ x between the x variables of the oscillators, and Fig. 4(b) shows the con-vergence of the phase difference φ to 0. We can see thatthe in-phase synchronized state is established faster inthe optimized case, and the results of the reduced phasemodel and direct numerical simulations agree well.For the FHN oscillators, the optimal linear filter canbe calculated from the time sequences of the limit-cyclesolution and PSF obtained numerically. Figure 5 showsthe synchronization of two coupled FHN oscillators, withand without linear filtering, where (cid:15) = 0 . P = 1, Q ≈ . φ = π/ φ ) − Γ( − φ ) of the phase coupling function Γ( φ ),(c) evolution of the difference ∆ x between the x vari-ables of the oscillators, and (d) convergence of the phasedifference φ toward 0. The linear stability is given by − Γ (cid:48) (0) ≈ .
844 for the case with the optimal filter andby − Γ (cid:48) (0) ≈ .
221 for the case without filtering, as shownby the straight lines in Fig. 5(b). The in-phase synchro-nized state is established faster in the optimized case,and the results of the reduced phase model and directnumerical simulations agree well. Because the FHN oscil-lator has the higher harmonic components in X ( ψ ) and Z ( ψ ), the optimal filter h ( τ ) can exploit these compo-nents, and hence the improvement in the linear stabilityis larger than that for the case with simple delay. no filteroptimal no filter (phase)no filter (DNS)optimal (phase)optimal (DNS) (a) (b) FIG. 4. Synchronization of two Stuart-Landau oscillators cou-pled with linear filtering. The results with the optimal filter-ing are compared with those without filtering. (a) Evolutionof the difference ∆ x in the x variables between the oscillators.(b) Evolution of the phase difference φ between the oscillators. (a) (b) no filteroptimal (c) no filteroptimal (d) no filter(phase)no filter(DNS)optimal(phase)optimal(DNS) FIG. 5. Synchronization of two FitzHugh-Nagumo oscillatorscoupled with linear filtering. In (b–d), the results with theoptimal filtering are compared with those without filtering.(a) Optimal linear filter h ( τ ). (b) Antisymmetric part of thephase coupling function Γ( φ ) − Γ( − φ ). (c) Evolution of thedifference ∆ x in the x variables between the oscillators. (d)Evolution of the phase difference φ between the oscillators. IV. NONLINEAR COUPLINGA. Mutual drive-response coupling
In this section, we consider the case of oscillators in-teracting through nonlinear coupling functionals. We as-sume that the coupling is of a drive-response type, i.e., itcan be written as a product of a response matrix of thedriven oscillator and a driving function that transformsthe signal from the other oscillator. The model is givenby ˙ X ( t ) = F ( X ( t )) + (cid:15) ˆA { X ( t )1 ( · ) } ˆ G { X ( t )2 ( · ) } , ˙ X ( t ) = F ( X ( t )) + (cid:15) ˆA { X ( t )2 ( · ) } ˆ G { X ( t )1 ( · ) } , (56)where the matrix ˆA : C → R N × N is a functional ofthe time sequence of each oscillator representing its re-sponse properties and ˆ G : C → R N is a functional thattransforms the time sequence of the other oscillator toa driving signal. It should be noted that we may alsoinclude self-coupling terms of the form (cid:15) ˆ I { X ( t ) ( · ) } , toeach equation, which allows the analysis, for example,of diffusive coupling that depends on the state differencebetween the oscillators. As explained in the AppendixA, it can be shown that inclusion of such self-couplingterms does not alter the results and the linear stabilityremains the same in the phase-reduction approximation.We thus analyze Eq. (56) hereafter.The coupling functionals in this case are given byˆ H { X ( t )1 , X ( t )2 } = ˆA { X ( t )1 ( · ) } ˆ G { X ( t )2 ( · ) } , ˆ H { X ( t )2 , X ( t )1 } = ˆA { X ( t )2 ( · ) } ˆ G { X ( t )1 ( · ) } , (57)and, as argued in Sec. IIB, at the lowest-order phase re-duction, these functionals can be expressed as ordinaryfunctions of the phase θ and θ as H ( θ , θ ) = A( θ ) G ( θ ) , H ( θ , θ ) = A( θ ) G ( θ ) , (58)where we introduced ordinary 2 π -periodic functions Aand G of θ and θ . Using these functions, the phasecoupling function is given byΓ( φ ) = (cid:68) Z ( ψ ) · A( ψ ) G ( ψ − φ ) (cid:69) ψ , (59)and the linear stability is characterized by − Γ (cid:48) (0) = (cid:68) Z ( ψ ) · A( ψ ) G (cid:48) ( ψ ) (cid:69) ψ = (cid:68) A † ( ψ ) Z ( ψ ) · G (cid:48) ( ψ ) (cid:69) ψ = − (cid:68) ddψ (cid:2) A † ( ψ ) Z ( ψ ) (cid:3) · G ( ψ ) (cid:69) ψ , (60)where the last expression is obtained by partial integra-tion using 2 π -periodicity of A( ψ ), Z ( ψ ), and G ( ψ ).Therefore, although we started from Eq. (56) with ageneral drive-response coupling that depends on the pasttime sequences of the oscillators, the linear stability canbe represented only by the present phase values of theoscillators at the lowest-order phase reduction. In thefollowing subsections, we consider the optimization of theresponse matrix A( ψ ) or the driving function G ( ψ ), rep-resented as functions of the phase ψ . B. Optimal response matrix
As for the first case, we optimize the response matrixA( ψ ) as a function of the phase ψ , assuming that thedriving functional ˆ G is given. We introduce a constraintthat the squared Frobenius norm of A( ψ ) averaged overone period of oscillation is fixed as (cid:68) (cid:107) A( ψ ) (cid:107) (cid:69) ψ = P , and consider an optimization problem:maximize − Γ (cid:48) (0) subject to (cid:68) (cid:107) A( ψ ) (cid:107) (cid:69) ψ = P, (61)where (cid:107) A (cid:107) = (cid:113)(cid:80) i,j A ij represents the Frobenius normof the matrix A = ( A ij ). By defining an objective func-tional, S { A , λ } = − Γ (cid:48) (0) + λ (cid:18)(cid:68) (cid:107) A( ψ ) (cid:107) (cid:69) ψ − P (cid:19) = (cid:68) Z ( ψ ) · A( ψ ) G (cid:48) ( ψ ) (cid:69) ψ + λ (cid:18)(cid:68) (cid:107) A( ψ ) (cid:107) (cid:69) ψ − P (cid:19) , (62)where λ is a Lagrange multiplier, and by taking the func-tional derivative with respect to each component, A ij , ofA, we obtain the extremum condition. In this case, δSδA ij ( ψ ) = 12 π Z i ( ψ ) G (cid:48) j ( ψ ) + λπ A ij ( ψ ) = 0 , (63)and we obtain A ij ( ψ ) = − λ Z i ( ψ ) G (cid:48) j ( ψ ) (64)i.e., A( ψ ) = − λ Z ( ψ ) G (cid:48) ( ψ ) † , (65)and the Lagrange multiplier is determined from the con-straint, (cid:68) (cid:107) A( ψ ) (cid:107) (cid:69) ψ = 14 λ (cid:68) (cid:107) Z ( ψ ) G (cid:48) ( ψ ) † (cid:107) (cid:69) ψ = P (66)as λ = − (cid:114) P (cid:68) (cid:107) Z ( ψ ) G (cid:48) ( ψ ) † (cid:107) (cid:69) ψ , (67)where the negative sign is chosen so that − Γ (cid:48) (0) >
0. Themaximum stability of the in-phase synchronized state is − Γ (cid:48) (0) = (cid:114) P (cid:68) (cid:107) Z ( ψ ) G (cid:48) ( ψ ) † (cid:107) (cid:69) ψ . (68) C. Optimal driving functional
We can also seek the function G ( ψ ) that provides theoptimal driving signal as a function of the phase ψ , as-suming that the response matrix ˆA is given. We constrainthe squared average of G ( ψ ) over one period of oscillationas (cid:68) | G ( ψ ) | (cid:69) ψ = P , and maximize the linear stability ofthe in-phase state:maximize − Γ (cid:48) (0) subject to (cid:68) | G ( ψ ) | (cid:69) ψ = P. (69)We define an objective functional, S { A , λ } = − Γ (cid:48) (0) + λ (cid:18)(cid:68) | G ( ψ ) | (cid:69) ψ − P (cid:19) = − (cid:68) ddψ (cid:2) A † ( ψ ) Z ( ψ ) (cid:3) · G ( ψ ) (cid:69) ψ + λ (cid:18)(cid:68) | G ( ψ ) | (cid:69) ψ − P (cid:19) , (70)where λ is a Lagrange multiplier. From the extremumcondition for S , we obtain δSδ G ( ψ ) = − π ddψ [A † ( ψ ) Z ( ψ )] + λπ G ( ψ ) = 0 (71)and the constraint on G . The optimal driving functionis given by G ( ψ ) = 12 λ ddψ [A † ( ψ ) Z ( ψ )] , (72)where the Lagrange multiplier λ should be chosen to sat-isfy the norm constraint,14 λ (cid:68) (cid:12)(cid:12)(cid:12)(cid:12) ddψ [A † ( ψ ) Z ( ψ )] (cid:12)(cid:12)(cid:12)(cid:12) (cid:69) ψ = P. (73)This yields λ = − (cid:115) P (cid:68) (cid:12)(cid:12)(cid:12)(cid:12) ddψ [A † ( ψ ) Z ( ψ )] (cid:12)(cid:12)(cid:12)(cid:12) (cid:69) ψ , (74)where the negative sign is taken to satisfy Γ (cid:48) (0) <
0. Themaximum stability is − Γ (cid:48) (0) = (cid:115) P (cid:68) (cid:12)(cid:12)(cid:12)(cid:12) ddψ [A † ( ψ ) Z ( ψ )] (cid:12)(cid:12)(cid:12)(cid:12) (cid:69) ψ . (75) D. Numerical examples
1. Optimal response matrix
As an example, we assume that the driving func-tional ˆ G { X ( t ) ( · ) } is simply given by ˆ G { X ( t ) ( · ) } = X ( t ),and seek the optimal response matrix A( ψ ) satisfying (cid:68) (cid:107) A( ψ ) (cid:107) (cid:69) ψ = P . For comparison, we also consider anidentity response matrix, A I = diag( (cid:112) P/ , (cid:112) P/ (cid:68) (cid:107) A I (cid:107) (cid:69) ψ = P . Note that both the x and y components are coupled, in contrast to the previ-ous section where only the x component is coupled. For the SL oscillator, the optimal response matrix canbe analytically expressed asA( ψ ) = (cid:114) P b × (cid:18) sin ψ ( b cos ψ + sin ψ ) − cos ψ ( b cos ψ + sin ψ )sin ψ ( b sin ψ − cos ψ ) cos ψ (cos ψ − b sin ψ ) (cid:19) , (76)and the phase coupling function is given by Γ( φ ) = − (cid:112) (1 + b ) P sin φ , which gives the optimal linear sta-bility − Γ (cid:48) (0) = (cid:112) (1 + b ) P . In contrast, for theidentity matrix A I , the phase coupling function isΓ( φ ) = − (cid:112) P/ b cos φ + sin φ ) and the linear stabilityis − Γ (cid:48) (0) = √ P / √
2. Thus, the linear stability is im-proved by a factor of (cid:112) b ). identityoptimal identity (phase)identity (DNS)optimal (phase)optimal (DNS) (a) (b) FIG. 6. Synchronization of Stuart-Landau oscillators withthe optimal response matrix. (a) Evolution of the difference∆ x between the x variables of the oscillators. (b) Evolutionof the phase difference φ between the oscillators. (a) (b) identityoptimal (c) identityoptimal (d) identity(phase)identity(DNS)optimal(phase)optimal(DNS) FIG. 7. Synchronization of FitzHugh-Nagumo oscillators withthe optimal response matrix. In (b–d), the results with theoptimal response matrix are compared with those with theidentity response matrix. (a) Four components of the opti-mal response matrix A( ψ ). (b) Antisymmetric part of thephase coupling function, Γ( φ ) − Γ( − φ ). (c) Evolution of thedifference ∆ x between the x variables of the oscillators. (d)Evolution of the phase difference φ between the oscillators. Figure 6 shows synchronization of two SL oscillatorsfor the cases with the optimal response matrix A( ψ ) and0with the identity matrix A I ( ψ ), where b = 1, (cid:15) = 0 . P = 2, and the initial phase difference is φ = π/
4. Fig-ure 6(a) shows the evolution of the difference ∆ x in the x variables between the two oscillators, and Fig. 6(b) showsthe convergence of the phase difference φ to 0. The in-phase synchronized state is more quickly established inthe optimized case, and the results of the reduced phasemodel and direct numerical simulations agree well.For the FHN oscillator, the optimal response matrixcan be calculated numerically. Figure 7 compares thesynchronization dynamics of two coupled FHN oscillatorswith the optimal and identity response matrices, where (cid:15) = 0 . P = 2, and the initial phase difference is φ = π/
4. Figure 7(a) shows four components of the opti-mal response matrix A( ψ ) for 0 ≤ ψ < π . It is notablethat the magnitude of A ( ψ ) is much larger than theother components, indicating that driving the y compo-nent of each oscillator by using the x component of theother oscillator is efficient in synchronizing the oscillatorsin this case. Figure 7(b) plots the antisymmetric part ofthe phase coupling functions for the optimal and iden-tity response matrices, which shows that a much higherstability is attained in the optimized case ( − Γ (cid:48) (0) ≈ . − Γ (cid:48) (0) ≈ . x between the two oscillators, and Fig. 7(d) shows theconvergence of the phase difference φ to 0. In order to usethe optimal response matrix, instantaneous phase valuesof the oscillators are necessary. In the direct numeri-cal simulations shown here, we approximately evaluatedthe phase value by linearly interpolating two consecutivecrossings times of the oscillator state at an appropriatePoincar´e section, and this value was used to generate thedriving signal. It can be seen from the figures that thein-phase synchronized state is established much faster inthe optimized case, and the results of the reduced phasemodel and direct numerical simulations agree well.
2. Optimal driving functional
For the numerical simulations, we assume that A( ψ )is simply given by an identity matrix, diag(1 , G ( ψ ) is then simply given as G ( ψ ) ∝ Z (cid:48) ( ψ ) from Eq. (72), with the norm constraint (cid:68) | G ( ψ ) | (cid:69) ψ = P . For the SL oscillator, the optimaldriving function is explicitly given by G ( ψ ) = (cid:114) P b (cid:18) cos ψ − b sin ψb cos ψ + sin ψ (cid:19) . (77)Figure 8 shows synchronization of two SL oscillatorscoupled through the optimal driving function, and cou-pled without transformation of the oscillator state, i.e.,ˆ G { X ( t ) ( · ) } = X ( t ), where b = 1, (cid:15) = 0 .
01, and P = 1,and the initial phase difference is φ = π/
4. Figure 8(a)shows the evolution of the difference ∆ x between the x no transformoptimal no transform (phase)no transform (DNS)optimal (phase)optimal (DNS) (a) (b) FIG. 8. Synchronization of Stuart-Landau oscillators coupledwith the optimal driving function. (a) Evolution of the dif-ference ∆ x between the x variables of the oscillators. (b)Evolution of the phase difference φ between the oscillators. variables of the two oscillators, and Fig. 8(b) shows theconvergence of the phase difference φ to 0. It is confirmedthat the linear stability of the in-phase synchronized stateis improved in the optimized case, and the results of thereduced phase model and direct numerical simulationsagree well.For the FHN oscillator, the norm of X ( ψ ) is (cid:68) | X ( ψ ) | (cid:69) ψ ≈ . P of G ( ψ )to this value. The optimal driving function can be calcu-lated from X ( ψ ) and Z ( ψ ) obtained numerically. Fig-ure 9 shows synchronization of two FHN oscillators cou-pled with the optimal driving function, as well as com-parison with the non-transformed case, where (cid:15) = 0 . P ≈ . φ = π/ G ( ψ ) for0 ≤ ψ < π , which is proportional to the derivative Z (cid:48) ( ψ ). Figure 9(b) plots the antisymmetric part of thephase coupling function for the optimal driving function G ( ψ ) with that without transformation, respectively, in-dicating a much higher linear stability in the optimizedcase ( − Γ (cid:48) (0) ≈ . − Γ (cid:48) (0) ≈ .
999 without optimization).Figure 9(c) shows a plot of the evolution of the dif-ference ∆ x between x variables of the oscillators, andFig. 9(d) shows the convergence of the phase difference φ to 0. Similar to the previous case with the optimal re-sponse matrix, instantaneous phase values of the oscilla-tors are approximately evaluated by linear interpolationand used to generate the optimal driving signal in thedirect numerical simulations. We can confirm that thein-phase synchronized state is established much faster inthe optimized case, and the results of the reduced phasemodel and direct numerical simulations agree well. V. DISCUSSION
We have shown that, by optimizing the mutual cou-pling between coupled oscillators, the linear stabilityof the in-phase synchronized state can be improvedand faster convergence to the synchronization can beachieved. We have shown that, even if we start from asystem of coupled oscillators with general coupling func-1 (a) (b) no transformoptimal (c) no transformoptimal (d) no transform(phase)no transform(DNS)optimal(phase)optimal(DNS)
FIG. 9. Synchronization of FitzHugh-Nagumo oscillators cou-pled with the optimal driving function. In (b–d), the re-sults with the optimal driving function are compared withthose without transformation. (a) Optimal driving function G ( ψ ) = ( G ( ψ ) , G ( ψ )). (b) Antisymmetric part of the phasecoupling function, Γ( φ ) − Γ( − φ ). (c) Evolution of the dif-ference ∆ x between the x variables of the oscillators. (d)Evolution of the phase difference φ . tionals that depend on the past time sequences of theoscillators, the system can be approximately reduced toa pair of simple ordinary differential equations that de-pend only on the present phase values of the oscillatorswithin the phase reduction theory, and the optimal cou-pling function between the oscillators can be obtained asa function of the phase values. Though we have consid-ered only the simplest cases where two oscillators withidentical properties are symmetrically coupled withoutnoise, the theory can also be extended to include hetero-geneity of the oscillators or noise.The linear coupling with time delay or linear filteringcan be realized without measuring the phase values of theoscillator, once the correlation functions of the PSF andthe limit-cycle orbit (or their derivatives) are obtained.The nonlinear coupling requires the measurement of thephase values of the oscillators, but can further improvethe linear stability of the synchronized state. We haveshown that a simple approximate evaluation of the phasevalues by a linear interpolation gives reasonable resultseven though it may yield a somewhat incorrect evaluationof the true phase values.It is interesting to compare the present analysis forstable synchronization between the two oscillators withthe optimization of driving signals for injection lockingof a single oscillator, which has been analyzed by Zlotnik et al. [24] and others (briefly explained in the AppendixB for a simple case). In Sec. IV C, we have obtained theoptimal driving functional. In particular, when A( ψ ) =K, where K is a constant matrix, the optimal driving signal is G ( ψ ) = 12 λ K † Z (cid:48) ( ψ ) (78)and the maximum stability is − Γ (cid:48) (0) = (cid:114) P (cid:68) | K † Z (cid:48) ( ψ ) | (cid:69) ψ . (79)This result coincides with the optimal injection signal forstable synchronization of a single oscillator, obtained byZlotnik et al. [24]. Thus, the optimal coupling betweenthe oscillators is realized by measuring the present phase ψ of the other oscillator and applying a driving signalthat is proportional to K † Z (cid:48) ( ψ ) to the oscillator.It is also interesting to note that we have obtainedsimilar expressions for the maximum stability in all ex-amples, − Γ (cid:48) (0) = (cid:114) P (cid:68) · · · (cid:69) ψ , where · · · depends on thequantity to be optimized. This is because we are essen-tially maximizing the inner product of the PSF with thederivative of the driving signal under a mean-square con-straint on the parameters or functions included in thedriving signal in all cases.The linear coupling schemes in Sec. III would beeasy to realize experimentally. The nonlinear couplingschemes in Sec. IV require the evaluation of the phase val-ues from the oscillators, but can yield an even higher lin-ear stability. These methods may be useful when higherstability of the in-phase synchronized state between oscil-lators is desirable in technical applications. It would alsobe interesting to study interactions between rhythmic el-ements e.g. in biological systems from the viewpoint ofsynchronization efficiency. ACKNOWLEDGMENTS
This work is financially supported by JSPS KAKENHIGrant Numbers JP16K13847, JP17H03279, 18K03471,JP18H032, and 18H06478.
Appendix A: Models with self coupling
We here show that the inclusion of self-coupling termsin the model does not alter the linear stability of thein-phase synchronized state under the phase-reductionapproximation. Suppose that we have additional self-coupling terms in the model as˙ X = F ( X ) + (cid:15) [ ˆ H { X ( t )1 ( · ) , X ( t )2 ( · ) } + ˆ I { X ( t )1 ( · ) } ] , ˙ X = F ( X ) + (cid:15) [ ˆ H { X ( t )2 ( · ) , X ( t )1 ( · ) } + ˆ I { X ( t )2 ( · ) } ] , (A1)where ˆ I { X ( t ) ( · ) } is a functional representing self cou-pling. A typical example is coupled oscillators with linear2diffusive coupling,˙ X = F ( X ) + (cid:15) ( X − X ) , ˙ X = F ( X ) + (cid:15) ( X − X ) , (A2)where we may take ˆ H { X ( t )1 ( · ) , X ( t )2 ( · ) } = X ( t ) andˆ I { X ( t )1 ( · ) } = − X ( t ). By phase reduction, we obtainthe phase equations˙ θ = ω + (cid:15) Z ( θ ) · [ H ( θ , θ ) + I ( θ )] , ˙ θ = ω + (cid:15) Z ( θ ) · [ H ( θ , θ ) + I ( θ )] , (A3)and the phase coupling function˜Γ( φ ) = (cid:68) Z ( ψ ) · [ H ( ψ, ψ − φ ) + I ( ψ )] (cid:69) ψ = Γ( φ ) + (cid:68) Z ( ψ ) · I ( ψ ) (cid:69) ψ , (A4)where Γ( φ ) is the phase coupling function for the casewithout the self-coupling term and the second term is aconstant. Thus, this model also has the in-phase syn-chronized state as a fixed point and its linear stability isequal to the case without the self-coupling term, − (cid:15) Γ (cid:48) (0) = − (cid:15) ˜Γ (cid:48) (0) . (A5) Appendix B: Optimal signal for injection locking
We here briefly review Zlotnik et al. ’s result [24] on theoptimal driving signal for injection locking for a simplecase. We consider a limit-cycle oscillator driven by aperiodic driving signal whose period is the same as thenatural period T of the oscillator,˙ X = F ( X ) + (cid:15) K f ( t ) , f ( t ) = f ( t + T ) , (B1)where X is the oscillator state, F ( X ) represents its dy-namics, and (cid:15) K f ( t ) is a weak periodic driving signal,where (cid:15) is a small positive parameter and a constant ma-trix K ∈ R N × N represents which components of X aredriven by f ( t ).By phase reduction, we obtain a phase equation˙ θ = ω + Z ( θ ) · K f ( t ) (B2)for the oscillator phase θ , where Z ( θ ) is the PSF. Defining θ − ωt = φ and averaging over one oscillation period yields˙ φ = Γ( φ ) . (B3) The phase coupling function Γ( ψ ) is expressed asΓ( φ ) = 12 π (cid:90) π Z ( φ + ψ ) · K f ( ψ/ω ) dψ = 12 π (cid:90) π Z ( φ + ψ ) · K ˜ f ( ψ ) dψ = (cid:68) Z ( φ + ψ ) · K ˜ f ( ψ ) (cid:69) ψ (B4)where we have defined ˜ f ( ψ ) = f ( ψ/ω ).By choosing the origin of the phase of the periodicsignal so that Γ(0) = 0 holds, the oscillator synchronizeswith the periodic at φ = 0, and the linear stability of thissynchronized state is given byΓ (cid:48) (0) = (cid:68) Z (cid:48) ( ψ ) · K ˜ f ( ψ ) (cid:69) ψ = (cid:68) K † Z (cid:48) ( ψ ) · ˜ f ( ψ ) (cid:69) ψ . (B5)We constrain the one-period average of ˜ f ( ψ ) as (cid:68) | ˜ f ( ψ ) | (cid:69) ψ = P, (B6)and consider an objective function S { ˜ f , λ } = − Γ (cid:48) (0) + λ ( | f | − P )= − (cid:68) K † Z (cid:48) ( ψ ) · ˜ f ( ψ ) (cid:69) ψ + λ (cid:18)(cid:68) | ˜ f ( ψ ) | (cid:69) ψ − P (cid:19) , (B7)where λ is a Lagrange multiplier. Extremum conditionsare given by δSδ ˜ f ( ψ ) = − π K † Z (cid:48) ( ψ ) + λπ ˜ f ( ψ ) = 0 , (B8) ∂S∂λ = (cid:68) | ˜ f ( ψ ) | (cid:69) ψ − P = 0 . (B9)The optimal driving signal is given by˜ f ( ψ ) = 12 λ K † Z (cid:48) ( ψ ) (B10)and the constraint is14 λ (cid:68) | K † Z (cid:48) ( ψ ) | (cid:69) ψ = P, (B11)which yields λ = (cid:114) P (cid:68) | K † Z (cid:48) ( ψ ) | (cid:69) ψ , (B12)where the negative sign should be taken in order thatΓ (cid:48) (0) <
0. Thus, the optimal driving signal is given by˜ f ( ψ ) = − (cid:118)(cid:117)(cid:117)(cid:116) P (cid:68) | K † Z (cid:48) ( ψ ) | (cid:69) ψ K † Z (cid:48) ( ψ ) (B13)and the maximum linear stability is given by − Γ (cid:48) (0) = − (cid:68) K † Z (cid:48) ( ψ ) · ˜ f ( ψ ) (cid:69) ψ = (cid:114) P (cid:68) | K † Z (cid:48) ( ψ ) | (cid:69) ψ . (B14)3 [1] A. T. Winfree, The Geometry of Biological Time (Springer, New York, 1980; Springer, Second Edition,New York, 2001).[2] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchroniza-tion: A Universal Concept in Nonlinear Sciences (Cam-bridge University Press, Cambridge, 2001).[3] S. H. Strogatz,
Sync: How Order Emerges from Chaos inthe Universe, Nature, and Daily Life (Hyperion Books,New York, 2003).[4] A. Sherman and J. Rinzel, Model for synchronization ofpancreatic beta-cells by gap junction coupling, Biophys.J. , 547 (1991).[5] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa,Spontaneous synchrony in power-grid networks, NaturePhys. , 191 (2013).[6] F. D¨orfler, F. Chertkov, and F. Bullo, Synchronization incomplex oscillator networks and smart grids, PNAS ,20052010 (2013).[7] S. H. Strogatz, Nonlinear Dynamics and Chaos, 2nd ed. (Westview Press, Boulder, 2015).[8] Y. Kuramoto,
Chemical Oscillations, Waves, and Tur-bulence (Springer, New York, 1984; Dover, New York,2003).[9] F. C. Hoppensteadt and E. M. Izhikevich,
Weakly Con-nected Neural Networks (Springer, New York, 1997).[10] G. B. Ermentrout and D. H. Terman,
MathematicalFoundations of Neuroscience (Springer, New York, 2010).[11] H. Nakao, Phase reduction approach to synchronizationof nonlinear oscillators, Contemp. Phys. , 188 (2016).[12] P. Ashwin, S. Coombes, and R. Nicks, Mathematicalframeworks for oscillatory network dynamics in neuro-science, J. Math. Neurosci. , 1 (2016).[13] V. Noviˇcenko and K. Pyragas, Phase reduction of weaklyperturbed limit cycle oscillations in time-delay systems,Physica D , 10901098 (2012).[14] K. Kotani, I. Yamaguchi, Y. Ogawa, Y. Jimbo, H. Nakao,and G. B. Ermentrout, Adjoint method provides phaseresponse functions for delay-induced oscillations, Phys.Rev. Lett. , 044101 (2012).[15] S. Shirasaka, W. Kurebayashi, and H. Nakao, Phase re-duction theory for hybrid nonlinear oscillators Phys. Rev.E , 012212 (2017).[16] Y. Park, K. M. Shaw, H. J. Chiel, and P. J. Thomas,The Infinitesimal Phase Response Curves of Oscillators inPiecewise Smooth Dynamical Systems, European Jour-nal of Applied Mathematics , 90540 (2018).[17] Y. Kawamura, H. Nakao, K. Arai, H. Kori, and Y. Ku-ramoto, Collective phase sensitivity, Phys. Rev. Lett. , 024101 (2008).[18] Y. Kawamura and H. Nakao, Collective phase descriptionof oscillatory convection, Chaos , 043129 (2013).[19] H. Nakao, T. Yanagita, and Y. Kawamura, Phase-reduction approach to synchronization of spatiotemporalrhythms in reaction-diffusion systems, Phys. Rev. X ,021032 (2014). [20] J. Moehlis, E. Shea-Brown, and H. Rabitz, Optimal in-puts for phase models of spiking neurons, J. Comput.Nonlin. Dyn. , 358 (2006).[21] T. Harada, H.-A. Tanaka, M. J. Hankins, and I. Z. Kiss,Optimal waveform for the entrainment of a weakly forcedoscillator, Phys. Rev. Lett. , 088301 (2010).[22] I. Dasanayake and J.-S. Li, Optimal design of minimum-power stimuli for phase models of neuron oscillators,Phys. Rev. E , 061916 (2011).[23] A. Zlotnik and J.-S. Li, Optimal entrainment of neuraloscillator ensembles, J. Neural Eng. , 046015 (2012).[24] A. Zlotnik, Y. Chen, I. Kiss, H.-A. Tanaka, and J.-S. Li, Optimal Waveform for Fast Entrainment of WeaklyForced Nonlinear Oscillators, Phys. Rev. Lett. ,024102 (2013).[25] A. Zlotnik, R. Nagao, I. Z. Kiss, and J.-S. Li, Phase-selective entrainment of nonlinear oscillator ensembles,Nature Comm. , 10788 (2016).[26] A. Pikovsky, Maximizing coherence of oscillations by ex-ternal locking, Phys. Rev. Lett. , 070602 (2015).[27] H.-A. Tanaka, Synchronization limit of weakly forcednonlinear oscillators, J. Phys. A: Math. Theor. ,402002 (2014).[28] H.-A. Tanaka, Optimal entrainment with smooth, pulse,and square signals in weakly forced nonlinear oscillators,Physica D , 1 (2014).[29] H.-A. Tanaka, I. Nishikawa, J. Kurths, Y. Chen, andI. Z. Kiss, Optimal synchronization of oscillatory chem-ical reactions with complex pulse, square, and smoothwaveforms signals maximizes Tsallis entropy, Europhys.Lett. , 50007 (2015).[30] D. Wilson, A. B. Holt, T. I. Netoff, and J. Moehlis, Op-timal entrainment of heterogeneous noisy neurons, Fron-tiers in neuroscience , 192 (2015).[31] K. Pyragas, A. P. Fedaraviˇcius, T. Pyragien´e, and P.A. Tass, Optimal waveform for entrainment of a spikingneuron with minimum stimulating charge, Phys. Rev. E98