Abrupt Desynchronization and Extensive Multistability in Globally Coupled Oscillator Simplices
aa r X i v : . [ n li n . AO ] M a y Abrupt Desynchronization and Extensive Multistability in Globally Coupled Oscillator Simplices
Per Sebastian Skardal ∗ and Alex Arenas Department of Mathematics, Trinity College, Hartford, CT 06106, USA Departament d’Enginyeria Informatica i Matem´atiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain
Collective behavior in large ensembles of dynamical units with non-pairwise interactions may play an im-portant role in several systems ranging from brain function to social networks. Despite recent work pointingto simplicial structure, i.e., higher-order interactions between three or more units at a time, their dynamicalcharacteristics remain poorly understood. Here we present an analysis of the collective dynamics of such asimplicial system, namely coupled phase oscillators with three-way interactions. The simplicial structure givesrise to a number of novel phenomena, most notably a continuum of abrupt desynchronization transitions withno abrupt synchronization transition counterpart, as well as extensive multistability whereby infinitely manystable partially synchronized states exist. Our analysis sheds light on the complexity that can arise in physicalsystems with simplicial interactions like the human brain and the role that simplicial interactions play in storinginformation.
PACS numbers: 05.45.Xt, 89.75.Hc
Research into the macroscopic dynamics of large ensem-bles of coupled oscillators have extended our understandingof natural and engineered systems ranging from cell cycles topower grids [1–5]. However, with few exceptions (including[6–8]), little attention has been paid to the synchronizationdynamics of coupled oscillator systems where interactions arenot pair-wise, but rather n -way, with n ≥ . Such interac-tions are called “simplicial”, where an n -simplex representsan interaction between n + 1 units, so -simplices describethree-way interactions, etc [9]. Recent advances suggest thatsimplicial interactions may be vital in general oscillator sys-tems [10–12] and may play an important role in brain dynam-ics [13–15] and other complex systems phenomena such as,the dynamics of collaborations [16] or social contagion [17].In particular, interactions in 2-simplices (named holes or cav-ities) are important because they can describe correlations inneuronal spiking activity (that can be mapped to phase oscil-lators [18]) in the brain [19] providing a missing link betweenstructure and function. In fact, coupled oscillator systems thatdisplay clustering and multi-branch entrainment have beenshown to be useful models for memory and information stor-age [20–26]. Despite these findings, the general collectivedynamics of coupled oscillator simplices and their utility instoring information are poorly understood.In this work we study large coupled oscillator simplicialcomplexes, considering the impact of -simplices, i.e., three-way interactions, on collective behavior. Specifically, we con-sider the 2- and 1-simplex multilayer system given by ˙ θ i = ω i + KN N X j =1 N X k =1 sin ( θ j + θ k − θ i ) , (1) φ i = ν i + κN N X j =1 sin( φ j − φ i ) + d sin( θ i − φ i ) . (2)The dynamics in the θ -layer are the natural generalizationof the classical Kuramoto model [27] with -simplex interac-tions (namely, coupling is sinusoidal and diffusive), where θ i represents the phase of oscillator i with i = 1 , . . . , N , ω i is itsnatural frequency which is assumed to be drawn from the dis-tribution g ( ω ) , and K is the global coupling strength. The θ -oscillators then drive a population of φ -oscillators in anotherlayer, for which a one-to-one correspondence exists follow-ing a multiplex structure [28]. Specifically, the φ -oscillatorsevolve subject to the classical Kuramoto model (with natu-ral frequencies ν i and coupling strength κ ), but with an ad-ditional driving term with strength d . While numerical in-vestigations of systems with non-pairwise interactions haveuncovered multistability and chaos [6–8], few analytical re-sults exist and their collective behavior remains largely unex-plored. Here we focus on large systems and obtain an ana-lytical description of the macroscopic dynamics using a par-tial dimensionality reduction obtained via a variation of theOtt-Antonsen ansatz [29, 30]. In particular, the 2-simplexmacroscopic dynamics are captured by a combination of twoorder parameters that capture the degree of synchronizationand asymmetry as oscillators organize into two distinct syn-chronized clusters. Here, clustering refers to multiple phase-locked groups [20, 25] at different locations on the circlerather than distinct groups of identically synchronized oscil-lators [31]. We uncover a novel phenomenon where a con-tinuum of abrupt desynchronization transitions emerge, eachat a different critical coupling strength, depending on theasymmetry of the system. Interestingly, no complementaryabrupt synchronization transitions occur [32, 33]. This con-tinuum stems from an extensive multistability whereby, forsufficiently strong coupling, an infinite number of distinct par-tially synchronized states are stable in addition to the incoher-ent state, which is stable for all finite coupling strengths. Thismultistability indicates the capability of storing a wide arrayof possible information as different oscillator arrangements.Serving as a minimal model for memory storage, the systemcaptures the critical properties of easily transitioning from aninformation storage state (i.e., synchronized) to the restingstate (i.e., incoherent) [34] via abrupt desynchronization. Thesystem then may return to another information storage statewith an appropriately chosen perturbation. Moreover, the 1-simplex layer can similarly store information provided thatthe driving strength from the 2-simplex layer is sufficientlystrong. The rich nonlinear dynamics that emerge in this rel-atively simple extension of pair-wise coupling to three-waycoupling highlights the complexity that may arise via simpli-cial interactions in systems like the human brain and the im-plications of these behaviors on information storage.Noting that the dynamics of the 1-simplex layer in Eq. (2)are driven by the dynamics of the 2-simplex layer in Eq. (1),we focus first on the θ -layer our analysis by introducing thegeneralized order parameters z q = N P Nj =1 e qiθ j for q = 1 and . Note that z is the classical Kuramoto order parameterwhile z typically measures clustering, which we will see inthis system. Using the polar decompositions z q = r q e iψ q werewrite Eq. (1) as ˙ θ i = ω i + Kr sin[2( ψ − θ i )] . (3)We then consider the continuum limit N → ∞ where thestate of the system can be described by a density function f ( θ, ω, t ) which describes the density of oscillator with phasebetween θ and θ + δθ and natural frequency between ω and ω + δω at time t . Because the number of oscillators in thesystem is conserved f must satisfy the continuity equation ∂ t f + ∂ θ ( f ˙ θ ) . Moreover, because each oscillator’snatural frequency is fixed and drawn from g ( ω ) the densityfunction f ( θ, ω, t ) may be expanded into a Fourier series ofthe form f ( θ, ω, t ) = g ( ω )2 π h P ∞ n =1 ˆ f n ( ω, t ) e inθ + c.c. i ,where ˆ f n ( ω, t ) is the n th Fourier coefficient and c.c. repre-sents the complex conjugate of the previous sum.We then consider the symmetric and asymmetric parts f s ( θ, ω, t ) and f a ( θ, ω, t ) , respectively, of f ( θ, ω, t ) whichsatisfy f ( θ, ω, t ) = f s ( θ, ω, t ) + f a ( θ, ω, t ) with symme-tries f s ( θ + π, ω, t ) = f s ( θ, ω, t ) and f a ( θ + π, ω, t ) = − f a ( θ, ω, t ) . Note that the linearity of the continuity equa-tion implies that if both f s and f a are solutions, thenso is f . While the asymmetric part f a does not al-low for dimensionality reduction, the symmetric part f s does. Noting that the Fourier series of f s is given by theeven terms of the Fourier series of f , i.e., f s ( θ, ω, t ) = g ( ω )2 π h P ∞ m =1 ˆ f m ( ω, t ) e imθ + c.c. i , we make the ansatzthat each even Fourier coefficient decays geometrically, i.e., ˆ f m ( ω, t ) = a m ( ω, t ) [29, 30]. Inserting this and Eq. (2) intothe continuity equation, we find that each subspace spannedby even terms e imθ collapse onto the same low-dimensionalmanifold characterized by the condition ∂ t a = − iωa + K (cid:0) z ∗ − z a (cid:1) . (4)Equation (4) describes the evolution of the complex func-tion a ( ω, t ) , and thereby f s , and depends on the order parame-ter z . Moreover, a ( ω, t ) can be linked to the order parameter z as follows. First, note that in the limit N → ∞ we have that z = RR f s ( θ, ω, t ) e iθ dθdω , and after inserting the Fourier series for f s this reduces to z = R g ( ω ) a ∗ ( ω, t ) dω . To fur-ther simplify the relationship we make the assumption that thefrequency distribution g ( ω ) is Lorentzian with mean ω andspread ∆ , i.e., g ( ω ) = ∆ / { π (cid:2) ( ω − ω ) + ∆ (cid:3) } , which hastwo simple poles in the complex plane at ω = ω ± i ∆ . Theintegral can then be evaluated using Cauchy’s Residue Theo-rem [35] by closing the integral contour with a semicircle ofinfinite radius in the lower-half plane and evaluating at the en-closed pole, yielding z = a ∗ ( ω − i ∆ , t ) . We then evaluateEq. (4) at ω = ω − i ∆ to obtain ˙ z = 2 iω z − z + K (cid:0) z − z ∗ z (cid:1) . (5)Next we introduce the rescaled parameters ˜ K = K/ ∆ and ˜ ω = ω / ∆ with rescaled time ˜ t = ∆ t , then use polar decom-positions, yielding (where we have dropped the ∼ -notation forconvenience) ˙ r = − r + Kr (1 − r ) cos(2 ψ − ψ ) , (6) ˙ ψ = 2 ω + Kr r r sin(2 ψ − ψ ) . (7)Equations (6) and (7) describe the dynamics of the even part f s , which falls onto a low dimensional manifold similar tothe Ott-Antonsen manifold and describes the evolution of z .However, these equations do not capture the asymmetric partof the dynamics, and moreover they depends on the asymmet-ric part via z ’s dependence on z . As we will see, this reflectsthe system’s dependence on asymmetry in oscillator arrange-ments between two clusters.To close the dynamics we apply a self-consistency anal-ysis to characterize the order parameter z . We first notethat, by entering the rotating frame θ θ + ω t we set ω = 0 so that ˙ ψ = ˙ ψ = 0 . Moreover, by rotating ini-tial conditions θ (0) θ (0) + ψ (0) we set ψ = ψ = 0 .Equation (2) then implies that oscillators that become phase-locked satisfy | ω i | ≤ Kr , in which case they relax to oneof the two stable fixed points θ i = θ ∗ ( ω i ) or θ ∗ ( ω i ) + π ,where θ ∗ ( ω ) = arcsin( ω/Kr ) / . These two fixed pointscorrespond to the two clusters that the phase-locked oscilla-tors organize into. Specifically, phase-locked oscillators start-ing near θ = 0 or π will end up at the fixed points θ ∗ ( ω ) or θ ∗ ( ω ) + π , respectively. The phase-locked population isdescribed by the density function f locked ( θ, ω ) = ηδ ( θ − θ ∗ ( ω )) + (1 − η ) δ ( θ − θ ∗ ( ω ) − π ) , (8)where the asymmetry parameter η describes the fraction ofphase-locked oscillators in the θ = 0 cluster. On the otherhand, oscillators satisfying | ω i | > Kr drift for all time andrelax to the stationary distribution f drift ( θ, ω ) = p ω − K r π [ ω + Kr sin(2 ψ − θ )] . (9)Next, the order parameter z is given by the integral z = coupling, K sy n c , r coupling, K sy n c , r (a) (b) = 0.8 = 1 decreasing K decreasing K = 1 = 0.8increasing K increasing K FIG. 1.
Synchronization profiles.
The order parameters (a) r and(b) r as a function of the coupling strength K for various asymme-try values η . Blue, red, green, orange, and purple circles representresults obtained from direct simulations of Eq. (2) with N = 10 oscillators for η = 1 , . , . , . , and . , respectively. RR f ( θ, ω, t ) e iθ dθdω , which after inserting the density f asdefined by Eqs. (8) and (9) reduces to r = (2 η − Z Kr − Kr s p − ( ω/Kr ) g ( ω ) dω, (10)where the contribution from the drifting oscillators vanishesdue to the symmetry of f drift . Returning to r , Eq. (6) impliesthat at steady state we have r = − p K r Kr . (11)Thus, the macroscopic steady-state is described by Eqs. (10)and (11).Interpreting these analytical results in the context of nu-merical simulations allows us to understand novel phenomenathat occur in the dynamics of Eq. (1). Beginning with simu-lations of a system with N = 10 oscillators whose naturalfrequencies are Lorentzian with ω = 0 and ∆ = 1 , we con-sider initial conditions of varying asymmetry, setting initialphases to θ i (0) = 0 and π with probabilities η and − η , re-spectively. We then begin simulations at K = 16 and afterreaching steady-state slowly decrease K to zero, then restoreit slowly to . In Figs. 1(a) and (b) we plot the resultingvalues for the order parameters r and r , respectively, for η = 1 (blue, top), . (red), . (green), . (orange), and . (purple, bottom). Arrows indicate the direction of K andresults corresponding to decreasing and increasing K are plot-ted as open and filled circles, respectively. As K is initiallydecreased solutions traverse different partially synchronizedstates that are determined by η , until each branch undergoesa discontinuous jump to the incoherent state at different criti-cal coupling strengths in abrupt desynchronization transitions.Importantly, this highlights both a rich extensive multistabil-ity (here five different branches are shown, but in the thermo-dynamic limit an infinite number of such states exist) and acontinuum of abrupt desynchronization transitions at differ- coupling, K sy n c , r coupling, K sy n c , r (a) = 1 = 1 (b) r r = 0.8 = 0.8 FIG. 2.
Synchronization branches.
For order parameters (a) r and(b) r , the synchronization branches predicted by Eqs. (10) and (11)plotted as solid curves for asymmetries η = 1 , . , . , . , and . as well as the minimum branches [predictions for which are givenin Eq. (15)]. Results obtained from direct simulations with N = 10 oscillators are plotted in blue, red, green, orange, and purple circlesfor each value of η and black circles represent minimum branches. ent coupling strengths. Partially synchronized branches arecharacterized by the asymmetry parameter η , indicating thatthese complex dynamics arise from different allocations ofphase-locked oscillators in the two clusters at θ = 0 and π .Next, as K is restored to its initial value of no spontaneoustransitions back to synchronization occur, indicating no abruptsynchronization to complement the abrupt desynchronizationtransitions. In Figs. 2(a) and (b) we show that our theory cap-tures these dynamics, plotting in solid curves the theoreticalpredictions of r and r , respectively, given by Eqs. (10) and(11). We use the same values of η as in Fig. 1 and overlaythe results from simulations in open circles, noting excellentagreement.To better understand the nature of the abrupt desynchro-nization and multistability phenomena, we investigate theminimal partially synchronized states, i.e., minimal non-zerovalues r min and r min allowable for different coupling strengths.We first plot numerical results for r min and r min in Figs. 2(a)and (b) using black circles. Interestingly, while r min appearsto decay monotonically with K , r min remains roughly con-stant. In fact, the minimal branch r min is constant and canbe leveraged to analytically describe the branches r min ( K ) and r min ( K ) as well as the corresponding asymmetry value η min ( K ) and critical coupling strengths K c ( η ) where theabrupt desynchronization occurs. We proceed by invertingEq. (11), obtaining Kr = 2 r / (1 − r ) , which can be in-serted into Eq. (10), yielding s r − r = √ K (2 η − × Z r / (1 − r ) − r / (1 − r ) vuut q − [ ω (1 − r ) / r ] g ( ω ) dω. (12)While Eq. (12) appears more complicated than Eq. (10), we asymmetry, ab r up t de sy n c ., K c theorysimulation phases, d i s t r i bu t i on , f () K = 6.7K = 6.6K = 6.5K = 6.4 (a) (b)
FIG. 3.
Abrupt desynchronization transition. (a) The critical cou-pling strength K c at which the abrupt desynchronization transitionoccurs as a function of the asymmetry η . The solid curve representsthe theoretical prediction given by Eq. (15) and black circles repre-sent observations from direct simulation with N = 10 oscillators.(b) Illustration of the distribution f ( θ ) of phases as the system passesthough the abrupt desynchronization transition. Here η = 0 . with K c ≈ . . Distributions for K = 6 . , . , . , and . are plot-ted in solid blue, dot-dashed red, dotted green, and dashed orangecurves. note that the coupling strength K has been entirely scaled outof the integral, appearing outside with (2 η − . We thereforeconclude that if the quantities √ K and η − cancel oneanother out, i.e., √ K (2 η − is constant, it follows that thesolution r in Eq. (12) is independent of K . We thereforepropose the ansatz √ K (2 η −
1) = const. and use the initialcondition η min ( K c (1)) = 1 , where K c (1) denotes the veryfirst coupling strength where a synchronized state is possiblewith η = 1 , yielding η min ( K ) = p K c (1)2 √ K + 12 . (13)Equation. (13) implies that along the minimum branch wehave that √ K (2 η −
1) = p K c (1) ≈ . , which can beused in Eq. (12) to compute the minimum branch of r , and inturn r via Eq. (11), yielding r min ( K ) ≈ . √ K and r min ( K ) ≈ . . (14)In Figs. 2(a) and (b) we plot these minimum branches in solidblack curves, which agree with the simulation results. Lastly,by inverting Eq. (13) we find the critical coupling strength K c as a function of η where the abrupt desynchronization transi-tion occurs, namely K c ( η ) ≈ . η − . (15)In Fig. 3(a) we plot the theoretical prediction of the abruptdesynchronization point K c ( η ) as a solid curve vs observa-tions from direct simulations as black circles, noting excellentagreement.The dynamics in the 2-simplex layer serves as a minimalmodel for memory and information storage, capturing a num- -1 coupling, d sy n c , = 0.8 = 1d = FIG. 4.
The order parameter ρ de-scribing synchronization of φ -oscillators in the 1-simplex layer as afunction of the driving strength d . Connected circles represent nu-merical simulations with oscillators. 2-simplex and 1-simplexcoupling are k = 12 and κ = 3 , and asymmetry values η = 1 , . , . , . , and . are colored blue, red, green, orange, and purple,respectively. dashed horizontal lines denote the values of r for each η and the dotted vertical line indicated d = κ . ber of critical properties. Each distinct synchronized state cor-responds to a specific piece of information, differentiated bythe clustering arrangement of the oscillators. Moreover, syn-chronized state the system can quickly and easily transitionto the resting state described by incoherence via the abruptdesynchronized transition. The microscopic properties of thisabrupt desynchronization transition are illustrated in Fig. 3(b),where for η = 0 . the distribution f ( θ ) of phases is plotted as K is decreased slowly through the critical value of K c ≈ . .Before the transition the distribution is asymmetrically clus-tered about θ = 0 and π and changes slowly until at K = K c all information is lost as the distribution becomes uniform,representing a resting state.Lastly, we shift our focus to the φ -layer dynamics, rewrit-ing Eq. (2) as ˙ φ i = ν i + κρ sin( ϕ − θ i ) + d sin( θ i − φ ) ,where ρ e iϕ = N P Nj =1 e iφ j . The dynamics of each φ i canbe written in terms of a single sine term, however, due to thefact that all θ i ’s are generically different, the target phase andamplitude are all distinct, making self-consistency analysesdifficult. However, the microscopic properties may still beunderstood by comparing κ and d . Namely, for κ ≫ d the dy-namics are effectively undriven, in which case oscillators thatbecome entrained form a single cluster. On the other hand,for d ≫ κ the driving term dominates, causing each φ i toentrain with θ i , thereby mirroring any cluster behavior. Thisis illustrated in Fig. 4, where we plot ρ as a function of d given κ = 3 and K = 12 for asymmetries η = 1 , . , . , . , and . in blue, red, green, orange, and purple circles,respectively. Dashed horizontal lines denote the r value foreach η , illustrating that ρ → r in the d ≫ κ regime and infact the 1-simplex layer can also store information. However,as d decreases into the d ≪ κ regime this property is lost asthe 1-simplex dynamics no longer reflect the structure in the2-simplex layer.The analysis presented here demonstrates how the exten-sion from pair-wise to more general simplicial (specifically, -simplex) interactions in coupled oscillator systems can giverise to a host of complex nonlinear phenomena including in-formation and memory storage. Moreover, these phenom-ena can be captured and described using analytical methods.In particular, we have characterized a continuum of abruptdesynchronization transitions that occur at different couplingstrengths without any abrupt synchronization transitions. Thiscontinuum stems from extensive multistability, whereby forsufficiently strong coupling an infinite number of partiallysynchronized states are stable. These different stable statesrepresent synchronized states organized via different asym-metries into two clusters of entrained oscillators. In additionto highlighting the possible complexity that may arise in cou-pled oscillator systems with simplicial interactions, we hy-pothesize that simplicial interactions may give rise to novelnonlinear phenomena in other complex systems. ∗ [email protected][1] S. H. Strogatz, Sync: the Emerging Science of Spontaneous Or-der (Hypernion, 2003).[2] A. Pikovsky, M. Rosenblum, and J. Kurths,