Development of structural correlations and synchronization from adaptive rewiring in networks of Kuramoto oscillators

Synchronization of non-identical oscillators coupled through complex networks is an important example of collective behavior, and it is interesting to ask how the structural organization of network interactions influences this process. Several studies have explored and uncovered optimal topologies for synchronization by making purposeful alterations to a network. On the other hand, the connectivity patterns of many natural systems are often not static, but are rather modulated over time according to their dynamics. However, this co-evolution and the extent to which the dynamics of the individual units can shape the organization of the network itself are less well understood. Here, we study initially randomly connected but locally adaptive networks of Kuramoto oscillators. In particular, the system employs a co-evolutionary rewiring strategy that depends only on the instantaneous, pairwise phase differences of neighboring oscillators, and that conserves the total number of edges, allowing the effects of local reorganization to be isolated. We find that a simple rule—which preserves connections between more out-of-phase oscillators while rewiring connections between more in-phase oscillators—can cause initially disordered networks to organize into more structured topologies that support enhanced synchronization dynamics. We examine how this process unfolds over time, finding a dependence on the intrinsic frequencies of the oscillators, the global coupling, and the network density, in terms of how the adaptive mechanism reorganizes the network and influences the dynamics. Importantly, for large enough coupling and after sufficient adaptation, the resulting networks exhibit interesting characteristics, including degree–frequency and frequency–neighbor frequency correlations. These properties have previously been associated with optimal synchronization or explosive transitions in which the networks were constructed using global information. On the contrary, by considering a time-dependent interplay between structure and dynamics, this work offers a mechanism through which emergent phenomena and organization can arise in complex systems utilizing local rules.

Electrical activity in populations of neurons, or the spread of information or disease across social networks are common examples of dynamical processes unfolding on networks.In many real-world systems, synchronization of dynamics, and its dependence on network architecture, is of critical interest.Importantly, network interactions can also be adaptive, such that connections rearrange over time according to the dynamical states of the nodes.A prevalent example of adaptation in connectivity takes place in biological neuronal networks, which reconfigure via local plasticity mechanisms that are informed by the dynamical relationships between directly coupled elements.Motivated by these types of systems, here we address the question of how both structured network topology and global synchronization can develop through a co-evolution of the underlying network connectivity and the dynamics.
Using the canonical Kuramoto model, we suggest a simple, adaptive rewiring scheme based on local phase information that can evolve initially unstructured random graphs towards topologies with specific organizational principles.In turn, we find that this process simultaneously enhances synchronization.The structure that arises in the co-evolving systems manifests as specific relationships between the natural frequencies of the oscillators and the network topology, properties that have previously been associated with optimal synchronization in which networks were purposefully constructed using complete information of the network and node properties.Yet here, without any global information, we consider a timedependent interplay between structure and dynamics, and offer a mechanism through which organization can emerge in complex systems utilizing local rules.
One of the most common and useful models for studying synchronization is the canonical Kuramoto model [42,43], which originally described the evolution of a population of N all-to-all coupled phase oscillators that were in general non-identical (see [44,45] for reviews).In recent years, this model has been extended to study systems with heterogeneous network topologies, in order to investigate how the architecture of complex connectivity affects the onset of synchronization in diverse oscillator populations [24,46].Such efforts have provided important insights into the nature of the synchronization transition in different graph models including those that display a scale-free degree distribution [47][48][49][50], those with small-world architecture [51], and those with community structure [52][53][54][55][56][57].
Within this body of work, particular attention has been paid to understanding what features of a network inhibit or enhance the ability to support collective dynamics.For the case of identical oscillators, this is often studied from the perspective of minimizing a ratio of eigenvalues that depend only on the structure of the network [58,59].However, the question of "optimal" networks for synchronization can be more interesting and complex when the oscillators' natural frequencies are heterogeneous, a characteristic of many real-world systems.In particular, for non-identical oscillators, a crucial consideration becomes how the structure of the network is intertwined with dynamical properties.For example, in the Kuramoto model, synchronization can be enhanced when there are specific types of correlations between node degrees and oscillator frequencies or between the natural frequencies of adjacent oscillators [60,61].Gómez-Gardeñes et al. [62] demonstrated that in scale-free networks, positive frequency-degree correlations can lead to a first-order, or explosive, transition to synchronization.More recently, discontinuous transitions have been found by imposing constraints on the minimal difference between connected nodes' natural frequencies [63,64].There has also been progress in analytical work determining network topologies that enhance synchronization.For example, it has been shown that optimal networks for synchronizing collections of non-identical oscillators exhibit particular relationships between Laplacian eigenvectors and oscillator frequencies [65][66][67].In addition, dimensionality reduction approaches [68] have recovered many previous numerical results, and have been used to derive analytical conditions for optimizing synchronization of Kuramoto oscillators in networks with attractive and repulsive interactions [69].
Notably, while the coupling structure between oscillators can drive their dynamics, network dynamics can also modulate structure.Specifically, in adaptive systems, the pattern of connectivity itself is continuously updated and modified in response to the dynamics that occur on top of it [70][71][72].Systems that display these processes can be observed across biological, ecological, social, and distribution networks [70,71], and collectively they can be characterized by topology and dynamical states that coevolve with one another.Kuramoto-like models in particular provide a useful framework in which to explore the effects of co-evolution and adaptation [73][74][75][76][77][78][79][80][81][82][83][84][85], and allow one to address questions such as (i) can and how might these systems organize themselves towards network configurations that enhance local or global synchronization?, and (ii) from an initially unstructured topology, can and how do different adaptive mechanisms lead to the emergence of certain architectural patterns or correlations between dynamical properties and network structure?
In addition to being adaptive, it is important to note that the evolution of many real-world networks is often governed by local rules, in which node dynamics update as a function of only neighboring node states, and in turn, the placement or weights of edges update primarily as a function of the states of the nodes they directly couple [70][71][72].This type of behavior is especially pertinent in biological systems, which typically co-evolve in the absence of global or top-down controllers of node states and/or network structure.A particularly salient example of this occurs in biological neuronal networks, where, under Hebbian plasticity rules, increases in synaptic weights occur when connected neurons exhibit correlated dynamics [86][87][88], while under anti-Hebbian plasticity rules, the opposite occurs [89,90].Other systems that obey such local adaptation are prevalent, and studied examples include models of reconfiguration of social networks under disease propagation [91] and opinion formation [92], or reorganization under feedback mechanisms in boolean models with applications to gene regulatory or neural networks [93,94].
In this study, we use the Kuramoto model to investigate how a simple, adaptive rewiring scheme can evolve initially unstructured random graphs towards ordered topologies, and also simultaneously lead to enhanced synchronization in the system.Importantly, the rule is in-formed by only local information of neighboring nodes' states at a given instant of time, and works by regularly breaking and randomly rewiring connections between more instantaneously phase-synchronized oscillators, while maintaining connections between more desynchronized pairs of oscillators.This process repeats continually over time, and can be thought of as a repulsive mechanism, or one that tends to represses assortativity (in terms of nodes connecting to other nodes with similar instantaneous states).We find that co-evolution of the network and dynamics can promote the degree of synchronization in the system, which occurs in tandem with the development of specific correlations between the topology and the natural frequencies of the oscillators.In previous work, these features have been imposed by purposeful selection, or have been shown to arise in work on optimizing synchronization.Here, however, the properties emerge from the interplay between network structure and dynamics.Focusing on the simplest situation of binary, undirected networks, we isolate the effects of adaptive reconfiguration alone, and thereby uncover a process through which heightened collective dynamics and organized network structure simultaneously arise in a local, unsupervised way.
The remainder of this paper is organized as follows.Sec.II states the formulation of Kuramoto dynamics on complex networks.In Sec.III we first briefly outline past work to motivate the specific mechanism studied here, and then detail the proposed co-evolutionary strategy.Sec.IV describes several interesting results of the adaptive process, and in Sec.V we discuss the implications of our findings and conclude.

II. THE KURAMOTO MODEL ON COMPLEX NETWORKS
Of the many models that exist for studying synchronization phenomena on complex networks, one of the most useful has been the paradigmatic Kuramoto model [42,43].It describes the dynamical evolution of a population of N phase oscillators coupled on a network according to the following equation: In this formulation, θ i is the instantaneous phase of the i th oscillator, ω i is its natural frequency, α is the overall coupling strength, and A is the N × N adjacency matrix describing the connectivity of the network.In this report, we consider binary, undirected networks, such that there is an edge between nodes i and j, 0 otherwise.
The natural frequencies {ω i } are distributed according to a probability density g(ω), which we will take to be symmetric and centered around a mean frequency of zero.
The overall amount of synchrony in the population at a given time t is typically quantified with the Kuramoto order parameter [42,43] which can be thought of as the centroid of the N phases on a unit circle in the complex plane [46].Here, ψ is the average phase of the population, and the modulus R, given by, quantifies the amount of phase coherence.When the oscillators' phases are uniformly spread, R ≈ 0 and the system exhibits low synchrony.On the other hand, when the phases become tightly clustered, R ≈ 1 and the system exhibits high levels of synchrony.
We can also use this order parameter, which ranges from R = 0 (complete incoherence) to R = 1 (complete phase synchronization), to monitor the global degree of synchrony in the system as a function of the coupling α.
In this case, one typically reports a time-averaged value computed on an interval of length T after several transient or relaxation time steps T R have been discarded.

III. MOTIVATION AND THE CO-EVOLUTIONARY MODEL
In this study, we consider a feedback process between dynamics and the restructuring of network topology that integrates different ideas and results from previous studies on adaptation or enhancing synchronization in networks of non-identical Kuramoto oscillators.To better frame and motivate our contributions, we briefly outline some past work on these topics below.

A. Inspiration from prior investigations
In general, it is expected that different adaptive strategies will lead to the emergence of different patterns in both the network topology and the dynamics, and a number of studies have explored these ideas using, for example, chaotic dynamics [74-76, 95, 96], models of neuronal dynamics [79,97,98], and non-identical Kuramoto oscillators [73,77,78,[80][81][82][83][84][85][99][100][101][102][103][104].Focusing on the latter of these classes, one early study found that dynamical rewiring to force links between nodes with more similar time-averaged frequencies creates strongly synchronized clusters of nodes, and the network reaches a smallworld configuration [80].More recently local, competitive adaptation mechanisms, which tend to strengthen (weaken) connections between more dynamically coherent (incoherent) oscillators have been shown to lead to the emergence of modular organization in Kuramoto networks [73,84,85,101], and positive feedback can simultaneously enhance synchronization and percolation in initially fragmented networks [105].In complementary efforts, Sendiña Nadal et al. [82] and Sendiña-Nadal et al.
[83] studied a growth process in which heterogeneous oscillators make connections to external pace-maker nodes so as to become locked with the pace-maker dynamics.When the attachment process is preferential and determined from differences in the dynamical states of the heterogeneous oscillators and pace-maker nodes, entrainment can occur simultaneously with the emergence of a power law degree distribution.In addition, adaptive processes that favor the strengthening of edges between more out of phase oscillators, can significantly improve global synchronization in non-identical Kuramoto networks [99,100].However, the resulting networks were not necessarily evolved under fixed total weight and the topology was not analyzed in depth for relationships between structure and dynamical properties as a function of the global coupling, both of which are aims of the present work.
Another line of inquiry revolves around the problem of optimizing synchronization of non-identical systems of Kuramoto oscillators in order to understand what structural properties of the network are important.For example, using an optimization procedure to maximize the order parameter R, Brede [60,61] uncovered key features that can enhance synchronization.These include the placement of more edges on oscillators with natural frequencies further from the mean (yielding positive correlations between degrees and frequency magnitudes), as well as the preferential attachment of oscillators with positive frequencies to other oscillators with negative frequencies (and vice-versa, yielding negative correlations between the intrinsic frequencies of adjacent oscillators.These findings have been corroborated in several other studies on optimizing networks of non-identical Kuramoto oscillators as well [106][107][108][109].Furthermore, forcing positive versus negative correlations between adjacent frequencies appears to change the critical coupling and exponents for the synchronization transition [110].Other work has considered the enhancement of both global and local phase synchronization [111,112], finding that local synchronization leads to an onset of collective dynamics at lower couplings -facilitated by clustering and the grouping together of nodes with similar frequencies -but makes the state of full synchronization more difficult to achieve.More recently, spectral analyses have been employed to show that synchronization is optimized under specific couplings of the natural frequencies to eigenvectors of the Laplacian, i.e., when the frequencies are maximally aligned with the largest eigenvector [65][66][67].Importantly, the previously reviewed works indeed uncover several crucial network features for optimizing collective dynamics.However, the question of if or via what type of mechanism such networks could be generated or evolved for using solely local adaptive strategies -and how such a process occurs over time and at different couplingsremains open. In this study, we numerically investigate the interplay between network structure and oscillator dynamics.We report on a local, state-dependent rewiring mechanism that can (i) evolve initially random and uncorrelated networks towards structured configurations with specific relationships between dynamical and topological properties and (ii) through a reciprocal process, simultaneously improve synchronization.While the ideas of dynamical self-organization from local rules and global network optimization strategies have been studied on separate fronts, here, we attempt to specifically consider them in tandem.
In what follows, we first state the initial setup and parameters of the system (Sec.III B) and then describe the co-evolutionary process in detail (Sec.III C).

B. Initial network construction
All simulations were carried out with a 4 th order Runge-Kutta method using a time step of ∆t = 0.02.Initial phases {θ i (0)} were distributed at random in the interval [−π, π], and the natural frequencies {ω i } were drawn at random from a uniform distribution in the range [−2, 2], denoted as {ω U }.In Appendix B), we also show results for the case of a Gaussian distribution with zero mean and unit standard deviation, denoted as {ω G }.
The initial network configurations are binary and undirected Erdös-Renyi (ER) random graphs of type G(N, L) (i.e., the network is drawn from the distribution of random graphs with N nodes and L edges), with corresponding average degree k = 2L/N .Importantly, this initialization produces networks without special topological characteristics and without relationships between network properties and dynamical properties.We will denote the initial network configurations as G o and the networks at the end of the adaptive process as G , and similarly we will use 'o' and ' ' to denote quantities computed on each network.For the remainder of the main text, we fix N = 100 and examine two different values of L chosen such that k = 12.5 or k = 25.(In the S.M. [113], we also explore the robustness of several results with respect to variations in system size, the network density, the initial network topology, and the presence of asymmetry in the frequency distribution).Reported measures correspond to ensemble averages over independent simulations using different instantiations of the initial network, initial phases, and node frequencies.FIG. 1.A schematic of the network update process.(a) At a time tm = TR + mT , node i is selected at random from all nodes in the network, and it breaks its connection to the node k, which (b) is the node it is instantaneously most in phase with.(c) A new edge (i, q) is then created between i and a randomly selected node q with phase θq

C. Mechanism of adaptive rewiring
We turn now to an explanation of the co-evolutionary scheme, whereby the network is restructured according to the dynamics.Motivated by previous literature [99,100], we take as our starting point a mechanism that aims to maintain connectivity between a given node and its neighbors with which it is more instantaneously out of phase, and rewire connections between a given node and its neighbor that it is most instantaneously synchronized with.After setting up the initial conditions and network G o as described in Sec.III B, the dynamics are first run for a relaxation interval T R .Then at regular times t m = T R + mT , where m is the number of attempted rewirings and T is an associated interval characterizing the adaptation process (which will be some multiple of the time step), rewiring is invoked as follows.A node i is chosen at random from the network and the quantity is computed for all j ∈ N i , where N i denotes the set of nodes directly connected to i.Note that this function is local in the sense that it depends only on the instantaneous phases of node i and of the nodes j ∈ N i that are neighbors of i at the current time t m .It takes on a value of 1 when θ i −θ j = ±π (maximal phase separation) and a value of 0 when θ i − θ j mod 2π = 0 (perfectly in-phase).The edge (i, k), where k is the node in N i minimizing f ik , is then broken and a new link (i, q) is formed between node i and a randomly selected node q that was not one of i s neighbors.This step corresponds to the breaking and rewiring of the link between node i and the neighbor k with which it is currently most in phase (Fig. 1).These dynamical update rules are repeated m times, resulting in an evolved network G .We use m = 2.5 × 10 4 , which allows us to observe interesting dynamical and structural changes in the system, while being within enough computational reason to allow for the exploration of several different network and natural frequency parameter changes (however, we also point out places where the capability to run more adaptation steps may offer addition insight).Finally, once rewiring has ceased, the dynamics continue to run atop G for another relaxation interval.For clarity, in the figures that follow, all quantities that change every time step will be labeled as "time, t", and all quantities that change only when the network is rewired will be labeled as "rewiring step, m".
Before continuing, we wish to point out some features of this mechanism.First, it applies to the case of binary, undirected connectivity, and maintains the density of the network throughout co-evolution.Setting these constraints disambiguates the role of rearrangements in network topology in shaping results from other factors such as increases in the total number of edges, heterogeneous weighting, or directionality.Furthermore, it allows for the cleanest comparison against much of the previous work on optimizing synchronization via rewiring.Second, we note that this adaptive process is a type of repulsive or suppressive strategy, but unlike [74-78, 99, 100], which consider weight plasticity, or [82,83] which consider a growth mechanism, we consider connectivity reorganization where edges in an initially random configuration are periodically pruned and rewired between the most instantaneously and locally in-phase oscillators, but remain in place between the more locally dissonant oscillators.The random rewiring after edge deletion introduces realistic stochasticity and allows for a sampling of the network, without breaking the locality condition for determining which edge is rewired [71].Though not a model of a specific system, this type of disassortative mechanism has real-world counterparts.For example, it is in the same vein as anti-Hebbian learning rules in neural systems, where synapses weaken between more dynamically correlated neurons and strengthen between more incoherent neurons [89,90].It may also mimic some strategies of opinion formation and influence on social networks, where people preferentially link with those of more different opinions from themselves [114].

A. Co-evolved networks exhibit enhanced global synchronization
Our first main result is that the adaptive rewiring mechanism is able to enhance global synchronization over a broad range of coupling values.This result holds for both frequency distributions (See Appendix B for the analysis with the normally distributed frequencies) and values of average degree.In addition, the result is quite robust to an order of magnitude difference in the structural reconfiguration interval, T .
We first discuss the time evolution of the order parameter R(t) for the co-evolving networks.Fig. 2 depicts examples of R(t) vs. time t for the case of uniformly distributed frequencies {ω U } and T = 0.2.The top row corresponds to networks with k = 25, at representative couplings (a) α = 0.065 and (b) α = 0.085, and the bottom row corresponds to networks with k = 12.5, at  In each case, the dynamics are first run atop the initially static ER network for several time steps.Adaptation begins at the time denoted by the first red line, and ends after several time steps at the second red line.We observe that at lower coupling values -where the dynamics on the initial network exhibit little coherence across time -the adaptive strategy is able to significantly increase R to an intermediate value during the rewiring stage, though the order parameter may still exhibit fluctuations.As α is increased, the order parameter in the non-adaptive regime sits at an intermediate average value, but once coevolution begins, the self-organizing network rearranges such that R again increases and reaches a value near 1.When rewiring ceases after several cycles, the resulting networks are able to maintain these states of heightened collective dynamics, and the global order parameter remains at an increased value from its initial location.However, though the time-averaged value R remains high for the networks with k = 12.5, we find that in some cases the order parameter still exhibits fluctuations, even after the adaptation period has ceased.This intuitively suggests that the co-evolved networks with higher average degree are more robust to the stochasticity in the rewiring and the exact placement of edges in the network, in terms of their ability to support a smooth, frequencysynchronized steady-state.Appendix B 1 contains additional figures of R(t) vs. time for the case of normally distributed frequencies.Also, in the S.M., we examine another measure of synchronization that quantifies how the number of locally-synchronized clusters [53] of oscillators changes as a function of time due to the network rearranging [113].
In order to obtain a more complete picture of the effect of adaptive rewiring, we performed a sweep over a comprehensive coupling range.For k = 12.5 we considered a range α ∈ [0, 0.4], and for k = 25, we considered α ∈ [0, 0.2]; in both cases, couplings were sampled at a resolution ∆α = 0.005.At each value of α, networks and initial conditions were initialized as described in Sec.III B. We then ran a set of simulations on the original, uncorrelated networks G o , and obtained a timeaveraged value of the global order parameter R (Eq. 5) over the last 1 × 10 4 time steps.A second set of simulations were then run with the same initial conditions, but under the co-evolutionary scheme (i.e. the network topology was allowed to co-evolve with the dynamics for m rewiring steps, after which a time-averaged order parameter was computed on the final adapted network G ).In each panel, the gray data points correspond to static ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the uniform distribution {ωU }, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.
We show the outcome of this analysis in Fig. 3, where each panel is a plot of R vs. α.Panels (a) and (b) correspond to networks with k = 12.5 and k = 25, respectively.Furthermore, we plot curves for two different values of the waiting time between structural changes to the network: T = 0.2 and T = 2.As can be seen from each of the curves, there is a broad range of α over which the rewiring mechanism for network restructuring leads to improvements (higher values) of the time-averaged global order parameter compared to the random networks.This enhancement does not begin immediately at α = 0, but once the coupling is high enough (still at a relatively low value), the order parameter for the adaptive networks G begins to increase at a much steeper rate.Once this begins, we find that R remains noticeably higher for the adaptive networks compared to the static networks G o , across all couplings beyond a certain point.These trends are robust for both rewiring intervals and average degree values, and as can been seen in Sec.B 1, the results hold for the case of normally distributed natural frequencies as well.(The S.M. [113] shows additional and qualitatively similar findings for simulations on slightly larger networks or those with lower mean degree, and for the case of a non-symmetric frequency distribution).It is also worth noting that -although the total number of times the network is allowed to rewire is limited by computational constraints -since there is no builtin condition for adaptation to cease and since there is a stochastic component to the network reconfiguration, the ability to run more adaptation steps may yield even further improved results.

B. Emerging structure and correlations between network topology and dynamical properties
Given that dynamical reconfiguration of the networkinformed by local information on the states of connected oscillators -can improve the overall amount of synchronization in the system, a second line of inquiry is understanding what properties of the evolved networks lend themselves to this capability.In particular, does the self-enacted rewiring mechanism generate the interesting topological features and correlations that arise with global optimization schemes, which in a cyclic process, then allow synchronization to occur?
We begin by investigating the networks G for certain relationships between their topology and the natural frequencies of the oscillators.Two properties in particular -degree-frequency correlations and frequencyneighbor frequency correlations -have been associated with networks optimized for synchronization [60,61,65] (see Sec.I and Sec.III).
Though the natural frequencies of each node remain fixed here, the way oscillators of different frequencies are coupled to one another changes over time from the initial, random configuration to the end of adaptation.In Fig. 4, we show examples of how two specific relationships manifest in an ER network G o and the corresponding adaptively rewired network G .The top and bottom two rows of Fig. 4 correspond to networks with k = 12.5 and k = 25, respectively, and in both cases, the first column corresponds to the ER network and the second column corresponds to the rewired network.For each node i, we first computed the offset ωi of i's intrinsic frequency from the mean of the population, where ωi = ω i − ω and ω is the mean intrinsic frequency over all nodes in the system.Panels (a,b) and (g,h) then show node degree k i vs. frequency offset ωi , and panels (d,e) and (j,k) show the average frequency offset of oscillator i s neighbors, ω Ni = j∈Ni ωj /k i , vs. frequency offset ωi .Note that the coupling values α are such that the original networks exhibited intermediate levels of synchrony, and the reor-ganized networks were able to entrain the population to a higher level of synchrony.
As expected, initially there is little correlation present between the topology of the network and the dynamical property of oscillator frequency, as illustrated by the lack of organization in the plots in the first column of Fig. 4.However, from observation of the second column, it is evident that when the co-evolutionary mechanism can enhance synchronization, it simultaneously leads to the emergence of very specific correlations between the network connectivity and the intrinsic frequencies of the oscillators.We first note the appearance of the marked vshaped curves characterizing the plots of k i vs. ωi (panels (b) and (h)), which signify that the node degree becomes positively correlated with the absolute value of the oscillator frequency offset.In other words, nodes with natural frequencies further from the mean natural frequency of the population gather proportionately more edges.
A second finding is that when increased global synchrony arises from a restructuring of the network, the final arrangement exhibits distinct relationships between the frequency offsets of a given oscillator and the frequency offsets of that central oscillator's direct neighbors on the network.The patterns in panel (e) and (k) showing ω Ni vs. ωi point to the fact that oscillators with positive natural frequency offsets tend to become connected to other oscillators with, on average, negative natural frequency offsets.Since the mean frequency ω ≈ 0 for the distributions we consider, this implies that oscillators with positive frequencies tend to become neighbored by oscillators with negative intrinsic frequencies, and vice-versa.In Sec.A we also define an additional measure of frequency-neighbor frequency organization as the correlation C ω, ω between natural frequency offsets and the sum of neighbor frequency offsets.Fig. 11 shows examples of this quantity for the same networks as those in Fig. 4.
To quantify these relationships and study how they evolve with the number of rewiring steps, we considered two summary statistics, following [60,61,110,111].The first of these measures is a simple (Pearson) correlation coefficient, C |ω|,k , to quantify the strength of the relationship between node degree and the magnitude of frequency offset.This measure increases steadily throughout co-evolution of the network and dynamics (Fig. 4c,i ).In addition, for each node i, we calculated the fraction of its neighbors f i that had natural frequency offsets of the opposite sign as compared to oscillator i s frequency offset, and then computed an average f = i f i /N over all nodes in the network.This metric also increases as the network is rewired, as observed in Fig. 4f,l.(See Fig. 11 for an example of the evolution of the additional measure,C ω, ω ).
We have thus far shown examples of emerging structural patterns that arise when a co-evolved network is clearly able to entrain the oscillators to a state of higher synchrony.However, the ability of this behavior to occur is also dependent on the coupling α.In order to better  understand the appearance of the topological and dynamical correlations and their dependence on the overall coupling, we computed the same measures (C |ω|,k and f ) as a function of α (see Figs. 5, 6, respectively).In each case, the measures were computed on the final networks G that exist at the end of the adaptation period.
For each combination of natural frequency distribution, average degree, and rewiring interval, we observe robust trends.In particular, C |ω|,k remains near zero at low coupling values, and then proceeds to quickly in-crease as a function of α until it saturates to a relatively steady value close to 1.In other words, as the overall coupling increases, the correlation between the node degrees and the magnitude of the frequency offsets becomes more apparent.The frequency-neighbor frequency relationship follows a similar evolution.Specifically, the mean fraction f (averaged over all nodes in the network) of an oscillator's neighbors that have frequency offsets of opposite sign relative to that of the central oscillator also grows with α, plateauing near a high value of  f ≈ 0.9.This points to a heightened mixing of oscillators with different intrinsic properties that arises with increased coupling.We also refer the reader to Fig. 12 for the analysis of the supplementary quantity, C ω, ω , which supports and further aids in the understanding of these findings.
There are three main features in these trends worth pointing out explicitly.First, for the rewiring intervals considered here, the general form of the curves appears to be relatively independent of the time T between a structural perturbation to the network.Second, the two quantities considered, C |ω|,k and f , both exhibit a clear change in behavior (signified by the onset of a rapid increase in value) at approximately the same coupling.Thus, in the proposed co-evolutionary scheme, the emergence of these relationships between the network topology and intrinsic properties of the dynamics seem to arise in tandem to one another, in that the appearance of one feature implies the development of the other.A final important observation is that the coupling at which these patterns begin to take shape is near the coupling at which the global order parameter begins to rise (compare to Fig. 3).(Fig. 12 shows that consistent behavior is found for C ω, ω as well).Thus, as found by Brede [60,61,110,111] in his work on optimization of synchronization of non-identical oscillators, enhanced collective behavior co-occurs with the materialization of specific relationships between network connectivity and the intrinsic frequencies.Here, we have shown that organized structure in the form of these correlations can emerge in networks from a simple, adaptive mechanism based on a combination of local state information and stochastic rewiring.We also again note that allowing the co-evolution to run longer may give heightened results, especially at low to intermediate couplings.

C. Analysis of the time-dependence of the adaptive mechanism
Our analysis thus far has considered the global order parameter and the correlations between the intrinsic frequencies and the network topology that arise in the adaptively rewired networks, and we have mainly focused on how these properties manifest after several rewirings of the network and how their strength depends on the overall coupling.While Figs. 2 and 4 do briefly examine the time-dependence of these properties as well, in order to gain a further understanding of the adaptation mechanism, we carry out a more in-depth exploration of how the co-evolutionary process unfolds over time at different couplings, and also how individual oscillators with different intrinsic dynamical properties (i.e.intrinsic frequencies) are affected by the rewiring mechanism.

Evolution of the instantaneous frequencies
To investigate the temporal evolution of the systemand how the adaptive scheme works from a more local standpoint -we first examine the instantaneous frequencies θi (t) as a function of time.Recall that the condition for a frequency-synchronized state corresponds to the instantaneous frequencies of all oscillators locking to the mean natural frequency of the population, ω .Thus, to better understand the mechanism of the co-evolution process, we specifically consider how oscillators of different intrinsic frequencies (in terms of how close or far ω i is to the average natural frequency ω ) evolve as a function of time, and how they may be differentially affected At very low coupling, adaptively rewiring the network appears to have negligible effect on the instantaneous frequencies.But as the coupling increases to an intermediate value near the transition to synchronization (see Fig. 3), we observe a change in the dynamics: oscillators with intrinsic frequencies near the center of the distribution start evolving with a frequency near the mean earliest, and then later in time, more outlying oscillators become incorporated into the coherent group.For the denser networks in particular, there are a few sampled coupling values for which the co-evolution results in partial synchronization of those oscillators around the center of the natural frequency distribution, but that coherent core does not extend to the most disparate oscillators.Thus, there does seem to be a dependence on the natural frequencies in terms of how the co-evolution develops and affects different oscillators over time.However, at even larger coupling, the transition in the dynamics becomes much faster, and the dependence on the intrinsic frequencies is less noticeable.The findings are similar for the case of normally distributed frequencies (Appendix B 3, Fig. 19) Before continuing, we note that the figures shown here (and throughout the rest of Sec.IV C) are for single realizations of the initial network topology and the intrinsic frequencies.For these particular instantiations, we have sampled coupling values that highlight different regimes -in regards to the behavior over time and the overall outcome -of the adaptive rewiring.However, it is important to state that over different realizations, we observe some variability in terms of how the co-evolution affects the dynamics and the network over time, and whether or not it is able to cause significant changes in the dynamics and the network structure at a given coupling.This is especially apparent for the sparser networks and low values of α.Though it is beyond the scope of the current work, it may be interesting in future work to investigate the degree of this variability and its dependence on the frequency distribution and properties of the network such as density.

Evolution of the correlations between topology and natural frequencies
We know from Sec. IV B that increases in the order parameter due to the adaptive mechanism are accompanied by the emergence of correlations between the network topology and the oscillator frequencies.Therefore, we now further examine how the network organization restructures over time as the system co-evolves.We focus in particular on the evolution of the degree of individual nodes across time, and also consider how the degree -natural frequency relationship proceeds as the network rearranges itself at different values of the global coupling.
For the same initial networks, natural frequencies, and coupling values, Fig. 8 shows the node degree k i vs. the rewiring step m and Fig. 9 shows the correlation C |ω|,k The natural frequencies were uniformly distributed in both cases.Below, we discuss the observations for each of the densities in turn.
For k = 12.5 and at low coupling (i.e.α = 0.04), the rewiring does not noticeably affect how edges are distributed on particular oscillators, and C |ω|,k fluctuates around zero.As the coupling increases, though, the coevolution begins to cause significant changes in the network (see panels for α = 0.114 and α = 0.0115).In particular, there is a short period of time in which edges become concentrated on oscillators with natural frequencies near the mean, whereas oscillators with intrinsic frequencies on the ends of the distribution remain with relatively low degrees.This is followed by the gradual spreading out of edges onto the more outlying oscillators with subsequent rewiring.We can quantify this behavior by considering the value of the correlation C |ω|,k .For α = 0.114 and α = 0.115 in this example, we see that a slight negative degree-frequency correlation develops briefly in the initial stages of the rewiring, followed by an increase in this quantity to a high, positive value.
For the case of k = 25, we again see that at very low coupling, the rewiring does not drive persistent changes in the system.Interestingly, though, for slightly larger coupling (i.e.α = 0.04 and α = 0.05 in this example) there is a regime during which edges consistently build up on oscillators with natural frequencies near the center of the distribution.This results in C |ω|,k becoming significantly negative due to the co-evolution of the network and dynamics, and unlike the situation in the sparser networks, this behavior can persist across the entire rewiring period (see, for example, α = 0.04).The concentration of edges on oscillators near the center of the distribution explains the dips in Fig. 5b, which shows C |ω|,k vs. α for the case of k = 25.We also note that near the value of the coupling where this dip occurs, the order parameter for the rewired networks begins to deviate in a positive direction from that of the ER networks (Fig. 3b), which is likely due to the formation of a locally synchronized cluster of oscillators with natural frequencies near the population average.We do not observe this type of organization for the sparser networks, suggesting an intricate dependence on the network density; this may be an interesting parameter to explore further in future work.As the coupling increases further to intermediate values (for example, α = 0.06 and α = 0.065), we find that there is first an increase in degree for oscillators with frequencies more closely surrounding the mean, which gives rise to a slightly negative degree-frequency correlation.But as the network continues to co-evolve, edges begin to localize on more outlying oscillators, and the degree-frequency correlation crosses zero and then starts to become positive.However, the most disparate oscillators may still not be able to gather enough edges, resulting in a positive correlation that is less than 1 (i.e., there is some scatter in the relationship).
As the coupling increases further, there is another shift in terms of how the adaptive mechanism affects the networks.For the example with k = 12.5, we observe some fluctuation in C |ω|,k when the rewiring begins, but eventually the period of marked negative degree-frequency correlation disappears and the edges rapidly become redistributed onto oscillators with the most disparate intrinsic frequencies.We also find that for the examples with k = 25, C |ω|,k almost immediately begins to rise -rather than decreasing first -at high coupling.In addition, the correlation in each case quickly saturates at a value close to 1, signifying a very strong degree-frequency relationship that extends to even the oscillators on the far edges of the natural frequency distribution.Thus, at large couplings the time evolution is similar for k = 12.5 and k = 25.If network co-evolution were allowed to increase for an even longer time, we may observe a plateau in C |ω|,k for more intermediate values of α, though not necessarily at a level corresponding to a near-perfect positive relationship.We repeat this analysis for the case of normally distributed frequencies {ω G } in Sec.B 3 of the Appendix and find qualitatively similar results.
Although a rigorous mathematical treatment of the coevolution is beyond the scope of this study, it is useful to postulate mechanisms that might be able to explain and complement at least some of the empirically-based descriptions discussed above.For example, one interesting observation is that, at low coupling and in the initial stages of rewiring, oscillators with more outlying natural frequencies tend to have lower degree than some oscillators with natural frequencies closer to the mean of the distribution.This seems to occur, at least to some extent, across both values of the network density and for both frequency distributions.In order to think about this phenomena, we first remark that changes in the structure of the network (and thus the dynamics) must be due to which edges are broken over time, since the location of a new edge is chosen at random.Second, oscillators with more disparate intrinsic frequencies will naturally want to rotate more rapidly than oscillators with more moderate intrinsic frequencies.Following this reasoning, we can then posit that the most outlying oscillators will have a greater chance of being most in phase with a focal node -which is selected at random -than oscillators with moderate intrinsic frequencies whose phases will tend to change less quickly and thus reduce the chance of being nearby the focal node.Therefore, this thought process suggests that when the selected focal node assesses its phase difference with its neighbors, the oscillators with natural frequencies most different from the mean will have a greater probability of becoming disconnected, and this would account for the observed lower degree of these nodes at low coupling and at the start of the adaptation period.Preliminary results show that making T smaller can sometimes slightly enhance or extend the region of couplings over which C |ω|,k goes negative, which is consistent with the proposed logic.Another observation, though, is that the degree-frequency correlation can change sign over time (go from slightly negative to positive), which does not seem directly obvious to explain from the previously outlined arguments.In order to make wholly concrete statements and to understand in more detail how the adaptation process gives rise to various results -such as the intricacies of the time-evolution of the network and also the dependence of that timeevolution on parameters like the global coupling, network density, and frequency distribution -will require a much more in-depth investigation.While outside the main contributions of this paper, it would be useful to formulate and carry out a more formal theoretical analysis in future work.

D. Spectral Analysis
As an additional connection to recent literature on optimizing synchronization in heterogeneous oscillator populations, we also carry out a spectral analysis, following the work of Skardal et al. [65].In a series of papers [65][66][67], the authors derived conditions that describe the deviation from full phase locking, and also uncovered conditions for promoting synchronization.By considering a linearized form of the dynamics valid in the highsynchrony regime, one of the main analytical results is that the global order parameter, R can in general be optimized by aligning the vector of intrinsic frequencies ω with the dominant eigenvector v N of the Laplacian matrix L (where the elements of the Laplacian are defined as L ij = δ ij k i − A ij ).Given the enhanced synchronization found here, it is interesting to ask whether the local rewiring rule that we study generates networks with similar spectral properties.To address this question, we examine the quantity where ω || ω|| is the (normalized) vector of natural frequency offsets, and v N ||v N || is the (normalized) dominant Laplacian eigenvector.In order to understand the dependence of this frequency-eigenvector alignment on the coupling strength, at each value of α we compute | on the final, co-evolved networks G , as well as on the original ER networks G o , for comparison.Fig. 10 (a) and (b) depict the results of this analysis for the two different values of the mean degree (see Fig. 22 for the corresponding analysis with normally distributed frequencies).At low coupling, | ω || ω|| , v N ||v N || | for the rewired networks remains low, near its initial value.Note that in this regime the order parameter is low as well (Fig. 3).As α increases, however, the frequency-eigenvector overlap begins to grow, and rapidly increases until a leveling out with further increase in the coupling; the plateau occurs close to the maximal value of one for both values of the mean degree.As with the relationships between natural frequency and node degree, and between natural frequency and neighbor natural frequency considered in the previous section, the increasing projection of the intrinsic frequencies onto the dominant Laplacian eigenvector is accompanied by an increase in R. Thus, the adaptive reconfiguration can cause persistent changes to the organization of the network that are consistent with the conditions predicted by Skardal et al. [65] for optimizing synchronization, and the result becomes more promi-nent at higher coupling.The investigation for the case of normally-distributed frequencies yields similar conclusions (see Appendix B 4, Fig. 22).This analysis provides a more analytical understanding of how the adaptive process is able to enhance the synchronization in the system.

V. DISCUSSION AND CONCLUSIONS
In this study, we examined co-evolution of network topology and Kuramoto phase-oscillator dynamics as a means to evolve initially unstructured networks towards organized architectures, and to simultaneously enhance synchronization in the system.In terms of the interplay between network topology and synchronization, it has been found that the presence of specific correlations between the structural layout of the network and the oscillator frequencies (which is a property of the dynamics) can greatly augment global synchronization.But these relationships usually arise through optimization strategies that utilize global information about the network or of node states as a function of time [60,61,[65][66][67][106][107][108][109][110][111][112].On the other hand, a different set of work has shown that adaptive strategies that suppress phase differences between Kuramoto oscillators [99,100] can lead to heightened synchronization.But in the latter case, the interactions between the topology of the adaptive networks and the node frequencies has not been explored or examined in depth.An interesting line of investigation is to therefore understand whether an adaptive rule for updating the structure of the network -based on local dynamical information -can shape the topological patterns and correlations with dynamical properties whose emergence simultaneously enhances synchronization.To this end, we studied a type of disassortative mechanism in which the edge between a randomly selected node and its most instantaneously synchronized neighbor is stochastically rewired, while all other edges are maintained.Co-evolution of the dynamics and network connectivity occurs through the repetition of this feedback process, whereby an initially random network continually reconfigures in response to the states of locally connected oscillators.
Through numerical simulation, we examined the timeevolution of this process and the dependence on the global coupling.We found that for a significant coupling range, the rewiring strategy was able to bring the system to a state of heightened collective behavior, as measured by the global order parameter.It is interesting to note that this eventual state of enhanced global coherency depended on the adaptive prevention of local synchronization, suggesting a trade-off between local and global dynamics.Other work on adaptive Kuramoto networks has shown that the opposite type of rule, i.e. a competitive strategy which strengthens connections between more in-phase oscillators at the expense of weakening connections elsewhere, can lead to modular organization and hence enhanced local rather than global synchrony [84,85,101].Perhaps most importantly, the enhancement of synchronization indeed co-occurred with the emergence of correlations in the network that have been shown to arise through optimization of the global order parameter.In particular, when the oscillators exhibited more coherent dynamics, the evolved networks tended to exhibit (1 ) positive correlations between node degrees and the magnitude of oscillators' difference in natural frequency from the mean of the population (i.e.magnitude of their frequency offset), and (2 ) the preference of connections between oscillators that have natural frequency offsets of opposite sign.We found that the emergence of these relationships and how the adaptive scheme reorganized the network topology over time depended on the global coupling parameter and the intrinsic frequencies, and -to some extent -the density of the network.We note that in addition to enhancing synchronization [60,61,110,111], the purposeful placement of these types of correlations has also been associated with first order, or explosive, synchronization transitions [62][63][64].Far fewer studies, however, have considered how these structural patterns and relationships might arise in a network from local rearrangements or adaptation (though see [103] for one example).
It is important to state that the results found in this study are in line with previous work that has examined adaptive schemes in which weights grow or shrink as a function of phase differences, which also find that rules that actively suppress differences in state are able to improve synchronization.However, there are some important distinctions between the strategies studied in [99] and the one studied here.In regards to the former, the network starts as completely disconnected, with no edges between any oscillators.The edge weights between all pairs of nodes are then allowed to increase, up to some bound.In this way, the total density of the network is allowed to change, connections can theoretically occur between all pairs of oscillators at a given time, and the individual edge weights can fluctuate.Other work has considered a situation in which the topology of the network is pre-defined and constant, while the weights can change [100].On the other hand, we wished to consider a situation in which the effect of topological organization alone could be isolated from other confounding features.We thus studied a case where the network begins connected, but in a random, disorganized arrangement, and then allowed the system to self-organize under the constraint of fixed total density and also binary and undirected connectivity.Together, these conditions mean that only a fraction of the nodes can be directly connected at a given time, and the goal is to understand if simple rearrangements in those connections, based on a local rule for determining the rewiring, can enhance synchronization.Thus, the path to a more coherent state is different here than in previous studies on co-evolutionary Kuramoto systems.
Of course, there are still some methodological considerations to make note of, as well as possibilities for future work.For example, we studied an adaptive mechanism that utilized only local information of a given node, and that preserved binary connectivity and the total number of edges, so as to isolate the effects of rearrangements in network topology from other factors.However, incorporating these constraints required that there be a random component to the rewiring process, and there was not a natural or self-employed stopping condition for network reconfiguration.One could thus further explore how results are affected by the length of time the system is allowed to co-evolve, and also how this relates to changes in other parameters, such as the size of the network, the mean degree, and the spread in the intrinsic frequencies.Another parameter that may be important to examine more in depth is the time-scale of the adaptation in the network.In addition, while some stochasticity or noise is likely a realistic feature of natural systems, it would be interesting to incorporate that into a related adaptation model that allows for weighted rather than just binary connectivity between network units.It is also important to note that, while meaningful insights can be gained from the empirical and observational type of analysis carried out in this study, in forthcoming work it will be useful to try and understand the origins of various results from a more fundamental and theoretical standpoint.Finally, we point out that continued investigation into the role of network topology and adaptation in shaping dynamics and structure may lead to a better understanding of the development and function of real-world networks, such as neuronal assemblies or large-scale brain structure and activity patterns.Indeed, there are many computational models of these systems in which this can be studied [6-9, 11, 79, 98, 115-118].In addition, for neural systems in particular, it is interesting to note that while synchronization is often a desired property, hyper-synchronization can also be detrimental, as is the case with epileptic seizures [119].Therefore, further examination of the tradeoff and transition between local and global synchrony in biologically motivated models [120,121] -and understanding how this might occur adaptively over time and influence the structure and function of a network -continues to be an exciting line of study.
In conclusion, understanding the concurrent influence of network architecture on the emergence of collective dynamics [2-4, 20, 21], and in turn, the effect of a dynamical process on reshaping or inducing network structure [70][71][72], is currently an active area of research across a broad set of disciplines, including the physical, social, and biological sciences.Along these lines, we have studied a dynamical rewiring scheme for networks of Kuramoto oscillators.The adaptive rule for the network was specifically inspired by previous work on optimizing networks for synchronization of heterogeneous oscillators [60,61,[65][66][67][106][107][108][109][110][111][112], from which we wished to uncover how similar dynamics and organization could occur through a co-evolution process.We found that a restructuring of the network, in which the effect of incoherence is suppressed by maintaining edges between disparate oscillators, while randomly rewiring edges between the most locally and instantaneously in-phase oscillators, led to the emergence of distinct topological patterns and correlations that concurrently enhance the system's ability to synchronize.This study thus sheds light on a mechanism for how enhanced synchronization and network structure might arise in a system that evolves and reconfigures according to local information alone, without knowledge of global connectivity or node states.

VI. SUPPLEMENTARY MATERIAL
Please see the Supplementary Material (S.M.) for an analysis of the time-evolution of locally synchronized clusters, and also for an examination of the robustness of some of the main results to simple variations in the network size, the initial network topology, and the frequency distribution.We also thank two anonymous reviewers whose comments greatly improved the quality of this work.The content is solely the responsibility of the authors and does not necessarily represent the official views of any of the funding agencies.
Appendix A: An additional measure to quantify relationships between network topology and oscillator frequencies In Sec.IV B, we used two measures to quantify emerging relationships between the intrinsic frequencies of the oscillators and the network topology: the correlation C |ω|,k between the magnitude of the oscillator frequency offset |ω i | and the node degree k i , and the mean fraction f of an oscillator's neighbors that have natural frequency offsets of opposite sign compared to the frequency offset of the central oscillator (averaged across all nodes in the system) [60,61].Here we define a related measure to further understand the organization that arises in the adaptively rewired networks.In particular, for each oscillator i, we assess the association between the frequency offset ωi and the sum of oscillator i s neighbors' frequency offsets j∈Ni ωj .We quantify this relationship in the system by considering the correlation C ω, ω between ωi and j∈Ni ωj across all nodes in the network.FIG.11.The relationship between the oscillator frequency offset ωi and the sum of neighbor frequency offsets j∈N i ωj can be used to further quantify the interplay between the network structure and the intrinsic frequencies of the oscillators that arises due to the co-evolutionary process.The top row shows an example of this relationship for a network with k = 12.5, and the bottom row shows examples for a network with k = 25; in both cases, the frequencies were drawn from the uniform distribution {ωU }.For each network density, the first column corresponds to an ER random graph that exhibits only intermediate levels of synchrony at the displayed coupling α (as measured by R ), and the second column corresponds to the adapted network, which exhibits a higher level of synchrony.(a,b); (d,e) The total neighbor frequency offset j∈N i ωj vs. frequency offset ωi.(c); (f ) The correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets j∈N i ωj vs. the number of rewiring steps m.Fig. 11 gives examples of this relationship for the same networks as those in Fig. 4 in the main text.The top and bottom rows correspond to networks with k = 12.5 and k = 25, respectively.Panels (a) and (d) show that in the ER networks, there is no clear relationship between frequency offset and the sum of neighbor frequency offsets, as expected.In the adaptively rewired networks -panels (b) and (e) -the order parameter has increased from its original value and there is a strong negative correlation between ωi and j∈Ni ωj .Not only do connections tend to form between oscillators that have frequency offsets of opposite sign (as measured by f ), but in addition, edges become distributed in the network such that the sum of each oscillator's neighbor frequency offsets proportionately cancels out each oscillator's own difference from the mean frequency of the population.Fig. 11c,f shows how the strength of this relationship -quantified by the correlation C ω, ω -evolves as the network is rewired.12.These plots depict the correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets j∈N i ωj as a function of the coupling.In each panel, the gray data points correspond to the initial, uncorrelated ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the uniform distribution {ωU }, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.
Fig. 12 shows C ω, ω as a function of the coupling α.As α increases, C ω, ω decreases to a strong negative value, and then remains approximately constant for larger values of the coupling.(Compare to Fig. 5 and Fig. 6, which show C |ω|,k vs. α and f vs. α, respectively).Finally, note that the strong decreases in C ω, ω occur in conjunction with increases in the order parameter (Fig. 3).

Appendix B: Analysis with normally distributed frequencies
The analysis in the main text was carried out using natural frequencies {ω U } drawn at random from the uniform distribution in the range [−2, 2].To demonstrate that the main results are not specific to a single choice of the frequency distribution, in this Appendix, we also consider the case of frequencies {ω G } drawn from a normal distribution with zero mean and unit standard deviation.The findings shown here are largely consistent with those reported in the main text.

Dependence of the order parameter on time and global coupling
Fig. 13 shows examples of the order parameter R(t) vs. time t at different values of the coupling and for the two different mean degrees k = 12.5 and k = 25.For these trials and values of the coupling, we observe that the order parameter increases due to the rewiring of the network, which takes place between the two red lines.(Fig. 2 from the main text shows examples for the case of uniformly distributed natural frequencies).
Fig. 14 shows examples of the time-averaged order parameter R vs. the coupling α.As with the case of uniformly distributed frequencies (Sec.IV A, Fig. 3), we find that the adapted networks exhibit heightened synchronization over a large coupling range.

Correlations between network topology and the intrinsic frequencies
Fig. 15 shows examples of node degree k i vs. intrinsic frequency offset ωi , average neighbor frequency offset ω Ni vs. node frequency offset ωi , and total neighbor frequency offset j∈Ni ωj vs. frequency offset ωi for ER networks and the corresponding co-evolved networks (see In each panel, the gray data points correspond to static ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the normal distribution {ωG}, and the mean degree of the networks were (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.
Sec. IV B and Appendix A for the definitions of these quantities).The top three and bottom three panels correspond to networks with k = 12.5 and k = 25, respectively, and the frequencies {ω G } were normally distributed.At the values of coupling used for these examples, we see that as the system reconfigures, the network begins to exhibit distinct patterns in terms of the organization of oscillators with different intrinsic frequencies.Each of the three metrics considered -C |ω|,k , f , and C ω, ω -exhibit a progression that allows for the eventual heightened degree of synchrony in the rewired network.These findings are similar to those discussed for the uniform frequency distribution (see Fig. 4 and Fig. 11 and the corresponding text in Sec.IV B).
We next show the evolution of the relationships between network topology and the intrinsic frequencies as a function of the coupling α.Figs. 16, 17, and 18 show C |ω|,k vs. α, f vs. α, and C ω, ω vs. α, respectively.The conclusions drawn for the case of normally distributed frequencies shown here are the same as those for the uniformly distributed frequencies examined in the main text.(See Figs. 5, 6, and 12 and the corresponding discussions in Sec.IV B and Appendix A).Briefly, each of the three measures exhibit transitions at similar values of the coupling, and the emergence of strong relationships between the network structure and oscillator frequencies arise near the coupling when the order parameter transitions from low to higher values.FIG. 17.These plots depict the mean fraction f (i.e.averaged over all nodes in the network) of an oscillator's neighbors that have natural frequency offsets of opposite sign compared to that of the central oscillator, as a function of the coupling.In each panel, the gray data points correspond to the initial, uncorrelated ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the normal distribution {ωG}, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.mally distributed and the same for both cases.Each row corresponds to one oscillator, and the rows are ordered by the quantity ωi = ω i − ω (i.e. the offset from the mean intrinsic frequency of the population).Adaptation of the network takes place between the two black lines.The results are qualitatively consistent with those described in the main text for the uniformly distributed frequencies (see the discussion in Sec.IV C and Fig. 7 for comparison).FIG.18.These plots depict the correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets j∈ mathcalN (i) ωj, as a function of the coupling.In each panel, the gray data points correspond to the initial, uncorrelated ER random graphs Go, and the orange and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the normal distribution {ωG}, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.
degree k i vs. the rewiring step m for the same set of networks, natural frequencies, and coupling values as in Fig. 20.In both figures, panel (a) corresponds to a network with k = 12.5, and panel (b) corresponds to a network with k = 25; the natural frequencies {ω G } are normally distributed and are the same in both cases.We observe similar types of behavior and regimes in terms of the evolution of these quantities as we did for the situation of uniformly distributed frequencies in the main text (see the discussion in Sec.IV C and Figs. 8 and 9 for comparison).One interesting result to point out for the examples using normally distributed frequencies shown here, is that the correlation C |ω|,k plateaus at a value near 0.5 (for k = 12.5) for several of the intermediate coupling values, before increasing further at higher coupling.As seen in Fig. 20, this is due to the oscillators with the most outlying natural frequencies (which can sometimes be more extreme for the normal distribution than for the uniform distribution) remaining with low degree.Though the results are qualitatively similar between the two frequency distributions, further work is needed to quantify and understand what may be subtle and intricate differences.

Spectral Analysis
In Sec.IV D we carried out a spectral-based analysis of the co-evolved networks inspired by [65].Fig. 22 shows the results of this analysis for the case of normally distributed frequencies {ω G }, which are consistent with those previously discussed in Sec.IV D using uniformly distributed frequencies (compare to Fig. 10).In short, for both values of the mean degree ( k = 12.ral frequency offsets and the normalized dominant Laplacian eigenvector of the rewired networks increases and then plateaus at a high value as the coupling increases (Fig. 22a,b).We refer the reader to Sec.IV D for a more detailed discussion of these findings.
FIG. 2. Ensemble averages of the number of synchronized clusters nsc(t) vs. time t, for various representative couplings α.In each case, the dynamics were first run atop an ER random graph, after which co-evolution of the network and dynamics took place between the two red lines.For these plots, the natural frequencies were drawn from the uniform distribution {ωU }, the mean degree of the network was k = 25, and each curve corresponds to an average over 10 different simulations.During the adaptation period, the network was continually rewired once every T = 0.2 time units.For high enough coupling, co-evolution results in a decrease in the number of synchronized clusters over time.
ulation of the dynamics, and for various values of the coupling α (the mean degree is k = 25 and the natural frequencies {ω U } are uniformly distributed).We see that the number of synchronized clusters is initially large, since the oscillators' phases are distributed at random at t = 0. Before the red line, the system evolves on a static ER network, and depending on the coupling strength, the number of components either fluctuates around a high value or -after a short initial transient -rapidly decreases to a lower value, denoting partial synchronization.Network adaptation begins at the time denoted by the red line, after which the behavior of n sc (t) is again determined by the global coupling.For very low α, the co-evolution process does not significantly affect the number of synchronized clusters.However, as the coupling increases, the network rewiring begins to have an impact on the dynamical clustering.At intermediate values of α, the number of synchronized components slowly decreases as a function of time.For higher coupling values, this transition takes place more rapidly.However, it appears that no matter the coupling regime, the transition from high to lower n sc does not occur abruptly at a single instance of time, but instead occurs over several time steps through the merging of clusters.In addition to the curves depicting individual examples of the number of  In each case, the dynamics were first run atop an ER random graph, after which co-evolution of the network and dynamics took place between the two red lines.For these plots, the natural frequencies were drawn from the normal distribution {ωG} and the mean degree of the network was k = 12.5.
During the adaptation period, the network was continually rewired once every T = 0.2 time units.For high enough coupling, co-evolution results in a decrease in the number of synchronized clusters over time.
synchronized clusters vs. time, Fig. S2 also shows averages of n sc (t) vs. t over 10 simulations using different instantiations of the initial conditions.Fig. S3 and Fig. S4 show n sc (t) vs. time t and the averages n sc (t) vs. t, respectively, for the case of mean degree k = 12.5 and normally distributed natural frequencies {ω G }.We find similar results to those for the denser networks.At low values of the coupling, network adaptation does not cause a consistent change in n sc .As the coupling increases, however, the number of components begins to decrease from its initial value, signifying the merging or growth of synchronized clusters due to rearrangements in the network.This decrease becomes more rapid at larger couplings, though does not occur discontinuously as a function of time.

II. VARIATION OF NETWORK SIZE, INITIAL NETWORK TOPOLOGY, AND FREQUENCY DISTRIBUTION
In the main text (and Appendix), we fixed the system size at N = 100 and initialized the coupling between oscillators to be ER random graphs with mean degree k = 25 or k = 12.5.We then examined two symmetric distributions (uniform {ω U } and normal {ω G }) for the FIG. 4. Ensemble averages of the number of synchronized clusters nsc(t) vs. time t, for various representative couplings α.In each case, the dynamics were first run atop an ER random graph, after which co-evolution of the network and dynamics took place between the two red lines.For these plots, the natural frequencies were drawn from the normal distribution {ωG}, the mean degree of the network was k = 12.5, and each curve corresponds to an average over 10 different simulations.During the adaptation period, the network was continually rewired once every T = 0.2 time units.For high enough coupling, co-evolution results in a decrease in the number of synchronized clusters over time.
intrinsic frequencies of the oscillators.In this section, we explore the robustness of some of the main results to simple variations in the network size, the initial network topology, and the frequency distribution.

A. Larger networks
We first investigate the effect of the co-evolutionary mechanism on the degree of synchronization in the system when the network is doubled in size to N = 200 (while also proportionately doubling the number of rewirings allowed for the adaptation process).Fig. S5a shows the order parameter R vs. coupling α for the case of uniformly distributed frequencies {ω U }, mean degree k = 25, and a rewiring time scale of T = 0.2 time units.We find qualitatively similar results to that of N = 100 (Fig. 3b of the main text), in that the curve corresponding to the adapted networks (blue) is shifted to the left compared to the curve corresponding to a set of static, random ER networks (gray).Thus, over a large range of global coupling, the co-evolved networks exhibit a heightened degree of collective dynamics.
In addition to the order parameter, we also assessed the co-evolved networks G for the emergence of correlations between their topology and the intrinsic frequencies of the oscillators (as done in Sec.IV B and Appendix A of the main text).The results of this analysis are shown in Fig. S5(b -d).As found previously for the case of N = 100 (Fig. 5b, 6b, and 12b), we find that the same set of relationships appear when the system size is doubled.Following the analysis described in the main text, for each node i, we computed its offset from the mean frequency ωi = ω i − ω , where ω is the mean natural frequency of the population.We then considered the correlations between |ω i | and node degree k i as a function of the coupling (Fig. S5b), and between ωi and the total neighbor frequency offset ωj (∀j ∈ N (i) as a function of the coupling (Fig. S5d ).In addition, for each oscillator i, we computed the fraction f i of i s neighbors that have a value of ω of opposite sign to that of ωi , and then computed the average f over all nodes in the network (Fig. S5c).As the coupling increases, there is an increasingly positive correlation C |ω|,k between node de-gree and the magnitude of the intrinsic frequency offset, and an increasingly negative correlation C ω, ω between the intrinsic frequency offset of each oscillator and the sum of its neighbors' frequency offsets.Furthermore, oscillators tend to become connected to other oscillators that have intrinsic frequency offsets on the opposite side of the mean.Although an in-depth examination into the effect of system size is beyond the scope of this study, it is important that the results are consistent for at least somewhat larger networks.We refer the reader to [2] for a detailed finite-size study and analysis of the critical coupling and critical exponents in oscillator populations with purposefully constructed correlations.

B. Smaller mean degree
This study has focused on networks with relatively large mean degree (yielding density values similar to that of large-scale anatomical brain networks, for example [3], but many other real-world systems -such as social networks [4] -are more sparsely connected.Here, we thus briefly examine whether various results hold for networks with a lower mean degree of k = 6.Fig. S6a shows the time-averaged order parameter R vs. coupling α for the static ER networks (gray curve) and the adapted networks (blue curve).The frequencies {ω U } were uniformly distributed the rewiring time scale was T = 0.2 time units.As for the denser networks considered in the main text (Fig. 3), here we also observe a significant enhancement of the order parameter in the co-evolved systems.
We also examined the locally adapted systems for the previously described relationships between the natural frequencies and the topological organization of the network.Fig. S6 panels (b-d) show C |ω|,k , f , and C ω, ω , respectively, as a function of the coupling.Again we find that as the coupling increases, there is an increasingly positive correlation between degree and the magnitude of intrinsic frequency offset, and an increasingly negative correlation between the intrinsic frequency offset of an oscillator and the sum of its neighbors' frequency offsets.Furthermore, oscillators tend to become connected to other oscillators with frequency offsets on the opposite side of the mean.Although the curves shown here for the case of k = 6 exhibit somewhat more fluctuations compared to the case of k = 25, the findings are in general robust and consistent with those found for networks with higher mean degree (compare to Figs. 5, 6, and 12 of the main text).

C. Initial network topology
The main interest and goal of this paper has been to define and explore a simple, local mechanism that evolves some initial network of coupled oscillators to a structured configuration that in turn supports enhanced collective dynamics.Given this objective, ER random graphs -which are considered to be relatively "disordered" -are a natural family of networks on which to begin the co-evolution and also to compare the adapted systems against.However, it may still be useful to examine if the results hold for a different initial network topology.To test this, we constructed a set of modular -small world (MSW) networks, which consist of a specified number of fully connected communities that are then linked together by placing edges at random between the modules.In contrast to the ER random graph model, these networks contain a significant amount of planted structure.For this analysis, we used networks of size N = 128, with N c = 8 modules, and a total number of edges such that the average degree was k = 20.Fig. S7 shows an example connection pattern.
Due to the nature of the adaptive mechanism, we expect that it will undo the seeded community structure and evolve this class of networks towards topologies with the same set of features we have discussed previously.To modular -small world network, test this expectation, we first simulated the Kuramoto dynamics on non-adaptive MSW networks, and then used the same set of MSW networks as the initial network configurations for the adaptive mechanism.Fig. S8a shows the time-averaged order parameter R vs. coupling α for the case of uniformly distributed frequencies {ω U }, and a rewiring time scale of T = 0.2 time units.The adaptively rewired networks (blue curve) show heightened global synchronization across a large range of coupling values compared to the static MSW networks (gray curve).Turning next to the correlations between the natural frequencies and the network structure, we find results consistent with those described in the main text.In particular, strong degree -frequency correlations and frequency -neighbor frequency mixing emerge in the selforganized systems (see Fig. S8b -d for plots of C |ω|,k , f , and C ω, ω vs. the coupling α).This investigation provides some evidence that the adaptive mechanism is agnostic to the initial network configuration from which it begins.

D. Natural frequency distribution
One can also explore in more depth how the distribution of natural frequencies affects certain results.In this study, we considered two standard distributions -uniform and normal -that are often used in studies on the Kuramoto model.Here, we also show some brief results for the case of natural frequencies {ω P } drawn from a power-law distribution in which P (ω) ∝ ω −γ , with γ = 3.The dynamics and network co-evolution start from ER random graphs with N = 200 nodes and mean degree k = 25.We find qualitatively similar results to those in the main text when using this skewed distribution for the intrinsic frequencies.Fig. S9 shows the time-averaged order parameter R vs. coupling α for the static ER networks (gray curve) and adaptively rewired networks (blue curve).The rewiring time scale was T = 0.2 time units.As with the two symmetric frequency distributions, we observe an increase in the order parameter for the networks evolved according to the local mechanism.We also examined the co-evolved networks G for the correlations between the natural frequencies of the oscillators and the topological organization of the network.Fig. S9b -d show plots of C |ω|,k , f , and C ω, ω vs. the coupling α.Consistent with the findings in the main text, we observe the emergence of strong degree-frequency correlations and frequencyneighbor frequency mixing patterns.
FIG.2.Examples of the global order parameter R(t) vs. time t for various representative couplings α.In each case, the dynamics were first run atop an initially ER random graph with average degree k , after which co-evolution of the network and dynamics took place between the two red lines.The natural frequencies were drawn from the uniform distribution {ωU }, and the mean degree k and coupling α used for each panel were (a) k = 25, α = 0.065, (b) k = 25, α = 0.085, (c) k = 12.5, α = 0.135, (d) k = 12.5, α = 0.22.During the adaptation period, the network was continually rewired once every T = 0.2 time units.The co-evolving networks exhibit enhanced collective dynamics, as observed by increases in the global order parameter.

2 FIG. 3 .
FIG.3.The time-averaged order parameter R vs. coupling α.In each panel, the gray data points correspond to static ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the uniform distribution {ωU }, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.

FIG. 4 .
FIG. 4. Relationships between the network structure and the intrinsic frequencies of the oscillators.The top two rows show examples for a network with k = 12.5, and the bottom two rows show examples for a network with k = 25; in both cases, the frequencies were drawn from the uniform distribution {ωU }.For each network density, the first column corresponds to an ER random graph that exhibits only intermediate levels of synchrony at the displayed coupling α (as measured by R ), and the second column corresponds to the adapted network, which exhibits a higher level of synchrony.These plots highlight key relationships that emerge from the co-evolutionary network update rule.(a,b); (g,h) Node degree ki vs. frequency offset ωi.(d,e); (j,k) Average neighbor frequency offset ω N i vs. frequency offset ωi.(c); (i) The correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi| vs. the number of rewiring steps m, and (f ); (l) the mean fraction f (i.e.averaged over all nodes in the network) of an oscillator's neighbors that have frequency offsets of opposite sign compared to that of the central oscillator vs. the number of rewiring steps m.

2 FIG. 5 .
FIG.5.These plots depict the correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi|, as a function of the coupling.In each panel, the gray data points correspond to the initial, uncorrelated ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the uniform distribution {ωU }, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.The dip observed in (b) at α ≈ 0.5 is due to the localization of edges on a cluster of oscillators with natural frequencies near the mean; this is examined further in Sec.IV C. All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.

FIG. 6 .
FIG.6.These plots depict the mean fraction f (averaged over all nodes in the network) of an oscillator's neighbors that have frequency offsets of opposite sign compared to that of the central oscillator, as a function of the coupling.In each panel, the gray data points correspond to the initial, uncorrelated ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the uniform distribution {ωU }, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.

FIG. 7 .
FIG. 7. Examples of the instantaneous frequencies θi(t) vs. time t, for various representative couplings α.The mean degree of the networks are (a) k = 12.5 and (b) k = 25, and the natural frequencies {ωU } were drawn from the uniform distribution.In all panels, each row corresponds to one oscillator, and the rows are ordered by the quantity ωi = ωi − ω (i.e. the offset from the mean intrinsic frequency of the population).For each coupling, the dynamics were first run atop an initially ER random graph, after which co-evolution of the network and dynamics took place between the two black lines.During the adaptation period, the network was continually rewired once every T = 0.2 time units.

FIG. 8 .
FIG. 8. Examples of the evolution of the node degree ki vs. the rewiring step m, for various representative couplings α.The mean degree of the networks are (a) k = 12.5 and (b) k = 25, and the natural frequencies {ωU } were drawn from the uniform distribution.In all panels, each row corresponds to one oscillator, and the rows are ordered by the quantity ωi = ωi − ω (i.e. the offset from the mean intrinsic frequency of the population).The network was continually rewired once every T = 0.2 time units.

14 FIG. 9 .
FIG. 9. Examples of the evolution of the correlation C | ω|,k between oscillator frequency offset ωi and the node degree ki,vs.the rewiring step m, for various representative couplings α.The mean degree of the networks are (a) k = 12.5 and (b) k = 25, and the natural frequencies {ωU } were drawn from the uniform distribution.In all panels, each row corresponds to one oscillator, and the rows are ordered by the quantity ωi = ωi − ω (i.e. the offset from the mean intrinsic frequency of the population).The network was continually rewired once every T = 0.2 time units.

2 FIG. 10 .
FIG. 10.Evolution of the overlap | ω|| ω|| , v N ||v N || | between the intrinsic frequencies ω and dominant Laplacian eigenvector vN as a function of the coupling α.In each panel, the gray data points correspond to the original, uncorrelated ER random graphs Go, and the blue and yellow curves correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The natural frequencies were drawn from the uniform distribution {ωU }, and the mean degree k of the networks are (a) k = 12.5 and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.
VII. ACKNOWLEDGMENTS L.P. acknowledges support from the National Science Foundation Graduate Research Fellowship Program.J.K. acknowledges support from the National Science Foundation Graduate Research Fellowship Program and NIH T32-EB020087, PD: Felix W. Wehrli.D.S.B. also acknowledges support from the John D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, and the National Science Foundation (BCS-1441502, CA-REER PHY-1554488, BCS-1631550, and CNS-1626008).
FIG.12.These plots depict the correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets j∈N i ωj as a function of the coupling.In each panel, the gray data points correspond to the initial, uncorrelated ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the uniform distribution {ωU }, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.

)
FIG.13.Examples of the global order parameter R(t) vs. time t, for various representative couplings α.In each case, the dynamics were first run atop an initially ER random graph with average degree k , after which co-evolution of the network and dynamics took place between the two red lines.The natural frequencies were drawn from a normal distribution {ωG}, and the mean degree k and coupling α used for each panel were (a) k = 25, α = 0.05, (b) k = 25, α = 0.075, (c) k = 12.5, α = 0.1, (d) k = 12.5, α = 0.2.During the adaptation period, the network was continually rewired once every T = 0.2 time units.The co-evolving networks exhibit enhanced collective dynamics, as observed by increases in the global order parameter.

2 FIG. 14 .
FIG.14.The time-averaged order parameter R vs. coupling α.In each panel, the gray data points correspond to static ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the normal distribution {ωG}, and the mean degree of the networks were (a) k = 12.5, and (b) k = 25.All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.

Fig. 19
Fig. 19 shows examples of θi (t) vs. t for several values of the coupling α around the point in which the dynamics transition from an incoherent state to a synchronized state.The top set of panels (a) are for a network with k = 12.5, and the bottom set of panels (b) are for a network with k = 25; the frequencies {ω G } were nor-

FIG. 15 .
FIG. 15.Relationships between the network structure and the intrinsic frequencies of the oscillators.The top three rows show examples for a network with k = 12.5, and the bottom three rows show examples for a network with k = 25; in both cases, the frequencies were drawn from the normal distribution {ωG}.For each network density, the first column corresponds to an ER random graph that exhibits only intermediate levels of synchrony at the displayed coupling α (as measured by R ), and the second column corresponds to the adapted network, which exhibits a higher level of synchrony.These plots highlight key relationships that emerge from the co-evolutionary network update rule.(a,b); (j,k) Node degree ki vs. frequency offset ωi.(d,e); (m,n) Average neighbor frequency offset ω N i vs. frequency offset ωi.(g,h); (p,q) Total neighbor frequency offset j∈N i ωj vs. frequency offset ωi.(c); (l) The correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi| vs. the number of rewiring steps m, and (f ); (o) The mean fraction f (i.e.averaged over all nodes in the network) of an oscillator's neighbors that have frequency offsets of opposite sign compared to that of the central oscillator vs. the number of rewiring steps m. (i); (r) The correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets j∈N i ωj vs. the number of rewiring steps m.

2 FIG. 16 .
FIG.16.These plots depict the correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi|, as a function of the coupling.In each panel, the gray data points correspond to the initial, uncorrelated ER random graphs Go, and the blue and yellow points correspond to the adapted networks G evolved under rewiring time scales of T = 0.2 and T = 2, respectively.The frequencies were drawn from the normal distribution {ωG}, and the mean degree of the networks are (a) k = 12.5, and (b) k = 25.The dip observed in (b) at α ≈ 0.4 is due to the localization of edges on a cluster of oscillators with natural frequencies near the mean of the distribution; this is examined further in Appendix B 3. All curves depict averages over 25 instantiations, and the lines between data points serve as guides for the eye.

Fig. 20
shows the node degree k i vs. the rewiring step m and Fig. 21 shows the correlation C |ω|,k between the absolute value of the frequency offset |ω i | and the node coupling, α coupling, α

FIG. 19 .
FIG.19.Examples of the instantaneous frequencies θi(t) vs. time t, for various representative couplings α.The mean degree of the networks are (a) k = 12.5 and (b) k = 25, and the natural frequencies {ωG} were drawn from the normal distribution.In all panels, each row corresponds to one oscillator, and the rows are ordered by the quantity ωi = ωi − ω (i.e. the offset from the mean intrinsic frequency of the population).For each coupling, the dynamics were first run atop an initially ER random graph, after which co-evolution of the network and dynamics took place between the two black lines.During the adaptation period, the network was continually rewired once every T = 0.2 time units.

FIG. 20 .FIG. 21 . 2 FIG. 22 .
FIG. 20.Examples of the evolution of the node degree ki vs. the rewiring step m, for various representative couplings α.The mean degree of the networks are (a) k = 12.5 and (b) k = 25, and the natural frequencies {ωG} were drawn from the normal distribution.In all panels, each row corresponds to one oscillator, and the rows are ordered by the quantity ωi = ωi − ω (i.e. the offset from the mean intrinsic frequency of the population).The network was continually rewired once every T = 0.2 time units.

FIG. 5 .
FIG. 5. Analysis of the synchronization and various correlations between the network topology and the intrinsic frequencies for networks of size N = 200.In each panel, the gray data points correspond to the initial ER networks Go with mean degree k = 25, and the blue points correspond to the adapted networks G evolved under a rewiring time scale of T = 0.2 time units.The natural frequencies were drawn from the uniform distribution {ωU }.(a) The time-averaged order parameter vs. coupling.(b) The correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi|, as a function of the coupling.(c) The mean fraction f of an oscillator's neighbors that have frequency offsets of opposite sign compared to the central oscillator, as a function of the coupling.(d) The correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets ωj, as a function of the coupling.All curves depict averages over 10 instantiations, and the lines between data points serve as guides for the eye.

FIG. 6 .
FIG.6.Analysis of the synchronization and various correlations between the network topology and the intrinsic frequencies for networks with sparser connectivity.In each panel, the gray data points correspond to the initial ER networks Go with mean degree k = 6, and the blue points correspond to the adapted networks G evolved under a rewiring time scale of T = 0.2 time units.The natural frequencies were drawn from the uniform distribution {ωU }.(a) The time-averaged order parameter vs. coupling.(b) The correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi|, as a function of the coupling.(c) The mean fraction f of an oscillator's neighbors that have frequency offsets of opposite sign compared to the central oscillator, as a function of the coupling.(d) The correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets ωj, as a function of the coupling.All curves depict averages over 10 instantiations, and the lines between data points serve as guides for the eye.

FIG. 7 .
FIG. 7.An example connection pattern of a modular -small world (MSW) network of size N = 128, with 8 communities and a mean degree k = 20.

FIG. 8 .
FIG.8.Analysis of the synchronization and various correlations between the network topology and the intrinsic frequencies when the initial network is modular-small world, rather than a random graph.In each panel, the gray data points correspond to the initial MSW networks Go, and the blue points correspond to the adapted networks G evolved under a rewiring time scale of T = 0.2 time units.The natural frequencies were drawn from the uniform distribution {ωU }.(a) The time-averaged order parameter vs. coupling.(b) The correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi|, as a function of the coupling.(c) The mean fraction f of an oscillator's neighbors that have frequency offsets of opposite sign compared to the central oscillator, as a function of the coupling.(d) The correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets ωj, as a function of the coupling.All curves depict averages over 10 instantiations, and the lines between data points serve as guides for the eye.

2 fFIG. 9 .
FIG.9.Analysis of the synchronization and various correlations between the network topology and the intrinsic frequencies when the frequencies are drawn from a power law distribution with exponent γ = 3.In each panel, the gray data points correspond to the initial ER networks Go, and the blue points correspond to the adapted networks G evolved under a rewiring time scale of T = 0.2 time units.The network size is N = 200 and their mean degree is k = 25.(a) The time-averaged order parameter vs. coupling.(b) The correlation C | ω|,k between node degree ki and the magnitude of the frequency offset |ωi|.(c) The mean fraction f of an oscillator's neighbors that have frequency offsets of opposite sign compared to the central oscillator, as a function of the coupling.(d) The correlation C ω, ω between oscillator frequency offset ωi and the sum of neighbor frequency offsets ωj, as a function of the coupling.All curves depict averages over 10 instantiations, and the lines between data points serve as guides for the eye.
13G.13.Examples of the global order parameter R(t) vs. time t, for various representative couplings α.In each case, the dynamics were first run atop an initially ER random graph with average degree k , after which co-evolution of the network and dynamics took place between the two red lines.The natural frequencies were drawn from a normal distribution {ωG}, and the mean degree k and coupling α used for each panel were