Symmetries and Cluster Synchronization in Multilayer Networks
F. Della Rossa, L. Pecora, K. Blaha, A. Shirin, I. Klickstein, F. Sorrentino
SSymmetries and Cluster Synchronization in Multilayer Networks
Fabio Della Rossa , Louis Pecora , Karen Blaha , Afroza Shirin , Isaac Klickstein , andFrancesco Sorrentino University of New Mexico, Albuquerque, NM, 87131, USA U.S. Naval Research Laboratory, Washington, DC 20375, USA Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy Corresponding author [email protected] 21, 2020
Abstract
Real-world systems in epidemiology, social sciences, power transportation, economics and engineeringare often described as multilayer networks. Here we first define and compute the symmetries of multilayernetworks, and then study the emergence of cluster synchronization in these networks. We distinguishbetween independent layer symmetries which occur in one layer and are independent of the other layersand dependent layer symmetries which involve nodes in different layers. We study stability of the clustersynchronous solution by decoupling the problem into a number of independent blocks and assessingstability of each block through a Master Stability Function. We see that blocks associated with dependentlayer symmetries have a different structure than the other blocks, which affects the stability of clustersassociated with these symmetries. Finally, we validate the theory in a fully analog experiment in whichseven electronic oscillators of three kinds are connected with two kinds of coupling.
Introduction
Real-world complex systems often contain heterogeneous components; these components may interact inmultiple ways via complex connectivity patterns, leading to complex dynamics. For example, the powergrid contains both a transmission and a communication network, and we must model both to understandphenomena such as cascade failures [1, 2] and blackouts duration [3, 4]. In this context, a recent paper [3]has shown that increased coupling between the power and the communication layers can be beneficial inreducing vulnerability of the system as a whole. Neurological systems can be modeled with varying levels ofcomplexity depending on the particular behavior of interest. Some approaches have a single kind of neuronswhich interact via multiple interaction layers [5, 6, 7, 8]; central pattern generator approaches [9, 10, 11]have multiple kinds of neurons (identified by kind, such as extensor, flexor, synapses generator, etc.) whichinteract via multiple interaction layers. The mathematical formalism introduced for these types of problemsis the multilayer network [12, 13, 14].We are interested in symmetries and synchronization of multilayer networks of coupled oscillators; syn-chronization is a collective behavior in which the dynamics of the network nodes converge on the sametime-evolution. The first study of synchronization in multilayer networks was presented in [8, 15] andmore recently in [16]. These papers studied complete synchronization (all nodes synchronize on the sametime-evolution) and consider diffusive coupling (the coupling matrices are Laplacian). Populations mayexhibit more complex forms of synchronization, such as clustered synchronization (CS), where clusters ofnodes exhibit synchronized dynamics but different clusters evolve on distinct time evolutions; many pa-pers consider CS in networks formed of nodes all of the same type and connections all of the same type[17, 18, 19, 20, 21, 22, 23, 24]. 1 a r X i v : . [ n li n . C D ] J u l he emergence of synchronized clusters in a network requires that the set of the network nodes is parti-tioned into subsets of nodes, called equitable clusters, such that all the nodes in the same cluster receive thesame total input from each one of the clusters, leading to the same temporal evolution [25, 26, 27]. Thesepartitions often arise from the network symmetries [28, 29]. A recent paper studied the effects of symmetryon collective dynamics in communities of coupled oscillators [30]. When equitable partitions arise from sym-metry, a rigorous, group theory-based framework exists to analyze the stability of cluster synchronization insimple networks [20, 21, 24]. However, cluster synchronization in multilayer networks has not been studiedother than in a recent paper which focused on experimental observations on a special multilayer networkcomposed of nodes all of the same type [31].This paper significantly advances the field of network dynamics by presenting one unified theory thataddresses the problem of cluster synchronization in arbitrary multilayer networks, where each layer is formedof homogeneous units, but different layers have different units. Our main contributions are twofold: we defineand compute the symmetries of multilayer networks and we study the emergence of cluster synchronization inthese networks analytically and experimentally. This involves analyzing stability of the cluster synchronoussolution, where the stability problem is decoupled into a number of lower dimensional blocks of equations anda validation of the theory in a fully analog experiment with three layers, each one formed of different typesof oscillators. We analyze CS in arbitrary multilayer networks, for which nodes are coupled through bothintra-layer connections (inside each layer) and inter-layer connections (between different layers.) We see thatCS patterns of synchronization in multilayer networks are determined by the symmetries of each individuallayer but also the particular pattern of interconnectivity between layers. We see that only nodes from thesame layer may be permuted among each other. However, other symmetries may involve simultaneous swapsof nodes in different layers. This leads to a classification of symmetries into Independent Layer Symmetries(ILS) which occur in one layer and are independent of swaps in other layers and Dependent Layer Symmetries(DLS) which require moving nodes in different layers. In what follows, we first present a general set ofdynamical equations for the time evolution of multilayer networks, we then define and compute the groupof symmetries of multilayer networks with different kinds of nodes (each corresponding to a different layer)and different kinds of connections, for which stability of the CS solution is analyzed, and finally we applythe theoretical framework to predict the emergence of CS in an experimental system. Results
Formulation and Dynamical Equations
Multilayer networks have different types of nodes interacting through different types of connections [12, 13].Previously defined subclasses of multilayer networks include multiplex [32, 33, 34] and multidimensional [35,36] networks. In a multiplex network, the same set of agents exists in all layers; for example, in a social system,each node is a person, each layer represents opinion on a topic, and links capture how social interactionsinfluence a person’s opinion on each topic. In a multidimensional network, different kinds of links connectthe same set of nodes, all of the same type. For a deeper discussion of particular multilayer networks andhow to reconcile each case with the general definition of a multilayer network, see [12, 13]. Sometimes theterm ‘multilayer network’ has been used in the literature to indicate generic networks formed of differenttypes of nodes and/or connections. For example, in [31] a ‘multilayer network’ has connections of differenttypes but nodes all of the same type. Here we consider the most general situation for which both nodesand connections can be of different types, with the case of Ref. [31] remaining a special case of multilayernetwork.Reference [37] introduced a mathematical model for the time evolution of a multiplex network. Followingprevious studies of synchronization in different instances of multilayer networks [8, 15, 37, 16, 31], next weprovide a general set of equations that describe the dynamics of a multilayer network. First we define thesets of nodes and of connections (or interactions) of a multilayer network. The nodes are arranged in sets { X α , α = 1 , . . . , M } , where each individual set corresponds to a given layer of the multilayer network. Theuncoupled dynamics of all the N α nodes in layer X α is the same: ˙ x αi = F α ( x αi ), i = 1 , . . . , N α , x α ∈ R n α .2he multilayer network has a total number of nodes equal to N = (cid:80) α N α . The set of connections correspondto either intra-layer interactions that connect nodes in the same layer or inter-layer interactions that connectnodes in different layers. The intra-layer interactions inside layer α are described by an adjacency matrix A αα , to which is associated a nonlinear coupling function H αα : R n α (cid:55)→ R n α and a coupling strength σ αα .The inter-layer interactions from layer β to layer α are described by an N α × N β adjacency matrix A αβ ,to which is associated a nonlinear coupling function H αβ : R n α (cid:55)→ R n β and a coupling strength σ αβ . Weassume throughout that all the couplings are undirected: A αβ,λ = A βα,λT , α, β = 1 , ..., M .We show an example of a multilayer network in Fig. 1 a . This network has two layers, labeled α and β , with intra-connectivity described by the matrices A αα and A ββ and inter-connectivity described by thematrix A αβ = ( A βα ) T , all shown in b . Fig. 1 c shows an independent layer symmetry. This symmetryinvolves only nodes in layer α ; we can swap α nodes 1 with 2 and α nodes 3 with 4 without affecting layer β . Fig. 1 d shows a dependent layer symmetry. This symmetry requires swapping nodes in both layer α and β ; when we swap α nodes 2 with 3 and α nodes 1 with 4, we must also swap β nodes 1 and 3. Note thatas dependent layer symmetries involve swapping nodes of different types, they are a characteristic feature ofmultilayer networks with different types of nodes.The dynamics of node i in layer α of the multilayer network is governed by the following set of differentialequations, ˙ x αi = F α ( x αi ) + σ αα N α (cid:88) j =1 A ααij H αα ( x αj ) (cid:124) (cid:123)(cid:122) (cid:125) intra-layer couplings + (cid:88) β (cid:54) = α σ αβ N β (cid:88) j =1 A αβij H αβ ( x βj ) (cid:124) (cid:123)(cid:122) (cid:125) inter-layer couplings , (1) i = 1 , ..., N α , α = 1 , ..., M . The underlying assumption here is that functions that are labeled differently aredifferent functions. By this we mean F α ( x ) (cid:54) = F β ( x ), for β (cid:54) = α .Next we introduce a vectorial notation. We combine the dynamical variables from each layer into a vectorof vectors x α = [ x α , x α , ..., x αN α ] (here and after, the comma stands for vertical stacking of vector and thesquare brackets stand for vector concatenation). Likewise, we can define, F α ( x α ) = [ F α ( x α ) , F α ( x α ) , ..., F α ( x αN α )] , (2) H αα ( x α ) = [ H αα ( x α ) , H αα ( x α ) , ..., H αα ( x αN α )] , (3) H αβ ( x β ) = [ H αβ ( x β ) , H αβ ( x β ) , ..., H αβ ( x βN β )] . (4)This notation suppresses the summations over the nodes in each layer. We obtain the set (one for each layer)of vector equations˙ x α = F α ( x α ) + σ αα A αα H αα ( x α ) + (cid:88) β (cid:54) = α σ αβ A αβ H αβ ( x β ) , α = 1 , ..., M. (5)Two nodes i and j are synchronized if x i ( t ) = x j ( t ) ∀ t . Synchronous motions define invariant manifoldsof Eqs. (5), i.e., x i = x j ⇒ ˙ x i = ˙ x j . In order to study cluster synchronization, we will look for symmetriesin the multilayer network; these are the node permutations that leave system (5) unchanged. The group ofsymmetries of the multilayer network is the set of permutations of the nodes that do not change Eqs. (5)for α = 1 , ..., M . Below we explain how to find this group of symmetries in arbitrary multilayer networks,and we investigate the relations between intralayer and interlayer connections in order to preserve certainsymmetries properties. Symmetries of Multilayer Networks
Symmetries of complex networks have received recent attention from the scientific community. These sym-metries influence network structural [28], spectral [38], and dynamical properties [39], including clustersynchronization [20, 21]. However, symmetries of multilayer networks have not been studied.3 αβ A αα = 0 1 0 111 0 1 010 1 0 111 0 1 011 1 1 10 , A αβ = 1 0 01 0 00 0 10 0 10 0 0 = ( A βα ) T A ββ = 0 1 01 0 10 1 0 αβ αβ ba dc Figure 1:
Illustration of Independent Layer Symmetries and Dependent Layer Symmetries.a
A simple two layer network with intra-connectivity described by the matrices A αα and A ββ and inter-connectivity described by the matrix A αβ = ( A βα ) T (in b ). c Example of an Independent Layer Symmetryfor the multilayer network in a . d Example of a Dependent Layer Symmetry for the multilayer network in a . Here we define the group of symmetries of general multilayer networks with both nodes of different typesand connections of different types. Although we will eventually use computational methods to find thesymmetry group of a particular network, it is best to understand what types of structures are imposed onthe symmetry group of a multilayered network in general. The study of the symmetry group provides insightinto the results of the numerical calculations and into whether what we find is general or just a property ofthe particular network we are analyzing. Next we develop some of the mathematical properties of the finalpermutation group structure of the multilayer systems. As we will see, the interplay between the symmetriesin each layer with the symmetries in other layers is subtle and leads to interesting structures in the finalgroup of the whole network.We first show that the multilayer structure imposes a block diagonal form on the permutations of thewhole group. Since each layer of the network contains a different kind of node, symmetries are not allowed tomove nodes from different layers. As a result, the symmetry group G of the multilayer network is represented4y block diagonal permutation matrices, i. e., each g ∈ G is of the form, g = g α . . . g β . . . g γ . . . ... ... ... . . . , (6)where g α is a permutation matrix for layer α , g β is a permutation matrix for layer β , and so on. In thissection we show that the symmetries of the group G of the multilayer network are formed of the symmetriesof the single-layers (i.e., g α ∈ G α , g β ∈ G β , ...) that are compatible with each other, where compatibility isdetermined by the inter-layer couplings. This means that if we perform a permutation on layer α , we mustperform a permutation on each other layer (which in some cases may be the identity permutation) by leavingthe structure of the overall network unaltered. As a result, compatibility links permutations from the groupsof symmetries of each layer, so that we preserve the symmetry of the whole multilayer network. The majorquestion we want to answer here is: how does compatibility relate permutations from different layers?Let us consider a multilayer system with two layers α and β . Eq. (5) becomes,˙ x α = F α ( x α ) + σ αα A αα H αα ( x α ) + σ αβ A αβ H αβ ( x β )˙ x β = F β ( x β ) + σ ββ A ββ H ββ ( x β ) + σ βα A βα H βα ( x α ) . (7)In Methods, we show that we can determine the symmetry group of a multilayer network with M > g ∈ G α and h ∈ G β .We know from Eq. (6) that a symmetry of this multilayer network is in the form (cid:18) g h (cid:19) , (8)so, which g ’s and h ’s are compatible? For the answer, we consider the inter-layer couplings terms, σ αβ A αβ H αβ ( xxx β ) and σ βα A βα H βα ( xxx α ) . (9)We say that symmetry-related nodes must be flow invariant. That is, the symmetries must guarantee thatsynchronized nodes have equal dynamical variables when we include inter-layer coupling.Application of the two permutations to the equations of the multilayer network results in g ( ˙ x α ) = F α ( g x α ) + σ αα A αα H αα ( g x α ) + σ αβ gA αβ H αβ ( x β ) h ( ˙ x β ) = F β ( h x β ) + σ ββ A ββ H ββ ( h x β ) + σ βα hA βα H βα ( x α ) , (10)where the g and h permutations can be taken into the arguments of F α , H αα and F β , H ββ , respectivelysince the functions operate sequentially on their vector arguments and we have used the properties that g commutes with A αα and h commutes with A ββ .To achieve flow invariance, we need g and h to act on x α and x β , respectively, in the arguments of theinterlayer coupling terms H αβ and H βα . Note that this does not require commutability ( gA αβ = A αβ g ) sincewe want to ‘exchange’ g for h and vice-versa (as the interlayer coupling matrices are generally not square,we do not expect commmutability). More generally, conjugacy relations are the requirements for symmetrycompatibility, i.e., gA αβ = A αβ h and hA βα = A βα g ; (11)which can be thought of as compatibility relations between permutations of the α and β layers. This meansthat the g ’s and the h ’s must be paired properly to satisfy Eq. (11). In general, not all the h ’s are compatiblewith all the g ’s, as the conjugacy relations restrict compatible permutations to subgroups of G α and G β . Thefinal group G of the multilayer network is determined by the structure of these particular subgroups.The permutations that fulfill the conjugacy relations (11) are defined by the following sets,5 α = { g ∈ G α | gA αβ = A αβ h and hA βα = A βα g for some h ∈ G β } (12)and H β = { h ∈ G β | hA βα = A βα g and gA αβ = A αβ h for some g ∈ G α } . (13)An element g ∈ H α is represented by an N α × N α matrix and an element h ∈ H β is represented by an N β × N β matrix. It is important to note that in general the matrix A αβ ∈ R N α × N β will have nontrivial leftand right null spaces. As a result, we may find more than one g that satisfies gA αβ = A αβ h for a given h and vice versa.In the Methods (Formal Proofs) we first prove that H α is a subgroup of G α and H β is a subgroup of G β .Then we show that for the case of undirected networks we only need either one of the following relationshipsto prove that H α and H β are subgroups, gA αβ = A αβ h or hA βα = A βα g. (14)We can now find the group of symmetries of the multilayer network G from the subgroups H α and H β .First, recall that the permutations of G have block structure of Eq. (8). While we use all permutations g ∈ H α and h ∈ H β to construct G , we can only pair the permutations that satisfy the conjugacy relation inEq. (11) (or the simpler version in Eq. (14)).We define an equivalence relation ∼ between the elements of H α : g ∼ g (cid:48) if gA αβ = g (cid:48) A αβ . Analogously, h ∼ h (cid:48) if hA βα = h (cid:48) A βα . One can verify that ∼ is an equivalence relation as it is reflexive, symmetric, andtransitive. We also see that if g ∼ g (cid:48) and g and h are conjugate, then so are g (cid:48) and h , indicating that ∼ defines a partition of each subgroup H α and H β into disjoint subsets (called equivalence classes) K αi and K βi , i = 1 , ..., K , respectively. Each subset K αi contains all the permutations g such that gA αβ is equal to agiven matrix M i ; correspondingly, each subset K βi contains all the permutations h such that hA βα is equalto the matrix M Ti . We can then construct the group of symmetries of the multilayer network G as follows, G = (cid:26)(cid:18) g h (cid:19) | g ∈ K αi and h ∈ K βi , for i = 1 , ..., K (cid:27) (15)Although the sets K αi and K βi have a one-to-one correspondence, the elements of each do not and often differin number, as shown in Methods for the example multilayer network in Fig. 1.Next, we present properties of the equivalence classes K αi , which then leads to the definition of Indepen-dent Layer Symmetries (ILS). Let K α be the equivalence class containing the identity. We can prove that K α is a normal subgroup of H α . Moreover, we can prove that all the K αi , i (cid:54) = 1, are left and right cosets of K α . Formal proofs for each statement are included in Methods (Formal Proofs.)The structure of the equivalence classes K αi gives insight into the symmetries; this facilitates the calcu-lations of the classes themselves and of the T matrices for each layer (see Methods for the structure of the T matrix). In particular the subgroups K α , K β , ... identify a special set of symmetries we will refer to asIndependent Layer Symmetries. The following general relationship between the symmetry operations from K α , K β ,... and the layer structures clarifies why we refer to the elements of K α , K β ,... as Independent LayerSymmetries. Given any g ∈ K α , then for each h ∈ K β , gA αβ = A αβ h = A αβ , (16)since the equivalence class subgroups are both associated with the identity operation. This means thatthe symmetry operations from the equivalence class subgroups do not affect the interlayer coupling and,hence, those operations (permutations) in one layer do not affect the types of allowed dynamics in the otherlayer(s). In particular, the loss of a symmetry associated with the ILS subgroup in the dynamics, i.e. asymmetry-breaking bifurcation, will not alter the possible clusters in other layers, provided no symmetriesin the other cosets ( K αi , i (cid:54) = 1) are also broken. We will show simple examples of this in the Methods Section6able 1: Symmetries of Real Multilayer Networks. For each network dataset, we include information onthe number of layers M , the number of nodes in each layer N α , α = 1 , ...M (the number of nodes N in alllayers for multiplex networks), the number of edges E , the order of the automorphism group |G| , and thecardinality of the largest orbital cluster max( |C k | ).Name M Number of nodes E |G| max( |C k | ) London Transport ** [40] N = 369 441 4 3 EU-AIR Transport [41] N = 450 3588 1475532 × ARXIV NETSCIENCE ** [42] N = 14489 59026 8349492 × PIERRE AUGER Coauthorship ** [42] N = 514 7153 19726 × CKM PHYSICIANS Social ∗ [43] N = 246 1551 120 2 BOS Genetic ∗ [44, 45] N = 321 325 2 × CANDIDA Genetic ∗ [44, 45] N = 367 397 2 × DANIORERIO Genetic ∗ [44, 45] N = 155 188 798 × US power-grid [46] N α = 4492; N β = 449 6994 5 × * Network treated as undirected. ** Weighted Network but treated as unweighted. and when presenting the experiment. Note that the stability of the clusters is a different issue, which willreceive separate consideration.To conclude, the ILS clusters are determined by the interplay between the intralayer couplings in eachlayer and the interlayer couplings between layers. In what follows we will refer to symmetries that are notILS, i.e. they involve simultaneously swapping nodes in different layers as Dependent Layer Symmetries(DLS).From knowledge of the group of symmetries of the multilayer network, the nodes in each layer α can bepartitioned into L α orbital clusters, C α , C α , ..., C αL α . Inside each layer, symmetries map into each other onlynodes in the same orbital cluster.In Table 1 we apply the techniques described in this section to compute the symmetries of severalmultilayer networks from datasets in the literature. For each network dataset, we include information onthe number of layers M , the number of nodes in each layer N α , α = 1 , ...M (the number of nodes N inall layers for multiplex networks), the number of edges E , the order of the automorphism group |G| , andthe cardinality of the largest orbital cluster max( |C k | ). The main underlying assumption is that nodes arehomogeneous inside each layer but nodes in different layers are not, so that nodes inside each layer canbe symmetric, but nodes from different layers cannot. While we understand this is an oversimplification,available datasets do not typically include information on the attributes of the individual nodes, and as aresult, our symmetry analysis is solely based on the multilayer network structure [37] and not on the specificnode attributes. An extensive description of each one of the datasets is included in the Supplementary Note1 and in the Supplementary Table 1. From Table 1 we see that several real multilayer networks possess verylarge numbers of symmetries.In Supplementary Note 2 we study the emergence of symmetries in artificially generated multilayernetworks, where each layer is a Scale Free (SF) network and nodes from different layers are randomlymatched to one another to obtain a multiplex network. We see that these networks typically do not displaysymmetries, except for the case that the power-law degree distribution exponents of the networks in each layerare low across the layers. It is interesting that some of the real multilayer networks we have analyzed displaymany more symmetries than these model multilayer networks. This observation motivated us to designa generating algorithm to construct multilayer networks with prescribed number of symmetries, which ispresented in Supplementary Note 3. 7 tability Analysis In what follows we study stability of the cluster-synchronous solution for a general multilayer networkdescribed by Eqs. (1). Given an orbital partition, we define the L α × L α intralayer quotient matrix Q αα such that for each pair of α -clusters ( C αu , C αv ), Q ααuv = (cid:88) j ∈C αv A ααij , i ∈ C αu , u, v = 1 , , . . . L α . (17)Analogously, we define the L α × L β interlayer quotient matrix Q αβ such that for each pair of clusters, thefirst cluster C αu from layer α and the second cluster C βv from layer β , Q αβuv = (cid:88) j ∈C βv A αβij , i ∈ C αu , u = 1 , , . . . L α , v = 1 , , . . . L β . (18)We can thus write the equations for the time evolution of the quotient multilayer network,˙ q αu = F α ( q αu ) + σ αα L α (cid:88) v =1 Q ααuv H αα ( q αv ) + (cid:88) β (cid:54) = α σ αβ L β (cid:88) v =1 Q αβuv H αβ ( q βv ) , (19) α = 1 , ..., M , u = 1 , ..., L α . Note that Eqs. (19) provides a mapping for each layer α from the nodecoordinates { x αi } , i = 1 , ..., N α to the quotient coordinates { q αu } , u = 1 , ..., L α , where x αi ( t ) ≡ q αu ( t ) if i ∈ C αu .We now linearize (1) about (19), δ ˙ x αi = DF α ( q αu ) δ x αi + (cid:88) β σ αβ N β (cid:88) j =1 A αβij DH αβ ( q βv ) δ x βj , i = 1 , ..., N α , (20)where again q βv is the time evolution of the quotient network node v that node j ∈ C βv maps to. Alsonote that in the above equation the summation in β runs over both intra-layer connections and inter-layerconnections.For each layer α , the above set of equations, can be written in vectorial form, by stacking together allthe individual perturbations applied to the vectors inside each layer, e.g., δ x α = [ δ x α , δ x α , ..., δ x αN α ], δ ˙ x α = (cid:32) L α (cid:88) u =1 E αu ⊗ DF α ( q αu ) (cid:33) δ x α + (cid:88) β σ αβ ( A αβ ⊗ I n α ) L β (cid:88) u =1 (cid:0) E βu ⊗ DH αβ ( q βu ) (cid:1) δ x β , (21) α = 1 , ..., M, where each indicator matrix E αu has dimension N α and is such that its diagonal entries areequal to one if node i of layer α is in cluster C αu and are equal to zero otherwise.We then stack all the layers one above the other to form the vector δ x = [ δ x α , δ x β , . . . ], and we rewrite(21) as δ ˙ x = (cid:80) u E αu ⊗ DF α ( q αu ) 0 · · · (cid:80) u E βu ⊗ DF β ( q βu ) · · · ... ... . . . + (cid:80) u A αα E αu ⊗ D ˆ H αα ( q αu ) (cid:80) u A αβ E βu ⊗ D ˆ H αβ ( q βu ) · · · (cid:80) u A βα E αu ⊗ D ˆ H βα ( q βu ) (cid:80) u A ββ E βu ⊗ D ˆ H ββ ( q βu ) · · · ... ... . . . δ x , (22)where D ˆ H αβ = σ αβ DH αβ . 8e are looking for a transformation that applied to Eq. (22), leaves the first terms on the right handside of (22) unchanged and decouples the second terms in independent blocks, independent of the DF andthe DH terms as they vary in time.From knowledge of the group of symmetries of the multilayer network G , we can compute the irreduciblerepresentations (IRRs) of G . For each layer α we can define an orthonormal transformation T α to the IRRcoordinate system (see [20]). We can then construct the following block diagonal orthonormal matrix, T = (cid:77) α T α , (23)that maps the entire multilayer network to the IRR coordinate system.A formal proof of the above particular structure of the matrix T can be found in Methods (Propertiesof the Equivalence Classes and structure of the matrix T ). Intuitively, the matrix T has a block diagonalstructure, where each block corresponds to a layer, because only nodes from the same layer can be swappedby a symmetry.Consider now the N -dimensional supra-adjacency matrix, A = A αα A αβ · · · A βα A ββ · · · ... ... . . . . (24)The transformed N × N block diagonal matrix B = T AT − is a direct sum ⊕ Rr =1 I d r ⊗ ˆ B r , where ˆ B r is a(generally complex) p r × p r matrix with p r the multiplicity of the r th IRR in the permutation representation, R the number of IRRs present and d r the dimension of the r th IRR, so that (cid:80) r d r p r = N . The matrix T contains information on which perturbations affecting different clusters get mapped to different IRRs[24]. There is one representation (labeled r = 1) which we call trivial and has dimension d = (cid:80) α L α . Allperturbations parallel to the synchronization manifold get mapped to this representation. Hence, the trivialrepresentation is associated with all the clusters C , ..., C L , C , ..., C L , ... However, it is possible that otherIRR representations are only associated with some of the clusters (not all of them).We can now define the (cid:80) α N α n α -dimensional orthonormal matrix ˜ T = (cid:76) α T α ⊗ I n α . Next, we will usethe matrix ˜ T to block-diagonalize Eq. (21). Applying the transformation η = ˜ T δ x to (22) we obtain,˙ η = (cid:80) u J αu ⊗ DF α ( q αu ) 0 · · · (cid:80) u J βu ⊗ DF β ( q βu ) · · · ... ... . . . + (cid:80) u B αα J αu ⊗ D ˆ H αα ( q αu ) (cid:80) u B αβ J βu ⊗ D ˆ H αβ ( q βu ) · · · (cid:80) u B βα J αu ⊗ D ˆ H βα ( q βu ) (cid:80) u B ββ J βu ⊗ D ˆ H ββ ( q βu ) · · · ... ... . . . η , (25)where each transformed indicator matrix J αu = T α E αu T αT and each block B αβ = T α A αβ T βT . The advantage of (25) over (22) lies in the block-diagonal structure of the matrix B = ⊕ Rr =1 I d r ⊗ ˆ B r . Theblock ˆ B is associated with motion along the synchronization manifold. The blocks ˆ B , ..., ˆ B R describe thedynamics transverse to the synchronization manifold. As a result, we have decoupled the dynamics alongthe synchronous manifold from that transverse to it [20]. Moreover, each transverse block r = 2 , .., R isassociated with either an individual cluster or a subset of intertwined clusters [20]. Thus the problem ofstudying the behavior of a perturbation away from the synchronous solution is typically reduced into manysmaller problems, which can be analyzed independently one from the other.Next we discuss how Dependent and Independent Layer Symmetries affect the stability analysis. Werecall that the matrix B = T AT − can be written as a direct sum of blocks ⊕ Rr =1 I d r ⊗ ˆ B r . As each rowof the matrix T is associated to a specific cluster [47], each one of the blocks ˆ B r corresponds to a set of9lusters which are identified by the rows of the matrix T . The trivial representation ( r = 1) is associatedwith all the clusters C , ..., C L , C , ..., C L , ... The corresponding rows of the matrix T are the eigenvectorsassociated to the eigenvalue 1 of the trivial representation; these have nonzero components all of the samesign. Now we look at the remaining ‘transverse’ blocks of the matrix B ( r > T are such that the sums of their entries is equal zero and are called symmetry breakings:each row, in fact, describes how a cluster, generated by a symmetry, may break into smaller ones. Thetransverse blocks can be divided into ILS blocks if the corresponding symmetry breakings are all from thesame layer and DLS blocks if the corresponding symmetry breakings are from different layers. We can thenwrite T = [ T T SYNC , T T ILS , T T DLS ] T and the transformed vector η = [ η T SYNC , η T ILS , η T DLS ] T .The ILS blocks are symmetry breakings generated by the ILS subgroup. From our definition of the ILSsubgroup, each K α is a normal subgroup of H α . Clifford theorem [48] states that each IRR of H α , whenrestricted on K α , is either itself an IRR of K α or breaks up into a direct sum of IRRs of K α of the samedimension. The rows of the change of coordinates T α ILS (for layer α ) to the IRR of K α , which are associatedwith the symmetry breaking perturbations of the ILS, are generated by the eigenvectors associated to theeigenvalue 1 of the projectors on the IRRs of K α . These rows generate invariant subspaces of minimaldimension, and must therefore be rows of T α .We now consider for simplicity the case of a network with two layers, labeled α and β . In general, thetransverse IRRs associated with the ILS subgroup in layer α will have structure,˙ η ILS = (cid:104)(cid:88) u J αu ⊗ DF α ( q αu ) + (cid:88) u B αα J αu ⊗ D ˆ H αα ( q αu ) (cid:105) η ILS . (26)Note the perturbation η ILS is independent of the dynamics on the β layer, i.e., of both DF β ( q βu ) and D ˆ H ββ ( q βu ). This indicates that the stability of ILS symmetries can be studied through a specific class ofthe Master Stability Function, which is the same as for single-layer networks [20]. The transverse IRRsassociated with the ILS subgroup in layer β will have an analogous structure as (26). On the other hand,the remaining transverse IRRs will have structure,˙ η DLS = (cid:18)(cid:20)(cid:80) u J αu ⊗ DF α ( q αu ) 00 (cid:80) u J βu ⊗ DF β ( q βu ) (cid:21) + (cid:20)(cid:80) u B αα J αu ⊗ D ˆ H αα ( q αu ) (cid:80) u B αβ J βu ⊗ D ˆ H αβ ( q βu ) (cid:80) u B βα J αu ⊗ D ˆ H βα ( q βu ) (cid:80) u B ββ J βu ⊗ D ˆ H ββ ( q βu ) (cid:21)(cid:19) η DLS . (27)Note that the perturbations η DLS depend on the dynamics of the systems in both layers α and β , throughthe mixed blocks appearing in Eq. (27). The structure of Eq. (27) is substantially different than that of Eq.(26), which may lead to dramatic effects in terms of stability.Consider now the multilalyer network in Fig. 1. The matrix T is equal to T = . . . . √ √ . − . . − . . − . − . . . . − . − . √ − √ SYNC (cid:27)
ILS breakings (cid:27)
DLS breakings (28)10he matrix B = T AT T is equal to: B = √ √ √ √ − √
20 0 0 0 0 0 √ SYNC block (cid:27)
ILS blocks (cid:27)
DLS block (29)As can be seen, the matrix B is the direct sum of one SYNC block, two independent ILS blocks (correspondingto breakings in the α layer which are independent of the β layer), and one DLS block (corresponding tobreakings in the α layer and in the β layer which are dependent on each other.) The transverse ILS andDLS blocks describe linear (local) stability of clusters, when the dynamics is linearized about the quotientnetwork time evolution. In what follows we will pay particular attention to stability of DLS blocks, whichare a specific feature of multilayer networks.To provide analytical insight, we first consider a two-layer network described by the following discrete-time dynamics, x αk +1 = rem( c α x αk + σA αα x αk + σA αβ x βk ) x βk +1 = rem( c β x βk + σA ββ x βk + σA βα x αk ) , (30)where the vectorial function rem( x ) returns a vector whose entries are the remainder of the integer divisionof the entries of the vector x by the scalar 1 and c α and c β are tunable layer-specific scalar parameters.Stability is described by the following set of equations, δ x αk +1 = c α δ x αk + σA αα δ x αk + σA αβ δ x βk δ x βk +1 = c β δ x βk + σA ββ δ x βk + σA βα δ x αk . (31)We now consider a generic DLS block of the form (cid:18) a bb a (cid:19) ( a = 0 and b = √ η k +1DLS = (cid:18) c α + σa σbσb c β + σa (cid:19) η k DLS , (32)with c α (cid:54) = c β and the case c α = c β corresponding to an ILS perturbation. The eigenvalues have the followingexpression, ρ ± = ¯ c + aσ ± ( b σ + δ c ) / , (33)where ¯ c = ( c β + c α ) / δ c = ( c β − c α ) / δ c = 0 corresponding to the case of an ILS.) From this equation wesee the larger (smaller) eigenvalue is an increasing (decreasing) function of δ c , and it follows that the bestcondition in terms of stability is achieved for δ c = 0. We thus conclude that DLS perturbations are moredifficult to stabilize than ILS perturbations.We present here a conjecture expanding on the above conclusion: that DLS clusters are generally moredifficult to stabilize compared to the case in which the systems in different layers are of the same type. Webase this on the following reasoning. Consider an ideal situation in which at first the parameters of thesystems in the different layers are identical and then they are increasingly perturbed to take on differentvalues in different layers. The Gershgorin circle’s theorem states that each eigenvalue of a matrix lies withinat least one of the Gershgorin discs centered at the entries on the main diagonal and having radius equalto the sum of the off diagonal entries. As perturbing the individual system’s parameters corresponds to11arying the centers of the Gershgorin discs (but not the radius), we can expect the eigenvalues to becomeincreasingly spread out as the systems are made increasingly different from one another, resulting in reducedstability.The particular network in Fig. 1 has one transverse irreducible representation (DLS) in the form of Eq.(32), with a = 0 and b = √
2, corresponding to simultaneous loss of synchronization between nodes (1 ,
3) inthe β layer and between the pairs of nodes (1 ,
2) and (3 ,
4) in the α layer. For this multilayer network, weconsider the dynamics of Eq. (30) and study the effects of changing the parameters c α and c β on stability.For the maps described by Eq. (30), as well as in the following for other kind of systems, we measure thepairwise synchronization error, defined as, E αij = < (cid:107) x αi − x αj (cid:107) >, (34)between nodes i and j in layer α = 1 , ..., M , where the symbol < ... > indicates a temporal average, computedafter the transient has elapsed.Panel a of Fig. 2 shows E α and E β vs the parameter δ c . δ c = 0 corresponds to identical systemsin the two layers, increasing values of δ c indicate the systems in the two layers are increasingly different.Synchronization is simultaneously lost in both layers for δ c (cid:39) .
35. Panel b is a plot of the eigenvalues ρ + and ρ − in Eq. (33) as a function of δ c , which shows that ρ + ( ρ − ) increases (decreases) with δ c . Loss ofstability occurs when either one of the two eigenvalues is either larger than 1 or smaller than −
1. It is thusexpected that stability may be lost for increasing values of δ c , i.e., as the systems in the two layers becomeincreasingly different. From Panel b of Fig. 2 we see that ρ + grows larger than 1 for δ c (cid:39) . M = 2 layers, A αα , A ββ , A αβ = ( A βα ) T corresponding to themultilayer network in Fig. 1, n α = n β = 2, x α = [ x α , y α ], x β = [ x β , y β ], F α ( x α ) = (cid:20) y α − x α − . y α [( x α ) − c α ] (cid:21) , F β ( x β ) = (cid:20) y β − x β − . y β [( x β ) − c β ] (cid:21) , (35)which both correspond to the dynamics of the Van der Pol oscillator, with layer-dependent parameters c α and c β . Moreover, we set H αβ ( x ) = (cid:18) − (cid:19) x , (36)for all pairs α, β = 1 , σ αβ = 0 .
15 for all pairs α, β = 1 ,
2. We set the layer-dependent parameters c α = (1 + δ c ) and c β = (1 − δ c ), so that δ c = 0 corresponds to identical systems in the two layers andincreasing values of δ c indicate the systems in the two layers are increasingly different. We then numericallyinvestigate stability of the DLS as the parameter δ c is increased. This can be seen in panel c of Fig. 2,which shows that synchronization between oscillators 1 and 3 from the α layer and synchronization betweenoscillators 1 and 3 from the β layer is simultaneously lost for δ c (cid:39) .
4. Panel d of Fig. 2 shows the numericallycomputed Maximum Lyapunov exponent of the DLS block (27). We also note that throughout the whole δ c interval considered, neither nodes 1 and 2, nor nodes 1 and 4 from the α layer ever synchronize (not shown).Overall, Fig. 2 shows that increased heterogeneity between the nodes in the two layers, can lead to loss ofDLS stability. Further numerical evidence of this is presented for the cases of the Lorenz oscillator and of theRoessler oscillator in Supplementary Note 4. Supplementary Note 5 investigates how varying the intra-layerand inter-layer coupling strengths affects CS stability. An experimental testbed for cluster synchronization
We apply the techniques from the previous sections on an experimental testbed circuit. We implement threedifferent kinds of electronic oscillators, using widely available, affordable components that can be assembledon breadboards. Simplicity, low-cost, ease of fabrication and availability of a large volume of previous studiesmake electronic circuits an ideal testbed for multilayer network studies [31]. The circuit is composed of threedifferent kinds/layers of nodes: one linear resonator, which we call the “jumper”, two FitzHugh-Nagumo(FHN) oscillators, and four Colpitts oscillators. The jumper is a linear resonator (i.e., without input, its12 . . . . . . . . δ c E α ( ) , E β ( ) c . . . . . δ c E α ( ) , E β () b . . . . . − δ c ρ + ( - ) , ρ − ( - ) d . . . . . − . − . − . − . . δ c M L E Figure 2:
Synchronization of DLS clusters. a
Discrete time maps. Synchronization errors E α (asterisks)and E β (diamonds) vs the parameter δ c . δ c = 0 corresponds to identical systems between the α and in the β layers. Increasing values of δ c are for increasingly different systems in the two layers. Synchronization islost in both layers for δ c (cid:39) . b Discrete time maps. Eigenvalues ρ + and ρ − in Eq. (33). As can beseen, ρ + ( ρ − ) increases (decreases) with δ c . Loss of stability occurs when either one of the two eigenvalues iseither larger than 1 or smaller than − ρ + becomes larger than 1 for δ c (cid:39) . c Van Der Pol oscillators.Synchronization errors E α (asterisks) and E β (diamonds) vs the parameter δ c . δ c = 0 corresponds toidentical systems between the α and in the β layers. Increasing values of δ c are for increasingly differentsystems in the two layers. Synchronization is lost in both layers for δ c (cid:39) . d Van Der Pol oscillators.Maximum Lyapunov Exponent of the DLS block (27) vs the parameter δ c .13scillations damp to zero) with two-dimensional governing equations; the FHN is a relaxational nonlinearoscillator with two-dimensional governing equations [49]; the Colpitts is a sinusoidal nonlinear oscillator withthree-dimensional governing equations [50]. All three kinds of nodes have similar uncoupled frequencies. Afull schematic of the experimental circuit is included in Figure 3. Figure 3 also states the measured parametervalues of each component.Figure 4 shows a simplified multilayer schematic of the circuit with the jumper as the α layer, the FHNsas the β layer, and the Colpitts as the γ layer. We couple the different oscillators in three ways. On theColpitts layer, we induce inductive coupling through mutual inductance, red in Fig. 4. Between the Colpittslayer and the FHN layer we introduce resistive coupling; this couples each FHN circuit with two Colpitts,blue in Fig. 4. Between the FHN and jumper layer, we induce inductive coupling between each FHN and aninductor in the jumper, orange in Fig. 4. FHN 1 Jumper FHN 2Colpitts 2Colpitts 1 Colpitts 4Colpitts 3 R R R R C F,1 L FJ + - L F R R R R C F,2 L FJ + - L F L J1 V L J2 R J C J x x R L,1 C ce,1 C be,1 R ee,1 V be,1 V ce,1 L C V cc V ee R L,2 R C R C C ce,2 C be,2 R ee,2 V be,2 V ce,2 L C V cc V ee x R L,3 C ce,3 C be,3 R ee,3 V be,3 V ce,3 L C V cc V ee R L,4 R C R C C ce,4 C be,4 R ee,4 V be,4 V ce,4 L C V cc V ee x -1.6 V+3.2 V25.3Ω27.4Ω 75.2Ω75.0 Ω 99.2 Ω 10.0 Ω10.2 nF 99.7 Ω115.9 Ω 1.99 kΩ1.98 kΩ 2.00 kΩ1.98 kΩ 117.0 Ω60.9 nF58.9 nF60.9 nF 1.17 nF 1.15 nF60.0 nF95 μH 95 μH 95 μH 3.0 V 95 μH150 μH 95 μH150 μH 27.6Ω25.5Ω74.9Ω75.7Ω59.7 nF60.5 nF61.4 nF60.2 nF 95 μH95 μH95 μH Transistor 2N2222 Figure 3:
Schematic of the experimental setup. R C = 2 . k Ω . We vary the magnetic coupling betweenColpitts 1 and 2 and Colpitts 3 and 4, k C , from − . − . x (red). We hold theseparation between the FHNs and the jumper (x in orange) constant such that the coefficient of mutualinductance, k F J , is maintained equal to 0.35.As shown in the above schematic and equations, we can describe the system with three layers composedas follows. Layer α includes the jumper alone, layer β includes 2 FHNs and the intra-layer connections, with G β = { () , (1 , } and layer γ includes 4 Colpitts and one intra-layer connection that connects nodes 1 with2, and nodes 3 with 4. The group of symmetries of this layer is G γ = { () , (1 , , (3 , , (1 , , , (1 , , , (1 , , } . (37)We have the inter-layer connections: A αβ = (cid:2) (cid:3) = ( A βα ) T , A βγ = (cid:20) (cid:21) = ( A γβ ) T . (38)14 F C C C C J R c k c k FJ α β γ Figure 4:
Multilayer representation of the experimental circuit.
Layer α contains the jumper node.Layer β contains two FHN oscillators. Layer γ contains four Colpitts oscillators with inductive intralayercoupling shown in red. The interlayer resistive coupling between the β and γ layers is shown in blue. Theinterlayer inductive coupling between the α and β layers is shown in orange; this coupling introduces avirtual FHN intralayer coupling which we discuss in the analysis of the system.All the permutations in the groups G α , G β , G γ are compatible, then H α = G α , H β = G β , H γ = G γ . The interlayer connections define the conjugate classes K α = (); K β = () , K β = (1 , K γ = { () , (1 , , (3 , , (1 , , } , K γ = { (1 , , , (1 , , } . (39)This set defines the group of symmetry of the multilayer network. There are three possible clustered patterns:(1a) layer β and layer γ fully synchronized, (1b) layer β fully synchronized and layer γ clustered synchronized(either 1 with 3 and 2 with 4 or, equivalently, 1 with 4 and 2 with 3), and (2) layer β not synchronizedand layer γ synchronized in clusters (1 with 2 and 3 with 4). The stability of each one of these clusteredpatterns is analyzed in detail in the Supplementary Note 6 and is here summarized in Figure 5. Plot a shows the two MLEs transverse to solution 1a. Solid and dotted lines refer to transverse blocks B and B b ,respectively (blocks defined in the Supplementary Note 6). Pattern 1a is stable when both the curves arenegative. b shows MLEs transverse to stable solutions on the synchronized manifold 1b. Blue curves referto the transverse MLE of solution of kind 1a, while the red curve refers to the transverse MLE of solution ofkind 1b. c shows MLEs transverse to stable solutions on the synchronized manifold 2. Blue curves refer tothe transverse MLE of solution of kind 1a, while the yellow curves refer to the transverse MLE of solutionof kind 2.The equivalence class subgroup K γ shows the ILS’s for the γ layer. The symmetry breaking of the(1 , ,
4) permutation for the [1 , , ,
4] cluster causes the cluster to break into two clusters, [1 , , [2 ,
4] or[1 , , [2 , β to remain synchronized in case (1b) above.Note, however, that breaking of the individual permutations like (1,2) by itself will cause a breaking of oneof the K γ permutations, so it is essential to check that the symmetry breakings of K γ are not inducing othersymmetry breakings in the same layer, but from other cosets. If they are not, then the other layers areunaffected by the loss of an ILS.With all three patterns, we can draw Figure 6, where we report the maximum of the two blue curves in15 B B M L E t r a n s v e r s e t o P a tt e r n kC p i t c h f o r k p i t c h f o r k p i t c h f o r k p i t c h f o r k p i t c h f o r k p i t c h f o r k p i t c h f o r k p i t c h f o r k p i t c h f o r k p i t c h f o r k -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4-4-3-2-101 10 M L E t r a n s v e r s e t o P a tt e r n b kC s u p o e r c r i t i c a l p i t c h f o r k s u b c r i t i c a l p i t c h f o r k s u b c r i t i c a l p i t c h f o r k s u p o e r c r i t i c a l p i t c h f o r k s u b c r i t i c a l p i t c h f o r k s u b c r i t i c a l p i t c h f o r k -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4-4-3-2-101 10 M L E t r a n s v e r s e t o P a tt e r n kC s u p o e r c r i t i c a l p i t c h f o r k s u p o e r c r i t i c a l p i t c h f o r k s a dd l e - n o d e s a dd l e - n o d e s u p o e r c r i t i c a l p i t c h f o r k s u p o e r c r i t i c a l p i t c h f o r k s a dd l e - n o d e s a dd l e - n o d e abc Figure 5:
Stability of different cluster synchronization patterns. a
MLEs transverse to solution 1a.Solid and dotted lines refer to different transverse blocks. When a line crosses zero it identifies a symmetrybreaking bifurcation in one of the other invariant manifold (red: bifurcations on CS manifold 1b - yellow:bifurcations on CS manifold 2). Pattern 1a is stable when both the curves are negative. b MLEs transverseto stable solutions on the synchronized manifold 1b. Blue curves refer to the transverse MLE of solution ofkind 1a, while the red curve refers to the transverse MLE of solution of kind 1b. Vertical lines indicate thebifurcations of the quotient system. c MLEs transverse to stable solutions on the synchronized manifold 2.Blue curves refer to the transverse MLE of solution of kind 1a, while the yellow curves refer to the transverseMLE of solution of kind 2. Vertical lines indicate the bifurcations of the quotient system. Colored linesrepresent symmetry breaking bifurcation (inferred from a ), while black lines are other bifurcations.Figure 5 a (blue), the red curve in Figure 5 b (red) and the yellow curves in Figure 5 c (yellow). We can thus16 M L E k c t ( s ) x • , β x • , γ ab β x γ x γ x γ x γ x β x Figure 6:
Onset of different cluster synchronization patterns as varying k C . We observe five kinds ofclustered behavior, which are denoted by the shading of the background. Shading with two-colored stripingindicates bistability between the two states represented by each color. a MLE transverse to each identifiedclustered solution. (Green) Pattern 1a; (Yellow) Pattern 2-( π/ π ); (Blue) Pattern1b; and (White) No pattern. b different stable solutions present in each of the structurally different regions(identified by a bifurcation of a stable solution both in one of the quotient networks or transverse to them).(Star) A stable solution of kind 1b occurs when the red line is negative; (Square) a modified version of kind1b occurs when the red line is positive.identify the parameter regions in which the different analyzed behaviors are present. Below we comment onFigure 6 and explain the sequence of bifurcations that occur in the system as the parameter k C is decreased(i.e., from right to left looking at the figure.) For large positive coupling, 1a is stable and it is the onlyattractor of the system. This region is shaded in green in Figure 6. At k C = 0 .
2, we observe the split of thetwo nodes in layer β (the FHNs, transverse direction toward 2 in the T matrix). The symmetry breaking(pitchfork) bifurcation is supercritical, and a solution of kind 2 is born. This region is shaded in yellow.At k C = 0 .
14 a second stable clustered solution of type 2 is born through a saddle-node bifurcation (viathe quotient network 2 dynamics). This solution is characterized by a phase separation of the two FHNoscillators of π , and a slightly higher oscillations frequency. It remains stable until k C is reduced to -0.28.17e thus have a region with two solutions of kind 2: the one created at k C = 0 . k C = 0 .
14 where the two FHNs are anti-phase (whose MLEs can be found looking at the yellowcurve from left to right). We call these two states pattern 2-( π/
2) and pattern 2-( π ), respectively. The statecreated at k C = 0 .
14 is shaded in lavender; the bistable region is striped lavender and yellow. At k C = 0 . π/
2) loses its stability giving rise to an unsychronous solution, shaded in white.This solution disappears through a saddle-node bifurcation at k C = 0. For k C ∈ [ − . , π ). The second (shaded in blue) has twopossible behaviors: when the red curve is negative, it is of kind 1b (see star for example), or when the redcurve is positive, it is a slightly non synchronized modification of pattern 1b (see square for example) inwhich the two FHNs remains locked, while the Colpitts are not perfectly synchronized. This last solution ispresent because the inductively coupled Colpitts pairs are antiphase; as a result, the FHN (which is coupledto both Colpitts) sees a very small net signal from the two Colpitts. Consequently, the Colpitts can verynearly synchronize even though they receive slightly different inputs. Finally, at k C = − .
28 the pattern 2-( π )solution undergoes a catastrophic symmetry breaking bifurcation, leaving 1b as the only possible attractor ofthe system (or its slight modification when the red curve is positive.) Experimental data showing qualitativeagreement with the numerical results in Figure 6 are included in the Supplemetary Note 7. Discussion
In this paper we study analytically, numerically, and experimentally the general problem of cluster syn-chronization (CS) in multilayer networks; these systems are composed of heterogeneous components thatmay interact in multiple ways. We first present a general set of equations that describes the dynamics of amultilayer network. We then define for the first time the group of symmetries of a multilayer network andexplain how to compute it. An analysis of several datasets of real multilayer networks shows they possesslarge number of symmetries.We investigate the stability of the CS patterns which correspond to the orbits of the multilayer network.The symmetries of each layer as well as the particular pattern of interlayer connections determine the CSpatterns in multilayer networks; we consequently describe how each layer’s symmetry group relates to thesymmetry groups of other layers through interlayer connections. The interplay between the intralayer andinterlayer couplings enables distinction between symmetries that affect the types of allowed dynamics in theother layers (DLS) and those that do not (ILS.) The latter form a normal subgroup of the full symmetrygroup. In particular, a symmetry-breaking bifurcation, associated with the ILS subgroup, breaks up aclustered pattern in one layer of the multilayer network without altering the possible clusters in the otherlayers.With an IRR change of coordinates we decouple the stability problem into several simpler (lower dimen-sional) problems. First, we decouple perturbations along directions parallel to the synchronization manifoldfrom those transverse to it; the latter determine stability of the clustered motion. Second, we decouple theequations for the transverse perturbations into several independent blocks; each block corresponds to thestability of either an individual cluster or a set of intertwined clusters [20]. When two or more clusters(which may belong to different layers) are intertwined, they are either all stable or all unstable. We seethat Dependent Layer Symmetries yield blocks of a different structure than those arising in the study ofsingle layer networks, which we show has a profound effect on the stability of the clusters involving thesesymmetries. In particular, we show analytically for a specific class of networks that DLS clusters are moredifficult to synchronize than in the case in which the systems in different layers are of the same type, which isalso confirmed numerically in simulations involving multilayer networks of Van der Pol, Lorenz, and Roessleroscillators (Supplementary Note 5.)We performe experiments with a fully analog multilayer network with seven electronic oscillators of threedifferent kinds coupled with two kinds of coupling. The testbed circuit has many features that occur innatural multilayer systems, like noise and parameter mismatches. The experiment is a good test of thetheoretical framework because it allow us to understand how theoretical predictions, made with simplifying18ssumptions, can be a guide to better understand real phenomena. In fact, the experimental results largelymatch the theoretical predictions; we observe all the predicted cluster states in the right part of the parameterspace. With respect to the experimental realization, the model includes several simplifications: we assumeall the nodes within the same layer are identical, we assume that all the connections of the same type areidentical, we use a simplified FHN model which neglects two resistors, we use ideal models for the operationalamplifiers and the transistors, and we assume no noise and neglect any stray inductance, capacitance, orresistance. We discuss why the experimental and theoretical results differ in our seven node electronic systemin Supplementary Note 8. A broader, multisystem analysis of the robustness of the method to heterogeneity,noise, and network nonideality would enhance the utility of the method.A variety of real world systems are multilayer networks that can exhibit clustering. Dynamical situations,to which our analysis may be relevant include: opinion dynamics and consensus among individuals interactingthrough different communication systems, for which clustering may show up on average, see e.g., [39], thedynamics of central pattern generators, small networks of similar neurons, which might show symmetriesand clustering since synchronization is part of their dynamics, see e.g., [10], and electronic networks andin general man-made systems formed of many identical subsystems or nodes, which may produce clustersynchronization. Recent work [51] has studied how network symmetries may affect synchronization modes ofpower grids and even suggested that symmetries may enhance complete synchronization in multilayer grids,characterized by the presence of different energy sources, such as power generators, wind turbines, solar, etc.Our work describes how clusters may arise in multilayer networks with a given structure; with further study,it may be possible to infer the structure of a multilayer network given the observation of clusters.A main limitation of our approach is its scalability with the number of nodes of a multilayer network.While the size of the network for which symmetries can be found is very large [52], the size for which thestability analysis can be performed is typically much smaller [23]. We also model the systems in each layeras exactly identical; this is not a characteristic of experimental systems (for a discussion of the discrepanciesbetween our experiment and our simulations see Supplementary Note 8.) A recent paper analyzed clustersynchronization in the presence of nearly identical systems [22], though the approach of [22] is not easilygeneralizable to the case of multilayer networks. Although our multilayer framework does not allow nodeheterogeneity and noise, we saw broad agreement between our predictions and experimental circuit behavior;this circuit contained slightly heterogeneous nodes with noise. Within this study, we cannot quantitativelystate when noise or heterogeneity will qualitatively impact our predictions, but the experiment demonstratesthat there is some tolerance. Further work will be needed to characterize these limitations.Finally, despite the generality of the theory proposed in this manuscript, several extensions are possible.An important direction for the research is to allow both intra-layer and inter-layer connections to be directed.Another direction is to study the formation of CS patterns that are not related to symmetries in multilayernetworks, see e.g. [24]. The case of group consensus which can be seen as a very special case (i.e., with lineardynamics) of the cluster synchronization problem considered here, has recently been studied in [53].
Data availability
All data supporting the findings of this study are available within the article and its supplementary infor-mation files and from the corresponding author upon reasonable request.
Methods
Formal Proofs
Theorem 1. H α is a subgroup of G α and H β is a subgroup of G β .Proof. To prove the theorem we must show that the identity element of G α is contained in H α , and whenever g and g are in H α , then so are h − and h h , so the elements of H α indeed form a group. We prove thetheorem for H α . The proof for H β is identical. 19dentity. Let e α be the identity of G α and e β be the identity of G β . We see that e a lpha is in H a lpha ,since e α A αβ = A αβ = A αβ e β . Inverse existence . For all g ∈ H α exist g − ∈ H α . By definition, if g ∈ H α then there exists h such that gA αβ = A αβ h . Consider the permutation g − , then gA αβ = A αβ h ⇒ g − gA αβ = g − A αβ h ⇒ A αβ h − h = g − A αβ h thus obtaining A αβ = g − A αβ h ⇒ g − A αβ = A αβ h − . Since G β is a group, h ∈ G β ⇒ h − ∈ G β , thusproving that g − ∈ H α . Closure . We need to prove that if g and g are in H α , then the product g g is also in H α . By definition,since g , g ∈ H α , there exist h , h such that g A αβ = A αβ h and g A αβ = A αβ h . Then g g A αβ = g A αβ h = A αβ h h . Since G β is a group, h , h ∈ G β ⇒ h h ∈ G β , thus proving that g g ∈ H α .A consequence of the theorem is that g ∈ H α ⇒ g ∈ G α : thus, to find H α , we have simply to take allthe permutations h ∈ G β and check which of the permutations in G α respect the conjugacy relation (11).Similarly, we can define the conjugate set Eq. (13), find its elements in the same way and show that it is asubgroup of G β .Since we are interested here only in undirected networks, we prove a corollary that simplifies the compu-tations of H α and H β . Corollary 1.
For the case of undirected networks we only need one of the following relationships to provethat H α and H β are subgroups, (i) gA αβ = A αβ h or (ii) hA βα = A βα g .Proof. Let’s choose to use relationship (i) to define H α and H β . This means eliminating requirement (ii)in the H α definition and replacing (ii) with (i) in the H β definition. We can use the same logic we used inthe previous theorem to show that these definitions of H α and H β also produce subgroups. The question is,are they the same as those of Eqs. (12) and (13)? We can show that this is the case using the fact that forundirected coupling the overall coupling matrix is symmetric. Hence, off-diagonal interlayer coupling blocksare related as in A βα = ( A αβ ) T . For permutations h − = h T , etc. Hence, since H α and H β are groups, h − and g − also obey the defining equation (i), so g − A αβ = A αβ h − ⇐⇒ g T A αβ = A αβ h T ⇐⇒ ( A αβ ) T g = h ( A αβ ) T ⇐⇒ ( A αβ ) T g = h ( A αβ ) T ⇐⇒ hA βα = A βα g (40)A consequence of this corollary is that for the case of undirected networks, if h ∈ G β is involved in aconjugacy relation for some g , then it is in H β . As a result, one does not need to search for both H α and H β , but populating H α automatically populates H β . Theorem 2. K α is a subgroup of H α .Proof. Identity . The identity e α is in K α by definition. Inverse . Let g ∈ K α then, e α A αβ = g − gA αβ = A αβ . But we also have gA αβ = A αβ , hence, g − A αβ = A αβ so g − ∈ K α . Closure . Let g and g ∈ K α then, g g A αβ = g A αβ = A αβ = ⇒ g g ∈ K α .Hence, K α is a subgroup of H α .The same reasoning holds for the β layer, etc. Next we state an additional property of the subgroups K α , K β , etc.: 20 orollary 2. K α is a normal subgroup. This means that if g ∈ K α , then ∀ g ∈ H α we have g − g g ∈ K α or, equivalently, g − K α g = K α .Proof. We know ∀ g ∈ H α , if gA αβ = A αβ h , then g − A αβ = A αβ h − . Hence, ∀ g ∈ K α , g − g gA αβ = g − g A αβ h = g − A αβ h = A αβ h − h = A αβ = ⇒ g − g g ∈ K α and K α is a normal subgroup.And, finally, we prove the following corollary. Corollary 3.
All the K αi , i (cid:54) = 1 , are left and right cosets of K α .Proof. If g ∈ K α and g ∈ K αi with gA αβ = A αβ h , then gg A αβ = gA αβ = A αβ h = ⇒ gg ∈ K αi .Similarly we can prove g g ∈ K αi . Since group products are unique, i.e. gg = gg = ⇒ g = g , we have g K α = K αi = K αi g . So , K αi are left and right cosets of K α . This also means all the K αi have the same numberof elements. Properties of the Equivalence Classes and structure of the matrix T Using the equivalence classes of each layer we can show the following is true: the T matrix for the entiremultilayer system is of a block diagonal form: T = (cid:18) T α T β (cid:19) (41)Can we find the matrices of each block independently? This may seem intuitively obvious, but it is notobvious that the arithmetic will work out. Consider the two layer system α and β as above. The group G of the full system is given by the union of direct products such as K αi × K βi . If n α is the number of elementsin K αi and n β is the number of elements in K βi , then the number of elements in an equivalence class directproduct is n α n β . And if K is the number of equivalence classes, then the number of terms in the wholedirect product is Kn α n β .In the calculation of the T matrix for the entire system we form the projection operators P ( l ) for each ofthe IRR labeled by l using the sums [20] P ( l ) = d ( l ) d (cid:88) C µ ( l ) C (cid:88) g ∈C g (42)where C is a conjugacy class, µ ( l ) C is the character of that class for the l th IRR, d ( l ) is the dimension of the l th IRR and d is the order (size) of the group. For the multilayer system d = Kn α n β . The elements of the β level group appear in the sum n α times and the elements of the α level group appear in the sum n β times.This will contribute factors to the sum for each layer. However, because d = Kn α n β is in the denominatorthe extra factors in the numerator will be cancelled in each layer leaving the correct dimension divisor foreach layer’s sum. Hence, we can find the T matrix for the whole system by finding the T matrices ( T α and T β ) for each layer independently and putting them into the block form to construct T . Computing the symmetry group of a simple multilayer network
Let us consider the simple example of a two-layer undirected network shown in Fig. 1. We show in detailthe calculations to determine the final group G of the full network.The group of symmetries of the two layers are (written in cyclic notation) G α = { () , (1 , , , (2 , , (1 , , , , (1 , , , , (1 , , (1 , , , (1 , , } , (43) G β = { () , (1 , } , (44)where, for example, (1,4)(2,3) means move node 1 to node 4 (and 4 to 1) and 3 to 2 (and 2 to 3). Also,(1,4,3,2) means move 1 to 4, 4 to 3, 3 to 2, and 2 to 1.21rom these we determine H α and H β using the results of Theorem 1. To find H α , we take all thepermutations h ∈ G β and we check which of the permutations in G α respect the conjugacy relation (11). Weobtain H α = { () , (1 , , , (1 , , , (1 , , } , (45) H β = { () , (1 , } , (46)Now we must define the ∼ relation. Applying the permutations in H α to A αβ we obtain() = (1 , , = = ⇒ () ∼ (1 , , K α = { () , (1 , , } (1 , , = (1 , , = = ⇒ (1 , , ∼ (1 , , K α = { (1 , , , (1 , , } while applying the permutation of H β to A αβ and recalling the conjugacy relation gA αβ = A αβ h we obtain () = = ⇒ K β = { () } (1 ,
3) = = ⇒ K β = { (1 , } . Combining the permutations from each pair of disjoint subsets as in Eq. (8), we obtain the full group, G = (cid:26)(cid:18) () 00 () (cid:19) , (cid:18) (1 , ,
4) 00 () (cid:19) , (cid:18) (1 , ,
2) 00 (1 , (cid:19) , (cid:18) (1 , ,
4) 00 (1 , (cid:19)(cid:27) . (47)Since K α and K β are the subgroups of H α and H β , respectively, this means that the separate clusters inthe α layer, [1 ,
2] and [3 , β dynamics or synchronous clusters (again, stability is a separate issue):[1 , , , → [1 ,
3] and [2 , , (48)or [1 , , , → [1 ,
4] and [2 , , (49)which would still allow [1,3] to be synchronous in the β layer. However, the bifurcation[1 , , → [1] , [2] , [3] , and [4] (50)would break a symmetry in the coset K α if the nodes 1 and 3 were synchronized. But if they are notsynchronized, then the system is operating in a subgroup of the original group and only the identity would22emain for the β layer. This means in this case 1 and 3 are each in their own (trivial) cluster already inthe β layer and that does not prevent either state [1 , ,
4] or [1] , [2] , [3] , [4] in the α layer. These cases canbe easily seen from Fig. 1 since node 1 in the β layer only depends on the sum of nodes 1 and 2 in the α layer. And similarly for the dependence of node 3 in the β layer on nodes 3 and 4 in the α layer. However,while these ILS relations are easy to see in this simple case, it would require the calculation of K α and K β for more complicated networks to find the ILS’s, associated clusters, and their allowed bifurcations. Symmetries of multiplex and multidimensional networks
Here we describe how to compute the group of symmetries of multiplex and multidimensional networks,which are special cases of the general multilayer network problem presented in the Results Section. We showthat the problems of computing the group of symmetries of a multiplex network and of a multidimensionalnetwork are closely related to each other. In what follows, we start by considering the problem of computingthe symmetry group of a multiplex network and mapping it to the problem of computing the symmetrygroup of a multidimensional network.Multiplex networks [32, 33, 34] are a particular class of multilayer networks; in multiplex networks,different features of the same set of agents are described in each different layer ( e.g., in a social system,each layer represents the opinion of a person on a different topic, and links capture how the different socialinteractions influence the person thinking on each topic). A multiplex network is thus formed of severallayers, with each layer containing the same number of nodes N . Moreover, since interlayer coupling onlyoccurs between node i in a given layer and the same node i in a different layer, A αβ = I N , where I N is the N -dimensional identity matrix.A general multiplex network is governed by the following set of equations [37, 54],˙ x αi = F α ( x αi ) + σ αα N (cid:88) j =1 A ααij H αα ( x αj ) + (cid:88) β (cid:54) = α σ αβ H αβ ( x βi ) , (51) i = 1 , ..., N , where the matrix A αα represents the intra-layer connectivity, the function F α represents theintrinsic dynamics of each node inside a layer and the functions H αα and H αβ represent the form of theintra- and inter-layer coupling, respectively.We reduce the set of Eqs. (51) to a much more compact form and use it to compute the symmetries of amultiplex network. Introducing x i = x αi x βi ... , F ( x i ) = F α ( x αi ) + (cid:80) λ (cid:54) = α σ αλ H αλ ( x λi ) F β ( x βi ) + (cid:80) λ (cid:54) = β σ βλ H βλ ( x λi )... , H λ ( x i ) = δ λα H αα ( x αi ) δ λβ H ββ ( x βi )... , σ λ = σ λλ A λ = A λλ (52)where λ ∈ { α, β, . . . } , δ αλ is the Kronecker delta (e.g., δ λα = 1 if λ = α , 0 otherwise) we can rewrite Eq.(51) as ˙ x = F ( x ) + Λ (cid:88) λ =1 σ λ A λ H λ ( x ) , (53)which is the equation of a multidimensional network [35, 36], i.e., a network with only one layer but differenttypes of interactions. Each interaction type is described by a different adjacency matrix A λ , λ = 1 , ..., Λ.Blaha et al. [31] recently presented a simple experimental realization of such a network.We have shown the equivalence between a multiplex network and a multidimensional network. Wenow discuss how to compute the symmetries of the multidimensional network (53), both analytically andcomputationally. The analytic approach gives insight into the origins of the final symmetry permutationgroup. The software approach allows a direct calculation of the full group, but without the insights into itsstructure.To analytically define the group of symmetries of the multidimensional network, we introduce G λ as thegroup of symmetries for each interaction type λ . Each element of the group G λ can be represented by a23ermutation matrix Π that commutes with A λ , Π A λ = A λ Π. Then, the symmetry group of the wholemultidimensional network G is given by G ∩ G ... ∩ G Λ , which is the largest common subgroup of the Λsymmetry groups {G λ } .Computationally, the group of symmetries of the multidimensional network G can be found using avail-able computational group theory tools, like GAP or SAGE [55, 56]. This software computes the group ofsymmetries of a labeled graph [57, 58]; a multidimensional network can easily be remapped to a labeledgraph. To remap the network, we define the labeled adjacency matrix A lab . Matrix entries A lab ij are definedby how pairs of nodes i and j interact. For each pair of nodes, one of three cases occurs: (i) there is nointeraction between i and j , (ii) there is one type of weighted interaction between i and j and (iii) thereare two or more types of weighted interactions between i and j . If there is no interaction between node i and node j (case i), then A lab ij = 0. We represent a single edge (case ii) with a 2-tuple ( τ, w ), where theinteger τ is the edge type and the real number w is the edge weight. We represent a multiple edge formedby q interactions (case iii) with a 2 q -tuple, ( τ , w , τ , w , , ..., τ q , w q ), τ < τ < ... < τ q , where each pair( τ i , w i ) represents an interaction type and the associated weight. We can partition the set of the networkedges into Z subsets of edges that are all represented by the same tuple. Then the entries of the matrix A lab are such that A lab ij = 0 if there is no interaction between nodes i and j and A lab ij = z , z = 1 , ..., Z if theedge between nodes i and j belongs to the subset z . This forces the software to consider only permutationsthat involve edges of the same type and of the same weight, i.e., permutation matrices Π that commute with A lab , Π A lab = A lab Π. As a result, the computed permutations Π that form G λ also commute with all theadjacency matrices of the multidimensional network A , . . . , A Λ : G is thus the largest common subgroup ofthe Λ symmetry groups. A network with three layers
In the Results Section we have considered multilayer networks with not more than two layers and only onetype of inter-layer coupling. Here we show how our procedure can be generalized to a multilayer netwotkwith M = 3 layers. Later we will consider the case of any number M of layers. For the case M = 3, thevector field is ˙ xxx α = F α ( xxx α ) + σ αα A αα H αα ( xxx α ) + σ αβ A αβ H αβ ( xxx β ) + σ αγ A αγ H αγ ( xxx γ )˙ xxx β = F β ( xxx β ) + σ ββ A ββ H ββ ( xxx β ) + σ βα A βα H βα ( xxx α ) + σ βγ A βγ H βγ ( xxx γ )˙ xxx γ = F γ ( xxx γ ) + σ γγ A γγ H γγ ( xxx γ ) + σ γα A γα H γα ( xxx α ) + σ γβ A γβ H γβ ( xxx β ) (54)If we apply three permutations g ∈ G α , h ∈ G β , and k ∈ G γ on the system we obtain g ˙ xxx α = F α ( gxxx α ) + σ αα A αα H αα ( gxxx α ) + σ αβ gA αβ H αβ ( xxx β ) + σ αγ gA αγ H αγ ( xxx γ ) h ˙ xxx β = F β ( hxxx β ) + σ ββ A ββ H ββ ( hxxx β ) + σ βα hA βα H βα ( xxx α ) + σ βγ hA βγ H βγ ( xxx γ ) k ˙ xxx γ = F γ ( kxxx γ ) + σ γγ A γγ H γγ ( kxxx γ ) + σ γα kA γα H γα ( xxx α ) + σ γβ kA γβ H γβ ( xxx β ) . (55)For a directed graph, we have 6 conjugacy relationships to be satisfied: gA αβ = A αβ h, gA αγ = A αγ k, hA βα = A βα g, hA βγ = A βγ k, kA γα = A γα g, and kA γβ = A γβ h. If the graph is undirected, then half are redundant ( e.g., gA αβ = A αβ h is the same relationship as hA βα = A βα g ). From these relationships we can define H α = (cid:8) g ∈ G α |∃ h ∈ G β and k ∈ G γ : gA αβ = A αβ h and gA αγ = A αγ k (cid:9) , H β = (cid:8) h ∈ G β |∃ g ∈ G α and k ∈ G γ : hA βα = A βα g and hA βγ = A βγ k (cid:9) , H γ = (cid:8) k ∈ G γ |∃ g ∈ G α and h ∈ G β : kA γα = A γα g and kA γβ = A γβ h (cid:9) . (56)Using the same reasoning as in Theorem 1, each subset in Eqs. (56) is a subgroup of its layer’s group. Theequivalence relationship ∼ is now defined by a set of two equations that must be satisfied, i.e., g ∼ g ⇐⇒ g A αβ = g A αβ and g A αγ = g A αγ K αij = { g ∈ H α : gA αβ = M i , gA αγ = M j } , K βij = { h ∈ H β : hA βα = M i , hA βγ = M j } , K γij = { k ∈ H γ : kA γα = M i , kA γβ = M j } . (57)There are K K disjoint sets for each subgroup H , where K is the number of different M i and K isthe number of different M j that can be obtained applying any of the compatible permutations. Finally, weform the final symmetry group of the multilayer network using the disjoint sets as in Eq. (15), namely, G = g h
00 0 k (cid:12)(cid:12)(cid:12) g ∈ K αij , h ∈ K βij , and k ∈ K γij for i = 1 , ..., K , j = 1 , ..., K . (58)As before with two layers the ILSs for this three-layer case will be in each of the layers’ equivalence classsubgroups, K αij , K βij , and K γij . Two interlayer coupling types
For the case of two or more coupling types between two layers, the extension is similar to the previoussection. The general equation in this case is˙ xxx α = F α ( xxx α ) + σ αα A αα H αα ( xxx α ) + σ αβ, A αβ, H αβ, ( xxx β ) + σ αβ, A αβ, H αβ, ( xxx β )˙ xxx β = F β ( xxx β ) + σ ββ A ββ H ββ ( xxx β ) + σ βα, A βα, H βα, ( xxx α ) + σ βα, A βα, H βα, ( xxx α ) . (59)Applying two permutations g ∈ G α and h ∈ G β to the system, we obtain 4 conjugacy relationships to besatisfied: gA αβ, = A αβ, h, gA αβ, = A αβ, h, hA βα, = A βα, g, hA βα, = A βα, g, , two of which are redundant for undirected graphs. The H groups are still defined by two of these relation-ships H α = (cid:8) g ∈ G α |∃ h ∈ G β : gA αβ, = A αβ, h and gA αβ, = A αβ, h (cid:9) , H β = (cid:8) h ∈ G β |∃ g ∈ G α : hA βα, = A βα, g and hA βα, = A βα, g (cid:9) , (60)and g ∼ g only when both produce the same left-hand-side for all conjugacy relationships. The disjoint K sets for each subgroup are generated as before (with two indices, since there are two constraints that definethe ∼ relationship) and the final symmetry group of the multilayer network using the disjoint sets is as inEq. (15). The ILSs are found as above. Any number of layers and inter-layer couplings
Here we present the general case of any number of layers and any number of inter-layer couplings. Extrap-olating from our discussion above, we can define the group of symmetries of a multilayer network with anynumber of layers and types of coupling. The number of conjugacy relationships to satisfy grows as the num-ber of unique types of inter-layer coupling grows; stated plainly, if we have A inter-layer adjacency matrices,we have A conjugacy relationships. In the case of three layers discussed earlier, there are six such inter-layercouplings; in the case of two layers and two types of interlayer couplings, there are four. The conjugacyrelationships must then be divided in the various layers to define the H groups, and two permutations ineach of the H groups have the equivalence relationship ∼ when they both give the same result in all theleft hand side of the H defining conjugacy relations. Consequently, the computational effort grows with thenumber of inter-layer couplings as O ( A ). 25 omputation of the symmetry group We have already discussed the computation of the symmetry group for the case of multiplex networks andmultidimensional networks. Here we briefly review the computation of the symmetry group for the case ofgeneral multilayer networks. In the presence of both nodes of different types and edges of different types,a symmetry is only allowed: (i) when the edges are interchangeable as discussed before for the case ofmultidimensional networks and (ii) when the nodes that are exchanged are of the same type. Availablecomputational group theory software [55, 56] allows us to handle both aspects, since they also managemulti-partite labeled graph [57, 58].To satisfy requirement (i), we provide the software with a labeled supra-adjacency matrix A lab for themultilayer network, which we define in a similar way as for the previously discussed case of the multidimen-sional network. First, we sequentially renumber all the nodes starting from layer α (so nodes 1 , . . . , N α arethe nodes in layer α , nodes N α + 1 , . . . , N α + N β are the nodes in layer β , and so on). Then we write thematrix A lab that describes the connections between the nodes of the multilayer network: A lab ij = 0 if there isno (intralayer or interlayer) connections between i and j . If they are connected, A lab ij is a tuple as describedearlier.To satisfy requirement (ii), we restrict the search of symmetries to permutations that only involve nodesof the same type (from the same layer). We thus provide the partition P = (cid:8) { , . . . , N α } , { N α + 1 , . . . , N α + N β } , . . . (cid:9) that describes how the renumbered nodes are split between the respective layers. This forces thesoftware to consider only permutations that involve nodes of the same type (within the same subset of thepartition.)Finding the group of symmetries of the multilayer network thus reduces to finding the group of symmetriesof the network described by the labeled supra-adjacency matrix A lab , Π A lab = A lab Π, with a predefinedpartition P of the nodes that are allowed to swap with one another.26 athematical model of the experimental system Applying Kirchhoff’s laws on the circuit, we can write down the following set of differential equations thatgovern the system’s dynamics (we color the coupling terms as in Figure 4):Jumper equations:( L J + L J ) ˙ I J = v J − V J − I J R J − k F J L F ,J ˙ I F, − k F J L F ,J ˙ I F, C J ˙ V J = I J FHN 1: L F, ˙ I F, = V F, − R , I F, − k F J L F ,J ˙ I J C F, ˙ V F, = − f [ V F, ] − I F, + 1 R C [( V ce, − V be, − V F, ) + ( V ce, − V be, − V F, )]FHN 2: L F, ˙ I F, = V F, − R , I F, − k F J L F ,J ˙ I J C F, ˙ V F, = − f [ V F, ] − I F, + 1 R C [( V ce, − V be, − V F, ) + ( V ce, − V be, − V F, )]Colpitts 1: C ce, ˙ V ce, = I L, − I c [ V be, ] + 1 R C ( V F, − ( V ce, − V be, )) C be, ˙ V be, = − ( V ee + V be, ) R ee, − I b [ V be, ] − I L, − R C [ V F, − ( V ce, − V be, )] L C ˙ I L, = V cc − V ce, + V be, − R L, I L, − k C L C ˙ I L, Colpitts 2: C ce, ˙ V ce, = I L, − I c [ V be, ] + 1 R C ( V F, − ( V ce, − V be, )) C be, ˙ V be, = − ( V ee + V be, ) R ee, − I b [ V be, ] − I L, − R C [ V F, − ( V ce, − V be, )] L C ˙ I L, = V cc − V ce, + V be, − R L, I L, − k C L C ˙ I L, Colpitts 3: C ce, ˙ V ce, = I L, − I c [ V be, ] + 1 R C ( V F, − ( V ce, − V be, )) C be, ˙ V be, = − ( V ee + V be, ) R ee, − I b [ V be, ] − I L, − R C [ V F, − ( V ce, − V be, )] L C ˙ I L, = V cc − V ce, + V be, − R L, I L, − k C L C ˙ I L, Colpitts 4: C ce, ˙ V ce, = I L, − I c [ V be, ] + 1 R C ( V F, − ( V ce, − V be, )) C be, ˙ V be, = − ( V ee + V be, ) R ee, − I b [ V be, ] − I L, − R C [ V F, − ( V ce, − V be, )] L C ˙ I L, = V cc − V ce, + V be, − R L, I L, − k C L C ˙ I L, The Colpitts nonlinearity arises from its transistor (2N2222) with gain: I b = (cid:40) , V be ≤ V thV be − V th R on , V be > V th , I c = β f I b . (61)27he FHN nonlinearity arises from the operational amplifier (AD844) which is arranged to provide a piecewiselinear cubic function [49] in the form: f ( V F,i ) = [ V F,i − h ( V F,i )] /R ,i , where h ( V F ) = − V R + , V F ≤ − R ,i V R + R ,i + R ,i (cid:16) R ,i R ,i + 1 (cid:17) V F , − R ,i V R + R ,i + R ,i < V F < R ,i V R + R ,i + R ,i V R + , V F ≥ R ,i V R + R ,i + R ,i . (62)Here, V R + = 6 . V is inferred from the FHN amplitude. Note that the expression h ( V F ) is independent of R ,i and R ,i if R ,i = R ,i .To study the cluster synchronization dynamics of this multilayer system, we vary one coupling and holdthe others constant. We vary the Colpitts intralayer coupling, k C , by finely varying the separation of theinductors of the Colpitts. We immobilize Colpitts layer nodes C and C , and move C and C with a caliper;this precisely and simultaneously changes the C to C and C to C separation. We impose the constantColpitts-FHN intralayer coupling, R C ,with a fixed resistor. We hold the FHN-jumper intralayer coupling, k F J , constant as well; we impose a constant separation between the inductors of the FHN and the jumperwith plastic clips.This experimental testbed differs fundamentally from the numerical examples we consider. Although wemodel our circuit as noise-free, electrical circuits are subject to several kinds of noise, including shot noise,burst noise, and thermal noise [59, 60]. Components such as transistors and operational amplifiers have noiseprofiles, with stronger noise at higher frequencies. Additionally, we neglect stray resistance, inductance,and capacitance which inevitably arises even in circuits which are carefully crafted to minimize such issues.Finally, small parameter mismatches lead the oscillators in each layer to be slightly heterogeneous; uncoupled,the nodes have slightly different frequencies and amplitudes.To write the system in the general form (1), we first assume identical oscillators on each layer, i.e., L J = L J = L F J ; L F ,J = L F ,J = L F J ; L F, = L F, = L F ; C F, = C F, = C F ; R , = R , = R ; C ce, = C ce, = C ce, = C ce, = C ce ; C be, = C be, = C be, = C be, = C be ; L C, = L C, = L C, = L C, = L C ; R ee, = R ee, = R ee, = R ee, = R ee . We can then write˙ x αi = F α ( x αi ) + σ αα N α (cid:88) j =1 A ααij H αα ( x αj ) (cid:124) (cid:123)(cid:122) (cid:125) intra-layer coupling + (cid:88) β (cid:54) = α σ αβ N β (cid:88) j =1 A αβij H αβ ( x βj ) . (cid:124) (cid:123)(cid:122) (cid:125) inter-layer coupling (63)Where x α = I J , x α = V J are the current and the voltage in the jumper; x βi, = I F,i , x βi, = V F,i are thecurrent and the voltage of FHN i ; x γi, = V ce,i , x γi, = V be,i , x γi, = I L,i are the voltages and the current ofColpitts i . Jumper layer: F α ( x α ) = (cid:18) ( R J x α + x α − v J ) / (cid:0) L F J ( k F J L F J − L F ) (cid:1) x α /C J (cid:19) , H αβ ( x βi ) = (cid:18) R x βi, − x βi, )0 (cid:19) , σ αβ = k F J L F − L F J k F J ) , A αβ = (cid:2) (cid:3) FHN layer F β ( x β ) = (cid:18) ( − L F ( R x β − x β ) − k F J L F v J ) / (cid:0) L F ( L F − k F J L F J ) (cid:1) ( − f ( x β ) − x β − x β /R C ) /C F (cid:19) , H ββ ( x βj ) = (cid:18) R x βj, − x βj, (cid:19) , σ ββ = L F J k F J L F ( L F J k F J − L F ) , A ββ = (cid:20) − − (cid:21) H βα ( x α ) = (cid:18) R J x α + x α (cid:19) , σ βα = k F J L F − L F J k F J ) = σ αβ , A βα = (cid:20) (cid:21) = ( A αβ ) T . βγ ( x β , x γ ) = (cid:18) x γ − x γ ) /C F (cid:19) , σ βγ = 1 R C , A βγ = (cid:20) (cid:21) . Note that the inductive coupling through the jumper generates a diffusive coupling in the β layer (as canbe seen in A ββ ). This is a consequence of the fact that the inductive coupling acts on the derivative ofthe current ˙ I F,i , see the orange terms in the Jumper equations. As a result, this induces a direct diffusivecoupling between the currents I F,i of the FHN circuits. This is represented in figure 4 by a dashed line inthe β layer. Colpitts layerF γ ( x γ ) = ( x γ − I C ( x γ ) − ( x γ − x γ ) /R C ) /C ce ( − ( V ee + x γ ) /R ee − I b ( x γ ) − x γ + ( x γ − x γ ) /R C )) /C be ( x γ − x γ − R L x γ + V CC (1 − k C )) / ( L (1 − k C )) , H γγ ( x γj ) = − ( x γj, − x γj, − R L x γj, ) , σ γγ = k C L (1 − k C ) , A γγ = H γβ ( x β , x γ ) = x β /C ce x β /C be , σ βγ = 1 R C = σ γβ , A γβ = ( A βγ ) T . Acknowledgements
This work was supported by the National Science Foundation through grant (Grant No. 1727948) and theOffice of Naval Research through the Naval Research Laboratory’s Basic Research Program. We thankB.Hunt for noting that the equivalence classes of each layer form cosets of a subgroup of that layer’s groupand M. Hossein-Zadeh and K. Huang for their support with the experimental activities.
Author contributions
F.D.R., L.P., and F.S. worked on the theory. K.B. performed the experiment with the help of F.D.R. A.S.worked on the analysis of real data. I.K. designed the algorithm for generating multilayer networks withsymmetries. F.S. supervised all aspects of the project. F.S. wrote the paper with the help of K.B. and L.P.
Competing interests
The authors declare no competing interests.
References [1] Buldyrev, S. V., Parshani, R., Paul, G., Stanley, H. E. & Havlin, S. Catastrophic cascade of failures ininterdependent networks.
Nature , 1025–1028 (2010).[2] Korkali, M., Veneman, J. G., Tivnan, B. F., Bagrow, J. P. & Hines, P. D. Reducing cascading failurerisk by increasing infrastructure network interdependence.
Scientific Reports , 44499 (2017).293] Rezai, A., Keshavarzi, P. & Moravej, Z. Key management issue in scada networks: A review. Engineeringscience and technology, an international journal , 354–363 (2017).[4] Rosato, V. et al. Modelling interdependent infrastructures using interacting dynamical models.
Inter-national Journal of Critical Infrastructures , 63–18 (2008).[5] Pereda, A. E. Electrical synapses and their functional interactions with chemical synapses. NatureReviews Neuroscience , 250–263 (2014).[6] Song, X., Wang, C., Ma, J. & Tang, J. Transition of electric activity of neurons induced by chemicaland electric autapses. Science China Technological Sciences , 1007–1014 (2015).[7] Adhikari, B. M., Prasad, A. & Dhamala, M. Time-delay-induced phase-transition to synchrony incoupled bursting neurons. Chaos: An Interdisciplinary Journal of Nonlinear Science , 023116 (2011).[8] Sorrentino, F. Synchronization of hypernetworks of coupled dynamical systems. New Journal of Physics , 033035 (2012).[9] Goulding, M. Circuits controlling vertebrate locomotion: moving in a new direction. Nature ReviewsNeuroscience , 507 (2009).[10] Lodi, M., Shilnikov, A. & Storace, M. Design of synthetic central pattern generators producing desiredquadruped gaits. IEEE Transactions on Circuits and Systems I , 1028–1039 (2018).[11] Danner, S. M., Wilshin, S. D., Shevtsova, N. A. & Rybak, I. A. Central control of interlimb coordinationand speed-dependent gait expression in quadrupeds. The Journal of Physiology , 6947–6967 (2016).[12] Kivela, M. et al.
Multilayer networks.
Journal of Complex Networks , 203–271 (2014).[13] Boccaletti, S. et al. The structure and dynamics of multilayer networks.
Physics Reports , 1–122(2014).[14] Taylor, D., Shai, S., Stanley, N. & Mucha, P. J. Enhanced detectability of community structure inmultilayer networks through layer aggregation.
Physical review letters , 228301 (2016).[15] Irving, D. & Sorrentino, F. Synchronization of a hypernetwork of coupled dynamical systems.
PhysicalReview E , 056102 (2012).[16] del Genio, C. I., G´omez-Garde˜nes, J., Bonamassa, I. & Boccaletti, S. Synchronization in networks withmultiple interaction layers. Science Advances , e1601679 (2016).[17] Belykh, V. N., Belykh, I. V. & Mosekilde, E. Cluster synchronization modes in an ensemble of coupledchaotic oscillators. Physical Review E , 036216 (2001).[18] Belykh, V. N., Osipov, G. V., Petrov, V. S., Suykens, J. A. & Vandewalle, J. Cluster synchronizationin oscillatory networks. Chaos: An Interdisciplinary Journal of Nonlinear Science , 037106 (2008).[19] Nicosia, V., Valencia, M., Chavez, M., D´ıaz-Guilera, A. & Latora, V. Remote synchronization revealsnetwork symmetries and functional modules. Physical Review Letters , 174102 (2013).[20] Pecora, L., Sorrentino, F., Hagerstrom, A., Murphy, T. & Roy, R. Cluster synchronization and isolateddesynchronization in complex networks with symmetries.
Nature Communications , 4079 (2014).[21] Sorrentino, F., Pecora, L. M., Hagerstrom, A. M., Murphy, T. E. & Roy, R. Complete characterizationof the stability of cluster synchronization in complex dynamical networks. Science Advances , e1501737(2016).[22] Sorrentino, F. & Pecora, L. Approximate cluster synchronization in networks with symmetries andparameter mismatches. Chaos: An Interdisciplinary Journal of Nonlinear Science , 094823 (2016).3023] Cho, Y. S., Nishikawa, T. & Motter, A. E. Stable chimeras and independently synchronizable clusters. Physical Review Letters , 084101 (2017).[24] Siddique, A. B., Pecora, L., Hart, J. D. & Sorrentino, F. Symmetry-and input-cluster synchronizationin networks.
Physical Review E , 042217 (2018).[25] Golubitsky, M. & Stewart, I. Rigid patterns of synchrony for equilibria and periodic cycles in networkdynamics. Chaos: An Interdisciplinary Journal of Nonlinear Science , 094803 (2016).[26] Belykh, I., Belykh, V., Nevidin, K. & Hasler, M. Persistent clusters in lattices of coupled nonidenticalchaotic systems. Chaos: An Interdisciplinary Journal of Nonlinear Science , 165–178 (2003).[27] Schaub, M. T. et al. Graph partitions and cluster synchronization in networks of oscillators.
Chaos:An Interdisciplinary Journal of Nonlinear Science , 094821 (2016).[28] MacArthur, B. D., S´anchez-Garc´ıa, R. J. & Anderson, J. W. Symmetry in complex networks. DiscreteApplied Mathematics (2008).[29] Klickstein, I. S. & Sorrentino, F. Generating graphs with symmetry.
IEEE Transactions on NetworkScience and Engineering (2018).[30] Skardal, P. S. Symmetry and symmetry breaking in coupled oscillator communities.
The EuropeanPhysical Journal B , 46 (2019).[31] Blaha, K. A. et al. Cluster synchronization in multilayer networks: A fully analog experiment with l coscillators with physically dissimilar coupling.
Physical Review Letters , 014101 (2019).[32] Verbrugge, L. M. Multiplexity in adult friendships.
Social Forces , 1286–1309 (1979).[33] Sol´a, L. et al. Eigenvector centrality of nodes in multiplex networks.
Chaos: An InterdisciplinaryJournal of Nonlinear Science , 033131 (2013).[34] Gomez, S. et al. Diffusion dynamics on multiplex networks.
Physical Review Letters , 028701 (2013).[35] Berlingerio, M., Coscia, M., Giannotti, F., Monreale, A. & Pedreschi, D. Multidimensional networks:foundations of structural analysis.
World Wide Web , 567–593 (2013).[36] Coscia, M., Rossetti, G., Pennacchioli, D., Ceccarelli, D. & Giannotti, F. You know because i know: amultidimensional network approach to human resources problem. In Proceedings of the 2013 IEEE/ACMInternational Conference on Advances in Social Networks Analysis and Mining , 434–441 (ACM, 2013).[37] Tang, L., Wu, X., L¨u, J., Lu, J.-a. & D’Souza, R. M. Master stability functions for complete, intralayer,and interlayer synchronization in multiplex networks of coupled r¨ossler oscillators.
Physical Review E , 012304 (2019).[38] MacArthur, B. D. & S´anchez-Garc´ıa, R. J. Spectral characteristics of network redundancy. PhysicalReview E , 026117 (2009).[39] Sorrentino, F., Siddique, A. B. & Pecora, L. M. Symmetries in the time-averaged dynamics of networks:Reducing unnecessary complexity through minimal network models. Chaos: An Interdisciplinary Jour-nal of Nonlinear Science , 011101 (2019).[40] De Domenico, M., Sol´e-Ribalta, A., G´omez, S. & Arenas, A. Navigability of interconnected networksunder random failures. Proceedings of the National Academy of Sciences , 8351–8356 (2014).[41] Cardillo, A. et al.
Emergence of network features from multiplexity.
Scientific reports , 1344 (2013).3142] De Domenico, M., Lancichinetti, A., Arenas, A. & Rosvall, M. Identifying modular flows on multilayernetworks reveals highly overlapping organization in interconnected systems. Physical Review X ,011027 (2015).[43] Coleman, J., Katz, E. & Menzel, H. The diffusion of an innovation among physicians. Sociometry ,253–270 (1957).[44] De Domenico, M., Nicosia, V., Arenas, A. & Latora, V. Structural reducibility of multilayer networks. Nature communications , 6864 (2015).[45] Stark, C. et al. Biogrid: a general repository for interaction datasets.
Nucleic acids research ,D535–D539 (2006).[46] Rossi, R. A. & Ahmed, N. K. The network data repository with interactive graph analytics andvisualization. In AAAI (2015). URL http://networkrepository.com .[47] Klickstein, I., Pecora, L. & Sorrentino, F. Symmetry induced group consensus.
Chaos: An Interdisci-plinary Journal of Nonlinear Science , 073101 (2019).[48] Clifford, A. H. Representations induced in an invariant subgroup. Annals of Mathematics
IEEE Transactionson Systems, Man, and Cybernetics
SMC-13 , 1010–1014 (1983).[50] Kennedy, M. Chaos in the colpitts oscillator.
IEEE Transactions on Circuits and Systems , 771–774(1994).[51] Ishizaki, T., Chakrabortty, A. & Imura, J.-I. Graph-theoretic analysis of power systems. Proceedings ofthe IEEE , 931–952 (2018).[52] MacArthur, B. D., S´anchez-Garc´ıa, R. J. & Anderson, J. W. Symmetry in complex networks.
DiscreteApplied Mathematics , 3525–3531 (2008).[53] Sorrentino, F., Pecora, L. M. & Trajkovic, L. Group consensus in multilayer networks.
IEEE Transac-tions on Network Science and Engineering (early access) (2020).[54] Leyva, I. et al.
Relay synchronization in multiplex networks.
Scientific Reports , 1–11 (2018).[55] Group, G. et al. Gap–groups, algorithms, and programming, version 4.4. 12 (2008). .[56] Stein, W. et al.
Sage: Software for algebra and geometry experimentation (2006). .[57] McKay, B. D. et al. Practical graph isomorphism (Department of Computer Science, Vanderbilt Uni-versity Tennessee, US, 1981).[58] McKay, B. D. & Piperno, A. Practical graph isomorphism, ii.