Avalanches impede synchronization of jammed oscillators
SSelf-organized avalanches in globally-coupled phase oscillators
William Gilpin ∗ Department of Applied Physics, Stanford University (Dated: June 14, 2019)Spontaneous desynchronization of coupled oscillators occurs in diverse systems spanning fromthe mammalian brain to alternating current power grids. Here, we model desynchronization us-ing a generalization of the classical Kuramoto model, which includes short-range phase repulsionamong individual oscillators. Surprisingly, we find that our model exhibits self-organized avalanchesat intermediate values of the repulsion strength, and that these avalanches have similar statisticalproperties and scaling exponents to cascades seen in real-world oscillator ensembles—including neu-ronal avalanches. We show that avalanches in our system arise due to a critical mechanism based oncompetition between the mean field and local interactions, which can be recreated using a cellularautomaton based on classical traffic models. We exactly solve our system in the many-oscillatorlimit, and find that criticality arises due to the system’s strong sensitivity to small-scale noise, herearising from local rearrangements among oscillators.
Globally-coupled phase oscillators represent a basic,tractable paradigm for understanding synchronization ina variety of systems—spanning from coordination amongneurons, to load-balancing on power grids, to interac-tions on social networks [1–3]. However, recent yearshave revealed a variety of exotic phenomena occurringeven in simple systems of phase oscillators, such chimerastates [4, 5], glass-like relaxation mediated by a ”volcano”transition [6], and oscillation death via broken rotationalsymmetry [7].Many such discoveries are framed in terms of theireffects on the synchronous state, which represents aglobally-stable solution of the underlying dynamicalequations. However, many real-world oscillator popula-tions, such as neuronal ensembles, exhibit a maximallydesynchronous state, in which the phases of individualoscillators become maximally spaced apart on the unitcircle [8, 9]. Such dynamics are achievable, in principle,by introducing negative couplings into standard phase os-cillator models [10, 11]; alternative approaches include in-troducing phase offsets or time delays in the interactionsamong oscillators, or introducing specific pairwise cou-plings among oscillators that embed them on a complexgraph [2, 12, 13]. The tunability of synchronization insuch models introduces the broader question of whetherunique phenomena can occur at the critical interme-diate point between the synchronous and maximally-desynchronous states.Here, we introduce a minimal generalization of clas-sical oscillator models, which allows a smooth transitionbetween full synchrony and maximal desynchrony. At in-termediate states, we observe a surprising phenomenonconsisting of avalanches, in which the system seeminglyapproaches synchrony, only to abruptly desynchronizeat quasi-random intervals. This phenomenon mirrorsavalanches seen in real-world oscillator networks, includ-ing neuronal ensembles, financial markets, and power ∗ [email protected] grids [14–16]. We trace these avalanches to a criticalmechanism driven by competition between large scale at-traction of oscillators to the mean field, and small-scalerearrangements, which we recreate using a cellular au-tomaton model inspired by self-organized criticality intraffic. Finally, we solve our model exactly in the many-oscillator limit, and show how the existence of avalanchesin the system arises due to the amplication of ”noise”provided by rearrangements.Our model consists of the classical Kuramoto modelof globally and attractively-coupled phase oscillators; wemodify this model by adding an additional term corre-spond to short-range repulsion among oscillators withproximal phases,˙ θ j = ω + 1 N N (cid:88) k =1 (cid:18) K sin( θ k − θ j ) − ddθ j V ( θ j , θ k ) (cid:19) (1)where the interaction potential introduces short-rangerepulsion. Here, due to its finite range and analytictractability, we choose a triangular step V ( θ j , θ k ) = V (1 − H ( | D ( θ j , θ k ) | − d ))(1 − (1 /d ) | D ( θ j , θ k ) | ) where d = 2 π/N and H denotes the unit step function. Theangular distance function D ( θ j , θ k ) corresponds to thesigned distance between two angles along the unit circle.When V = 0, Eq. 1 becomes the standard Kuramotomodel with global coupling and no phase offsets; thismodel exhibits a well-known and solvable phase transi-tion to synchrony as K increases [3]. This synchroniza-tion transition can be readily observed via the dynamicsof the first-order Daido order parameter R ( t ), where R ( t ) e i Φ( t ) ≡ N (cid:88) k =1 e iθ k ( t ) , which reaches a stable steady-state R = 1 when K issufficiently high. However, in the case of non-zero repul-sion V >
0, Eq. 1 admits another stable solution when K = 0, corresponding to equal spacing of the oscillatorsalong the unit circle, R = 0. Thus, the additional term a r X i v : . [ n li n . AO ] J un in Eq. 1 serves to stabilize the maximally desynchronizedstate via steric interaction.Unexpectedly, at intermediate values of K and V ,the order parameter undergoes irregular oscillations, dis-playing epochs of gradual synchronization punctuated byrapid desynchronization (Figure 1A). During synchro-nization, R exhibits consistent exponential growth at afixed rate; however, at varying intervals this growth halts,and the order parameter rapidly decreases to a value ap-proaching the maximally desynchronous state R = 0. In-spection of individual oscillators indicates that the dy-namics during synchronization largely match those ofthe classical Kuramoto model at the same value of K ;however, the desynchronization period occurs much morerapidly, and with a wide range of amplitudes. Becauseour model contains no explicit noise source, we hypothe-size that the onset of desynchronization arises from fine-scale details of the packing of oscillators along the unitcircle. All variation therefore originates from the ran-dom initial placement of the oscillators on the interval[0 , π ]—small changes in the initial conditions thereforehave a pronounced effect on the timing of desynchroniza-tion events at later times.These irregular oscillations are particularly strikinggiven the homogeneity of the oscillators, and the relativesimplicity of their interactions. Self-organized quasiperi-odicity has previously been reported in the phases of in-dividual phase oscillators with nonlinear coupling [17]; inthis system, when the coupling K is too low to producefull synchrony, the mean field fails to entrain individualoscillators—however, R ( t ) remains periodic. The behav-ior of Eq. 1 also differs from partial synchronization ob-served in chimera states, in which distinct synchronousand desynchronous subpopulations coexist within an en-semble of oscillators [2, 4, 5, 18]. While chimera statescan produce irregular oscillations of R ( t ), the form andcharacter of chimera states arises specifically from thepresence of a mixture of short- and long-ranged pairwisecouplings. In contrast, the number of short-range in-teractions experienced by an oscillator subject to Eq. 1varies continuously as other oscillators enter and exit itsrepulsive radius d .We therefore ask whether the observed cascades trulyarise due to a critical mechanism, or another source ofchaotic dynamics (such as multiscale correlated fluctu-ations) [19]. We observe that the power spectrum ofthe first-order parameter R ( t ) exhibits 1 /f decay (Fig-ure 1B); additionally, the higher-order Daido order pa-rameters R n ≡ (cid:80) Nj =1 e inθ j ( t ) , n > /f fluctuations in the order parameter oc-curring in the original ”sandpile” model of self-organizedcriticality [20, 21], as well as experimental studies of re-laxation in granular media [22] and frictional interfaces[23]. Additionally, the distribution of desynchroniza- frequency f B burst size s O r d e r p a r a m e t e r R ( t ) P S D ( f ) r -2 s -3/2 -1 -3 -4 -14 P ( s ) Time t A CD
FIG. 1. (A) Schematic of the two steady states of the repul-sive oscillator model. (B) Avalanches in the Daido(1) orderparameter, R ( t ), in ensembles of N = 15 (blue), N = 30(yellow), and N = 60 (turquoise) oscillators. (C) The powerspectral density, and (D) avalanche size distribution for longerversions of the simulations in the previous panel. Shadedranges are standard deviation within each bin of histogram.For this figure, K = 0 . V = 0 . tion event magnitudes, p ( S ), displays power-law behav-ior with critical exponent − /
2, while the distribution oftimes between events, p ( τ ) is exponential—both consis-tent with a self-organized critical mechanism, in whichavalanches of oscillator rearrangements punctuate grad-ual synchronization [19, 23] (Figure 1C, D). Moreover,the distribution of avalanche durations has a critical ex-ponent of − θ space). This clustering increases the steric pres-sure on synchronized oscillators, increasing the probabil-ity that two oscillators will overlap sufficiently to eitherexchange positions, or exert a joint force on a third oscil-lator, thereby triggering a cascade of oscillator rearrange-ments and transient desynchronization. We quantify thiseffect by computing the net repulsive force acting an os-cillator j , due the cumulative effect of other oscillatorsfalling within its interaction radius d , F j = − N (cid:88) k =1 ddθ j V ( θ j , θ k )For an oscillator with an equal number of neighborsfalling within d on its right and left sides, the net forcedue to the repulsive interactions is zero. Thus, F j is pro-portional to the gradient of the local density of oscilla-tors. Figure 2A shows (cid:80) Nj =1 |F j | at each timepoint dur-ing a typical avalanche build-up and release cycle. Duringsynchronization, brief peaks in the net force appear andthen relax, due to oscillators having sufficient space torearrange and maximize their spacing. However, as thedegree of synchronization increases, oscillators becomeless likely to have space to freely rearrange, causing re-pulsive forces to aggregate within the cluster. Eventually,a large-scale rearrangement occurs, triggering a cascadeof rearrangements that break apart the cluster—visibleas a sustained period of non-zero net repulsive forces.The gradual build-up and rapid release of force withinan oscillator cluster mirrors criticality observed in otherphysical process, such as stick-slip friction, earthquakes,and traffic [23]. In the supplementary material, we de-scribe how the continuous motion of oscillators along theunit circle, as well as the formation of ”jams” in regionsof high density (visible as synchronization epochs in thedynamics of R ( t )), motivates us model criticality in thesystem using a cellular automaton inspired by modelsof traffic flows (Supplementary Figure S3). We startwith the classic Nagel-Schreckenberg cellular automatonmodel for traffic on a periodic domain [27]. In this model, N ”cars” are distributed among M sites, represented by alength M string with site values 0 (unoccupied) or 1 (oc-cupied, without overlaps). In each timestep, cars at eachoccupied site move forward a number of units determinedby their current velocity, and then update their velocityusing various acceleration rules (such as random acceler-ation or braking). A central feature of the model is thata car’s velocity is bounded by the number of empty sitespreceding it, thus preventing passing—a constraint thattriggers the formation of large-scale jams at high densi-ties N/M →
1, and which allows the system to exhibitself-organized criticality when the density is sufficient forjams to transiently form and then break apart [28]. Here,we modify this model in two ways: (1) we impose an ac-celeration rule that depends, in part, on the distance ofcars from all other cars, and (2) if the net force on agiven car reaches a critical threshold, it can accelerate ata rate depending on the magnitude and direction of itsnet force (see supplementary material for a more detaileddescription of the model). We find that these two mod-ifications are sufficient to replicate the gradual build-upand rapid breakup of synchronized states that we observein the continuous-time dynamics of the oscillator array(Figure 2). Moreover, the statistical properties of the twomodels are comparable; the power spectral density of R for the traffic model (as calculated by assuming the carstravel on the unit circle) displays 1 /f scaling, as well asexponential waiting times between the events. This sim-ilarity underscores the central role of force build-up andrearrangement in determining the avalanche dynamics,despite the differences between the two models—namely,that only one jam can form at a given time, and only atan exponential rate, in the oscillator model.That rearrangements trigger cascades indicates therole of inter-oscillator repulsion in amplifying local eventsinto global cascades. Avalanches therefore are the re-sult of Eq. 1 becoming extremely sensitive to noise (hereprovided by rearrangements) at intermediate values of V /K . In order to better understand this critical mecha-nism, we next seek to solve the repulsive oscillator modelanalytically in the continuum (large N ) limit. In thislimit, the system comprises an ”oscillator fluid” with thedynamical equation ∂ t f + ∂ θ ( vf ) = 0 . (2) AB Time (1/ω) ( / N ) Σ N | F | d e n s i t y ππ . . . FIG. 2. (A) The summed net force |F j | vs. time during onesynchronization transition ( N = 60). The order parameter R ( t ) is overlaid in turquoise. (B) The density of oscillatorsalong the unit circle, where the motion of the mean field hasbeen subtracted from the individual phases. Color scale runsfrom 0 (black) to 0 . where f ( θ ) is the probability density of oscillators at po-sition θ along the unit circle. The force on this densityis given by a continuum analogue of Eq. 1, v = ω + K (cid:90) sin( θ (cid:48) − θ ) dθ (cid:48) − (1 /N )( V + ξ ( t )) ∂ θ f. (3)The first two term in this expression follow directly fromthe Kuramoto model [30]. The remaining term propor-tional to ∂ θ f follows from the continuum limit of therepulsive interaction (see supplementary material for aderivation of this term from Eq. 1). Similar density gra-dient terms appear in continuum versions of critical sand-pile models [29]; here, the term stabilizes the uniform dis-tribution f = 1 / (2 π ) , R = 0. The noise term ξ ( t ) triggerscascades; rearrangements are abstracted as uncorrelatedstochastic events, with amplitude proportional to the lo-cal oscillator density. We discuss noise in greater detailbelow.In order to solve for the dynamics of R ( t ) subject toEq. 2, we impose a form for f ( θ ) using the Ott-Antonsenansatz [30], which assumes that the density distributionhas the form of a Fourier series, f ( θ ) = 12 π (cid:32) ∞ (cid:88) n =1 a n e inθ + c.c. (cid:33) . Inserting this equation in Eq. 2 and reducing the re-sulting dynamical equation (see supplementary material)shows that the full dynamics are captured by those of theorder parameter,˙ R = − R ( t )2 π (cid:18) πK ( R ( t ) −
1) + (1 /N )( V + ξ ( t )) P ( R ( t )) (cid:19) ,P ( R ( t )) ≡ (cid:0) R ( t ) + R ( t ) + R ( t ) + R ( t ) (cid:1) (4)where R ( t ) e i Φ( t ) ≡ (cid:82) e iθ (cid:48) f ( θ (cid:48) ) dθ (cid:48) . In the zero-noise limit( ξ = 0), Eq. 4 admits two physically meaningful solu-tions, R = 0 , R ≈ (cid:115) − V πKN V πKN , (5)(see supplementary materials for a description of the ex-act nonzero solution, which requires the root of the high-degree polynomial P ( R )). These solutions exchange sta-bility via a transcritical bifurcation at V = πKN (Figure3A, dashed gray trace). Comparison of discrete oscilla-tor simulations for various N show general agreementwith the analytic model below the critical V /K ; be-low this point the simulations converge to a partially-synchronized state with amplitude predicted by Eq. 5(Figure 3A, colored markers). However, as V → πKN ,the dynamics of the discrete simulations undergo anabrupt transition to avalanche dynamics, with transitionpoint and average fluctuation amplitude both dependenton N . These fluctuations do not occur in the continuummodel without noise.To include fluctuations in the continuous model, weconsider the case of non-zero ξ ( t ) in Eq. 3. Analysisof the Ott-Antonsen ansatz and the resulting dynami-cal equations shows that the simplest way that noise canaffect the dynamics of R ( t ) is through a term of the form ξ ( t ) ∂ θ f , as appears in Eq. 3; other approaches, such asadditive or multiplicative noise in Eq. 3, do not affect thedynamics of R ( t ) (see supplementary material). Basedon the numerical simulations of discrete oscillators, wereason that the abrupt breakdown of the phase-lockedstate depends on specific features of rearrangements thatdepend on N , V , and K in a nonlinear way. We treatthis dependence phenomenologically, and impose a func-tional form for the noise term consisting of a sharp cut-off, implying that rearrangements are impossible below acritical force, ξ ( t ) = ξ N H (( V − V a ) /N ) η ( t ) , where ξ is the noise amplitude, and V a sets the minimumdegree of repulsion at which avalanches can occur. Theseare treated as fitting parameters, and η ( t ) is assumed tobe uncorrelated Langevin noise.Given a non-zero value of the noise term ξ ( t ), thedynamics of the order parameter R ( t ) under Eq. 5 canbe simulated as a Stratonovich process (see supplemen-tary material for methods). Comparative simulationsshow that the stochastic continuum model produces qual-itatively similar dynamics to the discrete simulationsover a range of values of V /K and N (Figure 3A,solid colored traces). The minimal V /K for the on-set of avalanches, as well as the average amplitude of R ( t ) during avalanches, both scale inversely with N —underscoring the role of finite-size effects in determiningthe influence of rearrangements on the dynamics. Abovethe critical V /K , the statistical properties of avalanchesin the continuous model show quantitative agreement Repulsion Strength V / (π N K) S t e a d y s t a t e a v e r a g e o r d e r 〈 R 〉 S t e a d y - s t a t e o r d e r p a r a m e t e r R N = 15 N = 30 N = 60DiscreteContinuum, noiseContinuum, no noise
FIG. 3.
Comparison of simulations and analyticmodel.
Results of the analytic continuum model in the ab-sence of noise (black dashed line), discrete simulations of Eq. 1for varying numbers of oscillators (colored markers), and sim-ulations of the continuum model (solid colored lines), withamplitude and noise threshold fitted to the discrete results.Empty circles indicate regions where the discrete and noisycontinuum models both undergo oscillations. with those observed in the discrete simulations, includ-ing an exponential waiting time distribution, a powerlaw size distribution with critical exponent − /
2, and apower spectrum with power law decay and critical ex-ponent − V a , ξ values characterizing rearrangementsin the discrete model at various N . Thus, while our noiseterm does not explicitly model small-scale details of rear-rangements, it shows how they essentially act as a sourceof random variation that is amplified by the dynamicalequation at intermediate values of V /K .We have shown that a minimal generalization of theKuramoto model can produce unexpectedly rich dynam-ics, mirroring those of avalanche models and other sys-tems exhibiting self-organized criticality. We anticipatepotential applications of our model to understandingavalanche-like dynamics seen in real-world oscillator net-works, ranging from power grids [16], to financial mar-kets [15], to neuronal circuits in the brain [14]. Sug-gestively, our observed scaling exponents match those ofexperimental studies of neuronal avalanches, which re-port scaling exponents of approximately 1 . − Supplementary material
For this arXiv submission, the supplementary materialhas been uploaded as an ancillary file. To obtain the SI,navigate to the abstract landing page, and then on theright-hand side select ”Other formats” under the ”Down-load” header. On the resulting page, click the link for”Download source” under the ”Source” header. A num-bered file will download. Either extract the file from theterminal, rename the file with the extension ”.gz”, andthen double-click the file to extract a directory containingthe PDF of supplementary material. [1] J. C. Gonz´alez-Avella, M. G. Cosenza, andM. San Miguel, Localized coherence in two interactingpopulations of social agents, Physica A: StatisticalMechanics and its Applications , 24 (2014).[2] M. J. Panaggio and D. M. Abrams, Chimera states: co-existence of coherence and incoherence in networks ofcoupled oscillators, Nonlinearity , R67 (2015).[3] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchroniza-tion: a universal concept in nonlinear sciences , Vol. 12(Cambridge university press, 2003).[4] Y. Kuramoto and D. Battogtokh, Coexistence of coher-ence and incoherence in nonlocally coupled phase oscilla-tors., Nonlinear Phenomena in Complex Systems , 380(2002).[5] D. M. Abrams and S. H. Strogatz, Chimera states forcoupled oscillators, Physical review letters , 174102(2004).[6] B. Ottino-L¨offler and S. H. Strogatz, Volcano transitionin a solvable model of frustrated oscillators, Physical re-view letters , 264102 (2018).[7] A. Zakharova, M. Kapeller, and E. Sch¨oll, Chimeradeath: Symmetry breaking in dynamical networks, Phys-ical review letters , 154101 (2014).[8] P. A. Tass, Phase resetting in medicine and biology:stochastic modelling and data analysis (Springer Science& Business Media, 2007).[9] O. V. Popovych, C. Hauptmann, and P. A. Tass, Ef-fective desynchronization by nonlinear delayed feedback,Physical review letters , 164102 (2005).[10] H. Hong and S. H. Strogatz, Conformists and contrariansin a kuramoto model with identical natural frequencies,Physical Review E , 046202 (2011).[11] L. Tsimring, N. Rulkov, M. Larsen, and M. Gabbay, Re-pulsive synchronization in an array of phase oscillators,Physical review letters , 014101 (2005).[12] V. Nicosia, M. Valencia, M. Chavez, A. D´ıaz-Guilera, andV. Latora, Remote synchronization reveals network sym-metries and functional modules, Physical review letters , 174102 (2013).[13] E. Omel’chenko, Y. L. Maistrenko, and P. A. Tass,Chimera states: The natural link between coherence andincoherence, Physical review letters , 044105 (2008). [14] J. M. Beggs and D. Plenz, Neuronal avalanches in neocor-tical circuits, Journal of neuroscience , 11167 (2003).[15] L. Bellenzier, J. V. Andersen, and G. Rotundo, Conta-gion in the world’s stock exchanges seen as a set of cou-pled oscillators, Economic Modelling , 224 (2016).[16] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa,Spontaneous synchrony in power-grid networks, NaturePhysics , 191 (2013).[17] M. Rosenblum and A. Pikovsky, Self-organized quasiperi-odicity in oscillator ensembles with global nonlinear cou-pling, Physical review letters , 064101 (2007).[18] E. A. Martens, S. Thutupalli, A. Fourri`ere, and O. Hal-latschek, Chimera states in mechanical oscillator net-works, Proceedings of the National Academy of Sciences , 10563 (2013).[19] G. Boffetta, V. Carbone, P. Giuliani, P. Veltri, andA. Vulpiani, Power laws in solar flares: self-organizedcriticality or turbulence?, Physical review letters , 4662(1999).[20] P. Bak, C. Tang, and K. Wiesenfeld, Self-organized crit-icality: An explanation of the 1/f noise, Physical reviewletters , 381 (1987).[21] H. J. Jensen, K. Christensen, and H. C. Fogedby, 1/fnoise, distribution of lifetimes, and a pile of sand, Phys-ical Review B , 7425 (1989).[22] G. A. Held, D. Solina, H. Solina, D. Keane, W. Haag,P. Horn, and G. Grinstein, Experimental study of critical-mass fluctuations in an evolving sandpile, Physical Re-view Letters , 1120 (1990).[23] H. J. S. Feder and J. Feder, Self-organized criticality in astick-slip process, Physical review letters , 2669 (1991).[24] N. Friedman, S. Ito, B. A. Brinkman, M. Shimono, R. L.DeVille, K. A. Dahmen, J. M. Beggs, and T. C. But-ler, Universal critical dynamics in high resolution neu-ronal avalanche data, Physical review letters , 208102(2012).[25] J. M. Palva, A. Zhigalov, J. Hirvonen, O. Korhonen,K. Linkenkaer-Hansen, and S. Palva, Neuronal long-range temporal correlations and avalanche dynamics arecorrelated with behavioral scaling laws, Proceedings ofthe National Academy of Sciences , 3585 (2013).[26] A. J. Fontenele, N. A. de Vasconcelos, T. Feliciano, L. A.Aguiar, C. Soares-Cunha, B. Coimbra, L. Dalla Porta, S. Ribeiro, A. J. Rodrigues, N. Sousa, et al. , Critical-ity between cortical states, Physical Review Letters ,208101 (2019).[27] K. Nagel and M. Schreckenberg, A cellular automatonmodel for freeway traffic, Journal de physique I , 2221(1992).[28] K. Nagel and M. Paczuski, Emergent traffic jams, Phys-ical Review E , 2909 (1995).[29] L. Gil and D. Sornette, Landau-ginzburg theory of self-organized criticality, Physical review letters , 3991(1996).[30] E. Ott and T. M. Antonsen, Low dimensional behavior oflarge systems of globally coupled oscillators, Chaos: AnInterdisciplinary Journal of Nonlinear Science , 037113(2008). [31] C. W. Lynn and D. S. Bassett, The physics of brainnetwork structure, function and control, Nature ReviewsPhysics , 1 (2019).[32] J. Cabral, E. Hugues, O. Sporns, and G. Deco, Role oflocal network oscillations in resting-state functional con-nectivity, Neuroimage , 130 (2011).[33] D. R. Chialvo, Emergent complex neural dynamics, Na-ture physics , 744 (2010).[34] R. Dodla and C. J. Wilson, Effect of phase responsecurve shape and synaptic driving force on synchroniza-tion of coupled neuronal oscillators, Neural computation , 1769 (2017).[35] A. Levina, J. M. Herrmann, and T. Geisel, Dynamicalsynapses causing self-organized criticality in neural net-works, Nature physics3