Feedback-induced desynchronization and oscillation quenching in a population of globally coupled oscillators
FFeedback-induced desynchronization and oscillation quenchingin a population of globally coupled oscillators
Ayumi Ozawa ∗ and Hiroshi Kori Department of Complexity Science and Engineering,The University of Tokyo, Chiba 277-8561, Japan (Dated: February 9, 2021)Motivated from a wide range of applications, various methods to control synchronizationin coupled oscillators have been proposed. Previous studies have demonstrated that globalfeedback typically induces three macroscopic behaviors: synchronization, desynchronization,and oscillation quenching. However, analyzing all of these transitions within a single the-oretical framework is difficult, and thus the feedback effect is only partially understood ineach framework. Herein, we analyze a model of globally coupled phase oscillators exposedto global feedback, which shows all of the typical macroscopic dynamical states. Analyticaltractability of the model enables us to obtain detailed phase diagrams where transitions andbistabilities between different macroscopic states are identified. Additionally, we proposestrategies to steer the oscillators into targeted states with minimal feedback strength. Ourstudy provides a useful overview of the effect of global feedback and is expected to serve asa benchmark when more sophisticated feedback needs to be designed.
I. INTRODUCTION
Synchronization is a self-organization phenomenon that occurs in interacting oscillatory ele-ments and is widely observed[1–5]. The resultant coherent oscillation is desirable in some systems.For example, synchronization is essential for the normal functioning of power grids [6]. Other bene-ficial effects of synchronization include the enhanced precision in biological oscillatory systems[4, 7–10], coordinated locomotion of animals and robots[11], and reduced congestion in models of trafficflow[12, 13].However, synchronization may also cause problems. At the Millennium Footbridge in London,the steps of the pedestrians were synchronized, and considerable lateral movement of the bridgewas observed[14]. Synchronization is also associated with neurological disorders. The local field po-tentials recorded from the brains of Parkinsonian patients and model animals often display markedoscillations, which are considered to be reflections of coherent neuronal activities[15]. Although ∗ [email protected] a r X i v : . [ n li n . AO ] F e b the mechanism of Parkinson’s disease is not well understood yet, exaggerated synchronization isone of the possible factors that induces related symptoms[16]. For some types of Parkinsonianpatients, deep brain stimulation (DBS), which involves electrical stimulation to particular regionsof the brain, may suppress the pathological collective oscillation and motor symptoms [15, 17].However, this treatment is sometimes accompanied by negative side effects[18, 19]. Thus, milderways of stimulation need to be developed.The wide range of desirable and undesirable synchronization phenomena has drawn substantialresearch attention to the control of synchronization. As a control implementation, global feed-back loops are known to be effective. In [20], the authors experimentally demonstrated that asynchronously oscillating state can be stabilized in a surface chemical reaction that inherentlyexhibits turbulent oscillatory dynamics. Global feedback may also desynchronize oscillator as-semblies. Several types of mean-field feedback that efficiently desynchronize oscillators have beenproposed [21–29]. Moreover, other behaviors may be realized through global feedback. It is provedthat in a particular class of phase oscillators, the oscillation death state, in which the mean-fieldoscillation terminates, can be stabilized via global mean-field feedback [27]. Oscillation death isknown to appear in various coupled oscillator systems [30].These extensive studies [20–29] elucidated that global feedback typically stabilizes the syn-chronous, asynchronous, or oscillation-death states. However, in the oscillator models and feed-back forms considered thus far, analytically treating all of these macroscopic states is difficult,and hence, a comprehensive phase diagram has not been obtained. With only partial knowledgebeing available on the phase diagram, feedback control may fail to realize a desired state owing tothe unexpected stabilization of other states. Detailed phase diagrams of an analytically tractablemodel would provide an insight into a general principle of feedback control and help design andtune the feedback scheme.Herein, we consider the Sakaguchi–Kuramoto model, which describes a population of noniden-tical phase oscillators with global coupling, as a coupled-oscillator model and incorporate a globalmean-field feedback loop into the system. We show that the system exhibits all three typicalmacroscopic states: the asynchronously oscillating state, the synchronously oscillating state, andthe oscillation-death state. By invoking the Ott–Antonsen ansatz[31, 32], we comprehensively per-form the existence and stability analysis of the macroscopic states, thus obtaining detailed phasediagrams in the space of feedback parameters for different coupling strengths. Our analysis eluci-dates (i) the dependency of the feedback effect on the parameters of the oscillator model; (ii) theoptimal feedback parameters for stabilizing the asynchronous state with minimal feedback strength;and (iii) the existence of the bistability between the oscillation-death state and one of the other twostates. The obtained phase diagrams are numerically validated. In addition, we propose a strategyfor realizing the asynchronous state in the bistable region. Finally, to support the robustness of ourresults, we numerically investigate another model that belongs to a more general class of oscillatormodels. II. MODEL
We consider globally coupled phase oscillators under global feedback, given as dθ i d ˜ t = ˜ ω i + ˜ KN N (cid:88) j =1 sin ( θ j − θ i + β ) + ˜ E sin( θ i + α ) F ( θ ) , (1)where θ i (˜ t ) and ˜ ω i ( i = 1 , . . . , N ) represent the phase and natural frequency of the i th oscillator,respectively; ˜ K ≥ E ≥ α and β are parameters, which are usually nonvanishing in real oscillator systems[33, 34]. Thefunction F ( θ ) = F ( θ , . . . , θ N ) describes a global feedback and sin( θ i + α ) is the phase sensitivityto the feedback. In particular, we consider F ( θ ) = 1 N N (cid:88) j =1 cos( θ j − δ ) , (2)= R cos(Θ − δ ) , (3)where δ is a parameter referred to as the phase offset in the feedback, and R = R ( t ) (0 ≤ R ≤ t ) (0 ≤ Θ < π ) are the order parameter and mean phase defined by r := Re iΘ = 1 N N (cid:88) j =1 e i θ j . (4)The R value indicates the synchronization level. The complex valued function r = r ( t ) is referredto as the complex order parameter. Equation (1) reduces to the Kuramoto–Sakaguchi model[35]in the absence of feedback, i.e., for ˜ E = 0. See Appendix A for the derivation of Eq. (1) froma general class of coupled limit-cycle oscillators. As discussed in Appendix A, F ( θ ) correspondsto a linear function of mean fields in the limit-cycle model introduced in Appendix. A, and theparameters ˜ E and δ can be tuned to arbitrary values when two output signals are observed fromindividual oscillators. Alternatively, it can be implemented when R ( t ) and Θ( t ) are inferred online.For analytical tractability, we assume ˜ ω i to be drawn from the Lorentzian distribution ˜ g (˜ ω ) = ˜ γπ ω − ω ) +˜ γ , where ˜ γ > ω >
0. Without loss of generality, we decrease the number ofparameters by introducing nondimensional quantities t = ω ˜ t, γ = ˜ γω , K = ˜ Kω , E = ˜ Eω , andfurther replacing θ i + α by θ i for i = 1 , . . . , N and δ + α by δ . The resultant equation is dθ i dt = ω i + KN N (cid:88) j =1 sin ( θ j − θ i + β ) + E sin θ i F ( θ ) , (5)or ˙ θ i = ω i + KR sin (Θ − θ i + β ) + ER cos(Θ − δ ) sin θ i , (6)where ω i is drawn from g ( ω ) = γπ ω − + γ . (7)Now, the mean frequency is set to unity. Equation (6) with Eqs. (2) and (7) is analyzed below.There are six parameters involved: N, K > , E > , β, δ , and γ .It is known that for N → ∞ , a certain class of oscillator assemblies including Eq. (6) has alow-dimensional manifold, on which a reduced dynamical equation can be obtained [31, 32]. Byfollowing [36], we obtain a closed equation for r on the manifold, given as˙ r = (cid:18) − γ + Ke i β (cid:19) r − Ke − i β | r | r − ER cos(Θ − δ )2 (cid:0) − r (cid:1) . (8)Here, we redefined r = Re iΘ as its continuous analogue r = (cid:90) (cid:90) ρ ( θ, ω, t ) e i θ dθdω, (9)where ρ ( θ, ω, t ) dθdω is the fraction of the oscillators with natural frequencies between ω and ω + dω and phases between θ and θ + dθ at time t . For finite but sufficiently large N , Eq. (8) is expectedto appropriately approximate the behavior of r ( t ) in Eq. (6) after a transient time.Henceforth, we assume | β | < π/
2; the coupling promotes synchronization.
III. EFFECT OF FEEDBACK ON THE MACROSCOPIC STATEA. Classification of macroscopic states
Dynamical properties of the system described by Eq. (8) in the absence of feedback ( E = 0)are evident. This system always has a global attractor for any parameter values within K ≥ − π < β < π . There are two global attractors: the steady solution r ( t ) = 0 for K ≤ K c andthe limit cycle solution r ( t ) = r e i ˜ ωt for K > K c , where K c = 2 γ/ cos β , r = [1 − K c /K ] / , and˜ ω = 1 + K sin β − γ tan β . The former corresponds to the asynchronous state, where the oscillatorsrotate with their natural frequencies. The latter corresponds to the synchronously oscillating state,where a subpopulation of the oscillators is phase-locked to the oscillating mean field. Figure 1 (a–c)shows typical dynamics of Eq. (6) for E = 0.In addition to these two attractors, a stable steady solution r ( t ) = r ∗ (cid:54) = 0 may arise in the pres-ence of the feedback ( E > r ( t ) = r ∗ ), let us consider Eq. (6). Inserting r ∗ = R ∗ e iΘ ∗ into Eq. (6), we obtain˙ θ i = ω i + A sin ( θ i + B ) , (10)where A = R ∗ [( E cos (Θ ∗ − δ ) − K cos (Θ ∗ + β )) + ( K sin (Θ ∗ + β )) ] / , (11)tan B = K sin (Θ ∗ + β ) E cos (Θ ∗ − δ ) − K cos (Θ ∗ + β ) . (12)A stable fixed point of Eq. (10) exists for | ω i | < | A | , which implies that the oscillators with | ω i | < | A | cease their oscillations. Thus, the solution r ( t ) = r ∗ (cid:54) = 0 corresponds to the oscillation-death state, where a subpopulation of the oscillators and the mean field quit oscillations. Figure1(d) and the blue dotted line in Fig. 1(a) exemplify the dynamics of the oscillators and the orderparameter, respectively, in the oscillation-death state. B. Bifurcation analysis and phase diagrams
The following bifurcation analysis of Eq. (8) enables us to obtain the phase diagrams shown inFig. 2.First, we analyze the bifurcation of the fixed point r = 0, i.e., the asynchronous state. Substi-tuting r = x + i y ( x, y ∈ R ) into Eq. (8) and linearizing the equations for dx/dt and dy/dt around( x, y ) = (0 , L = Λ − E cos δ − Ω − E sin δ Ω Λ , (13)where Λ = − γ + K cos β and Ω = 1 + K sin β . The quantity Ω represents the frequency of collectiveoscillation at its onset, i.e., at K = K c , in the absence of feedback. FIG. 1. (Color online) Three types of macroscopic state of the N = 100 oscillators described by Eq. (6). (a)Time series of the real part of the complex order parameter r for different parameter regimes. The black solid,red dashed, and blue dotted lines represent the solutions corresponding to the asynchronously oscillatingstate, the synchronously oscillating state, and the oscillation-death state, respectively. The coupling strengthand feedback parameters are set to K = 0 .
01 and E = δ = 0 (solid black line), K = 0 . E = δ = 0 (reddashed line), and K = 0 . E = 3 .
0, and δ = − π/ γ = 0 . β = 0. (b)–(d) Time series of the individual oscillators. The vertical axis is the index of the oscillator,and the color scale represents its phase. The natural frequency ω i is set to be monotonically increasing withthe oscillator index i : ω i = 1 + γ tan (cid:104) iπN − ( N +1) π N (cid:105) . The parameters for (b),(c), and (d) are identical tothose of the black solid, red dashed, and blue dotted lines in Fig. 1(a), respectively. The saddle-node, transcritical, pitchfork, and Hopf bifurcations are the codimension-one bifur-cations that alter the local stability and/or number of fixed points. The first three bifurcationsoccur when one of the eigenvalues of the stability matrix vanishes, which can be captured by itsnecessary condition | L | = 0. Hopf bifurcation occurs when both eigenvalues of L become purelyimaginary, i.e., Tr L = 0 and | L | > L = 0 for E under the condition | L | > E hopf = 4Λcos δ , (14)where δ satisfies − Λ + Ω + 2ΛΩ tan δ > . (15)The curve defined by Eqs. (14) and (15) is depicted by the orange dashed line in Fig.2(a)–(c). Because Tr L = 2Λ (1 − E/E hopf ), the asynchronous state is unstable for
E < E hopf (resp.
E > E hopf ) when Λ > < | L | = 0 yields another bifurcation curve: E pf = 2 (cid:0) Λ + Ω (cid:1) Λ cos δ − Ω sin δ . (16)On this curve, the pitchfork bifurcation occurs because Eq. (8) is invariant under the change of r → − r . See Appendix B for the transformation of Eq. (8) to a normal form and for a briefexplanation on why this bifurcation cannot be a transcritical or saddle-node bifurcation. Because | L | can be described as | L | = (cid:0) Λ + Ω (cid:1) (1 − E/E pf ), the asynchronous state is unstable for E > E pf when ΛΩ (cid:54) = 0. The stability of a fixed point changes through this bifurcation if the non-zeroeigenvalue of the stability matrix is negative; i.e., Tr L <
0, or2Λ < Λ + Ω Λ − Ω tan δ . (17)The curve defined by Eqs. (16) and (17) is higlighted using magenta solid line labeled as PFs inFig. 2. On this bifurcation curve, the asynchronous state loses its stability. Moreover, weaklynonlinear analysis performed in Appendix B implies that this bifurcation is subcritical in theparameter region considered in Fig. 2. This suggests the existence of the bistability between theasynchronous and oscillation-death states near the curve, which will be numerically confirmed inSec. III C.The pitchfork bifurcation occurs to an unstable fixed point if the non-zero eigenvalue of thestability matrix on the bifurcation curve is positive; i.e., Tr
L >
0. We thus obtain the curve forthis bifurcation, which is defined by Eq. (16) and Eq. (17) but with the opposite inequality. Theobtained curve is shown by the dot–dashed line labeled as PFu in Fig. 2. Because no stable stateemerges with this bifurcation, the curve is not relevant to the phase diagram.Next, we investigate the bifurcation of the fixed point r ∗ (cid:54) = 0, i.e., the oscillation-death state.For convenience, we express Eq. (8) using polar coordinates: dRdt = R (cid:26) − γ + (cid:0) − R (cid:1) (cid:20) K β − E δ − δ ) (cid:21)(cid:27) (18a) d Θ dt = 1 + (cid:0) R (cid:1) (cid:20) K β + E δ − sin ( δ − (cid:21) . (18b)A fixed point r ∗ = R ∗ e iΘ ∗ is given as a solution to dR/dt = 0 and d Θ /dt = 0; however, thisis difficult to solve explicitly. Nevertheless, we can still obtain the stability matrix by linearizingEq. (18) around ( R, Θ) = ( R ∗ , Θ ∗ ), as M = − R ∗ γ − R ∗ − E R ∗ (cid:0) − R ∗ (cid:1) sin( δ − ∗ ) − R ∗ R ∗ +1 E (cid:0) R ∗ + 1 (cid:1) cos( δ − ∗ ) . (19)Below, we show that the bifurcation curve can be obtained as a function of R ∗ and drawn in thephase diagram by varying R ∗ in the range of 0 < R ∗ < | M | = 0 and using the identity sin ( δ − ∗ ) + cos ( δ − ∗ ) = 1, we obtainsin ( δ − ∗ ) = ± (cid:112) ξ , (20a)cos ( δ − ∗ ) = ∓ ξ (cid:112) ξ , (20b)where ξ = (cid:0) R ∗ − (cid:1) γ (cid:0) R ∗ + 1 (cid:1) . (21)By substituting Eq. (20) into Eq. (18) and inserting dR/dt = d Θ /dt = 0, we obtain E = ± (cid:112) ξ + 1 (cid:0) (2 γ + KQ − cos β ) + γξ ( KQ + sin β + 2) (cid:1) ξ (cid:0) γ − KQ − cos β + γKQ +2 sin β (cid:1) , (22a)cos δ = ± ξ (2 γ + KQ − cos β ) (cid:0) γ − KQ − cos β + γKQ +2 sin β (cid:1)(cid:112) ξ + 1 Q − ((2 γ + KQ − cos β ) + γξ ( KQ + sin β + 2) ) ± ξ (cid:112) ξ + 1 , (22b)sin δ = ∓ ξ ( KQ + sin β + 2) (cid:0) γ − KQ − cos β + γKQ +2 sin β (cid:1)(cid:112) ξ + 1 Q + ((2 γ + KQ − cos β ) + γξ ( KQ + sin β + 2) ) ± (cid:112) ξ + 1 , (22c)where Q + = 1 + R ∗ and Q − = − R ∗ . Inserting 0 < R < r ∗ (cid:54) = 0 below the curve. Such a fixed point may not disappear unless abifurcation involving that point occurs, thus it should persist up to the parameter region of E = 0.However, this contradicts the fact that no isolated fixed point except r = 0 exists at E = 0. Thus,a saddle–node bifurcation occurs on the curve. We now know that a pair of fixed points arisethrough the bifurcation, but their stability remains unclear. Expressing Tr M as a function of R ∗ ,we numerically confirmed that Tr M < δ . To this end, for some parameter regions, we derive a necessary condition and a sufficientcondition for a nonzero fixed point ( R ∗ , Θ ∗ ) to exist. If a nonzero fixed point ( R ∗ , Θ ∗ ) exists, itmust satisfy d Θ /dt | Θ=Θ ∗ = 0, or, E = 4 + 2 K (1 + R ∗ ) sin β (1 + R ∗ )(sin( δ − ∗ ) − sin δ ) . (23)The numerator on the right-hand side is positive when Ω (cid:48) := 1 + K sin β >
0. In this case, a lowerbound E lower ≤ E is given by Eq. (23) with R ∗ = 1 and sin( δ − ∗ ) = 1, i.e., E lower = 2Ω (cid:48) − sin δ . (24)Thus, when Ω (cid:48) >
0, a necessary condition for the existence of a nonzero fixed point is given by E ≥ E lower . (25)Next, we address a sufficient condition for the existence of a nonzero fixed point. We restrictourselves to the parameter region where Ω (cid:48) > K > K c holds. Then, for the parameterregion under consideration, as shown in Appendix D, a nonzero fixed point exists if E satisfiesboth of the following inequalities: E > E lower (1 + η ) (26)and E < K cos β δ , (27)where η = K c /K (2 − K c /K ) (1 + K sin β ) . (28)We can find such E for β (cid:39) K when a value of δ is given. See Appendix Dfor details.0From Eq. (25) and (26), for Ω (cid:48) > β (cid:39)
0, and sufficiently large K , it follows that the saddle–node bifurcation should occur at E ∈ [ E lower , E lower (1 + η )]. As η → K → ∞ , the bifurcationcurve is well approximated by E lower for large K . In Fig. 2(c), we observe that E lower , which isdepicted as the approximate SN (aSN) curve, approximates the saddle–node bifurcation curve well,although the sufficient condition given by Eqs. (26) and (27) cannot be satisfied in a range of δ .Similarly to the saddle–node bifurcation curve, the Hopf bifurcation curve is obtained by im-posing Tr M = 0 and | M | >
0. The former condition yields E = (cid:118)(cid:117)(cid:117)(cid:116) (cid:2) Q − s + ( c + 2 γ (2 R ∗ + 1)) − γ R ∗ (cid:3) Q Q − s + 16 γ R ∗ Q Q − , (29a)cos δ = 2 (cid:0) c + 2 γ (cid:0) R ∗ + 1 (cid:1)(cid:1) EQ + Q − , (29b)sin δ = − sEQ + − (cid:115) − γ R ∗ E Q Q − , (29c)where (cid:16) c + 2 γ (cid:16) R ∗ + 1 (cid:17)(cid:17) − γ R ∗ + s Q − ≤ , (30) s = K (cid:16) R ∗ (cid:17) sin β + 2 , (31) c = K ( − R ∗ ) cos β. (32)Inserting 0 < R ∗ < | M | > >
0, a heteroclinic bifurcation curve extends from this point. The navysolid curve labeled by HC shows the heteroclinic bifurcation curve obtained using the softwareXPPAUT[38].Figure 2 indicates that the heteroclinic bifurcation curve ends up at a point on the saddle–nodebifurcation curve, where another codimension-two bifurcation should occur. The observation ofvector fields around this bifurcation point, shown in Fig. 3, reveals that the saddle–node bifurcationcurve is subdivided into two parts at this point, and the one labeled by SNIC in Fig. 2 correspondsto saddle–node bifurcation on an invariant circle. A very similar codimension-two bifurcation isreported in [39], where SN, SNIC, and homoclinic bifurcation curves meet. We obtain a heteroclinic1bifurcation curve rather than a homoclinic one because of the invariance of our system under r → − r .The three bifurcation curves HB, HC, and SNIC in Fig. 2 form the boundary of the stablesynchronous state, provided that there is no other bifurcation involving periodic solutions, such asthe saddle–node bifurcation of limit cycles. Then, the bistable region exists in the area surroundedby PFs, HC, and SN. Altogether, we obtain three qualitatively different phase diagrams dependingon K values, as shown in Fig. 2. Our extensive numerical analysis performed in Sec. III C andAppendix F indicates that the phase diagrams given in Fig. 2 are comprehensive.Other bifurcations may occur typically for Ω < <
0. The situation Ω < − π < β < K . Insuch a situation, in the absence of feedback, the mean field oscillates with a frequency oppositeto that of the typical natural frequency of individual oscillators, i.e., ω = 1. This implies thatsynchronized individual oscillators also have frequencies opposite to their natural ones. However,such a phenomenon is not commonly observed in limit-cycle oscillators and should be regarded asan artifact owing to the use of phase approximation for a case of strong coupling. We thereforeomit the case of Ω < C. Numerical verification
To verify the the analyses in Sec. III B, we simulated Eq. (6) for N = 2000. We set ω i = ω + γ tan (cid:104) iπN − ( N +1) π N (cid:105) with ω = 1, which converges to the Lorentzian distribution at N → ∞ [40].In Fig. 4(a) and its magnifications Fig. 4(b) and (c), we show the value of (cid:104) R (cid:105) for different feedbackparameters, where the angle brackets denote the long-time average. As the initial condition, weemploy the uniform state, i.e., θ i (0) = 2 π ( i − /N for i = 1 , · · · , N , in Fig. 4(a,b) and the fullysynchronized state, i.e., θ i (0) = 0 for i = 1 , · · · , N , in Fig. 4(c). The parameters are same as inFig. 2(b), and we draw the same bifurcation curves in Fig. 4(a–c) for comparison.In the black regions in Fig. 4, (cid:104) R (cid:105) (cid:39) FIG. 2. (Color online) Phase diagrams of the macroscopic state based on the stable solutions of Eq.(8) for(a) K = 0 . < K c , (b) K = 0 . > K c , and (c) K = 1 . (cid:29) K c , where K c = 2 γ/ cos β = 0 .
2. Other parametersare γ = 0 . β = 0. The black, green, and shaded regions correspond to the asynchronous, synchronouslyoscillating, and oscillation-death states, respectively. death states is predicted. To clarify which region of nonvanishing (cid:104) R (cid:105) in Fig. 2 corresponds to thesynchronously oscillating or oscillation-death states, we further measure (cid:104) ζ (cid:105) , where ζ = | r − (cid:104) r (cid:105)| .From the definition of R and ζ , (cid:104) R (cid:105) > (cid:104) ζ (cid:105) = 0 imply that the system is in the oscillation-deathstate, thus we confirm the predicted bistability as well as the existence of the stable oscillation-deathstate inside the SN and SNIC curves. IV. OPTIMAL FEEDBACK PARAMETERS
We consider
K > K c , or equivalently Λ := − γ + K cos β >
0, for which the system falls intothe synchronously oscillating state in the absence of the feedback, and determine the value ofthe phase offset δ that minimizes the required feedback strength E to suppress the synchronized3 FIG. 3. (Color online) Typical vector fields of Eq. (8) for different areas of the phase diagrams in Fig. 2.The vector fields on the complex planes are drawn with cyan arrows. Filled (resp. open) squares representstable (resp. unstable) spirals, and filled (resp. open) circles represent stable (resp. unstable) nodes. Theopen triangles represent saddles. Stable limit cycles are illustrated using black solid curves. The circle drawnwith the dashed line on each panel depicts the unit circle on the complex plane. The parameters for (a) and(b) are the same as those shown in Fig. 2(a) and (b), respectively. FIG. 4. (Color online) Simulation results of Eq. (6). Long-time average of (a–c) R and (d–f) ζ . Parametersare the same as those in Fig. 2(b), and the same bifurcation curves are drawn here. The parameter range in(b,c,e,f) is the same as that in the boxed area in (a,d). The initial condition is the uniform phase distributionin (a,b,d,e) and the fully synchronized state in (c,f). oscillations. We can achieve this by leading the system to (i) the asynchronous state and (ii) theoscillation-death state. Their optimal parameter sets ( δ, E ) are denoted by (i) ( δ async , E async ) and(ii) ( δ death , E death ).The point ( δ async , E async ) can be determined analytically. Because the asynchronous state isstable for E above the HB curve and below the PFs curve, ( δ async , E async ) is given by the minimum5of the HB curve. Further, δ = 0 provides the minimum when Eq. (15) holds for δ = 0, resulting in δ async = 0 , (33a) E async = 4Λ = − γ + 2 K cos β. (33b)This is the case when Λ ≤ | Ω | = | K sin β | , which typically arises for small K or large tan β .For Λ > | Ω | , the smallest | δ | that satisfies Eq. (15) provides the minimum; i.e., δ async = Ω | Ω | arcsin (cid:18) Λ − Ω Λ + Ω (cid:19) , (34a) E async = 2 (cid:0) Λ + Ω (cid:1) | Ω | . (34b)Although ( δ death , E death ) can only be numerically determined using Eqs. (22) and (29), anapproximate expression can be obtained from (24): δ death ≈ − π , (35a) E death ≈ K sin β. (35b)Figure 5 shows the parameter-dependency of E async given by Eqs. (33b) and (34b) and E death obtained numerically using Eqs. (22) and (29). The general tendency is well captured by Eqs. (33b)and (35b). The solid lines represent the parameter set at which E async = E death . Based onEqs. (33b) and (35b), we can roughly estimate that the asynchronous (oscillation-death) state canbe achieved with a smaller feedback strength when − γ + 2 K cos β is small (large) compared to1 + K sin β .When we desire to suppress the collective oscillation without causing oscillation death, we needto consider the bistability between the asynchronous and oscillation-death states. For example,see Fig. 2(c), where δ async = 0. Suppose that we increase E while fixing δ = 0. Then, theoscillation-death state will be obtained before the asynchronous state. By further increasing E ,we will eventually arrive at the HB curve, above which the asynchronous state is stable. However,because of the bistability, the oscillation-death state is expected to be sustained even in that region.Therefore, to realize the asynchronous state, we need to use a larger δ value at which we first arriveat the monostable region of the asynchronous state. Once the asynchronous state is realized, wecan vary δ to δ = 0 and decrease E to E async . To keep E as small as possible during the wholemanipulation, we should employ a δ value close to that of the intersection of the HB and SN curves.Using the HB curve given by (14) and aSN curve given by Eq. (24), the approximate intersection6 FIG. 5. (Color online) Parameter dependency of (a,b) E async and (c,d) E death . The solid lines representsthe parameter values at which E death = E async holds. (a,c) γ = 0 .
04. (b,d) γ = 0 .
4. On the left side of thedashed lines in (a) and (b), which depict K cos β = 2 γ , we have E async = 0 because the solution r = 0 isstable even without feedback. can be found as( δ, E ) ≈ (cid:32) π − α, Ω (cid:48) + 4Λ Ω (cid:48) (cid:33) , (36)where α = arcsin (cid:34) Ω (cid:48) (cid:112) Ω (cid:48) + 4Λ (cid:35) . (37)Using this δ value, we can efficiently steer the system into the asynchronous state.Note that in contrast to the case of the stabilization of the asynchronous state, the feedbackdoes not vanish when the oscillation-death state is achieved. Moreover, the minimum value of E does not necessarily imply that the intensity of the feedback | Ef ( r ) | is minimized as it also dependson r . Instead, this optimization does minimize the possible feedback intensityWe perform numerical simulations of Eq. (6) to verify whether a near-optimal feedback properlyworks. Figure 6 shows the time-series of the collective oscillation Re ( r ( t )) and the individual phases7 FIG. 6. (Color online) Alteration of the dynamics of the N = 100 oscillators by the feedback. In each figure,the upper panel displays the time series of Re ( r ) while the lower panel shows the phases of the oscillators.(a) The feedback parameters are set as E = 0 . ≈ E async and δ = 0 = δ async to stabilize the asynchronousstate with small E . The intrinsic parameters of the oscillators, namely K , β , and γ , are the same as thosein Fig. 2(b). (b) The feedback parameters are set as E = 1 . ≈ E death and δ = − π/ ≈ δ death to induce theoscillation-death state with small E in the population that has the same parameters as Fig. 2(c). θ i ( t ) before and after the onset of the feedback. In Fig. 6(a), the feedback with the parameters E = 0 . ≈ E async and δ = 0 = δ async is applied at t = 5050, as marked by the black arrow. Uponthe onset of the feedback, the population begins to be desynchronized, and r decreases with time.In Fig.6(b), the feedback parameters are chosen such that the oscillation-death state is inducedwith small feedback strength: E = 1 . ≈ E death and δ = − π/ ≈ δ death . The figure indicates thatthe oscillations of the individual oscillators and the mean field terminate immediately because ofthe feedback. V. INVESTIGATION ON THE ROBUSTNESS OF THE EFFECT OF FEEDBACKUSING A DIFFERENT MODEL OF OSCILLATORS
Our analyses thus far are based on the model given by Eq. (6). As presented in Appendix A,this model is derived from a general class of coupled oscillator models. However, we assumed thatcoupling, inhomogeneity, and feedback are sufficiently weak to employ averaging approximations8and that the functions contain only the first harmonics. Furthermore, we assumed the naturalfrequencies to obey Lorentzian distribution to obtain the reduced dynamical equation given inEq. (8). Here, to exemplify the robustness of the results to the violation of these assumptions, weprovide numerical results for a model with the form given by Eq. (A5). Specifically, we consider˙ θ i = 1 + µ i cos θ i + KN Z v ( θ i ) N (cid:88) j =1 V ( θ j ) + EZ f ( θ i ) f ( θ ) . (38)We adopt the pulse-like signal V ( θ ) = ν n (1 + cos θ ) n used in [41], where n is the parameter onthe sharpness of V , and ν n = 2 n ( n !) / (2 n )! normalizes (cid:82) π V ( θ ) dθ to be 2 π . We set n = 10. Thephase sensitivity functions Z v and Z f are chosen to weakly include the second Fourier mode: Z v ( θ ) = − sin θ + 0 . θ, (39) Z f ( θ ) = sin θ + 0 . θ. (40)Finally, µ i is drawn from Gaussian distribution with mean 0 and standard deviation 0 . (cid:104) R (cid:105) and (cid:104) ζ (cid:105) , which are shown inFig. 7(a) and (b), respectively. These figures qualitatively agree with Fig. 4(a) and (d), suggestingthe robustness of the results. VI. CONCLUSION AND DISCUSSION
Motivated from the wide range of the applications of synchronization control, we analyzed an in-homogeneous population of phase oscillators exposed to global feedback. Detailed phase diagramsof the collective state are obtained based on the bifurcation analysis of the macroscopic equationderived using the Ott–Antonsen theory. The diagram displayed three types of macroscopic states,namely, the synchronously oscillating state, the asynchronous state, and the oscillation-death state.Exact and approximate optimizations of the feedback parameters for steering a synchronously os-cillating population into the asynchronous or the oscillation-death state with minimum feedbackstrength are also presented. Although we assumed several conditions such as the weakness ofthe coupling and the feedback in the derivation of the model equation, the numerical investiga-tion in Sec. V demonstrates that our results do not change qualitatively even when some of theassumptions are violated to some extent.Herein, we focused on linear feedback F given by Eq. (2), and our extensive analysis revealed itsutility for synchronization control. Linear feedback can be regarded as a basic methodology, andour study is expected to serve as a benchmark when more sophisticated feedback is to be designed.9 FIG. 7. Long-time average of (a) R and (b) ζ in the system of pulse-coupled oscillators under the feedbackdescribed by Eq. (38). Initial conditions are given by θ i (0) = 2 π ( i − /N . A natural extension is to make the feedback function F nonlinear in R . Although it will not changethe linear stability of the asynchronous state, it may alter the stability and the existence of otherstates in addition to the amplitude of the collective oscillation[42] A class of nonlinear feedbackhas been proposed in [23] for reducing the amplitude of the collective oscillation, and furtherinvestigated in successive studies [42–46]. In [47, 48], the authors showed that a class of delayednonlinear feedback can stabilize complex dynamical states including a type of desynchronized stateand demonstrated its ability to control electrochemical oscillators. Furthermore, regarding DBS,smooth feedback may not fulfill safety requirements[49]. In [28], the authors compare two types0of feedback: a smooth feedback and a series of pulses that are amplified according to the smoothfeedback. They found that the pulses have a similar desynchronizing effect to the smooth feedback.The same approach might be applicable to the feedback studied in this article.Our theoretical results can be demonstrated in some experimental systems. Recently, techniquesfor inferring phase dynamics and their interactions from observed oscillatory signals have beendeveloped, and they have been utilized in experimental studies for predicting and controlling thedynamics of the oscillators [1, 34, 50–52]. These experiments have revealed various coupling andphase sensitivity functions; in some systems, the first Fourier component is dominant, while inother systems, some of the higher and the constant components are also prominent. In the formercase, our analytical results may be verified quantitatively by estimating the values of ω i , K, and β in Eq. (5) and implementing the global feedback loop. The feedback parameters δ and E can betuned if R and Θ can be inferred online or two outputs from individual oscillators are available, asdetailed in Appendix A.Finally, we remark on the limitation of the current study. it should be noted that our resultsare based on phase oscillator models. A qualitatively different phase diagram may be obtained forlimit-cycle oscillators whose amplitudes considerably deviate from that of the unperturbed periodicorbit. Therefore, to understand the effect of a large amplitude deviation on synchronization control,it is important to investigate models of limit-cycle oscillators and compare them with our results. ACKNOWLEDGMENTS
The authors thank Kei-Ichi Ueda for helpful discussions on numerical bifurcation analysis. H.K.acknowledges the financial support from JSPS KAKENHI Grant No. 18K11464.H.K. conceived the study and supervised the project. A.O. performed the analytical and nu-merical investigations. A.O. wrote the manuscript with support from H.K..
Appendix A: Derivation of our model given in Eq. (1)
We consider a general class of oscillator network given by˙ x i = u i ( x i ; p i , q i ) , (A1)where x i = ( x i , y i , . . . ) and u i ( i = 1 , . . . , N ) are the state and the local vector field of the i thoscillator, respectively, and p i and q i are parameters. The interactions and external forcing areassumed to be given through variations in the parameters as p i = p + ∆ p i and q i = q + ∆ q i ,1where p and q are the parameter values common for all the oscillators and ∆ p i and ∆ q i describeperturbations. When we consider global coupling and feedback, the perturbations may be given as∆ p i = K (cid:48) N N (cid:88) j =1 v ( x i , x j ) , (A2)∆ q i = ∆ q ≡ E (cid:48) f ( x , . . . , x N ) , (A3)where K (cid:48) and E (cid:48) are the strengths of coupling and feedback, respectively, and v and f are thefunctions describing coupling and feedback, respectively.We follow the standard procedure of the phase reduction to obtain the corresponding phasedescription to Eq. (A1) [5, 53]. Let us introduce ∆ u i ( x i ; p i , q i ) = u i ( x i ; p i , q i ) − u ( x i ; p i , q i ) for i = 1 , . . . , N , which describes inhomogeneity in inherent oscillator properties. We assume that theunperturbed system, i.e.,˙ x = u ( x ; p , q ) , (A4)has a stable limit cycle x ∗ ( t ). The phase of the unperturbed system is defined as a scalar field Φ( x )for the basin of attraction for the limit cycle such that the contour of Φ( x ) describe the isochronof the unperturbed system[5, 54]. Using this scalar field, the phase of the i th oscillator is definedas θ i = Φ( x i ( t )) ( i = 1 , . . . , N ).We assume that the orbital stability of the cycle x ∗ ( t ) in the unperturbed system given by Eq.(A4) is sufficiently higher than perturbation strength. Then, to the lowest order in perturbationsstrengths, each phase obeys˙ θ i = ω + Z ( θ i ) · U i ( θ i ) + K (cid:48) N N (cid:88) j =1 Z v ( θ i ) V ( θ i , θ j ) + E (cid:48) Z f ( θ i ) F ( θ ) , (A5)where ω is the natural frequency of the limit cycle x ∗ ; U i , V, F are the parametric representationsof ∆ u i , v, f on the limit-cycle x ∗ in terms of the phases, respectively; and Z , Z v , Z f are the phasesensitivity functions, which can be expressed in terms of the derivatives of u ( x ; p, q ) and Φ ( x ).All the functions are 2 π -periodic in each argument.When the magnitudes of the perturbation terms, i.e., the second to fourth terms of the right-hand side in Eq. (A5), are sufficiently small compared to ω , we may further simplify the equationusing an averaging approximation. The resultant equation is˙ θ i = ω + ∆ ω i + K (cid:48) N N (cid:88) j =1 Γ v ( θ i − θ j ) + E (cid:48) N N (cid:88) j =1 Γ f ( θ i − θ j ) , (A6)2where the constant ∆ ω i and the functions Γ v and Γ f can be expressed in terms of the functionsappearing in Eq. (A5). If the last term in Eq. (A5) is not very small compared to ω , we can stillaverage the other terms to obtain˙ θ i = ω + ∆ ω i + K (cid:48) N N (cid:88) j =1 Γ v ( θ i − θ j ) + E (cid:48) Z f ( θ i ) F ( θ ) . (A7)In Eq. (A7), we may consider larger E (cid:48) values than in Eq. (A6), which is the reason why we employthis type of model in this work.Our model given in Eq. (1) is a version of Eq. (A7), where we assume that only the first harmon-ics are present in all the functions appearing in Eq. (A7); i.e., Γ v , Z f and F . This assumption isvalid when we consider limit-cycle oscillators close to a Hopf bifurcation point, in which the phasesensitivity and the wave forms are well approximated by the functions with only the first harmonicsand with the constant and first harmonics terms, respectively. Let us further assume that fromeach oscillator we can observe two quantities, such as x j ( t ) and y j ( t ). If we denote the trajectoryof the unperturbed limit cycle x ∗ ( t ) by χ (Φ), i.e., χ (Φ( x ∗ ( t ))) = x ∗ ( t ), with its elements being χ = ( χ x , χ y , . . . ), the variation of χ x ( θ ) and χ y ( θ ) are almost sinusoidal near the Hopf bifurcationpoint. Therefore, the unperturbed waveforms can be denoted by χ x ( θ ) (cid:39) ¯ χ x + A x cos( θ − δ x ) and χ y ( θ ) (cid:39) ¯ χ y + A y cos( θ − δ y ), where ¯ χ x and ¯ χ y are the average of χ x ( θ ) and χ y ( θ ), respectively, A x and A y are the oscillation amplitudes, and δ x and δ y are the phase offsets in the waveforms. Wethen give the feedback function f as f = N (cid:88) j =1 (cid:20) a x j − ¯ x j A x + b y j − ¯ y j A y (cid:21) , (A8)where ¯ x j and ¯ y j denote the average values of x j and y j , respectively; and a and b are our controlparameters. In the lowest order phase description, Eq. (A8) results in F ( θ ) = N (cid:88) j =1 [ a cos( θ j − δ x ) + b cos( θ j − δ y )] . (A9)We can further transform (A9) to F ( θ ) = E N (cid:88) j =1 cos ( θ j − δ ) (A10)= E R cos (Θ − δ ) , (A11)where E = (cid:104) ( a cos δ x + b cos δ y ) + ( a sin δ x + b sin δ y ) (cid:105) / , (A12)tan δ = a sin δ x + b sin δ y a cos δ x + b cos δ y . (A13)3Therefore, we can give arbitrary E and δ values by appropriately assigning a and b values. Because˜ E in Eq. (1) is given by ˜ E = E (cid:48) E , we can also give an arbitrary ˜ E value. Appendix B: Classification of the zero-eigenvalue bifurcation at the origin
The codimention-one bifurcation involving zero-eigenvalue at the origin is limited to the pitch-fork bifurcation. One possible approach for this is to consider the facts that the saddle–nodebifurcation may not occur because the constant solution r = 0 may not vanish in Eq. (8) andthe pitchfork bifurcation rather than the transcritical bifurcation occurs because the symmetry ofEq. (8) implies that the emergence of a constant solution r = r ∗ (cid:54) = 0 must be accompanied withthe emergence of r = − r ∗ as well.An alternative way is to perform the center manifold reduction[37], which will clarify that thebifurcation is actually the pitchfork one and whether the bifurcation is super- or sub-critical. Letthe bifurcation parameter be µ = E − E pf , where E pf is given by Eq. (16). Inserting r = u + i v and E = E pf + µ into Eq. (8), we obtain ˙ u ˙ v = Λ − E pf cos δ − Ω − E pf sin δ Ω Λ uv + p ( u, v ) q ( u, v ) , (B1)where p ( u, v ) = − µ u cos δ + v sin δ ) + E pf + µ (cid:0) u − v (cid:1) ( u cos δ + v sin δ ) − K (cid:0) u + v (cid:1) ( u cos β + v sin β ) , (B2) q ( u, v ) = u v [( E pf + µ ) cos δ − γ − Λ] + uv [( E pf + µ ) sin δ + Ω − u (Ω −
1) + v ( − γ − Λ) . (B3)To reduce the system, let us transform the variables as follows: uv = − ΛΩ Λ sin δ +Ω cos δ Ω sin δ − Λ cos δ ˆ u ˆ v , (B4)which yields ˙ˆ u ˙ˆ v = λ ˆ u ˆ v + a µ ˆ u + a µ ˆ v + a ˆ u + O (cid:0) v , uv , u v (cid:1) bµ ˆ u + O (cid:0) µ ˆ v, ˆ v , ˆ u ˆ v , ˆ u ˆ v, ˆ u (cid:1) (B5)4where ˆ λ = cos δ (cid:0) Λ − Ω (cid:1) − δ Λ cos δ − Ω sin δ , (B6) a = (Ω sin δ − Λ cos δ ) − Λ ) cos δ + 2ΛΩ sin δ ] , (B7) a = Ω − Λ ) cos δ + 2ΛΩ sin δ ] , (B8) a = (cid:0) Λ + Ω (cid:1) (cid:8)(cid:0) − γ ΛΩ + Λ − Ω (cid:1) sin δ + (cid:2) γ (cid:0) Λ − Ω (cid:1) + 2ΛΩ (cid:3) cos δ (cid:9) Ω [(Ω − Λ ) cos δ + 2ΛΩ sin δ ] , (B9) b = (Λ cos δ − Ω sin δ ) − Ω ) cos δ − δ ) . (B10)The center manifold ˆ v = c (ˆ u, µ ) up to the second order is given by c (ˆ u, µ ) = − b ˆ λ µ ˆ u + O (cid:0) ˆ uµ , ˆ u (cid:1) . (B11)On this center manifold, the dynamics is reduced to˙ˆ u = a µ (cid:18) − a b ˆ λa µ (cid:19) ˆ u + a ˆ u + O (cid:0) µ ˆ u , µ u (cid:1) . (B12)Equation (B12) implies that the pitchfork bifurcation occurs at µ = 0. The origin changes itsstability through this bifurcation when ˆ λ <
0. It is supercritical for a < a > λ <
0, Ω >
0, and E pf is sufficientlysmall compared to Ω /γ . Because E ≥
0, and the numerator of E pf is positive, the denominator of E pf is also positive:Λ cos δ − Ω sin δ > , (B13)which implies that the denominator of ˆ λ is positive. Thus, if ˆ λ <
0, the numerator of ˆ λ is negative,i.e., (cid:0) Λ − Ω (cid:1) cos δ − δ < . (B14)Under this condition, a is positive if (cid:0) γ (cid:0) Λ − Ω (cid:1) + 2ΛΩ (cid:1) cos δ > (cid:0) γ ΛΩ − Λ + Ω (cid:1) sin δ (B15) ⇐⇒ Ω (Λ cos δ − Ω sin δ ) > − Λ (Λ sin δ + Ω cos δ ) + γ (cid:2) δ − (cid:0) Λ − Ω (cid:1) cos δ (cid:3) (B16)We first consider the case of Λ ≥
0. Then, the first term of the right-hand side of Eq. (B16) isnegative because Eqs. (B13) and (B14) yield0 ≤ Λ (Λ cos δ − Ω sin δ ) < Ω (Λ sin δ + Ω cos δ ) . (B17)5As for the second term, noting that | cos δ | , | sin δ | ≤
1, we obtain2ΛΩ sin δ − (cid:0) Λ − Ω (cid:1) cos δ ≤ + Ω ≤ (cid:0) Λ + Ω (cid:1) (B18)Equations (B16), (B17), and (B18) yield the following sufficient condition for a to be positive:Ω (Λ cos δ − Ω sin δ ) > γ (cid:0) Λ + Ω (cid:1) , (B19)which is equivalent to E pf = 2 (cid:0) Λ + Ω (cid:1) Λ cos δ − Ω sin δ < Ω γ . (B20)Therefore, for Λ ≥
0, the bifurcation is subcritical when E pf is sufficiently small compared to Ω /γ .We next consider the case of Λ <
0. If Λ sin δ + Ω cos δ < ≥
0. When Λ sin δ + Ω cos δ ≥
0, we evaluate the right-hand sideof Eq. (B16) as follows. As | Λ | < γ holds for Λ < − Λ (Λ sin δ + Ω cos δ ) < γ ( | Λ | + Ω) . (B21)This equation and Eqs. (B16) and (B18) yield the sufficient condition for the subcritical bifurcationΩ (Λ cos δ − Ω sin δ ) > γ (cid:0) Λ + Ω (cid:1) (cid:18) | Λ | + Ω2 (Λ + Ω ) (cid:19) , (B22)which is equivalent to E pf < Ω γ (cid:16) | Λ | +Ω2(Λ +Ω ) (cid:17) . (B23)From Eq.(B20) and Eq.(B23), it is shown that the pitchfork bifurcation is subcritical if E pf issufficiently small compared to Ω /γ .Near the parameter regions considered in Fig. 2, E pf is small enough, and hence the pitchforkbifurcation involving a stable fixed point is subcritical. Appendix C: Weakly nonlinear analysis of the Hopf bifurcation
We show below that the Hopf bifurcation at the origin is always supercritical. Substituting r = u + i v and the value of the feedback strength at the bifurcation point, given by Eq. (14), intoEq. (8), we have ˙ u ˙ v = − Λ −
2Λ tan δ − ΩΩ Λ uv + g ( u, v ) h ( u, v ) , (C1)6where g ( x, y ) = 2Λ (cid:0) x − y (cid:1) ( x + y tan δ ) − K (cid:0) x + y (cid:1) ( x cos β + y sin β ) , (C2) h ( x, y ) = K (cid:0) x + y (cid:1) ( x sin β − y cos β ) + 4Λ xy ( x + y tan δ ) (C3)To simplify the calculation, we change the coordinates to ˜ u ˜ v = Ω / ˜ ω Λ / ˜ ω uv , (C4)where ˜ ω = √− Λ + Ω + 2ΛΩ tan δ . Then, ˜ u and ˜ v obey ˙˜ u ˙˜ v = − ˜ ω ˜ ω ˜ u ˜ v + ˜ g (˜ u, ˜ v )˜ h (˜ u, ˜ v ) , (C5)where ˜ g (˜ u, ˜ v )˜ h (˜ u, ˜ v ) = Ω / ˜ ω Λ / ˜ ω g ( u (˜ u, ˜ v ) , v (˜ u, ˜ v )) h ( u (˜ u, ˜ v ) , v (˜ u, ˜ v )) . (C6)Then, z = u + i v obeys˙ z = i˜ ωz + ˜ g + i˜ h. (C7)Using a near-identity transformation, Equation (C7) is further reduced to the following normalform of the Hopf bifurcation [37]˙ w = i˜ ωw + d | w | ¯ w + O (cid:16) | w | (cid:17) , (C8)where w, d ∈ C , andRe ( d ) = − γ (cid:18) δ (cid:19) . (C9)The bifurcation is supercritical when Re ( d ) <
0. Noting that γ >
0, we may further showRe ( d ) < ΛΩ tan δ ≥
0, it is obvious that Re ( d ) < ΛΩ tan δ <
0, the following inequality holds:1 + ΛΩ tan δ > δ. (C10)Moreover, Eq. (15) implies that the Hopf bifurcation at the origin occurs only when − Λ + Ω + 2ΛΩ tan δ > . (C11)From Eqs. (C10) and (C11), we have1 + ΛΩ tan δ > (cid:18) ΛΩ (cid:19) > d ) <
0. Therefore, the Hopf bifurcation at the origin is supercritical for any param-eter values.7
Appendix D: Derivation of the sufficient condition for a nonzero fixed point
Here, under the conditions
K > K c and Ω (cid:48) = 1 + K sin β >
0, we derive Eqs. (26) and (27),i.e., a sufficient condition for Eq. (18) to have a nonzero fixed point. Note that an intersection ofthe two nullclines − γ + (cid:0) − R (cid:1) (cid:20) K β − E δ − δ ) (cid:21) = 0 (D1)and 1 + (cid:0) R (cid:1) (cid:20) K β + E δ − sin ( δ − (cid:21) = 0 (D2)gives a nonzero fixed point. Thus there exists a fixed point ( R ∗ , Θ ∗ ) such that R lower < R ∗ < R on the nullcline always satisfies R lower < R ∗ <
1, and (ii) the nullcline givenby Eq. (D2) passes through the region R lower < R ∗ < R -Θ plane.The first condition is equivalent to that the following inequality holds for any Θ: R lower2 < − γ K cos β − E [cos ( δ − δ ] < , (D3)which is satisfied when E < K cos β − γ (cid:0) − R lower2 (cid:1) − δ . (D4)Next we discuss the second condition. Equation (D2) yields R = − S (sin ( δ − − ,where S ( x ) = − K β + E x − sin δ ) (D5)is a monotonically increasing function of x . Because sin ( δ − , R lower < R ∗ < S ( − <
12 (D6)and S (1) > (cid:0) R (cid:1) − . (D7)When Ω (cid:48) >
0, the inequality given by (D6) holds for any E ≥
0. In contrast, inequality (D7) holdswhen
E > E lower (cid:34) − R lower2 (cid:0) R lower2 (cid:1) (1 + K sin β ) (cid:35) , (D8)8where E lower is given by Eq. (24).Therefore, if Eqs. (D4) and (D8) are satisfied, a nonzero fixed point ( R ∗ , Θ ∗ ) that satisfies R lower < R ∗ < R lower = (cid:113) − K c K in Eqs. (D8) and (D4).We can find a value of E that satisfies both of Eqs. (26) and (27) when E lower (1 + η ) < K cos β δ , (D9)which is equivalent to21 − sin δ (cid:20) K − K c + sin β (cid:21) − cos β δ < . (D10)Equation (D10) holds when β (cid:39) K is sufficiently large for a given value of δ . Appendix E: Analysis on the codimension-two bifurcation point
By imposing Tr L = | L | = 0, we obtain the values of the feedback parameters at which the Hopfbifurcation curve (14) and the pitchfork bifurcation curve (16) meet assin δ = Ω | Ω | Λ − Ω Λ + Ω , (E1)cos δ = 2Λ | Ω | Λ + Ω , (E2) E = 2 (cid:0) Λ + Ω (cid:1) | Ω | . (E3)Substituting r = u + i v together with Eq. (E1)–(E3) into Eq. (8), we have ˙ u ˙ v = − Λ − Λ Ω Ω Λ uv + g ( u, v ) h ( u, v ) , (E4)where g ( u, v ) = (Λ − γ ) u + (cid:18) −
2Ω + Λ Ω (cid:19) u v − ( γ + 3Λ) uv + (cid:18) − Λ Ω (cid:19) v , (E5) h ( u, v ) = (Ω − u − ( γ − u v + (cid:18) − − Ω + 2Λ Ω (cid:19) uv − ( γ + Λ) v . (E6)Next, we transform the linear part into the Jordan normal form by the following change of variables: ˇ u ˇ v = uv , (E7)9which yields ˙ˇ u ˙ˇ v = ˇ v + c ˇ u + c ˇ u ˇ v + c ˇ u ˇ v + c ˇ v ,d ˇ u + d ˇ u ˇ v + d ˇ u ˇ v + d ˇ v , (E8)where c = (cid:0) Λ + Ω (cid:1) (Λ − γ Ω) / Ω , (E9) c = (cid:0) γ ΛΩ − Λ Ω − − Ω − Ω (cid:1) / Ω , (E10) c = (3Λ − γ Ω) / Ω , (E11) c = ( − / Ω , (E12) d = (cid:0) Λ + Ω (cid:1) / Ω , (E13) d = (cid:0) Λ + Ω (cid:1) ( − γ Ω − / Ω , (E14) d = (cid:2)(cid:0) + Ω (cid:1) − (cid:0) − γ Λ + Λ + Ω (cid:1)(cid:3) / Ω , (E15) d = [ − Λ − Ω( γ − / Ω . (E16)Equation (E8) is then reduced to ˙ U ˙ V = V + d U + (3 c + d ) U V + O (cid:0) U , U V, U V , U V , U V , V (cid:1) (E17)by the following near-identity transformation ˇ u ˇ v = UV + (2 c + d ) U + ( c + d ) U V + c U V − c U + d U V + d U V . (E18)The signs of d and d (cid:48) := 3 c + d = − γ (cid:0) Λ + Ω (cid:1) / Ω determine the types of the codimension-one bifurcation that occurs near the codimension-two bifurcation point[37]. The heteroclinic bifur-cation as well as the Hopf and the pitchfork bifurcations occurs for d d (cid:48) <
0, while the bifurcationinvolving a pair of homoclinic orbit occurs for d d (cid:48) > >
0, andthe latter is the case for Ω < >
0. In practice, we numerically observed it for some parameter setseven when Ω > Appendix F: Numerical search for limit-cycle solutions of Eq. (8)
While the bifurcation analyses in Sec. III is comprehensive in regard to the local bifurcations,some global bifurcations such as saddle–node bifurcation of periodic orbits are not considered there.These global bifurcations might create or destroy stable limit cycles, affecting the stability regionof the synchronously oscillating state in Fig. 1.Thus, we numerically verified that the stability boundary of the synchronously oscillating stateconsists of the Hopf, SNIC, and heteroclinic bifurcation curves obtained in Sec. III as detailedbelow. For each parameter set, 100 initial conditions for Re ( r ) and Im ( r ) are drawn from uniformdistribution on the unit disk, and the type of the attractor to which each orbit converges is detected.In the black region in Fig. 8, at least one orbit converges to a limit cycle. The edge of the blackregion agrees with the bifurcation curves that are inferred to form the stability boundary of thecollective oscillation.The type of the attractor is classified according to the following criteria. (i) Let ∆ r t be the normof the variation of the orbit between time t and t + ∆ t , where ∆ t = 0 .
01. The orbit is considered tobe converged to a fixed point if the phase point is moving slowly (∆ r t , ∆ r t − ∆ t , ∆ r t − t < − )and is slowing exponentially ( | ∆ r t / ∆ r t − ∆ t − ∆ r t − ∆ t / ∆ r t − t | < − ). (ii) Let x max and x min bethe maximum and the minimum, respectively, of Re ( r ) during an interval (10000 < t < y n be the n th intersection of the line Re ( r ) = ( x max + x min ) / r ) >
0. The orbit is considered to be converged to a limit cycle when | y n − y n − | < − . With this criteria, we conclude that every orbit is converged to either a fixedpoint or a limit cycle. [1] A Pikovsky, M Rosenblum, and J Kurths. Synchronization: A Universal Concept in Nonlinear Sciences .2001.[2] Leon Glass. Synchronization and rhythmic processes in physiology.
Nature , 410(6825):277–284, 2001.[3] Alex Arenas, Albert D´ıaz-Guilera, Jurgen Kurths, Yamir Moreno, and Changsong Zhou. Synchroniza-tion in complex networks.
Phys. Rep. , 469(3):93 – 153, 2008.[4] A T Winfree.
The Geometry of Biological Time . Springer, New York, 2nd edition, 2001.[5] Y Kuramoto.
Chemical Oscillations, Waves, and Turbulence . Springer, New York, 1984.[6] Adilson E Motter, Seth A Myers, Marian Anghel, and Takashi Nishikawa. Spontaneous synchrony inpower-grid networks.
Nat. Phys. , 9(3):191–197, 2013.[7] J T Enright. Temporal precision in circadian systems: a reliable neuronal clock from unreliable com- FIG. 8. (Color online) The stable region of the limit cycle solution. The region where at least one stablelimit cycle exists is filled with black. The bifurcation curves obtained in Sec. III are also plotted. The Hopf,SNIC, and heteroclinic bifurcation curves agree with the boundary of the black region.ponents?
Science (80-. ). , 209(4464):1542–1545, sep 1980.[8] Erik D Herzog, Sara J Aton, Rika Numano, Yoshiyuki Sakaki, and Hajime Tei. Temporal Precision inthe Mammalian Circadian System: A Reliable Clock from Less Reliable Neurons.
J. Biol. Rhythms ,19(1):35–46, feb 2004.[9] Daniel J Needleman, Paul H E Tiesinga, and Terrence J Sejnowski. Collective enhancement of precisionin networks of coupled oscillators.
Phys. D Nonlinear Phenom. , 155(3-4):324–336, 2001.[10] Hiroshi Kori, Yoji Kawamura, and Naoki Masuda. Structure of cell networks critically determinesoscillation regularity.
J. Theor. Biol. , 297:61–72, 2012.[11] Auke Jan Ijspeert. Central pattern generators for locomotion control in animals and robots: A review.
Neural Networks , 21(4):642–653, 2008.[12] Stefan L¨ammer, Hiroshi Kori, Karsten Peters, and Dirk Helbing. Decentralised control of materialor traffic flows in networks using phase-synchronisation.
Phys. A Stat. Mech. its Appl. , 363(1):39–47,2006.[13] D R ALEKO and S Djahel. An IoT Enabled Traffic Light Controllers Synchronization Method forRoad Traffic Congestion Mitigation. In , pages 709–715, 2019.[14] Steven H Strogatz, Daniel M Abrams, Allan McRobie, Bruno Eckhardt, and Edward Ott. Crowdsynchrony on the Millennium Bridge.
Nature , 438(7064):43–44, 2005.[15] Constance Hammond, Hagai Bergman, and Peter Brown. Pathological synchronization in Parkinson’sdisease: networks, models and treatments.
Trends Neurosci. , 30(7):357–364, 2007.[16] Matthew M McGregor and Alexandra B Nelson. Circuit Mechanisms of Parkinson’s Disease.
Neuron ,101(6):1042–1056, 2019. [17] Melissa J Armstrong and Michael S Okun. Diagnosis and Treatment of Parkinson Disease: A Review. JAMA , 323(6):548–560, feb 2020.[18] Patricia Limousin and Tom Foltynie. Long-term outcomes of deep brain stimulation in Parkinsondisease.
Nat. Rev. Neurol. , 15(4):234–242, 2019.[19] Jeff M Bronstein, Michele Tagliati, Ron L Alterman, Andres M Lozano, Jens Volkmann, AlessandroStefani, Fay B Horak, Michael S Okun, Kelly D Foote, Paul Krack, Rajesh Pahwa, Jaimie M Henderson,Marwan I Hariz, Roy A Bakay, Ali Rezai, William J Marks Jr, Elena Moro, Jerrold L Vitek, Frances MWeaver, Robert E Gross, and Mahlon R DeLong. Deep Brain Stimulation for Parkinson Disease: AnExpert Consensus and Review of Key Issues.
Arch. Neurol. , 68(2):165, feb 2011.[20] Minseok Kim, Matthias Bertram, Michael Pollmann, Alexander von Oertzen, Alexander S Mikhailov,Harm Hinrich Rotermund, and Gerhard Ertl. Controlling Chemical Turbulence by Global DelayedFeedback: Pattern Formation in Catalytic CO Oxidation on Pt(110).
Science (80-. ). , 292(5520):1357– 1360, may 2001.[21] Michael Rosenblum and Arkady Pikovsky. Delayed feedback control of collective synchrony: An ap-proach to suppression of pathological brain rhythms.
Phys. Rev. E , 70(4):041904, oct 2004.[22] Michael G Rosenblum and Arkady S Pikovsky. Controlling Synchronization in an Ensemble of GloballyCoupled Oscillators.
Phys. Rev. Lett. , 92(11):114102, mar 2004.[23] Oleksandr V Popovych, Christian Hauptmann, and Peter A Tass. Effective Desynchronization byNonlinear Delayed Feedback.
Phys. Rev. Lett. , 94(16):164102, apr 2005.[24] Irmantas Ratas and Kestutis Pyragas. Controlling synchrony in oscillatory networks via an act-and-waitalgorithm.
Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. , 90(3):32914, sep 2014.[25] Natalia Tukhlina, Michael Rosenblum, Arkady Pikovsky, and J¨urgen Kurths. Feedback suppression ofneural synchrony by vanishing stimulation.
Phys. Rev. E , 75(1):011918, jan 2007.[26] Ming Luo, Yongjun Wu, and Jianhua Peng. Washout filter aided mean field feedback desynchronizationin an ensemble of globally coupled neural oscillators.
Biol. Cybern. , 101(3):241–246, 2009.[27] Alessio Franci, Antoine Chaillet, Elena Panteley, and Fran¸coise Lamnabhi-Lagarrigue. Desynchroniza-tion and inhibition of Kuramoto oscillators by scalar mean-field feedback.
Math. Control. Signals, Syst. ,24(1):169–217, 2012.[28] Oleksandr V Popovych, Borys Lysyansky, Michael Rosenblum, Arkady Pikovsky, and Peter ATass. Pulsatile desynchronizing delayed feedback for closed-loop deep brain stimulation.
PLoS One ,12(3):e0173363, mar 2017.[29] Shijie Zhou, Peng Ji, Qing Zhou, Jianfeng Feng, J¨urgen Kurths, and Wei Lin. Adaptive elimination ofsynchronization in coupled oscillator.
New J. Phys. , 19(8):083004, 2017.[30] Aneta Koseska, Evgeny Volkov, and J¨urgen Kurths. Oscillation quenching mechanisms: Amplitude vs.oscillation death.
Phys. Rep. , 531(4):173–199, 2013.[31] Edward Ott and Thomas M Antonsen. Low dimensional behavior of large systems of globally coupledoscillators.
Chaos An Interdiscip. J. Nonlinear Sci. , 18(3):037113, sep 2008. [32] Edward Ott and Thomas M Antonsen. Long time evolution of phase oscillator systems. Chaos AnInterdiscip. J. Nonlinear Sci. , 19(2):023117, may 2009.[33] Peter Ashwin, Stephen Coombes, and Rachel Nicks. Mathematical Frameworks for Oscillatory NetworkDynamics in Neuroscience.
J. Math. Neurosci. , 6(1):2, 2016.[34] Tomislav Stankovski, Tiago Pereira, Peter V. E. McClintock, and Aneta Stefanovska. Coupling func-tions: Universal insights into dynamical interaction mechanisms.
Rev. Mod. Phys. , 89(4):045001, nov2017.[35] Hidetsugu Sakaguchi and Yoshiki Kuramoto. A Soluble Active Rotater Model Showing Phase Transi-tions via Mutual Entertainment.
Prog. Theor. Phys. , 76(3):576–581, sep 1986.[36] Ken H Nagai and Hiroshi Kori. Noise-induced synchronization of a large population of globally couplednonidentical oscillators.
Phys. Rev. E , 81(6):065202, jun 2010.[37] J Guckenheimer and P Holmes.
Nonlinear Oscillations, Dynamical Systems, and Bifurcations of VectorFields . Springer-Verlag, 1983.[38] Bard Ermentrout.
Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT forresearchers and students . SIAM, 2002.[39] Lauren M Childs and Steven H Strogatz. Stability diagram for the forced Kuramoto model.
Chaos AnInterdiscip. J. Nonlinear Sci. , 18(4):43128, dec 2008.[40] H Daido. Scaling behaviour at the onset of mutual entrainment in a population of interacting oscillators.
J. Phys. A. Math. Gen. , 20(10):L629–L636, 1987.[41] Joel T Ariaratnam and Steven H Strogatz. Phase Diagram for the Winfree Model of Coupled NonlinearOscillators.
Phys. Rev. Lett. , 86(19):4278–4281, may 2001.[42] Denis S Goldobin and Arkady Pikovsky. Effects of Delayed Feedback on Kuramoto Transition.
Prog.Theor. Phys. Suppl. , 161:43–52, jan 2006.[43] Oleksandr V Popovych, Christian Hauptmann, and Peter A Tass. Control of Neuronal Synchrony byNonlinear Delayed Feedback.
Biol. Cybern. , 95(1):69–85, 2006.[44] Oleksandr V Popovych and Peter A Tass. Multisite Delayed Feedback for Electrical Brain Stimulation, 2018.[45] Xiaohan Zhang and Shenquan Liu. Nonlinear delayed feedback control of synchronization in an exci-tatory–inhibitory coupled neuronal network.
Nonlinear Dyn. , 96(4):2509–2522, 2019.[46] Oleksandr V Popovych, Christian Hauptmann, and Peter A Tass. Impact of Nonlinear Delayed Feed-back on Synchronized Oscillators.
J. Biol. Phys. , 34(3):267–279, 2008.[47] Istv´an Z Kiss, Craig G Rusin, Hiroshi Kori, and John L Hudson. Engineering Complex DynamicalStructures: Sequential Patterns and Desynchronization.
Science (80-. ). , 316(5833):1886 – 1889, jun2007.[48] Hiroshi Kori, Craig G Rusin, Istv´an Z Kiss, and John L Hudson. Synchronization engineering: The-oretical framework and application to dynamical clustering.
Chaos An Interdiscip. J. Nonlinear Sci. ,18(2):026111, jun 2008. [49] Daniel R Merrill, Marom Bikson, and John G R Jefferys. Electrical stimulation of excitable tissue:design of efficacious and safe protocols. J. Neurosci. Methods , 141(2):171–198, 2005.[50] Istv´an Z Kiss, Yumei Zhai, and John L Hudson. Predicting Mutual Entrainment of Oscillators withExperiment-Based Phase Models.
Phys. Rev. Lett. , 94(24):248301, jun 2005.[51] Hiroshi Kori, Yoshiki Kuramoto, Swati Jain, Istv´an Z Kiss, and John L Hudson. Clustering in globallycoupled oscillators near a Hopf bifurcation: Theory and experiments.
Phys. Rev. E , 89(6):62906, jun2014.[52] D Iatsenko, P V E McClintock, and A Stefanovska. Extraction of instantaneous frequencies from ridgesin time–frequency representations of signals.
Signal Processing , 125:290–303, 2016.[53] Yasuaki Kobayashi and Hiroshi Kori. Design principle of multi-cluster and desynchronized states inoscillatory media via nonlinear global feedback.
New J. Phys. , 11(3):33018, 2009.[54] Arthur T Winfree. Biological rhythms and the behavior of populations of coupled oscillators.