Bifurcations in the Kuramoto model on graphs

The Kuramoto model (KM) of coupled phase oscillators provides an important paradigm for studying collective dynamics and synchronization in ensembles of interacting dynamical systems. In its original form, the KM describes the dynamics of all-to-all coupled phase oscillators with randomly distributed intrinsic frequencies

Here, $θi:R→T:=R/2πZ$ is the phase of the oscillator $i∈[n]$⁠. The intrinsic frequencies $ωi,i∈[n],$ are independent identically distributed random variables. Throughout this paper, we assume that the density of the absolutely continuous probability distribution of $ω1$⁠, $g,$ is an even smooth function. The sum on the right-hand side of (1) models the interactions between oscillator $i$ and the rest of the network. Parameter $K$ controls the strength of the interactions. The sign of $K$ determines the type of interactions. The coupling is attractive if $K>0$⁠, and is repulsive otherwise. Sufficiently strong attractive coupling favors synchrony.

For small positive values of $K$⁠, the KM shows little coherence. The phases fill out the unit circle approximately uniformly [Fig. 1(a)]. This dynamical regime is called an incoherent state. It persists for positive values of $K$ smaller than the critical value $Kc=2[πg(0)]−1$⁠. For values of $K>Kc$⁠, the system undergoes a gradual build-up of coherence approaching complete synchronization as $K→∞$ [Fig. 1(b)]. Kuramoto identified the critical value $Kc$ and described the transition to synchrony, using the complex order parameter

Below, we will often use the polar form of the order parameter

Specifically, he showed that for $t≫1$ (cf. Ref. 20)

where

Here, $0≤r(t)≤1$ and $ψ(t)$ stand for the modulus and the argument of the order parameter, respectively, defined by the right-hand side of (2). The value of $r$ is interpreted as a measure of coherence in the system dynamics. Indeed, if all phases are independent random variables distributed uniformly on $T$ (complete incoherence), then $r=o(1)$ with probability one, by the Strong Law of Large Numbers. If, on the other hand, all phase variables are equal, then $r=1$⁠. Equations (4) and (5) suggest that the system undergoes a pitchfork bifurcation en route to synchronization. This bifurcation is clearly seen in the numerical experiments [see Fig. 1(c)]. Furthermore, Kuramoto’s scenario of the transition to synchrony was recently confirmed by rigorous mathematical analysis of the KM.^3,6,7

In this paper, we investigate bifurcations in the KM on a variety of graphs ranging from symmetric Cayley graphs to random small-world and power law graphs. To visualize the connectivity of a large graph, we use the pixel picture of a graph. This is a square black and white plot, where each black pixel stands for entry $1$ in the adjacency matrix of the graph (see Fig. 2). As we show below, the graph structure plays a role in the transition to synchrony. Furthermore, some graphs may force bifurcations to spatial patterns other than synchrony. To highlight the role of the network topology, we present numerical experiments of the KM on small-world graphs.²¹ These graphs are formed by replacing some of the edges of a regular graph with random edges [see Fig. 2(a)]. These results demonstrate that some graphs may exhibit bifurcations to spatial patterns other than synchrony.

Figure 3(a) shows that for the KM on small-world graphs like in that on complete graphs, the order parameter undergoes a smooth transition. Near the critical value $Kc+≈3.2$⁠, the asymptotic in time value of the order parameter starts to grow monotonically and approaches a value slightly greater than $0.5$ for increasing values of $K$⁠. This change in the order parameter corresponds to the transition to synchronization. Interestingly, in contrast to the original KM, the model on small-world graphs exhibits another transition at a negative value of $K$⁠, $Kc≈−27$ [Fig. 3(c)]. This time it grows monotonically for decreasing values of $K$ and approaches the value approximately equal to $0.0475$⁠. The corresponding state of the network is shown in Fig. 3(c). In the pattern shown in this figure, the oscillators are distributed randomly about a $2$-twisted state. A $q$-twisted state (⁠$q∈Z$⁠) is a linear function on the unit circle: $Tq(x)=2πqxnmod2π,$ $x∈[n],$ $q$-twisted states have been studied before as stable steady states in repulsively coupled KM on Cayley graphs.^17,22 By changing parameters controlling small-world connectivity, we find different $q$-twisted states bifurcating from the incoherent state for decreasing negative $K$ [see Figs. 4(a) and 4(b)]. These numerical experiments show that in the KM on small-world graphs, the incoherent state is stable for $K∈(Kc−,Kc+)$⁠. For values $K>Kc+,$ the system undergoes the transition to synchrony, while for $K<Kc−$⁠, it leads to random patterns localized around $q$-twisted states.

The loss of stability of the incoherent state in the KM on spatially structured networks was analyzed in Refs. 4 and 5. The linear stability analysis of the incoherent state yields the critical values $Kc±$⁠,⁴ and the analysis of the bifurcations at $Kc±$ explains the emerging spatial patterns.⁵ In the present work, using the insights from Refs. 4 and 5, we conduct numerical experiments illustrating the role of the network topology in shaping the transition from the incoherent state to coherent structures like synchronous and twisted states in the KM on graphs. In particular, we present numerical simulations showing for many distinct graphs the onset of synchronization takes place at the same critical value and the graph structure has only higher order effects. We identify the dominant structural properties of the graphs shaping the transition to synchronization in the KM.

The effective analysis of the KM on large graphs requires an analytically convenient description of graph sequences. To this end, we employ the approach developed in Refs. 14 and 15 inspired by the theory of graph limits¹¹ and, in particular, by W-random graphs.¹² We explain the graph models used in this paper in Sec. II. The W-random graph framework affords an analytically tractable mean field description of the KM on a broad class of dense and sparse graphs. The mean field partial differential equation for the KM on graphs is explained in Sec. III. There, we also explain the generalization of the order parameter suitable for the KM on graphs. After that we go over the main results of Refs. 4 and 5. In Sec. IV, we present our numerical results. We use carefully designed examples of graphs to highlight the structural properties affecting the onset of synchrony. For instance, we show that the KMs on Erdős-Renyi and Paley graphs have the same critical value $Kc$ marking the onset of synchronization. Furthermore, the mean field limits of the KM on these graphs coincide with that for the KM on weighted complete graph. We demonstrate the similarity in the transition to synchronization in the KM on bipartite graphs and on the family of stochastic block graphs interpolating a disconnected graph and a weighted complete graph. This time the mean field equations for the two models are different, but the critical value $Kc$ remains the same. Furthermore, we investigate the bifurcations in the KM on small-world and power law graphs illuminating the effects of these random connectivity patterns on the dynamics of the KM.

We begin with the description of graph models that will be used in this paper. Let $Γn=⟨V(Γn),E(Γn),An⟩$ be a weighted directed graph on $n$ nodes. $V(Γn)=[n]$ stands for the node set of $Γn$⁠. $An=(an,ij)$ is an $n×n$ weight matrix. The edge set

If $Γn$ is a simple graph, then $An$ is the adjacency matrix.

The KM on $Γn$ is defined as follows:

The sequence $αn≠0$ is only needed and ${Γn}$ is a sequence of sparse graphs. Without proper rescaling, the continuum limit (as $n→∞$⁠) of the KM on sparse graphs is trivial. The sequence ${αn}$ will be explained below, when we introduce sparse $W$-random graphs. Until then, one can assume $αn≡1$⁠.

We will be studying solutions of (6) in the limit as $n→∞$⁠. Clearly, such limiting behavior is possible only if the sequence of graphs ${Γn}$ is convergent in the appropriate sense. We want to define ${Γn}$ in such a way that it includes graphs of interest in applications and, at the same time, is convenient for deriving the continuum limit. To achieve this goal, we use the ideas from the theory of graph limits.¹¹ Specifically, we choose a square integrable function $W$ on a unit square $I2:=[0,1]2$⁠. $W$ is used to define the asymptotic behavior of ${Γn}$ as $n→∞$⁠. It is called a graphon in the theory of graph limits.¹¹ Let

We define $W0⊂W$ as a set of measurable symmetric functions on $I2$ with values in $I$⁠. We discretize the unit interval $I$⁠: $xn,j=j/n,j∈{0}∪[n]$ and denote $In,i:=(xn,i−1,xn,i],i∈[n].$

The following graph models are inspired by W-random graphs.^2,12 Let $Γn=⟨V(Γn),E(Γn),An=(an,ij)⟩$⁠.

• $(DD)$
  Weighted deterministic graph $Γn=H(n,W)$ is defined as follows:

  $an,ij=⟨W⟩In,i×In,j:=n2∫In,i×In,jW(x,y)dxdy,$

  (8)

  where $W∈W$⁠.
• $(RD)$
  Dense random graph $Γn=Hr(n,W)$⁠. Let $W∈W0$ and $an,ij,1≤i≤j≤n,$ be independent identically distributed (IID) random variables such that

  $Pan,ij=1=⟨W⟩In,i×In,jandPan,ij=0=1−Pan,ij=1.$

  (9)
• $(RS)$
  Sparse random graphs $Γn=Hr(n,W,{αn})$⁠. Let $W∈W$ be a nonnegative function and positive sequence $1≥αn↘0$ satisfy $nαn→∞$ as $n→∞$

  $Pi→j=αn⟨W~n⟩In,i×In,jandW~n(x,y):=min{αn−1,W(x,y)}.$

  (10)

We illustrate the random graph models in (RD) and (RS) with the following examples.

The graphon $W$ carries all the information needed for the derivation of the continuum limit for the KM on deterministic and random graphs (DD), (RD), and (RS). However, for the KM on random graphs (RD) and (RS) taking the continuum limit involves an additional step. Namely, we first approximate the model on random graph by that on the averaged deterministic weighted graph $Γn=H(n,W)$

where $W¯n,ij:=⟨W~n⟩In,i×In,j.$

The approximate model (13) is formally derived from the original KM (6) on a random graph $Γn=Hr(n,W,{αn})$ by averaging the right hand side of (6) over all possible realizations of $Γn$⁠. The justification of averaging is given in Ref. 8.

Solutions of the KM with distributed frequencies such as incoherent state and the solutions bifurcating from the incoherent state are best described in statistical terms. To this end, suppose $ρ(t,u,ω,x)$ is the conditional density of the random vector $(u,ω)$ given $ω$⁠, and parametrized by $(t,x)∈R+×I$⁠. Here, the spatial domain $I=[0,1]$ represents the continuum of the oscillators in the limit $n→∞$ (see Ref. 4 for precise meaning of the continuum limit). Then, $ρ$ satisfies the following initial value problem:

where initial density $ρ0(u)$ independent of $ω$ and $x$ for simplicity (see Ref. 8 for the treatment of a more general case). The velocity field $V$ is defined via

The density $ρ(t,⋅)$ approximates the distribution of the oscillators around $T$ at time $t∈[0,T]$⁠. Specifically, in the large $n$ limit the empirical measures defined on the Borel $σ$-algebra $B(G),$$G:=T×R$⁠:

converge weakly to the absolutely continuous measure

provided the initial data (7) converge weakly to $μ0$ [cf. Ref. 4 (Theorem 2.2), see also Ref. 8].

The mean field limit provides a powerful tool for studying the KM. It gives a simple analytically tractable description of complex dynamics in this model. In particular, the incoherent state in the mean field description corresponds to the stationary density $ρu=g(ω)/2π$⁠. The stability analysis of $ρu,$ a steady state solution of (15), yields the region of stability of the incoherent state and the critical values of $K$⁠, at which it loses stability. Furthermore, the analysis of the bifurcations at the critical values of $K$ explains the phase transitions in the KM (cf. Ref. 5). Importantly, we can trace the role of the network topology in the loss of stability of the incoherent state and emerging spatial patterns. In the remainder of this section, we informally review some of the results of Refs. 4 and 5, which will be used below.

The key ingredient in the analysis of the KM on graphs, which was not used in the analysis of the original KM, is the following kernel operator $W:L2(I)→L2(I)$⁠:

Recall that $W∈L2(I2)$ is a symmetric function. Thus, $W:L2(I2)→L2(I2)$ is a self-adjoint compact operator. The eigenvalues of $W$⁠, on one hand, carry the information about the structure of the graphs in the sequence ${Γn}$⁠; on the other hand, they appear in the stability analysis of the incoherent state. Thus, through the eigenvalues of $W$⁠, we can trace the relation between the structure of the network and the onset of synchronization in the KM on graphs.

Since $W$ is self-adjoint and compact, the eigenvalues of $W$ are real and bounded with the only possible accumulation point at $0$⁠. In all examples considered in this paper, the largest eigenvalue $μmax>0$ and the smallest $μmin≤0$⁠. The linear stability analysis in Ref. 4 shows that the incoherent state is stable for $K∈[Kc−,Kc+],$ where

Except for small-world graphs, for all other graphs considered below, the smallest eigenvalue $μmin=0$⁠. Thus, the incoherent state in the KM on these graphs is stable for $K≤Kc+$⁠. In this section, we comment on the bifurcation at $Kc+$ and postpone the discussion of the bifurcation at $Kc−$ until Sec. V, where we deal with the KM on small-world graphs.

The analysis of the bifurcations in the KM on graphs relies on the appropriate generalization of Kuramoto’s order parameter (cf. Ref. 4)

and its continuous counterpart

The value of the continuous order parameter (22)$h(t,x)$ carries the information about the (local) degree of coherence at point $x$ and time $t$⁠. It is adapted to a given network connectivity through the kernel $W$⁠.

The main challenge in the stability analysis of the incoherent state lies in the fact that for $K∈[Kc−,Kc+]$⁠, the linearized operator has continuous spectrum on the imaginary axis and no eigenvalues (cf. Ref. 4). To overcome this difficulty, in Ref. 5, the generalized spectral theory was used to identify generalized eigenvalues responsible for the instability at $Kc+$ and to construct the finite-dimensional center manifold. The explanation of these results is beyond the scope of this paper. An interested reader is referred to Ref. 5 for details. Here, we only present the main outcome of the bifurcation analysis. The center manifold reduction for the order parameter near $Kc+$ yields a stable branch of solutions

where $C(x)$ is defined through the eigenfunction of $W$ corresponding to $μmax$ (see Ref. 5 for details). Formula (23) generalizes the classical Kuramoto’s formula describing the pitchfork bifurcation in the all-to-all coupled model to the KM on graphs. The network structure enters into the description of the pitchfork bifurcation through the largest eigenvalue $μmax$ and the corresponding eigenspace.

In this section, we present numerical results elucidating some of the implications of the bifurcation analysis outlined in Sec. III.

Our first set of examples deals with the KM on Erdős-Rényi and Paley graphs (see Fig. 5). The former is a random graph, whose edges are selected at random from the set of all pairs of vertices with fixed probability $p$ (see Example II.1). Before giving the definition of the Paley graphs, we first recall the definition of a Cayley graph on an additive cyclic group $Zn=Z/nZ$⁠.

As all Cayley graphs, Paley graphs are highly symmetric [Fig. 5(b)]. However, Paley graphs also belong to pseudorandom graphs, which share many asymptotic properties with Erdős-Rényi graphs with $p=1/2$⁠.¹ In particular, both Erdős-Rényi and Paley graphs have constant graph limit $W1/2≡1/2$ as $n→∞$⁠. The mean field limit for the KM for Erdős-Rényi and Paley graphs has the same kernel constant $W=W1/2≡1/2$ [cf. (15)], exactly as in the case of the KM on weighted complete graphs. Since $μmax(W1/2)=1/2$⁠, the transition to synchronization takes place at

for all three networks (see Fig. 6).

Our second set of examples features another pair of distinct networks, which are not pseudorandom, and yet have the same synchronization threshold.

Here, we consider the KM on bipartite graphs $Bn,n=H(2n,Wbp)$ [see Fig. 7(a)]

Next, we consider the family of graphs interpolating between the disconnected graph with two equal components and a weighted complete graphs, $Γn,α=Hr(n,Wd,α)$⁠:

A simple calculation yields largest positive eigenvalues for each of these families of graphs

Thus, as in the first set of examples, the KM on $Bn,n$ and $Γn,α$ undergoes transition to synchronization at $Kc+=4[πg(0)]−1$ (see Fig. 8).

Next, we turn to the KM on power law graphs. To generate graphs with power law degree distribution, we use the method of sparse W-random graphs.² Specifically, $Γn=Hr(n,Wγ,n−β)$ with graphon $Wγ$ defined in (12). For graphs constructed using this method, it is known that the expected degree of node $i∈[n]$ is $O(n1+γ−βi−γ)$ and the edge density is $O(n−β)$ (cf. Ref. 9). Thus, $Γn=Hr(n,Wγ,n−β)$ is a family of sparse graphs with power law degree distribution. The mean field limit for the KM on $Γn$ is given by (15) with $W:=Wγ$⁠. The analysis of the integral operator $W$ [cf. (20)] with kernel $W:=Wγ$ shows that it has a single nonzero eigenvalue (cf. Ref. 13)

Thus, the synchronization threshold for the KM on power law graphs is

Note that as $α→1/2,$$Kc+$ tends to $0$⁠. Thus, the KM on power law graphs features remarkable synchronizability despite sparse connectivity. Numerics in Fig. 9 illustrate the onset of synchronization in power law networks.

Small-world graphs are obtained by random rewiring of regular $k$-nearest-neighbor networks.²¹ Consider a graph on $n$ nodes arranged around a circle, with each node connected to $k=⌊rn⌋$ neighbors from each side for some $r∈(0,0.5)$⁠. With probability $p∈[0,0.5]$⁠, each of these edges connecting $k$ neighbors from each side is then replaced by a random long-range (i.e., going outside the $k$-neighborhood) edge. The pixel picture of a representative small-world graph is shown in Fig. 2(a). A family of small-world graphs parametrized by $p$ interpolates between regular $k$-nearest-neighbor graph (⁠$p=0$⁠) and a fully random Erdős-Rényi graph (⁠$p=0.5$⁠). For intermediate values of $p,$ small-world graphs combine features of regular and random connectivity, just as seen in many real-world networks.²¹ 

It is convenient to interpret a small-world graph as a $W$-random graph $Sn,p,r=Hr(n,Wp,r)$ [cf. (11)]. The mean field limit of the KM on $Sn,p,r$ is then given by (15) with $W:=Wp,r$⁠. The corresponding kernel operator is given by

where $Kp,r(x)=(1−p)1|x|≤r(x)+p1|x|>r(x).$

Using (24), we recast the eigenvalue problem for $Wp,r$ as follows:

By taking the Fourier transform of both sides of (25), we find the eigenvalues of $Wp,r$

The corresponding eigenvectors are $wk=ei2πkx,k∈Z.$ Note that $μk=μ−k,$ since $Kp,r$ is even [cf. (26)]. Thus, the eigenspace corresponding to $μ=μ0$ is spanned by $w0=1$⁠. For $μ≠μ0,$ the eigenspace corresponding to $μ$ is spanned by

where $Iμ={k∈N:μk=μ}.$ The largest positive eigenvalue of $Wp,r$ is

Therefore, the onset of synchronization in the KM on small-world graphs takes place at

[see Fig. 3(b)].

Importantly, $Wp,r$ also has negative eigenvalues. Since $Wp,r$ is a compact operator, $0$ is the only accumulation point of the spectrum of $Wp,r$⁠. Thus, there is the smallest (negative) eigenvalue $μmin=min{μk:k∈N}$⁠. Let $q∈N$ be such that $μq=μmin$⁠. Assuming that the multiplicity of $μmin$ is $2$⁠, the center manifold reduction performed for the order parameter in Ref. 5 yields the following stable branch bifurcating from the incoherent state (⁠$h≡0$⁠) at $K=Kc−$⁠:

where $ϕ∈[0,1)$ depends on the initial data.

Equation (28) implies that at $K=Kc−$⁠, the KM on small-world graphs undergoes a pitchfork bifurcation. Unlike the bifurcation at $K=Kc+$ considered earlier, in the present case, the center manifold is two-dimensional. It is spanned by $w±q=e±i2πqx$⁠. In the remainder of this section, we analyze stable spatial patterns bifurcating from the incoherent state at $K=Kc−$⁠. To this end, we rewrite the velocity field (17) using the order parameter

The velocity field in the stationary regime is then given by

Using the polar form of the order parameter

from (30), we obtain

Since $ρ∞$ is a steady state solution of (15), it satisfies

From (32) and (33), we have

Here, $K=Kc−−κ$ and $δ$ stands for the Dirac delta function. Equation (34) describes partially phase locked solutions. They combine phase locked oscillators (type I)

and those whose distribution is given by the second line in (34) (type II). The oscillators of the first type form $q$-twisted states subject to noise $Y$⁠. Note that $Y$ is a function of $ω$ and, therefore, is random. Just near the bifurcation, where $0<κ≪1,$ the probability of $|ω|>KR∞(x,K)$ is high and, thus, most oscillators are of the second type [Figs. 10(a)]. However, as we move away from the bifurcations, we see the number of oscillators of the first type increases and the noisy twisted states become progressively more distinct [Figs. 10(b) and 10(c)]. Thus, the bifurcation analysis of the KM on small-world graphs identifies a family of stable twisted states (cf. Refs. 17 and 22). By changing parameters of small-world connectivity, one can control the winding number of the emerging twisted states [see Fig. 4(c)].

In this paper, we selected several representative families of graphs to illustrate the link between network topology and synchronization and pattern formation in the KM of coupled phase oscillators. In particular, we showed that the transition to synchronization in the KM on pseudorandom graphs (e.g., Erdős-Rényi, Paley, and complete graphs) starts at the same critical value of the coupling strength and proceeds in practically the same way. The bifurcation plots shown in Fig. 6 are very similar, although plots for the KM on Erdős-Rényi and Paley graphs show more variability. The almost identical bifurcation plots for these models are due to the fact that all three models result in the same mean field equation. This means that the empirical measures generated by the trajectories of these models are asymptotically the same limit as $n→∞$⁠. The differences seen in plots (a)–(c) of Fig. 6 are due to the higher order moments, which are not captured by the mean field limit. Other families of graphs considered in this paper include the bipartite graph and a family of stochastic block graphs interpolating between a graph with two disconnected components and Erdős-Rényi graph. The KM on all these graphs (including the disconnected graph) feature the same transition to synchronization. Finally, we studied the bifurcations in the KM on small-world and power law graphs due to their importance in applications. A remarkable feature of the KM on small-world graphs is the presence of the bifurcation leading to stable noisy twisted states. Twisted states are known as attractors in repulsively coupled KM with identical intrinsic frequencies.^16–18,22 For the KM with random frequencies from the Cauchy distribution, the transition from the incoherent state to random twisted states was shown in Ref. 19 using the Ott-Antonsen ansatz. In this paper, we identified the bifurcation in the repulsively coupled KM on small-world graphs with intrinsic frequencies from arbitrary absolutely continuous distribution with even sufficiently smooth density.

To perform numerical experiments presented in this paper, we had to overcome several challenges. Verification of the bifurcation scenarios in the KM required a large number of simulations of large systems of ordinary differential equations with random coefficients. Thus, the speed of computations was critical in this project. All computations were completed in MATLAB utilizing GPU computations for dramatic speed-up in the computational time (compared to CPU computations), see Fig. 11. Time steps were performed using Heun’s method with $Δt=0.01$⁠, which was sufficient for stable simulation results.

Each simulation was initialized with a random state vector $(un,1,…,un,n)T$ with each component independently chosen from a uniform distribution on $[0,2π)$ representing the phase of each oscillator. We additionally initialize a random vector of internal frequencies chosen from a normal distribution. Unless otherwise noted, in all simulations we take $n=4001$⁠, which is a prime $1$ modulo $4$⁠, so that the theory developed for Paley graphs applies. We observed that taking $Tfinal=20$ was sufficient for systems to exhibit synchronization or $q$-twisted states.

Computationally solving (26) for the minimal negative eigenvalue is accomplished by observing the trivial bound $μk>−1kπ$ and iteratively computing $μmin=μk0$ for $k0∈{1,…,M}$ for increasing values of $M$ until $μmin≤−1Mπ$⁠.

The work of the second author was supported in part by the National Science Foundation Division of Mathematical Sciences (NSF DMS) No. 1715161. The authors would like to thank Nicholas Battista for helpful discussions related to implementation of the numerical simulations.