What adaptive neuronal networks teach us about power grids
WWhat adaptive neuronal networks teach us about power grids
Rico Berner , , ∗ Serhiy Yanchuk , and Eckehard Sch¨oll Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, Hardenbergstr. 36, 10623 Berlin, Germany and Institut f¨ur Mathematik, Technische Universit¨at Berlin, Straße des 17. Juni 136, 10623 Berlin, Germany (Dated: June 12, 2020)Power grid networks, as well as neuronal networks with synaptic plasticity, describe real-worldsystems of tremendous importance for our daily life. The investigation of these seemingly unrelatedtypes of dynamical networks has attracted increasing attention over the last decade. In this Letter,we provide insight into the fundamental relation between these two types of networks. For this,we consider well-established models based on phase oscillators and show their intimate relation. Inparticular, we prove that phase oscillator models with inertia can be viewed as a particular classof adaptive networks. This relation holds even for more general classes of power grid models thatinclude voltage dynamics. As an immediate consequence of this relation, we find a novel type ofmulticluster state for phase oscillators with inertia. Moreover, the phenomenon of cascading linefailure in power grids is translated into an adaptive neuronal network.
Complex networks describe various processes in natureand technology, ranging from physics and neuroscienceto engineering and socioeconomic systems. Of particularinterest are adaptive networks, where the connectivitychanges in response to the internal dynamics. Such net-works can model, for instance, synaptic neuronal plas-ticity [1–4] or, generally, learning, memory, and devel-opment in neural circuits. Adaptive networks have beenreported for chemical [5], epidemic [6], biological, and so-cial systems [7]. A paradigmatic example of adaptivelycoupled phase oscillators has recently attracted much at-tention [8–17] and it appears to be useful for predictingand describing phenomena in more realistic and detailedmodels [18–21].A different class of network models describing powersystems as well as micro and macro power grids has beenanalyzed intensively [22–26]. It was shown that simplelow-dimensional models capture certain aspects of theshort-time dynamics of power grids very well [27–30]. Inparticular, the model of phase oscillators with inertia,also known as swing equation, has been widely used inworks on synchronization of complex networks and as aparadigm for the dynamics of modern power grids [31–47]. The phenomenon of converse symmetry breaking,predicted by this model, was demonstrated experimen-tally [48].Over the last years, studies on both types of models,phase oscillators with inertia and adaptively coupled os-cillators, revealed a plethora of common dynamical sce-narios including solitary states [41, 44, 45, 49], frequencyclusters [14, 15, 50, 51], chimera states [10, 12, 52], hys-teretic behavior and non-smooth synchronization transi-tions [9, 53–55]. Moreover, hybrid systems with phase dy-namics combining inertia with adaptive coupling weightshave been investigated, for instance, to account for achanging network topology due to line failures [56] orto include voltage dynamics [57].Despite the apparent qualitative similarities of the twotypes of models, so far nothing is known about their quantitative relationship. In particular, the followingquestion arises: is there a relation between phase oscilla-tor models with inertia and models of adaptively coupledphase oscillators?In this Letter, we show that dynamical power gridmodels have deep relations with adaptive networks. Inparticular, phase oscillators with inertia are a special sub-class of phase oscillators with adaptive couplings. Forthis, we introduce the so-called pseudo coupling matrix.In order to emphasize the strong implications of our find-ings, we provide two examples that illustrate our unifyingtheoretical approach. Such an approach permits to findmulticlusters which have not been known for oscillatorswith inertia. We also propose that the line failure ef-fect known for power grids might have its counterpart inadaptive neural networks with the presence of short termsynaptic depression being a short time equivalent to thefailure of a power line.We begin by introducing the two classes of models.The first class describes N adaptively coupled phase os-cillators and reads [10, 14]˙ φ i = ω i + N (cid:88) j =1 a ij κ ij f ( φ i − φ j ) , (1)˙ κ ij = − (cid:15) ( κ ij + g ( φ i − φ j )) , (2)where φ i ∈ [0 , π ) represents the phase of the i th oscil-lator ( i = 1 , . . . , N ), ω i is its natural frequency, and κ ij is the coupling weight of the connection from node j to i . Further, f and g are 2 π -periodic functions where f isthe coupling function and g is the adaptation rule, and (cid:15) (cid:28) a ij ∈ { , } of the adjacency matrix A . Note that thephase space of Eqs. (1)–(2) is ( N + N )-dimensional.The second class of models is given by N coupled phase a r X i v : . [ n li n . AO ] J un oscillators with inertia [41, 43] M ¨ φ i + γ ˙ φ i = P i + N (cid:88) j =1 a ij h ( φ i − φ j ) , (3)where M is the inertia coefficient, γ is the damping con-stant, P i is the power of the i th oscillator (related to thenatural frequency ω i = P i /γ ), h is the coupling function,and a ij is the adjacency matrix as defined in Eq. (1). Wenote that the phase space of (3) is 2 N -dimensional, i.e.,of lower dimension than that of Eqs. (1)–(2).In the following, we show that the class of phase oscil-lator models with inertia is a natural subclass of systemswith adaptive coupling weights where the weights denotethe power flows between the corresponding nodes. Wefirst write Eq. (3) in the form˙ φ i = ω i + ψ i , (4)˙ ψ i = − γM ψ i − γ N (cid:88) j =1 a ij h ( φ i − φ j ) . (5)where ψ i is the deviation of the instantaneous phase ve-locity from the natural frequency ω i . We observe thatthis is a system of N phase oscillators (4) augmented bythe adaptation (5) of the frequency deviation ψ i . Simi-lar systems with such a direct frequency adaptation havebeen studied in [58–61]. Note that the coupling betweenthe phase oscillators is realized in the frequency adapta-tion which is different from the classical Kuramoto sys-tem [62]. As we know from the theory of adaptivelycoupled phase oscillators [10, 14], an adaptation of thefrequency can also be achieved indirectly by a properadaptation of the coupling matrix.In order to introduce coupling weights into system (4)–(5), we express the frequency deviation ψ i as the sum ψ i = (cid:80) Nj =1 a ij χ ij of the dynamical power flows χ ij fromthe nodes j that are coupled with node i . The power flowsare governed by the equation ˙ χ ij = − (cid:15) ( χ ij + g ( φ i − φ j )),where g ( φ i − φ j ) ≡ − h ( φ i − φ j ) /γ are their stationaryvalues [43] and (cid:15) = γ/M . It is straightforward to checkthat ψ i , defined in such a way, satisfies the dynamicalequation (5).As a result, we have shown that system (4)–(5), andtherefore the swing equation (3) can be equivalently writ-ten as the following system of adaptively coupled phaseoscillators ˙ φ i = ω i + N (cid:88) j =1 a ij χ ij , (6)˙ χ ij = − (cid:15) ( χ ij + g ( φ i − φ j )) . (7)The obtained system corresponds to (1)–(2) with cou-pling weights χ ij and coupling function f ( φ i − φ j ) ≡ χ . Note that the base network topology a ij of the phase os-cillator system with inertia Eq. (3) is untouched by thetransformation. For a more rigorous and detailed expla-nation of the equivalence of the two systems, we refer tothe Supplemental Material [76].The obtained result suggests that the power grid modelis a specific realization of adaptive neuronal networks. In-deed, we proceed one step further and show that morecomplex models for synchronous machines like the swingequation with voltage dynamics [44, 57] can be repre-sented as adaptive network as well. The model reads M ¨ φ i + γ ˙ φ i = P i + N (cid:88) j =1 E i E j a ij h ( φ i − φ j ) (8) m i ˙ E i = − E i + E f,i + N (cid:88) j =1 a ij E j v ( φ i − φ j ) , (9)where the additional dynamical variable E i is the voltageamplitude. The functions h and v are 2 π -periodic, and m i and E f,i are machine parameters [44, 57]. All othervariables and parameters are defined as in (3). Re-writingthese equations as an adaptive network yields the samesystem (6)–(7) supplemented by Eq. (9) where g ( φ ) ≡− E i E j h ( φ ) /γ and (cid:15) = γ/M . Here it turns out that thevoltage dynamics can be interpreted as meta-plasticity ofthe coupling weights which is known from neuroscienceand technology [63–65]. For further details we refer tothe Supplemental Material [76].In the following, we provide two examples where theunification is used to find novel dynamical scenarios inphase oscillator models with inertia as well as in adap-tive networks by the transfer of known results from oneparadigm to the other. Mixed frequency cluster states in phase oscillator mod-els with inertia.
In the first example, we provide new in-sights into the emergence of multifrequency-cluster statesand report a novel type of multicluster state for the phaseoscillator model with inertia. In a multicluster state alloscillators split into M groups (called clusters) each ofwhich is characterized by a common cluster frequencyΩ µ . In particular, the temporal behavior of the i th os-cillator of the µ th cluster ( µ = 1 , . . . , M ) is given by φ µi ( t ) = Ω µ t + ρ µi + s µi ( t ) where ρ µi ∈ [0 , π ) and s µi are bounded functions describing different types of phaseclusters characterized by the phase relation within eachcluster [14]. Various types of multicluster states includ-ing the special subclass of solitary states have been ex-tensively described for adaptively coupled phase oscilla-tors [10, 15, 49]. For phase oscillator models with inertia,however, only one type of multicluster state is known sofar which is the in-phase multicluster [41, 51, 55]. More-over, little is known about the characteristic features thatstabilize the cluster.In Figure 1(a,c), we present a 4-cluster state of in-phasesynchronous clusters on a globally coupled network. As (a) (d)(c)(e) (f)(b) I nd e x i Index j φ j h ˙ φ j i Index j χ ij σ − σ FIG. 1. Hierarchical multicluster states in networks of cou-pled phase oscillators with inertia. The panels (a,b), (c,d)and (e,f) show the temporally averaged phase velocities (cid:104) ˙ φ j (cid:105) ,phase snapshots φ j ( t ) and the pseudo coupling matrices χ ij ( t ), respectively, at t = 10000. In (e) the oscillator in-dices are sorted in increasing order of their mean phase ve-locity. The states were found by numerical integration of (3)with identical oscillators P i = 0, h ( φ ) = − σγ sin( φ + α ), anduniform random initial conditions φ i (0) ∈ (0 , π ), ψ i (0) ∈ ( − . , . α is a phase-lag of the inter-action [66]. Parameters: (a,c,e) globally coupled network, M = 1, γ = 0 . σ = 0 . α = 0 . π ; (b,d,f) nonlocallycoupled ring network with coupling radius P = 40, M = 1, γ = 0 . σ = 0 . α = 0 . π ; N = 100. we know from the findings for adaptive networks, (hier-archical) multicluster states are built out of single clus-ter states whose frequency scales approximately with thenumber N µ of elements in the cluster. We find that in thezeroth-order expansion in γ the collective cluster frequen-cies are given by Ω µ ≈ − σN µ sin α , see also [14, 76]. Em-ploying the theory developed for adaptive networks [14],we find that multicluster states exist in the asymptoticlimit ( γ →
0) also for networks of phase oscillators withinertia if the cluster frequencies are sufficiently differentmeaning the clusters are hierarchical in size. Remarkably,the pseudo coupling matrix displayed in Fig. 1(e) showsthe characteristic block diagonal shape that is known foradaptive networks. In particular, the oscillators withineach cluster are more strongly connected than the oscil-lators between different clusters. I nd e x i (a) (b)Index j (c) Index j (d) I nd e x i FIG. 2. Temporal evolution of the pseudo coupling matrix χ ij for the phase oscillator model with inertia. The pseudocoupling matrix is shown at (a) t = 100, (b) t = 1750, (c) t = 5000, and (d) t = 10000 and the color code is chosen asin Fig. 1. Starting from an incoherent state in panel (a), thelargest cluster is formed first (b), and the other clusters arethen successively formed depending on their size. All param-eters and the initial condition as in Fig. 1(a,c,e). In order to extend the zoo of possible building blocksfor multiclusters, we consider a single splay-type clusterwhere φ j = 2 πkj/N with wavenumber k ∈ N . Splaystates are characterized by the vanishing local order pa-rameter R j = | (cid:80) Nk =1 a jk exp(i φ k ) | = 0 [67–69]. Takingthis latter property into account, we find that the clusterfrequency is Ω = 0. Hence, there is no way to scale thecluster frequency by scaling the cluster size. In order toovercome this, we use the findings on local splay stateson nonlocally coupled networks where each node is cou-pled to all nodes within a certain coupling range P , i.e., a ij = 1 if 0 < ( i − j ) mod N ≤ P and a ij = 0 other-wise. As it was shown in [49], splay states on nonlocallycoupled rings might lead to a non-vanishing local orderparameter, and hence, to scalable cluster frequencies.In Figure 1(b,d,f), we present a hierarchical mixed-type multicluster on a nonlocally coupled ring of phaseoscillators with inertia. It consists of one large splay clus-ter with wavenumber k = 2 and a small in-phase cluster.The emergence of such a multicluster state breaks the di-hedral symmetry of the nonlocally coupled ring network.This symmetry breaking causes a slight deviation fromthe ideal phase distribution of a splay state φ j = 2 πkj/N and of an in-phase state φ j = constant (Fig. 1d). We notethat to the best of our knowledge this type of multiclusterstate has not yet been reported in networks of coupledphase oscillators with inertia.Another novel observation for multicluster states innetworks of phase oscillators with inertia is their hierar-chical emergence. As reported in [10], the clusters emergein a temporal sequence from the largest to the smallestone. In Fig. 2, we show that this particular feature isalso found in phase oscillators with inertia.Summarizing the first example, we have shown thatthe findings for multicluster states for adaptively cou-pled phase oscillators can be transferred to networks ofphase oscillators with inertia. Note that the systems con-sidered above are homogeneous. However, heterogeneousreal-world networks can be treated with the methods es-tablished in this Letter as well. To show this, we haveanalyzed the dynamical characteristics of solitary statesin the German ultra-high voltage power grid network [70].Our pseudo coupling approach allows for an in-depth de-scription of the power flow for each line. It shows theemergence of high power fluctuations at the solitary nodeand the spreading of those fluctuations over the powergrid, see the Supplemental Material for more details [76].In the next example, we show how phenomena knownfrom power grid networks can be transferred to adaptivenetworks. Cascading line failures in adaptive networks:
In thesecond example, we show that the dynamical cascadingof line failures, which has been observed in power grids,may also occur in adaptive networks of phase oscilla-tors (1)–(2). We also propose a possible interpretationwith respect to neuronal networks.We use the following set-up that has been already em-ployed in the context of power grid networks [43]. Let usinterpret the power flow on a line in a power grid as the(localized) synaptic input F ij ( t ) = κ ij ( t ) f ( φ i ( t ) − φ j ( t ))from oscillator j to oscillator i . We say that a linefails if the corresponding synaptic input exceeds a cer-tain threshold K ∈ R , i.e., | F ij ( t ) | > K at some time t .Correspondingly, the link is cut off, i.e., a ij = a ji = 0. Apossible neuronal interpretation of such a temporal cut-off may be related to the presence of short term synapticplasticity[71, 72]. Indeed, when the signal between neu-rons or neuronal regions exceeds a certain critical level,then the corresponding connections can be affected byshort term activity-dependent depression. As a result,such an activity implies an effective cut-off, at least tem-porarily. Hence, allowing for a line to fail temporarily canbe regarded as including short term activity-dependentsynaptic depression into the adaptive network model.For the simulation presented in Fig. 3, we inte-grate (1)–(2) numerically with f ( φ ) = − φ and aHebbian-like adaptation rule g ( φ ) = − cos φ for a 5-node network. The natural frequencies are chosen as( ω , . . . , ω ) = ( − . , . , − . , − . , .
8) where we inter-pret node 2 and 5 as highly active neuronal units (hubs)due to their positive frequency. Note that in (1)–(2) allfrequencies are relative, and they can be considered asdeviations from some mean value. As initial conditionfor the simulation we choose a locally stable equilibriumthat solves 0 = ω i − (cid:80) j =1 a ij cos( φ i − φ j ) sin( φ i − φ j ), Time L i n e f a il u r e s | F i j ( t ) | (a)(b)FIG. 3. Cascading line failures in a network of 5 adap-tively coupled phase oscillators (1)–(2) with f ( φ ) = − φ , g ( φ ) = − cos φ , ( ω , . . . , ω ) = ( − . , . , − . , − . , . (cid:15) = 0 .
01. Panel (a) shows the connectivity structure ofthe nodes before (left inset) and after the cascading line fail-ures (right inset), and the number of line failures occurringduring the numerical simulation is presented vs time. Theline which is cut off due to an abrupt line failure at t = 0 . F ij during the cascading line failures. The colors of linescorrespond to the colors of the links displayed in (a). As ini-tial condition we choose a stable equilibrium of (1)–(2) withconnectivity structure as shown in the left inset of (a). Thedotted horizontal line indicates the threshold K = 1 .
04 abovewhich the lines fail. The dashed horizontal lines correspondto the stable equilibria after the link 2 − i = 1 , . . . ,
5. The corresponding initial coupling matrixtakes the form κ ij = cos( φ i − φ j ). Hence, F ij ( t ) = F ji ( t )at any time of the simulation due to the choice of thesymmetry preserving adaptation rule.We run the simulation first for 0 . t = 0 . | F ij ( t ) | , which showsthe dynamic occurrence of cascading line failures afterthe topological perturbation of the network. As a result,almost all lines fail. Moreover, the line failures lead tocomplete isolation of the hubs.It is important to note that the failure results from thedynamical properties of the network. In fact, there ex-ists a locally stable equilibrium for the initial networkwithout the connection between node 2 and 4. Thisequilibrium would meet the requirement | F ij | < K , seedashed lines in Fig. 3(b). However, this equilibrium isnot reached during the transient phase after the initialperturbation. Conclusions.
In conclusion, we find a striking relationbetween phase oscillators with inertia, which are widelyused for modeling power grids [31, 32, 34, 36, 38, 42, 43,45], and adaptive networks of phase oscillators, whichhave ubiquitous applications in physical, biological, so-cioeconomic or neuronal systems. The introduction ofthe pseudo coupling matrix allows us to split the totalinput from all nodes into node i into power flows. Thusthe frequency deviation ψ i = ˙ φ i − ω i in the phase os-cillator model with inertia corresponds to the adaptivelyadjusted total input which an oscillator receives. Thisgives insight into the concept of phase oscillator mod-els with inertia, which effectively takes into account thefeedback loop of self-adjusted coupling with all other os-cillators. Additionally, our novel theoretical frameworkallows for a generalization to swing equations with volt-age dynamics [57] and to a large class of second-orderconsensus models [73, 74, 76].Our first example shows that the theory of buildingblocks developed for adaptively coupled phase oscilla-tors can be transferred to explain the emergence of aplethora of known as well as novel multicluster states innetworks of coupled phase oscillators with inertia. In thesecond example, motivated by previous findings on thedynamical cascading of line failures in power grid net-works, we demonstrate an analogous effect in networksof adaptively coupled oscillators. While the implicationsof this effect have already been known for years in thecontext of power grids, cascading network failures havejust recently been proposed to be important for patho-logical neuronal states like Alzheimer disease [75]. Ourinsights might trigger new modeling approaches to obtaina better understanding of the function and the dysfunc-tion of the human brain.This work was supported by the German Re-search Foundation DFG, Project Nos. 411803875 and440145547. ∗ [email protected][1] H. Markram, J. L¨ubke, M. Frotscher, and B. Sakmann,Science , 213 (1997).[2] L. F. Abbott and S. Nelson, Nat. Neurosci. , 1178(2000).[3] N. Caporale and Y. Dan, Annu. Rev. Neurosci. , 25(2008).[4] C. Meisel and T. Gross, Phys. Rev. E , 061917 (2009).[5] S. Jain and S. Krishna, Proc. Natl. Acad. Sci. , 543(2001).[6] T. Gross, C. J. D. D’Lima, and B. Blasius, Phys. Rev.Lett. , 208701 (2006).[7] T. Gross and B. Blasius, J. R. Soc. Interface , 259(2008).[8] R. Guti´errez, A. Amann, S. Assenza, J. G´omez-Garde˜nes, V. Latora, and S. Boccaletti, Phys. Rev. Lett. , 234103 (2011).[9] X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Phys. Rev.Lett. , 038701 (2015).[10] D. V. Kasatkin, S. Yanchuk, E. Sch¨oll, and V. I. Neko-rkin, Phys. Rev. E , 062211 (2017).[11] M. M. Asl, A. Valizadeh, and P. A. Tass, Front. Physiol. , 1849 (2018).[12] D. V. Kasatkin and V. I. Nekorkin, Chaos , 093115(2018).[13] D. V. Kasatkin and V. I. Nekorkin, Eur. Phys. J. Spec.Top. , 1051 (2018).[14] R. Berner, E. Sch¨oll, and S. Yanchuk, SIAM J. Appl.Dyn. Syst. , 2227 (2019).[15] R. Berner, J. Fialkowski, D. V. Kasatkin, V. I. Nekorkin,S. Yanchuk, and E. Sch¨oll, Chaos , 103134 (2019).[16] R. Berner, J. Sawicki, and E. Sch¨oll, Phys. Rev. Lett. , 088301 (2020).[17] P. Feketa, A. Schaum, and T. Meurer, arXiv:1912.03922(2019).[18] O. V. Popovych, M. N. Xenakis, and P. A. Tass, PLoSONE , e0117205 (2015).[19] L. L¨ucken, O. V. Popovych, P. Tass, and S. Yanchuk,Phys. Rev. E , 032210 (2016).[20] S. Chakravartula, P. Indic, B. Sundaram, and T. Killing-back, PLoS ONE , e0178975 (2017).[21] V. R¨ohr, R. Berner, E. L. Lameu, O. V. Popovych, andS. Yanchuk, PLoS ONE , e0225094 (2019).[22] A. R. Bergen and D. J. Hill, IEEE T. Power Apparatusand Syst. , 25 (1981).[23] F. M. A. Salam, J. E. Marsden, and P. P. Varaiya, IEEETrans. Circuits Syst. , 673 (1984).[24] P. W. Sauer and M. A. Pai, Power system dynamics andstability , Vol. 101 (Upper Saddle River, NJ: Prentice hall,1998).[25] G. Filatrella, A. H. Nielsen, and N. F. Pedersen, Eur.Phys. J. B , 485 (2008).[26] J. Schiffer, D. Zonetti, R. Ortega, and A. M. Stankovic,Automatica , 135 (2016).[27] F. D¨orfler, M. Chertkov, and F. Bullo, Proc. Natl. Acad.Sci. U.S.A. , 2005 (2013).[28] T. Weckesser, H. Johannsson, and J. Ostergaard, in (IEEE, Piscataway, NJ, USA,2013).[29] T. Nishikawa and A. E. Motter, New J. Phys. , 015012(2015).[30] S. Auer, K. Kleis, P. Schultz, J. Kurths, and F. Hell-mann, Eur. Phys. J. Spec. Top. , 609 (2016).[31] F. D¨orfler and F. Bullo, SIAM J. Control Optim. ,1616 (2012).[32] M. Rohden, A. Sorge, M. Timme, and D. Witthaut,Phys. Rev. Lett. , 064101 (2012).[33] S. P. Cornelius, W. L. Kath, and A. E. Motter, Nat.Commun. , 1942 (2013).[34] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa,Nat. Phys. , 191 (2013).[35] F. D¨orfler and F. Bullo, Automatica , 1539 (2014).[36] P. J. Menck, J. Heitzig, J. Kurths, and H. J. Schellnhu-ber, Nat. Commun. , 3969 (2014).[37] F. A. Rodrigues, T. K. D. M. Peron, P. Ji, and J. Kurths,Phys. Rep. , 1 (2016).[38] D. Witthaut, M. Rohden, X. Zhang, S. Hallerberg, andM. Timme, Phys. Rev. Lett. , 138701 (2016). [39] S. Auer, F. Hellmann, M. Krause, and J. Kurths, Chaos , 127003 (2017).[40] V. Mehrmann, R. Morandin, S. Olmi, and E. Sch¨oll,Chaos , 101102 (2018).[41] P. Jaros, S. Brezetsky, R. Levchenko, D. Dudkowski,T. Kapitaniak, and Y. Maistrenko, Chaos , 011103(2018).[42] B. Sch¨afer, C. Beck, K. Aihara, D. Witthaut, andM. Timme, Nature Energy , 119 (2018).[43] B. Sch¨afer, D. Witthaut, M. Timme, and V. Latora,Nat. Commun. , 1975 (2018).[44] H. Taher, S. Olmi, and E. Sch¨oll, Phys. Rev. E ,062306 (2019).[45] F. Hellmann, P. Schultz, P. Jaros, R. Levchenko, T. Kap-itaniak, J. Kurths, and Y. Maistrenko, Nat. Commun. , 592 (2020).[46] C. H. Totz, S. Olmi, and E. Sch¨oll, arXiv preprintarXiv:1908.11649 (2019).[47] X. Zhang, C. Ma, and M. Timme, Chaos , 063111(2020).[48] F. Molnar, T. Nishikawa, and A. E. Motter, Nat. Phys. , 351 (2020).[49] R. Berner, A. Polanska, E. Sch¨oll, and S. Yanchuk, Eur.Phys. J. Spec. Top. (2020), arXiv: 1911.00320.[50] I. Belykh, B. N. Brister, and V. N. Belykh, Chaos ,094822 (2016).[51] L. Tumash, S. Olmi, and E. Sch¨oll, Chaos , 123105(2019).[52] S. Olmi, Chaos , 123125 (2015).[53] S. Olmi, A. Navas, S. Boccaletti, and A. Torcini, Phys.Rev. E , 042905 (2014).[54] J. Barr´e and D. M´etivier, Phys. Rev. Lett. , 214102(2016).[55] L. Tumash, S. Olmi, and E. Sch¨oll, Europhys. Lett. ,20001 (2018).[56] Y. Yang and A. E. Motter, Phys. Rev. Lett. , 248302(2017).[57] K. Schmietendorf, J. Peinke, R. Friedrich, andO. Kamps, Eur. Phys. J. Spec. Top. , 2577 (2014).[58] J. A. Acebr´on and R. Spigler, Phys. Rev. Lett. , 2229(1998).[59] J. A. Acebr´on, L. L. Bonilla, C. J. P´erez Vicente, F. Ri-tort, and R. Spigler, Rev. Mod. Phys. , 137 (2005).[60] D. Taylor, E. Ott, and J. G. Restrepo, Phys. Rev. E ,046214 (2010).[61] P. S. Skardal, D. Taylor, and J. G. Restrepo, Physica D , 27 (2013).[62] Y. Kuramoto, Chemical Oscillations, Waves and Turbu-lence (Springer-Verlag, Berlin, 1984).[63] W. C. Abraham and M. F. Bear, Trends Neurosci. ,126 (1996).[64] W. C. Abraham, Nat. Rev. Neurosci. , 387 (2008).[65] R. A. John, F. Liu, N. A. Chien, M. R. Kulkarni, C. Zhu,Q. D. Fu, A. Basu, Z. Liu, and N. Mathews, Adv. Mater. , 1800220 (2018).[66] H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys ,576 (1986).[67] J. G. Restrepo, E. Ott, and B. Hunt, Phys. Rev. E ,036151 (2005).[68] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, andC. Zhou, Phys. Rep. , 93 (2008).[69] M. Schr¨oder, M. Timme, and D. Witthaut, Chaos ,073119 (2017).[70] J. Egerer, Open Source Electricity Model for Ger- many (ELMOD-DE) , Tech. Rep. (Deutsches Institut f¨urWirtschaftsforschung (DIW), 2016).[71] M. H. Hennig, Front. Comput. Neurosci. , 45 (2013).[72] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski, Neuronal Dynamics: From single neurons to networksand models of cognition (Cambridge University Press,2014).[73] W. Yu, G. Chen, M. Cao, and J. Kurths, IEEE Trans.Sys. Man Cyb. , 881 (2010).[74] W. Ren, R. W. Beard, and E. M. Atkins, in Proceedingsof the 2005, American Control Conference (IEEE, 2005).[75] D. T. Jones, D. S. Knopman, J. L. Gunter, J. Graff-Radford, P. Vemuri, B. F. Boeve, R. C. Petersen, M. W.Weiner, and C. R. Jack, Jr., Brain , 547 (2016).[76] See Supplemental Material below for an analytical andnumerical analysis of a solitary state in the Germanultra-high voltage power grid network; a rigorous proofof the dynamical equivalence between the phase oscilla-tor network with inertia and the corresponding pseudoadaptive network model and between the swing equationwith voltage dynamics and the pseudo adaptive networkmodel with meta-plasticity; the pseudo adaptive couplingscheme for second-order consensus models; and an ap-proximation of the cluster frequencies in a multiclusterstate.
Supplemental Material on:What adaptive neuronal networks teach us about power grids
Rico Berner , , ∗ , Serhiy Yanchuk , and Eckehard Sch¨oll Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, Hardenbergstr. 36, 10623 Berlin, Germany Institut f¨ur Mathematik, Technische Universit¨at Berlin, Straße des 17. Juni 136, 10623 Berlin, Germany
I. SOLITARY STATES IN THE GERMANULTRA-HIGH VOLTAGE POWER GRIDNETWORK
In this section, we show that multifrequency-clusterstates, as discussed for Fig. 1 in the main text, may alsooccur in real world power grid networks. For the simula-tion, we consider the Kuramoto model with inertia givenby Eq. (3) of the main text where we set the couplingfunction as h ( φ ) = − σ sin φ . The network structure andthe power distribution are taken from the ELMOD-DEdata set provided in [70].In Figure S1, we provide a visualization of the Ger-man ultra-high voltage power grid. In order to de-termine the net power consumption/generation P i foreach node in Fig S1(a), the individual power generationand consumption for each unit are compared. We get P i = ( P total /C Total ) C off ,i − P off ,i where P Total = 36 GWand C Total ≈ .
343 GW are the off-peak power con-sumption and generation of the whole power grid net-work, respectively, and C off ,i and P off ,i are the off-peakpower consumption and generation for each individualunit, respectively. We further fix the parameters inmodel (3) of the main text as follows: M = Iω G with I = 40 × kg m and ω G = 2 π
50 Hz, γ = M a with a = 2 Hz. For details regarding the realistic modelingand the restrictions we refer the reader to [36, 44].In Figure S2 we show a solitary state obtained by thesimulation of model (3) with the parameters as describedabove for t = 600 and uniformly distributed random ini-tial conditions φ ∈ (0 , π ), ˙ φ ∈ ( − , µ ( µ = 1 , . . . , M ) where M isthe number of clusters. Within the cluster the oscilla-tors’ temporal behavior is the same up to some boundedvariations, i.e., φ µi ( t ) = Ω µ t + s µi ( t ) µ = 1 , . . . , M,i = 1 , . . . , N µ (S1)where the bounded function s µi ( t ) describes the i th os-cillator of the µ th cluster that has a total number of N µ oscillators. The temporal averages of the oscilla-tors’ phase velocities are obtained by neglecting the tran-sient period t ∈ [0 , ≈ − .
407 Hz. Similar results have been recently obtained in [44, 45]. Remarkably, the mean phase veloc-ities of the solitary nodes is very close to their naturalfrequency, see Fig. S2(b). This means that the solitarystates decouple on average from the mean field of theirneighborhood, i.e., ˙ φ Solitary = ω Solitary + (cid:80) j a ij χ ij withtemporal average (cid:104) (cid:80) j a ij χ ij (cid:105) small compared to ω Solitary .In order to shed light on further characteristics of thesolitary states, we consider the power flows, i.e, the el-ements of the pseudo coupling matrix χ ij introducedin (S8) of the main text. In Figure S2(c), we display thetemporal evolution of two typical elements of the cou-pling matrix. In both cases, the coupling between twonodes from the coherent cluster (orange) and between anode from the coherent cluster and a solitary node (blue),the average coupling value between the nodes is non zero.Both coupling weights vary periodically in time but withdifferent amplitudes. For the coupling between two nodesof the coherent cluster, the small variations stem fromthe small difference in their individual bounded tempo-ral dynamics which depends on their natural frequenciesand the individual topological neighborhoods. Due tothe realistic set-up the dynamical network is very het-erogeneous. In contrast to the case of two nodes of thecoherent clusters, the coupling between the solitary nodeand a node from the coherent cluster possesses a muchhigher variation in time and changes periodically. To un-derstand this observation, we derive an asymptotic ap-proximation for the dynamics of the solitary states.Using the approach similar to [14], we apply a mul-tiscale ansatz in (cid:15) = 1 /K (cid:28) φ repre-sents the phase of the solitary node with natural fre-quency ω and φ represents the phase of the coherentcluster with natural frequency Ω . The pseudo cou-pling weights between the two clusters are denoted by χ µν ( µ, ν = 1 , µ (cid:54) = ν ). The ansatz reads φ µ =Ω µ ( τ , τ , . . . )+ (cid:15) ( φ (1) ) µ + · · · and χ µν = χ (0) µν + (cid:15)χ (1) µν + · · · with τ p = (cid:15) p t , p ∈ N .Omitting technical details, in the first order approx-imation in (cid:15) we obtain φ = ω − ( K/ω ) cos( ωt ) and φ = Ω + ( K/ω ) cos( ωt ). Additional corrections to theoscillator frequencies appear in the third and higher or-ders of the expansion in (cid:15) and depend explicitly on ω .The latter fact is consistent with the numerical observa-tion in Fig. S2(b) that solitary nodes with a lower naturalfrequency may differ more strongly from their own nat-ural frequency than the solitary oscillators with a highernatural frequency. Moreover, the corrections explain a L o n g i t ud e ( d e g ) (a) Latitude (deg)(b) Net powers (MW) n FIG. S1. (a) Map of the German ultra-high voltage powergrid consisting of 95 net generators (green squares) and 343net consumers (red dots) connected by 662 bidirectional trans-mission lines (black lines). (b) The histogram shows the dis-tribution of the net power P i for each generator (green) andconsumer (red). The data displayed in panel (a) and (b) aretaken from the ELMOD-DE data set reported in Ref. [70] shift of average coupling weights away from zero sincesome of the coefficients χ ( p )0 do not vanish necessarily.From the approximation, we additionally obtain thatthe coupling weights between the solitary node and thecoherent cluster oscillates with the amplitude K/ω (upto the first order in (cid:15) ). Using the values obtained bythe numerical simulation, we obtain that χ , ≈ .
76 cos((16 . / π ) t ) which agrees with Fig. S2(c). Theshift, however, can be qualitatively explained only bycontributions of higher order expansions in (cid:15) as well asby the presence of the high degree of heterogeneity in the h ˙ φ i i ( s ) Index i Time χ i j (a)(c) 21410 7220 h ˙ φ i i ( s ) ω i = P i /γ (1/s)(b)FIG. S2. (a) Mean phase velocity (cid:104) ˙ φ i (cid:105) for each node inthe German power grid network presented in Fig. S1(a). (b)Mean phase velocity (cid:104) ˙ φ i (cid:105) vs. natural frequency ω i = P i /γ for each node in the German ultra-high voltage power grid.The dashed line shows the relation (cid:104) ˙ φ i (cid:105) = ω i . (c) Temporalevolution of two typical elements from the pseudo couplingmatrix χ introduced in (S8) of the main text. One elementcorresponds to a link between a solitary node and a nodefrom the synchronous cluster (blue) and the other elementcorresponds to a link between two nodes of the synchronouscluster (orange). Set-up: σ = 800 MW, t = 600, uniformlydistributed random initial conditions φ ∈ (0 , π ), ˙ φ ∈ ( − , real power grid set-up.In Figure S3, we provide a complete overview with re-spect to the pseudo coupling matrix for the results ob-tained in the simulation of the German ultra-high volt- (a) (b) L o n g i t ud e ( d e g ) (f)Latitude (deg)Latitude (deg)(c) (e)(d)Latitude (deg)Latitude (deg) L o n g i t ud e ( d e g ) FIG. S3. Panels (a) and (b) show a map of the German ultra-high voltage power grid as given in Fig. S1. Net generators andconsumers are displayed as green squares and red dots, respectively. The solitary nodes presented in Fig. S2 are displayed inblue. For each transmission line, we display in panel (a) and (b) the average pseudo coupling weight (cid:104) χ ij ( t ) (cid:105) and the maximaltemporal variation of χ ij ( t ), i.e., max( χ ij ( t )) − min( χ ij ( t )), respectively. The temporal evolution is evaluated over an averagingwindow of 100 time units. All values are coded by the line width in panel (a) and (b). In panel (b), we plot a dashed line if thetemporal variation is smaller than 0 .
05. Panels (c) and (d) provide blow-ups of panel (a) and (b) for the solitary node i = 235(red shading in all respective panels), respectively. Panels (e) and (f) provide blow-ups of panel (a) and (b) for the solitarynode i = 214 (blue shading in all respective panels), respectively. All parameters are as in Fig. S2. age power grid, see also Fig. S2. We present the averagecoupling weights as well as their temporal variations inFigure S3(a) and (b), respectively. As we know from thediscussion in the main text, the coupling weights corre-spond to the dynamics of the power flow of each trans-mission line. We further know from the asymptotic the-ory given above that the average value of the power flowbetween a solitary node and a node from the coherentcluster is small but not necessarily zero. This in fact issupported by Fig. S3(a),(c),(e). More insightful are the temporal variations of the power flow. Here, only a fewlines in Fig. S3(b) show significant temporal variations.In particular, these lines are between solitary nodes andthe coherent cluster, which is in agreement with the re-sults from the asymptotic approach and Fig. S2(c). Theblow-ups provided in Fig. S3(d,f) support the latter ob-servation by showing the highest values of the temporalvariation of the power flow for lines from and to the soli-tary nodes. Besides, Fig. S2(c) shows how far into thenetwork power fluctuations are spread in the presence ofsolitary states. It is visible that even between nodes ofthe coherent cluster high power fluctuations exist. Thesefluctuations would not be present if all oscillators weresynchronized. II. EQUIVALENCE OF TWO CLASSES OFPHASE OSCILLATOR MODELS
In the following, we show that the class of phase os-cillator models with inertia is a subclass of the phaseoscillator models with adaptive coupling weights. Forthis, we rewrite equation (3) as a system of first orderdifferential equations˙ φ i = ω i + ψ i , (S2)˙ ψ i = − γM ψ i − γ N (cid:88) j =1 a ij h ( φ i − φ j ) , (S3)where we set ω i = P i /γ .Consider the N dimensional pseudo coupling matrix χ with entries χ ij ∈ R . Let the instantaneous fre-quency deviation ψ i from the natural frequency ω i begiven by the row sum of the pseudo coupling matrix,i.e., ψ i = (cid:80) Nj =1 a ij χ ij . Then the equations (S2)–(S3) areequivalently written as the following system of adaptivelycoupled phase oscillators˙ φ i = ω i + N (cid:88) j =1 a ij χ ij , (S4)˙ χ ij = − (cid:15) ( χ ij + g ( φ i − φ j )) , (S5)with coupling function f ( φ i − φ j ) ≡ (cid:15) = γ/M and g ( φ i − φ j ) ≡ − h ( φ i − φ j ) /γ . See also the discussionof (6)–(7) in the main text.In the following, we show the dynamical equivalencebetween (S2)–(S3) and (S4)–(S5). We introduce the fol-lowing notation a i = ( a i , . . . , a iN ) , χ i = ( χ i , . . . , χ iN ) , and the N × N matrix B = a . . .
00 . . . . . . ...... . . . . . . 00 · · · a N . Define further the N variables ( ζ , . . . , ζ N ) T = B ( χ , . . . , χ N ). The kernel of the matrix B is given byker( B ) = χ ∈ R N : N (cid:88) j =1 a ij χ ij = 0 , ∀ i ∈ { , . . . , N } and is N − N dimensional in case a i (cid:54) = 0 for all i = 1 , . . . , N . Let the N − N variables ξ m ( m =1 , . . . , N − N ) define a coordinate system of ker( B ).Then the equations (S4)–(S5) can be written as˙ φ i = ω i + ζ i , (S6)˙ ζ i = − γM ζ i − γ N (cid:88) j =1 a ij h ( φ i − φ j ) , (S7)˙ ξ m = − γM ξ m + h m ( φ , . . . , φ N ) , (S8)where h m are functions determined by the particularchoice for the variables ξ m . The functions h m dependonly on the variables φ j . Thus Eq. (S8) describes thedynamics of N − N slave variables governed by the dy-namics of the 2 N dimensional system (S6)–(S7).We note that given ζ i , the N × N pseudo couplingmatrix χ is not uniquely defined. However, with theparticular choice above we are able to provide a phys-ical meaning of the coupling weights χ ij . For this, con-sider the power flows F ij from node j to node i given by F ij = g ( φ i − φ j ) [43]. Assume that each χ ij is driven bythe power flow from j to i and that asymptotically for t → ∞ χ ij → F ij in equilibrium. Then χ ij inherits themeaning of a dynamical power flow. III. SWING EQUATION WITH VOLTAGEDYNAMICS
In this section we show that the swing equation withvoltage dynamics can be regarded as a model of adap-tively coupled phase oscillators with meta-adaptivity [44,57]. The swing equation reads: M i ¨ φ i + γ ˙ φ i = P i + N (cid:88) j =1 a ij E i E j h ( φ i − φ j ) (S9) m i ˙ E i = − E i + E f,i + N (cid:88) j =1 a ij E j v ( φ i − φ j ) (S10)where the additional dynamical variable E i is the voltageamplitude, M i are the inertia constants, γ is the damp-ing constant, P i is the power of each oscillator and the2 π -periodic functions h and v are the coupling functions.The values for m i and E f,i take other machine parame-ters into account [44, 57]. All other variables and param-eters are defined as for (1)–(2) in the main text. We notethat the phase space of (S9)–(S10) is 3 N dimensional.By using the technique developed in the main text andsection II of this supplemental material, we may rewrite(S9)–(S10) as˙ φ i = ω i + N (cid:88) j =1 a ij χ ij , (S11)˙ χ ij = − M i ( γχ ij − E i E j h ( φ i − φ j )) , (S12) m i ˙ E i = − E i + E f,i + N (cid:88) j =1 a ij E j v ( φ i − φ j ) , (S13)where we introduce the coordinate changes χ ij → χ ij + P i /γ , E i → E i + E f,i and set ω i = P i /γ . Due to thevoltage dynamics (S10) the adaptation function g ( φ ) = E i ( t ) E j ( t ) h ( φ ) in (S12) possesses additional adaptivity.This kind of meta-adaptivity (meta-plasticity) has beenshown to be of importance in neuronal networks [63, 64]as well as for neuromorphic devices [65]. IV. SECOND-ORDER CONSENSUS MODELS
In this section we show that a second-order consen-sus model can be formally written as a dynamical net-work with adaptive complex coupling scheme. Consensusdescribes the result of a decision making process of au-tonomous mobile agents with positions x i and velocities v i . The decision making process is described by the con-sensus protocol that is given as a dynamical system on acomplex network structure. Consensus is achieved if theagents synchronize as time tends to infinity. Consensusmodels have a wide range of applications and are of par-ticular importance in social science and engineering [74].Let us consider the following second-order consensusmodel [73]˙ x i = v i , (S14)˙ v i = ρ N (cid:88) j =1 l ij v j + σ N (cid:88) j =1 a ij h ( x i − x j ) , (S15)where the dynamical variables x i , v i ∈ R d , a ij are theentries of the adjacency matrix of the network, l ij theentries of the Laplacian matrix of the network, i.e., l ij = a ij for i (cid:54) = j , l ii = − (cid:80) Nj =1 ,j (cid:54) = i a ij , and ρ, σ ∈ R arecoupling constants. Let us introduce the vector-valuedpseudo coupling matrix χ ij ∈ R d by v i = (cid:80) Nj =1 a ij χ ij .Using this the model (S14)–(S15) can be written as˙ x i = N (cid:88) j =1 a ij χ ij , (S16)˙ χ ij = − ρl ii χ ij + ρ N (cid:88) k =1 a jk χ jk + σ h ( x i − x j ) . (S17)By using the same arguments as given in section II,the dynamical equivalence between both models (S14)–(S15) and (S16)–(S17) can be proved. With this, we have shown that a second-order consensus model can bewritten as a dynamical network with a complex adap-tive coupling scheme rather than a fixed coupling ma-trix. Note that the elements of the complex dynamicalcoupling scheme χ ij are not uniquely defined, but mightbe chosen regarding their physical meaning. V. CLUSTER FREQUENCIES IN GLOBALLYCOUPLED PHASE OSCILLATOR MODELSWITH INERTIA
In this section, we give an approximation of the clusterfrequencies in multicluster states for N globally coupledphase oscillators with inertia. Let us consider the modelused in Fig. 1 of the main text. Then the correspondingadaptive network is as follows, (6)-(7) in the main text,˙ φ i = σ N (cid:88) j =1 χ ij , (S18)˙ χ ij = − γ ( χ ij + sin( φ i − φ j + α )) , (S19)with γ (cid:28) χ ij (cid:55)→ σχ ij . For simplicity, we assume that the whole set ofoscillators divides into two phase-locked groups of oscil-lators φ µi with µ = 1 , i = 1 , . . . , N µ , N + N = N . Inthe same way χ splits up into four blocks χ µνij describ-ing the coupling between ( µ (cid:54) = ν ) and within the clusters( µ = ν ). Note that this ansatz can be generalized to anynumber of clusters M . Assume further that the motion ofthe phase oscillators is approximately given by a commoncluster frequency, i.e., φ µi ≈ Ω µ t + ϑ µi with ϑ µi ∈ [0 , π ).Substituting this ansatz, we findΩ µ = − σN µ R ( φ µ ) sin( φ i − ψ ( φ µ ) + α ) + σ N ν (cid:88) j =1 χ µνij , where we have used that χ µµij = − sin( φ µi − φ µj + α ) withinthe cluster and the definition of the complex mean field Z ( φ µ ) and the real order parameter R ( φ µ ) for a vectorof phases φ µ = ( φ µ , . . . , φ µN µ ) T Z ( φ µ ) = N µ (cid:88) j =1 e i φ j = R ( φ µ ) e i ψ ( φ µ ) . (S20)It can be shown that if the relative motion between theclusters is sufficiently large, i.e., Ω − Ω (cid:29) γ , then thecoupling weights between the clusters χ ij , χ ij scale with γ and can be approximately neglected. For rigorous re-sults, we refer to Ref. [14]. Finally, we findΩ µ = − σN µ R ( φ µ ) sin( φ i − ψ ( φ µ ) + α )as the zeroth approximation in γ . Hence, in-phase clus-ters ( R ( φ µ ) = 1) have a collective frequency Ω µ = − σN µ sin( α ). Thus, the frequency difference can be con-trolled by scaling the size of the individual clusters. Incase of a splay configuration within each cluster, i.e., φ µi = 2 πk µ i/N µ = with wave number k µ ∈ N , the or-der parameter vanishes ( R ( φ µ ) = 0) and hence Ω µµ