PRACTICAL SYNCHRONIZATION OF GENERALIZED KURAMOTO SYSTEMS WITH AN INTRINSIC DYNAMICS

. We study the practical synchronization of the Kuramoto dynamics of units distributed over networks. The unit dynamics on the nodes of the network are governed by the interplay between their own intrinsic dynamics and Kuramoto coupling dynamics. We present two suﬃcient conditions for practical synchronization under homogeneous and heterogeneous forcing. For practical synchronization estimates, we employ the conﬁguration diameter as a Lyapunov functional, and derive a Gronwall-type diﬀerential inequality for this value.


1.
Introduction. Collective synchronized behavior in coupled oscillators often appears in many complex biological systems, such as groups of fireflies, neurons, and cardiac pacemaker cells [1,5,31,32]. The synchronization phenomenon arising from a pair of pendulum clocks hanging on the same bar was first reported in the literature by Huygens in 1665. However, its mathematical treatment was first investigated by two pioneers, Winfree [31] and Kuramoto [19], about forty years ago. Since then, Kuramoto's first-order model and its extension have been extensively studied in various disciplines [1,8,13,23,24,27]. The Kuramoto model has simple intrinsic dynamics governed by the natural frequency, so that the uncoupled Kuramoto oscillator's phase has linear dynamics on the unit circle. A natural questions is "If the intrinsic dynamics are heterogeneous and rather complicated, can we still expect some kind of synchrony among oscillators?" 788 SEUNG-YEAL HA, SE EUN NOH AND JINYEONG PARK Such situations can be easily found in several examples, e.g., the daily cycling of light and darkness affecting human sleep rhythms [30]. External fields can also model the external current applied to a neuron, so as to describe the collective properties of excitable systems with planar symmetry. For other physical devices, such as Josephson junctions, a periodic external force can model an oscillating current across the junctions.
The main purpose of this paper is to study the dynamics of Kuramoto units with heterogeneous intrinsic dynamics located on a symmetric and connected network G = (V, E, Ψ) where V = {u 1 , · · · , u N } and E ⊂ V × V are vertex and edge sets, respectively and Ψ = (ψ ij ) is an N × N matrix whose element ψ ij the capacities of the edge connecting from u j to u i . For a given network or graph G, we assume that Kuramoto oscillators are located at the nodes of the network V , and that they interact symmetrically through the interacting channels registered by the connection topology E and Ψ. Let ζ i = ζ i (t) be a quantifiable description of the state of unit at node i. In the absence of coupling, we assume that ζ i ∈ R is governed by its own dynamics:ζ (1) where we assume that the forcing function F i : R × R + → R is a C 1 -function. We now consider the case where the above decoupled dynamics (1) interact with each other through the connecting edges of the network. In this case, we assume that the dynamics of ζ i are governed by the following coupled non-autonomous system: where K is a positive coupling strength. The static interaction matrix Ψ = (ψ ij ) is assumed to be symmetric and connected in the sense that (i) ψ ij = ψ ji ≥ 0, 1 ≤ i, j ≤ N, (ii) For any (i, j) ∈ V × V , there is a shortest path from i to j, say i = k 0 → k 1 → · · · → k dij = j, (k l , k l+1 ) ∈ E, l = 0, 1, . . . , d ij − 1. ( Here, d ij denotes the distance between nodes i and j, i.e., the length of the shortest path from node i to node j. Note that for all-to-all coupling with ψ ij = 1 N and F i = Ω i , system (2) reduces to the Kuramoto model: Thus, we can view system (2) as a generalized Kuramoto model. Kuramoto-type models with external forcing terms has been addressed in the literature, e.g., [2,26,28], and can be used to model the sleep-wake cycle. The sleep-wake cycle and circadian rhythms are phase-locked to each other in the 24-hour period of outside world. Many biological experiments have shown that, in isolation from a 24-hour periodic environment such as the light-dark cycle, the various circadian rhythms, e.g., feeding, body temperature, and neuroendocrine variables, as well as the pattern of sleep and wakefulness, were maintained. However, a certain internal desynchronization phenomenon occurs, i.e., separate rhythmic variables oscillating with different periods. Many mathematical models have been deveploped to explain these phenomena, notably the Kuramoto model with an external periodic force In particular, Sakaguchi [28] showed numerically that the forced entrainment is not always achieved, and analytic studies of this feature have also been reported [2,9,26]. Our main interest is not restricted to periodic forcing terms. Instead, we consider general forcing terms without assuming the periodicity of F i in the first argument. In the presence of heterogeneous forcing terms in (2), in general the difference ζ i −ζ j andζ i −ζ j do not converge to a constant value asymptotically. Thus, we cannot use the concept of complete synchronization of employed in [11,13]. For this reason, we need to adopt a different notion of synchrony, namely "practical synchronization" roughly saying that ζ-differences as values in R can be made as small by taking large coupling strength K (See Definition 2.1). This clearly generalizes the concept of the complete synchronization in [11,13]. For a motivation of practical synchronization and its relevance to biological and engineering applications, we refer to Remark 1.
The main novelty of this paper is to provide two frameworks to ensure the practical synchronization of system (2) in terms of the forcing F i , coupling strength K, and initial configurations ζ 0 under homogeneous and heterogeneous forcing. (See Section 2.1) The rest of this paper is organized as follows. In Section 2, we present several a priori estimates for later sections and greenprovides the summary of framework and main results. In Section 3, we provide a complete synchronization for homogeneous forcing, and in Section 4, we study the practical synchronization of heterogeneous forcing. Finally, Section 5 summarizes our main results.

Preliminaries.
In this section, we study the concept of practical synchronization, and provide several basic estimates to be used in later sections. We first set the phase-diameter and energy as follows. For a configuration ζ = (ζ 1 , · · · , ζ N ) ∈ R N , we set We define the concept of complete synchronization and practical synchronization for Kuramoto oscillators as follows. Remark 1. 1. In previous literature, the practical synchronization appeared in chaotic systems [4,10,12,20,21,22] and in the first-order linear consensus model [17]. In Definition 2.1, we closely follows the stronger notion of practical synchronization from [17] saying that ζ-differences can be made arbitrary small by suitable controls. Note that in our system (2), the magnitude of control terms which are the sinusoidal couplings is dominated by the coupling strength K. In [4,10,12,20,21,22,29], the weaker concepts of practical synchronization comparing the Definition 2.1 were used to denote the uniform boundedness of phase differences in time but no restriction on the bound of the phase differences according to the coupling K. The numerical experiment of [29] shows that large coupling strength is necessary to obtain sufficiently small bound of phase diameter. 2. For the study of synchronization phenomena of Kuramoto oscillators with intrinsic dynamics, complete synchronization cannot occur in general. To see this we consider the following two-oscillator system: Numerical simulation result in Figure 1 clearly shows that the differences ζ 1 − ζ 2 andζ 1 −ζ 2 do not converges to zero so we cannot obtain complete synchronization in phase and frequency [9,11,13]. However, we can observe that the differences become smaller as coupling strength K is increased, in other words, this system is practically synchronized in the sense of Definition 2.1 Lemma 2.2. [14] Suppose that the network (V (G), E(G)) is connected and the set {ζ i } has zero mean: Then, we have where the constant L * is given by Here, E c denotes the complement of the edge set E in V × V and |E c | denotes its cardinality.
Lemma 2.3. Suppose that the phase configuration ζ = (ζ 1 , · · · , ζ N ) ∈ R N satisfies Then, we have In particular, if we have the additional zero sum condition where the constant C ∞ is defined by relation (5).   Proof. We use the following elementary inequality

SEUNG-YEAL HA, SE EUN NOH AND JINYEONG PARK
to obtain Here, the third inequality uses Lemma 2.2.
Then, the variance V(ζ) satisfies where D(F) is the diameter of the family of forcing terms Proof. We first note that the average ζ c and its perturbationζ i := ζ i − ζ c satisfẏ We multiply the second equation of (6) by 2ζ i , sum with respect to i, and divide by N to find Here, we used the symmetry of the network, ψ ij = ψ ji , and the trick i ↔ j. We now consider the two terms separately.

793
• (Estimate of I 11 ): There exists ζ * ij on the segment between ζ i and ζ j such that where we used the fact that ζ i − ζ j =ζ i −ζ j . Note that where we used (8) and • (Estimate of I 12 ): It follows from Lemma 2.3 that we have We combine (8)-(9) together with (7) to obtain the desired estimate.

Remark 2.
Note that for an all-to-all coupling with ψ ij = 1 N , we have N C ∞ = sin D 0 . Lemma 2.5. For T ∈ (0, ∞], let η = η(t) and ζ = ζ(t) be the corresponding solutions to systems (1) and (2) with the same initial data ζ 0 and satisfying the following a priori conditions: 1. The total energy for the decoupled system is bounded: 2. The phase-diameter is confined to the half-circle region: there exists D 0 ∈ (0, π) such that Then, the coupled solution

SEUNG-YEAL HA, SE EUN NOH AND JINYEONG PARK
Proof. We will show that the energy of the coupled system is smaller than that of the decoupled system. Then, by the framework of (10), the energy for the decoupled system is bounded, and we can derive the desired result. For the boundedness of E(ζ), we multiply (2) by 2ζ i and sum with respect to i to obtain where we used Lemma 2.3 to find We now consider the solution to the decoupled system η i with the same initial data: To obtain the time derivative of the energy functional of decoupled system E(η), we estimate d dt Then, by the comparison principle of ODEs, we have Then, the a priori condition (10) yields 2.1. The descriptions of main results. In this subsection, we summarize our two main results on synchronization phenomena for homogeneous and heterogeneous forcing. We first consider the the case where all forcing functions are the same, i.e., F i = F for all 1 ≤ i ≤ N . Therefore, the system (2) becomes For this homogeneous forcing, we have the following complete synchronization.
Theorem 2.6. Suppose that the forcing function F and the coupling strength K satisfy Then for any solution ζ = ζ(t) to the system (11), we have where C 0 is a positive constant defined by On the other hand, we consider heterogeneous forcing terms: There exists a pair i = j such that F i = F j .
We adopt the following framework A on the family F = {F 1 , · · · , F N } of forcing and network structure Ψ and the coupling strength K.
• (A3) The initial configuration satisfies the following boundedness condition: Theorem 2.7. Suppose that the frameworks A holds, then practical synchronization is achieved: 3. Complete synchronization: Homogeneous forcing. In this section, we consider the special case where all forcing terms are equal, so that each member of some given uncoupled system has identical dynamics: In this case, We recall that system (2) becomes The average phase ζ c and fluctuationsζ i = ζ i − ζ c satisfẏ Lemma 3.1. Suppose that the forcing function F and the coupling strength K satisfy Proof of claim. We split the proof into two parts. In Step A, we show that the set T is nonempty, and in Step B, we show that T * = ∞ using the differential inequality obtained in Lemma 2.4.
We are now ready to prove the theorem 2.6.
The proof of Theorem 2.6. We multiply the second equation in (14) by 2ζ i , sum with respect to i, and divide by N to obtain d dt where ζ * ij (t) is a point in the interval connecting ζ i (t) and ζ j (t). Remark 3. Note that for the homogeneous forcing case if the identical forcing function F is non-increasing, i.e., ∂F ∂ζ ≤ 0, we may allow K to be negative, as long as C 0 is positive. Although the negative coupling strength reduce the synchronization, the effect of non-increasing forcing function F accelerate the synchronization.

4.
Practical synchronization: Heterogeneous forcing. In this section, we present several sufficient conditions for practical synchronization in terms of initial configurations, parameters and forcing terms.

Bounded forcing.
We first consider the forcings F = {F 1 , · · · , F N } satisfying the following boundedness conditions: Note that the following forcings satisfy the boundedness condition (17): Proof of claim. We split the proof into two parts. In Step A, we show that the set T is nonempty, and in Step B, we show that T * = ∞ using the differential inequality obtained in Lemma 2.4.

SEUNG-YEAL HA, SE EUN NOH AND JINYEONG PARK
Note that the condition on K guarantees that the coefficient of the second term on the r.h.s. of (19) is positive.
Then, it follows from (19) that Y (t) satisfies By Gronwall's lemma, we then have This implies On the other hand, note that the condition on K and the initial configuration are equivalent to saying that the r.h.s. of the above relation is less than or equal to D0 √ 2 .
We are now ready to provide our second main theorem by combining the results of Lemmas 2.4 and 4.1.
The proof of Theorem 2.7. We repeat the argument presented in Lemma 4.1 to derive the estimate By letting t → ∞, we obtain .
2. For the Kuramoto model with F i = Ω i and ψ ij = 1 N , we have Thus, the conditions on Ω i , K and the initial configuration ζ 0 in Lemma 4.1 reduce to which are weaker than in [7]. 3. Note that for the linear stable dynamics we have Thus, Theorem 2.7 cannot be applied to this simple case where the decoupled system has bounded solutions only. However, if we check the proof of Lemma 4.1 more carefully, what we need is boundedness over the bounded phase space, not over the whole space R for ζ i , i.e., once the coupled system (2) has only bounded solutions that are confined to the compact state space N , then we can replace 4. For a linear-time varying multi-agent systems, the practical synchronization has been studied in [17].
Below, we will show that if the uncoupled system (1) has a bounded solution for a given initial configuration, then the solution to the coupled system (2) exhibits practical synchronization.
Corollary 1. Assume that the following conditions hold.
1. The initial configuration satisfies the following boundedness condition: 2. For a given family of forcing F = {F 1 , · · · , F N }, the decoupled systeṁ has the globally bounded solution: ζ ∈ N : compact subset of R N .
3. The family of forcing F, network structure, and coupling strength satisfy the conditions: Then, the practical synchronization holds: Proof. It follows from Lemma 2.5 that the state space for ζ is bounded, so we can use the modified diameterD(F): to apply the same argument as in Theorem 2.7. This completes the proof.
Remark 5. Corollary 1 covers the case where F i is given by the gradient field of the double well potential, i.e., In this case, the solution to the decoupled system is bounded.

4.2.
Unbounded forcing. In this subsection, we consider the case where the decoupled system has bounded and unbounded solutions at the same time, and see that the unbounded solutions can be turned into bounded solutions by coupling with bounded solutions. Consider a nonlinear system with linear intrinsic dynamics F i (ζ, t) = p i ζ. Then the system (2) with all-to-all coupling ψ ij = 1 N becomes

PRACTICAL SYNCHRONIZATION ESTIMATES OF KURAMOTO OSCILLATORS 801
where p i is a constant satisfying the negative sum condition: Note that system (26) has a trivial equilibrium solution ζ e : ζ e := (0, · · · , 0).
When the coupling is turned off, i.e., K = 0, the state ζ i can go to infinity or zero exponentially fast, depending on the sign of p i . If all p i are negative, then the uncoupled dynamics have a bounded solution, and this case can be covered by Corollary 1. Thus without loss of generality, we may assume that at least one of the p i is positive. In a near-equilibrium regime ζ ≈ ζ e = 0, the dynamics of the nonlinear system (26) can be studied via the linear system near ζ e : Before we present a uniform boundedness of ζ to the linear system (28) for sufficiently large K, we consider the following dynamics for two oscillators: The linear system (29) can be rewritten in matrix form as d dt By direct calculation, the matrix M 2 has two real eigenvalues: It is easy to see that λ − < 0 for sufficiently large K. On the other hand, note that Thus, if p 1 + p 2 < 0, K > 2p1p2 p1+p2 , then both ζ 1 and ζ 2 decay to zero so that we have practical synchronization. Before we proceed to the general case, we consider the two explicit examples corresponding to the case p 1 + p 2 > 0. In this case, we will not have practical synchronization.
Again, by direct calculation, we have Thus, we can conclude that with p 1 + p 2 > 0, system (29) cannot have a practical synchronization.
We now return to the linear system (28) associated with (26) rewritten in matrix form: where Note that the coefficient matrix M N is symmetric, so that all eigenvalues of M N are real. Below, we will show that, under the condition (27), if the coupling strength K is sufficiently large, then all eigenvalues of M N become negative so that the trivial equilibrium solution to (32) is exponentially stable. For this, we recall several lemmas in relation to eigenvalues of this linearized system.
Then, the following inequalities hold: Proof. For a proof, we refer to Theorem III.2.1 in [3].
Proof. The proof is given in Appendix A.
Now, we are ready to prove the negativity of eigenvalues of M N . For a matrix A, let σ(A) denote the set of eigenvalues of A, i.e., the spectrum of A. Proposition 1. Let ζ = (ζ 1 , · · · , ζ N ) be a global solution to system (28) with negative sum condition (27). Then, for a sufficiently large K, the solution ζ converges to zero exponentially fast, independent of the initial configuration ζ 0 .
Proof. For the desired estimate, it suffices to show that the eigenvalues for the coefficient matrix M N are negative: Without loss of generality, we may assume that p 1 ≥ p 2 ≥ . . . ≥ p N . Suppose that σ(M N ) = {µ 1 , . . . , µ N } is arranged in descending order, and set H N and P N as Then, it is easy to see that the spectra of the matrices H N and P N are given in descending order: It follows from Lemma 4.2 that we have We now set K to be sufficiently large satisfying Therefore, under the condition (34), µ 2 is negative, and hence µ i , i = 3, · · · , N are also negative: Now, we need to show that µ 1 is negative. Since the determinant of a matrix is the product of all eigenvalues, it is enough to show that In Lemma 4.3, for sufficiently large K, This yields that, for sufficiently large K, It follows from (35)-(36) that all eigenvalues of M N are negative for sufficiently large K. This implies that the solution to the linear system (28) decays to zero exponentially fast.

5.
Conclusion. The dynamics of an oscillatory system over a network are often influenced by internal and external forces, e.g., the daily cycling of light and darkness affecting human sleep rhythms. In this case, due to forcing effects, all relative distances and relative rate of changes between the states of units do not asymptotically approach a constant. Thus, we can at least expect that all states are confined to some bounded region. The question is whether we can control this bounded region by the strength of the coupling strength K. Indeed, we have shown that the diameter of the state set is of order O(K −1 ) for sufficiently large K s, i.e., "practical synchronization" occurs asymptotically for sufficiently large K. In particular, Theorem 2.7 implies that, if the forcing functions are C 1 and the state space is bounded, then practical synchronization can be achieved for a sufficiently large coupling strength and some restricted class of initial configurations.
Appendix A. Proof of Lemma 4.3. In this section, we provide the proof of Lemma 4.3 using elementary row operations and the Laplace expansion. We next estimate the M i separately.