Rich-club network topology to minimize synchronization cost due to phase difference among frequency-synchronized oscillators

# Rich-club network topology to minimize synchronization cost due to phase difference among frequency-synchronized oscillators

Takamitsu Watanabe 111takawatanabe-tky@umin.ac.jp Department of Physiology, School of Medicine, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
###### Abstract

Functions of some networks, such as power grids and large-scale brain networks, rely on not only frequency synchronization, but also phase synchronization. Nevertheless, even after the oscillators reach to frequency-synchronized status, phase difference among oscillators often shows non-zero constant values. Such phase difference potentially results in inefficient transfer of power or information among oscillators, and avoid proper and efficient functioning of the network. In the present study, we newly define synchronization cost by the phase difference among the frequency-synchronized oscillators, and investigate the optimal network structure with the minimum synchronization cost through rewiring-based optimization. By using the Kuramoto model, we demonstrate that the cost is minimized in a network topology with rich-club organization, which comprises the densely-connected center nodes and peripheral nodes connecting with the center module. We also show that the network topology is characterized by its bimodal degree distribution, which is quantified by Wolfson’s polarization index. Furthermore, we provide analytical interpretation on why the rich-club network topology is related to the small amount of synchronization cost.

###### pacs:
89.75.-k, 05.45.Xt, 88.80.H-

## I Introduction

As power grids Kundur1993Elgerd1982Book (); Dorfler2010ACC () and networks of bursting neurons Varela2001NatRevNeuro (); Fell2011NatRevNeuro (), functions of some complex networks of oscillators are based on not only synchronization of frequencies of the oscillators, but also synchronization of their phases. However, in general, frequency synchronization is more achievable than phase synchronization. Phase difference among frequency-synchronized oscillators often falls into a non-zero constant, and such non-zero phase difference would avoid proper and efficient functioning of the complex networks.

In power grids, alternating voltage of the power plants in the grids should be synchronized around certain specific frequency (e.g., 50 Hz in the most parts of Europe and 60 Hz in the north America) Kundur1993Elgerd1982Book (), and disruption of the frequency synchronization causes a blackout in a large area Kundur1993Elgerd1982Book (); HillChen2006IEE (). In addition, the phases of the voltages of the power plants are also required to be synchronized. As discussed in Appendix A, the difference in voltage phases among power plants inevitably causes power loss consumed as heat in power lines Kundur1993Elgerd1982Book (); Dorfler2010ACC (). In this sense, the phase difference in power grids can be regarded as synchronization cost. Considering recent increasing environmental awareness and soaring global demand of natural resources USEFC2008 (); SmartGrid2011IEEE (), it is necessary to reduce the synchronization cost due to the phase difference among frequency-synchronized power plants.

Synchronization in large-scale brain networks also requires less difference in phase of neuronal activity among different brain regions. A series of prior experimental researches have shown that various important functions in large-scale brain networks are based on not only frequency- but also phase- synchronizationVarela2001NatRevNeuro (); Fell2011NatRevNeuro (). A previous electrophysiological study showed that spike activity recorded from monkeys’ cortex exhibited phase synchronization in various frequency bands while the monkeys were conducting tasks that required integration of visual processing and motor responses Bressler1993Nature (). In another careful electrophysiology study, Roelfsema and his colleagues recorded local field potentials (LFP) in the cerebral cortex in cats, and revealed that phase synchronization between LFPs recorded in the visual and parietal cortices was increased only when the cats focused their attention on stimuli Roelfsema1997Nature (). Studies using electro-encephalogram to record human brain activity found that increase of phase synchronization in various frequency was associated with learning and perception of images Miltner1999Nature (); Hipp2011Neuron (). Furthermore, the frequency- and phase- synchronization are considered to occur in large-scale brain networks Bressler1993Nature (); Traub1996Nature () with zero time-lag Roelfsema1997Nature (); Womelsdorf2007Science (); Fell2011NatRevNeuro (). These researches suggest that phase synchronization in large-scale brain networks is related to crucial functions such as integration of multiple information Varela2001NatRevNeuro (), neural communication Fries2005TrendsCogSci (), and spike-timing-dependent plasticity Fell2011NatRevNeuro (). Actually, it is known that some types of the disruption of the synchronization cause dysfunctions of memory learning Stopfer1997Nature () or psychiatric disorders Tononi2000RBRR (). Considering these findings, it is to some extent reasonable to hypothesize that large-scale brain networks have a specific organization that minimizes phase differences among brain activity, and enables optimal phase synchronization in the entire networks.

These previous literatures indicate the importance to reduce phase difference among frequency-synchronized oscillators, which can be regarded as a type of cost that is spent during synchronization. However, little is examined about the optimal network topology that reduces this type of synchronization cost due to phase difference. Indeed, a series of previous literatures have investigated optimal network structures by introducing a different type of synchronization cost, which is needed for building or maintain of the optimal network infrastructure. A study regarded parameters based on coupling strength among oscillators as a cost, and demonstrated that homogeneous and uniform distribution of the coupling strength enhances the tendency of synchronization Motter2005PRE (). Another study employing the same definition of synchronization cost suggests that more heterogeneous network structures are required to simultaneously achieve both the maximum synchronizability and minimum synchronization cost based on coupling strength Nishikawa2006PhyD (). Another study revealed the optimal distribution of coupling strength that increases synchronizability among Kuramoto oscillators Brede2008EPJB (). However, these prior researches have not focused on synchronization cost due to phase difference. Consequently, despite a line of prior researches on optimal conditions for synchrony in networks Arenas2008RhysRep (), it remains unclear about optimal network structures that minimize synchronization cost due to phase difference among oscillators.

In the present study, therefore, we examine the optimal network topology to minimize the phase-difference-based synchronization cost. We take the following five steps:

1. First, we define the synchronization cost, , due to phase difference between frequency-synchronized phase oscillators and in the Kuramoto model. We adopt the Kuramoto model because the model has been used as approximation of various systems including power grids Buzna2009PRE (); Fioriti2009CIIS (); Dorfler2010ACC () and large-scale brain networks Breakspear2010FrontHN (); GomezGardenes2010PLosOne ().

2. Second, by using the definition, we numerically calculate the mean of synchronization cost, , in the entire network.

3. Third, by using a rewiring strategy Donetti05PRL (), we show that ”rich-club” network topology Zhou2004IEEEComm (); Colizza2006NatPhys (); GomezGardenes2010PLosOne (); Heuvel2011JNS (), which consists of densely inter-connected modules and peripheral low-degree nodes, is the optimal topology with the minimum synchronization cost.

4. Forth, we characterize the rich-club network topology by quantifying the bimodality of the degree distribution of the network.

5. Finally, we provide an analytical interpretation on why the rich-club network organization is associated with a small amount of the synchronization cost.

## Ii Method

### ii.1 Definition of synchronization cost

We first define synchronization cost due to phase difference in an unweighted and undirected network which is described by an adjacency matrix and consists of phase oscillators. is when oscillators and are connected, and is when they are not. According to the Kuramoto model, the phase of the oscillator , , is described as

 ˙θi=ωi−ϵ∑jAijsin(θi−θj), (1)

where is the natural frequency of the oscillator , and is a coupling strength. To reduce computational cost for the following rewiring-based optimization, we assume that is constant for any combination of oscillators. In this Kuramoto model, the synchronization cost between oscillators and , is defined based on the phase difference between the connected oscillators as follows:

 Sij=(θi−θj)2. (2)

Note that the synchronization cost is only defined after the network of the oscillators reaches to a state of frequency synchronization. As described in Appendix A, in power grids, can be regarded as an index that is proportional to power loss due to difference in voltage phase between power plants and .

### ii.2 Estimation of the mean synchronization cost

Based on the definition of , we numerically estimate for each edge in the following four steps for a given network:

1. We set normally-distributed for each node. It is because that previous studies on real networks such as power grids and brain networks have assumed that the natural frequencies of belonging oscillators are symmetrically fluctuating around the averaged frequency Filatrella2008EPJB (); Fioriti2009CIIS (); Fell2011NatRevNeuro ().

2. Based on the Kuramoto model described in Eq. (1), we numerically estimate frequency-synchronized status, where becomes a common constant value, , for any .

3. In the frequency-synchronized status, each oscillator has a different specific phase, . Based on the set of , we then evaluate for each edge.

4. As described in (i), the set of the natural frequency, , is fluctuating over time. Thus, the is also fluctuating over time. Therefore, we repeat the procedure (i)-(iii) 200 times with different sets of , and obtain 200 different sets of . Then, we average the over time, obtaining for each edge. Finally, we average the across edges, and obtain for the entire network.

### ii.3 Rewiring-based optimization

To search for the optimal network topology with the least , we adopt the rewiring method that previous studies used to explore the network topology with the largest synchronizability Donetti05PRL (). We apply the following rewiring-based optimization procedure to a given connected network with nodes and mean degree of : At each step, the number of rewired edges is randomly determined based on an exponential distribution. The set of edges to be rewired is also randomly chosen in a given network. After the rewiring, we estimate frequency-synchronized status and obtain . The attempted rewiring is rejected if the updated network is disconnected. Otherwise, the rewiring is accepted if , or with probability where is a temperature-like parameter and Donetti05PRL (). The initial rewiring is conducted at , and, after the first rewiring, is set as where is the largest in the first rewiring trials. After that, is decreased in every rewiring trials. This estimation process iterated until there is no change in more than 50 successive rewiring steps. We apply this rewiring-based optimization to three different initial networks: an Erdős-Rényi (ER) random model, Watts-Strogatz (WS) model Watts98Nature (), and a Barabási-Albert (BA) model Barabasi99Science () with and Donetti05PRL (). In all the cases, the coupling strength, , is set to be . Each set of the natural frequencies of Kuramoto oscillators, , is randomly chosen from the normal distribution with an average of and a standard deviation of .

During the optimization, we trace the standard order parameter, , and local synchoronizability, Gomez-Gardenes2007PRLPRE (), defined as follows:

 reiψ=1N∑jeiθj, (3) rlocal=12Nl∑i∑j∈Γi∣∣∣limΔt→∞1Δt∫tr+Δttrei[θi(t)−θj(t)]dt ∣∣∣, (4)

where is the total number of edges, is the set of neighbors of node , and is a global phase. Furthermore, after the optimization is completed, we compare the optimized networks from the different initial networks by estimating the following basic topological properties: mean of shortest path length Bunde1991Book (), mean of clustering coefficient Watts98Nature (), mean of betweenness centrality Newman2001PRE (), and degree correlation Newman2001PRE (). We conducted ten optimizations of ten different networks for each type of the initial networks, and averaged these basic properties.

• The shortest path length, , is defined as the shortest distance between two nodes and Bunde1991Book (). The averaged shortest path length, , is defined as the average value of over all the possible pairs of nodes in the network.

• The clustering coefficient for node , , measures the local group cohesiveness Watts98Nature (), which is defined as the ratio of the number of links between the neighbors of and the maximum number of such links. We define the mean clustering coefficient, , as an average value over all the nodes.

• The betweenness centrality for node , , is defined as the number of shortest paths between pairs of nodes that pass through a given node Newman2001PRE (). We define the mean betweenness centrality, , as an average value over all the nodes.

• The degree correlation for a network is defined as the Pearson assortativity coefficient of the degrees, Newman2001PRE (). The coefficient enables us to quantify the preference for high-degree nodes to attach to other high-degree nodes. Networks with this preference show large .

### ii.4 Estimation of rich-club coefficient

We estimate reich-club coefficient, , for both initial networks and optimized networks. According to the previous studies Zhou2004IEEEComm (); Colizza2006NatPhys (); GomezGardenes2010PLosOne (); Heuvel2011JNS (), the coefficient for each degree is calculated as

 Φ(k)=2E>kN>k(N>k−1), (5)

where represents the number of edges among nodes that have more than degrees. As in previous studies Colizza2006NatPhys (); GomezGardenes2010PLosOne (); Heuvel2011JNS (), we calculate normalized rich-club coefficients, through dividing the raw value, , by the mean of rich-club coefficients of 100 random networks (ER models), , as follows,

 Φnorm(k)=Φ(k)Φrandom(k). (6)

## Iii Results

### iii.1 Rewiring-based optimization

Fig. 1 shows representative results of the rewiring-based optimizations. was decreased from approximately to even when the initial network structure was different. In all the initial networks, reached to a stable status after approximate 400 steps of rewiring. Strikingly speaking, we cannot guarantee that the optimal network was found, but this result suggests that a reasonably robust approximation of the optimal topology was obtained in this method.

During the optimization, the standard order parameter, , were fluctuating just below 1 during the optimization (a small panel in Fig. 1 A). The local synchonizability, , showed the similar fluctuation below 1. Considering the previous studies on these parameters Gomez-Gardenes2007PRLPRE (), these behaviors of and are considered to be related to the amount of the coupling strength, . The previous studies Gomez-Gardenes2007PRLPRE () have demonstrated that, when the coupling strength is more than 0.2, both of and reach a plateau that is near to 1 regardless of network topology. In the present study, the coupling strength, , was set at 0.3, because global synchronization is necessary for the estimation of . This relatively large coupling strength could result in the saturation of the global and local order parameters, and during the optimization process.

Fig. 1 B shows that the optimized networks for the three different initial networks commonly exhibit a characteristic topology, which has a densely interconnected core nodes and peripheral nodes dangling the core module. The heterogeneous network features were also observed in the basic network properties in the optimized networks (Tab. 1). Compared with the initial networks, the optimized networks tended to show larger averaged values of the shortest path length, , betweenness centrality, , and degree correlation, . The averaged values of the clustering coefficients, , were smaller in the optimized networks. These results suggest that, through the optimization process, the network seems to enlarge its heterogeneity.

### iii.2 Rich-club organization

This heterogeneous network topology has been reported as ”rich-club” organization in a series of previous theoretical and experimental studies Zhou2004IEEEComm (); Colizza2006NatPhys (); GomezGardenes2010PLosOne (); Heuvel2011JNS (). The prior literatures have characterized the organization by estimating normalized rich-club coefficients, described in Eq. 6. If the network has rich-club organization, should be more than , and increase monotonically as increases.

Fig. 2 shows the comparison in between initial networks and optimized networks. To clarify the difference in the rich-club coefficients, we adopted as the initial networks larger networks than shown in Fig. 1 B (i.e., ). Before the optimization, was not always larger than (e.g., ER and WS models), and did not show monotonic increase along , which is consistent with a previous study Colizza2006NatPhys (). In contrast, in the optimized networks, the rich-club coefficients were larger than in almost all the range of , and monotonically increased as increased. These phenomena were observed commonly among the three different optimized networks that were derived from the three different initial networks. In addition to the appearance of the networks in Fig. 1 B, this estimation of supports the notion that the networks with rich-club organization has the minimum or a very small amount of synchronization cost due to phase difference among frequency-synchronized oscillators.

### iii.3 Bimodal Degree Distribution

As shown in Fig. 1 B, the rich-club network topology consists of high-degree nodes cluster and low-degree peripheral nodes. Therefore, we hypothesized that the topology can be characterized by a bimodal degree distribution. To test the hypothesis, we estimated Wolfson’s polarization index, Wolfson1997RIW (). The Wolfson’s index for degree distribution is defined as:

 (7)

where is the mean of the degree, , and denotes the median of the degree. and are the mean values of and of , respectively. represents Gini inequality index, which is defined as . This Wolfson’s polarization index shows the extent of the bimodality of the distribution. If the distribution is completely the same as a uniform distribution, the is . If the half of population has nothing and the other half shares everything, the reaches a maximum, . In the present case, a larger indicates that the network has a more bimodal and bipolarized degree distribution.

We estimated during the rewiring-based optimization. Because can be calculated more accurately for networks with more nodes, we used the BA model with and as an initial network for the optimization. As a result, in the course of the above-mentioned optimization, decreased during the rewiring-based optimization (Fig. 3A). Meanwhile, as decreases, increases (Fig. 3B). Actually, the degree distribution changed from a power-law distribution (Fig. 3C) to a bimodal distribution (Fig. 3D). This relation was also observed for different initial networks (e.g., ER model). This correlation supports the hypothesis that rich-club network with small can be characterized by its bimodal degree distribution.

### iii.4 Analytical Interpretation

We finally provide an analytical interpretation on why the rich-club network has less . Using mean-field approximation, the Eq. (1) in the frequency-synchronized status can be described as , where is defined in the Eq.  (3). Therefore, if is small enough, , and is described as

 (θi−θj)2=(ωi−Ω)2ϵ2ki2+(ωj−Ω)2ϵ2kj2−2(ωi−Ω)(ωj−Ω)ϵ2kikj, (8)

for a set of . is obtained as averaged across an enough number of sets of . As in the above-described numerical estimation, we assume that is distributed according to the normal distribution with a mean of and a standard deviation of , and that the synchronized frequency is always in every set of . Because we can also assume that is nearly equal to , is considered to be equal to , and is considered to be equal to zero. Consequently, we obtain the approximation of as follows:

 ~⟨Sij⟩=σ2ϵ2(1ki2+1kj2). (9)

This approximation was validated through comparison of with the real , as shown in Fig. 4 A (). The two parameters had a large negative value of Pearson’s correlation coefficient ().

This expression of gives qualitative explanation on relationship between rich-club network topology and a small value of : To achieve a small amount of , should be small. When the node has a high degree, it makes more contribution to a smaller to connect with a the node with a high degree. It is also case when the node has a small degree. As a result, for a small value of , high degree nodes should connect with other high degree nodes, and low degree nodes should not have edges with other low degree nodes, but with high degree nodes. Consequently, high degree nodes tend to be gathered and create a densely-connected core module, and low degree nodes tend to connect with high degree nodes in the core module. Overall, the optimized network with a small amount of is likely to be a rich-club network.

Note that it is difficult to further extend this approximation. If this approximation of is accurate enough, a simple calculation of Eq. 9 leads us to the proportional relationship between and . Given is a constant value as in the present study, should be proportional to . However, as shown in Fig. 4 B, we could not observe a linear relationship. This inaccurate approximation of may be caused by accumulation of the small difference between and . This result suggests that we cannot extend this approximation to represent only by .

## Iv Discussion

The present study introduced synchronization cost based on phase difference among frequency-synchronized oscillators. Using the rewiring-based optimization Donetti05PRL (); Gorochowski2010PRE (), we showed that the synchronization cost is minimized in a rich-club network topology. Furthermore, we demonstrated that the network topology can be characterized by the bimodality of its degree distribution. Finally, we provided analytical explanation on why the rich-club network topology is associated with a small amount of synchronization cost.

The concept of synchronization cost is not a novel idea of the present study. As described in Sec. I, a line of previous studies have investigated a different type of synchronization cost, which is based on coupling strength Motter2005PRE (); Nishikawa2006PhyD (); Brede2008EPJB (). Whereas the present synchronization cost due to phase difference can be regarded as dynamic cost per unit time, the cost based on coupling strength can be considered as static cost that is related to building and maintaining of network infrastructures. Interestingly, the optimal network topology with the least cost depends on which of the two types of synchronization cost we adopt. The optimal networks for the synchronization cost based on coupling strength often show more homogeneous properties Motter2005PRE () than those for the other synchronization cost. The homogeneity of networks is desired to enhance synchronizability Nishikawa2003PRL (); Donetti05PRL (). Therefore, it may be necessary to investigate what network structures balance these two types of synchronization cost.

The present synchronization cost in the present study can be another concept of load assigned to edge in a complex network. Previous studies used edge-betweenness as edge load Holme2002PRE (); Guimera2003PRE (), which is useful in various situations from human interaction Guimera2003PRE () to data transmission in computer networks Holme2002PRE (). However, because the edge betweenness does not consider synchrony in networks, its properties have evident difference from those of the synchronization cost. For example, as shown in Sec. III.4, the synchronization cost is lower between high degree nodes, and higher between low degree nodes. In contrast, the edge betweenness tends to be higher in edges bridging high degree nodes and be lower in edges bridging low degree nodes. These distinct properties suggest the possibility that the synchronization cost can be another concept of edge load.

The synchronization cost in the present study, , has a mathematical expression similar to that for local synchronizability, Gomez-Gardenes2007PRLPRE (). However, the two parameters focus on different phases of synchronization in complex networks. The local synchronizability enables us to quantify the local construction of the synchronization pattern. Therefore, it is useful to investigate properties of networks that are not yet fully synchronized. In contrast, the synchronization cost in the present study can be only estimated in fully-synchronized networks. Therefore, in the present study, we used a relatively large coupling strength, and achieved full synchronization throughout the optimizations. As a result, in the present study, the local synchronizability was always saturated.

Although the present study did not adopt models specific to any real networks, the findings may help understanding large-scale brain networks. Recently, a few of recent studies have reported the existence of the rich-club organization in the large-scale brain networks. A previous empirical study has demonstrated the existence of the rich-club organization in the large-scale human brain networks Heuvel2011JNS (). Another study has investigated the anatomical connectivity in the cerebral cortex of cats, and has showed that rich-club organization controls the dynamic transition of synchronization in the brain GomezGardenes2010PLosOne (). A recent review has suggested that the organization is a cost-effective network topology for the brain networks, which are required to be adapted to various cognitive functions Bullmore2012NatRevNeuro (). In addition to the context of cost-effectiveness, the rich-club network topology is robust to random attack Paul2006PhysicaA (). This previous study analytically and numerically demonstrated that networks robust to random attack have similar structures observed in the present study. The robust networks highly-interconnected hub modules and peripheral nodes (leaf nodes) that have a single edge. This network topology has bimodal degree distribution and shows rich-club organization. Considering these prior literatures, it is suggested that the rich-club organization is beneficial for the large-scale brain networks to efficiently and robustly maintain its wide range of functions based on synchronization.

###### Acknowledgements.
The author acknowledges the support from the Japan Society for the Promotion of Science (JSPS) Research Fellowship for Young Scientists (222882).

## References

• (1) O. I. Elgerd, Electric energy systems theory (McGraw-Hill Companies, 1982); P. S. Kundur, Power system stability and control (McGraw-Hill Companies, 1993).
• (2) F. Dörfler and F. Bullo, American Control Conference 2010, 930-937 (2010).
• (3) F. Varela, J.P. Lachaux, E. Rodriguez, and J. Martinerie, Nat Rev Neurosci. 2(4):229-39 (2001).
• (4) J. Fell, and N. Axmacher, Nat Rev Neurosci. 12(2):105-18 (2011).
• (5) D.J. Hill, and G. Chen, IEEE Int. Symposium on Circuits and Systems 711-715 (2006).
• (6) H. Gharavi, and R. Ghafurian, Proceedings of the IEEE 99, 917-921 (2011); G.W. Arnold, ibid. 99, 922-927 (2011).
• (7) Energy Future Coalition, Report of Smart Grid Working Group, (2008).
• (8) S.L. Bressler, R. Coppola, and R. Nakamura, Nature. 366(6451):153-6 (1993).
• (9) P.R. Roelfsema, A.K. Engel, P. König, and W. Singer, Nature. 385(6612):157-61 (1997).
• (10) W.H. Miltner, C. Braun, M. Arnold, H. Witte, E. Taub. Nature. 397(6718):434-6 (1999).
• (11) J.F. Hipp, A.K. Engel, and M. Siegel, Neuron. 69(2):387-96, (2011).
• (12) R.D. Traub, M.A. Whittington, I.M. Stanford, and J.G. Jefferys, Nature. 383(6601):621-4 (1996).
• (13) T. Womelsdorf, J-M. Schoffelen, R. Oostenveld, W. Singer, R. Desimone, A.K. Engel, and P. Fries, Science. 316(5831):1609-12 (2007).
• (14) P. Fries, Trends Cogn Sci (Regul Ed). 9(10):474-80 (2005).
• (15) M. Stopfer, S. Bhagavan, B.H. Smith, and G. Laurent, Nature. 390(6655):70-4 (1997).
• (16) G. Tononi, and G.M. Edelman, Brain Res Brain Res Rev. 31(2-3):391-400 (2000).
• (17) A.E. Motter, C. Zhou, and J. Kurths, Phys. Rev. E 71, 016116 (2005).
• (18) T. Nishikawa, and A.E. Motter, Physica D 224, 77-89 (2006).
• (19) M. Brede, Eur. Phys. J. B. 62, 87-94 (2008).
• (20) A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469, 93-153 (2008).
• (21) L. Buzna, S. Lozano, and A. Diaz-Guilera, Phys. Rev. E 80, 066120 (2009).
• (22) V. Fioriti, S. Ruzzante, E. Castorini, E. Marchei, and V. Rosato, Critical Information Infrastructure Security 14-23 (2009).
• (23) M. Breakspear, S. Heitmann, and A. Daffertshofer, Front. Hum. Neurosci. 4:190 (2010).
• (24) J. Gómez-Gardeñes, G. Zamora-López, Y. Moreno, and A. Arenas. PLoS ONE. 5(8):e12313 (2010).
• (25) L. Donetti, P.I. Hurtado, and M.A. Muñoz, Phys. Rev. Lett. 95, 188701 (2005).
• (26) S. Zhou, and R. Mondragon. IEEE Commun Lett. 8(3):180-2 (2004).
• (27) V. Colizza, A. Flammini, M. Serrano, and A. Vespignani, Nat Phys. 2(2):110-5 (2006).
• (28) M.P. van den Heuvel, and O. Sporns, Journal of Neuroscience. 31(44):15775-86 (2011).
• (29) G. Filatrella, A.H. Nielsen, and N.F. Pedersen, Eur. Phys. J. B. 61, 485-491 (2008).
• (30) D.J. Watts and S.H. Strogatz, Nature 393, 440 (1998).
• (31) A.L. Barabási and R. Albert, Science 286, 5439 (1999).
• (32) J. Gómez-Gardeñes, Y. Moreno, and A. Arenas, Phys. Rev. Lett. 98 034101 (2007), Phys. Rev. E. 75 066106 (2007).
• (33) A. Bunde, and S. Havlin, Fractals and Disordered System (Springer Verlag, 1991).
• (34) M.E.J. Newman, Phys. Rev. E. 64 016131 (2001).
• (35) M. Wolfson, Review of Income and Wealth 43, 4 401 (1997).
• (36) TE Gorochowski, M di Bernardo, CS Grierson. Phys Rev E. 81 056212 (2010).
• (37) T. Nishikawa, A.E. Motter, Y.C. Lai, and F.C. Hoppensteadt, Phys. Rev. Lett. 91, 014101 (2003).
• (38) P. Holme, Phys. Rev. E 66, 036119 (2002).
• (39) R. Guimerà, L. Danon, A. Diaz-Guilera, F. Giralt, and A. Arenas, Phys. Rev. E 68, 065103 (2003).
• (40) E. Bullmore, and O. Sporns, Nat Rev Neurosci. 336-349, (2012).
• (41) G. Paul, S. Sreenivasan, S. Havlin, and H.E. Stanley, Physica A. 370(2):854-62 (2006).

## Appendix A Definition of Synchronization Cost

In this section, we explain why power loss consumed in the electric line between power plants can be represented by square of difference in phase of voltage between the two power plants.

In the following model, as in previous studies Fioriti2009CIIS (); Dorfler2010ACC (), we do not consider the effect of the length of the power line on the power loss. To estimate the power loss in a typical power line shown in Fig. 5, we estimate active power flow ( and ), reactive power flow ( and ), and delayed reactive power flow ( and ) as follows Kundur1993Elgerd1982Book ():

 Pij=ViVjsin(θi−θj)Zij2/f0Lij+Vi2−ViVjcos(θi−θj)Zij2/Rij, (10) Qij=−ViVjsin(θi−θj)Zij2/Rij+Vi2−ViVjcos(θi−θj)Zij2/f0Lij, (11) Qci=f0Cij2Vi2, (12)

where represents synchronized angular frequency of alternating voltage, and . , , and are obtained by exchanging and . Using these power flows, the active power loss due to resistance, , is calculated as , whereas the reactive power loss due to inductance, , is estimated as as follows:

 Pijloss=RijZij2(−2ViVjcos(θi−θj)+Vi2+Vj2), (13) Qijloss=Zij2f0Lij(−2ViVjcos(θi−θj)+Vi2+Vj2)+f0Cij2(Vi2+Vj2). (14)

The total power loss is estimated as a combination of the active power loss and the reactive power loss Kundur1993Elgerd1982Book (). By using a second-order Taylor expansion, we regard the total power loss, , as , where and are constants (). Therefore, we define the synchronization cost, , for a power line between power plants and as

 Sij=(θi−θj)2. (15)

## Appendix B Power Grid as Kuramoto Model

In this section, we explain that, under several assumptions, we can approximate power grids by the first-order Kuramoto model of nonuniform oscillators.

As previous studies Buzna2009PRE (); Fioriti2009CIIS (); Dorfler2010ACC (), we model a power grid as follows: The structure of the power grid with power plants is represented as an unweighted and undirected adjacency matrix , where a node represents a power plant and an edge a power line. is when power plants and have a power line between them, and is when they do not. According to the previous studies Buzna2009PRE (); Fioriti2009CIIS (); Dorfler2010ACC (), the phase of the output voltage of the power plant , , is described as

 ˙θi=fiDi−∑jWijDiAijsin(θi−θj), (16)

where denotes a damping constant, is an amount of power transfer between power plants and , and represents the natural frequency of the output voltage from the power plant . To reduce computational cost for the following rewiring-based optimization, we assume that for any power line. Because is specific to power plant , we replace the value with . Consequently, the voltage phase of the power plants can be expressed in the Kuramoto model as

 ˙θi=ωi−ϵ∑jAijsin(θi−θj). (17)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters