Stability of entrainment of a continuum of coupled oscillators
aa r X i v : . [ n li n . AO ] S e p Stability of entrainment of a continuum of coupled oscillators
Jordan Snyder, a) Anatoly Zlotnik, and Aric Hagberg Department of Mathematics, University of California, Davis,CA 95616 Applied Mathematics and Plasma Physics (T-5), Theoretical Division,Los Alamos National Laboratory, Los Alamos, NM 87544
Complex natural and engineered systems are ubiquitous and their behavior is challengingto characterize and control. We examine the design of the entrainment process for anuncountably infinite collection of coupled phase oscillators that are all subject to the sameperiodic driving signal. In the absence of coupling, an appropriately designed input canresult in each oscillator attaining the frequency of the driving signal, with a phase offsetdetermined by its natural frequency. We consider a special case of interacting oscillators inwhich the coupling tends to destabilize the phase configuration to which the driving signalwould send the collection in the absence of coupling. In this setting we derive stabilityresults that characterize the trade-off between the effects of driving and coupling, andcompare these results to the well-known Kuramoto model of a collection of free-runningcoupled oscillators.Keywords: synchronization, control, nonlinear dynamics a) Electronic mail: [email protected] wo well-understood approaches can be applied to impose coherent behavior in a diversepopulation of dynamical systems: the “top-down” approach of applying a common driv-ing signal, and the “bottom-up” approach of imposing pairwise coupling. While these ap-proaches yield similar behaviors, their precise characteristics can put them in opposition. Inthis article we study a situation that highlights both the synergy and tension that can existbetween driving and coupling in collections of oscillators.I. INTRODUCTION A growing variety of collective dynamic behaviors in networked periodic phenomena are ob-served across disciplines, with many well-known examples such as circadian cycles, injection-locked semiconductors, and emerging applications including battery charging cycles and neuralinformation processing. Analysis, design, and control of interacting collections of rhythmic oroscillatory processes whose scale or complexity is beyond the scope of classical dynamical sys-tems theory requires new mathematical frameworks. Developing continuum approximations tooscillatory phenomena over large-scale networks will create paths towards practical solutions tomotivating applications.One basic class of phenomena that is of interest in chronobiology , electrochemistry ,neuroscience , and power grid engineering is the formation of coherent behavior in collectionsof interacting units. Such coherence can be imposed from outside, or can arise through the intrin-sic interactions themselves. There is a long history of studying the emergence of coherent motionin oscillators using phase model representations, dating back to Winfree and Kuramoto . Thecanonical examples that were posed in these early studies have been widely examined in subse-quent work , because they exhibit a rich phenomenology while admitting beautiful mathematicaldescriptions within an extensive range of analytical settings. While original studies focused onmutual entrainment , in which coherent motion arises purely from interactions between individualunits, more recent studies have investigated the effect of externally-imposed coherence in the formof external driving.A pioneering study of forced, coupled oscillations was performed by Sakaguchi , who consid-ered an infinite, heterogeneous, population of oscillators subject to global sinusoidal coupling anduniform sinusoidal forcing. By deriving a self-consistent equation for the order parameter measur-2ng phase alignment, Sakaguchi was able to predict transitions between regimes of incoherence,mutual entrainment, and forced entrainment. These transitions were subsequently investigated inmore detail by Antonsen and collaborators using a linear stability approach . The detailed na-ture of the bifurcations remained elusive, and suggested an underlying two-dimensional structurewhich had yet to be exploited. This two-dimensional structure was indeed discovered by Ott andAntonsen in their seminal work , which uncovered a particular low-dimensional manifold thatcaptures much of the asymptotic behavior of a wide family of models of coupled phase oscillators;in particular, the forced Kuramoto model was identified as a possible application of this dimensionreduction. Subsequent work has shown that, under certain mild conditions, this manifold is glob-ally attractive . This reduction represents an enormous simplification, because in many casesit permits closed-form evolution equations for the synchrony order parameter directly. Buildingon the framework established by Ott and Antonsen , Childs and Strogatz were able study thedynamics of the forced Kuramoto model on the two-dimensional attractive manifold, and found acomplete picture of a system’s bifurcation structure. It should be noted that studies of the effect offorcing in this context have almost exclusively considered a sinusoidal forcing function, due to theanalytical tractability of the resulting model. Recently, complementary work has been done thatexamines the role of random forcing applied to a population of sinusoidally-coupled oscillators .Beyond characterizing the phenomenology of natural and engineered complex oscillating sys-tems, emerging applications in neural systems, electrochemistry, and power grid engineering re-quire new capabilities to control and manipulate the behavior of such phenomena. Indeed, theability to control a system is the ultimate validation of our understanding of its behavior. For os-cillating systems, a general picture of frequency modulation by external forcing was first laid out in1946 by Adler , who derived equations describing the amount by which an external drive signalcan shift an oscillator’s frequency and amplitude. The idea of “injection locking” has since beenof major importance in many fields of engineering . One prediction made by Adler was thatan oscillator driven at a frequency different from its own may lock to the driving frequency, andexhibit a phase shift relative to the drive signal which is determined by its natural frequency. Forthe simple case considered originally, this function is sinusoidal. However, for a general forcingfunction and a general phase oscillator higher harmonics may be present, as seen in experimentsand derived analytically . The general framework of using periodic forcing signals to controlthe entrainment of nonlinear oscillating systems has been exploited to explore energy- and time-optimal control strategies for entrainment of one or more phase oscillators . The effect of3oupling on the efficacy of these control strategies remains unexplored.A challenge in specifying the forcing input to control collections of coupled oscillators is thatthey are underactuated; the entire collection of similar dynamical systems with possibly compli-cated individual behavior must be controlled using a small number of inputs. To overcome thischallenge we observe that the entrained or coherent state of a controlled collection of oscillatorsis characterized not only by synchronization to a forcing frequency, but also by the distribution ofsubsystems on the neighborhood of a nominal periodic orbit. For a finite collection, it is possibleto construct a forcing signal to achieve precise control of the relative phases of an ensemble ofstructurally similar oscillators with slight heterogeneity in frequencies . With the understandingthat such “phase-selective control” is possible for small, finite collections, we examine how themathematical framework can be extended to continuum systems. Further, we examine the effectof coupling between subsystems, which tends to drive phase differences to zero.In this paper we explore a continuum approximation of a very large collection of coupledoscillators subject to a common periodic (but non-sinusoidal) forcing, so that both coupling andforcing influence the collective behavior. Specifically, we consider a situation in which the forcingdrives individual phases to be maximally different (in a certain precise sense), while the couplingtends to align the phases. To quantify the trade-off between these two effects, we compute, asa function of the coupling strength, the asymptotic stability of a fixed point in which the phasesshow no global alignment. By finding the critical coupling strength above which this fixed pointis unstable, we demonstrate that mutual synchronization of entrained coupled oscillators occurs before mutual synchronization of unforced coupled oscillators, despite the imposed diversity ofphases. Moreover, numerical experiments confirm that the external forcing has facilitated phasealignment which is greater than that in the unforced case. Our results demonstrate that measuringonly phase alignment is bound to miss important information about the global organization of apopulation of oscillators. II. PRELIMINARIESA. Entrainment of oscillators
We next describe how a heterogeneous population of oscillators can be caused to move at asingle frequency by application of a suitable forcing function. Mathematical details can be found4n several standard references .As our basic model of an oscillator, we take a phase model , first popularized by Winfree . For i = , . . . , N , the i th oscillator is described by the ODE˙ ψ i = ω i + Z ( ψ i ) u , (1)where ψ i ∈ [ , π ) is the phase, “ ˙ ” denotes the derivative with respect to time, ω i ∈ R is thenatural frequency, u = u ( t ) is an external forcing, and Z ( ψ i ) is known as the phase response curve ,or PRC. The PRC determines the change in phase resulting from an infinitesimal external forceapplied at a given phase on the limit cycle . The equation (1) can be derived by consideringthe lowest-order approximation of the effect of an external force acting on a system near a stablelimit cycle , and in this sense is representative of a wide class of forced periodic motions.A standard approach to analyzing entrainment is to take u ( t ) = v ( Ω t ) , where v has period 2 π sothat Ω denotes the (angular) frequency of the driving signal. If Ω is not too far from ω i , we supposethat ψ i will behave as Ω t , plus a slowly-varying phase offset. We formalize this supposition bymaking the change of coordinates ψ i = Ω t + φ i , where φ i now represents the phase offset. In the φ i coordinate system, the dynamics now read˙ φ i = ∆ω i + Z ( φ i + Ω t ) v ( Ω t ) , (2)where we have introduced the frequency detuning ∆ω i ≡ ω i − Ω . Finally it is possible to approxi-mate (2) by the time-averaged system ˙ ϕ i = ∆ω i + Λ v ( ϕ i ) , (3)where we have introduced the interaction function Λ v ( ϕ ) = π π Z Z ( ϕ + θ ) v ( θ ) d θ , (4)in the sense that there exists a change of variables ϕ i = φ i + h ( ϕ i , φ i ) that maps solutions of (2) tothose of (3).If the frequency detunings { ∆ω i } are such that (3) has a stable fixed point solution for all i = , . . . , N , then the phases of all oscillators will be constant in the moving reference framewith frequency Ω . In other words, the entire population can be entrained by the driving signal u = v ( Ω t ) . 5espite having equal frequencies, the oscillators will, in general, have different phases, sincethe solutions to the fixed point equation ∆ω i + Λ v ( ϕ i ) = ∆ω i . This fact canbe exploited to design a forcing function that elicits frequency locking with a known distributionof phases, irrespective of initial conditions . B. The Kuramoto Model
To frame our study of phase coupling, we discuss some standard methods and results relatingto synchronization of phase oscillators. Kuramoto introduced a model of the form˙ ϕ i = ω i + KN N ∑ j = sin ( ϕ j − ϕ i ) , (5)which was derived as a “simplest” model for a collection of self-sustained linearly-coupledoscillators . Here { ϕ i } are the phases of N oscillators, { ω i } their natural frequencies (whichwe allow to take any real values), and K > ϕ i = ω i + KR sin ( Φ − ϕ i ) , (6)where we have used the synchrony R ∈ [ , ] , and the average phase Φ ∈ [ , π ) , first introducedby Kuramoto and defined by the formula Re i Φ = N N ∑ j = e i ϕ j . (7)In this sense, this form of coupling is mean-field in character, as each phase feels a force deter-mined by an average over the entire population.The key features of this model are1. The oscillators have differing intrinsic frequency: ω i = ω j ,2. The coupling tends to drive phases towards the mean (provided K >
0, which we assumethroughout).These two features are at odds with each other, and they undergo a trade-off at a critical valueof the coupling strength, K = K unf c (where we use the superscript ”unf” to emphasize that thisis the critical coupling strength in the unforced case). If K < K unf c the population of oscillators6oes not show global alignment towards any particular phase, while for K > K unf c , this situationbreaks down and a subset of the oscillators attains the same frequency and group together in phase,establishing a preferred direction and a nonzero value of the synchrony R .To make these statements precise it is useful to consider a mean-field approximation . Wesuppose that the population of oscillators is large enough that averaging over this population iswell approximated by averaging over a probability distribution that describes the behavior of atypical oscillator. General background on the mean field Kuramoto model can be found in variousreview articles .The main result we quote from the extensive body of literature on the Kuramoto model is that inthe limit of N → ∞ , if the oscillators’ natural frequencies are drawn at random from a probabilitydistribution having density g ( ω ) , unimodal and symmetric about zero, then the critical couplingstrength described above is given by K unf c = π g ( ) . (8)The expression (8) can be taken as a precise quantification of the trade-off between intrinsicdisorder ( g ( ) ) and coupling ( K ). The possibly surprising fact that K unf c depends only on the valueof g at the center of the distribution, and no other features of this distribution, is because the firstoscillators to synchronize are those whose natural frequencies lie at the center of the distribution.The rest of the density g then determines the growth of R with K > K unf c .In what follows, we will define a new model, show that it exhibits behavior that is qualitativelysimilar to that of the Kuramoto model, and find the location of the corresponding critical point.The expression (8) will serve as reference to interpret our results. III. MODEL FOR FORCING OF COUPLED OSCILLATIONSA. Finite N We now formulate a model of a population of oscillators that exhibits both frequency alignmentby broadcast forcing and phase alignment by attractive coupling. Many similar models have beendeveloped , and our present formulation aims to augment the rich existing literature.In general, we can consider each oscillator to respond to external forcing according to onephase response curve, and to respond to forcing from its neighboring oscillators according to7nother phase response curve. That is,˙ ψ i = ω i + Z e ( ψ i ) u ( t ) + KN N ∑ j = Z c ( ψ i ) f ( ψ j ) (9)where Z e is the PRC for external forcing, Z c is the PRC for coupling, and f () describes the forcean oscillator exerts on its neighbors as a function of its phase. The prefactor K / N allows us toadjust the coupling strength K in a way that allows comparison between different values of N .Assuming, as before, that u ( t ) = v ( Ω t ) with v having period 2 π , we move into a rotatingreference frame with frequency Ω and average over one period of the driving signal, obtaining theaveraged equations ˙ ϕ i = ∆ω i + Λ v ( ϕ i ) + KN N ∑ j = g ( ϕ j − ϕ i ) (10)where ϕ i , ∆ω i , and Λ v are defined as before (3) and g ( ∆ϕ ) = ( π ) − R π Z c ( θ + ∆ϕ ) f ( θ ) d θ .Clearly, many different systems may be defined in this form given appropriate choices for Z e , Z c , v , and f . In order to exhibit the qualitative features of phase dispersion caused by external forc-ing combined with phase alignment caused by coupling, while retaining tractability, we assumethat Z c and f are such that g ( ∆ϕ ) = sin ( ∆ϕ ) .Hence, we take a model of the form˙ ϕ i = ∆ω i + Λ v ( ϕ i ) + KN N ∑ j = sin ( ϕ j − ϕ i ) , (11)which can also be written ˙ ϕ i = ∆ω i + Λ v ( ϕ i ) + KR sin ( Φ − ϕ i ) , (12)with R and Φ defined as in (7).As a first step, we choose { ∆ω i } and Λ v such that all oscillators can be entrained individually,but the resulting phase offsets are as far as possible from alignment. This can be achieved bysetting ∆ω i = iN − , (13)and Λ v ( ϕ ) = − ϕπ , ϕ ∈ ( − π , π ] . (14)We refer to the function defined in (14) as the sawtooth interaction function , or just sawtooth , asit has a sawtooth shape when plotted on R . 8he standard unforced Kuramoto model with this choice of natural frequencies has been re-cently studied by Ottino-L¨offler and Strogatz , who found the asymptotic behavior of the lockingthreshold as N → ∞ , in agreement with results in the thermodynamic limit obtained earlier byPaz´o . These results will serve as a reference to put our findings in context. For now, we returnto the forced case.In the absence of coupling ( K = i th oscillator will be driven to a phase offset ϕ ∗ i definedby ∆ω i + Λ v ( ϕ ∗ i ) = = ⇒ ϕ ∗ i = π∆ω i = π iN − π . (15)A straightforward calculation shows that for this phase configuration, the synchrony is R =
0. Forthis reason, we refer to this fixed point as the desynchronized state . Another term used to describesuch a state is splay state . The point ϕ ∗ = ( ϕ ∗ i ) ∈ ( − π , π ] N is a fixed point of the dynamics (11)for any value of coupling strength K .In this respect, the situation is similar to the incoherent state discussed for the Kuramoto modelin Section II B, with the key difference that in this case, all oscillators have attained identicalfrequency locking to the forcing input. We proceed to study the asymptotic stability of this fixedpoint as a function of K , and obtain a critical coupling strength K c analogous to K unf c as defined in(8). B. The N → ∞ limit Next we introduce a thermodynamic limit of the model (11), and the fixed point correspondingto that defined in (15).We replace our population of oscillators, formerly a collection of N individual oscillators withnatural frequencies evenly spaced from − [ − , ] .Because our state of interest for finite N is such that each oscillator’s phase is fixed at a valuedetermined by its natural frequency, we describe the state of our infinite system by a function ϕ ( ω ) that gives the phase of any oscillator having natural frequency ω . As the system evolves the wholefunction ϕ ( ω ) will change in time, but for visual clarity we omit writing the time-dependenceexplicitly when discussing fixed points. This sort of formulation is used, for example, by Mirolloand Strogatz , except that oscillators are indexed by their “pinning phase” rather than their naturalfrequency. We describe this work in more detail in Section IV D.9o determine fixed points, we must establish the dynamics in the appropriate continuum setting.The intrinsic dynamics and effects of forcing remain the same, so we only need to concern our-selves with the coupling term. For finite N , we simply had an average over the population, and inthe infinite setting, we use a mean-field approach to say that averaging over the infinite populationis equivalent to averaging over the distribution of natural frequencies . Our infinite-dimensionaldynamics are ∂ t ϕ ( ω ) = ω + Λ v ( ϕ ( ω )) + K Z R g ( ω ′ ) sin ( ϕ ( ω ′ ) − ϕ ( ω )) d ω ′ , (16)where g is the density of the distribution of natural frequencies. These dynamics can be rewrittenin the form ∂ t ϕ ( ω ) = ω + Λ v ( ϕ ( ω )) + KR sin ( Φ − ϕ ( ω )) , (17)where R and Φ are the synchrony and average phase, defined for the infinite system as Re i Φ = Z R g ( ω ) e i ϕ ( ω ) d ω . (18)Using the sawtooth interaction function introduced above (see (14)), and g ( ω ) = / ω ∈ [ − , ] and 0 elsewhere, the fixed point condition for ϕ now reads0 = ω − ϕ ( ω ) π + K Z −
12 sin ( ϕ ( ω ′ ) − ϕ ( ω )) d ω ′ , (19)A straightforward calculation shows that the function ϕ ( ω ) = πω satisfies the condition (19).Note that this is precisely the infinite- N analog of the finite- N fixed point defined in (15). In whatfollows, we perform a linear stability analysis, finding the coupling strength K c at which this statebecomes unstable. IV. STABILITY ANALYSIS OF THE ENTRAINMENT PHASE DISTRIBUTIONA. Finite N We now analyze the stability of the fixed point ϕ ∗ = ( ϕ ∗ i ) ∈ ( − π , π ] N as defined in equation(15).Asymptotic stability of ϕ ∗ is controlled by the spectrum σ ( J ) of the Jacobian J of the right-hand side of (11) with respect to ϕ , evaluated at ϕ ∗ . If every element of σ ( J ) has negative real10art, then ϕ ∗ is an asymptotically stable fixed point, and if any element of σ ( J ) has positive realpart, then ϕ ∗ is unstable . The matrix elements of J are J i j = Λ ′ v ( ϕ ∗ i ) − KN ∑ k = i cos ( ϕ ∗ i − ϕ ∗ k ) ! δ i j + ( − δ i j ) KN cos (cid:0) ϕ ∗ i − ϕ ∗ j (cid:1) , (20)where δ i j is the Kronecker delta. A straightforward calculation (see Appendix A 1) shows that σ ( J ) = (cid:26) − π , − π + K (cid:27) , (21)so the desynchronized state has a critical point K c = / π and is linearly stable when K < / π , andlinearly unstable for K > / π . B. The N → ∞ limit Finally we will perform a linear stability analysis of the desynchronized fixed point of theinfinite- N model (16). For details of the calculation presented below, see Appendix A 2.To obtain a linearization of the dynamics near the fixed point ϕ ∗ ( ω ) = πω , we consider aninfinitesimal perturbation, ϕ ( ω ) = ϕ ∗ ( ω ) + εη ( ω ) , (22)where 0 < ε ≪ η : [ − , ] → R is a function which we take to be bounded and measurable.Inserting this form into (16) and collecting terms by order of ε yields O ( ε ) : ∂ t ϕ ∗ ( ω ) = ω − ϕ ∗ ( ω ) π + K Z −
12 sin ( ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω )) d ω ′ O ( ε ) : ∂ t η ( ω ) = − πη ( ω ) + K Z −
12 cos ( ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω )) η ( ω ′ ) d ω ′ . (23)As expected, the O ( ε ) equation holds by the fact that ϕ ∗ is a fixed point, and the O ( ε ) gives thetime evolution of small perturbations around ϕ ∗ .We can diagonalize the dynamics (23) by writing η as a Fourier series, η ( ω ) = ∑ k ∈ Z c k ( t ) e ik πω . (24)Inserting the form (24) into the O ( ε ) equation (23), we find that ∂ t c k = (cid:18) − π + δ | k | , K (cid:19) c k . (25)11ll Fourier components of the perturbation η except for the first decay exponentially in timewith a rate 1 / π , while the first Fourier component will grow or shrink with time, depending onthe sign of − / π + K /
2. Specifically, if K < / π , then the first Fourier mode also decays in time,while if K > / π , the first Fourier mode grows in time, and the fixed point ϕ ∗ is unstable. Hencewe have, as in the finite- N case, the critical coupling strength K c = / π . C. Interpretation
In both the finite- and infinite-dimensional versions of our model, we have found that nonzerosynchrony spontaneously develops as the coupling strength K exceeds K c = / π . We contrast thisresult with that for the corresponding unforced model,˙ ϕ i = ω i + KN N ∑ j = sin ( ϕ j − ϕ i ) , (26)where ω i = π i / N − π . While the standard result (8) does not directly apply in this case, sincethe uniform density is not unimodal, it has been established by Paz´o that the synchronizationtransition does in fact occur at K unf c = / π = / ( π g ( )) , which is twice the value at which theforced model begins to show nonzero synchrony. This result can be considered surprising, giventhat we have taken a forcing term, (14), that was designed specifically to drive the system to a stateof zero synchrony.The situation becomes clearer if we compare the desynchronized state present in the forcedmodel to the incoherent state in the unforced model. The desynchronized state, defined by ϕ ( ω ) = πω , has zero synchrony as measured by the order parameter R . However, it has the property thatevery oscillator moves at equal frequency. This is in contrast with the incoherent state of theunforced Kuramoto model, in which each oscillator moves at its own natural frequency. Hence, inthe sense of frequencies, the desynchronized state is far more organized than the incoherent state,although this fact is missed by the synchrony parameter R , which only measures instantaneousalignment of phases.To understand the role of frequency alignment in establishing phase alignment, it is instructiveto consider again the standard unforced Kuramoto model. As we have already quoted (8), thecritical coupling strength is K unf c = / π g ( ) , where g is the density of the distribution of naturalfrequencies. Intuitively, this expression captures the trade-off between disorder in the natural fre-quencies and the ordering influence of coupling; the tighter the distribution of natural frequencies,12he larger g ( ) , and the smaller K unf c . In other words, the coupling strength must be large enoughto overcome the diversity of natural frequencies in order to bring about a preferred phase.In the desynchronized state of the forced model, the oscillators move with a single frequency.Hence, there is no disorder to be overcome by the coupling. All that keeps the system in thedesynchronized state is the forcing, which appears as the eigenvalue − / π in the spectrum of theJacobian. The second eigenvalue, − / π + K /
2, directly captures the trade-off between the drivingand the coupling, showing that the stability of the entrained state is the only force that needs to becountered by coupling.
D. On the Relation to Previous Work
Finally, we discuss the relationship of the present model to previous work on models of glob-ally coupled oscillators subject to common forcing. The existing literature has focused almostexclusively on sinusoidal forcing , owing to the analytical progress that this assumptionallows.One such model was discussed by Ott and Antonsen as a possible application of the powerfuldimension reduction known as the Ott-Antonsen (OA) ansatz. While it is the case that the OAansatz can describe the fixed point that we consider, the dynamics away from the fixed point donot leave the OA manifold invariant, precisely because the sawtooth forcing function we consideris not sinusoidal.Another system much more closely similar to ours is the “random pinning” model studied byMirollo and Strogatz . The random pinning model consists of a system of N spins, with each onepinned by an anonymous driving force to a particular (randomly chosen) phase. In explicit terms,the dynamics are ˙ ϕ i = sin ( α i − ϕ i ) + KN N ∑ j = sin ( ϕ j − ϕ i ) , (27)where { α i } are random quantities sampled from the uniform distribution on the unit circle. Theonly difference between this equation and the one that we study is the term ω i − ϕ i / π is replaced bysin ( α i − ϕ i ) . It remains the case that in the absence of coupling, each oscillator evolves accordingto an autonomous ODE on the unit circle with one stable fixed point, and that the state in whicheach oscillator is at its individual fixed point has R ≈ α . Owing to the regularity of the sine function, it is possibleto obtain precise analytical results on the existence, number, and stability of fixed points. Ourformulation is not amenable to the same analysis, for the reason that the sawtooth forcing functionwe consider has infinitely many Fourier modes. V. NUMERICAL SIMULATIONS
Here we present some numerical studies of the dynamical system (11), which confirm thebifurcation at K = K c = / π and illuminate the system’s behavior away from the bifurcation point.To serve as reference, we also present data from numerical solution of the system in the absenceof forcing which is the Kuramoto model with evenly spaced natural frequencies.As we can see in Fig. 1, the synchrony R achieved at any value of the coupling strength K greaterthan 2 / π is greater in the forced case than in the unforced case. This confirms the conclusion thatentrainment by broadcast periodic forcing has brought the system closer to synchrony, as measuredby the order parameter R .In the numerical simulations above, we find a sharp increase in the steady-state value of R as afunction of coupling strength K . To obtain a deeper understanding of the nature of this transitionand of the R > N values, numerical continuation ofthe R > K using the numerical continuationsoftware AUTO .To perform numerical continuation with AUTO, it is first necessary to locate an attractor (inour case, a fixed point) on the branch of interest. For each N from 3 to 100, this was accomplishedby numerical integration of (11) until stationarity with K = .
7. This value of K was chosen asit is greater than K c = / π , and was observed to lead to an R > R × [ − π , π ] N ∋ ( K , ϕ ) , searching in the negative K direction from the user-supplied fixed point. AUTO equation and constants files, including initial fixed point locations for3 ≤ N ≤ N = . . . R > K = K s ( N ) < K c = / π . The unstable portion of this branch exists for all14 IG. 1. Synchronization R vs. coupling strength K . In the unforced case (dashed line) the synchronizationthreshold is K unf c = / π . When forcing is added to drive the system to a splay state of equally distributedphase angles it synchronizes at a lower coupling strength K c = / π . The data were generated from asimulation of N =
50 oscillators, starting from random initial conditions. For the forced case, integrationwas carried out until the system was determined to be at a fixed point. For the unforced case, integrationwas carried out until the system was determined to be in a statistically steady state. K ∈ [ K s ( N ) , K c ] , and meets the R = K = K c . Fig. 2 clearly shows, for N = , ,
20, the existence of a bistable region [ K s ( N ) , K c ] , implying that hysteresis is possible upon slow variation of K .Moreover, we find that the shape of the bifurcation diagram in the bistable region obeys a strongregularity across different values of N . In particular, the width of the bistable region, namely K c − K s ( N ) , follows a power-law scaling with N , with exponent − .
67. Additionally, the value of R at the saddle-node point, which we denote R s ( N ) , is observed to approach a value R c = / − .
29 (see Fig. 3). We expect, therefore, that theinfinite- N system will exhibit a jump bifurcation at K = / π with a height (as measured by R ) of1 /
2, but without a hysteresis loop. 15
IG. 2. Bifurcation diagrams for the finite- N system showing the bistable region as it depends on N . Solidlines indicate stable fixed points; dotted lines, unstable. Data generated using AUTO software . The situation is similar to that investigated by Paz´o , who found the locking threshold for the(unforced) Kuramoto model with evenly spaced natural frequencies. In contrast with the typicallyconsidered case in which the density g of the natural frequency distribution has g ′′ ( ) <
0, leadingto a continuous synchronization transition , the uniform distribution has g ′′ ( ) =
0, and the tran-sition is discontinuous. Precise results for the height of the jump, R unf c , and the scaling of R − R unf c for K > K unf c were derived using a self-consistent approach .Correspondingly, Paz´o found, in the finite- N system, a phenomenon of global frequency align-ment for K below the infinite- N critical point, K unf c . Specifically, it happens that as couplingstrength is increased, oscillators with nearby frequencies lock to each other, forming clumps,which then merge as K is further increased. The final merge occurs at K = K s ( N ) , which ap-proaches K unf c from below as N → ∞ , according to K unf c − K s ( N ) ∼ N − µ with µ ≈ /
2. We shouldnote that for finite N , the transition in the unforced case is not hysteretic, as it is in the forced case.16 I. CONCLUSIONS
We have explored, using an idealized model, the interplay between two ways in which a pop-ulation of phase oscillators may be caused to behave coherently: common periodic forcing andattractive coupling. Based on the synchrony order parameter, R , forcing and coupling can appearto be at odds; the forcing drives R towards zero while the coupling drives R towards one. How-ever, as we demonstrate both analytically and numerically, this view is inherently limited, sincefor K above K c = / π , the forced system exhibits greater phase alignment than the correspondingunforced system. An intuitive explanation for this mismatch is that the parameter R measures onlyphase alignment, and is prone to miss the necessary precondition of frequency alignment.Though we have gained considerable intuition from the results already obtained, there is morework to be done. First, we are still lacking analytical understanding of the upper branch of solu-tions, which would include an expression for the height of the jump and the scaling of R with K above the jump (see Fig. 2).Another set of questions involves the (in)feasibility of the sawtooth interaction function Λ v (de-fined as Λ v ( ϕ ) = − ϕ / π for ϕ ∈ ( − π , π ] (14)). A simple argument reveals that for any integrableforcing waveform v ∈ L ( , π ) , the corresponding interaction function Λ v will be continuous on S , a condition which the sawtooth does not satisfy. It remains unexplored to what degree theresults presented here may be approximated by interaction functions that approximate a sawtooth.One could investigate the scaling of dynamical properties with the energy of the input signal used.Techniques for analysis and control of entrainment processes can be used to examine and evenmanipulate numerous processes in biology . In addition to developing an initial mathematicalframework for characterizing stability of coherence phase structures in a continuum of interactingoscillators, our work presents a potential path towards addressing a compelling biological applica-tion. Specifically, although some disagreement about the nature and phenomenology of epilepsyexists in the neuroscience literature , studies in animal models have indicated that control of syn-chronization of neural dynamics can mitigate epileptiform activity . It is understood that neuralstimulation is an underactuated system because one or a few electrodes are used to control themean field of a very large collection of interacting neurons, which for practical purposes may beapproximated by a continuum . The ability to characterize the stability of phase decoherencein continuum models of general coupled oscillators could determine the possibility of develop-ing effective desynchronizing stimuli for treatment of epilepsy. The criterion that is derived in17ection IV C and validated by numerical experiments in Section V could in principle be testedexperimentally . ACKNOWLEDGMENTS
This work was supported by the Laboratory Directed Research and Development programthrough the Center for Nonlinear Studies at Los Alamos National Laboratory under Departmentof Energy Contract No. DE-AC52-06NA25396. We gratefully acknowledge support from the USArmy Research Office MURI award W911NF-13-1-0340 and Cooperative Agreement W911NF-09-2-0053. The authors thank the anonymous reviewers for the constructive comments, and J.S.acknowledges Raissa D’Souza for useful discussions.
Appendix A: Stability Analysis Calculations1. Spectrum of the Jacobian in Finite Dimensions
Here we calculate the spectrum of the N × N matrix J whose entries are J i j = Λ ′ v ( ϕ ∗ i ) − KN ∑ k = i cos ( ϕ ∗ i − ϕ ∗ k ) ! δ i j + ( − δ i j ) KN cos (cid:0) ϕ ∗ i − ϕ ∗ j (cid:1) , (A1)where δ i j is the Kronecker delta. By symmetry of the phase configuration, the sum in the diagonalterm is independent of i , and can be computed by noticing0 = N ∑ k = cos ( ϕ ∗ i − ϕ ∗ k )= cos ( ϕ ∗ i − ϕ ∗ i ) + ∑ k = i cos ( ϕ ∗ i − ϕ ∗ k ) , (A2)where the first equality in (A2) follows from symmetry ( R = ∑ k = i cos ( ϕ ∗ i − ϕ ∗ k ) = − cos ( ϕ ∗ i − ϕ ∗ i ) = − . (A3)This allows us to write the matrix entries in a simpler form, which will facilitate calculation ofeigenvalues, J i j = (cid:18) Λ ′ v ( ϕ ∗ i ) + KN (cid:19) δ i j + ( − δ i j ) KN cos (cid:0) ϕ ∗ i − ϕ ∗ j (cid:1) = − π δ i j + KN cos ( ϕ ∗ i − ϕ ∗ j ) . (A4)18ote that we have used Λ ′ v ( ϕ ∗ i ) = − / π for all i , and δ i j + ( − δ i j ) cos ( ϕ ∗ i − ϕ ∗ j ) = cos ( ϕ ∗ i − ϕ ∗ j ) for all i , j . We can write this in matrix form as J = K C − π I , (A5)where I is the identity matrix and C is the matrix with entries C i j = N − cos ( ϕ ∗ i − ϕ ∗ j ) . This formmakes it clear (see calculation starting at (A7)) that to find the eigenvalues of J for arbitrary valuesof K , it suffices to find the eigenvalues of C . To do this, we can write the action of C on an arbitraryvector x as ( C x ) i = N ∑ j = C i j x j = N − N ∑ j = cos ( ϕ ∗ i − ϕ ∗ j ) x j = cos ( ϕ ∗ i ) " N − ∑ j cos ( ϕ ∗ j ) x j + sin ( ϕ ∗ i ) " N − ∑ j sin ( ϕ ∗ j ) x j , where we have used the sum angle identity for cosine. The range of C is spanned by the vectors e = ( cos ( ϕ ∗ i )) Ni = and e = ( sin ( ϕ ∗ i )) Ni = . Each of these is in fact an eigenvector with eigenvalue1 /
2, which follows from ( C e ) i = cos ( ϕ ∗ i ) " N − ∑ j cos ( ϕ ∗ j ) cos ( ϕ ∗ j ) + sin ( ϕ ∗ i ) " N − ∑ j sin ( ϕ ∗ j ) cos ( ϕ ∗ j ) = cos ( ϕ ∗ i ) " N − ∑ j cos ( ϕ ∗ j ) = cos ( ϕ ∗ i ) " N − ∑ j + cos ( ϕ ∗ j ) =
12 cos ( ϕ ∗ i ) = ( e ) i , (A6)and similarly for e . Hence e and e are eigenvectors of C with eigenvalue 1 /
2, and all othereigenvalues of C are zero.Finally, we can find the eigenvalues of J for arbitrary K . Notice that λ ∈ σ ( J ) ⇐⇒ det ( J − λ I ) = ⇐⇒ det (cid:18) K C − (cid:18) π + λ (cid:19) I (cid:19) = ⇐⇒ det (cid:18) C − K − (cid:18) π + λ (cid:19) I (cid:19) = ⇐⇒ K − (cid:18) π + λ (cid:19) ∈ σ ( C ) . (A7)Hence the eigenvalues λ of J are of the form λ = − / π + K µ , for µ ∈ σ ( C ) = { , / } . In otherwords, σ ( J ) = (cid:26) − π , − π + K (cid:27) . (A8)19 . Linearization in Infinite Dimensions Here we present the details of linearizing the infinite-dimensional dynamics (16) at the desyn-chronized fixed point ϕ ∗ ( ω ) = πω .First, inserting the form ϕ ( ω ) = ϕ ∗ ( ω ) + εη ( ω ) into equation (16) yields ∂ t ( ϕ ∗ + εη ) = ω − ϕ ∗ ( ω ) + εη ( ω ) π + K Z −
12 sin (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) + ε (cid:0) η ( ω ′ ) − η ( ω ) (cid:1)(cid:1) d ω ′ . Next we expand the sine function in the integrand around the point ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) , and obtainsin ( ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) + ε ( η ( ω ′ ) − η ( ω ))) = sin ( ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ))+ ε cos ( ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ))( η ( ω ′ ) − η ( ω )) + O ( ε ) . From here we can read off the terms of order ε from each side of the equation, and get ∂ t ϕ ∗ = ω − ϕ ∗ ( ω ) π + K Z −
12 sin (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) (cid:1) d ω ′ , which clearly holds, as each side evaluates to zero for all ω ∈ [ − , ] .Next, we gather terms of order ε and obtain (dropping the ε factor from all terms) ∂ t η = − η ( ω ) π + K Z −
12 cos (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) (cid:1) (cid:2) η ( ω ′ ) − η ( ω ) (cid:3) d ω ′ . We can in fact simplify the integral above by noticing that Z −
12 cos (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) (cid:1) (cid:2) η ( ω ′ ) − η ( ω ) (cid:3) d ω ′ = Z −
12 cos (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) (cid:1) η ( ω ′ ) d ω ′ − η ( ω ) Z −
12 cos (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) (cid:1) d ω ′ = Z −
12 cos (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) (cid:1) η ( ω ′ ) d ω ′ , (A9)which follows from the symmetry of the phase configuration ϕ ∗ . We then arrive at the linearizeddynamics as presented in the main text, (23), which we repeat here for completeness, ∂ t η ( ω ) = − πη ( ω ) + K Z −
12 cos ( ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω )) η ( ω ′ ) d ω ′ . (A10)20inally, we demonstrate the diagonalization of the linearized dynamics (A10) in the Fourierbasis. As η is a function on [ − , ] , the appropriate Fourier basis is { e ik πω | k ∈ Z } , so we write η ( ω ) = ∑ k ∈ Z c k ( t ) e ik πω , (A11)with the understanding that η is real-valued and the coefficients { c k } will obey c k = c − k , wherethe bar denotes complex conjugate.Next, we use ϕ ∗ ( ω ) = πω and Euler’s formula to writecos (cid:0) ϕ ∗ ( ω ′ ) − ϕ ∗ ( ω ) (cid:1) = (cid:16) e i π ( ω ′ − ω ) + e − i π ( ω ′ − ω ) (cid:17) . (A12)Inserting (A12) and (A11) into (A10) gives ∂ t η ( ω ) = − πη ( ω ) + K Z − ∑ k ∈ Z c k ( t ) e ik πω ′ (cid:16) e i π ( ω ′ − ω ) + e − i π ( ω ′ − ω ) (cid:17) d ω ′ . (A13)The only terms of the sum that do not vanish in the integral are those with k = ±
1. For the k = ± Z − c ( t ) e i πω ′ e − i π ( ω ′ − ω ) d ω ′ = c ( t ) e i πω (A14)and likewise for k = −
1. This shows that the coupling term acts on η diagonally in the Fourierbasis. Equating Fourier coefficients on each side of (A13), we obtain k = ± ∂ t c k ( t ) = (cid:18) − π + K (cid:19) c k ( t ) (A15) k = ± ∂ t c k ( t ) = − π c k ( t ) . (A16) REFERENCESREFERENCES L. M. Prolo, J. S. Takahashi, and E. D. Herzog, “Circadian rhythm generation and entrainmentin astrocytes,” Journal of Neuroscience , 404–408 (2005). I. Z. Kiss, Y. Zhai, and J. Hudson, “Emerging coherence in a population of chemical oscillators,”Science , 1676–1678 (2002). 21
I. Z. Kiss, C. G. Rusin, H. Kori, and J. L. Hudson, “Engineering complex dynamical structures:sequential patterns and desynchronization,” Science , 1886–1889 (2007). P. J. Uhlhaas and W. Singer, “Neural synchrony in brain disorders: relevance for cognitive dys-functions and pathophysiology,” Neuron , 155–168 (2006). E. Rodriguez, N. George, J.-P. Lachaux, J. Martinerie, B. Renault, and F. J. Varela, “Perception’sshadow: long-distance synchronization of human brain activity,” Nature , 430–433 (1999). F. D¨orfler, M. Chertkov, and F. Bullo, “Synchronization in complex oscillator networks andsmart grids,” Proceedings of the National Academy of Sciences , 2005–2010 (2013). A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,”Journal of theoretical biology , 15–42 (1967). A. T. Winfree,
The geometry of biological time , Vol. 12 (Springer Science & Business Media,2001). Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators,” in
Interna-tional symposium on mathematical problems in theoretical physics (Springer, 1975) pp. 420–422. Y. Kuramoto,
Chemical oscillations, waves, and turbulence (Springer Science & Business Me-dia, 2012). E. Brown, P. Holmes, and J. Moehlis, “Globally coupled oscillator networks,” in
Perspectivesand Problems in Nolinear Science (Springer, 2003) pp. 183–215. H. Sakaguchi, “Cooperative Phenomena in Coupled Oscillator Systems under External Fields,”Progress of Theoretical Physics , 39–46 (1988). T. Antonsen Jr, R. Faghih, M. Girvan, E. Ott, and J. Platig, “External periodic driving of largesystems of globally coupled phase oscillators,” Chaos: An Interdisciplinary Journal of NonlinearScience , 037112 (2008). E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems of globally coupledoscillators,” Chaos: An Interdisciplinary Journal of Nonlinear Science , 37113 (2008). E. Ott and T. M. Antonsen, “Long time evolution of phase oscillator systems,”Chaos: An Interdisciplinary Journal of Nonlinear Science , 023117 (2009), arXiv:0902.2773. E. Ott, B. R. Hunt, and T. M. Antonsen, “Comment on Long timeevolution of phase oscillator systems [Chaos 19 , 023117 (2009)],”Chaos: An Interdisciplinary Journal of Nonlinear Science , 025112 (2011), arXiv:1005.3319. L. M. Childs and S. H. Strogatz, “Stability diagram for the forced Kuramoto model,” Chaos: An22nterdisciplinary Journal of Nonlinear Science , 43128 (2008). A. V. Pimenova, D. S. Goldobin, M. Rosenblum, and A. Pikovsky, “Interplay of couplingand common noise at the transition to synchrony in oscillator populations,” Scientific reports , 38518 (2016). R. Adler, “A Study of Locking Phenomena in Oscillators,”Proceedings of the IRE , 351–357 (1946). B. Razavi, “A study of injection pulling and locking in oscillators,”Proceedings of the IEEE 2003 Custom Integrated Circuits Conference, 2003. , 1415–1424 (2003). A. C. Barnes, R. C. Roberts, N. C. Tien, C. A. Zorman, and P. X. L. Feng, “Siliconcarbide (SiC) membrane nanomechanical resonators with multiple vibrational modes,”2011 16th International Solid-State Sensors, Actuators and Microsystems Conference, TRANSDUCERS’11 , 2614–2617 (2011). J. Hunter and J. Milton, “Amplitude and Frequency Dependence of Spike Timing: Implicationsfor Dynamic Regulation,” Journal of Neurophysiology , 387–394 (2003). A. Zlotnik and J.-S. Li, “Optimal Subharmonic Entrainment of Weakly Forced Nonlinear Oscil-lators,” SIAM Journal on Applied Dynamical Systems , 1654–1693 (2014). D. Efimov, P. Sacr´e, and R. Sepulchre, “Controlling the Phase of an Oscillator: A Phase Re-sponse Curve Approach,” in
Joint 48th Conference on Decision and Control (2009) pp. 7692–7697. A. Zlotnik and J.-S. Li, “Optimal entrainment of neural oscillator ensembles,”Journal of Neural Engineering , 046015 (2012), arXiv:1202.5080. A. Zlotnik, Y. Chen, I. Z. Kiss, H.-A. Tanaka, and J.-S. Li, “Optimal Waveform for Fast Entrain-ment of Weakly Forced Nonlinear Oscillators,” Physical Review Letters , 024102 (2013). A. Zlotnik, R. Nagao, I. Z. Kiss, and J.-S. Li, “Phase-selective entrainment of nonlinear oscilla-tor ensembles,” Nature communications (2016), 10.1038/ncomms10788. S. H. Strogatz,
Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chem-istry and Engineering (Westview Press, 2001). F. C. Hoppensteadt and E. M. Izhikevich,
Weakly connected neural networks , Vol. 126 (SpringerScience & Business Media, 2012). H. Nakao, “Phase reduction approach to synchronisation of nonlinear oscillators,” ContemporaryPhysics , 188–214 (2016). B. Ermentrout, “Type I membranes, phase resetting curves, and synchrony,” Neural computation , 979–1001 (1996). 23 M. A. Schwemmer and T. J. Lewis, “The theory of weakly coupled oscillators,” in
Phase re-sponse curves in neuroscience (Springer, 2012) pp. 3–31. S. H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in popula-tions of coupled oscillators,” Physica D: Nonlinear Phenomena , 1–20 (2000). J. A. Acebr´on, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler,“The Kuramoto model: A simple paradigm for synchronization phenomena,”Reviews of Modern Physics , 137–185 (2005), arXiv:0306625 [cond-mat]. R. E. Mirollo and S. H. Strogatz, “Jump Bifurcation and Hystere-sis in an Infinite-Dimensional Dynamical System of Coupled Spins,”SIAM Journal on Applied Mathematics , 108–124 (1990). B. Ottino-L¨offler and S. H. Strogatz, “Kuramoto model with uniformly spaced frequen-cies: Finite-N asymptotics of the locking threshold,” Physical Review E , 062220 (2016),arXiv:1512.02321. D. Paz´o, “Thermodynamic limit of the first-order phase transition in the Kuramoto model,”Physical Review E , 046211 (2005), arXiv:0509020 [nlin]. E. J. Doedel, T. F. Fairgrieve, B. Sandstede, A. R. Champneys, Y. A. Kuznetsov, and X. Wang,“Auto-07p: Continuation and bifurcation software for ordinary differential equations,” Tech.Rep. (2007). H.-A. Tanaka, “Entrainment limit of weakly forced nonlinear oscillators,” in
Mathematical Ap-proaches to Biological Systems (Springer, 2015) pp. 77–93. P. Jiruska, M. De Curtis, J. G. Jefferys, C. A. Schevon, S. J. Schiff, and K. Schindler, “Syn-chronization and desynchronization in epilepsy: controversies and hypotheses,” The Journal ofphysiology , 787–797 (2013). L. B. Good, S. Sabesan, S. T. Marsh, K. Tsakalis, D. Treiman, and L. Iasemidis, “Control ofsynchronization of brain dynamics leads to control of epileptic seizures in rodents,” Internationaljournal of neural systems , 173–196 (2009). S. Ching, E. N. Brown, and M. A. Kramer, “Distributed control in a mean-field cortical networkmodel: implications for seizure suppression,” Physical Review E , 021920 (2012).24 IG. 3. (upper) The extent of the stable R > K c = / π as a function of N . The data showapproximately a power law scaling N − . for N between 3 and 100. (lower) The difference between thevalue of synchrony R at the saddle-node point and a numerically estimated critical value of R c = / N . The data show approximately a power law scaling N − . for N between 3 and 100. Circlesrepresent data measured from AUTO simulation, solid line is a power law fit.between 3 and 100. Circlesrepresent data measured from AUTO simulation, solid line is a power law fit.