The mean field analysis for the Kuramoto model on graphs I. The mean field equation and transition point formulas

In his classical work on synchronization, Kuramoto derived the formula for the critical value of the coupling strength corresponding to the transition to synchrony in large ensembles of all-to-all coupled phase oscillators with randomly distributed intrinsic frequencies. We extend the Kuramoto's result to a large class of coupled systems on convergent families of deterministic and random graphs. Specifically, we identify the critical values of the coupling strength (transition points), between which the incoherent state is linearly stable and is unstable otherwise. We show that the transition points depend on the largest positive or/and smallest negative eigenvalue(s) of the kernel operator defined by the graph limit. This reveals the precise mechanism, by which the network topology controls transition to synchrony in the Kuramoto model on graphs. To illustrate the analysis with concrete examples, we derive the transition point formula for the coupled systems on Erd\H{o}s-R\'{e}nyi, small-world, and $k$-nearest-neighbor families of graphs. As a result of independent interest, we provide a rigorous justification for the mean field limit for the Kuramoto model on graphs. The latter is used in the derivation of the transition point formulas.


Introduction
Synchronization of coupled oscillators is a classical problem of nonlinear science with diverse applications in science and engineering [3,38]. Physical and technological applications of synchronization include power, sensor, and communication networks [9], mobile agents [35], electrical circuits [1], coupled lasers [18], and Josephson junctions [44], to name a few. In biological and social sciences, synchronization is studied in the context of flocking, opinion dynamics, and voting [14,30]. Synchronization plays a prominent role in physiology and in neurophysiology, in particular. It is important in the information processing in the brain [36] and in the mechanisms of several severe neurodegenerative diseases such as epilepsy [43] and Parkinsons Disease [22]. This list can be continued.
Identifying common principles underlying synchronization in such diverse models is a challenging task. In seventies Kuramoto found an elegant approach to this problem. Motivated by problems in statistical physics and biology, he reduced a system of weakly coupled limit cycle oscillators to the system of equations for the phase variables only 1 . The resultant equation is called the Kuramoto model (KM) [20]. Kuramoto's method applies directly to a broad class of models in natural science. Moreover, it provides a paradigm for studying synchronization. The analysis of the KM revealed one of the most striking results of the theory of synchronization. For a system of coupled oscillators with randomly distributed intrinsic frequencies, Kuramoto identified the critical value of the coupling strength, at which the gradual buildup of coherence begins. He introduced the order parameter, which describes the degree of coherence in a coupled system. Using the order parameter, Kuramoto predicted the bifurcation marking the onset of synchronization.
Kuramoto's analysis, while not mathematically rigorous, is based on the correct intuition for the transition to synchronization. His discovery initiated a line of fine research (see [40,41,39,5] and references therein). It was shown in [40,41] that the onset of synchronization corresponds to the loss of stability of the incoherent state, a steady state solution of the mean field equation. The latter is a nonlinear hyperbolic partial differential equation for the probability density function describing the distribution of phases on the unit circle at a given time. The bifurcation analysis of the mean field equation is complicated by the presence of the continuous spectrum of the linearized problem on the imaginary axis. To overcome this problem, in [5] Chiba developed an analytical method, which uses the theory of generalized functions and rigged Hilbert spaces [12].
The Kuramoto's result and the analysis in [40,41,5] deal with all-to-all coupled systems. Real world applications feature complex and often random connectivity patterns [34]. The goal of our work is to describe the onset of synchronization in the KM on graphs. To this end, we follow the approach in [27,28]. Specifically, we consider the KM on convergent families of simple and weighted graphs, for which we derive and rigorously justify the mean field limit. Our framework covers many random graphs widely used in applications, including Erdős-Rényi and small-world graphs. With the mean field equation in hand, we derive the transition point formulas for the critical values of the coupling strength, where the incoherent state loses stability. Thus, we identify the region of linear stability of the incoherent state. In the follow-up work, we will show that in this region the incoherent state is asymptotically stable (albeit in weak topology), analyze the bifurcations at the transition points, and develop the center manifold reduction.
The original KM with all-to-all coupling and random intrinsic frequencies has the following form: Here, θ i : R → S := R/2πZ, i ∈ [n] := {1, 2, . . . , n} is the phase of the oscillator i, whose intrinsic frequency ω i is drawn from the probability distribution with density g(ω), n is the number of oscillators, and K is the strength of coupling. The sum on the right-hand side of (1.1) describes the interactions of the oscillators in the network. The goal is to describe the distribution of θ i (t), i ∈ [n], for large times and n ≫ 1.
Since the intrinsic frequencies are random, for small values of the coupling strength K > 0, the dynamics of different oscillators in the network are practically uncorrelated. For increasing values of K > 0, however, the dynamics of the oscillators becomes more and more synchronized. To describe the degree of synchronization, Kuramoto used the complex order parameter: (1.2) Here, 0 ≤ r(t) ≤ 1 and ψ(t) stand for the modulus and the argument of the order parameter defined by the right-hand side of (1.2). Note that if all phases are independent uniform random variables and n ≫ 1 then with probability 1, r = o(1) by the Strong Law of Large Numbers. If, on the other hand, all phase variables are equal then r = 1. Thus, one can interpret the value of r as the measure of coherence in the system dynamics. Numerical experiments with the KM (with normally distributed frequencies ω i 's) reveal the phase transition at a certain critical value of the coupling strength K c > 0. Specifically, numerics suggest that for t ≫ 1 (cf. [39]) Assuming that g is a smooth even function that is decreasing on ω ∈ R + , Kuramoto derived the formula for the critical value K c = 2 πg(0) .
Furthermore, he formally showed that in the partially synchronized regime (K > K c ), the steady-state value for the order parameter is given by Recently, Chiba and Nishikawa [7] and Chiba [5] confirmed Kuramoto's heuristic analysis with the rigorous derivation of (1.3) and analyzed the bifurcation at K c .
In this paper, we initiate a mathematical investigation of the transition to coherence in the Kuramoto model on graphs. To this end, we consider the following model: is an n × n symmetric matrix of weights. Note that in the classical KM (1.1), every oscillator is coupled to every other oscillator in the network, i.e., the graph describing the interactions between the oscillators is the complete graph on n nodes. In the modified model (1.5), we supply the edges of the complete graph with the weights (W nij ). Using this framework, we can study the KM on a variety of deterministic and random (weighted) graphs. For instance, let W nij , 1 ≤ i < j ≤ n be independent Bernoulli random variables P(W nij = 1) = p, for some p ∈ (0, 1). Complete the definition of W n by setting W nji = W nij and W nii = 0, i ∈ [n]. With this choice of W n , (1.5) yields the KM on Erdős-Rényi random graph.
The family of Erdős-Rényi graphs parameterized by n is one example of a convergent family of random graphs [23]. The limiting behavior of such families is determined by a symmetric measurable function on the unit square W (x, y), called a graphon. In the case of the Erdős-Rényi graphs, the limiting graphon is the constant function W ≡ p. In this paper, we study the KM on convergent families of deterministic and random graphs. In each case the asymptotic properties of graphs are known through the limiting graphon W . The precise relation between the graphon W and the weight matrix W n will be explained below.
In studies of coupled systems on graphs, one of the main questions is the relation between the structure of the graph and network dynamics. For the problem at hand, this translates into the question of how the structure of the graph affects the transition to synchrony in the KM. For the KM on convergent families of graphs, in this paper, we derive the formulas for the critical values We derive (1.6) from the linear stability analysis of the mean field limit of (1.5). The latter is a partial differential equation for the probability density function corresponding to the distribution of the phase variables on the unit circle (see (2.10), (2.11)). For the classical KM, the mean field limit was derived by Strogatz and Mirollo in [40]. We derive the mean field limit for the KM on weighted graphs (1.5) and show that its solutions approximate probability distribution of the phase variables on finite time intervals for n ≫ 1. Here, we rely on the theory of Neunzert developed for the Vlasov equation [31,32] (see also [4,8,13]), which was also used by Lancellotti in his treatment of the mean field limit for the classical KM [21].
With the mean field limit in hand, we proceed to study transition to coherence in (1.5). As for the classical KM, the density of the uniform distribution is a steady state solution of the mean field limit. The linear stability analysis in Section 3 shows that the density of the uniform distribution is neutrally stable for K − c ≤ K ≤ K + c and is unstable otherwise. Thus, the critical values K ± c given in (1.6) mark the loss of stability of the incoherent state. The bifurcations at K ± c and the formula for the order parameter corresponding to (1.4) will be analyzed elsewhere using the techniques from [5,6].
Sections 4 and 5 deal with applications. In the former section we collect approximation results, which facilitate application of our results to a wider class of models. Further, in Section 5, we discuss the KM for several representative network topologies : Erdős-Rényi, small-world, k-nearest-neighbor graphs, and the weighted ring model. We conclude with a brief discussion of our results in Section 6.

The mean field limit
Throughout this paper, we will use a discretization of I = [0, 1] : which satisfies the following property 3) The weighted graph Γ n = G(W, X n ) on n nodes is defined as follows. The node and the edge sets of Γ n are V (Γ n ) = [n] and respectively. Each edge {i, j} ∈ E(Γ n ) is supplied with the weight W nij := W (ξ ni , ξ nj ).
On Γ n , we consider the KM of phase oscillatorṡ The phase variable θ ni : R → S := R/2πZ corresponds to the oscillator at node i ∈ [n]. Throughout this paper, we identify θ ∈ S with its value in the fundamental domain, i.e., θ ∈ [0, 2π). Further, we equip S with the distance The oscillators at the adjacent nodes interact through the coupling term on the right hand side of (2.5). The intrinsic frequencies ω 1 , ω 2 , . . . are IID RVs. Assume that ω 1 has absolutely continuous probability distribution with a continuous density g(ω). The initial condition are sampled independently from the conditional probability distributions with densitiesρ 0 θ|ω (θ, ω i , ξ ni ), i ∈ [n]. Here,ρ 0 θ|ω (θ, ω, ξ) is a nonnegative continuous function on G := S × R × I that is uniformly continuous in ξ, i.e., ∀ǫ > 0 ∃δ > 0 such that We want to show that the dynamics of (2.5) subject to the initial condition (2.7) can be described in terms of the probability density functionρ(t, θ, ω, x) satisfying the following Vlasov equation and the initial conditionρ By (2.9),ρ(0, θ, ω, x) is a probability density on (G, B(G)): (2.13) Here, B(G) stands for the Borel σ-algebra of G.
Below, we show that the solutions of the IVPs for (2.5) and (2.10), generate two families of Borel probability measures parametrized by t > 0. To this end, we introduce the following empirical measure To compare measures generated by the discrete and continuous systems, following [32], we use the bounded Lipschitz distance: where L is the set of functions and M stands for the space of Borel probability measures on G. Here, for P = (θ, ω, x) and P ′ = (θ ′ , ω ′ , x ′ ). The bounded Lipschitz distance metrizes the convergence of Borel probability measures on G [10, Theorem 11.3.3].
We are now in a position to formulate the main result of this section.
Proof. We rewrite (2.5) as followṡ subject to the initial condition (2.20) As before, we consider the empirical measure corresponding to the solutions of (2.19), (2.20) We need to show that µ n t and µ t are close for large n. This follows from the Neunzert's theory [32].
for some C > 0 independent from n.
For a given ν . ∈ C(0, T ; M), consider the following equation of characteristics: where P = (θ, ω, x) and Under our assumptions on W , (2.24) has a unique global solution, which depends continuously on initial data. Thus, (2.25) generates the flow T t,s : Following [31], we consider the fixed point equation: It is shown in [31] that under the conditions (I) and (II) given below, for any ν 0 ∈ M there is a unique solution of the fixed point equation (2.26) ν . ∈ C(0, T ; M). Moreover, for any two initial conditions ν 1,2 0 ∈ M, we have sup for some C > 0. By construction of T t,s and (2.21), the empirical measure µ n . satisfies the fixed point equation (2.26). By [31,Theorem 1], ν t , the solution of the (2.26), is an absolutely continuous measure with densityρ(t, ·) for every t ∈ [0, T ], provided ν 0 is absolutely continuous with densityρ(0, ·) (cf. (2.13)). Furthermore,ρ(t, P ) is a weak solution of the IVP for (2.10), (2.11), and (2.12). Therefore, since both the empirical measure µ n . and its continuous counterpart µ . (cf. (2.21) and (2.17)) satisfy the fixed point equation (2.26), we can use (2.27) to obtain (2.22). It remains to verify the following two conditions on the vector fieldṼ [ν . ], which guarantee the solvability of (2.26) and continuous dependence on initial data estimate (2.27) (cf. [32]): ](t, P ) is continuous in t and is globally Lipschitz continuous in P with Lipschitz constant 3 L 1 , which depends on W .
] is Lipschitz continuous in the following sense: for some L 2 > 0 and for all µ . , ν . ∈ C(R, M) and (t, For the Lipschitz continuous function W , it is straightforward to verify conditions (I) and (II). In particular, (I) follows from the weak continuity of µ t (cf. (2.23)) and Lipschitz continuity of W and sin x. The second condition is verified following the treatment of the mechanical system presented in [32] (see also [21]). We include the details of the verification of (II) for completeness.
Let P = (θ, ω, x) ∈ G be arbitrary but fixed and define which verifies the condition (II). where BL(G) stands for the space of bounded real-valued Lipschitz functions on G with the supremum norm. Since BL(G) is a separable space, let {f m } ∞ m=1 denote a dense set in BL(G). Using (2.7) and (2.14), we have RVs Y m,ni , i ∈ [n], are independent and uniformly bounded. Further, Because f m is Lipschitz and ρ 0 θ|ω is uniformly continuous in ξ (cf. (2.8)), the function is continuous on I.

Linear stability
In the previous section, we established that the Vlasov equation (2.10), approximates discrete system (2.5) for n ≫ 1. Next, we will use (2.10) to characterize the transition to synchrony for increasing K. To this end, in this section, we derive the linearized equation about the incoherent state, the steady state solution of the mean field equation. In the next section, we will study how the spectrum of the linearized equation changes with K.
In this section, we assume that probability density g is a continuous and even function monotonically decreasing on R + .

Linearization
First, settingρ(t, θ, ω, x) = ρ(t, θ, ω, x)g(ω), from (2.10) we derive the equation for ρ: By integrating (2.10) over S, one can see that and, thus, In addition, Let The expression in the curly brackets has the following form: Thus, where the linear operator G is defined by We arrive at the linearized equation:

Fourier transform
We expand Z into Fourier series whereẐ k stands for the Fourier transform of Ẑ Z k = 1 2π S Z(θ, ·)e ikθ dθ.
The linear stability of ρ u is, thus, determined by the time-asymptotic behavior ofẐ k , k ≥ 1. To derive the differential equations forẐ k , k ≥ 1, we apply the Fourier transform to (3.8): Using the definition of the Fourier transform and integration by parts, we have It remains to compute (G[Z]) k , k ≥ 1. To this, end we rewrite (3.12) Using (3.12), we compute where P[Z] := I R W (x, y)Z(t, λ, y)g(λ)dλdy. (3.14) The combination of (3.10), (3.11), and (3.13) yields the system of equations forẐ k , k ≥ 1:

Spectral analysis
In this section, we study (3.15), which decides the linear stability of the mixed state. We rewrite (3.15) as Equations (3.14) and (3.18) define linear operators P and T on the weighted Lebesgue space L 2 (X, gdωdx) with X := R × I. Proof. Consider the multiplication operator M iω : L 2 (X, gdωdx) → L 2 (X, gdωdx) defined by It is well known that M iω is closed and the (continuous) spectrum of M iω lies on the imaginary axis σ c (M iω ) = i · supp(g). Since W (x, y) is square-integrable by the assumption, P is a Hilbert-Schmidt operator on L 2 (X, gdωdx) and, therefore, is compact [46]. Then, the statement of the lemma follows from the perturbation theory for linear operators [17].
We define a Fredholm integral operator W on L 2 (I) by (3.20) Since W is compact, the set of eigenvalues σ p (W) is a bounded countable with the only accumulation point at the origin. Since W is symmetric, all eigenvalues are real numbers.

Lemma 3.2. The eigenvalues of T are given by
Proof. Suppose v ∈ L 2 (X, gdωdx) is an eigenvector of T corresponding to the eigenvalue λ: Then By multiplying both sides of (3.23) by g(ω) and integrating with respect to ω, we have where V := R v(ω, ·)g(ω)dω ∈ L 2 (I).
If 0 ∈ σ p (W) and V is a corresponding eigenvector, then Equation (3.23) yields Thus, ζ = 0 is not an eigenvalue of T .
The second equation of (3.29) gives (3.30) Since g is even, y = 0 is a solution of (3.30). Furthermore, since g is unimodal, there are no other solutions of (3.30). Thus, λ is real.
For the uniqueness, it is sufficient to show that the function is monotonically decreasing in x except at the jump point x = 0. If this were not true, there would exist two eigenvalues for some interval K 1 < K < K 2 , and two eigenvalues would collide and disappear at K = K 1 or K = K 2 . This is impossible, because D(λ) is holomorphic in λ.
To formulate the main result of this section, we will need the following notation: . Theorem 3.4 follows from Lemma 3.3. It shows that the incoherent state is linearly (neutrally) stable for K ∈ [K − c , K + c ] and is unstable otherwise. The critical values K ± c mark the loss of stability the incoherent state. The detailed analysis of the bifurcations at K ± c will be presented elsewhere. Remark 3.5. For the classical Kuramoto model on the complete graph, W ≡ 1 and ζ max (W) = 1. Thus, we recover the well-known Kuramoto's transition formula (1.3) [20]. In the general case, the transition points depend on the graph structure through the extreme eigenvalues of the kernel operator W. Remark 3.6. For nonnegative graphons W , ζ max (W) coincides with the spectral radius of the limiting kernel operator W (cf. [42]):

This can be seen from the variational characterization of the eigenvalues of a self-adjoint compact operator, which also implies
(W[f ], f ).

Approximation
Equation (2.5) may be viewed as a base model. To apply our results to a wider class of deterministic and random networks, we will need approximation results, which are collected in this section.

Deterministic networks
Consider the Kuramoto model on the weighted graphΓ n = [n], E(Γ n ),W whereW n = (W nij ) is a symmetric matrix.
Denote the corresponding empirical measure bỹ First, we show that if W n andW n are close, so are the solutions of the IVPs for (2.5) and (4.1) with the same initial conditions. To measure the proximity of W n = (W nij ) ∈ R n×n andW n = (W nij ) and the corresponding solutions θ n andθ n , we will use the following norms: where θ n = (θ n1 , θ n2 , . . . , θ nn ) and W n = (W nij ).
The following lemma will be used to extend our results for the KM (2.5) to other networks.

Random networks
We now turn to the KM on random graphs. To this end, we use W-random graphΓ n = G r (X n , W ), which we define next. As before, X n is a set of points (2.1),(2.2).Γ n is a graph on n nodes, i.e., V (Γ n ) = [n]. The edge set is defined as follows: The decision for each pair {i, j} is made independently from the decisions on other pairs.
The KM on the W-random graphΓ n = G r (X n , W ) has the following form: where e nij , 1 ≤ i ≤ j ≤ n are independent Bernoulli RVs: and e nij = e nji . Proof. The proof follows the lines of the proof of Lemma 4.1. As before, we set up the equation for φ ni :=θ ni − θ ni : As in (4.9), we have bound where we used 0 ≤ e nij ≤ 1.
Next, we turn to the first term on the right hand side of (4.15). For this, we will need the following definitions: and Z n = (Z n1 , Z n2 , . . . , Z nn ). With these definitions in hand, we estimate I 1 as follows: The combination of (4.15)-(4.17) yields, Our next goal is to estimate  To this end, we will use the following observations. Note that η nik and η nil are independent for k = l and where we used P(e nij = 1) = W nij .
Further, and c nikl η nik η nil (4.26) and, finally, We have six summation indices i, k, l, j, p, q ranging from 1 to n. Since E η nik = 0 for i, k ∈ [n], and RVs η nik and η njp are independent whenever {i, k} = {j, p}, the nonzero terms on the right-hand side of (4.27) fall into two groups: There are n 2 terms of type I and 3n 3 (n − 1) terms of type II. Thus, where we used (4.27), (4.22), (4.23), and the bound on |c nikl | in (4.25).
Next, for a given ǫ > 0 we denote and use Markov's inequality and (4.28) to obtain

Erdős-Rényi graphs
Let p ∈ (0, 1), X n be defined in (2.1), (2.2), and W p (x, y) ≡ p. Then Γ n,p = G r (X n , W p ) is a family of Erdős-Rényi random graphs. To apply the transition point formula (3.35), we need to compute the largest eigenvalue of the self-adjoint compact operator W p : L 2 (I) → L 2 (I) defined by Proof. Suppose λ ∈ R\{0} is an eigenvalue of W p and v ∈ L 2 (I) is the corresponding eigenvector. Then Since the right-hand side is not identically 0, v = constant = 0. By integrating both sides of (5.2), we find that ζ = p.
Thus, for the KM on Erdős-Rényi graphs, we have

Small-world network
Small-world (SW) graphs interpolate between regular nearest-neighbor graphs and completely random Erdős-Rényi graphs. They found widespread applications, because they combine features of regular symmetric graphs and random graphs, just as seen in many real-world networks [45].
Following [28,29], we construct SW graphs as W-random graphs [24]. To this end let X n be a set of points from (2.1) satisfying (2.2), and define W p,r : I 2 → I by where p, r ∈ (0, 1/2) are two parameters.
The justification of the mean field limit in Section 2 relies on the assumption that the graphon W is Lipschitz continuous. To apply Theorem 2.2 to the model at hand, we approximate the piecewise constant graphon W p,r by a Lipschitz continuous W ǫ p,r : Further, we approximate the SW graph Γ n = G r (X n , W p,r ) by Γ ǫ n = G r (X n , W ǫ p,r ). The approximation results in Section 4 guarantee that the empirical measures generated by the solutions of the IVPs for the KMs on Γ n and Γ ǫ n (with the same initial data) are O(ǫ) close with high probability for sufficiently large n. Thus, below we derive the transition point formula for the KM on (Γ ǫ n ) and take the limit as ǫ → 0 to obtain the critical values for the KM on SW graphs (Γ n ).
where K p,r is a 1−periodic function on R, whose restriction to the interval [−1/2, 1/2) is defined as follows The eigenvalue problem for W p,r can be rewritten as Thus, the eigenvalues of W p,r are given by the Fourier coefficients of K p,r (x) as for k ∈ Z. The corresponding eigenvectors v k = e i2πkx , k ∈ Z, form a complete orthonormal set in L 2 (I).

Coupled oscillators on a ring
A common in applications type of network connectivity can be described as follows. Consider n oscillators placed uniformly around a circle. They are labelled by integers from 1 to n in the counterclockwise direction.
To each potential edge {i, j} ∈ [n] 2 we assign a weight where G is a 1-periodic even measurable bounded function.
where r ∈ (0, 1/2) is a parameter. With this choice of G := G r , we obtain a k-nearest-neighbor model, in which each node is connected to k = ⌊rn⌋ from each side.
Example 5.5. Another representative example was used by Kuramoto and Battogtokh [19]. Here, let X n = {1/n, 2/n, . . . , 1} and the restriction of the 1-periodic even function G to [0, 1/2] is defined by G(x) := e −κx , where κ > 0 is a parameter. With this choice of G, we obtain the KM where the strength of interactions decreases exponentially with the distance between oscillators.
As in our treatment of the KM on SW graphs in §5.2, the integral operator W : L 2 (I) → L 2 (I) can be written as a convolution W[f ](x) = I G(x − y)f (y)dy, f ∈ L 2 (I). (5.12) It is easy to verify that the eigenvalues of W are given by the Fourier coefficients of G(x).
For instance, for the k-nearest-neighbor network in Example 5.4, by setting p = 0 in (5.8), for the network at hand we obtain ζ max (W) = r and, thus, Note that, in accord with our physical intuition, the transition point is inversely proportional to the coupling range.
If the explicit expression for the largest positive eigenvalue of W is not available, for a network with nonnegative graphon W , the transition point can be estimated using the variational characterization of the largest eigenvalue (3.36). Specifically, from (3.36), we have: , and, thus, K + c ≤ 2 πg(0) W L 1 (I 2 ) .

Discussion
In this work, we derived and rigorously justified the mean field equation for the KM on convergent families of graphs. Our theory covers a large class of coupled systems. In particular, it clarifies the mathematical meaning of the mean field equation used in the analysis of chimera states (see, e.g., [33]). Moreover, we show how to write the mean field equation for the KM on many common in applications random graphs including Erdős-Rényi and small-world graphs, for which it has not been known before.
We used the mean field equation to study synchronization in the KM on large deterministic and random graphs. We derived the transition point formulas for the critical values of the coupling strength, at which the incoherent state looses stability. The transition point formulas show explicit dependence of the stability boundaries of the incoherent state on the spectral properties of the limiting graphon. This reveals the precise mechanism by which the network topology affects the stability of the incoherent state and the onset of synchronization. In the follow-up work, we will show that the linear stability analysis of this paper can be extended to show nonlinear stability of the incoherent state albeit with respect to the weak topology. There we will also present the bifurcation analysis for the critical values K ± c . The analysis of the KM on small-world networks shows that, unlike in the original KM (1.1), on graphs the incoherent state may remain stable even for negative values of K, i.e., for repulsive coupling. In fact, the bifurcations at the left and right endpoints can be qualitatively different. In the small-world case, the center manifold at K − c is two-dimensional, whereas it is one-dimensional at K + c . These first findings indicate that the bifurcation structure of the KM on graphs (2.5) is richer than that of its classical counterpart (1.1) and motivates further investigations of this interesting problem. In the future, we also plan to extend our analysis to the KM on certain sparse graphs, including sparse power law networks considered in [16].