Phase Locking Induces Scale-Free Topologies in Networks of Coupled Oscillators

An initial unsynchronized ensemble of networking phase oscillators is further subjected to a growing process where a set of forcing oscillators, each one of them following the dynamics of a frequency pacemaker, are added to the pristine graph. Linking rules based on dynamical criteria are followed in the attachment process to force phase locking of the network with the external pacemaker. We show that the eventual locking occurs in correspondence to the arousal of a scale-free degree distribution in the original graph.


Introduction
Many biological, neural, chemical and technological systems find a suitable representation as growing networks of interconnected dynamical units [1][2][3]. It is nowadays established that such networks, as they are observed in their mature state, are characterized by specific topological features, as, for instance, relatively small characteristic distances between any two nodes, high clustering properties, and fat tailed shapes in the distribution of their connectivities. In particular, most of them are known to exhibit the so-called scale-free (SF) property [4] consisting in the fact that they feature a power-law distribution P(k),k 2c in the node connectivity k (degree). From the other side, it is observed that their complex functioning and organization is often associated with the adjustment of all (or a relevant portion of) their dynamical components into a collective synchronized motion.
It is therefore crucial to understand the intimate relationship between the topological structure displayed in the resulting graphs, and the mechanisms leading to the arousal of such synchronized behaviors. For instance, recent studies have shown that i) the ability of a graph to give rise to a synchronous behavior can be greatly enhanced by exploiting the topological structure emerging from independent statistically driven growth processes [5][6][7]; ii) proper topological mechanisms of rewiring/decoupling can enhance the arousal of a synchronized behavior [8]; iii) a dynamical evolution of the underlying topology of a graph is eventually able to stabilize a synchronous motion also in those cases in which synchronization would be prevented in static graph configurations [9][10][11].
However, the question of how node dynamics can shape the network has not been extensively addressed. The issue has only been considered within the game theory framework [12,13], where a not growing network of players is shaped in a decision game. In this Letter, we show that a growing process entirely guided by dynamical criteria to force frequency and phase locking of an original set of networking oscillators is able to fully reshape the connectivity of the graph, and that the entrainment process induces a scale-free degree distribution in the pristine network.

Methods
Without lack of generality, we exemplify our discussion with reference to an initial t = 0 graph G 0 of bi-directionally coupled Kuramoto phase oscillators [14,15] where i runs from 1 to n 1 , {v 0i } is the set of natural frequencies of the phase oscillators, uniformly distributed within the range 0.560.25, k i (t = 0) is the initial degree of (the number of connections pertinent to) the i th oscillator, d 1 is a coupling constant, and the set {a ij } are the elements of the n 1 6n 1 adjacency matrix A = (a ij ) describing the structure and topology of the wiring of connections in G 0 (a ij = 1 if oscillators i and j are connected, while a ij = 0 otherwise). Time integration is here performed by means of the Heun method with an integration step Dt in = 0.1. Though the main results of our study are independent on the initial choice of G 0 , from the time being we generate G 0 by means of the model introduced in Ref. [16]. Precisely, we start with a ring lattice of n 1 nodes, each one connected to its k 0i = 2m 0 nearest neighbors, and obtain G 0 by randomly rewiring each link with probability p = log n 1 /n 1 . The resulting G 0 approximates the so called small world property, featuring an average shortest path length / log n 1 , and an exponentially decaying degree distribution P 0 (k) that is well peaked around the mean value AEkae = 2m 0 . Furthermore, the coupling constant d 1 is selected so as the initial graph does not display a phase synchronized motion.
In our generic trials, the initial network G 0 is let evolve in its unsynchronized motion from t = 0 to t = 30 time units, and, at subsequent times t l = t 0 +lDt, a forcing network is grown on top of the evolution of G 0 . Namely, at each time interval Dt (which will be taken as a multiple of the integration step Dt in , i.e. Dt = sDt in ), a new node is added, forming m connections with nodes in G 0 . The added nodes are identical phase oscillators that follow the instantaneous phase of an external pacemaker _ w w p~vp , and all newly formed links correspond to unidirectional (forcing) interactions to nodes in G 0 . Notice that such forcing process is fully equivalent to considering a unique external forcing node (the pacemaker) which forms successive (and possibly multiple) connections with oscillators in G 0 in the attempt of locking their frequencies, and, as so, it can be taken as a representation of phenomena occurring in social science [the emergence of consensus driven by a leader opinion (media, press, fashion, publicity, etc..)] or in biological systems (the entrainment of circadian rhythms driven by the main circadian clock).
The key point is the selection mechanism through which the added nodes are linked to G 0 . We here consider a dynamical criterion fully driven to enhance phase locking: when the l th new node is attached to G 0 , it forms m connections preferentially with those nodes in G 0 whose instantaneous phase at time t l , w j (t l ), holds more closely a given phase condition. Specifically, we consider a generic parameter dM(0,2p) and establish the first of the m connections with that node j l whose actual value of the phase holds the condition min j~1,...,n1 d{Dh j mod2p with When m.1, we iteratively repeat the same condition excluding those nodes that already received a link at the same time step. In the following, we will set m = 1. However, the reported scenario is independent on the specific choice of m, the only difference being that for m = 1, the added nodes do not form additional cycles nor loops in the original graph G 0 . As for the parameter d, we will show that it will not affect qualitatively the reported scenario. The only constrain is that it cannot be taken equal to 0 nor to 2p, as these values correspond to the stable fixed point emerging during the locking of a single phase oscillator, and therefore these settings would determine a situation in which the first node in G 0 becoming locked with the pacemaker, would, from there on, attract all the rest of connections. Each newly formed connection is assigned a coupling coefficient d p and a coupling direction from the added node to the selected node in G 0 (different kinds of coupling interactions have been also studied, and a detailed report will be presented elsewhere [17]). The evolution of the graph is therefore described by where k i (t) is the time evolving degree of the i th node that accounts for new connections the i th node is receiving from added nodes, and the matrix B = (b ij ) is a size evolving matrix of n 1 6l(t) elements (with l(t)#n 2 ), whose entries b ij are equal to 1 if the j th added node formed a connection with the i th node in G 0 , and zero otherwise. Figure 1 reports the quantity R = AER(t)ae t vs. the parameter space (v p , d p ) for n 1 = 100, d 1 = 0.2, n 2 = 200, s = 1, and d = p (anti-phase coupling condition). Here AE…ae t denotes an average over time (performed after the growing process is finished), and R(t) stands for the phase synchronization order parameter [14] R t ð Þ~1 n 1 X n1 j~1 e iw j t ð Þ .

Results
From Figure 1 it is evident that the threshold for the setting of the phase locking of R(t) (R.1) crucially depends on the frequency of the external pacemaker v p . Specifically, as far as v p is close to v = 0.5 (the average frequency of the oscillators in G 0 ), the locking process occurs already for a relatively small value of d p , whereas, as v p deviates significantly from v the value of d p producing phase locking becomes larger and larger. A more quantitative description of the process can be gathered by inspection of Figure 2, where we report the time evolution of the mean frequency of the oscillators , and of R(t), for three pacemaker frequencies (v p = 0.1, 0.5 and 0.9) and two distinct coupling strengths (d p = 0.5 and d p = 5.5). It is seen that, while the low coupling regime is not associated to a phase locking of G 0 with the pacemaker, in the high coupling regime v(t) converges (after the growing process has ended) to the external forcing frequency and, at the same time R(t) converges to unity, and s v (t) vanishes.
The results reported in Figures 1 and 2 are generic, and the same qualitative scenario (though for different values of d p and v p ) characterizes the evolution of the system at different system sizes, for different initial topologies in the connectivity of G 0 , for different values of dM(0,2p) and for different values of m.
Let us now move to report the central point of this study, related to the investigation of the peculiar topological structures induced in the final degree distribution of G 0 as a result of the phase locking process. For instance, in Figure 3 we depict the final graphs obtained for n 1 = 100 and n 2 = 1000, with the nodes of the original graph (the added nodes) depicted in black (blue). The left and right networks correspond to very different situations. In the left network, the growing process is unable to produce locking of the phases of the oscillators in G 0 to the external pacemaker, and one sees that the distribution of blue attachments is rather homogeneous. At variance, the right graph corresponds to a case in which the forcing nodes eventually succeed to entrain the phases of the oscillators in G 0 , and eye inspection reveals a high heterogeneity of the blue attachments, with the simultaneous presence of few hubs (nodes with very high degree) and many seldomly connected nodes. In order to describe more quantitatively the situation, we perform large trials of numerical simulations with n 1 = 1000, n 2 = 10000, and d 1 = 0.2, and we monitor the time evolution of the degree distribution P t (k) of all nodes originally belonging to G 0 during the process of locking. In fact, we here measure the cumulative degree distribution P c t k ð Þ, given by P c t k ð Þ~P kmax k 0 P t k 0 ð Þ. This is because the summing process of the P(k) smooths the statistical fluctuations generally present in the tails of the distribution. As a generic property, it is important to remark that, if a power law scaling is observed in the behavior of P c (k) (i.e. if P c k ð Þ*k {c c ), this implies that also the degree distribution P(k) is characterized by a power law scaling P(k),k 2c , with c = 1+c c . Figure 4 reports how P c t k ð Þ evolves in two distinct situations: for a process that does not lead to any locking (left panel) and for a process that eventually leads to locking of G 0 to the frequency of the pacemaker (right panel). In the former case, while P c t k ð Þ deviates significantly from P c 0 k ð Þ in the course of time, it never assumes a power-law shape, while in the latter case the locking process (manifested by the evolution of R(t) to 1 in the inset) is accompanied by the convergence of P c t k ð Þ to an asymptotic distribution P c fin k ð Þ which features a power-law shape. The difference in the final distributions for the non locked and locked networks, and the convergence in this latter case of P c t k ð Þ to a SF distribution P c fin k ð Þ is a generic feature in the parameter    Figure 5(c)]. In all cases, solid (dashed) lines correspond to the locked (non locked) regime, obtained for high (low) values of d p , and solid red lines indicate the best power-law fits. Whenever the forcing nodes eventually induce the locking of G 0 to the pacemaker frequency, the final degree distribution displays a power-law (scale-free) behavior P(k),k 2c . The specific slope c of the power law scaling depends on the specific choice of the external frequency v p . In our trials, we always observed values of c in the range (2,3), in accordance to the values measured for most of the real world networks [1][2][3].

Discussion
The formation of a scale free degree distribution can be qualitatively understood as follows. Each node in G 0 can be in two different states: a state in which its phase is locked with that of the pacemaker (corresponding to Dh = 0) and an unlocked state, in which Dh(t) rotates in the unit circle clockwise or counterclockwise, depending on whether the difference between the frequency of the node and that of the pacemaker is positive or negative. Having set  d strictly different from 0, this implies that, once a node reaches the locked state, it cannot be further attached by the pacemaker. As a given node is connected to the pacemaker, it is reasonable to assume that its phase dynamics will be slightly attracted to that of the pacemaker. Initially, all nodes are not locked, and therefore the probability to find their phase differences in the proximity of the attachment condition will not depend on their specific initial frequency. When d p is small enough so as not to determine locking of the network, only the very few nodes whose original frequency was close to that of the pacemaker will reach the locked state, whereas all the remaining nodes will evolve unsynchronized, and for them we will have a sort of random attachment with equal probability leading to a final degree distribution which substantially differs from a scale free network (Figures 5(a)-(c), dashed lines). On the contrary, at high values of d p , nodes will be sequentially attracted to the locked state, starting again from those whose original frequency was close to that of the pacemaker. In this latter situation, the more nodes gets the locked state (where any further attachment is prevented), the higher is the probability that the remaining unlocked nodes will get subsequent connections, resulting in a kind of preferential attachment process in which the higher is the absolute value of the difference between the original frequency of a node and that of the pacemaker, the higher is the chance that this node will get connections during the growing process. This is confirmed in Fig. 5(d), where we report the final number of connections k i (t fin ) acquired by each node as a function of its initial frequency v 0i for v p = 0.5 and two values of d p . There, one clearly see that, while in the unlocked regime (upper plot) the distribution is almost uniform reflecting a random-like attachment, in the locked regime (lower plot), the distribution is strongly nonuniform, and promotes those nodes with higher frequency mismatch with the pacemaker. Notice that the few points in Figure 5(d) of high degree in proximity to v p correspond to oscillators whose initial frequency was, indeed, close to v p , but whose instantaneous frequency at the instant at which the growing process starts was moved away from v p during the initial, unforced, evolution of G 0 .
In conclusion, we have shown that a growing process entirely guided by dynamical criteria to promote and phase locking of an original network is able in a robust way to fully reshape the connectivity of the graph, and that the locking process is associated with the emergence of a scale-free degree distribution in the network connectivity. This fact can therefore give new hints on the fundamental processes that rule the growth of some of the real world networks, that ubiquitously feature such kind of connectivity distributions.