Cluster synchronization and isolated desynchronization in complex networks with symmetries

Abstract

Synchronization is of central importance in power distribution, telecommunication, neuronal and biological networks. Many networks are observed to produce patterns of synchronized clusters, but it has been difficult to predict these clusters or understand the conditions under which they form. Here we present a new framework and develop techniques for the analysis of network dynamics that shows the connection between network symmetries and cluster formation. The connection between symmetries and cluster synchronization is experimentally confirmed in the context of real networks with heterogeneities and noise using an electro-optic network. We experimentally observe and theoretically predict a surprising phenomenon in which some clusters lose synchrony without disturbing the others. Our analysis shows that such behaviour will occur in a wide variety of networks and node dynamics. The results could guide the design of new power grid systems or lead to new understanding of the dynamical behaviour of networks ranging from neural to social.

ARTICLE

Received 23 Sep 2013 | Accepted 2 May 2014 | Published 13 Jun 2014

Synchronization is of central importance in power distribution, telecommunication, neuronal and biological networks. Many networks are observed to produce patterns of synchronized clusters, but it has been difcult to predict these clusters or understand the conditions under which they form. Here we present a new framework and develop techniques for the analysis of network dynamics that shows the connection between network symmetries and cluster formation. The connection between symmetries and cluster synchronization is experimentally conrmed in the context of real networks with heterogeneities and noise using an electro-optic network. We experimentally observe and theoretically predict a surprising phenomenon in which some clusters lose synchrony without disturbing the others. Our analysis shows that such behaviour will occur in a wide variety of networks and node dynamics. The results could guide the design of new power grid systems or lead to new understanding of the dynamical behaviour of networks ranging from neural to social.

DOI: 10.1038/ncomms5079

Cluster synchronization and isolated desynchronization in complex networks with symmetries

Louis M. Pecora1, Francesco Sorrentino2, Aaron M. Hagerstrom3,4, Thomas E. Murphy4,5 & Rajarshi Roy4,6,7

1 US Naval Research Laboratory, Department of Material Science and Technology, Washington, District of Columbia 20375, USA. 2 Department of Mechanical Engineering, University of New Mexico, Albuquerque, New Mexico 87131, USA. 3 Department of Physics, University of Maryland, College Park, Maryland 20742, USA. 4 Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland 20742, USA. 5 Department of Electrical and Computer Engineering, University of Maryland, College Park, Maryland 20742, USA. 6 Department of Physics, University of Maryland, College Park, Maryland 20742, USA. 7 Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742, USA. Correspondence and requests for materials should be addressed to L.M.P. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5079

Synchronization in complex networks is essential to the proper functioning of a wide variety of natural and engineered systems, ranging from electric power grids to

neural networks1. Global synchronization, in which all nodes evolve in unison, is a well-studied effect, the conditions for which are often related to the network structure through the master stability function2. Equally important, and perhaps more ordinary, is partial or cluster synchronization (CS), in which patterns or sets of synchronized elements emerge3. Recent work on CS has been restricted to networks where the synchronization pattern is induced either by tailoring the network geometry or by the intentional introduction of heterogeneity in the time delays or node dynamics411. These anecdotal studies illustrate the interesting types of CS that can occur, and suggest a broader relationship between the network structure and synchronization patterns. Recent studies have begun to draw a connection between network symmetry and CS, although all have considered simple networks where the symmetries are apparent by inspection1214. More in-depth studies have been done involving bifurcation phenomena and synchronization in ring and point symmetry networks15,16. Finally, a form of CS can occur in situations where input and output couplings match in a cluster, but there is no symmetry16,17.

Here we address the more common case where the intrinsic network symmetries are neither intentionally produced nor easily discerned. We present a comprehensive treatment of CS, which uses the tools of computational group theory to reveal the hidden symmetries of networks and predict the patterns of synchronization that can arise. We use irreducible group representations to nd a block diagonalization of the variational equations that can predict the stability of the clusters. We further establish and observe a generic symmetry-breaking bifurcation termed isolated desynchronization, in which one or more clusters lose synchrony while the remaining clusters stay synchronized. The analytical results are conrmed through experimental measurements in a spatiotemporal electro-optic network. By statistically analysing the symmetries of several types of networks, as well as electric power distribution networks, we argue that symmetries, clusters and isolated desynchronization can occur in many types of complex networks and are often found concurrently. Throughout the text, we use the abbreviation ID for isolated desynchronization.

ResultsThe dynamical equations. A set of general dynamical equations to describe a network of N-coupled identical oscillators are

_

xit F xit

sXjAijHxj; i 1; ; N; 1

where xi is the n-dimensional state vector of the ith oscillator, F describes the dynamics of each oscillator, A is a coupling matrix that describes the connectivity of the network, s is the overall coupling strength and H is the output function of each oscillator. Equation (1) or its equivalent forms provide the dynamics for many networks of oscillators, including all those in refs 14,6 8,1114. This includes some cases of time delays in the coupling functions. As noted in ref. 1, it is only necessary for the form of the equations of motion or, more importantly, the variational equations to have the form of equation (1) near the synchronization manifolds. The form of equation (1) also applies to discrete time systems or more general coupling schemes18. In addition, the same form emerges in the more general case for the variational equations where the vector eld and the coupling combine into one function, viz. F(xi(t), {xj(t)}), where the argument in brackets is the input of all nodes connected to

node i, so long as the nodes are treated as having the same basic dynamics. See ref. 18 for more explanation. Generalizations of equation (1) have been studied in refs 5,9,10,19. The types of natural and man-made systems, which can be modelled by equations of the same form as equation (1), are large20,21. These include genetic networks, circadian networks, ecology, neuronal networks, cortical networks, consensus problems, opinion formation, power grids and concentration of metabolites in a cell to name a few.

The same analysis as presented here will apply to more general dynamics. Here, for simplicity, we take all nodes to have identical dynamics and be bidirectionally coupled to other nodes in the network by couplings of the same weight, that is, Aij is taken to be the (symmetric) adjacency matrix of 1s and 0s with the factor of s controlling the weight of the couplings. The vector eld F can contain self-feedback terms. Here we take the self-feedback to be identical for all nodes.

The cases of Laplacian coupling schemes are encompassed by our analysis as shown here. This is because the row sums are not affected by the symmetry operations. Laplacian coupling schemes are often used to tune the network specically for global synchronization, which is not our goal here. We think our scheme may be more representative of networks that form naturally (for example, neurons) where row sums will not necessarily appear in the feedback. The cases of directional (non-symmetric) coupling, variations of coupling weight and non-identical nodes can all be handled in much the same way as presented here, although it is not clear that those more general cases will have as many symmetries as our simpler case.

The symmetries of the network form a (mathematical) group

G. Each symmetry g of the group can be described by a permutation matrix Rg that re-orders the nodes in a way that leaves the dynamical equations unchanged (that is, each Rg commutes with A). The set of symmetries (or automorphisms)15,22 of a network can be quite large, even for small networks, but it can be calculated from A using widely available discrete algebra routines23,24. Figure 1a shows three graphs generated by randomly removing 6 edges from an otherwise fully connected 11-node network. Although the graphs appear similar and exhibit no obvious symmetries, the rst instance has no symmetries (other than the identity permutation), while the others have 32 and 5,760 symmetries, respectively. Thus, for even a moderate number of nodes (11), nding the symmetries can become impossible by inspection.

Once the symmetries are identied, the nodes of the network can be partitioned into M clusters by nding the orbits of the symmetry group: the disjoint sets of nodes that when all of the symmetry operations are applied permute among one another15. As equation (1) is essentially unchanged by the permutations, the dynamics of the nodes in each cluster can be equal, which is exact synchronization. Hence, there are M synchronized motions {s1,y, sM}, one for each cluster. In Fig. 1a, the nodes have been coloured to show the clusters. For the rst example, which has no symmetries, the network divides into M N trivial clusters with

one node in each. The other instances have ve and three clusters, respectively. Once the clusters are identied, equation (1) can be linearized about a state where synchronization is assumed among all of the nodes within each cluster. This linearized equation is the variational equation and it determines the stability of the clusters.

The variational equations. Equation (1) is expressed in the node coordinate system, where the subscripts i and j are identied with enumerated nodes of the network. Beyond identifying the symmetries and clusters, group theory also provides a powerful way to transform the variational equations to a new

2 NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5079 ARTICLE

a

" #

Zt;

4

XMm1Jm DH smt

where we have linearized about synchronized cluster states {s1,...,sM}, Z(t) is the vector of variations of all nodes transformed to the IRR coordinates, and DF and DH are the Jacobians of the nodes vector eld and coupling function, respectively. We note that this analysis holds for any node dynamics, steady state, periodic, chaotic and so on.

We can write the block diagonal B as a direct sum

Ll1Id

l

0 symmetries

(11 clusters)

32 symmetries (5 clusters)

5,760 symmetries (3 clusters)

b

Cl, where Cl is a (generally complex) pl pl matrix

with pl the multiplicity of the lth IRR in the permutation

representation {Rg}, L the number of IRRs present and d(l)

the dimension of the lth IRR, so that

A =

B =

A =

A =

PLl1 dlpl N (refs 25,26).For many transverse blocks, Cl is a scalar, that is, pl 1. However,

the trivial representation (l 1), which is associated with the

motion in the synchronization manifold has p1 M. Note that

the vector eld F can contain a self-feedback term bxi as in the experiment and other feedbacks are possible, for example, row sums of Aij, as long as those terms commute with the Rg. In all these cases, the B matrix will have the same structure.

Figure 1c shows the coupling matrix B in the IRR coordinate system for the three example networks. The upper-left block is an M M matrix that describes the dynamics within the synchro

nization manifold. The remaining diagonal blocks describe motion transverse to this manifold and so are associated with loss of synchronization. Thus, the diagonalization completely decouples the transverse variations from the synchronization block and partially decouples the variations among the transverse directions. In this way the stability of the synchronized clusters can be calculated using the separate, simpler, lower-dimensional ODEs of the transverse blocks to see if the non-synchronous transverse behaviour decays to zero.

An electro-optic experiment. Figure 2a shows the optical system used to study cluster synchronization. Light from a 1,550-nm light-emitting diode passes through a polarizing beamsplitter and quarter wave plate, so that it is circularly polarized when it reaches the spatial light modulator (SLM). The SLM surface imparts a programmable spatially dependent phase shift x between the polarization components of the reected signal, which is then imaged, through the polarizer, onto an infrared camera27. The relationship between the phase shift x applied by the SLM and the normalized intensity I recorded by the camera

is I(x) (1 cosx)/2. The resulting image is then fed back

through a computer to control the SLM. More experimental details are given in the Methods section.

The dynamical oscillators that form the network are realized as square patches of pixels selected from a 32 32 tiling of the SLM

array. Figure 2b shows an experimentally measured camera frame captured for 1 of the 11-node networks considered earlier in Fig. 1 (A full video is presented in the Supplementary Movie 1). The patches have been falsely coloured to show the cluster structure and the links of the network are overlaid to illustrate the connectivity. The phase shift of the ith region, xi, is updated iteratively according to:

xt 1i bI xti

s

c

B =

B =

Figure 1 | Three randomly generated networks with varying amounts of symmetry and associated coupling matrices. (a) Nodes of the same colour are in the same synchronization cluster. The colours show the maximal symmetry the network dynamics can have given the graph structure. (b) A graphic showing the structure of the adjacency matrices of each network (black squares are 1, white squares are 0). (c) Block diagonalization of the coupling matrices A for each network. Colours denote the cluster, as in a. The 2 2 transverse block for the 32 symmetry case

comes from one of the IRRs being present in the permutation matrices two times. The matrices for the 32 symmetry case are shown in the Methods section.

coordinate system in which the transformed coupling matrix B TAT 1 has a block diagonal form that matches the cluster

structure, as described below. The transformation matrix T is not a simple node re-ordering, nor is it an eigendecomposition of A. The process for computing T is non-trivial, and involves nding the irreducible representations (IRRs) of the symmetry group. We call this new coordinate system the IRR coordinate system. In the Methods section, we show the steps necessary to obtain the symmetries, clusters and the transformation T.

Once we have T, we can transform the variational equations as follows. Let Cm be the set of nodes in the mth cluster with

synchronous motion sm(t). Then, the original variational equations about the synchronized solutions are (in vectorial form and in the node coordinates),

d_

xt

XMm1Em DF smt sA

" #

dxt;

2

XMm1Em DH smt

where the Nn-dimensional vector dx(t) [dx1(t)T, dx2(t)T,...,

dxN(t)T]T and E(m) is an N-dimensional diagonal matrix such that

Emii

1; if i 2 Cm;

0; otherwise;

3

" #

i 1,y,N. Note that

XjAijI xtj

d

PMm1 Em IN, where IN is the N-dimensional identity matrix.

Applying T to equation (2), we arrive at the variational matrix equation shown in equation (4), where g(t) T # In dx(t), J(m) is

the transformed E(m) and B is the block diagonalization of the coupling matrix A,

_

Zt

XMm1Jm DF smt sB In

mod 2p 5

where b is the self-feedback strength and the offset d is introduced to suppress the trivial solution xi 0. Equation (5)

NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5079

a b

SLM

1 mm

LED

Camera

QWP

PBS

Lens

Computer

Figure 2 | Experimental conguration. (a) Light is reected from the SLM and passes though polarization optics thus, that the intensity of light falling on the camera is modulated according the phase shift introduced by the SLM. Coupling and feedback are implemented by a computer. (b) An imageof the SLM recorded by the camera in this conguration. Oscillators are shaded to show which cluster they belong to, and the connectivity of the network is indicated by superimposed grey lines. The phase shifts applied by the square regions are updated according to equation (5).

is a discrete time equivalent of equation (1). Depending on the values of b, s and d, equation (5) can show constant, periodic or chaotic dynamics. There are no experimentally imposed constraints on the adjacency matrix Aij, which makes this system an ideal platform to explore synchronization in complex networks.

Figure 3 plots the time-averaged root-mean square synchronization error for all four of the non-trivial clusters shown in Fig. 2b, as a function of the feedback strength b, holding s constant. We nd qualitatively similar results if s is chosen as a bifurcation parameter with b held constant. In Fig. 3ce, we plot the observed intracluster deviations xti xt for three specic

values of b indicated by the vertical lines in Fig. 3a,b, showing different degrees of partial synchronization that can occur, depending on the parameters.

Details of the calculations are given in the Methods section. Supplementary Movie 1 illustrates the experimentally recorded network behaviour for the case of b 0.72p, where the system

clearly partitions into four synchronized clusters plus one unsynchronized node.

Together, Fig. 3a,ce illustrate two examples of a bifurcation commonly seen in experiment and simulation: isolated desynchronization, where one or more clusters lose stability, while all others remain synchronized. At b 0.72p (Fig. 3c), all four of the

clusters synchronize. At b 1.4p (Fig. 3d), the magenta cluster,

which contains four nodes, has split into two smaller clusters of two nodes each, while the other two clusters remain synchronized.

Between b 0.72p and b 1.767p, two clusters, shown in

Fig. 1 as red and blue, respectively, undergo isolated desynchronization together. In Fig. 3a, the synchronization error curves for these two clusters are visually indistinguishable. The synchronization of these two clusters is intertwined: they will always either synchronize together or not at all. Although it is not obvious from a visual inspection of the network that the red and blue clusters should form at all, their intertwined synchronization properties can be understood intuitively by examining the connectivity of the network. Each of the two nodes in the blue cluster is coupled to exactly one node in the red cluster. If the blue cluster is not synchronized, the red cluster cannot synchronize, because its two nodes are receiving different input. The group analysis treats this automatically and yields a transverse 2 2 block in Fig. 1c.

The isolated desynchronization bifurcations we observe are predicted by computation of the maximum Lyapunov exponent (MLE) of the transverse blocks of equation (4), shown in Fig. 3. The region of stability of each cluster is predicted by a negative MLE. Although there are four clusters in this network, there are

only three MLEs: the two intertwined clusters are described by a two-dimensional block in the block-diagonalized coupling matrixB. These stability calculations reveal the same bifurcations as seen in the experiment.

Subgroup decomposition and isolated desynchronization. The existence of isolated desynchronizations in the network experiments raises several questions. As the network is connected why doesnt the desynchronization pull other clusters out of sync? What is the relation of ID to cluster structure and network symmetry? Is ID a phenomenon that is common to many networks? We provide answers to all these questions using geometric decomposition of a group, which was developed in refs 28,29. This technique enables a nite group to be written as a direct product of subgroups G H1 ... HV where V is the number of

subgroups and all the elements in one subgroup commute with all the elements in any other subgroup. This means that the set of nodes permuted by one subgroup is disjoint from the set of nodes permuted by any other subgroup. Then, each cluster (say, Cj) is

permuted only by one of the subgroups (say, Hk), but not by any

others. There can be several clusters permuted by one subgroup. This is the case of the red and blue clusters in the 32 symmetry network in Fig. 1, because the associated Hk cannot have a geo

metric decomposition, but may have a more structured decomposition such as a wreath product30.

We can show that the above decomposition guarantees that the nodes associated with different subgroups all receive the same total input from the other subgroups nodes. Hence, nodes of each cluster will not see the effects of individual behaviour of the other clusters associated with different subgroups. This enables the clusters to have the same synchronized dynamics even when another cluster desynchronizes. If that state is stable, we have ID.

To see this let Hk, a subgroup of G, permute only cluster Cm

and p be the permutation on the indices of nodes in Cm for one

permutation Rg, g A Hk. Assume xi is not in Cm; thus, it is not

permuted by Rg and recall that A commutes with all permutations in G, then we have (just concentrating on the terms from Cm),

Rg _

xt

i

_

s ARgHx

i

s

xit s RgAHx

i

XjAijH xpj

;

6

where p(l) is, in general, another node in Cm and the sums over

other clusters are unchanged. This shows that all nodes in Cm are

coupled into the ith node in the same way (the same Aij factor). Similarly, if we use a permution Rg0 on the cluster Cm

0 containing

4 NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5079 ARTICLE

1.5

a

c d e

1.0

0.8

0

0.8

x

2

2

x RMS

0.5

x

x

2

2

0

b

Cluster 1 Clusters 2 & 3 Cluster 4

c d e

Cluster transverse

Lyapunov exponent

2

2

[afii9811]

2

0 [afii9811] 0 30 0 30 0 30

2[afii9811]

[afii9811]

2

3[afii9811] 2

t t t

[afii9826]

Figure 3 | Experimental observation of isolated and intertwined desynchronization. (a) CS error as the self-feedback, b, is varied. For all cases considered, d 0.525 and s 0.67p. Colours indicate the cluster under consideration and are consistent with Fig. 1. (b) MLE calculated from simulation.

(ce) Synchronization error time traces for the four clusters, showing the isolated desychronization of the magenta cluster and the isolated desychronization of the intertwined blue and red clusters.

xi, we can show that all the nodes of Cm are coupled in the same

way to the nodes in Cm. Hence, nodes of Cm each receive the same

input sum from the nodes of Cm whether the nodes of Cm are

synchronized or not. This explains how the cluster Cm can

become desynchronized, but the nodes of Cm can still be

synchronizedthey all have the same input despite the Cm

desynchronization, thus making the Cm synchronous state ow

invariant. If it is also stable, this is the case of ID. This argument is easily generalized to the case when Hk permutes nodes of

several clusters as this will just add other similar sums to equation (6). The latter case explains the intertwined desynchronization in the experiment and is a more general form of ID.

Symmetries for different network topologies. How common is such an ID situation we outlined above? We have examined symmetry statistics for some classes of random and semi-random graph types, which suggest that when symmetries are present the opportunity for ID dynamics will be common, although the stability for such will depend on the dynamical systems of the network nodes. We emphasize that we are not attempting to compare symmetry statistics across different graph models, but only to generate networks that have different topologies (for example, random, tree-like and bipartite) to show common phenomena.

We examined 10,000 realizations of three random and semi-random networks each with 100 nodes: (1) randomly connected nodes (random graphs) similar to ErdsRenyi graphs31 generated as described above but with 100 nodes and 50 random edges deleted; (2) scale-free (tree) graphs following Barabasi and Albert32,33; and (3) random bipartite23 graphs using the RandomBipartite function in Sage with 90 nodes in 1 partition and 10 in the other.

Random graphs were generated by starting with 100 nodes completely connected and randomly deleting 50 edges. Scale-free Barabasi and Albert graphs were based on the original Barabasi and Albert preference algorithm32 using the Sage routine RandomBarabasiAlbert. These had 100 nodes with 99 edges and a tree structure. Bipartite networks were generated with the Sage routine RandomBipartite (n1, n2, p), which

generates a graph with two sets of nodes (n1 in the rst

partition and n2 in the second partition, where n1 n2 100) and

connections from nodes in the rst partition to nodes in the second partition are added with probability p. We used n1 10,

n2 90 and p 0.85. Ten thousand realizations of each graph

type were generated. We tested several 10,000 realizations and we see very little variation in statistics between realizations of the same class leading us to believe that we are sampling fairly and enough to trust our results. We also checked for equivalent (isomorphic) graphs to see how much repetition we had. The random systems yielded on average 1 equivalent pair per 10,000 realizations. The scale-free cases yielded o0.01% equivalent graphs. Apparently, we are not near the maximum number of inequivalent graphs for any of the classes. Even with just 100 realizations, the main trends in number of symmetries and other statistics are evident, although such small samples occasionally miss those symmetry cases that are not too common in the class.

Figure 4 shows the different topologies of the networks. Figure 5a shows the cumulative distribution of symmetries for each type of network. All have similar distributions overall, but on different scales of symmetries. Almost all graphs for each type have several non-trivial clusters and more than one subgroup. The median numbers of clusters for the random, Barabasi and Albert, and bipartite networks are 11, 19 and 15, respectively. The median numbers of subgroups are 8, 15 and 15, respectively. The per cent of cases where the number of subgroups is less than the number of clusters (intertwined cases) are 98.45, 98.48 and 0%, respectively. Thus, in these networks symmetries, clusters and subgroups are simultaneously present and, thus, so is the scenario for ID.

We also studied symmetries, clusters and subgroup decompositions in small-world graphs. Small-world graphs31,34 were generated by starting with a ring of nearest-neighbour connected nodes, then adding a xed number of edges to give the same number of edges as the random graphs in the text. We found we had to add many edges beyond the usual few used to generate the small-world effect, because adding only a few edges beyond the ring rarely resulted in any symmetries. As a result, the small-world examples approached being similar to the random graphs so we do not display their results, although the two systems each have symmetries that the other does not; hence, they appear to not be exactly identical.

NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5079

a

b

c

Figure 4 | Three graph topologies. Graph types analysed are (a) complete graph with random edges deleted, (b) BarabasiAlbert tree graphs and (c) bipartite graphs. The number of nodes shown is smaller (2025) than the actual number used in the calculations (100) so that the topological structure of each type is clear.

a b c

1.0

0.20

0.3

0.8

Random BA

BP

Random BA

BP

0.16

Probability

0.6

0.12

0.2

Random BA

BP

0.4

0.08

0.1

0.2

0.04

0.0

0.00

0.0

106 1026 1046 1066 1086

0 10 20 30

Number of symmetries

Number of nontrivial clusters

0 5 10 15 20 25

Number of subgroups

Figure 5 | Symmetry, cluster and subgroup statistics for three types of networks. The networks are random, Barabasi and Albert (BA in the gure), and the bipartite case (BP in the gure). The statistics are (a) the cumulative distribution of the number of symmetries (the dashed line is the median), (b) the counts of the number of non-trivial clusters and (c) the counts of the number of subgroups in the decomposition. The statistics were calculated using 10,000 realizations of each network.

Finally, we examined some man-made and natural networks. The man-made networks are the Nepal power grid35 and the Mesa del Sol electrical grid36. We show the Nepal grid, as its small size is easier to display in Fig. 6. In addition, shown is the block diagonalization of the coupling matrix. Here we treat the grid analogous to ref. 1 in which all power stations are identical with the same bidirectional coupling along each edge which ts the model in equation (1). In this model, global synchronization (the usual desired state) is a solution of the equations of motion. However, as we have pointed out, if symmetry induced clusters exist, CS is also a solution of the equations of motion. The Nepal network has 86,400 symmetries, three non-trivial clusters (plus two trivial ones) and three subgroups (one for each nontrivial cluster). This implies it is possible for this network to split into three sets of synchronized clusters. This also implies that depending on the exact dynamics, the parameters and the stability of global and cluster states, it is possible that the CS state may be a route to desynchronization of the global state in the Nepal grid or other grids with similar dynamics.

The Mesa del Sol electric grid has 4,096 symmetries, 20 nontrivial clusters and 10 subgroups. The network has three intertwined clusters, two with four clusters and one with ve clusters. This is shown in Fig. 7 a circle plot of the Mesa del Sol network, which, because of the network size (132 nodes), exposes the cluster structure much better. As the Nepal example indicates, the symmetry structure makes global synchronization, CS and ID a possibility in the Mesa del Sol network. We note that the foregoing analysis of symmetries considers only the topology of the network. The dynamics of real electrical power networks is better described by a complex admittance matrix, for which the symmetries may be different.

Many other networks were studied for symmetries in refs 28,29 for the purpose of nding motifs and redundancies, but not

dynamics. Those networks were Human B Cell Genetic Interactions, Caenorabhditis elegans genetic interactions, BioGRID data sets (Human, Saccharomyces cerevisiae, Drosphila and Mus musculus), the internet (Autonomous Systems Level), and the US Power Grid. All the networks had many symmetries ranging in number from on the order of 1013 to 1011,298, and could be decomposed into many subgroups (from 3 to more

Tibet

Nepal

50 km

N

B =

India

Figure 6 | Geographical diagram of the Nepal power grid network. Colours are used to indicate the computed cluster structure. The matrix (inset) shows the structure of the diagonalized coupling matrix, analogous to Fig. 1a. The diagonal colours indicate which cluster is associated with each column.

6 NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2014 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5079 ARTICLE

Xg2KRg 7

where K is a conjugacy class, that is, aK(l) the character of that class for the lth IRR,

d(l) is the dimension of the lth IRR and is the order (size) of the group. Remark: the trivial representation (all IRR matrices 1 and a(l) 1) is always present and is

associated with the synchronization manifold. All other IRRs are associated with transverse directions. Next we use singular value decomposition on P(l) to nd the basis for the projection subspace for the lth IRR. Finally, we construct T by stacking the row basis vectors of all the IRRs, which will form an N N matrix.

We show the results of this applied to the 32-symmetry case in Fig. 1:

A

0 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 1 1 0 1 1 1 1 0 0

XKalK

Figure 7 | Network and cluster structure of the Mesa del Sol electric grid. Colours are used to denote clusters. Nodes coloured white are trivial clusters, containing only one element.

2

66666666666666664

3

77777777777777775

than 50). The subgroups were 90% or more made up of basic factors (not intertwined) consisting of various orders n of the symmetric group Sn. Hence, viewed as dynamical networks, all could show ID in the right situations.

DiscussionThe phenomena of symmetry-induced CS and ID appear to be possible in many model, man-made and natural networks, at least when modelled as unweighted couplings and identical systems. This is because our work here and that of refs 28,29 show that many types of networks can have symmetries and refs 20,21 show that many man-made and natural networks have dynamics similar to or reducible to equation (1) or its generalizations mentioned above. We have shown that ID is explained generally as a manifestation of clusters and subgroup decompositions. Furthermore, computational group theory can greatly aid in identifying CS in complex networks where symmetries are not obvious or far too numerous for visual identication. It also enables explanation of types of desynchronization patterns and transformation of dynamic equations into more tractable forms. This leads to an encompassing of or overlap with other phenomena, which are usually presented as separate. This list includes (1) remote synchronization13 in which nodes not directly connected by edges can synchronize (this is just a version of CS), (2) some types of chimera states27,37, which can appear when the number of trivial clusters is large and the number of non-trivial clusters is small, but the clusters are big (see ref. 38 for some simple examples), (3) partial synchronization where only part of the network is synchronized (shown for some special cases in ref. 39).

We note that although we have concentrated mostly on the maximal symmetry case, we can also examine the cases of lower symmetry induced by bifurcations that break the original symmetry and the same group theory techniques will apply to those cases. Some of this is developed for simple situations (rings or simple networks with point group symmetry) in ref. 15, but we now have the ability to extend this to arbitrary complex networks.

Methods

Calculating the IRR transformation. Below are the steps necessary to determine the symmetries of the network, obtain the clusters, nd the IRRs and, the most crucial part, calculate the transformation T from the node coordinates to the IRR coordinates that will block diagonalized A, as A commutes with all symmetries of the group40.

Using the discrete algebra software, it is straightforward to determine the group of symmetries of A, extract the orbits that give the nodes in each cluster and extract the permutation matrices Rg, and use the character table of the group and the traces of the Rgs to determine which IRRs are present in the node-space representation of the group. Remark: this step is discussed in any book on representations of nite groups (for example, ref. 22). After this, we put each Rg into its appropriate conjugacy class. To generate the transformaion T from the group information, we have written a code41 on top of the discrete algebra software, which for each IRR present constructs the projection operator P(l) (ref. 22) from the node coordinates onto the subspace of that IRR, where l indexes the set of IRRs present. This is done by calculating,

Pl

dl

h

8

The non-trivial clusters are the nodes [1, 8], [2, 3, 7, 9], [4, 6], [5, 10] (the numbering of nodes matches the row and column numbers of A). The transformation matrix is

T

2

66 66666666666664

0 0 0 0 0 0 0 0 0 0 2 0 0 0

2

p 0

2

p 0 0 0 0 0

2

p 0 0 0 0 0 0

2

p 0 0 0

0 0 0 0

2

p 0 0 0 0

2

p 0

1

2

0 1 1 0 0 0 1 0 1 0 0 2

p 0 0 0 0 0 0

2

p 0 0 0

0 0 0 0

2

p 0 0 0 0

2

p 0

3

77777777777777775 9

p 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 0

0 0

2

p 0 0 0 0 0

0 0 0

2

p 0

2

2

p 0 0 0 0

2

p 0 0

0

2

p 0 0 0 0

and the block diagonal coupling matrix is

B

2

66666666666666664

3

77777777777777775 10

0

2

p

2

p 0 2 0 0 0 0 0 0

2

p 1 2 2 2

p 0 0 0 0 0 0 0 2 1 1 2

2

p 0 0 0 0 0 0

2 2

2

p 2 1 1 2

2

2

p 0 0 0 0 0 0

2

p 2

2

p 2

2

p 2 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0

0 0 0 0 0 1 1 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 2 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Experimental methods. The SLM (Boulder Nonlinear Systems P512-1550) has an active area of 7.68 7.68 mm2 and a resolution of 512 512 pixels. The camera

(Goodrich SU320KTSW-1.7RT/RS170) has an indium gallium arsenide focal plane array with 320 256 pixels and an area of 8 6.4 mm2. Using a lens

assembly, the SLM was imaged onto the camera with the magnication and orientation adjusted so that each 2 2 pixel area on the SLM is projected onto a

camera pixel. The cameras frame rate was B8 Hz and was synchronized with the SLMs refresh rate.

NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

& 2014 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5079

Each oscillator in equation (5) corresponds to a square patch of 16 16 pixels

on the SLM, which is imaged onto an 8 8 pixel region of the cameras focal plane

array. The intensity I(x) is the average of camera pixel values in this area. The

same phase shift x is applied by each of the SLM pixels in the patch. We employ a spatial calibration and lookup table to compensate for inhomogeneities in the SLM.

Data analysis methods. The phase shifts imparted by each oscillator xtj are recorded and used to determine the degree of synchronization, as shown in Fig. 3.

An individual oscillators deviation from synchronization at a given time as shown in Fig. 3ce is measured by Dxti xti xt, where xt denotes a spatial average of

the phases of all of the nodes within a particular cluster at time t. To quantify the degree of synchronization within a cluster as shown in Fig. 3a, we calculate the

root-mean square synchronization error DxRMS xti x

t

2

1=2

where / ST

indicates an average over a time interval T (here taken to be 500 iterations).

References

1. Motter, A. E., Myers, S. A., Anghel, M. & Nishikawa, T. Spontaneous synchrony in power-grid networks. Nat. Phys. 9, 191197 (2013).

2. Pecora, L. M. & Carroll, T. L. Master stability functions for synchronized coupled systems. Phys. Rev. Lett. 80, 21092112 (1998).

3. Zhou, C. & Kurths, J. Hierarchical synchronization in complex networks with heterogeneous degrees. Chaos 16, 015104 (2006).

4. Do, A.-L., Hofener, J. & Gross, T. Engineering mesoscale structures with distinct dynamical implications. New J. Phys. 14, 115022 (2012).

5. Dahms, T., Lehnert, J. & Schooll, E. Cluster and group synchronization in delay-coupled networks. Phys. Rev. E 86, 016202 (2012).

6. Fu, C., Deng, Z., Huang, L. & Wang, X. Topological control of synchronous patterns in systems of networked chaotic oscillators. Phys. Rev. E 87, 032909 (2013).

7. Kanter, I., Zigzag, M., Englert, A., Geissler, F. & Kinzel, W. Synchronization of unidirectional time delay chaotic networks and the greatest common divisor. Europhys. Lett. 93, 6003 (2011).

8. Rosin, D. P., Rontani, D., Gauthier, D. J. & Scholl, E. Control of synchronization patterns in neural-like boolean networks. Phys. Rev. Lett. 110, 104102 (2013).

9. Sorrentino, F. & Ott, E. Network synchronization of groups. Phys. Rev. E 76, 056114 (2007).

10. Williams, C. R. S. et al. Experimental observations of group synchrony in a system of chaotic optoelectronic oscillators. Phys. Rev. Lett. 110, 064104 (2013).

11. Belykh, V. N., Osipov, G. V., Petrov, V. S., Suykens, J. A. K. & Vandewalle, J. Cluster synchronization in oscillatory networks. Chaos 18, 037106 (2008).12. DHuys, O., Vicente, R., Erneux, T., Danckaert, J. & Fischer, I. Synchronization properties of network motifs: Inuence of coupling delay and symmetry. Chaos 18, 037116 (2008).

13. Nicosia, V., Valencia, M., Chavez, M., Diaz-Guilera, A. & Latora, V. Remote synchronization reveals network symmetries and functional modules. Phys. Rev. Lett. 110, 174102 (2013).

14. Russo, G. & Slotine, J.-J. E. Symmetries, stability, and control in nonlinear systems and networks. Phys. Rev. E 84, 041929 (2011).

15. Golubitsky, M., Stewart, I. & Schaeffer, D. G. Singularities and Groups in Bifurcation Theory Vol. II (Springer-Verlag, New York, NY, 1985).

16. Golubitsky, M. & Stewart, I. The Symmetry Perspective: From Equilibrium to Chaos in Phase Space and Physical Space (Berkhauser-Verlag, Basel, 2002).

17. Judd, K. Networked dynamical systems with linear coupling: synchronisation patterns, coherence and other behaviours. Chaos 23, 043112 (2013).

18. Fink, K. S., Johnson, G., Carroll, T., Mar, D. & Pecora, L. Three-oscillator systems as universal probes of coupled oscillator stability. Phys. Rev. E 61, 50805090 (2000).

19. Irving, D. & Sorrentino, F. Synchronization of dynamical hypernetworks: dimensionality reduction through simultaneous block-diagonalization of matrices. Phys. Rev. E 86, 056102 (2012).

20. Arenas, A., Diaz-Guilera, A. A., Kurths, J., Moreno, Y. & Zhou, C. Synchronization in complex networks. Phys. Rep. 469, 93153 (2008).

21. Newman, M. Networks, An Introduction Ch. 18 (Oxford University Press, 2011).

22. Tinkham, M. Group Theory and Quantum Mechanics (McGraw-Hill, New York, NY, 1964).

23. Stein, William. SAGE: Software for Algebra an Geometry Experimentation (http://www.sagemath.org/sage/; http://sage.scipy.org/

Web End =http://www.sagemath.org/sage/; http://sage.scipy.org/ , 2013).

24. The GAP Group. GAP: Groups, Algorithms, and Programming, Version 4-4 (http://www.gap-system.org

Web End =http://www.gap-system.org, 2005).

25. Vallentin, F. Symmetry in semidenite programs,. Linear Algebra Appl. 430,

360369 (2009).

26. Goodman, R. & Wallach, N. R. Representations and Invariants of the Classical Groups (Cam-bridge University Press, 1998).

27. Hagerstrom, A. M. et al. Experimental observation of chimeras in coupled-map lattices. Nat. Phys. 8, 658661 (2012).

28. MacArthur, B. D. & Sanchez-Garcia, R. J. Spectral characteristics of network redundancy. Phys. Rev. E 80, 026117 (2009).

29. MacArthur, B. D., Sanchez-Garcia, R. J. & Anderson, J. W. On automorphism groups of networks. Discrete Appl. Math. 156, 35253531 (2008).

30. Dionne, B., Golubitsky, M. & Stewart, I. Coupled cells with internal symmetry:I. wreath products. Nonlinearity 9, 559574 (1996).31. Bollobas, B. & Riordan, O. M. Handbook of Graphs and Networks (eds Bornholdt, S. & Schuster, H. G.) (Wiley-VCH, Weinheim, 2003).

32. Barabasi, A.-L. & Albert, R. Emergence of scaling in random networks. Science 286, 509512 (1999).

33. Albert, R. & Barabasi, A.-L. Statistical mechanics of complex networks. Rev. Mod. Phys. 74, 4896 (2002).

34. Watts, Duncan J. & Strogatz, Steven H. Collective dynamics of smallworld networks. Nature 393, 440442 (1998).

35. Nepal electricity authority annual report 2011. available at http://www.nea.org.np

Web End =www.nea.org.np (2011).

36. Abdollahy, S., Lavrova, O., Mammoli, A., Willard, S. & Arellano, B. Innovative Smart Grid Technologies (ISGT), 2012 IEEE PES 16 (IEEE, 2012).

37. Abrams, D. M. & Strogatz, S. H. Chimera states for coupled oscillators. Phys. Rev. Lett. 93, 174102 (2004).

38. Laing, C. R. The dynamics of chimera states in heterogeneous Kuramoto networks. Phys. D 238, 15691588 (2009).

39. Belykh, V. N., Belykh, I. V. & Hasler, M. Hierarchy and stability of partially synchronous oscillations of diffusively coupled dynamical systems. Phys. Rev. E 62, 63326345 (2000).

40. Sagan, B. E. The Symmetric Group (Wadsworth Brooks, Pacic Grove, CA, 1991).

41. Hagerstrom, A. Network Symmetries and Synchronization (https://sourceforge.net/projects/networksym/

Web End =https:// https://sourceforge.net/projects/networksym/

Web End =sourceforge.net/projects/networksym/ , 2014).

Acknowledgements

We acknowledge help and guidance with computational group theory from David Joyner of the US Naval Academy, and information and computer code from Ben D. MacArthur and Ruben J. Sanchez-Garcia, both of the University of Southampton, to address the problem of nding group decompositions. We thank Vivek Bhandari for providing us a copy of the Nepal Electricity Authority Annual Report35. We thank Shahin Abdollahy and Andrea Mammoli for providing us data on the Mesa Del Sol electric network36.

Author contributions

L.M.P. conceived project and group theory analysis. F.S. analysed stability and performed simulations. A.M.H. performed experiments and simulations. T.E.M. and R.R. contributed to experimental design and guidance. All authors contributed to writing the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Pecora, L. M. et al. Cluster synchronization and isolated desynchronization in complex networks with symmetries. Nat. Commun. 5:4079 doi: 10.1038/ncomms5079 (2014).

T

8 NATURE COMMUNICATIONS | 5:4079 | DOI: 10.1038/ncomms5079 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2014 Macmillan Publishers Limited. All rights reserved.