Richclub network topology to minimize synchronization cost due to phase difference among frequencysynchronized oscillators
Abstract
Functions of some networks, such as power grids and largescale brain networks, rely on not only frequency synchronization, but also phase synchronization. Nevertheless, even after the oscillators reach to frequencysynchronized status, phase difference among oscillators often shows nonzero 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 frequencysynchronized oscillators, and investigate the optimal network structure with the minimum synchronization cost through rewiringbased optimization. By using the Kuramoto model, we demonstrate that the cost is minimized in a network topology with richclub organization, which comprises the denselyconnected 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 richclub network topology is related to the small amount of synchronization cost.
pacs:
89.75.k, 05.45.Xt, 88.80.HI 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 frequencysynchronized oscillators often falls into a nonzero constant, and such nonzero 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 frequencysynchronized power plants.
Synchronization in largescale 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 largescale 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 electroencephalogram 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 largescale brain networks Bressler1993Nature (); Traub1996Nature () with zero timelag Roelfsema1997Nature (); Womelsdorf2007Science (); Fell2011NatRevNeuro (). These researches suggest that phase synchronization in largescale brain networks is related to crucial functions such as integration of multiple information Varela2001NatRevNeuro (), neural communication Fries2005TrendsCogSci (), and spiketimingdependent 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 largescale 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 frequencysynchronized 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 phasedifferencebased synchronization cost. We take the following five steps:

First, we define the synchronization cost, , due to phase difference between frequencysynchronized 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 largescale brain networks Breakspear2010FrontHN (); GomezGardenes2010PLosOne ().

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

Third, by using a rewiring strategy Donetti05PRL (), we show that ”richclub” network topology Zhou2004IEEEComm (); Colizza2006NatPhys (); GomezGardenes2010PLosOne (); Heuvel2011JNS (), which consists of densely interconnected modules and peripheral lowdegree nodes, is the optimal topology with the minimum synchronization cost.

Forth, we characterize the richclub network topology by quantifying the bimodality of the degree distribution of the network.

Finally, we provide an analytical interpretation on why the richclub 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
(1) 
where is the natural frequency of the oscillator , and is a coupling strength. To reduce computational cost for the following rewiringbased 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:
(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:

We set normallydistributed 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 ().

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

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

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 Rewiringbased 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 rewiringbased 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 frequencysynchronized 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 temperaturelike 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 rewiringbased optimization to three different initial networks: an ErdősRényi (ER) random model, WattsStrogatz (WS) model Watts98Nature (), and a BarabásiAlbert (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, GomezGardenes2007PRLPRE (), defined as follows:
(3)  
(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 highdegree nodes to attach to other highdegree nodes. Networks with this preference show large .
ii.4 Estimation of richclub coefficient
We estimate reichclub 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
(5) 
where represents the number of edges among nodes that have more than degrees. As in previous studies Colizza2006NatPhys (); GomezGardenes2010PLosOne (); Heuvel2011JNS (), we calculate normalized richclub coefficients, through dividing the raw value, , by the mean of richclub coefficients of 100 random networks (ER models), , as follows,
(6) 
Iii Results
iii.1 Rewiringbased optimization
Fig. 1 shows representative results of the rewiringbased 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 GomezGardenes2007PRLPRE (), these behaviors of and are considered to be related to the amount of the coupling strength, . The previous studies GomezGardenes2007PRLPRE () 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 Richclub organization
This heterogeneous network topology has been reported as ”richclub” organization in a series of previous theoretical and experimental studies Zhou2004IEEEComm (); Colizza2006NatPhys (); GomezGardenes2010PLosOne (); Heuvel2011JNS (). The prior literatures have characterized the organization by estimating normalized richclub coefficients, described in Eq. 6. If the network has richclub 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 richclub 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 richclub 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 richclub organization has the minimum or a very small amount of synchronization cost due to phase difference among frequencysynchronized oscillators.
iii.3 Bimodal Degree Distribution
As shown in Fig. 1 B, the richclub network topology consists of highdegree nodes cluster and lowdegree 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 rewiringbased 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 abovementioned optimization, decreased during the rewiringbased optimization (Fig. 3A). Meanwhile, as decreases, increases (Fig. 3B). Actually, the degree distribution changed from a powerlaw 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 richclub network with small can be characterized by its bimodal degree distribution.
iii.4 Analytical Interpretation
We finally provide an analytical interpretation on why the richclub network has less . Using meanfield approximation, the Eq. (1) in the frequencysynchronized status can be described as , where is defined in the Eq. (3). Therefore, if is small enough, , and is described as
(8) 
for a set of . is obtained as averaged across an enough number of sets of . As in the abovedescribed 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:
(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 richclub 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 denselyconnected 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 richclub 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 frequencysynchronized oscillators. Using the rewiringbased optimization Donetti05PRL (); Gorochowski2010PRE (), we showed that the synchronization cost is minimized in a richclub 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 richclub 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 edgebetweenness 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, GomezGardenes2007PRLPRE (). 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 fullysynchronized 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 largescale brain networks. Recently, a few of recent studies have reported the existence of the richclub organization in the largescale brain networks. A previous empirical study has demonstrated the existence of the richclub organization in the largescale human brain networks Heuvel2011JNS (). Another study has investigated the anatomical connectivity in the cerebral cortex of cats, and has showed that richclub organization controls the dynamic transition of synchronization in the brain GomezGardenes2010PLosOne (). A recent review has suggested that the organization is a costeffective network topology for the brain networks, which are required to be adapted to various cognitive functions Bullmore2012NatRevNeuro (). In addition to the context of costeffectiveness, the richclub 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 highlyinterconnected hub modules and peripheral nodes (leaf nodes) that have a single edge. This network topology has bimodal degree distribution and shows richclub organization. Considering these prior literatures, it is suggested that the richclub organization is beneficial for the largescale 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 (McGrawHill Companies, 1982); P. S. Kundur, Power system stability and control (McGrawHill Companies, 1993).
 (2) F. Dörfler and F. Bullo, American Control Conference 2010, 930937 (2010).
 (3) F. Varela, J.P. Lachaux, E. Rodriguez, and J. Martinerie, Nat Rev Neurosci. 2(4):22939 (2001).
 (4) J. Fell, and N. Axmacher, Nat Rev Neurosci. 12(2):10518 (2011).
 (5) D.J. Hill, and G. Chen, IEEE Int. Symposium on Circuits and Systems 711715 (2006).
 (6) H. Gharavi, and R. Ghafurian, Proceedings of the IEEE 99, 917921 (2011); G.W. Arnold, ibid. 99, 922927 (2011).
 (7) Energy Future Coalition, Report of Smart Grid Working Group, (2008).
 (8) S.L. Bressler, R. Coppola, and R. Nakamura, Nature. 366(6451):1536 (1993).
 (9) P.R. Roelfsema, A.K. Engel, P. König, and W. Singer, Nature. 385(6612):15761 (1997).
 (10) W.H. Miltner, C. Braun, M. Arnold, H. Witte, E. Taub. Nature. 397(6718):4346 (1999).
 (11) J.F. Hipp, A.K. Engel, and M. Siegel, Neuron. 69(2):38796, (2011).
 (12) R.D. Traub, M.A. Whittington, I.M. Stanford, and J.G. Jefferys, Nature. 383(6601):6214 (1996).
 (13) T. Womelsdorf, JM. Schoffelen, R. Oostenveld, W. Singer, R. Desimone, A.K. Engel, and P. Fries, Science. 316(5831):160912 (2007).
 (14) P. Fries, Trends Cogn Sci (Regul Ed). 9(10):47480 (2005).
 (15) M. Stopfer, S. Bhagavan, B.H. Smith, and G. Laurent, Nature. 390(6655):704 (1997).
 (16) G. Tononi, and G.M. Edelman, Brain Res Brain Res Rev. 31(23):391400 (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, 7789 (2006).
 (19) M. Brede, Eur. Phys. J. B. 62, 8794 (2008).
 (20) A. Arenas, A. DiazGuilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469, 93153 (2008).
 (21) L. Buzna, S. Lozano, and A. DiazGuilera, Phys. Rev. E 80, 066120 (2009).
 (22) V. Fioriti, S. Ruzzante, E. Castorini, E. Marchei, and V. Rosato, Critical Information Infrastructure Security 1423 (2009).
 (23) M. Breakspear, S. Heitmann, and A. Daffertshofer, Front. Hum. Neurosci. 4:190 (2010).
 (24) J. GómezGardeñes, G. ZamoraLó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):1802 (2004).
 (27) V. Colizza, A. Flammini, M. Serrano, and A. Vespignani, Nat Phys. 2(2):1105 (2006).
 (28) M.P. van den Heuvel, and O. Sporns, Journal of Neuroscience. 31(44):1577586 (2011).
 (29) G. Filatrella, A.H. Nielsen, and N.F. Pedersen, Eur. Phys. J. B. 61, 485491 (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ómezGardeñ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. DiazGuilera, F. Giralt, and A. Arenas, Phys. Rev. E 68, 065103 (2003).
 (40) E. Bullmore, and O. Sporns, Nat Rev Neurosci. 336349, (2012).
 (41) G. Paul, S. Sreenivasan, S. Havlin, and H.E. Stanley, Physica A. 370(2):85462 (2006).


initial  optimized  initial  optimized  initial  optimized  initial  optimized  
From ER models  1.9  3.5 0.025  0.21  0.13 0.010  44  122.2 3.7  0.034  0.25 0.12 
From BA models  1.8  3.3 0.021  0.33  0.15 0.012  48  113 2.8  0.15  0.28 0.021 
From WS models  2.3  3.1 0.017  0.62  0.16 0.010  63  105 2.4  0.047  0.25 0.011 

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 ():
(10)  
(11)  
(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:
(13)  
(14) 
The total power loss is estimated as a combination of the active power loss and the reactive power loss Kundur1993Elgerd1982Book (). By using a secondorder 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
(15) 
Appendix B Power Grid as Kuramoto Model
In this section, we explain that, under several assumptions, we can approximate power grids by the firstorder 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
(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 rewiringbased 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
(17) 