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

TL;DR

This paper investigates how community structure and distributed time delays in coupled oscillators lead to complex, non-monotonic synchronization transitions, including asymmetric suppression and competition among communities.

Contribution

It introduces a detailed analysis of the combined effects of community structure and time delays on synchronization, revealing non-monotonic and asymmetric transition phenomena.

Findings

Increasing coupling can both promote and inhibit synchronization.

Asymmetric suppression occurs where communities compete to dominate synchronization.

Analytical descriptions confirm numerical simulation results.

Abstract

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.

Equations

$${\dot{\theta }}_i^{\sigma }={\omega }_i^{\sigma }+\sum _{{\sigma }^{\prime }=1}^C\frac{K^{\sigma {\sigma }^{\prime }}}{N_{{\sigma }^{\prime }}}\sum _{j=1}^{N_{{\sigma }^{\prime }}}\sin [{\theta }_j^{{\sigma }^{\prime }}(t-{\tau }_{ij}^{\sigma {\sigma }^{\prime }})-{\theta }_i^{\sigma }(t)],$$

$${\dot{z}}_{\sigma }$$

$$+$$

$${\dot{z}}_1$$

$${\dot{z}}_2$$

$$T{\dot{w}}_1=z_1-w_1,T{\dot{w}}_2=z_2-w_2.$$

$${\left(\frac{K}{2}\right)}^2=-\frac{(1+{\omega }^2{T}^2{)}^2}{1-{\omega }^2{T}^2}[(\omega -{\Omega }_1)(\omega -{\Omega }_2)+{\kappa }^2],$$

$$\kappa (2\omega -{\Omega }_1-{\Omega }_2)=\frac{2\omega T[{\kappa }^2-(\omega -{\Omega }_1)(\omega -{\Omega }_2)]}{1-{\omega }^2{T}^2}.$$

$$K_1$$

$$K_2$$

$$K_3$$

$${\left(\frac{K}{2}\right)}^2=\frac{{\Omega }^2}{(1+r^2{)}^2}+{\left(\frac{r^2}{(1-r^2)}-\kappa \right)}^2.$$

$$2(1+{\Omega }^2{T}^2)\sqrt{\frac{\kappa }{T}}<K<2\Omega .$$

$$K^{\ast }$$

$$\displaystyle\pm\sqrt{n^{2}(T\kappa+1)^{2}-2n(T\kappa-1)^{2}+(T\kappa-1)^{2}}\bigg{]},$$

$$K^{†}$$

$$K$$

$$\displaystyle+n^{2}(1-r^{2})\sqrt{\frac{(\kappa-r^{2}(1+\kappa))^{2}}{n^{2}(1-r^{2})^{2}}+\frac{(2n-1)\Omega^{2}}{n^{4}(1+r^{2})^{2}}}\bigg{]}.$$

$$z_{\sigma }=\iint f_{\sigma }(\theta ,{\omega }^{\sigma },t)d\theta d{\omega }^{\sigma },$$

$$w_i^{\sigma {\sigma }^{\prime }}=\int z_{{\sigma }^{\prime }}(t-\tau )h_{\sigma {\sigma }^{\prime }}(\tau )d\tau \equiv w_{{\sigma }^{\prime }},$$

$$0$$

$$=\frac{\partial }{\partial t}{f}_{\sigma }+\frac{\partial }{\partial \theta }\left[f_{\sigma }\left({\omega }^{\sigma }+\frac{K^{\sigma \sigma }}{2i}\left(z_{\sigma }{e}^{-i{\theta }^{\sigma }}-z_{\sigma }^{\ast }{e}^{i{\theta }^{\sigma }}\right)+\sum _{{\sigma }^{\prime }\ne \sigma }^C\frac{K^{\sigma {\sigma }^{\prime }}}{2i}\left(w_{{\sigma }^{\prime }}{e}^{-i{\theta }^{\sigma }}-w_{{\sigma }^{\prime }}^{\ast }{e}^{i{\theta }^{\sigma }}\right)\right)\right].$$

$$f_{\sigma }(\theta ,{\omega }^{\sigma },t)=\frac{g_{\sigma }({\omega }^{\sigma })}{2\pi }\left[1+\sum _{n=1}^{\infty }{\hat{f}}_{\sigma ,n}({\omega }^{\sigma },t)e^{in\theta }+\text{c.c.}\right],$$

$$f_{\sigma }(\theta ,{\omega }^{\sigma },t)=\frac{g_{\sigma }({\omega }^{\sigma })}{2\pi }\left[1+\sum _{n=1}^{\infty }{a}_{\sigma }^n({\omega }^{\sigma },t)e^{in\theta }+\text{c.c.}\right].$$

$${\partial }_t{a}_{\sigma }=-i{\omega }^{\sigma }{a}_{\sigma }+\frac{k}{2}\left(z_{\sigma }^{\ast }-z_{\sigma }{a}_{\sigma }^2\right)+\frac{K}{C}\sum _{{\sigma }^{\prime }\ne \sigma }\left(w_{{\sigma }^{\prime }}^{\ast }-w_{{\sigma }^{\prime }}{a}_{\sigma }^2\right),$$

$$z_{\sigma }=\int g_{\sigma }({\omega }^{\sigma })a_{\sigma }^{\ast }({\omega }^{\sigma },t)d{\omega }^{\sigma }.$$

$${\dot{z}}_{\sigma }$$

$$(1+T_{\sigma }s){\overline{w}}_{\sigma }(s)={\overline{z}}_{\sigma }(s),$$

$$T_{\sigma }{\dot{w}}_{\sigma }$$

$${\dot{\theta }}_i^{\sigma }$$

$${\dot{w}}_{\sigma }$$

$$\displaystyle J=\left(\begin{array}[]{cccc}\kappa+i\Omega&0&0&\frac{K}{2}\\ 0&\kappa-i\Omega&\frac{K}{2}&0\\ \frac{1}{T}&0&-\frac{1}{T}&0\\ 0&\frac{1}{T}&0&-\frac{1}{T}\\ \end{array}\right)$$

$$a_0$$

$$a_1$$

$$a_2$$

$$a_3$$

$$a_4$$

$$T>0,$$

Full text

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 Strogatz2003 ; Pikovsky2003 . Examples from both natural and engineered systems include cardiac pacemaker dynamics Glass1988 , power grids Skardal2015SciAdv , and Josephson junctions Wiesenfeld1998PRE . The Kuramoto phase oscillator model and its variants have proven particularly useful in building an understanding of collective behavior Kuramoto1984 , and a large body of literature has identified features that give rise to rich nonlinear behavior, including external forcing Childs2008Chaos , multimodal frequency distributions Martens2009PRE , and mixed sign coupling Hong2011PRL .

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) Prindle2012Nature ; Scott2016ACS , interconnected regional or national power grid networks Nishikawa2015NJP ; Kim2015NJP , microgrids Lai2016IEEE , and synchronization of neuronal oscillations from different brain regions Deco2011Frontiers ; Sporns2016ARP . 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 Lee2009PRL ; Skardal2012PRE ; Yeung1999PRL ; Skardal2018IJBC ; Montbrio2004PRE ; Barreto2008PRE ; Pietras2016PRE , their combined effect on the collective dynamics remains, with a few exceptions (e.g., Ref. Wang2012Comm ), relatively unexplored. Here we address this gap and uncover a number of novel nonlinear behaviors.

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\geq 2$ communities of coupled phase oscillators governed by

[TABLE]

where $\theta_{i}^{\sigma}$ represents the phase of oscillator $i$ in community $\sigma$, $\omega_{i}^{\sigma}$ is its natural frequency, which is assumed to be drawn from the distribution $g_{\sigma}(\omega)$, $K^{\sigma\sigma^{\prime}}$ is the coupling strength between oscillators in communities $\sigma$ and $\sigma^{\prime}$, $\tau_{ij}^{\sigma\sigma}$ is the time delay between oscillators $i$ and $j$ in communities $\sigma$ and $\sigma^{\prime}$, which is assumed to be drawn from the distribution $h_{\sigma\sigma^{\prime}}(\tau)$, and $N_{\sigma}$ is the number of oscillators in community $\sigma$.

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^{\sigma\sigma}$ and $K=CK^{\sigma\sigma^{\prime}}/2$ (for $\sigma^{\prime}\neq\sigma$) denoting intra- and inter-community coupling. Next, we assume that within each community time delays are zero, i.e., $h_{\sigma\sigma}(\tau)=\delta(\tau)$, but between different communities the distribution $h_{\sigma\sigma^{\prime}}$ is exponential with mean $T_{\sigma^{\prime}}$, namely, $h_{\sigma\sigma^{\prime}}(\tau)\propto e^{-\tau/T_{\sigma^{\prime}}}$. We also consider the case where all communities are of the same size, i.e., $N_{\sigma}=N$ for all $\sigma$ and we assume that frequency distributions are Lorentzian with comunity-specific mean $\Omega_{\sigma}$ and width $\Delta$, i.e., $g_{\sigma}(\omega)=\Delta/\{\pi[\Delta^{2}+(\omega-\Omega_{\sigma})^{2}]\}$. Next, seeking a description for the local order parameters $z_{\sigma}=N^{-1}\sum_{j=1}^{N}e^{i\theta_{j}^{\sigma}}$ describing the degree of synchronization within each community, we apply the dimensionality reduction technique of Ott and Antonsen Ott2008Chaos ; Ott2009Chaos to obtain the following system of reduced equations

[TABLE]

where $w_{\sigma}$ represents a time-delayed order parameter for community $\sigma$, with $\sigma=1,2,\dots,C$ and the width parameter $\Delta$ of the natural frequency distributions have been scaled out. (Details of the dimensionality reduction are provided in Appendix A).

II.1 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 $\kappa=(k-2)/2$, noting that for $K=0$ each community undergoes a transition to synchronization at $k=2$. Therefore $\kappa>0$ ($\kappa<0$) indicates that the isolated communities would be synchronized (incoherent). Equations (2)-(3) become

[TABLE]

Equations (4)–(6) display some remarkable dynamics, which we illustrate in Fig. 2 for parameters $\kappa=0.02$, $\Omega_{1}=-1$, $\Omega_{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 [math] to $4$, then slowly decreasing it back to [math]. (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 $\psi_{1},\psi_{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}\approx 0$ or $r_{2}\approx 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 $\Omega_{\sigma}$ and have the same spread $\Delta=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” Gardenes2011PRL ) 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 local 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\times 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}=\tilde{z}_{1}e^{i\omega t}$, $z_{2}=\tilde{z}_{2}e^{i\omega t}$, $w_{1}=\tilde{w}_{1}e^{i\omega t}$, $w_{2}=\tilde{w}_{2}e^{i\omega t}$, we obtain the following equations for $\omega$ and $K$:

[TABLE]

In what follows we will focus on the case $0<\kappa\ll 1$. (The case of general $\kappa$ is treated exactly for the symmetric case $\Omega_{1}=-\Omega_{2}$ below). To balance Eq. (8) when $\kappa\ll 1$ there are three options, to first order in $\kappa$: $\omega=\Omega_{1}+\omega_{1}\kappa$, $\omega=\Omega_{2}+\omega_{2}\kappa$, and $\omega=\omega_{3}\kappa$, where $\omega_{i}$ are constants to be determined. Inserting these in Eqs. (7)-(8) we find the corresponding values of $K$, to leading order in $\kappa$:

[TABLE]

These values of $K$ are real and positive if the frequencies $\Omega_{1}$ and $\Omega_{2}$ have opposite sign. In general, depending on the values of $\Omega_{1}$, $\Omega_{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 $\kappa=0.02$, $\Omega_{1}=-1$, and $\Omega_{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 $K_{1}<K<K_{2},K_{3}$ or $K_{2}<K<K_{1},K_{3}$, and incoherence when $K_{1},K_{2}<K<K_{3}$. (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 $\text{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}\propto\kappa^{1/2}r_{1,2}$ (see Appendix D), so that for $\kappa\ll 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 delay $T_{c}=\sqrt{\sqrt{\Omega_{2}}-\sqrt{-\Omega_{1}}}/\sqrt{\sqrt{-\Omega_{1}}\Omega_{2}^{2}-\sqrt{\Omega_{2}}\Omega_{1}^{2}}$, at which point no asymmetric suppression state exists and the local synchronization state gives way directly to incoherence. Assuming that $|\Omega_{2}|>|\Omega_{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 symmetric case, $\Omega_{2}=-\Omega_{1}=\Omega$. Searching for phase-locked solutions of Eqs. (4)–(6) of the form $z_{1}=re^{i\psi}$, $z_{2}=re^{i\psi+\alpha}$, with $r$, $\psi$, $\alpha$ constants and $w_{1,2}=z_{1,2}$, gives the following implicit expression for the synchronized branch:

[TABLE]

By taking the limit $r\to 0^{+}$ we see that this branch begins at $K=K_{3}$, as expected. Whether this bifurcation is subcritical or supercritical depends on whether $K^{\prime\prime}(0)=-4\frac{\kappa+\Omega^{2}}{\sqrt{\kappa^{2}+\Omega^{2}}}$ is negative or positive, respectively, indicating that bistability exists when $\kappa>-\Omega^{2}$. Note that an isolated community (i.e., at $K=0$) is incoherent (synchronized) for $\kappa<0$ ($\kappa>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 phase-locked, 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\leq\kappa\ll 1$, these conditions reduce to

[TABLE]

In Fig. 4(b) we plot these bifurcation curves in the stability diagram in $(K,k)$ space for the symmetric case, using $\Omega=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.

II.2 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 $\Omega=\Omega_{\text{even}}=-\Omega_{\text{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_{\sigma}$ versus coupling $K$ for $\kappa=-0.05$, $T=1$, and $\Omega=2$. Since $\kappa<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

[TABLE]

where $n=C/2$. We note that the $\pm$ 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\to\infty$. 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^{\dagger}$ branch intersects the $K^{*}$ branch before the minimum is reached. In the $C\to\infty$ limit the bifurcation comes solely from the $K^{*}$ branch (using the $+$ sign), yielding $K_{\infty}^{*}=4\kappa(T^{2}\Omega^{2}+(T\kappa-1)^{2})/(T\kappa-1)$. Moreover, the symmetry $\Omega_{\text{even}}=\Omega_{\text{odd}}$ allows us to calculate 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$)

[TABLE]

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 Prindle2012Nature ; Scott2016ACS , power grids Nishikawa2015NJP ; Kim2015NJP and brain dynamics Deco2011Frontiers ; Sporns2016ARP , 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_{\sigma}\to\infty$ 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 $\sigma$ can be described by the density function $f_{\sigma}(\theta,\omega,t)$ that describe the fraction of oscillators in each community $\sigma$ with phase between $\theta$ and $\theta+d\theta$ and frequency between $\omega$ and $\omega+d\omega$ at time $t$. We seek to describe the macroscopic dynamics as described by the local order parameters $z_{\sigma}=N_{\sigma}^{-1}\sum_{j=1}^{N_{\sigma}}e^{i\theta_{j}^{\sigma}}$. We also introduce the set of time delayed order parameters $w_{i}^{\sigma\sigma^{\prime}}=N_{\sigma^{\prime}}^{-1}\sum_{j=1}^{N_{\sigma^{\prime}}}e^{i\theta^{\sigma^{\prime}}(t-\tau_{ij}^{\sigma\sigma^{\prime}})}$. In the continuum limit the order parameters can be written

[TABLE]

and

[TABLE]

where we used the fact that the time delays $\tau_{ij}^{\sigma\sigma^{\prime}}$ in the definition of $w_{i}^{\sigma\sigma^{\prime}}$ are chosen from the distribution $h_{\sigma\sigma^{\prime}}(\tau)\propto s^{-\tau/T_{\sigma^{\prime}}}$, which is independent of $i$ and $\sigma$.

Next, due to the conservation of oscillators, the density functions each satisfy their own continuity equation

[TABLE]

Since the natural frequencies remain fixed and are drawn from the distributions $g_{\sigma}(\omega^{\sigma})$, the density functions can be written in a Fourier series of the form

[TABLE]

where c.c. stands for the complex conjugate of the previous term. Following the dimensionality reduction technique of Ott and Antonsen Ott2008Chaos ; Ott2009Chaos we propose the ansatz that the Fourier coefficients decay geometrically, i.e.,

[TABLE]

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

[TABLE]

where we have used $K^{\sigma\sigma}=k$ and $K^{\sigma\sigma^{\prime}}=2K/C$ (for $\sigma\neq\sigma^{\prime}$).

Next, the function $a_{\sigma}$ can be tied directly to the local order parameter by inserting Eq. (21) into Eq. (17), yielding

[TABLE]

Assuming that the frequency distributions are Lorentzian, i.e., $g_{\sigma}(\omega^{\sigma})=\Delta/\{\pi[\Delta^{2}+(\omega^{\sigma}-\Omega_{\sigma})^{2}]\}$, the integral can be evaluated using Cauchy’s residue theorem, yielding $z_{\sigma}=a_{\sigma}^{*}(\Omega_{\sigma}-i\Delta,t)$. Evaluating Eq. (22) at $\omega^{\sigma}=\Omega_{\sigma}-i\Delta$ and taking a complex conjugate yields

[TABLE]

Lastly, to close the dynamics we convert Eq. (18) to a differential equation for $w_{\sigma}$. Taking a Fourier transform, we have that

[TABLE]

and converting back to the time domain yields

[TABLE]

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. Lee2011Chaos . 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 Lee2009PRL , so we focus on validating its interplay with the community structure. The hybrid approach consists in considering the equations

[TABLE]

where $z_{\sigma}=r_{\sigma}e^{i\psi_{\sigma}}$, $w_{\sigma}=\rho_{\sigma}e^{i\phi_{\sigma}}$, 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 $\sigma^{\prime}$, $w_{\sigma^{\prime}}$, is effectively a delayed version of the order parameter from community $\sigma^{\prime}$. Both Eqs. (27)-(28) and Eqs. (1) have the same low-dimensional description, Eqs. (24)-(26). For a further discussion of the hybrid approach, see Ref. Lee2011Chaos .

We begin by considering the asymmetric two community case with mean frequencies $\Omega_{1}=-1$ and $\Omega_{2}=2$ and $\kappa=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. (27)-(28) with $N=2\times 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 $(\Omega_{1},\Omega_{2})=(-2,2)$ and $(\Omega_{1},\Omega_{2},\Omega_{3},\Omega_{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\times 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, $\sigma=1,2\dots,2n.$

C.1 Two communities with symmetric frequencies: exact analysis

First, we consider the case of two communities, $n=1$, with symmetric mean frequencies, i.e., $\Omega_{1}=-\Omega_{2}\equiv\Omega$. The Jacobian associated with perturbations from the incoherent solution $z_{1}=z_{2}=w_{1}=w_{2}=0$ of Eqs. (4)-(6) is

[TABLE]

The characteristic equation is given by $a_{0}+a_{1}\lambda+a_{2}\lambda^{2}+a_{3}\lambda^{3}+a_{4}\lambda^{4}=0$, with

[TABLE]

Applying the Routh-Hurwitz criterion Hurwitz1895 , 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_{1}^{2}a_{4}+a_{0}a_{3}^{2}$ are satisfied. These inequalities simplify to

[TABLE]

When $T\to 0$ with $\kappa>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 $\kappa$, the inequalities reduce to $K_{1,2}<K<K_{3}$, with

[TABLE]

To leading order in $\kappa$, these inequalities reduce to Eq. (13),

[TABLE]

and agree with Eqs. (8)-(10) when $\Omega_{1}=-\Omega_{2}$.

C.2 Even communities with symmetric frequencies

Here we consider the case of an even number $C=2n$ of communities with mean frequencies $\Omega_{\sigma}=\Omega$, $\sigma=1,2,\dots,n$, $\Omega_{\sigma}=-\Omega$, $\sigma=n+1,\dots,2n$, which corresponds to the case shown in Fig. (4). Linearizing Eqs. (8)-(10) about the incoherent state we obtain

[TABLE]

Anticipating that communities with the same frequency will remain synchronized, we introduce the variables

[TABLE]

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. (34)-(36) we obtain the decoupled systems

[TABLE]

and

[TABLE]

First we focus on the perturbations that leave communities with the same mean frequencies synchronized, i.e., $z_{\pm}$ and $w_{\pm}$. We determine the location of potential bifurcations by finding the values of $K$ where the growth rate associated to Eqs. (38)-(41) is purely imaginary. We set $z_{\pm}(t)=e^{i\omega t}\tilde{z}_{\pm}$, $w_{\pm}(t)=e^{i\omega t}\tilde{w}_{\pm}$ and obtain

[TABLE]

Eliminating $\tilde{z}_{\pm}$ and $\tilde{w}_{\pm}$ leads to the complex equation for $K$ and $\omega$

[TABLE]

Separating the equation into real and imaginary parts and solving the algebraic equations, we find the critical values

[TABLE]

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\alpha t}\tilde{x}_{j}$, $y_{j}(t)=e^{i\alpha t}\tilde{y}_{j}$ in Eqs. (42)-(45) and eliminating $\tilde{x}_{j}$, $\tilde{y}_{j}$, we obtain the critical values

[TABLE]

From the exact form of the eigenvalues associated to Eqs. (42)-(45), one can check that their real part becomes negative for $K>K_{T}$. For $\kappa>0$, the incoherent state is already incoherent [the cusp in the curves in Fig. 4(b) occurs at $\kappa=(\sqrt{2n-1}-n)/(n+\sqrt{2n-1})<0$]. For the region of interest, $\kappa<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. (46), (47) 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 (24)-(26) such that communities with the same frequency have the same order parameters, and set $z_{\sigma}=w_{\sigma}=r$ for $\sigma=1,2\dots,n$, $z_{\sigma}=w_{\sigma}=re^{i\alpha}$ for $\sigma=n+1,\dots,2n$, where $r$ and $\alpha$ are constants. This Ansatz gives

[TABLE]

Separating the equation into real and imaginary parts and solving for $K$ gives Eq. (16) in the main text,

[TABLE]

Appendix D Two communities with asymmetric frequencies: Asymmetric suppression

In this section we consider the case of two communities with asymmetric frequencies, i.e., $\Omega_{2}\neq-\Omega_{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.

D.1 Critical coupling constants

Linearizing Eqs. (24)-(26) about the incoherent state, and again defining $k/2=1+\kappa$ and setting $\Delta=1$ we obtain

[TABLE]

We find the location of bifurcations by looking for purely imaginary growth rates. To this end, we set $\delta z_{1,2}=e^{i\omega t}\tilde{z}_{1,2}$, $\delta w_{1,2}=e^{i\omega t}\tilde{w}_{1,2}$, and obtain

[TABLE]

Eliminating $\tilde{z}_{1,2}$ and $\tilde{w}_{1,2}$ we obtain a complex equation for $K$ and $\omega$, with its real and imaginary parts giving

[TABLE]

Eliminating $K$ from these equations we obtain a cubic equation for $\omega$,

[TABLE]

In what follows we will consider the limit $\kappa\ll 1$ and argue that the three solutions to the cubic are, to leading order, $\omega\approx\Omega_{1}$, $\omega\approx\Omega_{2}$, and $\omega\approx\kappa(\Omega_{1}+\Omega_{2})/(2T\Omega_{1}\Omega_{2})$. When $\kappa\ll 1$, the leading order terms in (56) are

[TABLE]

Since the left hand side is of order $\kappa$, the right hand side needs to be of order $\kappa$, which can be accomplished by the three choices $\omega=\Omega_{1}+\omega_{1}\kappa+\mathcal{O}(\kappa^{2})$, $\omega=\Omega_{2}+\omega_{2}\kappa+\mathcal{O}(\kappa^{2})$, or $\omega=\omega_{3}\kappa+\mathcal{O}(\kappa^{2})$, where $\omega_{i}$ are to be determined.

First consider the case $\omega=\Omega_{1,2}+\omega_{1,2}\kappa+\mathcal{O}(\kappa^{2})$. Inserting this in (57) we obtain

[TABLE]

which, when inserted in Eq. (54) gives, to leading order in $\kappa$,

[TABLE]

Now consider the case $\omega=\omega_{3}\kappa+\mathcal{O}(\kappa^{2})$. We obtain

[TABLE]

and inserting $\omega=\omega_{3}\kappa$ in Eq. (54) gives, to leading order,

[TABLE]

For $K=0$, the communities are decoupled, and since there is no intra-community delay, each community synchronizes at $k=2$ (corresponding to $\kappa=0$). Therefore the incoherent state is unstable for small enough positive $\kappa$, and loses stability when $K$ first increases beyond the smallest of $K_{1}$, $K_{2}$, and $K_{3}$.

D.2 Asymmetric suppression

As discussed in the main text, for small $\kappa$ 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. (50)-(53) we obtain

[TABLE]

At $K=K_{1}$, $\omega=\Omega_{1}+\omega_{1}\kappa$, and we get to leading order in $\kappa$

[TABLE]

Similarly, at $K=K_{2}$ we get $|\tilde{z}_{1}|\sim\kappa^{1/2}|\tilde{z}_{2}|$. For the parameters used in Figure 1 (a) in the main text, we obtain $|\tilde{z}_{1}|\approx 0.053|\tilde{z}_{2}|$ at $K=K_{2}$, and for those used in Figure 1 (b) we obtain $|\tilde{z}_{1}|\approx 0.11|\tilde{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.