Birth and stabilization of phase clusters by multiplexing of adaptive networks
BBirth and stabilization of phase clusters by multiplexing of adaptive networks
Rico Berner , , ∗ Jakub Sawicki , 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, Hardenbergstr. 36, 10623 Berlin, Germany (Dated: October 1, 2019)We propose a concept to generate and stabilize diverse partial synchronization patterns (phaseclusters) in adaptive networks which are widespread in neuro- and social sciences, as well as biology,engineering, and other disciplines. We show by theoretical analysis and computer simulations thatmultiplexing in a multi-layer network with symmetry can induce various stable phase cluster statesin a situation where they are not stable or do not even exist in the single layer. Further, we developa method for the analysis of Laplacian matrices of multiplex networks which allows for insight intothe spectral structure of these networks enabling a reduction to the stability problem of single layers.We employ the multiplex decomposition to provide analytic results for the stability of the multilayerpatterns. As local dynamics we use the paradigmatic Kuramoto phase oscillator, which is a simplegeneric model and has been successfully applied in the modeling of synchronization phenomena ina wide range of natural and technological systems.
Complex networks are an ubiquitous paradigm in na-ture and technology, with a wide field of applicationsranging from physics, chemistry, biology, neuroscience,to engineering and socio-economic systems. Of particu-lar interest are adaptive networks, where the connectivitychanges in time, for instance, the synaptic connectionsbetween neurons are adapted depending on the relativetiming of neuronal spiking [1–5]. Similarly, chemical sys-tems have been reported [6], where the reaction ratesadapt dynamically depending on the variables of the sys-tem. Activity-dependent plasticity is also common inepidemics [7] and in biological or social systems [8]. Syn-chronization is an important feature of the dynamics innetworks of coupled nonlinear oscillators [9–13]. Vari-ous synchronization patterns are known, like cluster syn-chronization where the network splits into groups of syn-chronous elements [14], or partial synchronization pat-terns like chimera states where the system splits into co-existing domains of coherent (synchronized) and incoher-ent (desynchronized) states [15–17]. These patterns werealso explored in adaptive networks [18–33]. Furthermore,adapting the network topology has also successfully beenused to control cluster synchronization in delay-couplednetworks [34].Another focus of recent research in network science aremultilayer networks, which are systems interconnectedthrough different types of links [35–38]. A prominentexample are social networks which can be described asgroups of people with different patterns of contacts orinteractions between them [39–41]. Other applicationsare communication, supply, and transportation networks,for instance power grids, subway networks, or airtrafficnetworks [42]. In neuroscience, multilayer networks rep-resent for instance neurons in different areas of the brain,neurons connected either by a chemical link or by an elec-trical synapsis, or the modular connectivity structure ofbrain regions [43–51]. A special case of multilayer net-works are multiplex topologies, where each layer contains the same set of nodes, and only pairwise connections be-tween corresponding nodes from neighbouring layers ex-ist [52–71].In spite of the lively interest in the topic of adaptivenetworks, little is known about the interplay of adap-tively coupled groups of networks [25, 72, 73]. Suchadaptive multilayer or multiplex networks appear nat-urally in neuronal networks, e.g., in interacting neuronpopulations with plastic synapses but different plasticityrules within each population [74, 75], or affected by dif-ferent mechanisms of plasticity [76], or the transport ofmetabolic resources [77]. Beyond brain networks, coex-isting forms of (meta)plasticity are investigated in neuro-inspired devices to develop artificially intelligent learningcircuitry [78].In this Letter we show that a plethora of novel pat-terns can be generated by multiplexing adaptive net-works. In particular, partial synchronization patternslike phase clusters and more complex cluster states whichare unstable in the corresponding monoplex network canbe stabilized, or even states which do not exist in thesingle-layer case for the parameters chosen, can be bornby multiplexing. Thus our aim is to provide fundamentalinsight into the combined action of adaptivity and multi-plex topologies. Hereby we elucidate the delicate balanceof adaptation and multiplexing which is a feature of manyreal-world networks even beyond neuroscience [79–82].As local dynamics we use the paradigmatic Kuramotophase oscillator model, which is a simple generic modeland has been successfully applied in the modeling of syn-chronization phenomena in a wide range of natural andtechnological systems [13].A general multiplex network with L layers each con-sisting of N identical adaptively coupled phase oscillators a r X i v : . [ n li n . AO ] S e p is described by˙ φ µi = ω − N N (cid:88) j =1 κ µij sin( φ µi − φ µj + α µµ ) − L (cid:88) ν =1 ,ν (cid:54) = µ σ µν sin( φ µi − φ νi + α µν ) , (1)˙ κ µij = − (cid:15) (cid:0) κ µij + sin( φ µi − φ µj + β µ ) (cid:1) , where φ µi ∈ [0 , π ) represents the phase of the i th oscil-lator ( i = 1 , . . . , N ) in the µ th layer ( µ = 1 , . . . , L ), and ω is the natural frequency. The interaction between thephase oscillators within each layer is described by thecoupling matrix elements κ µij ∈ [ − , κ µij are determined adaptively, whereasthe inter-layer coupling weights σ µν ≥ α µν are the phase lags of the interaction [83].The adaptation rate 0 < (cid:15) (cid:28) β µ of the adaptation functionsin( φ µi − φ µj + β µ ) describes different rules that may oc-cur in neuronal networks. For instance, for β µ = − π/
2, aHebbian-like rule [84–86] is obtained where the coupling κ ij is increasing between any two systems with close-byphases, i.e., φ i − φ j ≈ β = 0, the link κ ij willbe strengthened if the i th oscillator is advancing the j th .Such a relationship is typical for spike-timing dependentplasticity in neuroscience [3, 5, 88, 89]. If β = π/
2, ananti-Hebbian-like rule emerges.Eq. (S1) has been widely used as a paradigmaticmodel for adaptive networks [18–30]. It generalizes theKuramoto-Sakaguchi model with fixed coupling topol-ogy [90–94].Let us note important properties of the model.First, ω can be set to zero without loss of general-ity due to the shift-symmetry of Eq. (S1), i.e. , con-sidering the co-rotating frame φ → φ + ωt . More-over, due to the existence of the attracting region G ≡ (cid:8)(cid:0) φ µi , κ µij (cid:1) : φ µi ∈ (0 , π ] , | κ µij | ≤ , i, j = 1 , . . . , N,µ = 1 , . . . , L } , one can restrict the range of the couplingweights to the interval − ≤ κ ij ≤ α , β , φ , κ ) (cid:55)→ ( − α , π − β , − φ , κ ) , ( α µµ , β µ , φ µi , κ µij ) (cid:55)→ ( α µµ + π, β µ + π, φ µi , − κ µij ) , where α , β , φ , κ abbreviate the whole set of variables andparameters, it is sufficient to analyze the system withinthe parameter region α ∈ [0 , π/ α µµ ∈ [0 , π ) ( µ (cid:54) = 1), α µν ∈ [0 , π ) ( µ (cid:54) = ν ) and β µ ∈ [ − π, π ).For a single layer, Eq. (S1) has been studied numeri-cally and analytically [18–22, 24]. In particular, it wasshown that starting from uniformly distributed random initial condition φ i ∈ [0 , π ), κ ij ∈ [ − ,
1] the system canreach different frequency multi-cluster states with hier-archical structure depending on the parameters α and β . The frequency multi-clusters in turn consist of severalone-clusters which determine the existence and stabilityof the former [23]. Therefore, these one-cluster states(with identical frequency, but different phase distribu-tions) constitute the building blocks of adaptively cou-pled phase oscillators.Before we consider multiple layers, we will introducenotions and patterns of the single-layer case, which willbe needed for the multiplex network. Each solution ofEq. (S1) for L = 1 or L = 2 is called a monoplex orduplex state, respectively.It is known that already the monoplex system (S1)possesses a huge variety of dynamical states such asmulti-clusters with respect to frequency synchronization,chaotic attractors, and chimera-like states [20–24]. Inthis article we will focus on the simplest patterns, namelyone-cluster states, and their generalization to the multi-plex case. The reasons for this are threefold. First, fromthe analytical point of view, the one-cluster states in themonoplex network are very well understood. Second, theone-clusters are building blocks for more complex dynam-ical states such as multi-cluster states. Therefore, theyare essential for the understanding of more complex dy-namical patterns. And third, stable chimera-like statesas they were studied in [24, 25] exist close to the bordersof the stability regions for one-clusters, so the existenceand stability of one-clusters may pave the way for ob-serving those hybrid patterns.In general, one-cluster states are given by equilibriarelative to a co-rotating frame [22] φ µi = Ω t + a µi ,κ µij = − sin( a µi − a µj + β µ ) , (2)with collective frequency Ω and relative phases a µi .Hence the second moment order parameter R ( a µ ) = N (cid:12)(cid:12)(cid:12)(cid:80) Nj =1 e i2 a µj (cid:12)(cid:12)(cid:12) with a µ ≡ ( a µ , . . . , a µN ) T can be used asa characteristic measure. In the case of monoplex sys-tems ( L = 1), three types of solutions exist (see Fig. 1)which are characterized by corresponding frequencies Ωas a function of ( α , β ) [22]: (a) Ω = cos( α − β ) / R ( a ) = 0 (Splay state), (b) Ω = sin α sin β if R ( a ) = 1 with a i ∈ { , π } (Antipodal state), (c)Ω = cos( α − β ) / − R ( a ) cos( ψ Q ) / < R ( a ) < a i ∈ { , π, ψ Q , ψ Q + π } (Double antipodal state)with ψ Q being the unique solution (modulo 2 π ) of(1 − q ) sin( ψ Q − α − β ) = q sin( ψ Q + α + β ) , where q = Q/N and Q ∈ { , . . . , N − } denotes thenumber of phase shifts a i ∈ { , π } . Here, splay statesare defined in a more general sense by R ( a ) = 0, whichincludes the states a i = 2 πi/N usually referred to assplay state [95]. κ ij − ψ (a) (b) (c) FIG. 1. Illustration of the three types of monoplex one-clusterstates of Eq. (2) ( L = 1) for an ensemble of 10 oscillators(green circles) with frequencies Ω (upper panels) and couplingstructure with weights κ ij (lower panels): One-cluster (a) ofsplay type ( R ( a ) = 0), (b) of antipodal type ( R ( a ) = 1),and (c) of double antipodal type with Q = 7. Parameters: α = 0 . π , β = 0 . π Let us now consider these one-cluster states in multi-plex structures. Therefore, we introduce the notion of lifted one-cluster states, where in each layer µ = 1 , . . . , L the state ( φ µi ( t ) , κ µij ( t )) is a monoplex one-cluster, i.e.,the phases a µi of the oscillators are of splay, antipodal,or double antipodal type. It can be shown [102] that induplex systems ( L = 2) the phase difference of oscillatorsbetween the layers ∆ a ≡ a i − a i takes only two valuesand solves ∆Ω = σ sin(∆ a + α ) + σ sin(∆ a − α ),where ∆Ω ≡ Ω( α , β ) − Ω( α , β ) is given above for thethree different one-cluster states (splay, antipodal, dou-ble antipodal). These states are shown in Fig. 2: Panels(a),(b),(d) display lifted states of splay, antipodal, andsplay type, respectively. The phase distributions in bothlayers are the same but shifted by the constant value ∆ a in agreement with the above equation. In contrast to thelifted states, Fig. 2(c) shows another possible one-clusterfor the duplex network. Due to the interaction of the twolayers we can find a phase distribution which is of doubleantipodal type in each layer but not a lifted state. Thismeans that these states are born by the duplex set-up.Moreover, in contrast to the other examples the phasedistribution between the layers does not agree, ψ (cid:54) = ψ .For the monoplex case, it has been shown that double an-tipodal states are unstable for any set of parameters [23].Hence, finding stable double antipodal states which in-teract through the duplex structure is unexpected.For more insight into the birth of phase-locked statesby multiplexing, Fig. 3 displays the emergence of doubleantipodal states in a parameter regime where they donot exist in single-layer networks. They are character-ized by the second moment order parameter R . It isremarkable that the new double antipodal state can befound for a wide range of the inter-layer coupling strengthlarger than a certain critical value σ c , and is clearly differ-ent from those of the monoplex. Below the critical value σ c , the double antipodal states are no longer stable, andmore complex temporal dynamics occurs which causes Layer 1 Layer 2 Layer 1 Layer 2 π π ∆ a ∆ aπ πψ ψ
10 30 50 10 30 50 10 30 50 10 30 50 j j φ µ j i φ µ j i (a)1030501030502 π π π π (b)(c) (d) FIG. 2. Different duplex states of Eq. (2) ( L = 2) for anensemble of 50 oscillators in each layer with color-coded cou-pling weights κ µij (upper panels, color code as in Fig.1), phases φ µj (lower panels): Duplex one-cluster states (a) of liftedsplay type ( R ( a µ ) = 0) for α / = 0 . π , σ / = 0 . R ( a µ ) = 1) for α = 0 . π , α = 0 . π , σ / = 0 .
62; (c) of double antidodal type (nota lifted state) for α / = 0 . π , σ / = 0 .
28; (d) of liftedsplay type for α = 0 . π , α = 0 . π , σ / = 0 .
8, and (cid:15) = 0 .
01. In the lower panels phase differences between thetwo layers are indicated by ∆ a ≡ a i − a i , and between thetwo new antipodal states (c) by ψ , ψ . temporal changes in R . This leads to non-vanishingtemporal variance indicated by the error bars in Fig. 3.In the following we show how the dynamics in a neigh-borhood of theses states can be lifted as well, i.e., we in-vestigate their local stability. Everything is exemplifiedfor antipodal states but can be generalized to lifted splayand lifted double antipodal states in a straightforwardmanner. To study the dynamics around the one-clusterstates described by Eq. (2), we linearize Eq. (S1) aroundthese states:˙ δφ µi = 1 N N (cid:88) j =1 (cid:2) sin(∆ a + β µ ) cos(∆ a + α µµ )∆ µµij δφ − sin(∆ a + α µµ ) δκ µij (cid:3) − M (cid:88) ν =1 σ µν cos(∆ a + α µν )∆ µνij δφ, ˙ δκ µij = − (cid:15) (cid:0) δκ µij + cos(∆ a + β µ )∆ µµij δφ (cid:1) (3)where ∆ µνij δφ ≡ δφ µi − δφ νj .In duplex networks, the coupling structure is given bya 2 × M with the N × N unity matrix I N : M = (cid:18) A m · I N n · I N B (cid:19) . (4)If A and B are diagonalizable N × N matrices whichcommute ( m, n ∈ R , n (cid:54) = 0), the following relation for φ µ j i R π π . . . .
200 0 . . π ψ ψ π j h R ( φ ) ih R ( φ ) i σ c σ FIG. 3. Birth of double antipodal state in a duplex network( N = 12) for a wide range of inter-layer coupling strength σ = σ = σ . The solid lines are the temporal averagesfor the second moment order parameter R of the individuallayers (layer 1: black, layer 2: red). The error bars for σ < σ c denote the standard deviation of the temporal evolution of R . The dashed horizontal lines represent the unique valuesof R for the double antipodal state in a monoplex network.The plot was obtained by adiabatic continuation of a duplexdouble antipodal state (see inset) in both directions startingfrom σ = 0 .
5. Parameters: α / = 0 . π , α / = 0 . β = 0 . π , β = − . π , and (cid:15) = 0 . the characteristic polynomial can be proven [102] usingSchur’s decomposition [96, 97]:det ( mn · I N − ( D A − λ I N )( D B − λ I N )) = 0 , (5)where D A and D B are the diagonal matrices correspond-ing to A and B , respectively. Note that Eq. (5) not onlysimplifies the calculation for the eigenvalues in the caseof a duplex structure, moreover, it is a general result onlinear dynamical systems on duplex networks. Therefore,this result is important for the investigation of stabilityand symmetry in multiplex networks.In the case of a duplex antipodal one-cluster stateEq. (S1) with a i ∈ { , π } and a i = a i − ∆ a ,Eq. (3) can be brought to the form (4) and pos-sesses the following set of Lyapunov exponents S = {− (cid:15), ( λ i, , λ i, , λ i, , λ i, ) i =1 ,...,N } where λ i, ,..., are thesolutions of polynomials containing the eigenvalues of themonoplex system [102].Thus, the stability analysis of the duplex system canbe reduced to that of the monoplex case. We are able toanalyze the stabilizing and destabilizing features of a du-plex network numerically and analytically. To illustratethe effect of multiplexing, the interaction between twoclusters of antipodal type is presented in Fig. 4. The sta-bility of these states is determined by integrating Eq. (2)numerically starting with a slightly perturbed lifted an-tipodal state. The states are stable if the numerical tra-jectory is approaching the lifted antipodal state. Oth-erwise, the state is considered as unstable. The blackcontour lines in Fig. 4 show the borders of the stabil-ity regions in dependence of the coupling strength σ ,as calculated from the Lyapunov exponents.The bordersare in remarkable agreement with the numerical results. σ σ = 0 . σ = 0 β /π β /π α / π (a) (b) FIG. 4. Regions of stability (blue) and instability (white)of the lifted antipodal state in the ( α , β ) parameter planefor different values of interlayer coupling (indicated by differ-ent blue shading) σ , where regions of stronger coupling σ (lighter blue) include such of weaker σ (darker blue). Stabil-ity regions for single-layer antipodal clusters are indicated byred hatched areas. The inter-layer coupling is considered as(a) unidirectional ( σ = 0) and (b) bidirectional ( σ = σ ).Parameters: α = 0 . π , β = − . π , α = 0, α = 0 . π ,and (cid:15) = 0 . We investigate two different situations in Fig. 4: Inboth panels the parameters for the first layer α , β arechosen such that the antipodal state is stable withoutinter-layer coupling. The stability of the duplex antipo-dal states is displayed in the ( α , β ) parameter planefor several values of the inter-layer coupling σ . To com-pare the effects of the duplex network with the mono-layer case, the stability regions for monoplex antipodalsstates are displayed, as well, as red hatched areas. Theyare markedly different. In Figure 4(a), the two layers areconnected unidirectionally ( σ = 0). It can be seen thatwith increasing inter-layer coupling weight σ the regionof stability for the lifted antipodal state also grows. Al-ready for small values of the inter-layer couplings σ , astabilizing effect of the duplex network can be noticed.For σ = 0 . σ similar to the unidirectional case. However, the regionsof stability grow at different rates in dependence on σ and non-monotonically with respect to the parameters α , β . Comparing the size of the stability region forboth cases, one can see that for small values of σ theregion for bidirectional coupling is larger. In turn, forhigher inter-layer coupling, the regions for the unidirec-tional case are larger. It is worth noting that in Fig. 4the stability regions for smaller values of σ are alwayscontained in the region for larger values of σ .In conclusion, we have proposed a concept to inducediverse partial synchronization patterns (phase clusters)in adaptively coupled phase oscillator networks. Whileadaptive networks have recently attracted a lot of atten-tion in the fields of neuro- and social sciences, biology,engineering, and other disciplines, and multilayer net-works are a paradigm for real-world complex networks,little has been known about the interplay of multilayerstructures and adaptivity. We have aimed to fill this gapwithin a rigorous framework of theoretical analysis andcomputer simulations. We have shown that multiplex-ing in a multi-layer with symmetry can induce variousstable phase cluster states like splay states, antipodalstates, and double antipodal states, in a situation wherethey are not stable or do not even exist in the singlelayer. Further, we have developed a novel method foranalysis of Laplacian matrices of duplex networks whichallows for insight into the spectral structure of these net-works, and can easily be generalized to more than twolayers [102]. This new approach of multiplex decompo-sition has a broad range of applications to physical, bio-logical, socio-economic, and technological systems, rang-ing from plasticity in neurodynamics or the dynamicsof linear diffusive systems [98, 99] to generalizations ofthe master stability approach [100, 101] for adaptive net-works [102]. We have used the multiplex decompositionto provide analytic results for the stability of lifted statesin the multilayer system. As local dynamics we haveused the paradigmatic Kuramoto phase oscillator model,supplemented by adaptivity of the link strengths with aphase lag parameter which can model a whole range ofadaptivity rules from Hebbian via spike-timing depen-dent plasticity to anti-Hebbian.This work was supported by the German ResearchFoundation DFG (Projects SCHO 307/15-1 and YA225/3-1 and Projektnummer - 163436311 - SFB 910). Wethank Serhiy Yanchuk for insightful discussions. ∗ [email protected] † [email protected][1] H. Markram, J. L¨ubke, 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] L. L¨ucken, O. Popovych, P. Tass, and S. Yanchuk, Phys.Rev. E , 032210 (2016).[6] S. Jain and S. Krishna, Proc. Natl. Acad. Sci. , 543(2001).[7] T. Gross, C. J. D. D’Lima, and B. Blasius, Phys. Rev.Lett. , 208701 (2006).[8] T. Gross and B. Blasius, J. R. Soc. Interface , 259(2008).[9] A. Pikovsky, M. G. Rosenblum, and J. Kurths, Syn- chronization: a universal concept in nonlinear sciences (Cambridge University Press, Cambridge, 2001).[10] S. H. Strogatz, Nature , 268 (2001).[11] R. Albert and A. L. Barab´asi, Rev. Mod. Phys. , 47(2002).[12] M. E. J. Newman, SIAM Review , 167 (2003).[13] S. Boccaletti, A. N. Pisarchik, C. I. del Genio, andA. Amann, Synchronization: From Coupled Systems toComplex Networks (Cambridge University Press, Cam-bridge, 2018).[14] T. Dahms, J. Lehnert, and E. Sch¨oll, Phys. Rev. E ,016202 (2012).[15] Y. Kuramoto and D. Battogtokh, Nonlin. Phen. in Com-plex Sys. , 380 (2002).[16] D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. ,174102 (2004).[17] M. J. Panaggio and D. M. Abrams, Nonlinearity , R67(2015).[18] T. Aoki and T. Aoyagi, Phys. Rev. Lett. , 034101(2009).[19] T. Aoki and T. Aoyagi, Phys. Rev. E , 066109 (2011).[20] V. I. Nekorkin and D. V. Kasatkin, AIP Conf. Proc. , 210010 (2016).[21] D. V. Kasatkin and V. I. Nekorkin, Radiophysics andQuantum Electronics , 877 (2016).[22] R. Berner, E. Sch¨oll, and S. Yanchuk, “Multi-clusters innetworks of adaptively coupled phase oscillators,” SIAMJ. Appl. Dyn. Syst., in print (2019), arXiv:1809.00573.[23] R. Berner, J. Fialkowski, D. V. Kasatkin, V. Nekorkin,E. Sch¨oll, and S. Yanchuk, “Hierarchical frequency clus-ters in adaptive networks of phase oscillators,” Chaos,in print (2019), arXiv:1904.06927.[24] D. V. Kasatkin, S. Yanchuk, E. Sch¨oll, and V. I. Neko-rkin, Phys. Rev. E , 062211 (2017).[25] D. V. Kasatkin and V. I. Nekorkin, Chaos , 1054(2018).[26] A. Gushchin, E. Mallada, and A. Tang, in Informa-tion Theory and Applications Workshop ITA 2015, SanDiego, CA, USA (IEEE, 2015) pp. 291–300.[27] C. B. Picallo and H. Riecke, Phys. Rev. E , 036206(2011).[28] L. Timms and L. Q. English, Phys. Rev. E , 032906(2014).[29] Q. Ren and J. Zhao, Phys. Rev. E , 016207 (2007).[30] V. Avalos-Gayt´an, J. A. Almendral, I. Leyva, F. Battis-ton, V. Nicosia, V. Latora, and S. Boccaletti, Phys. Rev.E , 042301 (2018).[31] L. Papadopoulos, J. Z. Kim, J. Kurths, and D. S. Bas-sett, Chaos An Interdiscip. J. Nonlinear Sci. , 073115(2017).[32] D. V. Kasatkin and V. I. Nekorkin, Eur. Phys. J. Spec.Top. , 1051 (2018).[33] D. V. Kasatkin, V. V. Klinshov, and V. I. Nekorkin,Phys. Rev. E , 1 (2019).[34] J. Lehnert, P. H¨ovel, A. A. Selivanov, A. L. Fradkov, andE. Sch¨oll, Phys. Rev. E , 042914 (2014).[35] S. Boccaletti, G. Bianconi, R. Criado, C. I. del Ge-nio, J. G´omez-Garde˜nes, M. Romance, I. Sendi˜na Nadal,Z. Wang, and M. Zanin, Phys. Rep. , 1 (2014).[36] M. De Domenico, A. Sol´e-Ribalta, E. Cozzo, M. Kivel¨a,Y. Moreno, M. A. Porter, S. G´omez, and A. Arenas,Phys. Rev. X , 041022 (2013).[37] M. De Domenico, V. Nicosia, A. Arenas, and V. Latora,Nat. Commun. , 6864 (2015). [38] M. Kivel¨a, A. Arenas, M. Barth´elemy, J. P. Gleeson,Y. Moreno, and M. A. Porter, J. Complex Networks , 203 (2014).[39] M. Girvan and M. E. J. Newman, Proc. Natl. Acad. Sci.USA , 7821 (2002).[40] R. Amato, N. E. Kouvaris, M. S. Miguel, and A. D´ıaz-Guilera, New J. Phys. , 123019 (2017).[41] R. Amato, A. D´ıaz-Guilera, and K.-K. Kleineberg, Sci.Rep. , 7087 (2017).[42] A. Cardillo, M. Zanin, J. G`omez Garde˜nes, M. Romance,A. Garcia del Amo, and S. Boccaletti, Eur. Phys. J. ST , 23 (2013).[43] D. Meunier, R. Lambiotte, and E. T. Bullmore, Front.Neurosci. , 200 (2010).[44] B. Bentley, R. Branicky, C. L. Barnes, Y. L. Chew,E. Yemini, E. T. Bullmore, P. E. V´etes, and W. R.Schafer, PLOS Comput. Biol. , 1 (2016).[45] F. Battiston, V. Nicosia, M. Chavez, and V. Latora,Chaos , 047404 (2017).[46] M. Vaiana and S. F. Muldoon, J. Nonlinear Sci. , 1(2018).[47] L. Ramlow, J. Sawicki, A. Zakharova, J. Hlinka, J. C.Claussen, and E. Sch¨oll, EPL , 50007 (2019).[48] A. Ashourvan, Q. K. Telesford, T. Verstynen, J. M. Vet-tel, and D. S. Bassett, PLoS One , e0215520 (2019).[49] C. Zhou, L. Zemanov´a, G. Zamora, C. C. Hilgetag, andJ. Kurths, Phys. Rev. Lett. , 238103 (2006).[50] C. Zhou, L. Zemanov´a, G. Zamora-L´opez, C. C. Hilgetag,and J. Kurths, New J. Phys. , 178 (2007).[51] R. Wang, P. Lin, M. Liu, Y. Wu, T. Zhou, and C. Zhou,Phys. Rev. Lett. , 38301 (2019).[52] X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Phys. Rev.Lett. , 038701 (2015).[53] V. A. Maksimenko, V. V. Makarov, B. K. Bera,D. Ghosh, S. K. Dana, M. V. Goremyko, N. S. Frolov,A. A. Koronovskii, and A. E. Hramov, Phys. Rev. E ,052205 (2016).[54] S. Jalan and A. Singh, Europhys. Lett. , 30002(2016).[55] S. Ghosh, A. Kumar, A. Zakharova, and S. Jalan, Eu-rophys. Lett. , 60005 (2016).[56] I. Leyva, R. Sevilla-Escoboza, I. Sendi˜na-Nadal,R. Guti´errez, J. M. Buld´u, and S. Boccaletti, Sci. Rep. , 45475 (2017).[57] R. G. Andrzejak, G. Ruzzene, and I. Malvestio, Chaos , 053114 (2017).[58] S. Ghosh, A. Zakharova, and S. Jalan, Chaos SolitonsFractals , 56 (2018).[59] N. Semenova and A. Zakharova, Chaos , 051104(2018).[60] M. Mikhaylenko, L. Ramlow, S. Jalan, and A. Za-kharova, Chaos , 023122 (2019).[61] J. Sawicki, I. Omelchenko, A. Zakharova, and E. Sch¨oll,Phys. Rev. E , 062224 (2018).[62] I. Omelchenko, T. H¨ulser, A. Zakharova, and E. Sch¨oll,Front. Appl. Math. Stat. , 67 (2019).[63] E. Rybalova, T. Vadivasova, G. Strelkova, V. An-ishchenko, and A. Zakharova, Chaos , 033134 (2019).[64] D. Nikitin, I. Omelchenko, A. Zakharova, M. Avetyan,A. L. Fradkov, and E. Sch¨oll, Phil. Trans. R. Soc. A , 20180128 (2019).[65] K. A. Blaha, K. Huang, F. Della Rossa, L. M. Pecora,M. Hossein-Zadeh, and F. Sorrentino, Phys. Rev. Lett. , 014101 (2019). [66] R. Sevilla-Escoboza, I. Sendi˜na-Nadal, I. Leyva,R. Guti´errez, J. M. Buld´u, and S. Boccaletti, ChaosAn Interdiscip. J. Nonlinear Sci. , 065304 (2016).[67] R. J. Requejo and A. D´ıaz-Guilera, Phys. Rev. E ,022301 (2016).[68] E. Pitsik, V. Makarov, D. Kirsanov, N. Frolov, M. Gore-myko, X. Li, Z. Wang, A. Hramov, and S. Boccaletti,New J. Phys. , 075004 (2018).[69] I. Leyva, I. Sendi˜na-Nadal, R. Sevilla-Escoboza, V. P.Vera-Avila, P. Chholak, and S. Boccaletti, Sci. Rep. ,8629 (2018).[70] S. Jalan, A. Kumar, and I. Leyva, Chaos An Interdiscip.J. Nonlinear Sci. , 041102 (2019).[71] N. S. Frolov, V. A. Maksimenko, V. V. Makarov, D. V.Kirsanov, A. E. Hramov, and J. Kurths, Phys. Rev. E , 022320 (2018).[72] O. V. Maslennikov and V. I. Nekorkin, Chaos , 121101(2018).[73] M. V. Goremyko, D. V. Kirsanov, V. O. Nedaivozov,V. V. Makarov, and A. E. Hramov, in Dyn. Fluctu-ations Biomed. Photonics XIV , Vol. 10063, edited byV. V. Tuchin, K. V. Larin, M. J. Leahy, and R. K.Wang (2017) p. 100631C.[74] A. Citri and R. C. Malenka, Neuropsychopharmacology , 18 (2008).[75] E. Edelmann, E. Cepeda-Prado, and V. Leßmann, Front.Synaptic Neurosci. , 032805 (2017).[76] F. Zenke, E. J. Agnes, and W. Gerstner, Nat. Commun. , 6922 (2015).[77] Y. S. Virkar, W. L. Shew, J. G. Restrepo, and E. Ott,Phys. Rev. E , 042310 (2016).[78] R. A. John, F. Liu, N. A. Chien, M. R. Kulkarni, C. Zhu,Q. Fu, A. Basu, Z. Liu, and N. Mathews, Adv. Mater. , 1800220 (2018).[79] T. Gross, C. J. D. D’Lima, and B. Blasius, Phys. Rev.Lett. , 208701 (2006).[80] S. Shai and S. Dobson, Phys. Rev. E , 042812 (2013).[81] L. Wardil and C. Hauert, Sci. Rep. , 5725 (2015).[82] P. Klimek, M. Diakonova, V. M. Egu´ıluz, M. San Miguel,and S. Thurner, New J. Phys. , 083045 (2016).[83] H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys ,576 (1986).[84] F. C. Hoppensteadt and E. M. Izhikevich, Biol. Cybern. , 129 (1996).[85] P. Seliger, S. C. Young, and L. S. Tsimring, Phys. Rev.E , 041906 (2002).[86] T. Aoki, Neural Networks , 11 (2015).[87] D. Hebb, The Organization of Behavior: A Neuropsycho-logical Theory , new edition ed. (Wiley, New York, 1949).[88] Y. Maistrenko, B. Lysyansky, C. Hauptmann,O. Burylko, and P. A. Tass, Phys. Rev. E ,066207 (2007).[89] O. Popovych, S. Yanchuk, and P. Tass, Sci. Rep. , 2926(2013).[90] Y. Kuramoto, Chemical Oscillations, Waves and Turbu-lence (Springer-Verlag, Berlin, 1984).[91] S. H. Strogatz, Physica D , 1 (2000).[92] J. A. Acebr´on, L. L. Bonilla, C. J. P. Vicente, F. Ritort,and R. Spigler, Rev. Mod. Phys. , 137 (2005).[93] A. Pikovsky and M. G. Rosenblum, Phys. Rev. Lett. ,264103 (2008).[94] O. E. Omel’chenko and M. Wolfrum, Phys. Rev. Lett. , 164101 (2012).[95] C. U. Choe, T. Dahms, P. H¨ovel, and E. Sch¨oll, Phys. Rev. E , 025205(R) (2010).[96] S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004).[97] J. Liesen and V. Mehrmann,
Linear Algebra , SpringerUndergraduate Mathematics Series (Springer Interna-tional Publishing, Cham, 2015).[98] S. G´omez, A. D´ıaz-Guilera, J. G´omez-Garde˜nes, C. J.P´erez Vicente, Y. Moreno, and A. Arenas, Phys. Rev.Lett. , 028701 (2013).[99] A. Sol´e-Ribalta, M. De Domenico, N. E. Kouvaris,A. D´ıaz-Guilera, S. G´omez, and A. Arenas, Phys. Rev.E , 032807 (2013).[100] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. ,2109 (1998).[101] L. Tang, X. Wu, J. L¨u, J. Lu, and R. M. D’Souza, Phys.Rev. E (2019), 012304.[102] See Supplemental Material (below) for phase differencesof one-cluster duplex states, the characteristic polyno-mial of duplex matrices, a general form of the multiplexdecomposition, a generalization of the master stabilityapproach to multiplex networks, an application to diffu-sively coupled linear systems, and for the derivation ofthe Lyapunov exponents of duplex antipodal states. Supplemental Material onBirth and stabilization of phase clusters by multiplexing of adaptive networks
Rico Berner , , ∗ , Jakub Sawicki , 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, Hardenbergstr. 36, 10623 Berlin, Germany
MODEL
We consider a multiplex network with L layers eachconsisting of N identical adaptively coupled phase oscil-lators dφ µi dt = ω − N N (cid:88) j =1 κ µij sin( φ µi − φ µj + α µµ ) (S1) − L (cid:88) ν =1 ,ν (cid:54) = µ σ µν sin( φ µi − φ νi + α µν ) ,dκ µij dt = − (cid:15) (cid:0) κ µij + sin( φ µi − φ νj + β µ ) (cid:1) , (S2)where φ µi ∈ [0 , π ) represents the phase of the i th os-cillator ( i = 1 , . . . , N ) in the µ th layer ( µ = 1 , . . . , L )and ω is the natural frequency. The interaction betweenthe phase oscillators within each layer is described by thecoupling matrix κ µij ∈ [ − , σ µν ≥ α µν can be considered as a phase lag ofthe interaction [S1].In order to include also more general dynamics onstatic networks we also consider the following dynami-cal equations with diffusive coupling [S2 ? , S3]˙ x ki = f ( x ki ( t )) − σ N (cid:88) j =1 a kij H ( x ki − x kj ) − ρ L (cid:88) l =1 m kl G ( x ki − x li ) , (S3)where x ki is a d -dimensional state vector of the i th nodein the k th layer, i.e. x ki ∈ C d , f ∈ C ( C d , C d ) describesthe local dynamics, the functions H, G ∈ C ( C d , C d ) de-termine the intra- and inter-layer coupling scheme, re-spectively. The parameter σ ∈ R is the intra-layer cou-pling strength, ρ ∈ R and is the inter-layer couplingstrength. The connectivities are given by the matrix el-ements a kij ∈ { , } of the N × N intra-layer adjacencymatrix A k and the L × L (inter)-layer adjacency matrix M with elements m kl ∈ { , } . EXISTENCE OF DUPLEX EQUILIBRIA INADAPTIVE NETWORKS
Suppose we have two one-cluster states where eachis of either splay, antipodal, or double antipodal typewhich form a duplex one-cluster (see Eq. (2) of the maintext) and φ µi = Ω( α µµ , β µ ) t + ω µ t + a µi ( µ = 1 , α µµ , β µ ) is given by Eq. (3), of the main text, and thecoupling weights are given by κ µij = − sin( a µi − a µj + β µ ).We verify by directly inserting that φ µi and κ µij solveEq. (S2). For the given ansatz, Eq. (S1) readsΩ( α , β ) + ω = 12 cos( α − β ) − (cid:60) (cid:16) e − i(2 a i + α + β ) Z ( a ) (cid:17) − σ sin(∆Ω t + ∆ ωt + a i − a i + α ) , andΩ( α , β ) + ω = 12 cos( α − β ) − (cid:60) (cid:16) e − i(2 a i + α + β ) Z ( a ) (cid:17) + σ sin(∆Ω t + ∆ ωt + a i − a i − α ) . where ∆Ω = Ω( α , β ) − Ω( α , β ) and ∆ ω = ω − ω ,respectively. Thus, φ = ( φ , φ ) is a duplex equilibriumif ∆Ω + ∆ ω = 0which is equivalent to∆Ω = σ sin( a i − a i + α )+ σ sin( a i − a i − α )for all i = 1 , . . . , N . Note that ∆Ω is not necessarily zeroeven if the phase-lag parameters for both layers agree.They can still differ in the type of one-cluster state. Theformer equation can be written as∆Ω C = sin( a i − a i + ν ) (S4)with sin( ν ) = 1 C (cid:0) σ sin( α ) − σ sin( α ) (cid:1) , (S5)cos( ν ) = 1 C (cid:0) σ cos( α ) + σ cos( α ) (cid:1) , where C = (cid:112) ( σ ) + ( σ ) + 2 σ σ cos( α + α ) . Whenever ( σ ) + ( σ ) + 2 σ σ cos( α + α ) ≥ σ ) + ( σ ) + 2 σ σ cos( α + α ) ≥ ∆Ω , (S6)Eq. (S4) has the two solutions a i − a i = arcsin(∆Ω /C ) − ν and a i − a i = π − arcsin(∆Ω /C ) − ν . Considering theinverse function arcsin : [ − , → [ − π/ , π/
2] applied toEq. (S5) determines ν to be either ν (cid:48) or π − ν (cid:48) , where ν (cid:48) := arcsin(sin( ν )) and sin( ν ) as given in (S5). Thesecond equation for cos( ν ) then fixes ν to take one of thevalues.The condition (S6) is a relation between all parametersof the system which has to be fulfilled for the existence ofduplex relative equilibria. Note that for any given inter-layer coupling σ (cid:54) = 0 and α + α (cid:54) = ± π/ ± π/ σ < ∞ suchthat the lifted one-clusters exist. In case of unidirectionalcoupling, i.e., σ = 0, the condition gives the minimumweight σ ≥ ∆Ω. MULTIPLEX NETWORKS AND THEIRDECOMPOSITION
In this section, we provide important tools and theo-rems to find the spectrum of multiplex networks.
Theorem 1.
Let R be a commutative subring of C N × N and let M ∈ R L × L . Then, det C M = det C (det R M ) . The proof can be found in Ref ? . This rather abstractresult allows for a very nice decomposition for pairwisecommuting matrices and yields a useful tool to study thelocal dynamics in multiplex systems. Proposition 2.
Let M ∈ C N × N be a unitary diago-nalizable matrix with M = U D M U H where U , U H and D M are a unitary, its adjoint and a diagonal matrix, re-spectively. Let further D M be the set of simultaneouslydiagonalizable matrices to M , i.e., the set of all matriceswhich commute pairwise and with M . Then, det A · · · A L ... . . . ... A L · · · A LL =det (cid:32) (cid:88) σ ∈ S L (cid:34) sgn ( σ ) L (cid:89) µ =1 D A µ,σ ( µ ) (cid:35)(cid:33) (S7) where A µν ∈ D M for µ, ν = 1 , . . . , L and S L is the set ofall permutations of the numbers , . . . , L . Proof. Consider any
A, B ∈ D M , then they are simulta-neously diagonalizable with M and hence A = D A U H and B = U D B U H with the same U . Thus, all A µν canbe diagonalized with the same U . Since U is unitary, i.e. (det U ) = 1, we finddet A · · · A L ... . . . ... A L · · · A LL = det D A · · · D A L ... . . . ... D A L · · · D A LL by applying the block diagonal matrices diag( U, · · · , U )and diag( U H , · · · , U H ) from the left and right, respec-tively. The set of diagonal matrices with usual matrixmultiplication and addition form a commutative sub-ring of C N × N . Applying Theorem 1 and using the well-known determinant representation of Leibniz, the expres-sion (S7) follows. Remark . The set D M consists of all matrices whichcommute with M and all the other elements of D M . Inparticular, the identity matrix I N ∈ D M for any M ∈ C N × N .In the following, we apply the last result to a duplexand triplex system and connect the local dynamics on theone-layer network to the multiplex case. We specify ourconsideration by defining two special multiplex systems. Definition 4.
Suppose
A, B, C ∈ C N × N and m ij ∈ C ( i, j = 1 , . . . , N × N block matrix M (2) = (cid:18) A m I m I B (cid:19) (S8)and the 3 N × N block matrix M (3) = A m I m I m I B m I m I m I C (S9)are called (complex) duplex and triplex network, respec-tively.Suppose we know how to diagonalize the individuallayer topologies. The next result shows how the eigenval-ues of the individual layers are connected to eigenvaluesof the multiplex system. This will be done for the duplexand triplex network. For the proof of our following state-ment, we provide two different proofs. The first approachmakes use of Schur’s decomposition [S6, S7] which will beused, later on, in order to derive the characteristic equa-tions. In particular, any m × m matrix M = (cid:18) A BC D (cid:19) inthe 2 × M = (cid:18) I p BD − I q (cid:19) (cid:18) A − BD − C D (cid:19) (cid:18) I p D − C I q (cid:19) . (S10)With this, a simplified form of the determinant of a 2 × M = (cid:18) A BC D (cid:19) is derived, namelydet( M ) = det( A − BD − C ) · det( D ) , (S11)An extension of the first approach to any number of lay-ers in the network can be found by induction but is verytechnical, see [ ? ]. The second approach uses Proposi-tion 2 which allows for a straightforward extension to anynumber of layers in a multiplex network. Proposition 5.
Suppose
A, B, C ∈ C N × N , they com-mute pairwise, and are diagonalizable with diagonal ma-trices D A , D B , D C and unitary matrix U . Then, theeigenvalues µ for the multiplex networks M (2) and M (3) can be found by solving the N quadratic µ − (( d A ) i + ( d B ) i ) µ + ( d A ) i ( d B ) i − m m = 0(S12) and cubic polynomial equations µ + a ,i µ + a ,i µ + a ,i = 0 , (S13) respectively, with a ,i = − (( d A ) i + ( d B ) i + ( d C ) i ) a ,i = ( d A ) i ( d B ) i + ( d A ) i ( d C ) i + ( d B ) i ( d D ) i − m m − m m − m m a ,i = m m ( d C ) i + m m ( d B ) i + m m ( d A ) i − ( d A ) i ( d B ) i ( d C ) i − m m m − m m m and ( d A ) i , ( d B ) i , and ( d C ) i the respective diagonal ele-ments of D A , D B , and D C .Proof. Since
A, B, C are diagonalizable and com-mute, Proposition 2 can be applied to both matrices M (2) , M (3) . Anyhow, for the matrix M (2) we will provideanother proof using Schur’s decomposition.The determinant is an antisymmetric multi-linearform. Thus, we can writedet (cid:16) M (2) − µ I N (cid:17) = det (cid:18) A − µ I N m · I N m · I N B − µ I N m (cid:19) = ( − N det (cid:18) m · I N A − µ I N B − µ I N m · I N . (cid:19) By assumption A and B are both diagonalizable withrespect to the unitary transformation matrix U , and soare A − µ I and B − µ I . This allows us to writedet (cid:18) m · I N A − µ I N B − µ I N m · I N . (cid:19) = det (cid:18) m I N D A − µ I N D B − µ I N m I N (cid:19) we apply the block diagonal matrices diag( U, · · · , U ) anddiag( U H , · · · , U H ) from the left and right, respectively.Now, using Schur’s decomposition (S10) the determinantcan written asdet (cid:18) m I N D A − µ I N D B − µ I N m I N (cid:19) = n N det (cid:18) m − n ( D A − µ I N ) ( D B − µ I N ) (cid:19) = det ( m m I N − ( D A − µ I N ) ( D B − µ I N )) . The last expression together with det (cid:0) M (2) − µ I N (cid:1) = 0yields the N quadratic equations (S12).Using that ( A − µ I ) , ( B − µ I ) , ( C − µ I ) commute pair-wise, Proposition 2 can be applied. We finddet (cid:16) M (3) − µ I N (cid:17) =det (( D A − µ I N ) [( D B − µ I N )( D C − µ I N ) − m m I N ] − m [ m ( D C − µ I N ) − m m I N ]+ m [ m m I N − m ( D B − µ I N )])The last expression together with det (cid:0) M (3) − µ I N (cid:1) = 0yields the N cubic equations (S13).Let us briefly discuss some special cases for both theduplex and triplex network. Consider a duplex networkwith master and slave layer, i.e. , either m = 0 or m =0. Then, the quadratic equations (S12) yield( µ − ( d A ) i ) ( µ − ( d B ) i ) = 0 . (S14)As shown in Proposition 2, the eigenvalues for specialtriplex networks can be found by solving cubic equations.For the solution even closed forms exist. Despite this, theexplicit form of the solutions is rather tedious, in general.However, if we consider A = B = C and a ring-likeinter-layer connection between the networks, i.e. , m = m = m = 0, then equation (S13) has the followingsolutions for all j = 1 , . . . , Nµ = − ( d A ) j + ( m m m ) / ,µ = − ( d A ) j + 12 i(i + √
3) ( m m m ) / ,µ = − ( d A ) j −
12 (i + √
3) ( m m m ) / , where i denotes the imaginary unit. In analogy to equa-tion (S14), a decoupling for the eigenvalues can be found.Consider three pairwise commuting matrices A, B, C ,and the structure between the layers is a directed chain,i.e., m = m = m = m = 0 , then( µ − ( d A ) i ) ( µ − ( d B ) i ) ( µ − ( d C ) i ) = 0 . (S15)In the following three section, different applications willbe briefly discussed. The master stability approach for multiplexnetworks
In Ref. [S9], the master stability function for dynam-ical systems on multiplex networks was introduced. Intheir article, the authors considered a systems similarto (S3) but with Laplacian instead of adjacency matri-ces. Since they assumed the coupling function H and G to be linear, both systems are equivalent. Consider nowthe synchronous solution which solves ˙ s = f ( s ). Themaster stability function corresponding to (S3) can bederived from the following variational equations˙ ξ = (cid:2) I NL ⊗ Df ( s ) − σ (cid:0) L intra ⊗ H (cid:1) − ρ (cid:0) L inter ⊗ G (cid:1)(cid:3) ξ , all necessary details can be found in Ref. [S9]. Here, ξ = x − I L ⊗ I N ⊗ s and x ∈ C L · N · d is the system statevector where all individual nodal states are stacked oneach other ordered by the layer and node index. Further,the intra-layer Laplacian is defined is defined as L intra = L (cid:77) l =1 L l = L . . . L L where L k = (cid:80) Nj =1 a k j . . . (cid:80) Nj =1 a kNj − A k . The inter-layer Laplacian is defined as L inter = L I ⊗ I N where L I = (cid:80) Ll =1 m l . . . (cid:80) Ll =1 m Ll − M, where m kl with m ll = 0 are the entries of the L × L matrix M . In Ref [S9], it is shown that if L intra and L inter commute, a master stability equation for system (S3) canbe found which reads˙ y = [ Df ( s ) − αH − βG ] y , where α = σλ , β = ρµ , λ and µ are the (complex)eigenvalues of L intra and L inter , respectively. Note that if H = G the master stability equation can be reduced to˙ y = [ Df ( s ) − γH ] y , (S16)with γ = α + β and y ∈ C d . Equation (S16) is calledthe master stability equation for the composite systemwhere a single supra-Laplacian matrix σ L intra + ρ L inter described the network topology [S10, S11].Direct evaluation shows that (cid:2) L intra , L inter (cid:3) = 0 isequivalent to m kl L l − m lk L k = 0 for all l, k = 1 , . . . , L . Thus, to require (cid:2) L intra , L inter (cid:3) = 0 yields a linear de-pendence between the individual layer topologies. UsingProposition 4 for composite system, the master stabilityfunction can be used under much milder conditions. Proposition 6.
Let us consider the multiplex dynamicalsystem (S3) , with linear function H = G . Suppose thatfurther all L l , l = 1 , . . . , L , commute pairwise. Then, themaster stability equation is given by ˙ y = [ Df ( s ) − µH ] y , (S17) and µ = µ ( λ i , . . . , λ li ) being a non-linear function, map-ping the i th eigenvalues λ li ( i = 1 , . . . , N ) of the layertopologies L l , to the ”master parameter” µ in (S17) . Thenon linear mapping is given as the formal solution to the L th order polynomial equation det (cid:32) (cid:88) σ ∈ S L (cid:34) sgn ( σ ) L (cid:89) l =1 D (cid:16) L supra lσ ( l ) − µ I N δ lσ ( l ) (cid:17) (cid:35)(cid:33) = 0 , (S18) where L supra = σ L intra + ρ L inter is a L × L block matrixdivided into N × N matrices which we individually referto with L supra kl , k, l = 1 , . . . , L , and δ kl is the Kroneckersymbol.Proof. Using H = G , the variation equation for (S3) onthe synchronous solution s ( t ) is given by˙ ξ = [ I NL ⊗ Df ( s ) − L supra ⊗ H ] ξ , By assumption all L l , l = 1 , . . . , L , commute pairwise and L inter is a L × L block matrix consisting of N × N identitymatrices multiplied by scalar. Hence, Proposition 2 canbe applied to ( L supra − µ I L · N ) in order to diagonalize L supra .With this Proposition, we have a very powerful tool inorder to investigate not only the influence of the multi-plex structure network on the stability of the synchronousstate but also the impact of different layer topologieswhich are not necessarily linearly dependent. As an ex-ample, we consider a duplex systems with [ L , L ] = 0and L supra = (cid:18) σL + ρm I − ρm I − ρm I σL + ρm I (cid:19) . Knowing the eigenvalues for L and L , the master pa-rameter µ is determined using Eq. (S12). The matrices L and L are Laplacian matrices which have at least onezero eigenvalue, corresponding to the so-called Goldstonemode ˆ1 = (1 , . . . N − times · · · , T . As a result there aretwo parameters µ = 0 and µ = ρ ( m + m ). The firstvalue corresponds to the Goldstone mode (cid:0) ˆ1 , ˆ1 (cid:1) T . Thesecond parameter is exclusively induced by the duplexstructure and completely independent from the individ-ual layer topologies. Further, we find that in case of m = − m another zero parameter µ exists which cor-responds to the eigenmode (cid:0) ˆ1 , − ˆ1 (cid:1) T . The other eigen-values of the supra-Laplacian matrix L supra are given bythe non-linear mappings µ ( λ , λ ) = σ ( λ + λ ) + ρ ( m + m )2 ± (cid:113) (cid:0) σ ( λ − λ ) + ρ ( m − m ) (cid:1) + 4 ρ m m (S19)for which we have formally solved equation (S12) withrespect to µ . Let us consider two special cases.First, we assume that there is no connection from thesecond to the first layer. We have a master-slave set-upwhich means that m = 0. With this, the master pa-rameter is µ ( λ , λ ) = σλ and µ ( λ , λ ) = σλ + ρm .Remarkably, in this set-up the stability of a synchronousstate in a duplex network is reduced to the pure one-layersystem. The stability in the duplex system is determinedby the spectrum of the individual layer topologies whereonly in the second layer the spectrum is shifted due tothe interaction.The second case starts from the consideration in [S9].In particular, we consider (cid:2) L intra , L inter (cid:3) = 0 which leadsto a pairwise linear dependence of all individual layertopologies, and hence λ li = λ i for all l = 1 , . . . , L and i = 1 , . . . , N with λ li ∈ C . Taking this into account,the equation for the master parameter yields µ ( λ , λ ) = σλ + ρ (cid:0) m + m (cid:1) and µ ( λ , λ ) = σλ . In order tosee that the master parameter agrees with the set for theparameter γ in (S16), we determine the eigenvalues of L intra and L inter individually. Since L (1) and L (2) are lin-early dependent, the set of eigenvalues for L intra consistsof the eigenvalues of L with double multiplicity. Us-ing Proposition 5, we find that L inter has eigenvalues 0and ( m + m ) each with multiplicity N . Since bothLaplacian matrices commute, there exists a common setof eigenvectors. We find γ = ρλ (1) for the eigenvector(1 , T ⊗ v and γ = ρλ (1) + σ ( m + m ) for the eigen-vector ( m , − m ) T ⊗ v where L v = λ v . Here, againeverything boils down to a pure one-layer set-up with anadditional shift. Analytic treatment of diffusive dynamics onmultiplex networks
In the previous section we considered the dynamics oflinear systems as they are given by the variational equa-tion (S16). In the context of diffusive systems on complex networks, recently linear diffusive processes were consid-ered in order to study the dynamics on social as wellas transport networks [S12, S13]. Compared with equa-tion (S3) in [S13], the authors investigate a duplex system( L = 2) with f = 0, σ = D k ( k = 1 , H, G = I , and ρm kl = D x . Additionally, they considere weighted net-works for which instead of a kij ∈ { , } coupling weights w kij ( k = 1 ,
2) were taken. Note that the former resultson multiplex matrices still hold true for weighted connec-tions.Let us assume that the super-Laplacian matrix is givenby L supra = (cid:18) D L + D x I − D x I − D x I D L + D x I (cid:19) . with [ L , L ] = 0. Due to the structure of the supra-Laplacian we are allowed to apply Proposition 6. In ac-cordance with [S13] and our former findings in Sec. , wehave the two eigenvalues µ = 0 and µ = 2 D x corre-sponding to the Goldstone mode (ˆ1 , ˆ1) T and the vector(ˆ1 , − ˆ1) T , respectively. The other eigenvalues are givenas solution to the equations (S12) and read µ i = D λ i + D λ i + 2 D x ± (cid:113) (cid:0) D λ i − D λ i (cid:1) + 4( D x ) , where λ i and λ i are the eigenvalues of the Laplacianmatrices L and L , respectively. Regions of synchronization in adaptive multiplexnetworks
For an arbitrary duplex equilibrium of the form φ µi =Ω t + a µi with a k = (0 , πN k, . . . , ( N − πN k ) T and a k = a k − ∆ a we start with the linearized system (4) of themain text. This can also be written in the block matrixform ˙ δφ ˙ δφ ˙ δκ ˙ δκ = A m I N B m I N A B C − (cid:15) I N C − (cid:15) I N δφ δφ δκ δκ with ( δφ µ , δκ µ ) T = ( δφ µ , . . . , δφ µN , δκ µ , . . . , δκ µ N ,δκ µ , . . . , δκ µNN ) T , the matrices A µ , B µ , and C µ followfrom system (4) of the main text, and m , m ∈ R . Withthe help of Schur’s decomposition the characteristic equa-tion for the linearized system takes the form( λ + (cid:15) ) N − N ) = 0det (cid:18) ( λ I N − A ) ( λ + (cid:15) ) − B C − ( λ + (cid:15) ) m I N − ( λ + (cid:15) ) m I N ( λ I N − A ) ( λ + (cid:15) ) − B C (cid:19) = 0 . (S20)The second equation has the block matrix form whichis required from Proposition 5. All blocks can be di-agonalized and commute since they all possess a cyclicstructure; compare Lemma 4.1 of [S14]. Thus, we are allowed to apply Proposition 5 which we use in order todiagonalize the matrix in Eq. (S20). For the diagonalizedmatrix we find the following equations for the diagonalelements µ i ( λ + (cid:15) ) m m − ( p i ( λ ; α , β , α , σ ) − µ i )( p i ( λ ; α , β , α , σ ) − µ i ) = 0 (S21)where i = 1 , . . . , N , p µi ( λ ; α µµ , β µ ) is a second order poly-nomial in λ which depends continuously on α and β aswell as functionally on the type of the one-cluster state.For every i ∈ { , . . . , N } , these equations will give us twoeigenvalues µ i, and µ i, for the matrix in Eq. (S20) de-pending on λ and the system parameters. Thus, we canwrite Eq. (S21) as( µ i − µ i, ( λ ; α , β , σ )) ( µ i − µ i, ( λ ; α , β , σ )) = 0where α , β , σ represent all system parameter chosen for(S1)–(S2). In order to find the eigenvalue λ of the lin- earized system (S20) one of the eigenvalues µ has tovanish. This means that we have to find λ such thatEq. (S21) equals µ i ( µ i − µ i, ( λ ; α , β , σ )) = 0which is equivalent to finding λ such that the followingquartic equation is solved p i ( λ ; α , β , α , σ ) p i ( λ ; α , β , α , σ ) − ( λ + (cid:15) ) m m = 0 . (S22)Note that here the diagonal elements of A are slightlydifferent from those in Prop. 4.2 of [S14] but they donot affect the result, i.e., the diagonal element equals ρ i ( α , β ) − m ( α ). The same holds true for A . Thus, with the two possible eigenvalues ρ i, , ( α µµ , β µ ) for themonoplex system from Corr. 4.3 of [S14] one finds thefollowing quartic equation which give the Lyapunov ex-ponents for the lifted duplex one-cluster (cid:2)(cid:0) λ − ρ i, ( α , β ) (cid:1) · (cid:0) λ − ρ i, ( α , β ) (cid:1) + m ( λ + (cid:15) ) (cid:3) × (cid:2)(cid:0) λ − ρ i, ( α , β ) (cid:1) (cid:0) λ − ρ i, ( α , β ) (cid:1) + m ( λ + (cid:15) ) (cid:3) − ( λ + (cid:15) ) m m = 0 . (S23)In case of an antipodal monoplex one-cluster given byEq. (2) of the main text with a i ∈ { , π } , then the setof Lyapunov exponents S for Eq. (4)in the main text isgiven by In the case of a duplex antipodal one-cluster state givenby Eq. (2) of the main text with a i ∈ { , π } and a i = a i − ∆ a , Eq. (3) possesses the following set of Lyapunovexponents S Duplex = {− (cid:15), ( λ i, , λ i, , λ i, , λ i, ) i =1 ,...,N } where λ i, ,..., solve the following N quartic equations( λ + (cid:15) ) m m − (cid:2)(cid:0) λ − ρ i, (cid:1) · (cid:0) λ − ρ i, (cid:1) + m ( λ + (cid:15) ) (cid:3) × (cid:2)(cid:0) λ − ρ i, (cid:1) (cid:0) λ − ρ i, (cid:1) + m ( λ + (cid:15) ) (cid:3) = 0 , (S24)with m ≡ σ cos(∆ a + α ), m ≡ σ cos(∆ a − α ) andthe eigenvalues ρ µi, , ≡ ρ i, , ( α µµ , β µ ) for the monoplexsystem S Monoplex = (cid:110) (0) , ( − (cid:15) ) ( N − N +1 , ( ρ ) N − , ( ρ ) N − (cid:9) (S25)where ρ and ρ solve ρ + (cid:0) (cid:15) − cos( α ) sin( β ) (cid:1) ρ − (cid:15) sin( α + β ) = 0. Here, the multiplicities for eacheigenvalue are given as lower case labels. The proof andthe results for other clusters can be found in [S14, S15].Using the Lyapunov exponents of the duplex antipodalclusters, the stability for duplex antipodal states can befound. The results are presented in Fig. 4 (in the maintext). ∗ [email protected] † [email protected][S1] H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys ,576 (1986). [S2] A. Keane, T. Dahms, J. Lehnert, S. Suryanarayana,P. H¨ovel, and E. Sch¨oll, Eur. Phys. J. B , 407 (2012).[S3] J. Lehnert, Controlling Synchronization Patterns inComplex Networks , Springer Theses (Springer Interna-tional Publishing, Cham, 2016).[S4] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. ,2109 (1998).[S5] J. R. Silvester, Math. Gaz. , 460 (2000).[S6] S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004).[S7] J. Liesen and V. Mehrmann,
Linear Algebra , SpringerUndergraduate Mathematics Series (Springer Interna-tional Publishing, Cham, 2015).[S8] P. D. Powell, arXiv: 1112.4379 (2011).[S9] L. Tang, X. Wu, J. L¨u, J.-a. Lu, and R. M. D’Souza,Phys. Rev. E , 012304 (2019).[S10] M. De Domenico, A. Sol´e-Ribalta, E. Cozzo, M. Kivel¨a,Y. Moreno, M. A. Porter, S. G´omez, and A. Arenas,Phys. Rev. X , 041022 (2013).[S11] M. Kivela, A. Arenas, M. Barthelemy, J. P. Gleeson,Y. Moreno, and M. A. Porter, J. Complex Networks ,203 (2014).[S12] M. Barth´elemy, Phys. Rep. , 1 (2011).[S13] S. G´omez, A. D´ıaz-Guilera, J. G´omez-Garde˜nes, C. J.P´erez-Vicente, Y. Moreno, and A. Arenas, Phys. Rev.Lett.110