# Heterogeneity induces emergent functional networks for synchronization

###### Abstract

We study the evolution of heterogeneous networks of oscillators subject to a state-dependent interconnection rule. We find that heterogeneity in the node dynamics is key in organizing the architecture of the functional emerging networks. We demonstrate that increasing heterogeneity among the nodes in state-dependent networks of phase oscillators causes a differentiation in the activation probabilities of the links. This, in turn, yields the formation of hubs associated to nodes with larger distances from the average frequency of the ensemble. Our generic local evolutionary strategy can be used to solve a wide range of synchronization and control problems.

###### pacs:

05.45.Xt,05.65.+b,05.10.-a,87.23.Kg^{†}

^{†}preprint: APS/123-QED

## I Introduction

Evolution is a fundamental force driving the organization and structure of natural systems. It is based on two key ingredients: variation and natural selection Darwin (1951). The first ensures the necessary mutation and recombination generating new species while the second determines the survival of the fittest to perform a certain function. Networks in Nature have been subject to the same powerful mechanisms that ultimately determined their structure, properties and functionality. The resulting networks have heterogenous topological structures, which researchers have been interested in together with their effects on dynamical processes Tan et al. (2014). Examples include epidemic spreading, opinion formation, and synchronization Boccaletti et al. (2006); Barrat et al. (2008). Often there is also heterogeneity in the nodes of a network. For example, in social networks, individuals have different personalities, which will have great impacts on their social relationships; or, in manufacturing, industrial products are slightly different from one other, affecting their impact and market shares. The relationship between the heterogeneity of the nodes and the structural properties of a network is little understood, particularly when the network evolution is state-dependent.

Here we suggest that heterogeneity in the nodes is a driving force behind the evolution of the network structure that determines its properties and function. To test this ansatz we take as a representative example the problem of evolving the network structure to achieve synchronization of coupled oscillators. This is one of the best understood and most widely studied type of collective behavior on networks Kuramoto (1984); Strogatz (2000); Acebrón et al. (2005); Arenas et al. (2008); Barrat et al. (2008).

So far, optimal network structures for synchronization have been studied mainly by using Monte Carlo methods Donetti et al. (2005); Rad et al. (2008); Yanagita and Mikhailov (2010); Gorochowski et al. (2010) or gradient-based learning strategies Tanaka and Aoyagi (2008); Skardal et al. (2014). These are based on the use of some objective function for synchronization (as for example the order parameter) which is used to find the optimal network whose structural properties are then surveyed. The Monte Carlo approach is a generic and powerful strategy but it is typically time-consuming, and increasingly cumbersome to apply to large-scale networks. Gradient-based methods assume some constraints to derive the evolution rule of the coupling strengths and the rules are often not local, in the sense that some global information on the entire network is used. Also, it has been shown that adaptive networks can yield the emergence of modular and scale-free structures, while enhancing synchronization Gutiérrez et al. (2011).

In this paper, we propose the use of an evolutionary strategy to find a functional structure for synchronization in a network of heterogeneous oscillators. In so doing we will show that heterogeneity in the nodes is instrumental in determining the properties of the resulting network. The goal of the strategy is to identify, over all possible unweighted network configurations, the structure with a minimal number of links, which guarantees frequency synchronization of its nodes. While the fundamental aim of our study is similar to that of the literature Donetti et al. (2005); Rad et al. (2008); Yanagita and Mikhailov (2010); Gorochowski et al. (2010); Tanaka and Aoyagi (2008); Skardal et al. (2014); Gutiérrez et al. (2011), the approach we propose is completely different. Indeed, our strategy uses adaptive schemes which are completely local and do not rely on any global synchronization measure. Moreover such schemes are deployed in a novel evolutionary manner.

## Ii Problem statement

We start by considering a network of general nonlinear coupled oscillators

(1) |

where is the -dimensional state of the -th oscillator, denotes its dynamics (note that oscillators can be slightly different from each other due to both parameters and model mismatches), is a generic coupling function and are time-varying coupling gains determining the strength of the coupling between neighboring oscillators.

We model the evolutionary pressures to reach synchronization by considering state-dependent second-order nonlinear dynamics for the gains dependent upon a double well potential . The gain dynamics are given by

(2) |

in which is a generic increasing function such that Note that this is a very general adaptive network equation relying on a decentralized, local, state-dependent interconnection rule. This system can be systematically reduced, under a standard technique Kuramoto (1984), to the network of adaptively coupled phase oscillators:

(3) | |||

(4) |

in which is the phase of the -th generic oscillator, is a generic -periodic function. We set the overall coupling strength to a unitary value, since it can be absorbed into a parameter defining the heterogeneity of the natural frequencies by rescaling time, i.e. by setting . In this paper we analyze, for the sake of clarity, the simplest case

(5) | ||||

(6) |

The effectiveness of edge snapping strategies to achieve synchronization has been discussed in DeLellis et al. (2010) and further details are given in Appendix A. Under such a forcing the dynamics of (starting from zero initial conditions and ), will either converge towards 0 (link is not present) or towards 1 (link is activated)

The differences in the natural frequencies of the oscillators originate from the heterogeneity of the node dynamics in weakly coupled nonlinear oscillators Kuramoto (1984). In what follows, these natural frequencies are selected deterministically from a Gaussian distribution with zero mean and standard deviation equal to . Therefore, the parameter can be used to “tune” the level of heterogeneity among nodes.

We note here that when the number of nodes is not so large, such as or , the natural frequencies sampled from a distribution can be biased. To avoid the effect of the biased sampling, we deterministically select the natural frequencies of the oscillators, similarly to Yanagita and Mikhailov (2010), as the -tuple satisfying the constraints:

where is the probability density function of a given distribution. It should be noted that for a large network, we performed our simulation taking the natural frequencies randomly from a distribution and the obtained results are qualitatively the same.

Next, we investigate how the evolution of the network is affected by tuning the heterogeneity in the nodes. To this aim we use the edge snapping strategy described above in a novel evolutionary manner (see Fig. 1) as explained in the next section.

## Iii Evolutionary Edge-Snapping

The evolutionary Edge-Snapping technique is based on two fundamental steps: one implementing the variation ingredient of evolution, the other its selection mechanism.

To implement the variation ingredient of evolution, a set of unweighted networks is generated using equations (3) and (4) starting the process from different sets of initial conditions. We consider a set of initial conditions randomly selected using a Latin Hypercube strategy McKay et al. (1979) in the range , . To obtain the “fitness” of each link, we next compute the probability of each link being activated as the fraction between the number of generated networks where that link is present, say , and the total number of trials, e.g. . This yields a stochastic matrix whose elements are the probabilities of activation of every possible link among nodes.

The selection rule is obtained by selecting only those links whose activation probability is above a certain critical threshold value , i.e. such that . We choose so as to guarantee that the resulting network is connected and has the smallest number of links. We shall term such a network as the minimal edge-snapping (ES) network.

The variation step of our evolutionary strategy relies on the generation of a set of unweighted network using equations (3) and (4) and starting the process from a different set of initial condition. With the aim of choosing a reasonable value for the number of trials , we plot in Fig. 2(a) the standard deviation of the link activation probability as a function of . As can be noted, the differentiation in the is quite constant as varies from to . Thus we select in all of our simulations. Indeed this guarantees a good degree of variation with the least computational cost. Finally, Fig. 2(b) confirms that the dynamical and structural properties of the emerging ES minimal structure do not show significant fluctuations when the value of is increased.

Note that the state space of the initial phases of many oscillators is a high-dimensional space (i.e. the aggregate N-dimensional state space obtained collecting the phase of each oscillator in the stack vector ). To obtain effective samplings from that space, we adopted a Latin Hypercube Sampling (LHS) strategy first proposed in McKay et al. (1979). LHS is a statistical method for generating a sample of plausible collections of parameter values from a multidimensional distribution. Specifically, let denote a variate random variable with probability density function for . Then the range space of each of the components of is partitioned in disjoint intervals of size . Taking the Cartesian product of these intervals yields cells each of probability size . Each cell can be labeled by a set of coordinates where is the interval number of component represented in cell . A LHS is obtained from a random selection of the cells , with the condition that for each the set is a permutation of integers . As a result, one random observation is made in each cell. The main advantage of the LHS strategy is that it does not require more samples for more dimension of the range space . This is the main reason why we use LHS in our method.

To measure the synchronization performance of a ES network, we consider an ensemble of phase oscillators connected by that network and evaluate Kuramoto order parameter as .

## Iv Emergence of Minimal networks

We first test our strategy by applying it to a small size network with and (Fig. 3). We obtain the matrix visualized in Fig. 2(a). In Fig. 2(b), as the threshold value is increased, the number of edges, , rapidly decreases while the value of the order parameter remains near unity.

In the figure, the normalized number of edges, which is divided by maximum links between nodes, i.e. , is plotted. Also above a certain threshold the network becomes disconnected. Therefore we choose obtaining the minimal ES network depicted in Fig. 2(c) which is characterized by edges and . We compare the minimal ES structure with the optimal network structure shown in Fig. 2(d) obtained from an exhaustive search and a Monte Carlo based method Gorochowski et al. (2010) maximizing the value of with the constraint that the total number of edges is equal to 7. We notice that the two networks share the same links.

Next, we study how heterogeneity induces functional structural properties of the network. Figure 4 shows the matrix as a function of the heterogeneity parameter when . We see that as is increased a differentiation becomes more and more apparent in the distribution of the link activation probabilities with edges between oscillators with relatively different frequencies becoming more likely to occur in the minimal ES structure.

Fig. 5(a) shows the standard deviation of the link activation probabilities as a linear function of in a larger network of oscillators. The structural properties of the emerging network are therefore induced by the node heterogeneity. This is confirmed in Fig. 5(b) where the maximum and minimum values of the node degree , corresponding to each minimal ES network, is plotted as a function of . The behaviors of the maximum value of (red dashed line) and the minimum of (black solid line) show an abrupt transition when passing from to . This suggests that the differentiation in the degree distribution of the minimal ES network becomes remarkable when heterogeneity in the nodes is increased from zero (identical oscillators) to a value greater than zero (non-identical oscillators).

The structural properties of the emergent minimal ES network are highlighted in Fig. 5(c)-(f) for a network of highly heterogeneous oscillators (). The activation probability of each link is plotted in Fig. 5(c) against the value of the difference between the natural frequencies of the oscillators at the endpoints. Links connecting more distant nodes tend to be activated with a higher likelihood confirming that differentiation among links is induced by heterogeneity in the nodes. Also, as shown in Fig. 5(d), hubs tend to be associated with oscillators whose frequencies are farther away from the average. The functional advantage of the emerging network is shown in Fig. 5(e). Indeed, we observe that the order parameter of the minimal ES structure is close to its maximal value for an all-to-all network of the same size, even if the number of links in the minimal ES network is remarkably lower than that in an all-to-all configuration. For the sake of comparison, the values of for a randomly generated network of the same number of edges is also depicted in Fig. 5(e). The sudden dip of is due to the graph becoming disconnected beyond that critical value of the threshold .

Notice that, as shown in Fig. 5(f), as the coupling strength is varied, the order parameter of the phase oscillators interconnected by the minimal ES network exhibits a sudden hysteretic change, associated to a discontinuous phase transition, whereas the system with a unimodal frequency distribution undergoes a continuous phase transition Kuramoto (1984). This discontinuous phase transition, also known as “explosive synchronisation”, has been studied in the literature Pazó (2005); Gómez-Gardeñes et al. (2011); Leyva et al. (2013); Zhang et al. (2013), also in the case of adaptive networks Hui-Jun et al. (2011); Zhang et al. (2015), revealing that the correlation between natural frequencies and the node degree, as shown in Fig. 4(d), can induce this phenomenon. Here, we wish to emphasise that the proposed evolutionary strategy, which functionally organizes the network structure for synchronization, changes the type of phase transition that would be generically observed otherwise, inducing explosive synchronisation.

Our results clearly show the role of node heterogeneity in inducing functional structures using an evolutionary strategy for network synchronization. In particular, differences in the node dynamics do influence the evolution of the network determining a differentiation in the link activation probabilities that is instrumental to obtain minimal structures with relatively high values of the order parameter. Also, hubs tend to emerge there where the distance from the average natural frequency is highest. Further simulations also confirmed that a similar structure of the emergent network can be induced by using a power-law rather than a normal distribution when selecting the heterogeneous natural frequencies of the oscillators (data not shown).

It is notable that the presence of hubs seems to characterize the emergent networks for synchronization when the nodes are heterogeneous as opposed to more homogenous structures, such as entangled networks, which have been suggested to be optimal structures in the homogeneous case Donetti et al. (2005). This is also confirmed in the case of Monte Carlo based optimal networks in Yanagita and Mikhailov (2010) where the presence of links between nodes with more distant frequencies is shown to be more likely and in the recent paper Skardal et al. (2014) based on the use of gradient-based methods. Here we obtain a further confirmation of these observations but via a generic local evolutionary strategy that is state-dependent and can be applied to a wider range of network synchronization and control problems.

## V Conclusions

Our results suggest that heterogeneity is the driving force determining the evolution of state-dependent functional networks. This can explain the structural properties detected in natural networks such as neural interconnections in the brain, gene regulatory networks or ecological networks where the states of the nodes typically affects the evolution of their interconnections Gross and Blasius (2008); Gorochowski et al. (2012); Avalos-Gaytán et al. (2012); Belykh et al. (2014). It can also be used in Dynamical Systems and Control theory to design state-dependent evolutionary strategies able to induce a desired collective behavior in a network of interest.

## Acknowledgments

This work was partially supported by JSPS KAKENHI Grants No. 24120708, No. 24740266, No. 25115719, and No. 26520206. FS would like to acknowledge support from the Network of Excellence MASTRI Materiali e Strutture Intelligenti (POR Campania FSE 2007/2013).

## Appendix A Description of the Edge Snapping Method

Edge Snapping DeLellis et al. (2010) is an adaptive strategy for the evolution of an unweighted network. Time-varying coupling gains are assigned to all pair of nodes and , with a second order dynamics affected by a double well potential defined as:

where is a damping coefficient, and are the states of the nodes at the endpoints of the edge . The driving force is a generic increasing function such that .

The gains’ dynamics mimics the damped motion of a particle in a one-dimensional space subject to a double-well potential as schematically outlined in Fig. 6. Indeed, in Fig. 6(a) the initial forcing is strong enough to drive the mass particle from the equilibrium at 0 (edge turned off) to the well associated to the equilibrium at 1 (edge activated). On the contrary, in Fig. 6(b) the forcing input to the edge snapping dynamics is not able to move the particle from the equilibrium at the origin. As a result of these dynamics, each coupling gains converges to either one of the equilibrium points, or .

Note that the dynamics of the gains is interdependent on the dynamics of the nodes’ state (in this letter, the dynamics of the states is given by a coupled oscillator dynamics among nodes). The resulting unweighted network is the outcome of the co-evolving dynamics of the nodes and the state-dependent network. This strategy is based on a distributed adaptive nonlinear approach and is therefore a generic decentralized approach relying only on a nonlinear potential to drive edge adaptation as explained in DeLellis et al. (2010).

In addition, this edge-snapping strategy has two parameters and , which can be used to control network evolution. Indeed the barrier of the potential between the two wells, acts as a constraint. As explained above, if the driving force is not strong enough, the edge, after a transient, will remain in the well corresponding to the absence of link. The height of the barrier can be tuned varying the parameter in the expression of the potential . The higher the barrier , the stronger the constraint.

## Appendix B Netevo

NetEvo is a computational framework designed to help understand the evolution of dynamical complex networks Gorochowski et al. (2012). It provides flexible tools for the simulation of dynamical processes on networks and methods for the evolution of underlying topological structures. To bring together simulation and evolution in a coherent way, the framework uses the idea of a supervisor, illustrated in Fig. 7. Evolution of the system is performed by the supervisor which can be viewed as a form of optimiser. This takes as input an initial topology, simulated output from the system and user defined constraints, and aims to return an optimal or enhanced topology. Changes to the system are assessed by using the performance measure -R (the opposite of the order parameter), with smaller values representing an improved performance. By default, NetEvo provides a supervisor that uses a Simulated Annealing meta-heuristic to search for near optimal configurations. This method has been shown to perform well for a wide range of problems with an unknown prior structure.

We tuned NetEvo to find an optimal structure, given an initial condition (the same used in the procedure for finding the minimal structure). Simulated annealing tends to avoid local minima (or maxima), so we could start the optimization from any random connected network structure. However, we decided to start ”near” the minimal structure, to facilitate the optimization (by near, we mean a structure obtained from the minimal structure, after rewiring about of its edges).

We note here, that it was necessary to run NetEvo several times (i.e. times), because of local maxima traps that the algorithm could not avoid. Finally, the optimal (or sub-optimal) structure is selected as the network that maximize (starting from ), among each of the results.

## References

- Darwin (1951) C. Darwin, The Origin of Species by Means of Natural Selection, or the Preservation of Favoured Races in the Struggle for Life (Oxford Univ. Press, 1951).
- Tan et al. (2014) S. Tan, J. Lu, G. Chen, and D. Hill, Circuits and Systems Magazine, IEEE 14, 36 (2014).
- Boccaletti et al. (2006) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Physics Reports 424, 175 (2006).
- Barrat et al. (2008) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, 2008).
- Kuramoto (1984) Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer, New York, 1984).
- Strogatz (2000) S. H. Strogatz, Physica D: Nonlinear Phenomena 143, 1 (2000).
- Acebrón et al. (2005) J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
- Arenas et al. (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics Reports 469, 93 (2008).
- Donetti et al. (2005) L. Donetti, P. I. Hurtado, and M. A. Muñoz, Phys. Rev. Lett. 95, 188701 (2005).
- Rad et al. (2008) A. A. Rad, M. Jalili, and M. Hasler, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, (2008).
- Yanagita and Mikhailov (2010) T. Yanagita and A. S. Mikhailov, Phys. Rev. E 81, 056204 (2010).
- Gorochowski et al. (2010) T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, Phys. Rev. E 81, 056212 (2010).
- Tanaka and Aoyagi (2008) T. Tanaka and T. Aoyagi, Physical Review E 78, 046210 (2008).
- Skardal et al. (2014) P. S. Skardal, D. Taylor, and J. Sun, Phys. Rev. Lett. 113, 144101 (2014).
- Gutiérrez et al. (2011) R. Gutiérrez, A. Amann, S. Assenza, J. Gómez-Gardeñes, V. Latora, and S. Boccaletti, Phys. Rev. Lett. 107 (2011).
- DeLellis et al. (2010) P. DeLellis, M. di Bernardo, F. Garofalo, and M. Porfiri, Circuits and Systems I: Regular Papers, IEEE Transactions on 57, 2132 (2010).
- McKay et al. (1979) M. D. McKay, R. J. Beckman, and W. J. Conover, Technometrics 21, 239 (1979).
- Pazó (2005) D. Pazó, Phys. Rev. E 72, 046211 (2005).
- Gómez-Gardeñes et al. (2011) J. Gómez-Gardeñes, S. Gómez, A. Arenas, and Y. Moreno, Phys. Rev. Lett. 106, 128701 (2011).
- Leyva et al. (2013) I. Leyva, A. Navas, I. Sendiña-Nadal, J. A. Almendral, J. M. Buldú, M. Zanin, D. Papo, and S. Boccaletti, Sci. Rep. 3 (2013).
- Zhang et al. (2013) X. Zhang, X. Hu, J. Kurths, and Z. Liu, Phys. Rev. E 88, 010802 (2013).
- Hui-Jun et al. (2011) J. Hui-Jun, W. Hao, and H. Zhong-Huai, Chinese Physics Letters 28 (2011).
- Zhang et al. (2015) X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Phys. Rev. Lett. 114, 038701 (2015).
- Gross and Blasius (2008) T. Gross and B. Blasius, J. R. Soc. Interface 5, 259 (2008).
- Gorochowski et al. (2012) T. E. Gorochowski, M. di Bernardo, and C. S. Grierson, Complexity 17, 18 (2012).
- Avalos-Gaytán et al. (2012) V. Avalos-Gaytán, J. A. Almendral, D. Papo, S. E. Schaeffer, and S. Boccaletti, Phys. Rev. E 86, 015101 (2012).
- Belykh et al. (2014) I. Belykh, M. di Bernardo, J. Kurths, and M. Porfiri, Physica D: Nonlinear Phenomena 267, 1 (2014).