Synchronization of Weakly Coupled Oscillators: Coupling, Delay and Topology
SSynchronization of Weakly Coupled Oscillators: Coupling, Delayand Topology
Enrique Mallada and Ao TangCornell University, Ithaca, NY 14853
Abstract
There are three key factors of a system of coupled oscillators that characterize the interactionamong them: coupling (how to affect) , delay (when to affect) and topology (whom to affect) . Foreach of them, the existing work has mainly focused on special cases. With new angles and tools,this paper makes progress in relaxing some assumptions of these factors. There are three mainresults in this paper. First, by using results from algebraic graph theory, a sufficient conditionis obtained which can be used to check equilibrium stability. This condition works for arbitrarytopology. It generalizes existing results and also leads to a sufficient condition on the couplingfunction with which the system is guaranteed to reach synchronization. Second, it is known thatidentical oscillators with sin() coupling functions are guaranteed to synchronize in phase on acomplete graph. Using our results, we demonstrate that for many cases certain structures insteadof exact shape of the coupling function such as symmetry and concavity are the keys for globalsynchronization. Finally, the effect of heterogenous delays is investigated. We develop a newframework by constructing a non-delayed phase model that approximates the original one in thecontinuum limit. We further derive how its stability properties depend on the delay distribution.In particular, we show that heterogeneity, i.e. wider delay distribution, can help reach in-phasesynchronization. a r X i v : . [ m a t h . O C ] M a r . INTRODUCTION The system of coupled oscillators has been widely studied in different disciplines rangingfrom biology [1–5] and chemistry [6, 7] to engineering [8, 9] and physics [10, 11]. The possiblebehavior of such a system can be complex. For example, the intrinsic symmetry of thenetwork can produce multiple limit cycles or equilibria with relatively fixed phases (phase-locked trajectories) [12], which in many cases can be stable [13]. Also, the heterogeneity inthe natural oscillation frequency can lead to incoherence [14] or even chaos [15].One particular interesting question is whether the coupled oscillators will synchronize inphase in the long run [16–20]. Besides its clear theoretical value, it also has rich applicationsin practice.In essence, there are three key factors of a system of coupled oscillators that characterizethe interaction among oscillators: coupling , delay and topology . For each of them, theexisting work has mainly focused on special cases as explained below. In this paper, furtherresearch will be discussed on each of these three factors: • Topology ( whom to affect , Section III B): Current results either restrict to completegraph or a ring topology for analytical tractability [19], study local stability of topol-ogy independent solutions over time varying graph [21–23], or introduce dynamic con-trollers to achieve synchronization for time-varying uniformly connected graphs [24,25]. We develop a graph based sufficient condition which can be used to check equi-librium stability for any fixed topology. It also leads to a family of coupling functionwith which the system is guaranteed to reach global phase consensus for arbitraryundirected connected graph using only physically meaningful state variables. • Coupling ( how to affect , Section III C): The classical Kuramoto model [14] assumes asin() coupling function. Our study hints that certain symmetry and convexity struc-tures are enough to guarantee global synchronization. • Delay ( when to affect , Section IV): Existing work generally assumes zero delay amongoscillators or require them to be bounded up to a constant fraction of the period [26].This is clearly not satisfactory especially if the oscillating frequencies are high. Wedevelop a new framework to study unbounded delays by constructing a non-delayedphase model that is equivalent to the original one. Using this result, we show thatwider delay distribution can help reach synchronization.In this paper we study weakly coupled oscillators, which can be either pulse-coupled orphase-coupled. Although most of the results are presented for phase-coupled oscillators, theycan be readily extended for pulse-coupled oscillators (see, e.g., [27, 28]). It is worth noticingthat results in Section III are independent of the strength of the coupling and therefore theweak coupling assumption is not necessary there. Preliminary versions of this work has beenpresented in [29] and [30].The paper is organized as follows. We describe pulse-coupled and phase-coupled oscil-lator models, as well as their common weak coupling approximation, in Section II. Usingsome facts from algebraic graph theory and potential dynamics in Section III A, we presentthe negative cut instability theorem in Section III B 1 to check whether an equilibrium isunstable. This then leads to Proposition 1 in Section III B 2 which identifies a class ofcoupling functions with which the system always synchronizes in phase. It is well known2hat the Kuramoto model produces global synchronization over a complete graph. In Sec-tion III C, we demonstrate that a large class of coupling functions, in which the Kuramotomodel is a special case, guarantee the instability of most of the limit cycles in a completegraph network. Section IV is devoted to the discussion of the effect of delay. An equivalentnon-delayed phase model is constructed whose coupling function is the convolution of theoriginal coupling function and the delay distribution. Using this approach, it is shown thatsometimes more heterogeneous delays among oscillators can help reach synchronization. Weconclude the paper in Section V.
II. COUPLED OSCILLATORS
We consider two different models of coupled oscillators studied in the literature. Thedifference between the models arises in the way that the oscillators interact between eachother, and their dynamics can be quite different. However, when the interactions are weak(weak coupling), both systems behave similarly and share the same approximation. Thisallows us to study them under a common framework.Each oscillator is represented by a phase θ i in the unit circle S which in the absence ofcoupling moves with constant speed ˙ θ i = ω. Here, S represents the unit circle, or equivalentlythe interval [0 , π ] with 0 and 2 π identified (0 ≡ π ), and ω = πT denotes the naturalfrequency of the oscillation. A. Pulse-coupled Oscillators
In this model the interaction between oscillators is performed by pulses. An oscillator j sends out a pulse whenever it crosses zero ( θ j = 0). When oscillator i receives a pulse, it willchange its position from θ i to θ i + εκ ij ( θ i ). The function κ ij represents how other oscillators’actions affect oscillator i and the scalar ε > δ satisfying δ ( t ) = 0 ∀ t (cid:54) = 0, δ (0) = + ∞ ,and (cid:82) δ ( s ) ds = 1. The coupled dynamics is represented by˙ θ i ( t ) = ω + εω (cid:88) j ∈N i κ ij ( θ i ( t )) δ ( θ j ( t − η ij )) , (1)where η ij > i and j ( η ij = η ji ), and N i isthe set of i ’s neighbors. The factor of ω in the sum is needed to keep the size of the jumpwithin εκ ij ( θ i ). This is because θ j ( t ) behaves like ωt when crosses zero and therefore thejump produced by δ ( θ j ( t )) is of size (cid:82) δ ( θ j ( t )) dt = ω − [28].The coupling function κ ij can be classified based on the qualitative effect it produces inthe absence of delay. After one period, if the net effect of the mutual jumps brings a pairof oscillators closer, we call it attractive coupling. If the oscillators are brought furtherapart, it is considered to be repulsive coupling. The former can be achieved for instance if κ ij ( θ ) ≤ θ ∈ [0 , π ) and κ ij ( θ ) ≥ θ ∈ [ π, π ). See Figure 1 for an illustration of anattractive coupling κ ij and its effect on the relative phases.This pulse-like interaction between oscillators was first introduced by Peskin [2] in 1975as a model of the pacemaker cells of the heart, although the canonic form did not appear inthe literature until 1999 [28]. In general, when the number of oscillators is large, there are3 θ / π κ st Jump θ / π κ nd Jump ij ijjij i
FIG. 1. Pulse-coupled oscillators with attractive coupling. several different limit cycles besides the in-phase synchronization and many of them can bestable [13].The question of whether this system can collectively achieve in-phase synchronization wasanswered for the complete graph case and zero delay by Mirollo and Strogatz in 1990 [20].They showed if κ ij ( θ ) is strictly increasing on (0 , π ) with a discontinuity in 0 (which resem-bles attractive coupling), then for almost every initial condition, the system can synchronizein phase in the long run.The two main assumptions of [20] are all to all comunication and zero delay. Whetherin-phase synchronization can be achieved for arbitrary graphs has been an open problem formore than twenty years. On the other hand, when delay among oscillators is introduced theanalysis becomes intractable. Even for the case of two oscillators, the number of possibilitiesto be considered is large [31, 32]. B. Phase-coupled Oscillators
In the model of phase-coupled oscillators, the interaction between neighboring oscillators i and j ∈ N i is modeled by change of the oscillating speeds. Although in general the speedchange can be a function of both phases ( θ i , θ j ), we concentrate on the case where the speedchange is a function of the phase differences f ij ( φ j ( t − η ij ) − φ i ( t )). Thus, since the net speedchange of oscillator i amounts to the sum of the effects of its neighbors, the full dynamics isdescribed by ˙ φ i ( t ) = ω + ε (cid:88) j ∈N i f ij ( φ j ( t − η ij ) − φ i ( t )) . (2)The function f ij is usually called coupling function, and as before η ij represents delay and N i is the set of neighbors of i . 4 θ / π f Attractive Coupling θ / π f Repulsive Coupling ijij
FIG. 2. Phase-coupled oscillators with attractive and repulsive coupling.
A similar definition for attractive and repulsive couplings can be done in this model. Wesay that the coupling function f ij is attractive if, without delays, the change in speeds bringsoscillators closer, and repulsive if they are brought apart. Figure 2 shows typical attractiveand repulsive coupling functions where arrows represent the speed change produced by theother oscillator; if the pointing direction is counter clockwise, the oscillator speeds up, andotherwise it slows down.When f ij = KN sin(), K > N i = N \{ i } , i.e.,complete graph topology) and no delay ( η ij = 0), see e.g. [19, 34], or to some regions of thestate space [26]. C. Weak Coupling Approximation
We now concentrate in the regime in which the coupling strength of both models is weak,i.e. 1 (cid:29) ε >
0. For pulse-coupled oscillators, this implies that the effect of the jumpsoriginated by each neighbor can be approximated by their average [27]. For phase-coupledoscillators, it implies that to the first order φ i ( t − η ij ) is well approximated by φ i ( t ) − ωη ij .The effect of these approximations allows us to completely capture the behavior of bothsystems using the following equation where we assume that every oscillator has the samenatural frequency ω and only keep track of the relative difference using˙ φ i = ε (cid:88) j ∈N i f ij ( φ j − φ i − ψ ij ) . (3)5or pulse-coupled oscillators, the coupling function is given by f ij ( θ ) = ω π κ ij ( − θ ) , (4)and the phase lag ψ ij = ωη ij represents the distance that oscillator i ’s phase can travelalong the unit circle during the delay time η ij . Equation (4) also shows that the attrac-tive/repulsive coupling classification of both models are in fact equivalent, since in order toproduce the same effect κ ij and f ij should be mirrored, as illustrated in Figure 1 and Figure2. Equation (3) captures the relative change of the phases and therefore any solution to (3)can be immediately translated to either (1) or (2) by adding ωt . For example, if φ ∗ is anequilibrium of (3), by adding ωt , we obtain a limit cycle in the previous models. Besidesthe delay interpretation for ψ ij , (3) is also known as a system of coupled oscillators with frustration , see e.g. [35].From now on we will concentrate on (3) with the understanding that any convergenceresult derived will be immediately true for the original models in the weak coupling limit.We are interested in the attracting properties of phase-locked invariant orbits within T N ,which can be represented by φ ( t ) = ω ∗ t N + φ ∗ , where N = (1 , . . . , T ∈ T N , and φ ∗ and ω ∗ are solutions to ω ∗ = ε (cid:88) j ∈N i f ij ( φ ∗ j − φ ∗ i − ψ ij ), ∀ i . (5)Whenever the system reaches one of these orbits, we say that it is synchronized or phase-locked. If furthermore, all the elements of φ ∗ are equal, we say the system is synchronized in-phase or that it is in-phase locked. It is easy to check that for a given equilibrium φ ∗ of(3), any solution of the form φ ∗ + λ N , with λ ∈ R , is also an equilibrium that identifies thesame limit cycle. Therefore, two equilibria φ , ∗ and φ , ∗ will be considered to be equivalent,if both identifies the same orbit, or equivalently, if both belongs to the same connected setof equilibria E φ ∗ := { φ ∈ T N | φ = φ ∗ + λ N , λ ∈ R } . (6) III. EFFECT OF TOPOLOGY AND COUPLING
In this section we concentrate on the class of coupling function f ij that are symmetric( f ij = f ji ∀ ij ), odd ( f ij ( − θ ) = − f ij ( θ )) and continuously differentiable. We also assumethat there is no delay within the network, i.e. ψ ij = 0 ∀ ij . Thus, (3) reduces to˙ φ i = ε (cid:88) j ∈N i f ij ( φ j − φ i ) . (7)In the rest of this section we progressively show how with some extra conditions on f ij wecan guarantee in-phase synchronization for arbitrary undirected graphs. Since we know thatthe network can have many other phase-locked trajectories besides the in-phase one, ourtarget is an almost global stability result [36], meaning that the set of initial conditionsthat does not eventually lock in-phase has zero measure. Latter we show how most of thephase-locked solution that appear on a complete graph are unstable under some generalconditions on the structure of the coupling function.6 . Preliminaries We now introduce some prerequisites used in our later analysis.
1. Algebraic Graph Theory
We start by reviewing basic definitions and properties from graph theory [37, 38] that areused in the paper. Let G be the connectivity graph that describes the coupling configuration.We V ( G ) and E ( G ) to denote the set of vertices ( i or j ) and undirected edges ( e ) of G . Anundirected graph G can be directed by giving a specific orientation σ to the elements in theset E ( G ). That is, for any given edge e ∈ E ( G ), we designate one of the vertices to be the head and the other to be the tail giving G σ .Although in the definitions that follow we need to give the graph G a given orientation σ ,the underlying connectivity graph of the system is assumed to be undirected . This is nota problem as the properties used in this paper are independent of a particular orientation σ and therefore are properties of the undirected graph G . Thus, to simplify notation we dropthe superscript σ from G σ with the understanding that G is now an induced directed graphwith some fixed, but arbitrarily chosen, orientation.We use P = ( V − , V + ) to denote a partition of the vertex set V ( G ) such that V ( G ) = V − ∪ V + and V − ∩ V + = ∅ . The cut C ( P ) associated with P , or equivalently C ( V − , V + ),is defined as C ( P ) := { ij ∈ E ( G ) | i ∈ V − , j ∈ V + , or vice versa. } . Each partition can beassociated with a vector column c P where c P ( e ) = 1 if e goes form V − to V + , c P ( e ) = − e goes form V + to V − and c P ( e ) = 0 if e stays within either set.There are several matrices associated with the oriented graph G that embed informationabout its topology. However, the one with most significance to this work is the orientedincidence matrix B ∈ R | V ( G ) |×| E ( G ) | where B ( i, e ) = 1 if i is the head of e , B ( i, e ) = − i is the tail of e and B ( i, e ) = 0 otherwise.
2. Potential Dynamics
We now describe how our assumptions on f ij not only simplifies considerably the dynam-ics, but also allows us to use the graph theory properties introduced in Section III A 1 togain a deeper understanding of (3).While f ij being continuously differentiable is standard in order to study local stability andsufficient to apply LaSalle’s invariance principle [39], the symmetry and odd assumptionshave a stronger effect on the dynamics.For example, under these assumptions the system (7) can be compactly rewritten in avector form as ˙ φ = − εBF ( B T φ ) (8)where B is the adjacency matrix defined in Section III A 1 and the map F : E ( G ) → E ( G ) is F ( y ) = ( f ij ( y ij )) ij ∈ E ( G ) . This new representation has several properties. First, from the properties of B one caneasily show that (5) can only hold with ω ∗ = 0 for arbitrary graphs [16] (since N ω ∗ = ω ∗ TN N = − ε TN BF ( B T φ ) = 0), which implies that every phase-locked solution is an7quilibrium of (7) and that every limit cycle of the original system (3) can be representedby some E ∗ φ on (7).However, the most interesting consequence of (8) comes from interpreting F ( y ) as thegradient of a potential function W ( y ) = (cid:88) ij ∈ E ( G ) (cid:90) y ij f ij ( s ) ds. Then, by defining V ( φ ) = ( W ◦ B T )( φ ) = W ( B T φ ), (8) becomes a gradient descent law for V ( φ ), i.e., ˙ φ = − εBF ( B T φ ) = − εB ∇ W ( B T φ ) = − ε ∇ V ( φ ) , where in the last step above we used the property ∇ ( W ◦ B T )( φ ) = B ∇ W ( B T φ ). Thismakes V ( φ ) a natural Lyapunov function candidate since˙ V ( φ ) = (cid:104)∇ V ( φ ) , ˙ φ (cid:105) = − ε |∇ V ( φ ) | = − ε (cid:12)(cid:12)(cid:12) ˙ φ (cid:12)(cid:12)(cid:12) ≤ . (9)Furthermore, since the trajectories of (8) are constrained into the N -dimensional torus T N , which is compact, V ( φ ) satisfies the hipotesis of LaSalle’s invariance principle (Theorem4.4 [39]), i.e. there is a compact positively invariant set, T N and a function V : T N → R thatdecreases along the trajectories φ ( t ). Therefore, for every initial condition, the trajectoryconverges to the largest invariant set M within { ˙ V ≡ } which is the equilibria set E = { φ ∈ T N | ˙ φ ≡ } = (cid:83) φ ∗ E φ ∗ . Remark 1.
The fact that symmetric and odd coupling induces potential dynamics is wellknow in the physics community [40]. However, it has been also rediscovered in the controlcommunity [17] for the specific case of sine coupling. Clearly, this is not enough to showalmost global stability, since it is possible to have other stable phase-locked equilibrium setsbesides the in-phase one. However, if we are able show that all the non-in-phase equilibriaare unstable, then almost global stability follows. That is the focus of the next section.
B. Negative Cut Instability Condition
We now present the main results of this section. Our technique can be viewed as ageneralization of [19]. By means of algebraic graph theory, we provide a better stabilityanalysis of the equilibria under a more general framework. We further use the new stabilityresults to characterize f ij that guarantees almost global stability.
1. Local Stability Analysis
In this section we develop the graph theory based tools to characterize the stability of eachequilibrium. We will show that given an equilibrium φ ∗ of the system (8), with connectivity8raph G and f ij as described in this section. If there exists a cut C ( P ) such that the sum (cid:88) ij ∈ C ( P ) f (cid:48) ij ( φ ∗ j − φ ∗ i ) < , (10)the equilibrium φ ∗ is unstable .Consider first an equilibrium point φ ∗ . Then, the first order approximation of (8) around φ ∗ is δ ˙ φ = − εB (cid:20) ∂∂y F ( B T φ ∗ ) (cid:21) B T δφ, were δφ = φ − φ ∗ is the incremental phase variable, and ∂∂y F ( B T φ ∗ ) ∈ R | E ( G ) |×| E ( G ) | is theJacobian of F ( y ) evaluated at B T φ ∗ , i.e., ∂∂y F ( B T φ ∗ ) = diag (cid:0) { f (cid:48) ij ( φ ∗ j − φ ∗ i ) } ij ∈ E ( G ) (cid:1) . Now let A = − εB (cid:104) ∂∂y F ( B T φ ∗ ) (cid:105) B T and consider the linear system δ ˙ φ = Aδφ.
Althoughit is possible to numerically calculate the eigenvalues of A given φ ∗ to study the stability,here we use the special structure of A to provide a sufficient condition for instability thathas nice graph theoretical interpretations.Since A is symmetric, it is straight forward to check that A has at least one positiveeigenvalue, i.e. φ ∗ is unstable, if and only if x T Ax >
0. Now, given any partition P =( V − , V + ), consider the associated vector c P , define x P such that x i = if i ∈ V + and x i = − if i ∈ V − . Then it follow from the definition of B that c P = B T x P which impliesthat − ε x TP Ax P = c TP (cid:20) ∂∂y F ( B T φ ∗ ) (cid:21) c P = (cid:88) ij ∈ C ( P ) f (cid:48) ij ( φ ∗ j − φ ∗ i ) . Therefore, when condition (10) holds, A = − εBDB T has at least one eigenvalue whosereal part is positive . Remark 2.
Equation (10) provides a sufficient condition for instability; it is not clearwhat happens when (10) does not hold. However, it gives a graph-theoretical interpretationthat can be used to provide stability results for general topologies. That is, if the minimumcut cost is negative, the equilibrium is unstable.
Remark 3.
Since the weights of the graph f (cid:48) ij ( φ ∗ j − φ ∗ i ) are functions of the phase difference, (10) holds for any equilibria of the form φ ∗ + λ N . Thus, the result holds for the whole set E φ ∗ defined in (6) . When (10) is specialized to P = ( { i } , V ( G ) \{ i } ) and f ij ( θ ) = sin( θ ), it reduces to theinstability condition in Lemma 2.3 of [19]; i.e., (cid:88) j ∈N i cos ( φ ∗ j − φ ∗ i ) < . (11)However, (10) has a broader applicability spectrum as the following example shows. Example 1.
Consider a six oscillators network as in Figure 3, where each node is linkedwith its four closest neighbors and f ij ( θ ) = sin( θ ) . Then, by symmetry, it is easy to verify hat φ ∗ = (cid:20) , π , π , π, π , π (cid:21) T (12) is an equilibrium of (7) . FIG. 3. The network of six oscillators (Example 4) t φ i π FIG. 4. Unstable equilibrium φ ∗ . Initial condition φ = φ ∗ + δφ We first study the stability of φ ∗ using (11) as in [19]. By substituting (12) in cos( φ ∗ j − φ ∗ i ) ∀ ij ∈ E ( G ) we find that the edge weights can only take two values: cos( φ ∗ j − φ ∗ i ) = (cid:40) cos( π ) = , if j = i ± π ) = − , if j = i ± Then, since any cut that isolates one node from the rest (like C = C ( { } , V ( G ) \{ } ) inFigure 3) will always have two edges of each type, their sum is zero . Therefore, (11) cannotbe used to determine stability. f we now use condition (10) instead, we are allowed to explore a wider variety ofcuts that can potentially have smaller costs. In fact, if instead of C we sum over C = C ( { , , } , { , , } ) , we obtain, (cid:88) ij ∈ C cos( φ ∗ j − φ ∗ i ) = − < , which implies that φ ∗ is unstable.Figure 4 verifies the equilibrium instability. By starting with an initial condition φ = φ ∗ + δφ close to the equilibrium φ ∗ , we can see how the system slowly starts to move awayfrom φ ∗ towards a stable equilibrium set.Furthermore, we can study the whole family of non-isolated equilibria given by φ ∗ = (cid:20) ε , π ε , π ε , π + ε , π ε , π ε (cid:21) T (13) where ε , ε , ε ∈ R , which due to Remark 3, we can reduce (13) to φ ∗ = (cid:20) , π λ , π λ , π, π λ , π λ (cid:21) T (14) with λ = ε − ε and λ = ε − ε . −4 −2 0 2 4−4−2024−6−5−4−3−2−1 λ λ FIG. 5. Minimum cut value C ∗ ( λ , λ ) showing that the equilibria (13) are unstable Instead of focusing on only one cut, here we compute the minimum cut value (10) overthe 31 possible cuts, i.e. C ∗ ( λ , λ ) := min P (cid:80) ij ∈ C ( P ) f (cid:48) ij ( φ j ( λ , λ ) ∗ − φ ∗ i ( λ , λ )) . Figure5 show the value of C ∗ ( λ , λ ) for λ i ∈ [ − π, π ] . Since C ∗ ( λ , λ ) is π -periodic on eachvariable and its value is negative for every λ , λ ∈ [ − π, π ] , the family of equilibria (14) (and consequently (13) ) is unstable. . Almost Global Stability Condition (10) also provides insight on which class of coupling functions can potentiallygive us almost global convergence to the in-phase equilibrium set E N . If it is possible tofind some f ij with f (cid:48) ij (0) >
0, such that for any non-in-phase equilibrium φ ∗ , there is a cut C with (cid:80) ij ∈ C f (cid:48) ij ( φ ∗ j − φ ∗ i ) <
0, then the in-phase equilibrium set will be almost globallystable [13]. The main difficulty is that for general f ij and arbitrary network G , it is not easyto locate every phase-locked equilibria and thus, it is not simple to know in what region ofthe domain of f ij the slope should be negative.We now concentrate on the one-parameter family of functions F b , with b ∈ (0 , π ), suchthat f ij ∈ F b whenever f ij is symmetric , odd , continuously differentiable and • f (cid:48) ij ( θ ; b ) > ∀ θ ∈ (0 , b ) ∪ (2 π − b, π ), and • f (cid:48) ij ( θ ; b ) < ∀ θ ∈ ( b, π − b ).See Figure 2 for an illustration with b = π . Also note that this definition implies that if f ij ( θ ; b ) ∈ F b , the coupling is attractive and f ij ( θ ; b ) > ∀ θ ∈ (0 , π ). This last property willbe used later. We also assume the graph G to be connected .In order to obtain almost global stability we need b to be small. However, since theequilibria position is not known a priori, it is not clear how small b should be or if there isany b > I be a compact connected subset of S and let l ( I ) be its length, e.g., if I = S then l ( I ) = 2 π . For any S ⊂ V ( G ) and φ ∈ T N , define d ( φ, S ) as the length of the smallestinterval I such that φ i ∈ I ∀ i ∈ S , i.e. d ( φ, S ) = l ( I ∗ ) = min I : φ i ∈ I , ∀ i ∈ S l ( I ) . Using this metric, together with the aid of Proposition 2.6 of [16] we can identify twovery insightful properties of the family F b whenever the graph G is connected. Claim 1. If φ ∗ is an equilibrium point of (8) with d ( φ ∗ , V ( G )) ≤ π , then either φ ∗ is anin-phase equilibrium, i.e. φ ∗ = λ N for λ ∈ R , or has a cut C with f (cid:48) ij ( φ ∗ j − φ ∗ i ) < ∀ ij ∈ C .Proof. Since d ( φ ∗ , V ( G )) ≤ π , all the phases are contained in a half circle and for theoscillator with smallest phase i , all the phase differences ( φ ∗ j − φ ∗ i ) ∈ [0 , π ]. However,since f ij ( · ; b ) ∈ F b implies f ij ( θ ; b ) ≥ ∀ θ ∈ [0 , π ] with equality only for θ ∈ { , π } ,˙ φ ∗ i = (cid:80) j ∈N i f ij ( φ ∗ j − φ ∗ i ) = 0 can only hold if φ ∗ j − φ ∗ i ∈ { , π } ∀ j ∈ N i . Now let V − = { i ∈ V ( G ) : d ( φ ∗ , { i, i } ) = 0 } and V + = V ( G ) \ V − . If V − = V ( G ), then φ ∗ is anin-phase equilibrium. Otherwise, ∀ ij ∈ C ( V − , V + ), f ij ( φ ∗ j − φ ∗ i ) = f ij ( π ) < b that guarantees the instabilityof the non-in-phase equilibria. Claim 2.
Consider f ij ( · ; b ) ∈ F b ∀ ij ∈ E ( G ) and arbitrary connected graph G . Then for any b ≤ πN − and non-in-phase equilibrium φ ∗ , there is a cut C with f (cid:48) ij ( φ ∗ j − φ ∗ i ; b ) < , ∀ ij ∈ C roof. Suppose there is a non-in-phase equilibrium φ ∗ for which no such cut C exists. Let V − = { i } and V +0 = V ( G ) \{ i } be a partition of V ( G ) for some arbitrary node i .Since such C does not exists, there exists some edge i j ∈ C ( V − , V +0 ), with j ∈ V +0 ,such that f (cid:48) i j ( φ ∗ j − φ ∗ i ; b ) ≥
0. Move j from one side to the other of the partition bydefining V − := V − ∪ { j } and V +1 := V +0 \{ j } . Now since f (cid:48) i j ( φ ∗ j − φ ∗ i ; b ) ≥
0, then d ( φ ∗ , V − ) ≤ b. In other words, both phases should be within a distance smaller than b .Now repeat the argument k times. At the k th iteration, given V − k − , V + k − , again we can findsome i k − ∈ V − k − , j k ∈ V + k − such that i k − j k ∈ C ( V − k − , V + k − ) and f (cid:48) i k − j k ( φ ∗ j k − φ ∗ i k − ; b ) ≥ d ( φ ∗ , { i k − , j k } ) ≤ b , d ( φ ∗ , V − k ) ≤ b + d ( φ ∗ , V − k − ) . Thus by solving the recursion we get: d ( φ ∗ , V − k ) ≤ kb. After N − V − N − = V ( G ) and d ( φ ∗ , V ( G )) ≤ ( N − b . Therefore,since b ≤ πN − , we obtain d ( φ ∗ , V ( G )) ≤ ( N − πN − π. Then, by Claim 1 φ ∗ is either an in-phase equilibrium or there is a cut C with f (cid:48) ij ( φ ∗ j − φ ∗ i ) < ∀ ij ∈ C . Either case gives a contradiction to assuming that φ ∗ is a non-in-phase equilibriumand C does not exists. Therefore, for any non-in-phase φ ∗ and b ≤ πN − , we can always finda cut C with f ij ( φ ∗ j − φ ∗ i ; b ) < ∀ ij ∈ C .Claim 2 allows us to use our cut condition (10) on every non-in-phase equilibrium. Thus,since (8) is a potential dynamics (c.f. Section III A 2), from every initial condition the systemconverges to the set of equilibria E . But when b ≤ πN − the only stable equilibrium set inside E is the in-phase set E N . Thus, E N set is globally asymptotically stable. We summarizedthis result in the following Proposition. Proposition 1 (Almost global stability) . Consider f ij ( θ ; b ) ∈ F b and an arbitrary con-nected graph G . Then, if b ≤ πN − , the in-phase equilibrium set E N is almost globallyasymptotically stable . This result provides a sufficient condition for almost global asymptotic stability to the in-phase equilibrium set E N . Although found independently, the same condition was proposedfor a specific piecewise linear f ij in [41]. Here we extend [41] in many aspects. For example,instead of assuming equal coupling for every edge, our condition describes a large family ofcoupling functions F b where each f ij can be taken independently from F b . Also, in [41] theconstruction of f ij ( θ ) assumes a discontinuity on the derivative at θ = b . This can pose aproblem if the equilibrium φ ∗ happens to have phase differences φ ∗ j − φ ∗ i = b . Here we donot have such problem as f ij is continuously differentiable.The condition b ≤ πN − implies that, when N is large, f ij should be decreasing in most ofit domain. Using (4) this implies that κ ij should be increasing within the region ( b, π − b ),which is similar to the condition on [20] and equivalent when b →
0. Thus, Proposition 1confirms the conjecture of [20] by extending their result to arbitrary topologies and a morerealistic continuous κ ij for the system (1) in the weak coupling limit.13 . Complete Graph Topology with a Class of Coupling Functions In this subsection we investigate how conservative the value of b found in Section III B 2is for the complete graph topology. We are motivated by the results of [19] where it isshown that f ( θ ) = sin( θ ) ( b = π ) with complete graph topology ensures almost globalsynchronization.Since for general f it is not easy to characterize all the possible equilibria of the system,we study the stability of the equilibria that appear due to the equivalence of (8) with respectto the action group S N × T , where S N is the group of permutations of the N coordinatesand T = [0 , π ) represents the group action of phase shift of all the coordinates, i.e. theaction of δ ∈ T is φ i (cid:55)→ φ i + δ ∀ i . We refer the readers to [12] and [16] for a detailed studyof the effect of this property.These equilibria are characterized by the isotropy subgroups Γ of S N × T that keep themfixed, i.e., γφ ∗ = φ ∗ ∀ γ ∈ Γ. In [12] it was shown that this isotropy subgroup takes the formof ( S k × S k × · · · × S k lB − ) m (cid:111) Z m where k i and m are positive integers such that ( k + k + · · · + k l B − ) m = N , S j is thepermutation subgroup of S N of j -many coordinates and Z m is the cyclic group with action φ i (cid:55)→ φ i + πm . The semiproduct (cid:111) represents the fact that Z m does not commute with theother subgroups.In other words, each equilibria with isotropy ( S k × S k × · · · × S k lB − ) m (cid:111) Z m is conformedby l B shifted constellations C l ( l ∈ { , , . . . l B − } ) of m evenly distributed blocks, with k l oscillators per block. We use δ l to denote the phase shift between constellation C and C l .See Figure 6 for examples these types of equilibria. d d k k k m=4 km=8 FIG. 6. Equilibria with isotropy ( S k × S k × S k ) (cid:111) Z (left) and ( S k ) (cid:111) Z (right) Here we will show that under mild assumptions on f and for b = π most of the equilibriafound with these characteristics are unstable. We first study all the equilibria with m even.In this case there is a special property that can be exploited.14hat is, when f ∈ F π such that f is even around π , we have g m ( δ ) := m − (cid:88) j =0 f ( 2 πm j + δ ) (15)= m/ − (cid:88) j =0 f ( 2 πm j + δ ) + f ( π + 2 πm j + δ )= m/ − (cid:88) j =0 f ( 2 πm j + δ ) + f (( 3 π πm j + δ ) − π m/ − (cid:88) j =0 f ( 2 πm j + δ ) + f ( − ( 2 πm j + δ ))= m/ − (cid:88) j =0 f ( 2 πm j + δ ) − f ( 2 πm j + δ ) = 0where the third step comes from f being even around π/ π -periodic, and the fourthfrom f being odd.Having g m ( δ ) = 0 is the key to prove the instability of every equilibria with even m .It essentially states that the aggregate effect of one constellation C l on any oscillator j ∈ V ( G ) \ C l is zero when m is even, and therefore any perturbation that maintains C l has nulleffect on j . This is shown in the next proposition. FIG. 7. Cut of Proposition 2, the red block represents one possible set V Proposition 2 (Instability for even m ) . Given an equilibrium φ ∗ with isotropy ( S k × S k ×· · · × S k lB ) m (cid:111) Z m and f ∈ F π even around π . Then, if m is even, φ ∗ is unstable.Proof. We will show the instability of φ ∗ by finding a cut of the network satisfying (10).Let V ⊂ V ( G ) be the set of nodes within one of the blocks of the constellation C andconsider the partition induced by V , i.e. P = ( V , V ( G ) \ V ). Due to the structure of φ ∗ ,1510) becomes (cid:88) ij ∈ C ( P ) f (cid:48) ( φ ∗ j − φ ∗ i ) = − k f (cid:48) (0) + l B (cid:88) l =1 k l g (cid:48) m ( δ l ) , where g (cid:48) m ( δ ) is the derivative of g m and δ l is the phase shift between the C and C l . Finally,since by assumptions g m ( δ ) ≡ ∀ δ then it follows that g (cid:48) m ( δ ) ≡ (cid:88) ij ∈ C ( P ) f (cid:48) ij ( φ ∗ j − φ ∗ i ) = − k f (cid:48) (0) < . Therefore, by (10), φ ∗ is unstable.The natural question that arises is whether similar results can be obtained for m odd.The main difficulty in this case is that g m ( δ ) = 0 does not hold since we no longer evaluate f at points with phase difference equal to π such that they cancel each other. Therefore, anextra monotonicity condition needs to be added in order to partially answer this question.These conditions and their effects are summarized in the following claims. Claim 3 (Monotonicity) . Given f ∈ F π such that f is strictly concave for θ ∈ [0 , π ] , then f (cid:48) ( θ ) − f (cid:48) ( θ − φ ) < , ≤ θ − φ < θ ≤ π (16) f (cid:48) ( θ ) − f (cid:48) ( θ + φ ) < , − π ≤ θ < θ + φ ≤ Proof.
The proof is a direct consequence of the strict concavity of f . Since f ( θ ) is strictlyconcave then basic convex analysis shows that f (cid:48) ( θ ) is strictly decreasing within [0 , π ]. There-fore, the inequality (16) follows directly from the fact that θ ∈ [0 , π ], θ − φ ∈ [0 , π ] and θ − φ < θ . To show (17) it is enough to notice that since f is odd ( f ∈ F π ), f is strictlyconvex in [ π, π ]. The rest of the proof is analogous to (16). Claim 4 ( f (cid:48) Concavity) . Given f ∈ F π such that f (cid:48) is strictly concave for θ ∈ [ − π , π ] .Then for all m ≥ , f (cid:48) ( πm ) ≥ f (cid:48) (0) .Proof. Since f (cid:48) ( θ ) is concave for θ ∈ [ − π, π ] then it follows f (cid:48) ( πm ) = f (cid:48) ( λ m − λ m ) π > λ m f (cid:48) (0) + (1 − λ m ) f (cid:48) ( π > λ m f (cid:48) (0)where λ m = m − m . Thus, for m ≥ λ m ≥ and f (cid:48) ( πm ) > f (cid:48) (0)as desired.Now we show the instability of any equilibria with isotropy ( S k × S k × · · · × S k lB ) m (cid:111) Z m for m odd and greater or equal to 7. Proposition 3 (Instability for m ≥ . Suppose f ∈ F π with f concave in [0 , π ] and f (cid:48) concave in [ − π , π ] , then for all m = 2 k + 1 with k ≥ the equilibria φ ∗ s with isotropy ( S k × S k × · · · × S k lB ) m (cid:111) Z m are unstable. IG. 8. Cut used in Proposition 3. The dots in red represent all the oscillators of some maximalset S with d ( φ ∗ , S ) < πm The proof of Proposition 3 also uses our cut condition to show instability, but with adifferent cut induced by the partition P = ( S, V ( G ) \ S ) of V ( G ) where S is set to the amaximal subset of V ( G ) such that d ( φ, S ) < πm , see Figure 8 for an illustration of P . Noticethat any of these partitions will include all the oscillators of two consecutive blocks of everyconstellation. The details of the proofs are rather technical and are relegated to AppendixA. IV. EFFECT OF DELAY
Once delay is introduced to the system of coupled oscillators, the problem becomes fun-damentally harder. For example, for pulse-coupled oscillators, the reception of a pulse nolonger gives accurate information about the relative phase difference ∆ φ ij = φ j − φ i betweenthe two interacting oscillators. Before, at the exact moment when i received a pulse from j , φ j was zero and the phase difference was estimated locally by i as ∆ φ ij = − φ i . How-ever, now when i receives the pulse, the difference becomes ∆ φ ij = − φ i − ψ ij . Therefore,the delay propagation acts as an error introduced to the phase difference measurement andunless some information is known about this error, it is not possible to predict the behavior.Moreover, as we will see later, slight changes in the distribution can produce nonintuitivebehaviors.Even though it may not be satisfactory for some applications, many existing works chooseto ignore delay. (see for e.g., [18–20]). That is mainly for analytical tractability. On theother hand, when delay is included [26] the studies concentrate on finding bounds on delaythat maintain stability.In this section we study how delay can change the stability in a network of weakly coupledoscillators. A new framework to study these systems with delay will be set up by constructingan equivalent non-delayed system that has the same behavior as the original one in thecontinuum limit. We then further use this result to show that large heterogeneous delaycan help reach synchronization, which is a bit counterintuitive and significantly generalizesprevious related studies [28, 42, 43]. We will assume complete graph to simplify notationand exposition although the results can be extended for a boarder class of densely connected17etworks.The contribution of this section is two fold. First, it improves the understanding of theeffect of delays in networks of coupled oscillators. And second, it opens new possibilities ofusing delay based mechanisms to increase the region of attraction of the in-phase equilib-rium set. We shall build on existing arguments such as mean field approximation [33] andLyapunov stability theory [17, 19] while looking at the problem from a different perspective. A. Mean Field Approximation
Consider the case where the coupling between oscillators is all to all and identical ( N i = N \{ i } , ∀ i ∈ N and f ij = f ∀ i, j ). And assume the phase lags ψ ij are randomly andindependently chosen from the same distribution with probability density g ( ψ ). By letting N → + ∞ and ε → εN =: ¯ ε a constant, (3) becomes v ( φ, t ) := ω + ¯ ε (cid:90) π − π (cid:90) + ∞ f ( σ − φ − ψ ) g ( ψ ) ρ ( σ, t ) dψdσ, (18)where ρ ( φ, t ) is a time-variant normalized phase distribution that keeps track of the fractionof oscillators with phase φ at time t , and v ( φ, t ) is the velocity field that expresses the netforce that the whole population applies to a given oscillator with phase φ at time t . Sincethe number of oscillators is preserved at any time, the evolution of ρ ( φ, t ) is governed by thecontinuity equation ∂ρ∂t + ∂∂φ ( ρv ) = 0 (19)with the boundary conditions ρ (0 , t ) ≡ ρ (2 π, t ). Equations (18)-(19) are not analyticallysolvable in general. Here we propose a new perspective that is inspired by the followingobservation.Consider the non-delayed system of the form˙ φ i = ω + ε (cid:88) j ∈N i H ( φ j − φ i ) , (20)where H ( θ ) = f ∗ g ( θ ) = (cid:90) + ∞ f ( θ − ψ ) g ( ψ ) dψ (21)is the convolution between f and g .By the same reasoning of (18) it is easy to see that the limiting velocity field of (20) is v H ( φ, t ) = ω + ¯ ε (cid:90) π H ( σ − φ ) ρ ( σ, t ) dσ = ω + ¯ ε (cid:90) π (cid:18)(cid:90) + ∞ f (( σ − φ ) − ψ ) g ( ψ ) dψ (cid:19) ρ ( σ, t ) dσ = ω + ¯ ε (cid:90) π (cid:90) + ∞ f ( σ − φ − ψ ) g ( ψ ) ρ ( σ, t ) dψdσ = v ( φ, t ) 18here in the first and the third steps we used (21) and (18) respectively. Therefore, (3) and(20) have the same continuum limit.Remark 4. Although (20) is quite different from (3) , both systems behave exactly the samein the continuum limit. Therefore, as N grows, (20) starts to become a good approximationof (3) and therefore can be analyzed to understand the behavior of (3) . θ / π θ / π θ / π g ∗ f f ∗ g FIG. 9. Effect of delay in coupling shape
Figure 9 shows how, the underlying delay (in this case the delay distribution) determineswhat type of coupling (attractive or repulsive) produces synchronization. The original func-tion f produces repulsive coupling, whereas the corresponding H is attractive. In fact, as wewill soon see, the distribution of delay not only can qualitatively affect the type of couplingbut also can change the stability of certain phase-locked limit cycles.We now study two example to illustrate how this new approximation can provide signif-icant information about performance and stability of the original system. We also providenumerical simulations to verify our predictions. B. Kuramoto Oscillators
We start by studying an example in the literature [44] to demonstrate how we can use theprevious equivalent non-delayed formulation to provide a better understanding of systemsof coupled oscillators with delay. When f ( θ ) = K sin( θ ), H ( θ ) can be easily calculated: H ( θ ) = (cid:90) + ∞ K sin( θ − ψ ) g ( ψ ) dψ = K (cid:90) + ∞ (cid:61) [ e i ( θ − ψ ) g ( ψ )] dψ = K (cid:61) [ e iθ (cid:90) + ∞ e − iψ g ( ψ ) dψ ]= K (cid:61) (cid:2) e iθ Ce − iξ (cid:3) = KC sin( θ − ξ )where (cid:61) is the imaginary part of a complex number, i.e. (cid:61) [ a + ib ] = b . The values of C > ξ are calculated using the identity Ce iξ = (cid:90) + ∞ e iψ g ( ψ ) dψ. This complex number, usually called “order parameter”, provides a measure of how thephase-lags are distributed within the unit circle. It can also be interpreted as the center ofmass of the lags ψ ij ’s when they are thought of as points ( e iψ ij ) within the unit circle S .19hus, when C ≈
1, the ψ ij ’s are mostly concentrated around ξ . When C ≈
0, the delay isdistributed such that (cid:80) ij e iψ ij ≈ φ i = ω + εKC (cid:88) j ∈N i sin( φ j − φ i − ξ ) . (22)Here we see how the distribution of g ( ψ ) has a direct effect on the dynamics. For ex-ample, when the delays are heterogeneous enough such that C ≈
0, the coupling term dis-appears and therefore makes synchronization impossible. A complete study of the systemunder the context of superconducting Josephson arrays was performed [44] for the completegraph topology. There the authors characterized the condition for in-phase synchroniza-tion in terms of K and Ce iξ . More precisely, when KCe iξ is on the right half of the plane( KC cos( ξ ) > KCe iξ is on theleft half of the plane ( KC cos( ξ ) < N (cid:80) Nl =1 e iφ l , becomes zero. ψ g Lags Distributions −1 −0.5 0 0.5 1−1−0.500.51
Order Parameters
FIG. 10. Delay distributions and their order parameter Ce iξ We now provide simulation results to illustrate how (22) becomes a good approximationof the original system when N is large enough. We simulate the original repulsive ( K < N = 5 , ,
50. Figure 11 shows that when N issmall, the phases’ order parameter of the original system (in red/blue) draw a trajectorywhich is completely different with respect to its approximation (in green). However, as N grows, in both cases the trajectories become closer and closer. Since K <
0, the trajectoryof the system with wider distribution ( C cos ξ <
0) drives the order parameter towards theboundary of the circle, i.e., heterogeneous delay leads to homogeneous phase . C. Effect of Heterogeneity
We now explain a more subtle effect that heterogeneity can produce. Consider the systemin (20) where H odd and continuously differentiable. Then, from Section III, all of the20 IG. 11. Repulsive sine coupling with heterogeneous delays oscillators eventually end up running at the same speed ω with fixed phase difference suchthat the sum (cid:80) i ∈N i H ( φ j − φ i ) cancels ∀ i . Moreover, we can apply (10) to assess the stabilityof these orbits. Therefore, if we can find a cut C of the network such that (cid:80) ij ∈ C H (cid:48) ( φ ∗ j − φ ∗ i ) < , the phase-locked solution will be unstable.Although this condition is for non-delayed phase-coupled oscillators, the result of thissection allows us to translate it for systems with delay. Since H is the convolution of thecoupling function f and the delay distribution function g , we can obtain H (cid:48) ( φ ∗ j − φ ∗ i ) < f (cid:48) ( φ ∗ j − φ ∗ i ) >
0. This usually occurs when the convolution widens the regionwith a negative slope of H . See Figure 9 for an illustration of this phenomenon.Figures 12 and 13 show two simulation setups of 45 oscillators pulse-coupled all to all.The initial state is close to a phase locked configuration formed of three equidistant clustersof 15 oscillators each. The shape of the coupling function f and the phase lags distributionsare shown in part a. We used (4) to implement the corresponding pulse-coupled system ( ?? ).While f is maintained unchanged between both simulations, the distribution g does change.Thus, the corresponding H = f ∗ g changes as it can be seen in part b; the blue, red, andgreen dots correspond to the speed change induced in an oscillator within the blue clusterby oscillators of each cluster. Since all clusters have the same number of oscillators, the neteffect is zero. In part c the time evolution of oscillators’ phases relative to the phase of ablue cluster oscillator are shown. Although the initial conditions are exactly the same, thewider delay distribution on Figure 13 produces negative slope on the red and green pointsof part b, which destabilizes the clusters and drives oscillators toward in-phase synchrony.Finally, we simulate the same scenario as in Figures 12 and 13 but now changing N andthe standard deviation, i.e. the delay distribution width. Figure 14 shows the computationof the synchronization probability vs. standard deviation. The dashed line denotes theminimum value that destabilizes the equivalent system. As N grows, the distribution shapebecomes closer to a step, which is the expected shape in the limit. It is quite surprising thatas soon as the equilibrium is within the region of H with negative slope, the equilibriumbecomes unstable as the theory predicts. 21 θ / π θ / π g a b f f ∗ g FIG. 12. Pulse-coupled oscillators with delay: Stable equilibrium
V. CONCLUSION
This paper analyzes the dynamics of identical weakly coupled oscillators while relaxingseveral classical assumptions on coupling, delay and topology. Our results provide globalsynchronization guarantee for a wide range of scenarios. There are many directions thatcan be taken to further this study. For example, for different topologies, to guarantee globalin-phase synchronization, how does the requirement on coupling functions change? Anotherspecific question is to complete the proof in Section III C for the cases when m = 1 , , Acknowledgments:
The authors thank Steven H. Strogatz of Cornell for useful discussions.
Appendix A: Proof of Proposition 3
As in Proposition 2 we will use our cut condition to show the instability of φ ∗ . Thus,we define a partition P = ( S, V ( G ) \ S ) of V ( G ) by taking S to be a maximal subset of V ( G ) such that d ( φ, S ) < πm , see Figure 8 for an illustration of P . Notice that any of thesepartitions will include all the oscillators of two consecutive blocks of every constellation.Instead of evaluating the total sum of the weights in the cut we will show that the sum22 θ / π θ / π g ba f f ∗ g FIG. 13. Pulse-coupled oscillators with delay: Unstable equilibrium of edge weights of the links connecting the nodes of one constellation in S with the nodesof a possibly different constellation in V ( G ) \ S is negative. In other words, we will focus onshowing (cid:88) ij ∈K l l f (cid:48) ( φ ∗ j − φ ∗ i ) < K l l = { ij : i ∈ C l ∩ S, j ∈ C l ∩ V ( G ) \ S } .Given any subset of integers J , we define g Jm ( δ ) = g m ( δ ) − (cid:88) j ∈ J f ( 2 πm j + δ ) . IG. 14. Pulse-coupled oscillators with delay: Synchronization probability
Then, we can rewrite (A1) as (cid:88) ij ∈K l l f (cid:48) ( φ ∗ j − φ ∗ i ) ==( g { , } m ) (cid:48) ( δ l l ) + ( g {− , } m ) (cid:48) ( δ l l )=2 g (cid:48) m ( δ l l ) − f (cid:48) ( δ l l + 2 πm ) − f (cid:48) ( δ l l ) − f (cid:48) ( δ l l − πm ) (A2)where δ l l ∈ [0 , πm ] is the phase shift between the two constellations. Then, if we can showthat for all δ ∈ [0 , πm ] (A2) is less than zero then for any values of l and l we will have(A1) satisfied.Since f is odd and even around π , f (cid:48) is even and odd around π and g (cid:48) m ( δ ) can be rewritten24s g (cid:48) m ( δ ) = f (cid:48) ( δ )+ (cid:88) ≤| j |≤(cid:98) k (cid:99) (cid:26) f (cid:48) ( δ + 2 πm j ) − f (cid:48) ( δ − sgn ( j ) πm + 2 πm j ) (cid:27) − (cid:104) f (cid:48) ( δ + πm k ) + f (cid:48) ( δ − πm k ) (cid:105) [ k odd] where [ k odd] is the indicator function of the event [ k odd], the sum is over all the integers j with 1 ≤ | j | ≤ (cid:98) k (cid:99) and k = m − The last term only appears when k is odd and in fact it is easy to show that it is alwaysnegative as the following calculation shows: − f (cid:48) ( δ + πm k ) − f (cid:48) ( δ − πm k ) == − f (cid:48) ( πm k + δ ) − f (cid:48) ( πm k − δ )= − f (cid:48) ( π − π m + δ ) − f (cid:48) ( π − π m − δ )= f (cid:48) ( π − δ + π m ) − f (cid:48) ( π − δ − π m )= f (cid:48) ( θ ) − f (cid:48) ( θ − φ ) < f (cid:48) being even, in step two we used k = m − and in stepthree we use f (cid:48) being odd around π . The last step comes from substituting θ = π − δ + π m , φ = πm and apply Claim 3, since for m ≥ ≤ θ − φ < θ ≤ π .Then it remains the show that the terms of the form f (cid:48) ( δ + πm j ) − f (cid:48) ( δ − sgn ( j ) πm + πm j )are negative for all j s.t. 1 ≤ | j | ≤ (cid:98) k (cid:99) . This is indeed true when j is positive since for all δ ∈ [0 , πm ] we get 0 ≤ δ − πm + 2 πm j < δ + 2 πm j ≤ π, for 1 ≤ j ≤ (cid:98) k (cid:99) and thus we can apply again Claim 3.When j is negative there is one exception in which Claim 3 cannot be used since − π ≤ δ + 2 πm j < δ + 2 πm j + πm ≤ , ∀ δ ∈ [0 , πm ]only holds for −(cid:98) k (cid:99) ≤ j ≤ −
2. Thus the term corresponding to j = − j = ± g (cid:48) m is strictly upper boundedfor all δ ∈ [0 , πm ] by g (cid:48) m ( δ ) Mathematical aspects of heart physiology . New York, NY, USA: Courant Instituteof Mathematical Sciences, New York University, 1975.[3] P. Achermann and H. Kunz, “Modeling circadian rhythm generation in the suprachiasmaticnucleus with locally coupled self-sustained oscillators: Phase shifts and phase response curves,” Journal of Biological Rhythms , vol. 14, no. 6, pp. 460 – 468, 1999.[4] J. Garcia-Ojalvo, M. B. Elowitz, S. H. Strogatz, and C. S. Peskin, “Modeling a syntheticmulticellular clock: Repressilators coupled by quorum sensing,” Proceedings of the NationalAcademy of Sciences of the United States of America , vol. 101, no. 30, pp. 10955–10960, 2004.[5] S. Yamaguchi, H. Isejima, T. Matsuo, R. Okura, K. Yagita, M. Kobayashi, and H. Okamura, Synchronization of cellular clocks in the suprachiasmatic nucleus,” Science , vol. 32, pp. 1408–1412, Nov. 2003.[6] I. Z. Kiss, Y. Zhai, and J. L. Hudson, “Emerging coherence in a population of chemicaloscillators,” Science , vol. 296, pp. 1676–1678, Nov. 2002.[7] R. York and R. Compton, “Quasi-optical power combining using mutually synchronized os-cillator arrays,” IEEE Transactions on Microwave Theory and Techniques , vol. 39, pp. 1000–1009, Jun. 1991.[8] Y.-W. Hong and A. Scaglione, “A scalable synchronization protocol for large scale sensornetworks and its applications,” IEEE Journal on Selected Areas in Communications , vol. 23,pp. 1085–1099, May. 2005.[9] G. Werner-Allen, G. Tewari, A. Patel, M. Welsh, and R. Nagpal, “Firefly-inspired sensor net-work synchronicity with realistic radio effects,” in SenSys: Proceedings of the 3rd InternationalConference on Embedded Networked Sensor Systems , (New York, NY, USA), pp. 142–153,ACM, 2005.[10] S. A. Marvel and S. H. Strogatz, “Invariant submanifold for series arrays of josephson junc-tions,” Chaos , vol. 19, p. 013132, Mar. 2009.[11] P. C. Bressloff and S. Coombes, “Travelling waves in chains of pulse-coupled integrate-and-fireoscillators with distributed delays,” Phys. D , vol. 130, no. 3-4, pp. 232–254, 1999.[12] P. Ashwin and J. W. Swift, “The dynamics of n weakly coupled identical oscillators,” J.Nonlinear Sci , vol. 2, no. 1, pp. 69–108, 1992.[13] G. B. Ermentrout, “Stable periodic solutions to discrete and continuum arrays of weaklycoupled nonlinear oscillators,” SIAM J. Appl. Math. , vol. 52, no. 6, pp. 1665–1687, 1992.[14] Y. Kuramoto, “International symposium on mathematical problems in theoretical physics,”in Lecture notes in Physics , vol. 39, p. 420, Springer, 1975.[15] O. V. Popovych, Y. L. Maistrenko, and P. A. Tass, “Phase chaos in coupled oscillators,” Phys.Rev. E , vol. 71, p. 065201, Jun. 2005.[16] E. Brown, P. Holmes, and J. Moehlis, “Globally coupled oscillator networks,” in Perspec-tives and Problems in Nonlinear Science: A Celebratory Volume in Honor of Larry Sirovich ,pp. 183–215, Springer, 2003.[17] A. Jadbabaie, N. Motee, and M. Barahona, “On the stability of the kuramoto model ofcoupled nonlinear oscillators,” in Proceedings of the American Control Conference. , vol. 5,pp. 4296–4301, June 30, 2004.[18] D. Lucarelli and I.-J. Wang, “Decentralized synchronization protocols with nearest neighborcommunication,” in Proceedings of the 2nd International Conference on Embedded NetworkedSensor Systems , 2004.[19] P. Monz´on and F. Paganini, “Global considerations on the kuramoto model of sinusoidallycoupled oscillators,” in Proceedings of the 44th IEEE Conference on Decision and Control,and European Control Conference , (Sevilla, Spain), pp. 3923–3928, Dec. 2005.[20] R. E. Mirollo and S. H. Strogatz, “Synchronization of pulse-coupled biological oscillators.,” SIAM J. Appl. Math , vol. 50, pp. 1645–1662, 1990.[21] L. Moreau, “Stability of continuous-time distributed consensus algorithms,” in Decision andControl, 2004. CDC. 43rd IEEE Conference on , vol. 4, pp. 3998 – 4003 Vol.4, dec. 2004.[22] W. Ren and R. Beard, “Consensus of information under dynamically changing interactiontopologies,” in American Control Conference, 2004. Proceedings of the 2004 , vol. 6, pp. 4939–4944 vol.6, 30 2004-july 2 2004.[23] G.-B. Stan and R. Sepulchre, “Dissipativity characterization of a class of oscillators and net- orks of oscillators,” in Decision and Control, 2003. Proceedings. 42nd IEEE Conference on ,vol. 4, pp. 4169 – 4173 vol.4, dec. 2003.[24] L. Scardovi, A. Sarlette, and R. Sepulchre, “Synchronization and balancing on the n-torus,” Systems & Control Letters , vol. 56, no. 5, pp. 335 – 341, 2007.[25] R. Sepulchre, D. A. Paley, and N. E. Leonard, “Stabilization of planar collective motion withlimited communication,” Automatic Control, IEEE Transactions on , vol. 53, no. 3, pp. 706–719, 2008.[26] A. Papachristodoulou and A. Jadbabaie, “Synchonization in oscillator networks with heteroge-neous delays, switching topologies and nonlinear dynamics,” in Proceedings of the 45th IEEEConference on Decision and Control , pp. 4307–4312, Dec. 2006.[27] E. M. Izhikevich, “Phase models with explicit time delays,” Phys. Rev. E , vol. 58, pp. 905–908,Jul. 1998.[28] E. M. Izhikevich, “Weakly pulse-coupled oscillators, fm interactions, synchronization, andoscillatory associative memory,” IEEE Transactions on Neural Networks , vol. 10, pp. 508–526, May. 1999.[29] E. Mallada and A. Tang, “Synchronization of phase-coupled oscillators with arbitrary topol-ogy,” in Proceedings of American Control Conference , 2010.[30] E. Mallada and A. Tang, “Weakly pulse-coupled oscillators: Heterogeneous delays lead tohomogeneous phase,” in IEEE Conference on Decision and Control , 2010.[31] U. Ernst, K. Pawelzik, and T. Geisel, “Synchronization induced by temporal delays in pulse-coupled oscillators,” Phys. Rev. Lett. , vol. 74, pp. 1570–1573, Feb. 1995.[32] U. Ernst, K. Pawelzik, and T. Geisel, “Delay-induced multistable synchronization of biologicaloscillators,” Phys. Rev. E , vol. 57, pp. 2150–2162, Feb. 1998.[33] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence . Berin Heidelberg New YorkTokyo: Springer-Verlag, 1984.[34] R. Sepulchre, D. A. Paley, and N. E. Leonard, “Stabilization of planar collective motion: All-to-all communication,” Automatic Control, IEEE Transactions on , vol. 52, no. 5, pp. 811–824,2007.[35] H. Daido, “Quasientrainment and slow relaxation in a population of oscillators with randomand frustrated interactions,” Phys. Rev. Lett. , vol. 68, pp. 1073–1076, Feb 1992.[36] A. Rantzer, “A dual to lyapunov’s stability theorem,” Systems and Control Letters , vol. 42,pp. 161–168, 2001.[37] B. Bollobas, Modern Graph Theory . New York: Springer, 1998.[38] C. Godsil and G. Royle, Algebraic Graph Theory . New York: Springer, 2001.[39] H. K. Khalil, Nonlinear systems; 3rd ed. Prentice-Hall, 1996.[40] F. C. Hoppensteadt and E. M. Izhikevich, “Weakly connected neural networks,” 1997.[41] A. Sarlette, Geometry and symmetries in coordination control . PhD thesis, University of Li`ege,Belgium, January 2009.[42] C. van Vreeswijk, L. Abbott, and G. B. Ermentrout, “When inhibition not excitation syn-chronizes neural firing,” Journal of Computational Neuroscience , vol. 1, no. 4, pp. 313–321,1994.[43] W. Gerstner, “Rapid phase locking in systems of pulse-coupled oscillators with delays,” Phys.Rev. Lett. , vol. 76, pp. 1755–1758, Mar 1996.[44] S. Watanabe and S. H. Strogatz, “Constants of motion for superconducting josephson arrays,” Phys. D , vol. 74, no. 3-4, pp. 197–253, 1994., vol. 74, no. 3-4, pp. 197–253, 1994.