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

Ric h-club net w ork top ology to minimize sync hronization cost due to phase d ifference among frequency-sync hron ized oscil lators T ak amitsu W atanab e ∗ 1 1 Dep artment of Physiolo gy, Scho ol of Me dicine, The Universi ty of T okyo, 7-3-1 Hongo, Bunkyo-ku, T okyo 113-8656, Jap an Abstract F unctions of some net w ork s , such a s p ow er grids a n d large-scale brain net w orks, rely on not only fr e- quency syn c hronization, but also phase synchronizati on. Neve rtheless, eve n after the oscillators reac h to frequency-sync hr onized stat us, phase difference among oscillators often sho ws non -zero constan t v alues. Such phase difference p oten tially results in inefficien t transfer of p ow er or inform ation among oscillato r s, and a vo id pr op er and efficient functioning of the net work. In the p resen t study , we newly define synchronizati on cost b y the phase difference among the frequen cy-sync hr on ized oscillators, and in vestig ate the optimal net wo rk structure with the minimum sync h ronization cost through rewiring- based optimization. By using the K uramoto mo del, we demonstrate that the cost is minimized in a n et work top olog y with r ic h-club organization, whic h compr ises the densely-connected cen ter no des and p eripheral no des conn ecting with the cente r mo dule. W e also sho w that the net work top ology is c haracterized by its bimo d al degree distribu tion, whic h is q u an tified by W olfson’s p olarizatio n index. F urthermore, w e p ro vid e analytical inte rpretation on w h y th e ric h -club net work top ology is related to the small amount of syn c hr on ization cost. ∗ tak aw atanab e-tk y @umin.ac.jp 1 I. INTR O DUCTION As pow er grids [1, 2] and net works of bursting neurons [3, 4], functions of some complex net w o rks of oscillators a r e based on not only sync hronizatio n o f frequencies of the oscillators, but also sync hronization of their phases. Ho we v er, in general, frequency sync hronization is more a c hiev able than phase sync hro nizatio n. Phase difference among frequency-sync hronized oscillators often falls in t o a non-zero constan t, and suc h non- zero phase difference would av oid prop er and efficie n t functioning of the complex net w or ks. In p o w er grids, alternating v oltage of the p o w er plan ts in the grids should b e sync hronized around certain sp ecific frequency (e.g., 50 Hz in the most parts of Europ e and 60 Hz in the north America) [1], and disruption of the frequency sync hronization causes a blac kout in a large area [1, 5]. In addition, the phases of the v oltages of the p ow er plan t s ar e a lso r equired t o b e sync hronized. As discussed in App endix A, the difference in v oltage phases among p ow er pla nts inevitably causes p ow er lo ss consumed as heat in p ow er lines [1, 2 ]. In this sense, the phase difference in p o wer g rids can b e regarded as sync hronization cost. Considering recen t increasing en vironmen tal a warene ss and soaring global demand of natural resources [6, 7], it is necessary to reduce the sync hronization cost due t o the phase difference among frequency-sync hronized p o w er pla n ts. Sync hronization in large- scale brain net w orks also requires less difference in phase of neu- ronal activit y among differen t brain regions. A series of prior exp erimen tal res ear ches ha ve sho wn that v arious imp orta n t functions in large-scale bra in net works a re based on not only frequency- but also phase- sync hronization[3 , 4 ]. A previous electrophy siolog ical study sho wed that spik e activity recorded from monk eys’ cortex exhibited phase sync hronization in v arious frequency bands while t he monk eys were conducting ta sks that required in tegratio n o f visual pro cessing and motor resp onses [8]. In another careful electroph ysiology study , Ro elfsema and his colleagues recorded lo cal field p ot entials (LFP) in the cerebral cortex in cats, a nd re- v ealed that phase sync hronization b et w een LFPs recorded in the visual and parietal cort ices w as increased only when t he cats fo cused their atten tion on stim uli [9]. Studies using electro- encephalogram to record h uman brain a ctivity fo und that increase o f phase sync hronization in v arious frequency w as asso ciated with learning and p erception of images [10, 1 1 ]. F urthermore, the frequency- and phase- sync hronization are considered to o ccur in large-scale brain netw orks [8, 12] with zero time-lag [4, 9, 1 3 ]. These researc hes suggest that phase sync hronization in 2 large-scale brain net works is related to crucial functions suc h as integration of m ultiple infor- mation [3], neural comm unication [14 ], a nd spik e-timing- dep enden t plasticity [4]. Actually , it is kno wn that some t yp es of the disruption of the sync hronization cause dysfunctions of mem- ory learning [15] or psyc hiatric disorders [16]. Considering these findings, it is to some exten t reasonable to h yp othesize that large-scale brain net w orks ha ve a sp ecific orga nization that min- imizes phase differences amo ng brain activity , and enables optimal phase sy nchronization in the en tir e net w o r ks. These previous literatures indicate the imp ortance to reduce phase difference among frequency-syn c hronized oscillators, whic h can b e regarded as a type of cost that is sp en t during sync hronization. H o wev er, little is examined ab out the optimal net work to p ology that reduces this type o f sync hronizatio n cost due to phase difference. Inde ed, a series of previous literatures ha ve in v estigated optimal net w ork structures by intro ducing a different t yp e of sync hronization cost, whic h is needed for building or maintain of the optimal net w o r k infrastructure. A study regarded par ameters ba sed o n coupling strength among oscillators as a cost, and demonstrated that homogeneous and unifo rm distribution of the coupling strength enhances the tendency of sync hronization [17]. Another study emplo ying the same definition of sync hro nizatio n cost suggests that more heterogeneous net w ork structures are required to sim ultaneously ac hiev e b oth t he maxim um sync hronizability and minim um sync hronization cost ba sed o n coupling strength [18]. Another study rev ealed the optimal distribution of coupling strength that in- creases sync hronizabilit y a mong Kuramoto oscillators [19]. Ho we v er, these prior research es ha ve not fo cused on sync hronization cost due to phase difference. Consequen tly , despite a line of prior researc hes on optimal conditions for sync hrony in net w o r ks [20], it remains unclear ab out optimal netw ork structures that minimize sync hronization cost due to phase difference among oscillators. In the presen t study , therefore, w e examine the optimal net work top olog y to minimize the phase-difference-based sync hronizatio n cost. W e ta ke the follo wing fiv e steps: (i) First, w e define the sync hronization cost, S ij , due to phase difference b etw een frequency- sync hronized phase oscillators i and j in the Kuramoto mo del. W e adopt the Kuramoto mo del b ecause t he mo del has b een used as appro ximation of v arious systems including p o w er grids [2, 21, 22] and large- scale brain net works [23, 2 4]. (ii) Second, by using the definition, w e nu merically calculate the mean of sync hronizatio n 3 cost, h S i , in the entire net work. (iii) Third, b y using a rewiring strategy [25], w e show that ”ric h- club” net w ork top ology [24, 26 –28], which consists of densely in ter-connected mo dules and p eripheral lo w- degree no des, is the optimal top ology with the minimum sync hronization cost. (iv) F orth, w e c har a cterize the ric h-club net work top ology by quan tif ying the bimo dality of the degree distribution of the netw ork. (v) Finally , we pro vide an analytical in terpretation on wh y the ric h- club netw ork organization is asso ciated with a small a moun t of the sync hronization cost. I I. METHO D A. Definition of sync hronization cost W e first define sync hronization cost due to phase difference in an unw eigh ted and undirected net w o rk whic h is describ ed b y an adjacency matrix A and consists of N phase oscillators. A ij is 1 when oscillators i and j are connected, a nd A ij is 0 when they are not. According to the Kuramoto mo del, the phase of the oscillator i , θ i , is describ ed as ˙ θ i = ω i − ǫ X j A ij sin( θ i − θ j ) , (1) where ω i is the natural frequency o f the oscillator i , and ǫ is a coupling strength. T o reduce computational cost for the following rewiring-ba sed optimization, w e assume t ha t ǫ is constant for a ny com binatio n of oscillators. In t his Kura mo t o mo del, the sync hro nization cost b et w een oscillators i and j , S ij is defined based on the phase difference b et we en the connected oscillators as follows : S ij = ( θ i − θ j ) 2 . (2) Note that the sync hronizatio n cost is only defined a fter the net w ork of the oscillators reach es to a state o f frequency sync hro nizatio n. As describ ed in App endix A, in p ow er grids, S ij can b e regarded as a n index tha t is prop ortio nal to p ow er loss due to difference in v oltag e phase b et wee n p ow er plants i and j . 4 B. Estimation of the mean sync hronization cost Based on t he definition of S ij , w e nume r ically es t imate S ij for eac h edge in the fo llowing four steps fo r a giv en net work: (i) W e set normally-distributed { ω i } for eac h no de. It is b ecause t ha t previous studies on real net w o r ks suc h as p ow er grids and brain net works hav e assumed that the natura l frequencies of b elonging oscillators a re symmetrically fluctuating a round the av eraged frequency [4, 2 2, 29]. (ii) Based on the Kuramoto mo del describ ed in Eq. (1), w e num erically estimate f r equency - sync hronized status, where ˙ θ i b ecomes a common constan t v alue, Ω, for any i . (iii) In the frequency-sync hronized status, eac h o scillator ha s a differen t sp ecific phase, θ i . Based on the set of { θ i } , w e then ev aluate S ij for eac h edge. (iv) As describ ed in (i), the set o f the natura l frequency , { ω i } , is fluctuating ov er time. Th us, the S ij is also fluctuating o ve r time. Therefore, w e rep eat the pro cedure (i)-(iii) 200 times with different sets of { ω i } , and obtain 2 0 0 differen t sets o f { S ij } . Then, w e a verage the { S ij } o ver time, obtaining h S ij i for each edge. F inally , w e a v erag e the h S ij i across edges, and obtain h S i for t he en tire net work. C. Rewiring-based optimization T o searc h for the optimal net w ork top ology with the least h S i , we adopt the rewiring metho d that previous studies used to explore the netw ork top olog y with the largest sync hronizabilit y [25]. W e apply the followin g rewiring-based optimization pro cedure to a giv en connected net- w ork with N no des a nd mean degree of h k i : A t eac h step, the n umber of rewired edges is ran- domly determined based on an exp o nential distribution. The set of edges to b e rewired is also randomly chosen in a give n net work. After the rewiring, w e estimate f requency-sy nchronized status and obtain h S i updated . The attempted r ewiring is rejected if the up da ted net work is disconnected. Otherwise, t he rewiring is accepted if ∆ h S i = h S i updated − h S i initial < 0, or with probabilit y p = min(1 , [1 − (1 − q )∆ h S i /T ] 1 / (1 − q ) ) where T is a temp erature-lik e parameter and q = − 3 [25]. The initial rewiring is conducted at T = ∞ , and, af ter the first N rewiring, T is set as (1 − q )(∆ h S i ) max where (∆ h S i ) max is the larg est ∆ h S i in the first N rewiring trials. After 5 that, T is decreased 10% in ev ery 10 rewiring trials. This estimation pro cess iterated until there is no c hange in more than 50 successiv e rewiring steps. W e apply this rewiring-based optimiza- tion to three differen t initial net works: an Erd˝ os-R ´ en yi (ER) random mo del, W atts-Strogatz (WS) mo del [30], and a Barab´ asi-Alb ert (BA) mo del [31] with N = 50 and h k i = 4 [25]. In all the cases, the coupling strength, ǫ , is set to b e 0 . 3. Eac h set of the natural frequencies of Kuramoto oscillators, { ω i } , is randomly chose n from the normal distribution with an av erage of 100 π and a standar d deviation of 1. During the optimization, w e tr a ce the standard order parameter, r , and lo cal sync horoniz- abilit y , r lo cal [32], defined a s follo ws: r e iψ = 1 N X j e iθ j , (3) r lo cal = 1 2 N l X i X j ∈ Γ i      lim ∆ t →∞ 1 ∆ t Z t r +∆ t t r e i [ θ i ( t ) − θ j ( t )] dt      , (4) where N l is the tota l nu m b er of edges, Γ i is the set of neigh b ors of no de i , a nd ψ is a globa l phase. F urthermore, after t he optimizatio n is completed, w e compare the optimized net w orks from the differen t initial net w orks b y estimating the followin g basic top olo gical prop erties: mean of shortest path length [33], mean of clustering co efficien t [30], mean of b etw eennes s cen trality [34], and degree correlation [34]. W e conducted ten optimizations of ten differen t netw orks fo r eac h t yp e of the initial net w orks, and av eraged t hese basic prop erties. • The shortest path length, ℓ ij , is defined as the shortest distance b etw een tw o no des i and j [33 ]. The a ve r a ged shortest path length, h ℓ i , is defined as the a verage v alue o f ℓ ij o ver all the p o ssible pairs o f no des in the net w ork. • The clustering co efficien t for no de i , C i , measures the lo cal group cohesiv eness [30], whic h is defined as the ratio of the n umber of links b etw een the neighbors of i and the maxim um n umber of suc h links. W e define the mean clustering co efficien t, h C i , as a n a verage v alue o ver all the no des. • The b et w eenness cen t ralit y for no de i , b i , is defined as t he n umber of shortest paths b e- t we en pairs of no des that pass through a giv en no de [34 ]. W e define the mean b etw eennes s cen trality , h b i , a s a n av erage v alue o v er all the no des. • The degree corr elat io n fo r a netw ork is defined a s the P earson a ssortativit y co efficien t of the degrees, r assortativ e [34]. The co efficien t enables us to quantify the preference fo r 6 high-degree no des to attac h to o t her high-degree no des. Net works with this preference sho w large r assortativ e . D. Estimation of ric h-club co efficien t W e estimate reic h-club co efficien t, Φ( k ), for b ot h initial net works and optimized net works . According to the previous studies [24, 26 – 28], the co efficien t for eac h degree k is calculated as Φ( k ) = 2 E >k N >k ( N >k − 1) , (5) where E >k represen ts the n umber of edges among N >k no des that ha v e more than k degrees. As in previous studies [24, 27, 28], w e calculate normalized ric h-club co efficien ts, Φ norm ( k ) through dividing the raw v alue, Φ( k ) , by the mean of ric h- club co efficien t s of 100 random net works (ER mo dels), Φ random ( k ), as fo llo ws, Φ norm ( k ) = Φ( k ) Φ random ( k ) . (6) I I I . RESUL TS A. Rewiring-based optimization Fig. 1 sho ws represen tativ e results of the r ewiring-based optimizations. h S i w as decreased from approxim ately 6 . 5 × 10 − 2 to 4 . 5 × 10 − 2 ev en when the initial net work structure was differen t. In all the initial netw orks, h S i reac hed to a stable status after approximate 400 steps of rewiring. Strikingly speaking, w e cannot guarante e that the optimal net w ork w a s found, but this result suggests that a reasonably robust a ppro ximation of the optimal top ology w a s obtained in this metho d. During the optimization, the standard order parameter, r , w ere fluctuating just b elow 1 during the optimization (a small panel in Fig. 1 A). The lo cal sync honizabilit y , r lo cal , sho w ed the similar fluctuation b elow 1. Considering the previous studies o n these parameters [32], these b eha viors of r and r lo cal are considered to b e related to the amoun t of the coupling strength, ǫ . The previous studies [32] ha ve demonstrated that, when the coupling strength is more than 0.2, b o t h of r and r lo cal reac h a plateau that is near to 1 regardless of net w ork top ology . In the presen t study , the coupling strength, ǫ , w as se t at 0.3, b ecause global sync hronizatio n is 7 necessary for the estimation of h S i . This relativ ely large coupling strength could result in the saturation of the global and lo cal order parameters, r and r lo cal during the optimization pro cess. Fig. 1 B sho ws that the optimized netw orks for the three differen t initial netw orks commonly exhibit a c ha r a cteristic top olog y , which has a densely inte rconnected core no des a nd p eripheral no des dangling the core mo dule. The heterogeneous net w ork features w ere also observ ed in the basic net w ork prop erties in the o pt imized net w orks (T ab. I). Compared with the initial net w o rks, the o ptimized netw orks tended to show larger a v eraged v alues of the shortest path length, h ℓ i , b et wee nness cen trality , h b i , and degree correlation, r assortativ e . The a veraged v alues of the clustering co efficien ts, h C i , we r e smaller in the optimized net works . These results suggest that, through the optimization pro cess, the netw ork seems to enlarge its heterogeneit y . B. Ric h-club organization This heterog eneous net work top ology has b een rep orted as ”rich-club” organization in a series of previous theoretical and exp erimen tal studies [24, 2 6–28]. The prio r literatures hav e c haracterized the organizatio n b y estimating no rmalized r ic h-club co efficien ts, Φ norm ( k ) de- scrib ed in Eq. 6. If t he net work has ric h- club organization, Φ norm ( k ) should b e more than 1, and increase monotonically as k increases. Fig. 2 sho ws the comparison in Φ norm ( k ) b etw een initial net works and optimized net w orks. T o clarify the difference in the ric h-club co efficien t s, we adopted as the initial net w orks larger net w o rks than sho wn in Fig . 1 B ( i.e., N = 100 , h k i = 4). Before the optimizatio n, Φ norm ( k ) w as not alw ays larg er than 1 ( e.g., ER and WS mo dels), and did not sho w monotonic increase along k , whic h is consisten t with a previous study [2 7]. In contrast, in t he optimized netw orks, the rich-club co efficien ts w ere larger than 1 in almost all the range of k , and monotonically increased as k increased. These phenomena were o bserv ed commonly among the three differen t optimized net w orks that w ere deriv ed from the three differen t initial net works. In addition to the app earance of the netw orks in Fig. 1 B, this estimation of Φ norm ( k ) supp o r ts the notio n that the net w o rks with rich-club organization has the minim um or a v ery small amount of sync hronization cost due to phase difference among frequency-sy nc hronized oscillators. 8 C. Bimo dal Degree Distribution As sho wn in Fig. 1 B, the ric h-club net w ork top ology cons ists o f high- degree no des cluster and lo w-degree p eripheral no des. Therefore, we h yp othesized that the top olog y can b e characterized b y a bimo dal degree distribution. T o test the hy p othesis, w e estimated W olfson’s p olarization index, ˆ P [35]. The W o lfson’s index for degree distribution is defined a s: ˆ P = 2 h k i m (2 ( h k i 2 − h k i 1 ) − G ) , (7) where h k i is the mean of the degree, k i , and m denotes the median of the degree. h k i 1 and h k i 2 are the mean v alues of { k i   k i < m } and of { k i   k i ≥ m } , resp ectiv ely . G represen ts Gini inequalit y index, whic h is defined as G = 1 2 h k i P N i =1 P N j =1   k i − k j   . This W olfson’s p olarization index sho ws the exten t o f the bimo dality of the distribution. If the distribution is completely the same as a uniform distribution, the ˆ P is 0. If the ha lf of p opulation has nothing and the other half shares ev erything, the ˆ P reac hes a maxim um, 0 . 2 5. In the presen t case, a larger ˆ P indicates that t he net w or k has a more bimo dal and bip olarized degree distribution. W e estimated ˆ P during the rewiring-based optimization. Because ˆ P can b e calculated more accurately for net works with more no des, w e used the BA mo del with N = 100 and h k i = 4 as a n initia l net work for the optimizatio n. As a result, in the course of the ab ov e-men tioned optimization, h S i decreased during the rewiring-based optimization (Fig. 3A). Mean while, as h S i decreases, ˆ P increases (Fig. 3B). Actually , the degree distribution changed from a p ow er- la w distribution (Fig . 3C) to a bimo dal distribution (Fig. 3D). This relation w as also observ ed for differen t initial net works (e.g., ER mo del). This correlation supp orts the hypothesis that ric h-club net w ork with small h S i can b e characterized b y its bimo dal degree distribution. D. Analytical Interpretation W e finally pro vide an a nalytical interpretation on wh y the rich-club net w ork has less h S i . Using mean-field approximation, the Eq. (1) in the frequency-sync hronized status can b e de- scrib ed as Ω = ω i − ǫk i sin( θ i − ψ ), where ψ is defined in the Eq. (3). Therefore, if | θ i − ψ | is small enough, θ i − ψ = 1 ǫk i ( ω i − Ω), and ( θ i − θ j ) 2 is describ ed a s ( θ i − θ j ) 2 = ( ω i − Ω) 2 ǫ 2 k i 2 + ( ω j − Ω) 2 ǫ 2 k j 2 − 2( ω i − Ω)( ω j − Ω) ǫ 2 k i k j , (8) 9 for a set of { ω i } . h S ij i is obtained as av eraged ( θ i − θ j ) 2 across an enough n um b er of sets of { ω i } . As in the a b o v e-describ ed n umerical estimation, w e assume that { ω i } is distributed according to the no r mal distribution with a mean of ω 0 and a standard deviation of σ , a nd that the sync hronized frequency is alw ay s Ω in ev ery set of { ω i } . Because we can also assume that ω 0 is nearly equal t o Ω, h ( ω i − Ω) 2 i is considered to b e equal to σ 2 , and h ( ω i − Ω)( ω j − Ω) i is considered to b e equal to zero. Consequen tly , w e obtain the a ppro ximation of h S ij i as follows: ˜ h S ij i = σ 2 ǫ 2 ( 1 k i 2 + 1 k j 2 ) . (9) This a ppro ximation w as v alidated through comparison of ˜ h S ij i with the r eal h S ij i , a s show n in Fig. 4 A ( σ = 1 , ǫ = 0 . 3 ) . The t w o parameters had a large negativ e v alue of P earson’s correlation co efficien t ( − 0 . 88 ). This expression of ˜ h S ij i giv es qualitativ e explanation on relationship b et w een ric h-club net- w ork top ology and a small v alue of h S i : T o ac hiev e a small amoun t of h S i , ˜ h S ij i should be small. When the no de i has a hig h degree, it mak es more con tr ibutio n to a smaller ˜ h S ij i to connect with a the no de j with a high degree. It is also case when the no de i has a small degree. As a result, for a small v alue o f ˜ h S ij i , high degree no des should connect with other high degree no des, and low degree no des should not hav e edges with other low degree no des, but with high degree no des. Consequen tly , high degree no des tend to b e g a thered and create a densely-connec ted core mo dule, and low degree no des tend to connect with high degree no des in the core mo dule. Ov erall, the optimized netw ork with a small amount o f h S i is lik ely to b e a ric h- club net work. Note that it is difficult to further extend this approximation. If this approximation of h S ij i is accurate enough, a simple calculatio n of Eq. 9 leads us to the prop o rtional relationship b et wee n h S i and 1 h k i  1 k  . Given h k i is a constan t v alue as in the presen t study , h S i should be prop ortional to  1 k  . Ho w ev er, as sho wn in F ig. 4 B, w e could not observ e a linear relationship. This inaccurate approximation of h S i ma y b e caused b y a ccumulation of the small difference b et wee n h S ij i and ˜ h S ij i . This result suggests that w e cannot extend this appro ximation to represen t h S i only by  1 k  . 10 IV. DISCUSSION The presen t study introduced sync hronization cost based on phase difference a mong frequency-syn c hronized oscillators. Using the rewiring-based optimization [25, 36], w e show ed that the sync hronization cost is minimized in a ric h-club net w o rk top ology . F urthermore, we demonstrated that the netw ork to p ology can b e ch aracterized by the bimo dality o f its degree distribution. Finally , we pro vided analytical explanation on wh y the ric h- club net work top olog y is asso ciated with a small a moun t of sync hronization cost. The concept of sync hronization cost is not a no v el idea of the presen t study . As describ ed in Sec. I, a line o f previous studies hav e inv estigated a differen t t yp e of sync hronization cost, whic h is based on coupling strength [17–19]. Whereas the presen t sync hronization cost due to phase difference can b e regarded as dynamic cost p er unit time, t he cost based o n coupling strength can b e considered as static cost that is related to building a nd main taining of net- w ork infrastructures. Inte restingly , the optimal net w o rk to p ology with the least cost dep ends on whic h of t he tw o t yp es of sync hronization cost w e adopt. The optimal net w orks for the sync hronization cost based on coupling strength often sho w more homogeneous prop erties [17] than those f or the other sync hronization cost. The homogeneity of netw orks is desired to en- hance sync hronizability [25, 37]. Therefore, it ma y b e necessary to inv estigate what netw ork structures balance these t wo types of sync hronization cost. The pr esen t sync hro nization cost in the presen t study can b e a no t her concept of loa d as- signed to edge in a complex net work. Previous studies used edge-b et w eenness a s edge load [38, 39], whic h is useful in v arious situatio ns from hum an interaction [3 9] to data tr a nsmission in computer net w o rks [38 ]. Ho we v er, b ecause the edge b et we enness do es not consider sync hron y in net w orks, its prop erties ha v e eviden t difference from those of the sync hro nizatio n cost. F o r example, as sho wn in Sec. I I I D, the sync hronization cost is low er b et wee n high degree no des, and hig her b etw een lo w degree no des. In con trast, the edge betw eennes s tends to b e higher in edges bridging high degree no des and b e lo wer in edges bridging lo w degree no des. These distinct prop erties suggest the p ossibilit y that the sync hro nization cost can b e another concept of edge load. The sync hronization cost in the presen t study , S ij , has a mathematical expression similar to that for lo cal sync hronizability , r lo cal [32]. Ho wev er, the t wo parameters fo cus on differen t phases of sync hronization in complex netw orks. The lo cal sync hronizability enables us to quan tify 11 the lo cal construction o f the sync hronization pattern. The refore, it is useful to inv estigate prop erties of netw orks that are not yet fully sync hronized. In con trast, the sync hronizatio n cost in the presen t study can b e only estimated in fully-sync hronized netw orks. Therefore, in the presen t study , w e used a relative ly la rge coupling strength, a nd achie ved full sync hronization throughout the optimizations. As a result, in the presen t study , the lo cal sync hronizability w a s alw ays saturated. Although the presen t study did not ado pt mo dels sp ecific to an y real net w orks, the findings ma y help understanding large-scale brain net works. Recen tly , a few o f recen t studies hav e re- p orted the existence of the ric h-club org a nization in the large-scale brain netw orks. A previous empirical study has demonstrated the existe nce of the rich-club org anization in the large-scale h uman brain net w orks [28]. Another study has inv estigated the anatomical connectivit y in the cerebral cortex of cats, a nd has sho w ed that ric h-club orga nization con trols the dynamic transition of sync hronization in the brain [24]. A recen t review has suggested that the orga- nization is a cost- effectiv e netw ork top olog y fo r the brain net works, whic h are required to b e adapted to v arious cognitive functions [40]. In addition to the con text of cost-effectiv eness, the ric h-club netw ork top ology is robust to ra ndo m attack [41]. This previous study ana lytically and n umerically demonstrated that net works robust to ra ndom a ttac k hav e similar structures observ ed in t he presen t study . The ro bust net w orks highly-in terconnected h ub mo dules and p eripheral no des ( le af no des) t hat hav e a single edge. This net work top ology has bimo dal degree distribution and show s ric h- club organization. Considering these prior literatures, it is suggested that the rich-club o rganization is b eneficial for the large-scale brain net works to efficien tly and robustly maintain its wide range of functions based on sync hronization. Ac kno w ledgmen ts The author ac kno wledges the support from the Japan So ciet y for the Promotion of Science (JSPS) Researc h F ellow ship for Y oung Scien tists (222 8 82). [1] O. I. Elgerd, Ele ctric ener gy systems th e ory (McGra w-Hill Companies, 1982) ; P . S. Kundu r, Power system stability and c ontr ol (McGra w-Hill C omp anies, 1993). [2] F. D¨ orfler and F. Bullo, American Con trol Conference 2010, 930-937 (2010). 12 [3] F. V arela, J.P . Lac h aux, E. Ro driguez, and J. Martinerie, Nat Rev Neurosci. 2(4) :229-39 (2001 ). [4] J. F ell, and N. Axmac her, Nat Rev Neurosci. 12(2) :105-18 (2011). [5] D.J. Hill, and G. C hen, IEEE Int. Symp osium on Circuits and S ystems 711-715 (2006). [6] H. Gharavi, and R. Gh afu rian, Pro ceedings of the IE EE 99 , 917-92 1 (2011); G.W. Arnold, ibid . 99 , 922-927 (2011). [7] Energy F u tu re Coalition, Rep ort of S mart Grid W orkin g Group , (2008). [8] S.L. Bressler, R. C op p ola, and R. Nak am ur a, Nature. 366(6451) :153 -6 (199 3). [9] P .R. Ro elfsema, A.K. Engel, P . K¨ onig, and W. Singer, Nature. 385(6612) :157-61 (1997) . [10] W.H. Miltner, C. Braun, M. Arnold, H. Witte, E. T aub. 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. T rau b , M.A. Whittington, I .M. Stanford, and J.G. Jefferys, Nature. 383(6601) :621-4 (1996 ). [13] T. W omelsdorf, J-M. Sc h offelen, R. Oostenv eld, W. S inger, R. Desimone, A.K. Engel, and P . F ries, S cience. 316(5831) :1609 -12 (2007). [14] P . F ries, T rends Cogn Sci (Regul Ed). 9(10) :474-80 (2005) . [15] M. Stopfer, S. Bhaga v an, B.H. S mith, and G. Lauren t, Nature. 390(6655) :70-4 (1997) . [16] G. T ononi, and G.M. Edelman, Brain Res Brain Res Rev. 31(2-3) :391-400 (2000). [17] A.E. Motter, C. Zhou, an d J. Ku rths, Phys. Rev. E 71 , 016116 (2005 ). [18] T. Nishik a wa, 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 . Kurth s, Y. Moreno, and C. Zh ou, 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. Ruzzan te, E. Castorini, E. Marc hei, and V. Rosato, Critical Information Infrastruc- ture Security 14-23 (2009) . [23] M. Breaksp ear, S. Heitmann, and A. Daffertshofer, F ront. Hum. Neurosci. 4 :190 (2010). [24] J. G´ o mez-Garde ˜ n es, G. Zamora-L´ opez, Y. Moreno, and A. Arenas. PLoS O NE. 5(8) :e1231 3 (2010 ). [25] L. Donetti, P .I. Hurtado, and M.A. Mu ˜ noz, Phys. Rev. Lett. 95 , 188701 (2005). [26] S. Zhou, and R. Mondr agon. IEEE Commun Lett. 8(3) :180-2 (2004). [27] V. Colizza, A. Flammini, M. Serrano, and A. V esp ignani, Nat Phys. 2(2) :110-5 (2006). [28] M.P . v an den Heuv el, and O. Sp orns, Journal of Neuroscience. 31(44) :15775 -86 (2011) . 13 [29] G. Filatrella, A.H. Nielsen, and N.F. Pedersen, Eur. Phys. J. B. 61 , 485-49 1 (2008). [30] D.J. W atts and S .H. Strogatz, Nature 393 , 440 (1998). [31] A.L. Barab´ asi and R. Alb ert, Science 286 , 5439 (1999). [32] J. G´ o mez-Garde ˜ nes, Y. Moreno, and A. Arenas, Ph ys. Rev. Lett. 98 034101 (2007), Phys. Rev. E. 75 066106 (2007 ). [33] A. Bund e, and S. Ha vlin, F r actals and Disor der e d System (Spr inger V e rlag, 1991). [34] M.E.J. Newman, Phys. Rev. E. 64 016131 (2001 ). [35] M. W olfson, Review of Income and W ealth 43 , 4 401 (1997). [36] TE Goro c h o wski, M d i Bernard o, CS Grierson. Phys Rev E. 81 056212 (2010). [37] T. Nishik a w a, A.E. Motte r, Y.C. Lai, and F.C. Hopp en steadt, Phys. Rev. Lett. 91 , 014101 (2003). [38] P . Holme, Ph ys. Rev. E 66 , 036119 (2002). [39] R. Guimer` a, L. Danon, A. Diaz-Guilera, F. Giralt, and A. Aren as, Phys. Rev. E 68 , 065103 (2003 ). [40] E. Bullmore, and O. S p orns , Nat Rev Neurosci. 336-349, (2012) . [41] G. Pa ul, S . Sreeniv asan, S. Ha vlin, and H.E. Stanley , Physica A. 370(2) :854-62 (2006). 14 6 5 100 iteration From an ER model A From a BA model From a WS model 0 200 300 400 7 WS Model ER Model BA Model -2 ×10 B iteration 0 r 0.94 0.97 200 400 FIG. 1: (Color online) A . Main pan el: Ch ange of the syn c hronization cost, h S i , dur in g rewiring-based optimization. De s pite differen t initial net works (ER, BA, and WS m o dels) with N = 50 , h k i = 4, th e sync h ronization cost con v erge to a similar amount of h S i . Su b panel: Ch ange of the standard order parameter, r , d uring the optimization. In con tr ast to h S i , the stand ard order parameter do es not show notable c h an ge, just flu ctuating b elo w 1. Th e lin e shows the c hange of r wh en the initial n et wo r k is the BA mo del. In cases of the other t w o initial net w orks , the similar fluctuations were observ ed. B . Net wo r k top ology optimized from differen t initial n et wo rks. O ptimized netw orks are similar to eac h other. They ha ve ric h -club n et wo rk top ology , w hic h consists of a densely-connected core mo du le and p eriph eral lo w degree no d es conn ecting with the core. The color in th e no d es represent the degree of the no des: dark er n o des ha ve more edges. 4 8 12 16 20 4 0 1 2 3 from a BA model from a ER model from a WS model 16 4 8 12 3 0 1 2 16 4 8 12 4 0 1 2 3 initial network optimized network FIG. 2: Difference in ric h-club co efficien ts b etw een the initial n etw orks (circles and dashed lines) and optimized n et wo rks (m ultiple marks and solid lines). While the normalized ric h -club coefficient s Φ norm ( k ) d o not show monotonic increase in the initia l n et wo r ks, those in th e optimize d netw orks monotonically in cr ease. These results suggest that the rewiring-based optimization c hanges the in itial net works to n etw orks with rich-cl ub organizat ion. T o clarify the difference b et w een b efore and after optimization, we adopted larger n et wo rks ( N = 100 , h k i = 4) th an in Fig. 1 ( N = 50 , h k i = 4). iteration # of Nodes # of Nodes A B C D 40 0 10 20 30 60 0 20 40 20 0 5 10 15 0 10 20 0.05 0.10 0.15 0.128 0.120 0.124 0 200 400 0.120 0.124 0.128 FIG. 3: As the s ync hr onization cost, h S i , decreases in the r ewir ing-based optimization (panel A ), the W olfson’s p olarization ind ex, ˆ P , increases (panel B ). Th is relation su ggests that the r ic h-club net work with a sm all amount of the syn c hr on ization cost can b e c haracterized by bimo dal degree distr ib ution, whic h is quan tified by W olfson’s p olarizatio n index. Ind eed, the degree d istr ibution c hanged from a p o wer-la w distribution (panel C ) to a bimod al distribution (panel D ). T o clarify the difference b et ween b efore and after optimization, we adopted larger net works ( N = 100 , h k i = 4) th an in Fig. 1 ( N = 50 , h k i = 4). 0.2 0.2 A B 0.6 0.6 1.0 1.0 1.4 1.4 1.8 1.8 -3 ×10 -1 ×10 -3 ×10 0.132 0.136 0.140 1.16 1.20 1.24 1.28 FIG. 4: A T h e analytical approximat ion of sync hronization cost, ˜ h S ij i , is predictiv e of the real sync h ronization cost, h S ij i , in the BA mo del with N = 200 , h N i = 10. B S imple calc ulations u sing ˜ h S ij i suggest a p ositive linear relati onship b et ween h S i and  1 k  . Ho we v er, ther e is not a strong correlation b et w een them, wh ic h s u ggests a limitation of the approximati on. T ABLE I: Basic top ological p rop erties of th e optimized net wo r ks. Despite d ifferen t initial net w orks, the optimized n et wo rks had similar net work top ological prop erties. The v alues for the initial n etw orks represent a v eraged v alues across ten estimations, whereas the v alues for the optimized net works sho w the mean ± s.d. across the ten estimations. h ℓ i h C i h b i r assortativ e initial optimized initial optimized initial optimized initial optimized F rom E R mo dels 1.9 3.5 ± 0.025 0.21 0.13 ± 0.010 44 122.2 ± 3.7 0.034 0.25 ± 0.12 F rom BA mo dels 1 .8 3.3 ± 0.021 0.33 0.15 ± 0.012 48 113 ± 2.8 -0.15 0.28 ± 0.021 F rom WS mo dels 2.3 3.1 ± 0.017 0.6 2 0.16 ± 0.010 63 105 ± 2.4 0.047 0.25 ± 0.011 App endix A: Definition of Sync hronization Cost In this section, w e explain why p ow er loss consumed in the electric line b et we en p ow er plants can b e represen ted by square of difference in phase of v oltage b et w een the tw o p ow er plan ts. In the following mo del, a s in previous studies [2, 22], we do not consider the effect of the length of the p o w er line on the p o we r loss. T o estimate the p ow er loss in a typic al p ow er line sho wn in Fig. 5 , w e estimate active p ow er flo w ( P ij and P j i ), r eactiv e p o w er flo w ( Q ij and Q j i ), and dela ye d reactiv e p o wer flo w ( Q ci and Q cj ) as follo ws [1]: P ij = V i V j sin( θ i − θ j ) Z ij 2 /f 0 L ij + V i 2 − V i V j cos( θ i − θ j ) Z ij 2 /R ij , (A1) Q ij = − V i V j sin( θ i − θ j ) Z ij 2 /R ij + V i 2 − V i V j cos( θ i − θ j ) Z ij 2 /f 0 L ij , (A2) Q ci = f 0 C ij 2 V i 2 , (A3) where f 0 represen ts sync hronized angular fr equency of alternat ing v oltage, and Z ij 2 = R ij 2 + ( f 0 L ij ) 2 . P j i , Q j i , and Q cj are obtained b y exc hanging i and j . Using these p o w er flows, the activ e p o w er loss due to resistance, P ij loss , is calculated as P ij + P j i , whereas the r eactiv e p o we r loss due to inductance, Q ij loss , is estimated as Q ij + Q j i + Q ci + Q cj as follows : P ij loss = R ij Z ij 2  − 2 V i V j cos( θ i − θ j ) + V i 2 + V j 2  , (A4) Q ij loss = Z ij 2 f 0 L ij  − 2 V i V j cos( θ i − θ j ) + V i 2 + V j 2  + f 0 C ij 2 ( V i 2 + V j 2 ) . (A5) The total p ow er loss is estimated as a com bination of the activ e p o w er loss and the reactive p o w er loss [1]. By using a second-order T aylor expansion, w e regard the tot a l p o w er loss, P ij loss + Q ij loss , as a 0 + a 1 ( θ i − θ j ) 2 , where a 0 and a 1 are constants ( a 1 > 0). Therefore, w e define the sync hronization cost, S ij , for a p o we r line b et wee n pow er plan ts i and j as S ij = ( θ i − θ j ) 2 . (A6) App endix B: Po w er Grid as Kuramoto Mo del In this section, we explain tha t , under sev eral assumptions, w e can approximate p ow er grids b y t he first-order Kuramoto mo del o f non uniform oscillators. As previous studies [2, 21, 22 ], we mo del a p ow er grid as f ollo ws: The structure of the p ow er grid with N pow er pla n ts is represen ted as a n un we igh ted and undirected adjacency matrix A , where a no de represen ts a p ow er plant and an edge a p o w er line. A ij is 1 when p ow er plan ts i and j hav e a p o w er line b et w een them, and A ij is 0 when they do not. According to the previous studies [2, 21 , 22], the phase of the output voltage of the p ow er plant i , θ i , is described as ˙ θ i = f i D i − X j W ij D i A ij sin( θ i − θ j ) , (B1) where D i denotes a damping constant, W ij is an amount of p ow er transfer b et we en p o w er plan ts i and j , and f i represen ts the natural fr equency of t he output v o lt age from the p ow er plan t i . T o reduce computat io nal cost for the fo llowing rewiring- ba sed optimization, w e assume that W ij /D i = ǫ for any p ow er line. Because f i /D i is sp ecific to p o w er plant i , we replace the v alue with ω i . Consequen tly , the v oltage phase o f the p o wer plan t s can b e expresse d in the Kuramot o mo del as ˙ θ i = ω i − ǫ X j A ij sin( θ i − θ j ) . (B2) A B Power Plant i Power Plant i Power Plant j Power Plant j FIG. 5: P anel A s ho ws a typica l p o we r line b et ween p o w er p lan ts i and j . V i cos θ i and V j cos θ j indicate the volta ge of the output fr om the p o w er plan ts. R ij , L ij , and C ij indicate resistance, inductance, and condu ctance b et w een the p ow er p lan ts. As sh o w in p anel B , S is defin ed in every p o w er line based on the p hase difference of the volt ages b et wee n the connecting p o we r plants.