Competitive Suppression of Synchronization and Non-Monotonic Transitions in Oscillator Communities with Distributed Time Delay

Community structure and interaction delays are common features of ensembles of network coupled oscillators, but their combined effect on the emergence of synchronization has not been studied in detail. We study the transitions between macroscopic states in coupled oscillator systems with community structure and time delays. We show that the combination of these two properties gives rise to non-monotonic transitions, whereby increasing the global coupling strength can both inhibit and promote synchronization, yielding both desynchronization and synchronization transitions. For relatively wide parameter choices we also observe asymmetric suppression of synchronization, where communities compete to suppress one another's synchronization properties until one or more win, totally suppressing the others to effective incoherence. Using the ansatz of Ott and Antonsen we provide analytical descriptions for these transitions that confirm numerical simulations.


I. INTRODUCTION
Understanding the emergence of collective behavior in ensembles of interacting dynamical systems remains an important area of research in the nonlinear dynamics community due to synchronization's central role in a wide range of phenomena [1,2]. Examples from both natural and engineered systems include cardiac pacemaker dynamics [3], power grids [4], and Josephson junctions [5]. The Kuramoto phase oscillator model and its variants have proven particularly useful in building an understanding of collective behavior [6], and a large body of literature has identified features that give rise to rich nonlinear behavior, including external forcing [7], multimodal frequency distributions [8], and mixed sign coupling [9].
In many applications, oscillators are organized into well defined communities, where oscillators in the same community are more strongly coupled with each other than with oscillators in other communities. Examples include engineered communities of synchronizing bacteria (so-called biopixels) [10,11], interconnected regional or national power grid networks [12,13], microgrids [14], and synchronization of neuronal oscillations from different brain regions [15,16]. In addition to the topological effects that community structure has on synchronization dynamics, time delays in the transmission of information from one community to another inevitably exist when network structure is related to an underlying geometry, as in the examples listed above. While the effects of community structure and time delay on the synchronization dynamics have been studied separately in the context of the Kuramoto model [17][18][19][20][21][22][23], their combined effect on the collective dynamics remains, with a few exceptions (e.g., Ref. [24]), relatively unexplored. Here we address this gap and uncover a number of novel nonlinear behaviors. * juanga@colorado.edu; The authors contributed equally to this work. † persebastian.skardal@trincoll.edu ; The authors contributed equally to this work. The most notable phenomenon that arises from the combination of time delays and community structure is a number of non-monotonic synchronization transitions, where increasing the global coupling strength can either inhibit or promote synchronization. In the case of two communities (illustrated in Fig. 1) this manifests first in a desynchronization transition where locally synchronized states give way to incoherence, followed by a (subcritical) synchronization transition where incoherence gives way to global synchronization. For more than two communities a third transition occurs with incoherence giving way to local synchronization, which is then followed by the transitions described above. In addition to these non-monotonic transitions, when system the communities' parameters are chosen to be asymmetric we find that the oscillator communities compete and asymmetrically suppress one another until one or more "win out", totally suppressing the others' synchronization properties. The role played by each community, i.e., which wins out and is able to remain synchronized longer, depends nonlinearly on the time delay.

II. OSCILLATOR COMMUNITIES WITH DISTRIBUTED TIME DELAY
We consider here a system of C ≥ 2 communities of coupled phase oscillators governed bẏ where θ σ i represents the phase of oscillator i in community σ, ω σ i is its natural frequency, which is assumed to be drawn from the distribution g σ (ω), K σσ ′ is the coupling strength between oscillators in communities σ and σ ′ , τ σσ ij is the time delay between oscillators i and j in communities σ and σ ′ , which is assumed to be drawn from the distribution h σσ ′ (τ ), and N σ is the number of oscillators in community σ.
To begin our analysis of the dynamics of Eq. (1) we make a few simplifying parameter choices. First we allow for two coupling strengths: k = K σσ and K = CK σσ ′ /2 (for σ ′ = σ) denoting intra-and inter-community coupling. Next, we assume that within each community time delays are zero, i.e., h σσ (τ ) = δ(τ ), but between different communities the distribution h σσ ′ is exponential with mean T σ ′ , namely, h σσ ′ (τ ) ∝ e −τ /T σ ′ . We also consider the case where all communities are of the same size, i.e., N σ = N for all σ and we assume that frequency distributions are Lorentzian with comunity-specific mean Ω σ and width ∆, i.e., g σ (ω) = ∆/{π[∆ 2 +(ω−Ω σ ) 2 ]}. Next, seeking a description for the local order parameters z σ = N −1 N j=1 e iθ σ j describing the degree of synchronization within each community, we apply the dimensionality reduction technique of Ott and Antonsen [25,26] to obtain the following system of reduced equationṡ where w σ represents a time-delayed order parameter for community σ, with σ = 1, 2, . . . , C and the width parameter ∆ of the natural frequency distributions have been scaled out. (Details of the dimensionality reduction are provided in Appendix A).

A. Two communities
To illustrate the rich dynamics introduced by the interplay between time delay and hierarchical community structure in the simplest setting, we consider the case of two communities with T 1 = T 2 = T , illustrated in Fig. 1. We also define κ = (k − 2)/2, noting that for K = 0 each community undergoes a transition to synchronization at k = 2. Therefore κ > 0 (κ < 0) indicates that the isolated communities would be synchronized (incoherent). Equations (2)-(3) becomė Equations (4)-(6) display some remarkable dynamics, which we illustrate in Fig. 2 for parameters κ = 0.02, Ω 1 = −1, Ω 2 = 2, and T = 1, T = 0.1 for panels (a) and (b), respectively. We plot the time-averaged values of r 1 = |z 1 | (solid blue) and r 2 = |z 2 | (dashed red) obtained from first slowly increasing K from 0 to 4, then slowly decreasing it back to 0. (For some parameters, r 1 and r 2 have oscillatory or chaotic dynamics, which are not the focus of this paper.) For sufficiently small K both communities are partially synchronized with r 1 , r 2 > 0 but are not synchronized with one another as the angles ψ 1 , ψ 2 do not phase-lock (not shown). We call this state local synchronization. As K is increased we observe that one of the communities nearly desynchronizes (r 1 ≈ 0 or r 2 ≈ 0) while the other remains synchronized, a state we call asymmetric suppression. (Although complete incoherence is not reached because of the pulling effect of the synchronized community, the transition is easy to see.) More specifically, for T = 1 community 1 undergoes this transition, whereas for T = 0.1 it is community 2 that undergoes this transition. Since the oscillator communities differ only in their mean frequencies Ω σ and have the same spread ∆ = 1 (which typically dictates local synchronization properties) the identity of the desynchronizing community depends on the characteristic time delay rather than any community-specific properties.
In particular, while both communities suppress one another's synchronization properties, one eventually wins out, remaining synchronized for larger values of K. Further increasing K yields complete incoherence, i.e., both r 1 , r 2 = 0, followed by a subcritical (sometimes called "explosive" [27]) transition to global synchronization. When K is decreased the system undergoes a similarly explosive desynchronization transition to incoherence at a different coupling strength (highlighting the existence of a hysteresis loop), followed by a return to lo- cal synchronization through the asymmetric suppression state. We note that the results of the reduced equations presented here are in excellent agreement with those obtained from direct simulations of a microscopic system based on Eq. (1) using communities of size N = 2 × 10 4 , which are presented in Appendix B, and similarly for results presented below.
To better understand this sequence of bifurcations, we perform a linear stability analysis of the incoherent state z 1 = z 2 = w 1 = w 2 = 0. Details can be found in Appendix C. Linearizing Eqs. (4)-(6) and looking for values of K where solutions have a purely imaginary growth rate, z 1 =z 1 e iωt , z 2 =z 2 e iωt , w 1 =w 1 e iωt , w 2 =w 2 e iωt , we obtain the following equations for ω and K: In what follows we will focus on the case 0 < κ ≪ 1. (The case of general κ is treated exactly for the symmetric case Ω 1 = −Ω 2 below). To balance Eq. (8) when κ ≪ 1 there are three options, to first order in κ: ω = Ω 1 +ω 1 κ, ω = Ω 2 +ω 2 κ, and ω = ω 3 κ, where ω i are constants to be determined. Inserting these in Eqs. (7)-(8) we find the corresponding values of K, to leading order in κ: These values of K are real and positive if the frequencies Ω 1 and Ω 2 have opposite sign. In general, depending on the values of Ω 1 , Ω 2 , and T , one can have any ordering of K 1 , K 2 , and K 3 , leading to different bifurcation structures. In Fig. 3 we obtain the stability diagram for the system by plotting K 1 (solid blue), K 2 (dashed red), and K 3 (dot-dashed black) as a function of T for κ = 0.02, Ω 1 = −1, and Ω 2 = 2.
The sequence of bifurcations shows regions of stability, with local synchronization when K < K 1 , K 2 , K 3 , global synchronization when K > K 3 , asymmetric suppression where (We note that the region of stability for global synchronization typically stretches below K 3 , which we will see below.) Regions of stability are labeled in Fig. 3 and a zoomed-in view of the small-T region is shown in panel (b). We also indicate the time delays T = 1 and T = 0.1 used in Figs. 2(a) and (b) using vertical dotted lines. Lastly, the interplay between K 1 , K 2 , and K 3 illuminates the asymmetric suppression state as follows. First, asymmetric suppression states are only attainable for time delays T where min(K 1 , K 2 ) < K 3 . They are then realized when K surpasses either K 1 or K 2 , but not both. The mode that becomes unstable at K = K 1 (K = K 2 ) is localized in community 1 (2). More precisely, at K = K 1,2 the mode with imaginary growth rate satisfies r 2,1 ∝ κ 1/2 r 1,2 (see Appendix D), so that for κ ≪ 1 the synchronized mode is localized in one community. When K 1 < K < K 2 or K 2 < K < K 1 , we find that the values of r 1 , r 2 saturate at values consistent with the mode localization predicted by the linear stability analysis. In particular, when K 1 < K < K 2 community 1 reaches near incoherence, whereas when K 2 < K < K 1 community 2 reaches near incoherence. Moreover, communities 1 and 2 swap roles in the asymmetric suppression state when K 1 and K 2 intersect at the critical time de- , at which point no asymmetric suppression state exists and the local synchronization state gives way directly to incoherence. Assuming that |Ω 2 | > |Ω 1 | this shows that for sufficiently small T it is community 1 that desynchronizes first, while for larger T it is community 2 that desynchronizes first.
To gain further insight into the hysteretic nature of the transition to global synchronization we now focus on the sym-metric case, Ω 2 = −Ω 1 = Ω. Searching for phase-locked solutions of Eqs. (4)-(6) of the form z 1 = re iψ , z 2 = re iψ+α , with r, ψ, α constants and w 1,2 = z 1,2 , gives the following implicit expression for the synchronized branch: By taking the limit r → 0 + we see that this branch begins at K = K 3 , as expected. Whether this bifurcation is subcritical or supercritical depends on whether is negative or positive, respectively, indicating that bistability exists when κ > −Ω 2 . Note that an isolated community (i.e., at K = 0) is incoherent (synchronized) for κ < 0 (κ > 0), and thus the transition at K 3 is subcritical if (but not only if) there is a synchronized branch at K = 0. In Fig. 4(a) we plot the value of r obtained from Eq. (12) together with the locally synchronized and incoherent states, with solid and dashed curves indicating stable and unstable branches. We note that both stable and unstable branches describe globally synchronized states, i.e., states where the communities are phaselocked, but with the unstable branch displaying smaller degree of local synchronization. Conditions for the stability of the incoherent state can be determined using the Routh-Hurwitz criterion on the Jacobian associated with the linearization of Eqs. (4)-(6) (see Appendix C 1). When 0 ≤ κ ≪ 1, these conditions reduce to In Fig. 4(b) we plot these bifurcation curves in the stability diagram in (K, k) space for the symmetric case, using Ω = 2 and T = 1. Along with the bifurcations indicated by Eq. (13) we also indicate the lower bound (in K) for the globally synchronized state obtained from Eq. (12). In addition to incoherence, local synchronization, and global synchronization, we identify bistable regions between incoherence and global synchronization and between local and global synchronization where the globally synchronized branch folds over either of the other states.

B. Many communities
Lastly, we consider larger systems comprised of many (C > 2) communities that display even richer dynamics. For simplicity we let C be even and consider the case where Ω = Ω even = −Ω odd . (We note that C > 2 yields a system that is qualitatively distinct from the C = 2 case because of the inclusion of time-delayed interactions between even communities and time-delayed interactions between odd communities.) Even for the C = 4 case the resulting dynamics are more complicated than the C = 2 case (presented above), most notably due to a richer sequence of non-monotonic synchronization transitions. We illustrate this in Fig. 5(a), where we plot the local order parameters r = r σ versus coupling K for κ = −0.05, T = 1, and Ω = 2. Since κ < 0 the systems begins in the incoherent state when K = 0, but as K is increased the system undergoes a first bifurcation to local synchronization followed by a second bifurcation back to incoherence. This is then followed by a third subcritical bifurcation to global synchronization. Decreasing K highlights another hysteresis loop, however the transitions are again more complicated, first with a bifurcation from global synchronization to local synchronization (via explosive desynchronization) then a second bifurcation back to the incoherent state.
A linear stability analysis of the incoherent state (See Appendix C 2) yields the critical values where n = C/2. We note that the ± choice in Eq. (14) corresponds to two different branches of the same curve. The combination of these critical values describes the bifurcation involving the incoherent state, and allows us to sketch the stability diagram for the incoherent state in Fig. 5(b) illustrating the C = 4, 16, and 64 cases along with the limiting case C → ∞.
The lower-left regions of the diagram represent the respective regions of stability for the incoherent state and highlight the potential for non-monotonic transitions as K is increased for fixed k. In particular, for the C = 4 and 16 cases the upper portion of the bifurcation curve (given by K * ) decreases, then increases, giving rise to non-monotonic transitions similar to that observed in Fig. 5(a). This phenomenon does not persist for arbitrarily large numbers of communities, however, as can be seen for the C = 64 case, where the K † branch intersects the K * branch before the minimum is reached. In the C → ∞ limit the bifurcation comes solely from the K * branch (using the + sign), yielding K * ∞ = 4κ(T 2 Ω 2 +(T κ−1) 2 )/(T κ−1). Moreover, the symmetry Ω even = Ω odd allows us to calcu-late the globally synchronized branch analytically using similar techniques as those used in the two community case (see Appendix C 2), resulting in the implicit equation (defining K in terms of r) The globally synchronized branch plotted in Fig. 5(a) is given by Eq. (16), with solid and dashed curves indicating stable and unstable branches.

III. DISCUSSION
Despite the possible presence of both time delays and community structure in a number of real-world systems with synchronization properties, e.g., bacteria [10,11], power grids [12,13] and brain dynamics [15,16], the collective dynamics that emerge from their combination in heterogeneous systems has to date remained relatively unexplored. In this work we have demonstrated that the combination of these two important properties in coupled oscillator systems gives rise to a rich landscape of dynamical phenomena that does not arise from either of these property in isolation. Using both numerical simulations and analytical techniques we have shown that such systems often go through non-monotonic sequences of synchronization transitions, whereby increasing the coupling strength can first inhibit synchronization, and then promote it. In the two community case this manifests in a first bifurcation from local synchronization to incoherence, then a second (subcritical) bifurcation from incoherence to global synchronization. In the presence of more than two communities we demonstrated that this sequence is more complicated, with an initial bifurcation from incoherence to local synchronization, followed by the sequence described above. Moreover, when community-wise parameters are chosen to be asymmetric, the communities asymmetrically suppress one another's synchronization properties, leading one community to win as it forces the other to near incoherence. Interestingly, the roles of the two communities, in terms of which one remains synchronized while the other is pushed to near incoherence, unexpectedly reverse depending on the time delay. The novel phenomena observed in this work demonstrates the rich dynamics that can emerge from different dynamical and structural properties in oscillator systems.

Appendix A: Derivation of the Low Dimensional Equations
First we present the derivation of the low dimension equations, namely Eqs. (2) and (3), from the full high dimensional system given by Eq. (1). We begin by considering the continuum limit N σ → ∞ taken in such a way that the fraction of oscillator in each community remains constant. In this limit the macroscopic system state of each community σ can be described by the density function f σ (θ, ω, t) that describe the fraction of oscillators in each community σ with phase between θ and θ + dθ and frequency between ω and ω + dω at time t. We seek to describe the macroscopic dynamics as described by the local order parameters z σ = N −1 σ Nσ j=1 e iθ σ j . We also introduce the set of time delayed order parameters In the continuum limit the order parameters can be written and where we used the fact that the time delays τ σσ ′ ij in the definition of w σσ ′ i are chosen from the distribution h σσ ′ (τ ) ∝ s −τ /T σ ′ , which is independent of i and σ.
Next, due to the conservation of oscillators, the density functions each satisfy their own continuity equation Since the natural frequencies remain fixed and are drawn from the distributions g σ (ω σ ), the density functions can be written in a Fourier series of the form where c.c. stands for the complex conjugate of the previous term. Following the dimensionality reduction technique of Ott and Antonsen [25,26] we propose the ansatz that the Fourier coefficients decay geometrically, i.e., a n σ (ω σ , t)e inθ + c.c. .

(A5)
Remarkably, the dynamics defined by the infinite collection of Fourier coefficients collapse onto the low dimensional manifold describing the evolution of a via the differential equation where we have used K σσ = k and K σσ ′ = 2K/C (for σ = σ ′ ). Next, the function a σ can be tied directly to the local order parameter by inserting Eq. (A5) into Eq. (A1), yielding Assuming that the frequency distributions are Lorentzian, i.e., }, the integral can be evaluated using Cauchy's residue theorem, yielding z σ = a * σ (Ω σ − i∆, t). Evaluating Eq. (A6) at ω σ = Ω σ − i∆ and taking a complex conjugate yieldṡ Lastly, to close the dynamics we convert Eq. (A2) to a differential equation for w σ . Taking a Fourier transform, we have that and converting back to the time domain yields thereby closing the dynamics.

Appendix B: Numerical Validation of the Low Dimensional Equations
Next we present simulations validating the low dimensional equations. Since it is computationally infeasible for us to account for both a distribution of time delays and community structure, we use the hybrid approach proposed in Ref. [28]. In this approach, the microscopic dynamics of the oscillators is taken into account, but the effect of the time delay is accounted via the techniques presented above. We note that the treatment of distributed time delays via the Ott-Antonsen Ansatz has been validated before [17], so we focus on validating its interplay with the community structure. The hybrid approach consists in considering the equationṡ where z σ = r σ e iψσ , w σ = ρ σ e iφσ , and we have taken advantage of the technique used above to deal with the time delays. Note that the time delays enter in the evolution of oscillator i by the fact that the input it receives from community σ ′ , w σ ′ , is effectively a delayed version of the order parameter from community σ ′ . Both Eqs. (B1)-(B2) and Eqs. (1) have the same low-dimensional description, Eqs. (A8)-(A10). For a further discussion of the hybrid approach, see Ref. [28].
We begin by considering the asymmetric two community case with mean frequencies Ω 1 = −1 and Ω 2 = 2 and κ = 0.02. In Fig. 6 we plot the steady-state values of r 1 and r 2 in solid blue and dashed red curves, obtained in the same way to produce Fig. 1, with time delays T = 1 and 0.1 in panels (a) and (b). We then overlay the results from direct simulations of Eqs. (B1)-(B2) with N = 2 × 10 4 oscillators in each community in blue circles and red crosses. We note that the results from simulations are in excellent agreement.
Next we consider the symmetric case with two and four communities with (Ω 1 , Ω 2 ) = (−2, 2) and (Ω 1 , Ω 2 , Ω 3 , Ω 4 ) = (−2, 2, −2, 2), respectively, and T = 1. In Fig. 7(a) and (b) we plot the steady-state values of r 1 (which are indistinguishable from r 2 , r 3 , and r 4 ) in blue for the two and four community cases, respectively. Solid and dashed curves indicating stable and unstable branches, obtained in the same way to produce Fig. 3(a) and Fig. 4(a). We then overlay the results from direct simulations with N = 2 × 10 4 oscillators in each community in blue circles, again noting excellent agreement.

Appendix C: Linear Stability Analysis of the Incoherent State
Here we present linear stability analyses of the incoherent state for different mean frequency configurations. We restrict our attention to an even number of communities, σ = 1, 2 . . . , 2n.

Two communities with symmetric frequencies: exact analysis
First, we consider the case of two communities, n = 1, with symmetric mean frequencies, i.e., Ω 1 = −Ω 2 ≡ Ω. The Jacobian associated with perturbations from the incoherent solution z 1 = z 2 = w 1 = w 2 = 0 of Eqs. (4)- (6) is The characteristic equation is given by a 0 + a 1 λ + a 2 λ 2 + a 3 λ 3 + a 4 λ 4 = 0, with Applying the Routh-Hurwitz criterion [29], we obtain necessary and sufficient conditions for linear stability. The eigenvalues have all negative real part, corresponding to linear stability of the incoherent state, if the Routh-Hurwitz inequalities a i > 0, a 3 a 2 > a 4 a 1 , and a 1 a 2 a 3 > a 2 1 a 4 + a 0 a 2 3 are satis-fied. These inequalities simplify to T > 0, κT < 1, When T → 0 with κ > 0 the inequalities are not satisfied, showing that time delay is necessary to stabilize the incoherent state. For a given positive time delay, and sufficiently small κ, the inequalities reduce to K 1,2 < K < K 3 , with To leading order in κ, these inequalities reduce to Eq. (13), and agree with Eqs. (8)- (10) when Ω 1 = −Ω 2 .
Perturbations z + and z − are our main interest. We will show that the perturbations x j , y j decay, implying that the only relevant synchronization modes are those in which communities with equal frequencies behave identically. Inserting the new variables in Eqs. (C1)-(C3) we obtain the decoupled systemṡ Tẏ j = x j − y j , j = 2, . . . , n, (C10) First we focus on the perturbations that leave communities with the same mean frequencies synchronized, i.e., z ± and w ± . We determine the location of potential bifurcations by finding the values of K where the growth rate associated to Eqs. (C5)-(C8) is purely imaginary. We set z ± (t) = e iωtz ± , w ± (t) = e iωtw ± and obtain Eliminatingz ± andw ± leads to the complex equation for K and ω .
Separating the equation into real and imaginary parts and solv-ing the algebraic equations, we find the critical values To complete the analysis, we need to consider the evolution of perturbations transversal to the manifold where communities with the same frequency are synchronized, i.e., x j , y j . Setting x j (t) = e iαtx j , y j (t) = e iαtỹ j in Eqs. (C9)-(C12) and eliminatingx j ,ỹ j , we obtain the critical values From the exact form of the eigenvalues associated to Eqs. (C9)-(C12), one can check that their real part becomes negative for K > K T . For κ > 0, the incoherent state is already incoherent [the cusp in the curves in Fig. 4 For the region of interest, κ < 0, we obtain then K T < 0, which means that the only relevant modes of desynchronization are those in which communities with the same mean frequency remain synchronized, which yield the critical values in Eqs. (C13), (C14) plotted in Fig. 4(b).
In addition to the linear stability of the incoherent state, we present here the derivation of Equation (16) for the globally synchronized branch. We look for steady state solutions of (A8)-(A10) such that communities with the same frequency have the same order parameters, and set z σ = w σ = r for σ = 1, 2 . . . , n, z σ = w σ = re iα for σ = n + 1, . . . , 2n, where r and α are constants. This Ansatz gives Separating the equation into real and imaginary parts and solving for K gives Eq. (16) in the main text, Appendix D: Two communities with asymmetric frequencies: Asymmetric suppression In this section we consider the case of two communities with asymmetric frequencies, i.e., Ω 2 = −Ω 1 . A linear stability analysis of the incoherent state allows us to find the critical coupling constants and the modes of instability that give rise to the asymmetrically suppressed states.
First consider the case ω = Ω 1,2 +ω 1,2 κ+O(κ 2 ). Inserting this in (D8) we obtain which, when inserted in Eq. (D5) gives, to leading order in κ, Now consider the case ω = ω 3 κ + O(κ 2 ). We obtain and inserting ω = ω 3 κ in Eq. (D5) gives, to leading order, For K = 0, the communities are decoupled, and since there is no intra-community delay, each community synchronizes at k = 2 (corresponding to κ = 0). Therefore the incoherent state is unstable for small enough positive κ, and loses stability when K first increases beyond the smallest of K 1 , K 2 , and K 3 .

Asymmetric suppression
As discussed in the main text, for small κ we find that when K 1 < K < K 2 < K 3 community 2 synchronizes while community 1 remains almost completely incoherent. (An analogous situation occurs if K 2 < K < K 1 < K 3 .) This can be understood from the linear stability analysis by determining the form of the modes of instability at the bifurcations that occur at K = K 1 and K = K 2 . From Eqs. (D1)-(D4) we obtainz At K = K 1 , ω = Ω 1 + ω 1 κ, and we get to leading order in κ Similarly, at K = K 2 we get |z 1 | ∼ κ 1/2 |z 2 |. For the parameters used in Figure 1 (a) in the main text, we obtain |z 1 | ≈ 0.053|z 2 | at K = K 2 , and for those used in Figure  1 (b) we obtain |z 1 | ≈ 0.11|z 2 | at K = K 1 . For the case of Figure 1(a), where K 1 < K < K 2 < K 0 , the mode that desynchronizes as K is decreased past K 2 is therefore localized mostly in community 2. We find that this linear analysis still describes the behavior of the two communities in the interval K 1 < K < K 2 , where community 2 remains synchronized while community 1 has a very small degree of synchronization.