Path-dependent Dynamics Induced by Rewiring Networks of Inertial Oscillators
William Qian, Lia Papadopoulos, Zhixin Lu, Keith Wiley, Fabio Pasqualetti, Danielle S. Bassett
PPath-dependent Dynamics Induced by Rewiring Networks of Inertial Oscillators
William Qian, Lia Papadopoulos, Zhixin Lu, Keith Wiley, Fabio Pasqualetti, and Danielle S. Bassett
1, 2, 3, 4, 5, 6, 7 Department of Physics & Astronomy, College of Arts & Sciences, University of Pennsylvania, Philadelphia, PA 19104 USA Department of Bioengineering, School of Engineering & Applied Science, University of Pennsylvania, Philadelphia, PA 19104 USA Department of Mechanical Engineering, University of California, Riverside, CA 92521 USA Department of Neurology, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA Department of Psychiatry, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA Santa Fe Institute, Santa Fe, NM 87501 USA To whom correspondence should be addressed: [email protected] (Dated: September 29, 2020)In networks of coupled oscillators, it is of interest to understand how interaction topology affectssynchronization. Many studies have gained key insights into this question by studying the classicKuramoto oscillator model on static networks. However, new questions arise when network structureis time-varying or when the oscillator system is multistable, the latter of which can occur when aninertial term is added to the Kuramoto model. While the consequences of evolving topology andmultistability on collective behavior have been examined separately, real-world systems such asgene regulatory networks and the brain can exhibit these properties simultaneously. How doesthe rewiring of network connectivity affect synchronization in systems with multistability, wheredifferent paths of network evolution may differentially impact system dynamics? To address thisquestion, we study the effects of time-evolving network topology on coupled Kuramoto oscillatorswith inertia. We show that hysteretic synchronization behavior occurs when the network density ofcoupled inertial oscillators is slowly varied as the dynamics evolve. Moreover, we find that certainfixed-density rewiring schemes induce significant changes to the level of global synchrony, and thatthese changes remain after the network returns to its initial configuration and are robust to awide range of network perturbations. Our findings suggest that the specific progression of networktopology, in addition to its initial or final static structure, can play a considerable role in modulatingthe collective behavior of systems evolving on complex networks.
I. INTRODUCTION
Understanding the emergence of collective behaviorsin systems of dynamical units coupled through complexnetworks remains an important goal in the study of dy-namical systems. The synchronization of coupled oscil-lators is a key example of such behavior [1], and compu-tational models have proven effective in gaining insightinto a number of real-world systems where this phenom-ena occurs, including the synchronization of power grids,the flashing of fireflies, and the dynamics of neuronal net-works [2–6]. More generally, a number of past studieshave focused on the question of how distinct behaviorsof coupled oscillators arise from distinct network topolo-gies, assuming that a given topology remains fixed for agiven system [7–9]. Yet, in many systems, network orga-nization is not static, but rather evolves over time; socialnetworks, neuronal networks, and biological regulatorynetworks are all examples of systems whose interactiontopology can change with time [6, 10–12]. The existenceof such time-evolving network systems motivates an in-vestigation of how specific pathways of network evolutionalter the dynamical behaviors of coupled oscillators.To date, studies of the evolution of coupled oscilla-tors on temporally-evolving networks have often used theKuramoto model. This model is an established systemfor studying synchronization behavior, widely used forits simplicity and analytical tractability. For Kuramotooscillators in the limit of fast network rewiring, prior work has shown that switching between different cou-pling topologies has the same effect as allowing oscilla-tor dynamics to evolve on a network with weights av-eraged over the different switching topologies [13]. Incontrast, another study investigated the effects of net-work connectivity that co-evolves with Kuramoto oscil-lator dynamics, and showed that an adaptive rewiringscheme where oscillators re-route links away from theirneighbors with which they are most in-phase can resultin network topologies that enhance synchronization [14].Other work has demonstrated that networks of phase-lagged Kuramoto oscillators with a biologically-inspiredHebbian learning rule gives rise to unique spatiotempo-ral modes of oscillation [15]. These studies are illustra-tive of the breadth of the field [16–21], which collectivelydemonstrates that the dynamics of Kuramoto oscillatorsdepend appreciably on the type of reconfiguration thatthe coupling network undergoes.In the presence of a bimodal natural frequency distri-bution, phase lags, or frequency-degree correlations, theKuramoto model can exhibit a variety of complex behav-iors, such as hysteretic transitions as a function of cou-pling strength [22–25]. However, under the most genericconditions, the Kuramoto model does not exhibit path-dependent dynamics. Specifically, when adiabatically in-creasing and then decreasing the coupling strength of aKuramoto oscillator population, the value of the orderparameter is typically identical along the forward andbackward transitions. Indeed, for non-negative values ofcoupling, systems of Kuramoto oscillators are monostable a r X i v : . [ n li n . AO ] S e p [26, 27], indicating that the path of network evolutiontaken during oscillator dynamics does not affect levelsof synchrony once a coupling pattern has been fixed. Inother words, history has no effect on the dynamics ofstandard Kuramoto oscillators once transient effects arediscarded.In contrast, systems of second-order Kuramoto oscil-lators are known to be sensitive to history. In particu-lar, the introduction of an inertial term to the Kuramotomodel has been shown to result in highly multistable dy-namics in certain parameter regimes [28, 29]. Unlike inthe standard Kuramoto model, adiabatically tuning thecoupling strength of inertial Kuramoto oscillators resultsin hysteretic synchronization transitions [30]. Specifi-cally, slowly increasing and then decreasing the couplingstrength of inertial Kuramoto oscillators creates a hys-teresis loop in the order parameter, indicating that iner-tial oscillator dynamics can depend significantly on priorconditions. Furthermore, it has been analytically proventhat any nonzero amount of inertia can induce these hys-teresis loops by turning a supercritical bifurcation intoa subcritical bifurcation [31]. These behaviors and thestudies unearthing them collectively suggest that path-dependent dynamics may arise from time-varying con-nectivity in networks of inertial Kuramoto oscillators.Like the standard Kuramoto model, the inertial Ku-ramoto model has also proven insightful for understand-ing real-world systems. The inertial Kuramoto modelwas first introduced to explain synchronization patternsin groups of fireflies [32]. It has since been used exten-sively to study the stability of power grids and the syn-chronization of Josephson junctions [33]. One interpreta-tion of the inertial term is that it extends the Kuramotomodel beyond the simplified and completely overdampedregime, where system dynamics behave analogously tocoupled units oscillating in an extremely viscous medium.The inclusion of inertia allows for both underdamped andoverdamped dynamics, depending on the value of the in-ertial constant. In the context of neuroscience, one sourceof biological support for the inclusion of an inertial termin equations for neural dynamics comes in the form of in-ertia being analogous to inductance [34]. Forms of induc-tance have been observed experimentally in squid axons,and the inclusion of inductive effects in models of neuronshas been shown to allow for richer modulation of tempo-ral dynamics [35, 36]. Other studies have used inertialphase oscillators as a simplified model for the dynamicsof a neuron with an axon and dendrite [37, 38], findingthat incorporating inertia to model dendritic dynamicscan alter responses to stimulation.Given the sensitivity of the inertial Kuramoto model tohistory, we hypothesize that network rewiring alone caninduce path-dependent behavior in networks of inertialKuramoto oscillators. This hypothesis, combined withthe relevance of the inertial Kuramoto model to real-world systems, prompts us to investigate how specificnetwork evolution pathways affect the collective behav-ior of inertial Kuramoto oscillators. In particular, some important questions arise: how does the process of net-work evolution impact the dynamics of a population ofinertial oscillators? In addition, which paths of networkevolution are effective in synchronizing or desynchroniz-ing inertial oscillators, and to what extent do these effectspersist after further network rewiring? Previous studieshave shown that for standard Kuramoto oscillators withstatic connectivity, modular network structures promotelocal synchrony but hinder global synchronization [9, 39].In addition, prior work suggests that high alignment be-tween the eigenvectors of the network Laplacian and theoscillators’ natural frequencies can enhance synchroniza-tion in the standard Kuramoto model [8]. These resultsprompt an analysis of whether modular networks andsynchrony-optimized networks behave similarly for thecase of inertial oscillators at each node, and how pathsof network evolution towards or away from these specialtopologies may affect system dynamics.The remainder of the paper is organized as follows.Section II introduces the inertial Kuramoto model oncomplex networks. In Section III, we examine how thecollective dynamics of inertial oscillators are affected bydifferent network evolution schemes, and we analyze therobustness of effects induced by network rewiring. Weconclude in Section IV with a discussion of our findingsas well as of possible areas for further study. II. THE INERTIAL KURAMOTO MODEL
A system of N inertial Kuramoto oscillators evolvesaccording to the equation m ¨ θ i + ˙ θ i = ω i + α N (cid:88) j =1 A ij sin( θ j − θ i ) , (1)where θ i represents the instantaneous phase of the i thoscillator, ω i is the natural frequency of the oscillator, α is the coupling strength, m is the inertial constant, and A is an N × N unweighted, undirected adjacency matrixrepresenting network connectivity. Note that in the over-damped limit m →
0, the original first-order Kuramotomodel is recovered.The instantaneous level of global synchrony in a popu-lation of oscillators is usually quantified by the modulusof the complex order parameter R ( t ) = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) j =1 e iθ j ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (2)which takes on values ranging from 0 to 1, with highervalues indicating higher levels of phase synchronization.We also introduce the time-averaged order parameter (cid:104) R (cid:105) = 1 T (cid:90) T R + TT R R ( t ) dt, (3) FIG. 1: A schematic of the network rewiring process.
Initially,the network connectivity evolves from G to G f through a series ofunweighted and undirected intermediate graphs (dashed arrow). Thenetwork connectivity then returns to G through the same series ofintermediate graphs (solid arrow). Throughout this rewiring process,oscillator dynamics evolve atop the time-varying connectivity. where T R represents a discarded transient period, and T is the length of the interval over which the order param-eter is averaged. III. SIMULATIONS AND REWIRINGPROCEDURES
We used N = 100 oscillators in all simulations. Ini-tial phases { θ i (0) } were selected at random from [ − π, π ],while initial frequencies { ˙ θ i (0) } and natural frequencies { ω i } were both selected at random from a uniform distri-bution in the interval [ − , G = ( V, E ) and G f = ( V, E f ), we gen-erated a sequence of intermediate graphs { G , G , . . . G f } that determined how network topology would vary overtime. Specifically, let S del = E \ E f denote the edgesin G but not in G f , and let S add = E f \ E be theedges in G f but not in G . We generate the i + 1 st intermediate graph G i +1 from G i by randomly remov-ing ≈ | S del | /f edges in S del ∩ E i from G i , and ran-domly adding ≈ | S add | /f edges in S add ∩ E i to G i , where G = ( V, E ) represents the complement graph of G .After generating the sequence of graphs { G , . . . , G f } ,we carried out a two-step process (Fig. 1). First, we simu-lated the time-evolution of inertial oscillator dynamics asnetwork connectivity evolved from G to G f through theseries of intermediate networks. Then, we continued thetime-evolution of inertial oscillator dynamics as networkconnectivity evolved from G f back towards G throughthe same series of intermediate graphs. For our simula-tions, we use f = 50 transition graphs, and the networkrewiring occurs every l = 5 × time-steps at ∆ t = 0 . T R = ( l × ∆ t ). Additional initial and final transients ofsimulations are also explicitly shown where appropriate. III.A. Varying Network Density
We first demonstrate that hysteretic synchronizationbehavior occurs while increasing and then decreasingnetwork density as oscillator dynamics evolve atop thetime-varying network structure. We generate graphs G through G f by starting with a random Erd˝os–R´enyigraph with an average degree of (cid:104) k (cid:105) = 10. We next addedges uniformly at random until a graph G f with an av-erage degree of (cid:104) k (cid:105) = 40 is reached. Then, we apply therewiring procedure (see Sec. III) to generate all inter-mediate graphs between G and G f . Starting with G we allow the dynamics of the oscillators to run atop thegraph with random initial conditions { θ i (0) } and { ˙ θ i (0) } .Next, we switch the network topology to G , a slightlydenser graph, using the final states of the oscillators fromrunning atop G as the initial conditions for G . This se-quential process is repeated until G f is reached, and isthen continued in reverse until the network returns to G . We hold m and α constant throughout the process.To determine if the presence of inertia gives rise topath-dependent behavior, we allow Kuramoto oscillatordynamics to evolve with and without inertia while wevary network density in the manner described above. Tocompare the two situations, we set the coupling valuesfor the non-inertial and inertial system such that the ini-tial level of synchrony is relatively low and comparablebetween the two cases. As expected, in the absence ofinertia, we find that oscillator dynamics evolve in a re-versible manner throughout the G → G f → G rewiringprocess, suggesting that dynamics are identical for thesame network structures regardless of the network evo-lution pathway taken to reach those structures (Fig. 2a).However, this reversibility is not observed in the presenceof inertia (Fig. 2b), where we instead observe asymmetrictrajectories of both phase synchronization and frequencyentrainment as a function of time.Given the irreversibility of collective dynamics in theinertial case, we hypothesized that a hysteresis loop of thetime-averaged order parameter should form as the net-work density is slowly increased and then decreased backto its initial value. We indeed observe this phenomenawhen inertia is present (Fig. 3b), but not for the stan-dard Kuramoto system (Fig. 3a). Note that this findingis consistent with prior work reporting hysteretic behav-ior in the second-order Kuramoto model while tuningthe coupling strength but holding network connectivityfixed [30]. Indeed, for Erd˝os–R´enyi random networks, itis intuitive that increasing and then decreasing networkdensity should have an effect similar to that of increasingand then decreasing the coupling strength.Our observations thus far leave unanswered the ques-tion of how varying network density in the manner we de-scribe affects oscillator dynamics when both inertia and(a) \ I y - 5 ~ - - - - - - - - - - - - - ~ cc .... (]J I QJ E ro L.. & L.. QJ -c L.. time, t (b) ~i I II I ·' \ I y -5 ,____ ________________ ___, cc .... (]J I QJ E ro L.. & L.. QJ -c L.. time, t FIG. 2: Oscillator dynamics when varying network density, with and without inertia.
Single-instance examples of all oscillators’instantaneous frequencies { ˙ θ i } (top panels) and the global order parameter R ( t ) (bottom panels) as a function of time as network density is firstincreased ( (cid:104) k (cid:105) = 10 →
40) and then decreased ( (cid:104) k (cid:105) = 40 →
10) as the dynamics run. (a)
The process with no inertia ( m = 0, α = 0 . (b) Theprocess with inertia ( m = 2, α = 0 . G to G f and from G f back to G are colored blue and red, respectively. Forvisual clarity, these examples were produced on a faster rewiring timescale than described in the main text ( l = 10 ). Parameters have beenchosen such that minimum and maximum levels of synchrony are comparable in the two cases. (a) __....._ cc: --......- 0.8 ... L QJ -'-' E o.6 ro L ~ L QJ E o.o....___,__-----.-------------.---------,---------,-------,-J
20 30 40 50
Graph Index (b) __......,., cc: ...._... ... L OJ "'-' OJ E ro L ro a.. L OJ "'C L ◄ ◄ ◄ [> [> ◄[> ◄[> [> ◄ ◄ ◄ [> [> [> [> [>[>
20 30 40 Graph Index
FIG. 3: Inertia causes hysteresis in the order parameter when varying network density.
The time-averaged order parameter (cid:104) R (cid:105) asnetwork density is first increased ( (cid:104) k (cid:105) = 10 →
40) and then decreased ( (cid:104) k (cid:105) = 40 → (a) The process with no inertia ( m = 0, α = 0 . (b) The process with inertia ( m = 2, α = 0 . (a) has been chosen to ensure that minimum and maximum levels ofsynchrony are comparable in the two cases. In both panels, the order parameter is plotted against the graph index of the intermediate graphs(rather than the average degree) to emphasize that (cid:104) R (cid:105) -values lying along the same vertical line were obtained from identical oscillatornetwork-connectivity. All curves depict averages over 25 instantiations. strong network coupling are present. To investigate thiscase, we increased the global coupling strength α for boththe inertial and non-inertial system such that the ini-tial synchrony level would be intermediately-valued andagain approximately the same for the two conditions. Inthis manner, the oscillators would initially exhibit par-tially synchronized dynamics in both cases. At high cou-pling, the order parameter for the model without iner-tia continues to exhibit reversible behavior (Fig. 4a). Incontrast, transitions from low-density networks towardsand away from high-density networks create a significantseparation between the forward and backward order pa-rameter curves when inertia is present (Fig. 4). However,the form of the irreversibility at high coupling is quali-tatively different than that observed with moderate cou-pling (Fig. 3b). As opposed to the case with moderatecoupling, no closed hysteresis loop is formed. Rather, athigh coupling, levels of synchrony remain markedly in-creased even after the original lowest-density network isrecovered. The shape of this trajectory suggests that,when the parameters and initial network connectivity of inertial oscillators allow for partially synchronized oscilla-tor dynamics, network evolution towards and away frommore synchronizable network structures may irreversiblyincrease the levels of global synchrony. III.B. Constant Network Density
The changes we report in the rewiring processes de-scribed above could be a function of both the changingtopology and the changing network density. To isolatethe effects of changing topology it is therefore of inter-est to consider rewiring processes that maintain the net-work density. This case is also particularly relevant toreal-world network systems wherein there often exists acost associated with the development and maintenance ofnetwork connections. For example, the energy consumedby synapses in mammalian brains places metabolic con-straints on brain development [40, 41]. A question thenarises whether there exist fixed-density rewiring schemesthat also produce significant separation between the for-(a) (b) __......_ cc: ......_~ 0 8 -- .... L.. QJ -1-J ~ ro L.. ~ ~ L.. QJ ~ [> [> [> [> [> [> [> [> [> [>[> [> [>
0 10 20 30
Graph Index
40 50
FIG. 4: Varying network density with strong coupling.
The time-averaged order parameter (cid:104) R (cid:105) as network density is first increased( (cid:104) k (cid:105) = 10 →
40) and then decreased ( (cid:104) k (cid:105) = 40 →
10) with strong coupling. (a)
The process with no inertia ( m = 0, α = 0 . (b) The processwith inertia ( m = 2, α = 0 . (a) has been chosen to ensure that the minimum and maximum levels of synchrony are comparable in the two cases.All curves depict an average over 25 instantiations. ward and backward order parameter curves. To answerthis question, it is useful to consider network evolutionpathways toward or away from topologies known to behighly synchronizable in the standard Kuramoto model.Prior work has demonstrated that networks of standard,first-order Kuramoto oscillators with optimal alignmentbetween the network Laplacian’s eigenvectors and the os-cillators’ natural frequencies are highly synchronizable[8]. To describe this alignment, let λ j and v j representthe j -th largest eigenvalue and its corresponding eigen-vector of the network Laplacian L ij = δ ij k i − A ij , where k i is the degree of node i . Following Ref. [8], in thestrongly synchronized regime, minimizing the synchronyalignment function J ( ω , L ) = 1 N N (cid:88) j =2 λ − j (cid:104) v j , ω (cid:105) , (4)serves to maximize the global order parameter R in thestandard Kuramoto model.Applying this approach, we generated synchrony-aligned networks of a given average degree (cid:104) k (cid:105) via ahill-climbing algorithm with the procedure described inRef. [8] (See Supplementary Material). Then, usingErd˝os–R´enyi graphs for G and synchrony-aligned graphsfor G f , we simulated the network evolution pathway de-fined by G → G f → G while maintaining a fixed net-work density ( (cid:104) k (cid:105) = 20).We begin by considering a situation of relatively highglobal coupling (Fig. 5). For the standard Kuramotomodel, rewiring towards synchrony-aligned networks in-creases the order parameter substantially, and as ex-pected, the synchrony level returns to its initial valuealong the same path as the network returns back to theoriginal Erd˝os–R´enyi graph (Fig. 5a). When inertia is in-corporated (and the coupling strength adjusted to obtaina similar level of initial synchrony), we again find thatnetwork evolution towards synchrony-aligned graphs en-hances the order parameter, and the system nears perfect synchrony at G f = G ∗ (Fig. 5b). Moreover, in contrastto the non-inertial case, the transition from Erd˝os–R´enyigraphs toward and away from synchrony-aligned graphscreates a significant separation between the forward andbackward order parameter curves. However, this rewiringdoes not elicit a closed hysteresis loop. Similar to the caseof varying network density with strong coupling (Fig. 4)for the second-order Kuramoto model, we find that thesteady-state level of synchrony is maintained at a signif-icantly higher value even after the system returns to theoriginal Erd˝os–R´enyi graph.We next quantify how the steady-state synchronizationgap R ssf − R ss changes over a swath of the inertia-couplingparameter space (0 . ≤ α ≤ .
5, 1 . ≤ m ≤ . R ss and R ssf represent the initial and final time-averaged order parameters, respectively, after discard-ing a long transient period (Fig. 6a). At high cou-pling and low inertia, there is little steady-state sepa-ration between the initial steady-state order parameter R ss and the final steady-state order-parameter R ssf . Thisbehavior is expected because oscillators with high cou-pling and low inertia reach close-to-perfect synchrony onthe initial Erd˝os–R´enyi connection topology (Fig. 12a);rewiring towards synchrony-aligned networks can there-fore only induce a small enhancement of the order pa-rameter (Fig. 6b). As detailed further in the followingparagraph, low coupling and high inertia also result innegligible steady-state separations R ssf − R ss (e.g., theparameter combination denoted by the green triangle inFig. 6d). In contrast, in the regime of moderate couplingand moderate inertia, the network evolution process hasa clear sustained effect on the system’s collective dynam-ics as reflected in the steady-state synchronization gap R ssf − R ss .To dig deeper into the behavior of the system, we nextconsider the fact that for some parameter combinations,the Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyinetwork evolution could induce a hysteresis loop but nota steady-state synchrony gap. To assess this more nu-(a) _........._ cc ----.-- 0.8 ... L OJ -1-' ~ ro L ~ L OJ "E ~
0 10 20 30 40 50
Graph Index (b) * ___....._ cc: ------- 0.8 ... L.. QJ -1-1 E o.6 C> ro L.. & L.. QJ °E o.o~-----------r----------.-------,----------,----------,--,
0 10 20 30 40 50
Graph Index
FIG. 5: Synchrony gains through network evolution at constant density and intermediate coupling.
The time-averaged orderparameter (cid:104) R (cid:105) as the network connectivity of inertial oscillators evolves from an Erd˝os–R´enyi graph towards a synchrony-aligned graph atconstant density ( (cid:104) k (cid:105) = 20), followed by reversal along the same set of intermediate networks until the original connectivity graph is recovered. (a) The process with no inertia ( m = 0, α = 0 . (b) The process with inertia ( m = 2, α = 0 . ( a) and (b) . All curves depict averagesover 25 instantiations. (a) Ο * Δ (b) Ο cc ---.--- 0.8 L.. OJ °E o.0~ ------------,-----------.------~--------r----------,-.
0 10 20 30 40 50
Graph Index (c) Ο * Δ (d) Δ FIG. 6: Effects of fixed-density rewiring from Erd˝os–R´enyi to synchrony-aligned networks and back, as a function of thecoupling strength and inertia.
We show the behavior of the order parameter (cid:104) R (cid:105) as rewiring occurs from Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi networks while maintaining a constant density of (cid:104) k (cid:105) = 20. (a) Gain in the level of global synchrony as a result of the rewiringprocess from Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi over a portion of the α - m parameter space. (b) The time-averaged orderparameter (cid:104) R (cid:105) throughout the rewiring process with parameters m = 1 . α = 0 . (c) Area between the forward (Erd˝os–R´enyi → synchrony-aligned) and backward (synchrony-aligned → Erd˝os–R´enyi) order parameter (cid:104) R (cid:105) curves over a portion of the α - m parameter space.Area is defined such that it is positive when the backward order parameter curve is above the forward curve. Note that we discardedcontributions from initial and final transients. (d) The time-averaged order parameter (cid:104) R (cid:105) throughout the rewiring process with parameters m = 2 . α = 0 .
24. The inset reports results from the same process but with the standard Kuramoto model, where the coupling strength hasbeen chosen such that initial levels of synchrony are comparable to the inertial case ( α = 0 . anced behavior, we calculated the area between the for-ward and backward order parameter (cid:104) R (cid:105) curves resultingfrom Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyinetwork evolution, over the same inertia-coupling param-eter space (Fig. 6c). For this analysis, the area is definedsuch that it is positive when the backward order parame-ter curve is above the forward curve (and we again ignorecontributions from initial and final transient periods).Interestingly, we observe a regime at low coupling andhigh inertia where no steady-state synchrony gap is pro-duced, but hysteresis loops of negative area are formed(Fig. 6d). That is, the order parameter actually decreases upon rewiring towards synchrony-aligned networks, andthen increases back to its initial value along the reversenetwork evolution pathway.This type of dynamical trajectory could be a naturalconsequence of the fact that the derivation of the syn-chrony alignment function used to produce synchrony-aligned networks employs the approximation of thestrong synchrony regime [8]. This fact in turn sug-gests that synchrony-aligned networks could be inef-fective in promoting synchronization when synchroniz-ability is already low as a result of parameter choices.However, it is also possible that the ineffectiveness ofsynchrony-aligned networks (and the emergence of hys-teresis loops characterized by negative area) in some pa-rameter regimes is a consequence of inertia rather thaninitial synchrony levels alone. Indeed, oscillators coupledthrough synchrony-aligned networks exhibit strong sen-sitivity to initial synchrony levels when inertia is present(see Supplementary Figs. 13 and 14). To probe this possi-bility further, we assessed the Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi rewiring process using the stan-dard Kuaramoto model, with a coupling strength chosento make initial synchrony levels comparable to that ofthe main panel in Fig. 6d. Consistent with the idea thatinertia is responsible for the ineffectiveness of synchrony-aligned networks in some parameter regimes, we findthat standard Kuramoto oscillators with similar levelsof initial synchrony still synchronize well when they arerewired towards a synchrony-aligned topology (Fig. 6dinset).
III.C. Further Network Perturbation
We next sought to quantify the robustness of increasesin synchronization due to network rewiring. We be-gan by taking the final, high-synchrony states of iner-tial oscillators obtained after rewiring towards and awayfrom synchrony-aligned graphs at intermediate coupling(Fig. 5b), and using these final states as initial conditionsfor a set of new simulations. The initial topologies G of these new simulations were the original Erd˝os–R´enyigraphs used and the parameters remained fixed at m = 2, α = 0 .
3. While holding network density constant, wethen rewired oscillator connectivity towards and awayfrom one of four final network topologies G f : 1) other Erd˝os–R´enyi graphs, 2) synchrony-misaligned graphs, 3)random modular graphs, or 4) frequency modular graphs(see Supplementary Materials for details on graph con-struction). These final network structures were chosen soas to assess the level of topological perturbation neededto effectively desynchronize systems of inertial oscillatorsin a high-synchrony state induced by a particular networkevolution history.We hypothesized that further network evolution to-ward and away from other Erd˝os–R´enyi graphs wouldhave little effect on global synchrony, and would pre-serve most of the prior synchrony gains. In contrast,we expected that networks with modular organizationmay effectively erase global synchrony gains resultingfrom a specific path of network evolution. In particular,we conjectured that frequency modular graphs—graphscreated by assigning oscillators with similar natural fre-quencies to the same module—would be most effective inperturbing rewiring-induced gains in global synchrony.Consistent with intuition, we found that rewiring to-wards other Erd˝os–R´enyi graphs had little-to-no effecton levels of synchrony (Fig. 7a). In every numericalexperiment, the system remained at the enhanced syn-chrony level acquired under Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi network evolution, with littledeviation throughout both the forward and backwardrewiring trajectories. This finding suggests that gains insynchrony due to network rewiring through synchrony-aligned graphs are quite robust to further random net-work perturbations.Still, it remains unclear as to which topologies mightbe able to desynchronize inertial oscillators with syn-chrony gains resulting from a particular network evolu-tion history. To probe this question further, one naturalidea is to use synchrony- misaligned graphs for the G f network structures (blue curves in Fig. 5a). Such net-works are constructed by maximizing (rather than min-imizing) the synchrony alignment function (Eq. 4), andthus should theoretically be quite difficult to synchro-nize. We found that rewiring trajectories towards thesynchrony-misaligned graphs induced partial desynchro-nization of the oscillators. Interestingly, though, we ob-served clear irreversibility in the order parameter as werewired from the synchrony-misaligned graphs back tothe original Erd˝os–R´enyi graphs. Specifically, the systemdid not fully return to the baseline synchrony level ob-tained with random initial conditions (green dotted linein Fig. 7a).For our final analysis, we wished to investigate whethernetworks that promote local synchrony can effectivelyreset gains in global synchrony resulting from networkhistory. To do so, we considered modular networks,which have topologies known to favor local synchronyover global synchrony. Rewiring towards both randommodular and frequency modular G f graphs greatly re-duced the global synchrony of the oscillators (see for-ward trajectories in Fig. 7b). However, only evolu-tion towards frequency modular graphs gave rise to ef-(a) -----~~► ~ ............. L.. QJ ~ ... " ... - Other ER
Synchrony-Mi sa I ig ned o.o_l__---0 +-------
Graph Index (b) __....._ cc: ---.-- 0.8
L.. OJ ~ Modular Frequency Modular o.ol___---------+---0-------10--2r-o--3~0--4--ro-------=r:so
Graph Index
FIG. 7: The robustness of gains in global synchrony from network evolution through synchrony-aligned graphs.
We considered thefinal states of oscillators from Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi shown in Fig. 5b. We then used these final states as initialconditions for a simulation that began at the same Erd˝os–R´enyi G connectivity and evolved toward and away from one of four final networktopologies G f . (a) Results of simulations in which G f was set to be another Erd˝os–R´enyi graph (orange), or in which G f was set to be asynchrony-misaligned graph (blue). (b) Results of simulations in which G f was set to be a modular graph (orange), or in which G f was set to bea frequency modular graph (blue). Parameters m = 2, α = 0 . G connectivity and random initial conditions; that is, the order parameter prior to anynetwork rewiring. All curves depict averages over 25 instantiations. fects that remained even after Erd˝os–R´enyi G con-nectivity was recovered (see backward trajectories inFig. 7b). This behavior might occur because, in addi-tion to discouraging global synchronization, frequencymodular graphs are more prone to allowing oscillatorsto evolve onto the cluster synchronization manifold, re-setting much of the history of global synchrony (see Sup-plementary Fig. 15). Note also that rewiring towardsfrequency-modular graphs yields a slightly lower level ofglobal synchrony than rewiring to modular graphs, whichmay also play a role in determining the final level of syn-chrony after rewiring back to the Erd˝os–R´enyi networks.In sum, our results indicate that the extent to whichenhanced synchrony is maintained after further networkrewiring depends on more than just how effectively the G f topology reduces global synchrony during its pres-ence. In particular, global synchrony while G f topologywas present was higher for synchrony-misaligned graphsthan for random modular graphs, but as the system re-turned to its initial topology along the backward transi-tion, the synchrony-misaligned pathway ultimately led tomore sustained desynchronization. This pattern of find-ings suggests that the specific pathway of network evo-lution can play a key role in modulating the collectivedynamics of coupled oscillators, beyond just the imme-diate effects that different network structures have onsynchrony. IV. DISCUSSION
In this paper, we investigated how various routes ofnetwork evolution affect systems of coupled oscillators.Networks of standard Kuramoto oscillators are monos-table [26, 27]. Therefore, under generic conditions theyare unaffected by network rewiring processes; the post-transient global synchrony of oscillators at a given time are a function of just the network structure present atthat time. However, it may not be the case that thispath-independent behavior persists in systems of inher-ently multistable oscillators. To probe the question ofwhether multistability in the dynamics of individual os-cillators gives rise to path-dependent behavior under net-work rewiring, we used the inertial Kuramoto model,which adds an inertial term to the standard Kuramotomodel and consequently exhibits sensitivity to initial con-ditions [28, 29]. Prior work has shown that systems ofinertial oscillators can exhibit hysteretic synchronizationtransitions as the coupling strength is increased and thendecreased [30]. For networked oscillators, this tuning ofthe coupling strength can be regarded as a global scalingof the strength of connections that leaves the networktopology intact.However, in many systems, it is the network organiza-tion itself—i.e., where edges exist or do not exist—that isdynamic, rather than the overall strength of each connec-tion [10–12]. In these cases, it then becomes interesting toask whether networks of inertial oscillators exhibit path-dependent dynamics induced by changes in network or-ganization alone. To answer this question, we developeda network rewiring procedure that isolates the effects ofnetwork evolution history. Specifically, we evolved thenetwork connectivity of systems of coupled inertial (andnon-inertial) oscillators towards a pre-specified final net-work structure, and then we reversed the rewiring processalong the same path. In this way, any path-dependentsynchronization behavior would be reflected as asymme-tries of the order parameter between the forward andbackward network rewiring trajectories.We first investigated the effects of slowly increasingthe network density of random graph topology and thenreversing the evolution until the original graph was recov-ered. For oscillators with moderate inertia and coupling,we found that this density-varying process could inducehysteretic synchronization behavior, with oscillators pre-ferring to stay in a more globally synchronized state formore of the backward rewiring process than the forwardrewiring process. This finding is in line with the afore-mentioned work on hysteretic behavior of inertial Ku-ramoto oscillators with tuning of the coupling strength[30]; it is natural to expect that increasing the densityof connections in non-sparse random networks will yieldsimilar effects to that of globally increasing the strengthof connections between oscillators. Expanding upon thisunderstanding, the formation of a hysteresis loop as net-work density is varied demonstrates the existence of situ-ations where identical network connectivity patterns giverise to different oscillator dynamics due to differences innetwork evolution history alone. Going a step further,we then analyzed the case of varying oscillator networkdensity at high coupling, where we uncovered a qualita-tively unique form of path-dependent behavior in whichnetwork rewiring results in irreversible gains in globalsynchrony.To isolate the role of network topology from density-driven effects in path-dependent behavior, we next stud-ied how inertial oscillators behave when their networktopology is rewired at constant density. Specifically, wegenerated networks known to be highly synchronizablefor the standard Kuramoto model, and analyzed the ef-fects of constant-density rewiring of initially randomly-coupled networks of inertial oscillators toward and thenaway from these synchrony-aligned networks. Notably,we found that even when density is held constantthroughout the network rewiring process, the dynam-ics of inertial oscillators can depend on previous networkevolution history, and that this path-dependence is sig-nificant for a considerable portion of the inertia-couplingparameter space. In addition, gains in synchrony dueto constant-density rewiring were robust to a number ofsubsequent network perturbations, with near-completereversal of effects occurring only in the extreme caseof further rewiring towards highly cluster-synchronizablenetworks. Collectively, these results demonstrate thatvariations in topology alone can drive path-dependentdynamics of inertial oscillators, and that the resultingeffects typically persist even with further network evolu-tion.
Opportunities for extensions and expansions.
Ourresults prompt a number of interesting directions for fur-ther investigation, particularly related to the nature ofthe graphs studied, the nature of the rewiring process,and expansions to other formalisms and models. Firstwe note that we considered networks with binary connec-tivity only. While a reasonable place to begin, it wouldbe interesting in the future to study inertial oscillatorsevolving on time-varying weighted networks as well. Inparticular, the weights could be made to vary continu-ously in time [5, 42, 43], possibly making stability analy-sis of oscillator dynamics evolving along different networkevolution pathways more analytically tractable. Such anextension could increase the relevance of our observations to real world systems, which are typically characterizedby a high density of edges whose weights can vary overseveral orders of magnitude [44–47].Second, we studied the effects of rewiring an initial net-work topology towards a final network topology, wherethe necessary edges to be moved were rewired in a ran-dom order. A future study could consider whether theorder in which edges are rewired plays a significant role inthe development of path-dependent behavior. Prior workhas shown that enforcing certain relationships betweenpairwise differences in the natural frequencies of oscilla-tors and their coupling patterns may promote more com-plex oscillator dynamics, such as explosive synchroniza-tion [48, 49]. Therefore, it is possible that first placingedges between oscillators with the least similar naturalfrequencies may greatly affect how synchronization devel-ops in both the forward and backward rewiring processes.Moreover, while the rewiring process we used was conve-nient for illustrating potential path-dependent behavior,it was controlled in the sense that only the edges that ul-timately needed to be added or removed to arrive at somefinal network structure were rewired, and each relevantedge was only altered once. It would be interesting to seehow inertial oscillators behave when the rewiring processoccurs in a more organic manner, such as by allowingfor all edges to be added and pruned repeatedly whilestill maintaining that some pre-determined final networkstructure is eventually reached.A third possible area for future study is to considerhow path-dependence arises in a system of inertial oscilla-tors adhering to an adaptive rewiring scheme, where thestates of the oscillators themselves inform the networkrewiring process [14, 19, 20, 43, 50–52]. In particular, in-vestigating path-dependence in systems of inertial oscil-lators under Hebbian or anti-Hebbian adaptive rewiring[53–55] may be insightful for understanding whether net-work evolution path-dependence plays a significant rolein the development of neuronal networks. It is well knownthat the human brain undergoes a variety of structuralchanges during development [56–58], not only in synap-tic density but also in topological characteristics such asdegree heterogeneity, clustering, and modularity [59, 60].Such changes complicate any inferences drawn from theexisting levels of synchronization, which may be both afunction of the current network topology and a functionof the network’s developmental history.Another interesting avenue for future work would beto examine how multilayer oscillator networks behave inthe presence of inertia and network rewiring. Prior workhas shown that multilayer oscillator networks exhibit avariety of rich behaviors [61]. For example, explosivesynchronization may occur in two-layer networks withadaptive coupling [62, 63]. Such behavior has even beenshown to occur when an oscillator network is coupled toan entirely different dynamical process, such as a nutrienttransport layer [64]. Future work could fruitfully investi-gate how multilayer oscillator networks evolve in the pres-ence of inertia, and examine the role of path-dependence0in multilayer oscillator networks with dynamic connec-tivity.Finally, it is worth noting that other variants of theKuramoto model can also exhibit multistability, suchas the Kuramoto-Sakaguchi model with time-delayedcoupling [22–25]. It would thus also be insightful toinvestigate the interplay between multistability arisingfrom these alternative means and the network evolutionof oscillator connectivity.
Conclusion.
Discerning the effects of dynamic networkorganization on the collective behavior of coupled dy-namical subunits remains an important area of study,with implications for a number of physical and biolog-ical systems [3–5, 10–12]. To understand whether os-cillators coupled through time-varying networks can besignificantly affected by the history of the coupling net-work, or if only the final network structures obtainedfrom the network evolution process are relevant, we havestudied the impact of various network rewiring pathwayson systems of inertial Kuramoto oscillators. We haveobserved many situations throughout this study wheremarkedly different oscillator dynamics and levels of syn-chrony can arise from the same network structure solelydue to differences in prior connectivity patterns. In par-ticular, we have found that network rewiring can drivehysteretic synchronization behavior of inertial oscillatorsvia increasing and decreasing network density, and thateven changes in topology alone from constant-densitynetwork rewiring can induce path-dependent behavior.These findings demonstrate that beyond the overdampedlimit, pathways of network evolution themselves play animportant role in regulating the behavior of networks ofcoupled subunits, and require consideration when study-ing the dynamics of systems evolving over complex net-works, and when devising strategies for their control.
V. ACKNOWLEDGMENTS
WQ acknowledges support from the Vagelos programat the University of Pennsylvania. FP and DSB ac-knowledge support from the National Science Founda-tion, through a collaborative grant funding mechanism(IIS-1926757). LP, ZL, KW, and DSB also acknowl-edge further support from the Paul G. Allen FamilyFoundation, the National Science Foundation (PHY15-54488, DMR-1420530), and the Army Research Office(W911NF-16-1-0474, W011MF-191-244). The content issolely the responsibility of the authors and does not nec-essarily represent the official views of any of the fundingagencies.
VI. CITATION DIVERSITY STATEMENT
Recent work in several fields of science has identi-fied a bias in citation practices such that papers from women and other minorities are under-cited relative tothe number of such papers in the field [65–70]. Here wesought to proactively consider choosing references thatreflect the diversity of the field in thought, form of con-tribution, gender, and other factors. We obtained pre-dicted gender of the first and last author of each ref-erence by using databases that store the probability ofa name being carried by a woman [65, 71]. By thismeasure (and excluding self-citations to the first andlast authors of our current paper), our references con-tain 9.4% woman(first)/woman(last), 7.8% man/woman,15.6% woman/man, 56.2% man/man, and 10.94% un-known categorization. This method is limited in that a)names, pronouns, and social media profiles used to con-struct the databases may not, in every case, be indicativeof gender identity and b) it cannot account for intersex,non-binary, or transgender people. We look forward tofuture work that could help us to better understand howto support equitable practices in science.1
VII. SUPPLEMENTARY MATERIAL
Generating synchrony-aligned networks : We gen-erated synchrony-aligned networks of a given average de-gree (cid:104) k (cid:105) via a hill-climbing algorithm [8] as follows: start-ing with an Erd˝os–R´enyi graph G with average degree (cid:104) k (cid:105) , an edge is deleted at random and replaced with anedge between two randomly chosen nonadjacent nodes,producing G (cid:48) . If J ( ω , L (cid:48) ) < J ( ω , L ), then G (cid:48) is accepted;otherwise, the original graph G is retained. This processwas repeated for 2 × iterations until a synchrony-aligned graph G (cid:63) was obtained. Synchrony-misalignednetworks were generated following the same procedure,but instead sought to maximize J ( ω , L ). Generating modular networks : Random modulargraphs were generated by randomly assigning each os-cillator to one of five modules, and then preferentiallyadding edges within modules so that 90% of edges wereintra-modular. Frequency modular graphs were gener-ated similarly, but instead assigned modules to oscil-lators based on the similarity of their natural frequen-cies. This assignment was done by dividing the interval[ ω min , ω max ] evenly into five intervals of the same length,where ω min and ω max correspond to the minimum andmaximum natural frequency, respectively. Each oscilla-tor was then assigned to a module corresponding to theone of five intervals that its natural frequency fell within. Assessing cluster synchrony:
We quantify the time-averaged pairwise synchrony between oscillators i and j in Figs. 15c and 15d by (cid:104) R i,j (cid:105) = 1 T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) T R + TT R
12 ( e iθ j ( t ) + e iθ i ( t ) ) dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (5)These values were computed while the modular orfrequency modular G f topology was present duringErd˝os–R´enyi → G f → Erd˝os–R´enyi network evolution.As in Fig. 7b, the initial conditions were obtained fromthe high-synchrony states resulting from Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi evolution (Fig. 5b).2(a)(b)
FIG. 8: Oscillator dynamics while varying network density.
Single-instance examples of the phases and frequencies of the oscillators vs.time t as network density is first increased ( (cid:104) k (cid:105) = 10 →
40) and then decreased ( (cid:104) k (cid:105) = 40 → l = 1 × ), with every tenth phase and frequency shown. The conditions for each simulation were chosen to make minimum andmaximum levels of synchrony comparable: (a) without inertia ( m = 0, α = 0 .
15) (b) with inertia ( m = 2, α = 0 . FIG. 9: Oscillator dynamics while undergoing constant-density rewiring.
Single-instance examples of the phases and frequencies of theoscillators vs. time t as oscillators undergo Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi rewiring at constant density ( (cid:104) k (cid:105) = 20), with andwithout inertia. The red dotted line indicates the beginning of reverse network evolution. The oscillators are sorted in ascending order of naturalfrequencies, with each row corresponding to the time series of a particular oscillator. For visual clarity, these examples were produced on a fasterrewiring timescale than described in the main text ( l = 1 × ), with every tenth phase and frequency shown. The conditions for each simulationwere chosen to make minimum and maximum levels of synchrony comparable: (a) without inertia ( m = 0, α = 0 .
15) (b) with inertia ( m = 2, α = 0 . ....-..... ex: ._____. ◄ ◄ ◄ ◄ ◄ ◄ ◄ ◄ ◄ / [> [> J> [> ◄ [> [> L.. [> (]) ◄ J> ""CO2 ◄ c::t> L.. . ~
40 60
80 100
Graph Index (b) ------- cc: ------ -. L.. QJ -1--1 QJ E ro L.. ro a_ L.. QJ -c L.. ◄ ~ • • ~ I I ' • ◄ • ' ◄ ,).Jflltl!llllilfflhlU!IIl~111II I ' - o.o-----~--------.-----------,---~---o
200 300
Graph Index
FIG. 10: Varying network density at different timescales.
The time-averaged order parameter (cid:104) R (cid:105) as network density is first increased( (cid:104) k (cid:105) = 10 →
40) and then decreased ( (cid:104) k (cid:105) = 40 →
10) at different values of l and f , where l × f is kept constant to fix the total integration time.Both figures were produced with parameters m = 2, α = 0 . l = 2 . × and f = 100. (b) Here, l = 5 × and f = 500. All curves depict averages over 25 instantiations. (a) .............. cc: --___.- 0.8 -. L.. OJ -'-■-J ~ c> ro L.. ~ L.. OJ ~ I> I> I> o.o~------+---------r---------.---------r---------.---
0 20 40 60 80 100
Graph Index (b) ex: ____, 0.8
L.. OJ °E -
0 100 200 300 400 500
Graph Index
FIG. 11: Constant-density rewiring at different timescales.
The time-averaged order parameter (cid:104) R (cid:105) as Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi rewiring occurs at constant density ( (cid:104) k (cid:105) = 20) for different values of l and f , where l × f is kept constant to fix the totalintegration time. Both figures were produced with parameters m = 2, α = 0 . l = 2 . × and f = 100. (b) Here, l = 5 × and f = 500. All curves depict averages over 25 instantiations. (a) (b) FIG. 12: Supplementary analysis of constant-density rewiring over the inertia-coupling parameter space. (a)
The initialsteady-state order parameter R ss prior to Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi rewiring. (b)
The final steady-state orderparameter R ssf after Erd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi rewiring.
Initial Order Parameter, R (0) R ss (b) Initial Order Parameter, R (0) R ss FIG. 13: The effect of the initial order parameter on networks of inertial oscillators with static Erd˝os–R´enyi connectivity.
Thesteady-state order parameter R ss resulting from static Erd˝os–R´enyi connectivity for 100 different sets of initial phases, with parameters (a) α = 0 . m = 2, and parameters (b) α = 0 .
24 and m = 2 .
3, corresponding to the parameters of Figs. 5b and 6d, respectively. Initial phaseswere sampled so that various levels of initial synchrony are represented. (a)
Initial Order Parameter, R (0) R ss (b) Initial Order Parameter, R (0) R ss FIG. 14: The effect of the initial order parameter on networks of inertial oscillators with static synchrony-aligned connectivity.
The steady-state order parameter R ss resulting from static synchrony-aligned connectivity for 100 different sets of initial phases, with parameters (a) α = 0 . m = 2, and parameters (b) α = 0 .
24 and m = 2 .
3, corresponding to the parameters of Figs. 5b and 6d, respectively. Initialphases were sampled so that various levels of initial synchrony are represented. i (c)(b) i (d) FIG. 15: Modular networks and their synchronizability during network evolution.
Here we show examples of modular networks usedfor G f in Fig. 7b. (a, b) Examples of a random modular network and a frequency modular network used. Natural frequencies of the oscillatorsare shown. (c,d)
The pairwise synchrony (cid:104) R i,j (cid:105) corresponding to these examples is shown for all pairs of vertices while the G f topology ( (c) random modular, or (d) , frequency modular) is present during Erd˝os–R´enyi → G f → Erd˝os–R´enyi network evolution. Rows and columns aresorted by modules, with modules delimited on the axes. Initial conditions were obtained from the high-synchrony states resulting fromErd˝os–R´enyi → synchrony-aligned → Erd˝os–R´enyi evolution (Fig. 5b). [1] Alex Arenas, Albert D´ıaz-Guilera, Jurgen Kurths, YamirMoreno, and Changsong Zhou. Synchronization in com-plex networks. Physics reports , 469(3):93–153, 2008.[2] Rapha¨el Sarfati, Julie Hayes, ´Elie Sarfati, and Orit Pe-leg. Spatiotemporal reconstruction of emergent flash syn-chronization in firefly swarms via stereoscopic 360-degreecameras. bioRxiv , 2020.[3] Adilson E. Motter, Seth A. Myers, Marian Anghel, andTakashi Nishikawa. Spontaneous synchrony in power-gridnetworks.
Nature Physics , 9(3):191–197, Oct 2013.[4] John Buck. Synchronous rhythmic flashing of fireflies. ii.
The Quarterly Review of Biology , 63(3):265–289, 1988.[5] D. Cumin and C.p. Unsworth. Generalising the ku-ramoto model for the study of neuronal synchronisa-tion in the brain.
Physica D: Nonlinear Phenomena ,226(2):181–196, 2007.[6] Rabiya Noori, Daniel Park, John D. Griffiths, SonyaBells, Paul W. Frankland, Donald Mabbott, and JeremieLefebvre. Activity-dependent myelination: A glial mech-anism of oscillatory self-organization in large-scale brainnetworks.
Proceedings of the National Academy of Sci-ences , 117(24):13227–13237, Jun 2020.[7] Francisco A. Rodrigues, Thomas K. DM. Peron, PengJi, and J¨urgen Kurths. The kuramoto model in complexnetworks.
Physics Reports , 610:1–98, Jan 2016.[8] Per Sebastian Skardal, Dane Taylor, and Jie Sun. Opti-mal synchronization of complex networks.
Physical Re-view Letters , 113(14), 2014.[9] Per Sebastian Skardal and Juan G. Restrepo. Hierarchi-cal synchrony of phase oscillators in modular networks.
Physical Review E , 85(1):016208, Jan 2012.[10] Guillaume Laurent, Jari Saram¨aki, and M´arton Kar-sai. From calls to communities: a model for time-varying social networks.
The European Physical JournalB , 88(11):301, Nov 2015.[11] Vince D. Calhoun, Robyn Miller, Godfrey Pearlson, andTulay Adalı. The chronnectome: Time-varying connec-tivity networks as the next frontier in fmri data discovery.
Neuron , 84(2):262–274, Oct 2014.[12] Sophie L`ebre, Jennifer Becq, Fr´ed´eric Devaux,Michael PH Stumpf, and Ga¨elle Lelandais. Statisticalinference of the time-varying structure of gene-regulationnetworks.
BMC Systems Biology , 4(1):130, Sep 2010.[13] Marco Faggian, Francesco Ginelli, Fernando Rosas, andZoran Levnaji´c. Synchronization in time-varying randomnetworks with vanishing connectivity.
Scientific Reports ,9(1), 2019.[14] Lia Papadopoulos, Jason Z. Kim, J¨urgen Kurths, andDanielle S. Bassett. Development of structural corre-lations and synchronization from adaptive rewiring innetworks of kuramoto oscillators.
Chaos: An Interdisci-plinary Journal of Nonlinear Science , 27(7):073115, 2017.[15] L. Timms and L. Q. English. Synchronization in phase-coupled kuramoto oscillator networks with axonal delayand synaptic plasticity.
Physical Review E , 89(3):032906,Mar 2014.[16] Wu-Jie Yuan and Changsong Zhou. Interplay betweenstructure and dynamics in adaptive complex networks:Emergence and amplification of modularity by adaptivedynamics.
Physical Review E , 84(1):016116, Jul 2011.[17] Wu-Jie Yuan, Jian-Fang Zhou, Qun Li, De-Bao Chen, and Zhen Wang. Spontaneous scale-free structure inadaptive networks with synchronously dynamical linking.
Physical Review E , 88(2):022818, Aug 2013.[18] Changsong Zhou and J¨urgen Kurths. Dynamical weightsand enhanced synchronization in adaptive complex net-works.
Physical Review Letters , 96(16):164102, Apr 2006.[19] Takaaki Aoki and Toshio Aoyagi. Co-evolution of phasesand connection strengths in a network of phase oscilla-tors.
Physical Review Letters , 102(3):034101, Jan 2009.[20] Takaaki Aoki and Toshio Aoyagi. Self-organized networkof phase oscillators coupled by activity-dependent inter-actions.
Physical Review E , 84(6):066109, Dec 2011.[21] Simone Baldi, Tian Tao, and Elias B. Kosmatopoulos.Adaptive hybrid synchronisation in uncertain kuramotonetworks with limited information.
IET Control Theory& Applications , 13(9):1229–1238, Jun 2019.[22] Diego Paz´o and Ernest Montbri´o. Existence of hysteresisin the kuramoto model with bimodal frequency distribu-tions.
Physical Review E , 80(4), 2009.[23] David M´etivier and Shamik Gupta. Bifurcations inthe time-delayed kuramoto model of coupled oscilla-tors: Exact results.
Journal of Statistical Physics ,176(2):279–298, 2019.[24] B. C. Coutinho, A. V. Goltsev, S. N. Dorogovtsev, andJ. F. F. Mendes. Kuramoto model with frequency-degreecorrelations on complex networks.
Physical Review E ,87(3), Apr 2013.[25] M. K. Stephen Yeung and Steven H. Strogatz. Time delayin the kuramoto model of coupled oscillators.
PhysicalReview Letters , 82(3):648–651, Jan 1999.[26] Shadisadat Esmaeili, Darka Labavi´c, Michel Pleimling,and Hildegard Meyer-Ortmanns. Breaking of time-translation invariance in kuramoto dynamics with multi-ple time scales.
EPL (Europhysics Letters) , 118(4):40006,may 2017.[27] Darka Labavi´c and Hildegard Meyer-Ortmanns. Long-period clocks from short-period oscillators.
Chaos:An Interdisciplinary Journal of Nonlinear Science ,27(8):083103, Aug 2017.[28] Simona Olmi. Chimera states in coupled kuramoto oscil-lators with inertia.
Chaos: An Interdisciplinary Journalof Nonlinear Science , 25(12):123125, 2015.[29] Patrycja Jaros, Serhiy Brezetsky, Roman Levchenko,Dawid Dudkowski, Tomasz Kapitaniak, and YuriMaistrenko. Solitary states for coupled oscillators withinertia.
Chaos: An Interdisciplinary Journal of Nonlin-ear Science , 28(1):011103, Jan 2018.[30] Simona Olmi, Adrian Navas, Stefano Boccaletti, andAlessandro Torcini. Hysteretic transitions in the ku-ramoto model with inertia.
Physical Review E , 90(4),Jun 2014.[31] J. Barr´e and D. M´etivier. Bifurcations and singularitiesfor coupled oscillators with inertia and frustration.
Phys-ical Review Letters , 117(21), 2016.[32] Hisa-Aki Tanaka, Allan J. Lichtenberg, and ShinichiOishi. First order phase transition resulting from finiteinertia in coupled oscillator systems.
Physical ReviewLetters , 78(11):2104–2107, 1997.[33] Kurt Wiesenfeld, Pere Colet, and Steven Strogatz. Fre-quency locking in josephson arrays: Connection with thekuramoto model.
Physical Review E , 57(2):1563–1569, Neural Networks , 53:165–172,May 2014.[35] A. Mauro, F. Conti, F. Dodge, and R. Schor. Sub-threshold behavior and phenomenological impedance ofthe squid giant axon.
The Journal of General Physiology ,55(4):497–523, Jan 1970.[36] Christof Koch. Cable theory in neurons with active, lin-earized membranes.
Biological Cybernetics , 50(1):15–33,1984.[37] Kevin Dolan, Milan Majtanik, and Peter A Tass. Phaseresetting and transient desynchronization in networks ofglobally coupled phase oscillators with inertia.
PhysicaD: Nonlinear Phenomena , 211(1-2):128–138, 2005.[38] Milan Majtanik, Kevin Dolan, and Peter A Tass. Desyn-chronization in networks of globally coupled neurons withdendritic dynamics.
Journal of biological physics , 32(3-4):307–333, 2006.[39] E. Oh, K. Rho, H. Hong, and B. Kahng. Modular syn-chronization in complex networks.
Physical Review E ,72(4):047101, Oct 2005.[40] Jan Karbowski. Approximate invariance of metabolicenergy per synapse during development in mammalianbrains.
PLoS ONE , 7(3):e33425, Mar 2012.[41] K. Fonseca-Azevedo and S. Herculano-Houzel. Metabolicconstraint imposes tradeoff between body size and num-ber of brain neurons in human evolution.
Proceedings ofthe National Academy of Sciences , 109(45):18571–18576,Nov 2012.[42] Rachel Leander, Suzanne Lenhart, and Vladimir Pro-topopescu. Controlling synchrony in a network of ku-ramoto oscillators with time-varying coupling.
PhysicaD: Nonlinear Phenomena , 301–302:36–47, May 2015.[43] Spase Petkoski and Aneta Stefanovska. Kuramotomodel with time-varying parameters.
Physical ReviewE , 86(4):046212, Oct 2012.[44] M. A. Serrano, M. Boguna, and A. Vespignani. Extract-ing the multiscale backbone of complex weighted net-works.
Proceedings of the National Academy of Sciences ,106(16):6483–6488, Apr 2009.[45] Farzad V. Farahani, Waldemar Karwowski, and Nic-hole R. Lighthall. Application of graph theory for iden-tifying connectivity patterns in human brain networks:A systematic review.
Frontiers in Neuroscience , 13:585,Jun 2019.[46] Emily M. Jin, Michelle Girvan, and M. E. J. Newman.Structure of growing social networks.
Physical Review E ,64(4):046132, Sep 2001.[47] Akrati Saxena and S. R. S. Iyengar. Evolving models formeso-scale structures. In , page 1–8. IEEE, Jan 2016.[48] I. Leyva, A. Navas, I. Sendi˜na-Nadal, J. A. Almendral,J. M. Buld´u, M. Zanin, D. Papo, and S. Boccaletti.Explosive transitions to synchronization in networks ofphase oscillators.
Scientific Reports , 3(1):1281, Dec 2013.[49] I. Leyva, I. Sendi˜na-Nadal, J. A. Almendral, A. Navas,S. Olmi, and S. Boccaletti. Explosive synchroniza-tion in weighted complex networks.
Physical Review E ,88(4):042808, Oct 2013.[50] Thilo Gross and Bernd Blasius. Adaptive coevolution-ary networks: a review.
Journal of The Royal Society Interface , 5(20):259–271, Mar 2008.[51]
Adaptive networks: theory, models and applications . Un-derstanding complex systems. Springer, 2009.[52] Jun-Fang Zhu, Ming Zhao, Wenwu Yu, ChangsongZhou, and Bing-Hong Wang. Better synchronizabilityin generalized adaptive networks.
Physical Review E ,81(2):026201, Feb 2010.[53] Ritwik K. Niyogi and L. Q. English. Learning-rate-dependent clustering and self-development in a net-work of coupled phase oscillators.
Physical Review E ,80(6):066213, Dec 2009.[54] Jared C. Bronski, Yizhang He, Xinye Li, Yue Liu,Danielle Rae Sponseller, and Seth Wolbert. The stabilityof fixed points for a kuramoto model with hebbian inter-actions.
Chaos: An Interdisciplinary Journal of Nonlin-ear Science , 27(5):053110, May 2017.[55] Per Sebastian Skardal, Dane Taylor, and Juan G. Re-strepo. Complex macroscopic behavior in systems ofphase oscillators with adaptive coupling.
Physica D:Nonlinear Phenomena , 267:27–35, Jan 2014.[56] Evelyn Tang, C Giusti, G L Baum, S Gu, E Pollock,A E Kahn, D R Roalf, T M Moore, K Ruparel, R C Gur,R E Gur, Theodore D Satterthwaite, and Danielle S Bas-sett. Developmental increases in white matter networkcontrollability support a growing diversity of brain dy-namics.
Nat Commun , 8(1):1252, 2017.[57] Eli J Cornblath, E Tang, G L Baum, T M Moore, A Ade-bimpe, D R Roalf, R C Gur, R E Gur, F Pasqualetti, T DSatterthwaite, and Danielle S Bassett. Sex differences innetwork controllability as a predictor of executive func-tion in youth.
Neuroimage , 188:122–134, 2019.[58] Graham L Baum, Z Cui, D R Roalf, R Ciric, R F Bet-zel, B Larsen, M Cieslak, P A Cook, C H Xia, T MMoore, K Ruparel, D J Oathes, A F Alexander-Bloch,R T Shinohara, A Raznahan, R E Gur, R C Gur, D SBassett, and Theodore D Satterthwaite. Developmentof structure-function coupling in human brain networksduring youth.
Proc Natl Acad Sci U S A , 117(1):771–778,2020.[59] Budhachandra S. Khundrakpam, Andrew Reid, JensBrauer, Felix Carbonell, John Lewis, Stephanie Ameis,Sherif Karama, Junki Lee, Zhang Chen, Samir Das, andet al. Developmental changes in organization of struc-tural brain networks.
Cerebral Cortex , 23(9):2072–2085,Sep 2013.[60] Wei Gao, Sarael Alcauter, J. Keith Smith, John H.Gilmore, and Weili Lin. Development of human brain cor-tical network architecture during infancy.
Brain Struc-ture and Function , 220(2):1173–1186, Mar 2015.[61] Ginestra Bianconi.
Multilayer Networks: Structure andFunction . Oxford University Press, 2018.[62] Xiyun Zhang, Stefano Boccaletti, Shuguang Guan, andZonghua Liu. Explosive synchronization in adaptiveand multilayer networks.
Physical Review Letters ,114(3):038701, Jan 2015.[63] Raissa M. D’Souza, Jesus G´omez-Garde˜nes, Jan Nagler,and Alex Arenas. Explosive phenomena in complex net-works.
Advances in Physics , 68(3):123–223, 2019.[64] Vincenzo Nicosia, Per Sebastian Skardal, Alex Arenas,and Vito Latora. Collective phenomena emerging fromthe interactions between dynamical processes in multi-plex networks.
Physical Review Letters , 118(13):138302,Mar 2017.[65] Jordan D. Dworkin, Kristin A. Linn, Erin G. Teich, Perry Zurn, Russell T. Shinohara, and Danielle S. Bassett. Theextent and drivers of gender imbalance in neurosciencereference lists.
Nature Neuroscience , 2020.[66] Daniel Maliniak, Ryan Powers, and Barbara F Walter.The gender citation gap in international relations.
Inter-national Organization , 67(4):889–922, 2013.[67] Neven Caplar, Sandro Tacchella, and Simon Birrer.Quantitative evaluation of gender bias in astronomicalpublications from citation counts.
Nature Astronomy ,1(6):0141, 2017.[68] Paula Chakravartty, Rachel Kuo, Victoria Grubbs, andCharlton McIlwain.
Journalof Communication , 68(2):254–266, 2018. [69] Yannik Thiem, Kris F. Sealey, Amy E. Ferrer, Adriel M.Trott, and Rebecca Kennison. Just Ideas? The Statusand Future of Publication Ethics in Philosophy: A WhitePaper. Technical report, 2018.[70] Michelle L Dion, Jane Lawrence Sumner, andSara McLaughlin Mitchell. Gendered citation patternsacross political science and social science methodologyfields.