Stability of synchronized oscillations in networks of phase-oscillators

We derive simple conditions for the stability or instability of the synchronized oscillation of a class of networks of coupled phase-oscillators, which includes many of the systems used in neural modelling.


introduction
The general equations for a system of n coupled phase-oscillators are given by (1) θ ′ i = g i (θ 1 , θ 2 , ..., θ n ), 1 ≤ i ≤ n, where θ i denotes the phase of the i-th oscillator, and the functions g i encapsulate both the internal dynamics of the i-th oscillator and its coupling to the other oscillators. Since the θ i 's denote phases, we assume that the functions g i are 2π periodic with respect to all variables.
The study of the dynamics of networks of phase-oscillators is a formidable mathematical problem, and it is therefore important to isolate special structures of such systems which facilitate the study of at least some aspects of their dynamics. Of course, in order for the work to be relevant to modelling, these structures must be natural in the sense that they reflect some modelling assumptions. In this work we identify a special structure for systems of the form (1) which makes the study of synchronization much more tractable. Moreover this structure is general enough to encompass many of the systems which have been used in neural modelling in the 'frequency domain' approach [8,9].
A synchronized oscillation of the system (1), is a solution for which θ i (t) = θ(t) for all 1 ≤ i ≤ n, with (2) θ(t + T ) = θ(t) + 2π ∀t ∈ R We shall say that the system is synchronized if there exists a stable synchronized oscillation (see section 3 for the definition of stability). The study of synchronization is an important part in the quest to understand the dynamics of networks of coupled oscillators [11,13]. We derive simple conditions for the stability/instability of synchronization for a class of systems. As we shall point out, these conditions generalize several stability conditions derived in previous works for more restricted classes of systems.
We first present a special case of the structure that we identify in this work, before turning to the general case. Assume that the functions g i have the special form: Where h(α) is a 2π-periodic function describing the internal dynamics of the oscillators, f (α, β) is a function which is 2π-periodic in both variables and describes the coupling of pairs of oscillators, and c ij are constants describing the strength of the influence of oscillator j on oscillator i. Let us note the modelling assumptions implied by such a form for the functions g i : (i) The influence of all the oscillators on oscillator i is the sum of terms each one of which represents the influence of one of the other oscillators. This is known in the literature as the assumption of conventional synaptic connections.
(ii) The type of coupling among pairs of oscillators, represented by the function f , is identical for each pair of oscillators, and only the strength c ij of the coupling is allowed to vary with i and j (including the possibility that c ij = 0, so that there is no influence of oscillator j on oscillator i). A more general model, still satisfying the assumption of conventional synaptic connections, would have f depending on the indices i, j, but such a model falls outside the scope of this study and does not admit the pleasing analytical results proved here.
(iii) The internal dynamics of the oscillators, represented by the function h, are identical. Models of the form (3) are already of great interest from a modelling perspective [1,2,4,5,9,10], but our results in fact apply to more general models in which the functions g i have the form: The function S(φ, y), 2π-periodic in the first variable, can be thought of as a cut-off function limiting the input to each oscillator (see [7,8] for such examples), so its dependence on y (which we do not restrict) might be sigmoidal or hump-shaped . We note that this structure relaxes the 'conventional synaptic connections' assumption (i). Systems of the form (3) are special cases of systems of the form (4) in which S(φ, y) = h(φ) + y.

Existence of a synchronized oscillation
We first deal with the rather trivial matter of existence of a synchronized oscillation for systems with the structure (4).
Condition (5) guarantees that each of the equations (1) reduce to (8) upon substituting θ 1 = θ 2 = ... = θ n = θ, while condition (6) implies that each solution of (8) satisfies (2) for an appropriate T > 0. It is easy to see that (6) is a necessary condition for the existence of a solution of (8) satisfying (2), hence for the existence of a synchronized oscillation. We note that the period T of the synchronized oscillation is given by the formula We make some remarks about the condition (5): (i) Condition (5) says that the sum of strengths of the couplings of each oscillator to all the others does not depend on the oscillator, and it is easy to see that if this condition is violated a synchronized oscillation will usually not exist, because such a function would have to satisfy two independent scalar differential equations. (ii) Condition (5) can be considered as an idealization -in a 'general' network we do not expect the sum of strengths of the connections to each oscillator to be exactly constant. As a consequence, in a 'general' network we do not expect a synchronized oscillation to exist at all. However, there are two very important classes of systems for which condition will naturally hold: symmetric networks and regular networks, as we explain below. Moreover, we expect many natural 'large' networks to satisfy (5) in an approximate way, leading to approximate synchronization -see below.
(iii) A symmetry of a network is a permutation τ ∈ S n with the property that c ij = c τ (i)τ (j) for all 1 ≤ i, j ≤ n. The set of symmetries of the network forms a subgroup G of S n . It is easy to check that if the subgroup G is transitive, that is if for any 1 ≤ i, j ≤ n, there exists a τ ∈ G such that τ (i) = j, then (5) holds. Thus if we call a network symmetric whenever its symmetry group is transitive, we have that any symmetric network satisfies (5). A simple case is a cyclically-symmetric ring of oscillators, in which the symmetry group is generated by the permutation τ (i) = i + 1 (mod n). For more on symmetry in networks of oscillators see [6]. (iv) We note the trivial fact that a linear combination of two matrices satisfying (5) also satisfies (5), so that our condition holds, for example, if we have two independent sets of connections, each of them with a different transitive symmetry group.
(v) Consider the case of a network in which the connections have either strength s or strength 0 (i.e., no connection). Such a network can be represented by a directed graph in which the oscillators are the nodes and an arrow from node j to node i indicates that c ij = s. If this graph is k-regular in the sense that each node has k incoming arrows, then condition (5) is satisfied. This includes 'random' networks generated by asigning k oscillators to each oscillator. (vi) If we relax the assumption of regularity and generate a random network by a (perhaps more natural) process of setting c ij = s with probability p and c ij = 0 with probability 1 − p, we can still expect that when the number n of oscillators is large, the network will satisfy (5) in an approximate way. Indeed by the central limit theorem the sums in (5) will be approximately pns, with an error of O( 1 √ n ). Thus such a network can be considered as a small perturbation of a network satisfying (5) with c = pns. Under such circumstances, we expect the perturbation to change the synchronized oscillation into a nearby entrained oscillation, that is a solution satifying θ i (t + T ) = θ i (t) + 2π for all 1 ≤ i ≤ n. Moreover we expect the functions θ i to be close to one another, though not identical. We thus have the phenomenon of 'near-synchronization', and we also expect that the stability of the 'near-synchronized' oscillation will be the same as that of the synchronized solution of the system satisfying (5) off which we perturb. These heuristic remarks need to be backed by a careful analysis, but we leave this to future work.
We formulate the corollary of theorem 1 applied to systems of the form (3):

stability of the synchronized oscillation
We now turn to the more interesting issue of stability of the synchronized oscillation.
We recall some standard definitions. Let d denote the euclidean distance in R n . A solution θ(t) = (θ 1 (t), ..., θ n (t)) of (1) is said to be Liapunov stable if for any ǫ > 0 there is a δ > 0 such that whenever t 0 ∈ R,θ 0 ∈ R n and d(θ 0 , θ(t 0 )) < δ, we have d(θ(t), θ(t)) < ǫ for all t > 0, whereθ(t) denotes the solution of (1) with initial conditionθ (0) We note that we are dealing here with local stability, so that the fact that the synchronized oscillation is stable does not imply that synchronization will be achieved from any initial state, but only that there is some open set of initial conditions that will lead to synchronization. The study of global synchronization properties of systems of phase-oscillators is a very intriguing problem.
In general the stability of a synchronized oscillation is to be determined by examining the time-periodic linear system obtained by linearizing (1) around the synchronized solution: which in general can only be done by numerical computation. The crucial point is that the special structure (4) entails a special structure for the linearized system (12), which makes it possible to proceed much further with the analytic investigation than in the general case. Differentiating (4) we have Thus substituting θ 1 = θ 2 = ... = θ n = θ (14) we obtain where I is the identity matrix, C is the connection matrix C ij = c ij , and the scalar functions a(θ), b(θ) are given by: Let us assume for the moment that the connection matrix C is diagonalizable: where λ k , 1 ≤ k ≤ n are the eigenvalues of C. We note that (5) implies that (1, 1, ..., 1) is an eigenvector of C with eigenvalue c, so that we may assume Let F (t) denote the n × n matrix-valued function which is the fundamental solution of (12), so that F (0) = I and (20) We recall from the stability theory for periodic solutions of autonomous systems that the synchronized oscillation is stable if all the eigenvalues of F (T ) (T denoting the period of θ), known as the characteristic multipliers, are in the interior of the unit disk in the complex plane, except for one eigenvalue which is 1 (related to time-translation invariance). If at least one of the eigenvalues is in the exterior of the unit disk, the synchronized oscillation is unstable (see, e.g., [12], Ch. V, theorem 8.4).
Thus µ 1 = 1 is the eigenvalue of F (T ) corresponding to time-translations, and the sufficient condition for stability of the synchronized solution is that all the other eigenvalues, given by (23), are inside the unit disc. This is equivalent to the condition that the numbers σ i have negative real parts for 2 ≤ i ≤ n.
Recalling the definition (17) of b(θ), we finally obtain the expression: We note now that although (27) was derived under the assumption that C is diagonalizable, it continues to be valid for any matrix C satisfying (5), because the set of diagonalizable matrices is dense and both sides of (27) are continuous in C. We thus have: (5) and (6) are satisfied. Let cf (θ, θ)) .

(i) If
then the synchronized oscillation of (7) is stable.
We formulate the stability result for the special case of conventional synaptic coupling: Theorem 4. Assume (5) and (10) are satisfied. Let .
The calculations involved in determining the stability of the synchronized oscillation using theorems 3 or 4 split into two parts: calculation of the integral χ, and study of the spectrum of the connection matrix C. The following reformulation of theorem 3 highlights the independence of these two aspects: Theorem 5. Assume (5), (6) hold, and let χ be defined by (28).
(i) If then the synchronized solution of any system of the form (7) with χ = 0 is unstable.
then the synchronized solution of (7) is stable if χ < 0 and unstable if χ > 0.
then the synchronized solution of (7) is stable if χ > 0 and unstable if χ < 0.
We note that χ depends on the connection matrix only through the 'overall strength' c.
We assume (5) holds. In order to satisfy (10) we must have ω−c sin(θ) cos(θ) > 0 for all θ, which is easily seen to be equivalent to the requirement Thus if (5) and 35) hold, there exists a synchronized oscillation of (34) computing χ according to (30) we obtain and we see that χ > 0 whenever (35) holds. Thus from theorem 4 we obtain that the synchronized solution is stable iff the real parts of the eigenvalues of C (other than c itself) are all smaller than c, or, in other words, when all the eigenvalues of C − cI have negative real parts.

Connection matrices and the stability of synchronized oscillations
The spectrum of the connection matrix C is, of course, totally independent of the nonlinearities h, f, S. This fact allows us, for some connection matrices, to conclude that the synchronized solution is unstable, independently of the internal dynamics of the oscillators and of the coupling functions, as shown by part (i) of theorem 5.
In this section we consider the role of the connection matrix through some examples, using some well-known results from linear algebra. 4.1. Two-oscillator networks. In the case n = 2 it is easy to see that if we assume no self-coupling, (c 11 = c 22 = 0) the general system (7) is in fact no more general than the particular case (11), and in this case we may also absorb h into the f , so we are dealing with the system: and we have λ 1 = c = 1, λ 2 = −1 Thus (33) holds and from theorem 5 we obtain: Theorem 6. For a system of two coupled phase-oscillators of the form (36) for which f (θ, θ) > 0 ∀θ ∈ R, let: The synchronized solution is stable if χ > 0 and unstable if χ < 0.

4.2.
Three-oscillator networks. We assume that there is no self-coupling, so that the connection matrix can be written in the form By computing the characteristic polynomial of C−cI, we find that the product of its two eigenvalues (besides the eigenvalue 0) is where (38) ∆ = 3c 2 − (c 12 + c 23 + c 31 )c + (c 12 c 23 + c 12 c 31 + c 23 c 31 ), and that their sum is Assume ∆ < 0. Then λ 2 ,λ 3 must be real (because otherwise they are conjugate which would imply that the product (λ 2 −c)(λ 3 −c) is non-negative), and (λ 2 − c), (λ 3 − c) are of opposite signs, hence part (i) of theorem 5 implies that the synchronized oscillation is unstable. Now assume ∆ > 0. If λ 2 ,λ 3 are non-real, hence conjugate, then we have 4.3. Larger networks. In principle, for larger networks, with theorem 5, one can obtain algebraic conditions involving the elements of the connection matrix and χ for determining the stability/instability of the synchronized solution, by applying the Routh-Hurwitz conditions to the characteristic polynomial of C − cI. Obviously these conditions become more complicated as the size of the network gets larger, and we shall not proceed along these lines but rather consider some specific classes of connection matrices for which simple results are possible.

4.4.
Non-negative connection matrices. A very important case is that in which all the entries of the connection matrix are non-negative. In this case we can apply the Perron-Frobenius theorem to obtain stability results. In fact to apply the Perron-Frobenius theorem we need to assume also that the matrix C is irreducible, which means that there does not exist a partition of the indices {1, ..., n} into two disjoint sets A, B such that i ∈ A, j ∈ B implies c ij = 0. In other words, we assume that we cannot partition the oscillators into disjoint sets A, B so that oscillators in B have no influence on oscillators in A (we note that a simple sufficient condition for C to be irreducible is that c ij = 0 for all i = j). The Perron-Frobenius theorem states that if C is nonnegative and irreducible then it has a unique eigenvalue λ which is positive and has an eigenvector whose components are positive. Moreover λ is a simple eigenvalue and all other eigenvalues of C are less than or equal in modulus to λ.
Since in our case we already know of the eigenvalue λ 1 = c which has an eigenvector (1, 1, ..., 1) with positive components, the uniqueness part of the Perron-Frobenius theorem implies that all other eigenvalues satisfy λ i = c and |λ i | ≤ c. These two facts imply that Re(λ i ) < c, so that (33) holds. We thus conclude: Theorem 8. If the connection matrix C satisfies (5) and is non-negative and irreducible, and (6) holds, then the synchronized solution of (7) is stable if χ > 0 and unstable if χ < 0, where χ is defined by (28).
Of course an analogous theorem is true for the case of non-positive matrices, where the roles of the conditions χ > 0, χ < 0 are reversed.
A special case of theorem 8 was proved in [4] (section 2.5), with χ given in a somewhat less explicit form. There it was assumed that the system is of the form (11) with h a constant and that f is of the special product form f (α, β) = f 1 (α)f 2 (β). The unneccesary assumption that C is a symmetric matrix was made.
We note that theorem 8 includes the case of global uniform coupling given by c ij = s for all i, j (or for all i = j). This case has been treated, for the special case of (11) in which the coupling f (α, β) depends only β, in [5], and for the special case in which the coupling is of the form f (α, β) = f 0 (β − α) + f 1 (α)f 2 (β) and h is constant, in [2]. The stability results for synchronized oscillations in these works are corollaries of theorem 8.

4.5.
Cyclically symmetric systems. We now consider the case of a ring of oscillators with cyclic symmetry, which means that the effect of the i-th oscillator on the j-th oscillator is the same as that of the i + k-th oscillator j+k-th oscillator (where i+k, j+k are considered modulo n). In mathematical terms this means that the connection matrix C is circulant [3], that is, one having the form with d ∈ R n some vector. As we noted before, a circulant matrix automatically satisfies (5).
We have the following formula for the eigenvalues of a circulant matrix, where ρ = e −2πi/n : The formula (41), together with the identity 1 − cos(x) = 2 sin 2 ( x 2 ), give us the following stability result for cyclically symmetric matrices: Theorem 9. Assume (6) is satisfied, and assume that the network is cyclically symmetric, its connection matrix given by (39). Let χ be defined by (28), and define, for 1 ≤ k ≤ n: γ k = n m=1 d m sin 2 π (k − 1)(m − 1) n .
(i) If some of the γ k 's are postitive and some are negative, the synchronized oscillation of (7) is unstable.
(ii) If χγ k > 0 for all 1 ≤ k ≤ n then the synchronized oscillation is stable.
(iii) If χγ k < 0 for all 1 ≤ k ≤ n then the synchronized oscillation is unstable.
4.6. Two populations of interacting oscillators. Let us consider 2n oscillators, subdivided into two groups, each of size n, with the pairs of oscillators within each group coupled with strength p, and any pair of oscillators belonging to different groups are coupled with strength q.
For example in the case n = 3 the connection matrix is 0 p p q q q p 0 p q q q p p 0 q q q q q q 0 p p q q q p 0 p q q q p p 0 Such a matrix satisfies (5) with c = (n−1)p+nq. Let us denote by e k ∈ R 2n (1 ≤ k ≤ 2n) vector with 1 in the k-th place and 0 elesewhere. It is is easy to verify that the eigenvalues of C are c with eigenvector (1, 1, 1..., 1), (n−1)p−nq with eigenvector (1, 1, ..1, −1, −1, ... − 1), and −p with the eigenvectors v i , w i (2 ≤ i ≤ n), given by v i = e 1 − e i , 2 ≤ i ≤ n, w i = e n+1 − e n+i , 2 ≤ i ≤ n.