Symmetry-breaking mechanism for the formation of cluster chimera patterns
Malbor Asllani, Bram A. Siebert, Alex Arenas, James P. Gleeson
SSymmetry-breaking mechanism for the formation of robust and stable chimerapatterns
Malbor Asllani , Bram A. Siebert , Alex Arenas , James P. Gleeson MACSI, Department of Mathematics and Statistics,University of Limerick, Limerick V94 T9PX, Ireland and Departament d’Enginyeria Inform`atica i Matem`atiques,Universitat Rovira i Virgili, 43007 Tarragona, Catalonia, Spain
Synchronization of the dynamics of the individual entities that constitute complex systems hasbeen observed and systematically studied for more than four decades now. An exotic behaviorthat has triggered scientists’ curiosity for a long time is the emergence of chimera states, long-lasting dynamical patterns of coexistence of synchronous and asynchronous clusters of nodes innetworked systems. Although many features have been understood so far, many questions remainunanswered. In particular, the stability of such states for finite networks of oscillators and their highsensitivity to the initial conditions have not yet been tackled exhaustively. Here we show that theformation of clustered synchronization can follow a global symmetry-breaking mechanism, which isalso responsible for selecting the nature of the final patterns. Through this novel method, we showthat it is possible to reconcile the pattern formation process with cluster synchronization, pavingthe way to a robust mechanism for explaining the formation of stable chimera states in networks ofcoupled dynamical systems.
The full understanding of the dynamics of complex sys-tems is probably one of the biggest challenge of physicsthis century. The emergence of collective behaviour andthe coexistence of order and disorder are fundamental keyplayers in this process. A very well–known phenomenonof collective behavior is that of synchronization, where aset of coupled oscillators spontaneously synchronize theirphases. This phenomenon has been widely studied in re-cent years and has been applied to explain the coordi-nation of the beats of the heart myocytes and the firingof brain neurons [1–5], the simultaneous flashing of malefireflies during the mating season [6] and the synchroniza-tion of the alternating current in power grids [7]. The-oretical studies have identified particular types of syn-chronization, among which it is worth mentioning the re-mote synchronization of two or more oscillators throughasynchronous ones [8], the synchronization of chaotic os-cillators and the phenomenon of cluster synchronizationwhen the set of oscillators is divided into subsets of syn-chronized ones [9] or the emergence of exotic states whereidentically coupled oscillators synchronise in clusters [10].The latter class also includes an essential and fasci-nating phenomenon of coexistence of order and disorderknown as chimera states [11–13]. Chimeras, named fromGreek mythology, are particular states where clusters ofcoherent and incoherent oscillators coexist. Since theirdiscovery, chimera states have aroused scientists’ curios-ity about the mechanisms that give birth to them, aboutthe effective lifespan of these states [ ? ], and their sen-sitivity to the initial conditions that make their experi-mental verification difficult to achieve [11, 14].Outstanding theoretical results have been obtained inthe study of chimeras, and promising methods have been proposed to reveal the mystery surrounding these exoticstates, such as whether they will exist, their stability,and specially, their origin. For example, it has beenfound that chimeras are stable in infinite-size networksand that they can be transient in finite ones [15]. So far,however, there is not an exhaustive answer to the sta-bility of chimera states [12, 16], or a robust mechanismthat explains how chimeras are generated [17]. An impor-tant advance in the area came from the group-theoreticalcharacterization of network symmetries that correspondsto the identification of partitions of nodes that can bestable or unstable [17–20]. Using this approach the au-thors are able to design the emergence of chimera states.Moreover, experimental results [21] prove the existence ofcluster synchronization and chimera states for the case ofcoupled Belousov-Zhabotinsky oscillators [22].Here, we aim at analyzing the chimera states, propos-ing an alternative path to the aforementioned cluster par-tition. The method we propose rigorously guarantees theformation of stable chimera patterns and simultaneouslyprocures a robust way for their emergence. More ex-plicitly, we will use a pattern formation approach basedon a global symmetry-breaking mechanism [22–24] thatyields (stable) cluster synchronization patterns as a re-sult of this procedure. Based on this principle, we willsuccessively generate stable chimera oscillations that arerobust to initial conditions. The great advantage of thismethod is that it has been extensively studied startingfrom the early 1970’s, and many groundbreaking rigorousanalytical results have been obtained regarding the pre-diction of the shape and stability of the final nonlinearpatterns. Furthermore, the pattern formation process isrobust to the choice of the initial conditions. These threeaspects are crucial because they will pave the way andfacilitate the process of obtaining stable chimera states. a r X i v : . [ n li n . AO ] F e b This way, instead of tackling the problem frontally tounderstand how to decompose the chimera in groups oflocally stable and unstable nodes, we follow a shortcutby rigorously proving the emergence of simultaneouslycoherent and incoherent oscillatory nodes.We start by considering the general setting of a net-work of coupled dynamical systems through the equa-tions ˙ x i = F ( x i ( t )) + σ (cid:12) N X i =1 W ij G ( x j ( t )) , (1)where each of the N individuals, identified as nodes ofthe network, are represented by M -dimensional vectorstate variables x i [3]. The local dynamics of the i -thnode is given by the nonlinear function F ( x i ), σ is thevector of the coupling strength, W ij are the weightedentries of the adjacency matrix and G ( x j ) is the cou-pling function. Notice that in the above equation, (cid:12) represents the element-wise product, and we have gener-alised the global coupling term σ into a vector, a gen-eralisation that will prove useful in the following. Westart by considering the same synchronised equilibriumstate for all the nodes x ∗ ( t ), which is stable in the ab-sence of coupling. Then by perturbing this state by asmall random term x i ( t ) = x ∗ ( t ) + ξ i ( t ) we can obtainat the first perturbative order the variational equations˙ ξ i = DF ( x ∗ ( t )) ξ i ( t ) + σ (cid:12) P Ni =1 W ij DG ( x ∗ ( t )) ξ j ( t ).Now starting from here and assuming that the linearizedspatial operator DG ( x ∗ ( t )) is diagonalizable, we canproceed further to obtain the Master Stability Function(MSF) [3, 25], which relates the stability of the uniformlysynchronized state x ∗ to the topology of the connectionnetwork. So, assuming that there is a linearly indepen-dent basis of eigenfunctions(-vectors) of the spatial oper-ator, we can obtain the diagonalized decoupled equations˙ ζ α = (cid:0) DF ( x ∗ ( t )) α + Λ ( α ) σ (cid:1) (cid:12) ζ α ( t ) (2)where by Λ ( α ) we have denoted the α -th eigenvalue of thespatial operator. To concretize our analysis and withoutloss of generality, throughout this text, we will consider adiffusive coupling through the Laplacian matrix L whoseentries are L ij = W ij − k i δ ij with k i = P Nj =1 W ij .The formalism described so far allows the local analy-sis of the linear stability of the originally uniform state.Depending on the nature of the latter, eq. (2) might benon-autonomous when the uniform state is a limit cycleor autonomous when it is a fixed point. The traditionalapproach in the study of synchronization phenomenon isbased on the former case when the system is already in asynchronized limit cycle manifold. In this case, the aimis to keep as stable as possible such a state since an even-tual instability would cause the desynchronization of thewhole system. As already anticipated, we will adopt herean alternative approach that to the best of our knowledge has not been used so far in the vast literature related tothe synchronization dynamics. The first major differenceis that we start from an unstable uniform equilibrium of afixed point, so x ∗ i ( t ) = x ∗ for all the nodes. Starting froma uniform stationary state will allow the use of the spec-tral analysis instead of the classical Lyapunov exponentused for non-autonomous systems [3, 25]. This makespossible the necessary analytical treatment for findingthe instability conditions and the description of the or-bits’ behavior in the initial linear regime. Since a homo-geneous stationary equilibrium lacks the oscillatory be-haviour that characterises synchronization dynamics, wewill recover final time-varying patterns by consideringa nonlinear function constituted by 3 variables, knownin the mathematical biology literature as species [24] F ( x i ) = (cid:2) f ( x i ) , g ( y i ) , h ( z i ) (cid:3) . For each node we writethe 3-dimensional state vector x i = [ u i , v i , z i ] with a littleabuse of notation. In other words, we have ˙ u i = f ( u i , v i , z i ) + σ u P Nj =1 L ij u j ˙ v i = g ( u i , v i , z i ) + σ v P Nj =1 L ij v j ˙ z i = h ( u i , v i , z i ) + σ z P Nj =1 L ij z j . (3)It is well known that a 3-species system is a sufficientrequirement to have an oscillatory instability and, as aconsequence, nonlinear oscillatory patterns [23, 24]. Thisis the key factor in producing the global oscillations ofour system, as is shown in Fig. 1 panels a ) − a ). Asan alternative to this model, oscillatory patterns can beobtained also for a 2-species model in a directed net-work where the oscillations are as a consequence of thecomplex spectrum of the Laplacian matrix [26, 27] (seethe Supplemental Material (SM) for a detailed descrip-tion). We should, however, emphasise that in general, ina complex network, the oscillations appear to be asyn-chronous following the shape of a traveling front in thelimit case of a regular lattice [23, 28] (see Fig. 1 c ) and d )). To understand the rationale behind such behavior,we should look for the shape of the eigenvectors, whichdescribe the linear evolution of the orbits in a short timeand also determine the general shape of the final nonlin-ear pattern. In fact in order to solve the linear (decou-pled) variational equations (2) we look for solutions ofthe form ζ i = P Nα =1 c α e γ α Φ ( α ) i where γ α and Φ ( α ) rep-resent the α -th growth mode and Laplacian eigenvectorrespectively. So, in the initial regime, once there is atleast one value of α for which the real part of γ α , knownas the dispersion relation, is positive, as is the case forFig. 1 a ), then the entries of vectors of the 3 species willbe scaled by the entries of the eigenvector Φ ( α ) . Also, in a3-species model, the imaginary part is different from zeroallowing an oscillatory instability to occur [23]. In regularmulti-dimensional lattices, such an eigenvector representsa discrete term of multi-dimensional Fourier series. Thisis why in such discrete domains, the final patterns will bea combination, depending on the number of the unstable(Fourier) modes, of time-varying sinusoidal shaped (trav-eling waves) [23, 24, 28]. Furthermore, such final patternsare not transient, but it can be shown that there existsa stability region (known also as “stability balloon”) [23]of parameters for which they are asymptotically stable.Having briefly introduced the mechanism responsiblefor shaping the stable nonlinear pattern, we next describehow to construct the coupling network to obtain clus-ters of coherent and incoherent nodes [ ? ]. We startconsidering the set of eigenvectors (cid:2) Φ (0) , Φ (1) , . . . , Φ ( N ) (cid:3) of the Laplacian of a connected and undirected networksorted in decreasing order of the respective modes (eigen-values). From the algebraic connectivity, we know thatthe network Laplacian will have a single zero eigenvalueΛ (0) = 0 whose eigenvector is uniformly distributed,for simplicity here assumed to be the scaled unit vec-tor Φ (0) = a = [ a, a, . . . , a ] T where a is the scalingconstant. Let us remark here that from eq. (2) thenull mode Λ (0) = 0 represents the a-spatial case mean-ing that the nodes should be stable once uncoupled. Sowe should look for the eigenvectors corresponding to thenonnull modes to destabilize our system globally. To ob-tain eigenvectors with clusters of entries that are close toeach other, in total analogy with the equilibrium state,we proceed by construction. We start by considering net-works with more than one null mode, for example, C sucheigenvalues Λ (0) = · · · = C Λ (0) = 0. This means thatour network has C disconnected components, here refer-ring to as disconnected clusters. The Laplacian matrixand its eigenvectors in this case are respectively a blockmatrix and block vectors L = diag [ L , . . . , n L , . . . , C L ]and Φ ( α ) = (cid:2) Φ ( α ) , . . . , n Φ ( α ) , . . . , C Φ ( α ) (cid:3) . Here, we havedivided the N nodes into C disconnected clusters withthe same number of nodes. Observe that the eigen-vector of the n -th null mode can include different so-lutions, but here we chose purposely the one that hasconstant entries for the blocks of each cluster n Φ (0) = (cid:2) n a . . . , n a | {z } N/C . . . , n a i . . . , n a i | {z } N/C . . . n a C . . . , n a C | {z } N/C (cid:3) . Insteadfor the eigenvectors of the non null modes the only possi-ble solution is n Φ ( α =0) = (cid:2) . . . , | {z } N/C . . . , n Φ ( α ) | {z } N/C . . . . . . , | {z } N/C (cid:3) where n Φ ( α ) is the α -th eigenvector of the n -th discon-nected cluster. We invite the interested reader to referto the Supplemental Material for a rigorous derivationof this assertion. For the next step, we weakly connectthe C disconnected clusters through a minimal numberof links. From a purely perturbative point of view, thismeans that the former spectrum will be weakly deformed.The first indication of this deformation is that due to theconnectedness, the graph Laplacian will have only onezero eigenvalue. All the remaining C − C − C − C − a )is known as the Fiedler eigenvalue, and it is an indicatorof the algebraic connectivity of the network [31]. Besidesthe modular eigenvalues, we also have the non-modularones, the remaining N − C eigenvalues. Following theperturbative analysis above, explained in more detail inthe SM, it is not possible that their eigenvectors havesegregated nonnull entries per module. This means thatif a non-modular eigenvalue would be the only unstablemode then the pattern of such setting would be irregularfor a given module and almost not expressed for the restof them. The linear stability analysis allows the exten-sion of pattern selection beyond the single unstable modecase. In fact, if two or more modes are unstable and nearthe bifurcation threshold, they will compete with eachother in the linear regime until the nonlinearites stabilisethe pattern shaping. It is reasonable to expect that thepattern outcome will reflect the shape of the eigenvec-tors of the unstable modes. To show this behaviour, inFig. 2 we have taken into consideration only two unsta-ble eigenvalues, a modular and a non-modular one. Thechoice has been such that in the dispersion relation rep-resentation, the two modes have comparable real parts.In this way we expect that the two different behaviourswill be present to similar levels. In fact, as can be seenfrom the pattern snapshot taken at a given time, Fig. 2 b ), the second module is non-homogeneous in the con-centration (of the first species in this case), and clearlyfrom panel c ) this state constitutes a chimera. Althoughthe results presented is both figures has been validated a ) a ) a )d ) d )
25 50 75 100 12500.20.40.60.81
EigenvectorPattern
III IV VIII d )
100 200 300 400 500255075100125
IIIIIIIVV
25 50 75 100 12500.20.40.60.81
EigenvectorPattern
I II III IV V b )
100 200 300 400 500255075100125
IIIIIIIVV c )
100 200 300 400 500255075100125
25 50 75 100 12500.20.40.60.81
EigenvectorPattern b )c )b )c ) FIG. 1.
Cluster synchronisation in modular networks.
The 3-species model used to generate the patterns follows thatof Zhabotinsky [29, 30] defined as f ( u, v, z ) = − c uv + c − c u/ ( q + u ) with q = 0 . g ( u, v, z ) = c uv − c v + c , and h ( u, v, z ) = c u − c v with general parameters c = c = 28 . c = c = 15 . c = c = 1, c = 25 . c = 3 . σ u = σ v = 0and only σ z is varying according to the specific case. In the columns from left to right are shown the networks in increasingmodularity, from a ring lattice (left) with σ z = 30 to a strong modular one (right) σ z = 59 passing through a intermediatemodularity (center) σ z = 4 . a ) − a ) The dispersion relation for the different setting of clusterisation of the network wherea single unstable mode has been selected. In the inset are shown both the real and the imaginary part of the dispersion relation. b ) − b ) Comparison between the eigenvectors of the corresponding unstable modes vs. the final pattern. Notice the changein the shape of the eigenvector from a Fourier (discretised) eigenfunction to a cluster segregated one. c ) − c ) Time evolutionof patterns of the first species, from a travelling wave in the ring lattice to a cluster synchronised one for the strong modularnetwork. Also the small size of the pattern is due to the tiny value of real part of the unstable eigenvalue, purposely selectednear the bifurcation threshold. (An example of a larger pattern can be found in SI.) d ) − d ) Graphical representation ofsnapshots of the oscillatory patterns. for the case of modules of the same size, in the SM we have shown that our method is robust even for the caseof networks with modules of considerably different sizes.We note that stationary Turing patterns [32] has beenrecently studied also in the context of modular networks[33].To conclude, in this paper we have studied the problemof emergence of synchronised and robust chimera statesthat follow the modular structure of the underlying net-work. To understand the mechanism that generates thesestates we have developed a theoretical framework basedon symmetry-breaking principles. Due to a global in-stability the system departs from an homogeneous stateand the perturbations associated to each node will beshaped according to the eigenvectors corresponding tothe unstable modes of the Laplacian matrix. This lin-ear regime behaviour is saturated by the nonlinear termsin the final stage, inheriting the eigenvector characteris-tics in the spatially extended pattern. In order to obtainstable chimera states we start with a setting of synchro-nised clusters. We show that is possible to reconstructnetworked structures that have particular eigenvectorswith certain features: they are organised in groups ofentries close to each other. Based in a perturbative ap-proach, we prove that these networks have C − C isthe number of their clusters. On the other hand the re- maining N − C non equilibrium eigenstates, here callednon-modular ones, present irregular shapes. At this pointwe could select (at least) two unstables eigenvalues, onemodular and one non-modular in the Master StabilityFunction such that due to their competition during thepattern formation process the final result was a chimerastate. The pattern formation mechanism presents ma-jor advantages, first, it assures the stability of the finalchimera states, and second, the generation process is ro-bust to the choice of initial conditions. These two aspectsare the main weak points of the classical methods for gen-erating chimera states. In this work we have particularlystudied the case when the coupled dynamical systemsare fixed points instead of the usual nonlinear oscillators.The main reason for that was because in this way it waspossible to make use of the solid classical results alreadyexisting in the theory of pattern formation. Neverthe-less we believe that the results shown here can be easilyextended also to the case of coupled oscillators. We areconfident that our results will shed light on a better un-derstanding of the generation mechanisms and stabilityof chimera states and also fill the theoretical gap of amathematical explanation for recent experimental obser-vations in this domain [21]. Acknowledgements:
B.A.S. acknowledges fundingfrom the Irish Research Council under Grant No.GOIPG/2018/3026. The work of J.P.G. and M.A. ispartly funded by Science Foundation Ireland (Grants No.16/IA/4470, No. 16/RC/3918, No. 12/RC/2289 P2, andNo. 18/CRT/6049) and cofunded under the EuropeanRegional Development Fund. A.A. acknowledges finan-cial support from Spanish MINECO (grant PGC2018-094754-B-C21), Generalitat de Catalunya (grant No.2017SGR-896), Universitat Rovira i Virgili (grantNo. 2019PFR-URV-B2-41), Generalitat de CatalunyaICREA Academia, and the James S. McDonnell Foun-dation (grant [1] A. Pikovsky, M. Rosenblum, J. Kurths,
Synchronization:A universal concept in nonlinear sciences , (CambridgeUniversity Press, 2001).[2] S. H. Strogatz & I. Stewart, Coupled oscillators and bio-logical synchronization.
Sci. Am. (6), 102–109 (1993).[3] A. Arenas, A. A. Diaz-Guilera, J. Kurths, Y. Moreno, C.Zhou, Synchronization in complex networks.
Phys. Rep. , 93–153 (2008).[4] D. A. Wiley, S. H. Strogatz, and M. Girvan, The size ofthe sync basin.
Chaos , 015103 (2006).[5] T. Chouzouris, I. Omelchenko, A. Zakharova, J. Hlinka,P. Jiruska, E. Sch¨oll, Chimera states in brain networks:Empirical neural vs. modular fractal connectivity. Chaos , 045112 (2018).[6] J. Buck & E. Buck, Synchronous fireflies. Sci. Am. (5), 74–85 (1976).[7] A. E. Motter, S. A. Myers, M. Anghel, T. Nishikawa, Spontaneous synchrony in power-grid networks.
Nat.Phys. , 191–197 (2013).[8] V. Nicosia, M. Valencia, M. Chavez, A. D´ıaz-Guilera, andV. Latora, Remote synchronization reveals network sym-metries and functional modules. Phys. Rev. Lett. ,174102 (2013).[9] V. N. Belykh, G. V. Osipov, V. S. Petrov, J. A. K.Suykens and J. Vandewalle, Cluster synchronization inoscillatory networks.
Chaos , 037106 (2008).[10] M. H. Matheny, et. al. , Exotic states in a simple networkof nanoelectromechanical oscillators. Science , 1057(2019).[11] Y. Kuramoto & D. Battogtokh, Coexistence of coherenceand incoherence in nonlocally coupled phase oscillators.
Nonlinear Phenom. Complex Syst. , 380–5 (2002).[12] D. M. Abrams & S. H. Strogatz, Chimera states for cou-pled oscillators. Phys. Rev. Lett. , 174102 (2004).[13] E. Sch¨oll, A. Zakharova, R. G. Andrzejak, Chimera statesin complex networks. Front. Appl. Math. Stat , 62(2019).[14] M. J. Panaggio & D. M. Abrams, Chimera states: co-existence of coherence and incoherence in networks ofcoupled oscillators. Nonlinearity , R67–R87 (2015).[15] M. Wolfrum & O. E. Omel’chenko, Chimera states arechaotic transients. Phys. Rev. E , 015201(R) (2011).[16] D. M. Abrams, R. Mirollo, S. H. Strogatz, D. A. Wiley,Solvable model for chimera states of coupled oscillators. Phys. Rev. Lett. , 084103 (2008).[17] Y. S. Cho, T. Nishikawa, A. E. Motter, Stable chimerasand independently synchronizable clusters.
Phys. Rev.Lett. , 084101 (2017).[18] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E.Murphy, R. Roy Cluster synchronization and isolateddesynchronization in complex networks with symmetries.
Nat. Commun. , 4079 (2014).
100 200 300 400 500255075100125
IIIIIIVVI
25 50 75 100 125-1-0.500.51
Mod EigenvectorNon-Mod EigenvectorPattern
I II III IV V a)c) c )b)d) FIG. 2.
Chimera states in modular networks a)
The dispersion relation shows that in order to obtain chimera states,(at least) two modes a modular and a non-modular should be simultaneously unstable. In the inset are shown both the realand the imaginary part of the dispersion relation. b) A comparison between the eigenvectors of the selected unstable modular(red stars) and non-modular (green squares) eigenvalues are shown here in comparison with the final pattern (blue circles)taken at a given fixed time. Notice how the irregularity of the non-modular eigenvector at the second module influences alsothe final pattern. c) The chimera states evolution of the first species is shown as a function of time where the second moduleis incoherent compared to all the remaining four ones. d) A graphical representation of the snapshots of the chimera statesin the 5 modules network. The model is the same as in Fig. 1 and the parameters used in this case are c = c = 27 . c = c = 15 . c = c = 1 c = 25 . c = 3 . σ u = σ v = 0 and σ z = 5 . Sci. Adv. , e1501737 (2016).[20] F. D. Rossa, L. M. Pecora, K. Blaha, A. Shirin, I. Klick-stein, F. Sorrentino, Symmetries and cluster synchro-nization in multilayer networks. Nat. Commun. , 1–17(2020).[21] M. R. Tinsley, S. Nkomo, T. L. Showalter, Chimera andphase-cluster states in populations of coupled chemicaloscillators. Nat. Phys. , 662–665 (2012).[22] G. Nicolis & I. Prigogine, Self-organization in nonequi-librium systems: From dissipative structures to orderthrough fluctuations , (J. Wiley & Sons, 1977).[23] M. Cross & H. Greenside,
Pattern formation and dy-namics in nonequilibrium systems , (Cambridge Univer-sity Press, 2009).[24] J. D. Murray,
Mathematical biology II: Spatial models andbiomedical applications , (Springer-Verlag, 2001).[25] L. M. Pecora, T. L. Carroll, Master stability functionsfor synchronized coupled systems.
Phys. Rev. Lett. ,2109–2112 (1998). [26] M. Asllani, J. D. Challenger, F. S. Pavone, L. Sacconi,D. Fanelli, The theory of pattern formation on directednetworks. Nat. Commun. (1), 1–9 (2014).[27] M. Asllani, T. Carletti, D. Fanelli, P. K. Maini, A univer-sal route to pattern formation in multicellular systems. Eur. Phys. J. B , 135 (2020).[28] H. G. Othmer, L. E. Scriven, Instability and dynamicpattern in cellular networks. J. Theor. Biol. (3), 507–537 (1971).[29] A. M. Zhabotinsky, M. Dolnik, I. R. Epstein, Pattern for-mation arising from wave instability in a simple reaction-diffusion system. J. Chem. Phys. , 10306 (1995).[30] M. Asllani, T. Biancalani, D. Fanelli, A. J. McKane, Thelinear noise approximation for reaction-diffusion systemson networks.
Eur. Phys. J. B , 476 (2013).[31] M. E. J. Newman, Networks , (Oxford University Press,2018).[32] A. M. Turing, The chemical basis of morphogenesis.
Phil.Trans. R. Soc. Lond. B , 37 (1952).[33] B. Siebert, C. L. Hall, J. P. Gleeson, M. Asllani, Roleof modularity in self-organization dynamics in biologicalnetworks.
Phys. Rev. E , 052306 (2020). upplemental MaterialSymmetry-breaking mechanism for the formation of robust and stable chimerapatterns
Malbor Asllani , Bram A. Siebert , Alex Arenas , James P. Gleeson MACSI, Department of Mathematics and Statistics,University of Limerick, Limerick V94 T9PX, Ireland and Departament d’Enginyeria Inform`atica i Matem`atiques,Universitat Rovira i Virgili, 43007 Tarragona, Catalonia, Spain
I. LAPLACIAN EIGENVECTORS OF C WEAKLY CONNECTED CLUSTERS NETWORK
Starting from a perturbative setting is possibleto show that the eigenvector of the n -th nullmode can select only one of the possible solu-tions, the one that has constant entries of alternat-ing signs for the blocks of each cluster n Φ (0) = (cid:2) n a . . . , n a | {z } N/C . . . , n a i . . . , n a i | {z } N/C . . . n a C . . . , n a C | {z } N/C (cid:3) T . In fact,the most intuitive eigenvector of the null mode is thecanonical one n Φ (0) = (cid:2) . . . | {z } N/C . . . , n . . . n | {z } N/C . . . , . . . | {z } N/C (cid:3) T .However, this choice is not possible since it does not re-spect the perturbative approach once the clusters areweakly connected. To prove that we consider the sim-ple case of 2 clusters, but the same easily follows for thecase of C clusters. The perturbative eigenvalues problemfor the canonical eigenvector is now written as L z }| {(cid:20) L E E
21 2 L (cid:21) (cid:20) + (cid:15) (cid:15) (cid:21) = λ (cid:15) (cid:20) + (cid:15) (cid:15) (cid:21) . where E , E are the intra-cluster weak coupling, (cid:15) , (cid:15) the perturbation of the uncoupled eigenvectors and λ (cid:15) the pertubation of the null mode. From this relation weobtain the following equalities ( L (cid:15) + E (cid:15) = λ (cid:15) + λ (cid:15) (cid:15) : O ( (cid:15) ) = O ( (cid:15) ) E + E (cid:15) + L (cid:15) = λ (cid:15) (cid:15) : O ( (cid:15) ) ? = O ( (cid:15) ) . Clearly since the last relation is not possible this choicefor the eigenvector is not permitted.On the contrary, for the eigenvectors of the nonnull modes the only possible solution is n Φ ( α =0) = (cid:2) . . . , | {z } N/C . . . , n Φ ( α ) | {z } N/C . . . . . . , | {z } N/C (cid:3) T where n Φ ( α ) is the α -theigenvector of the n -th disconnected cluster. In fact, if forthe 2 cluster case we consider as candidate (cid:2) n Φ ( α ) , (cid:3) T from the same reasoning as above we have L z }| {(cid:20) L E E
21 2 L (cid:21) (cid:20) n Φ ( α ) + (cid:15) + (cid:15) (cid:21) = (cid:16) Λ ( α ) + λ (cid:15) (cid:17) (cid:20) n Φ ( α ) + (cid:15) + (cid:15) (cid:21) , so ( L (cid:15) + E + E (cid:15) = λ (cid:15)n Φ ( α ) + Λ ( α ) (cid:15) + λ (cid:15) (cid:15) E n Φ ( α ) + E (cid:15) + L (cid:15) = Λ ( α ) + Λ ( α ) (cid:15) + λ (cid:15) + λ (cid:15) (cid:15) . This is equivalent to ( O ( (cid:15) ) = O ( (cid:15) ) O ( (cid:15) ) ? = O ( (cid:15) )which again cannot be true. This concludes that theonly possible solution for the Laplacian eigenvectors of anetwork with disconnected clusters is defined as above. II. CLUSTER SYNCHRONISED ANDCHIMERA STATES IN MODULAR DIRECTEDNETWORKS
In this section we will consider the case when the clus-ter synchronisation and chimera states emerge due to themodularity of directed networks. A major contribute ofdirected networks is due to the fact that the spectrum ofthe Laplacian matrix is complex which cause a Turing-Hopf instability in the linear regime yielding this wayoscillatory patterns. Let us remind here that such pat-terns are not otherwise possible for a 2-species reaction-diffusion system. Also the directedness of the networkedstructure contributes in the increment of chances foremergence of patterns due to a higher dispersion rela-tion compared to the symmetric network or the contin-uous domain. For further details the interested readershould refer to [1, 2]. It is important also to remind that,in general, a directed networked structure might be alsostrongly non-normal as observed in many real-world net-works [3]. In this case the pattern formation process fol-lows different mechanisms induced by the non-normalityof the Laplacian matrix [4]. However, we decided to skipsuch scenario that goes beyond the scope of this paperand focus instead to the strictly normal case. a r X i v : . [ n li n . AO ] F e b We start by considering a 2-species activator-inhibitormodel of a reaction-diffusion system. The concentrationfor the two species are denoted u i , for the activator and v i , for the inhibitor and where the index i refers to oneof the N cells. The corresponding equations take thefollowing form: ˙ u i = f ( u i , v i ) + D u P Nj =1 L ij u j ˙ v i = g ( u i , v i ) + D v P Nj =1 L ij v j where D u and D v are the diffusion constants and f ( · , · ), g ( · , · ) are the two nonlinear functions representing the in-teraction of the species inside each node i . The spatialinteractions on the other side, are represented by the dif-fusion operator with entries L ij which are a reflection ofthe topology of the networked support.To proceed with the stability analysis, we linearise thereaction-diffusion system (1) around the uniform steadystate ( u ∗ , v ∗ ): (cid:18) ˙ δ u ˙ δ v (cid:19) = (cid:18) f u I + D u L f v I g u I g v + D v L (cid:19) (cid:18) δ u δ v (cid:19) = (cid:0) J + D (cid:1) (cid:18) δ u δ v (cid:19) where here δ u , δ v are the vectors of the perturbations.The Jacobian matrix for the reaction terms, evaluated atthe equilibrium point is given by J = (cid:18) f u I f v I g u I g v I (cid:19) , andthe diffusion one for both species is D = (cid:18) D u L OO D v L (cid:19) ,where O is the N × N zero-valued matrix. We considerthe basis of the eigenvectors Φ αi , with α = 1 , . . . , N , ofthe Laplacian operator L to the corresponding eigenval-ues Λ α . Due to the asymmetry of the adjacency matrix A , the spectrum Λ α is, in principle, complex. In particu-lar is has been shown in [1, 2] that for balanced directednetworks where the numbers of incoming edges is equalto the outgoing ones k outi = k ini , such a basis alwaysexists. This will be also the case in the following discus-sion. In order to solve this linearised system, we expandthe perturbations on the basis of the eigenvectors, i.e. δu i = P Nα =1 b α Φ αi and δv i = P Nα =1 c α Φ αi . Followingthis, it is possible to reduces the 2 N × N system to a2 × α = 1 , . . . , N :det (cid:0) J α − λ α I (cid:1) = det (cid:18) f u + D u Λ α − λ α f v g u g v + D v Λ α − λ α (cid:19) = 0where J α is the 2 × J α ≡ J + DΛ α where D = diag ( D u , D v ), Λ α = diag (Λ α , Λ α ). Ifthe real part of λ α , known as the dispersion relation, hasat least a mode Λ α =0 for which takes positive values, thenthe steady state becomes unstable. Mathematically thisis formalised as: λ α = 12 [(tr J α ) Re + γ ] + 12 [(tr J α ) Im + δ ] ι (1) where γ = s A + √ A + B δ = sgn( B ) s − A + √ A + B A = [(tr J α ) Re ] − [(tr J α ) Im ] − [(det J α ) Re ] B = 2 (tr J α ) Re (tr J α ) Im − [(det J α ) Im ] . In the above expression we have denoted by sgn( · ) thesign function and () Re , () Im indicate, the real and imagi-nary parts of the respective function. To proceed furtherwith the calculations, we make use of the definition of asquare root of a complex number. By taking z = a + bι where ι = √− √ z = ± r a + | z | b ) r − a + | z | ι ! . The instability sets in when | (tr J α ) Re | ≤ γ , a conditionthat can be expressed by the following inequality [1]: S (Λ ( α ) Re )[Λ ( α ) Im ] ≤ − S (Λ ( α ) Re ) , (2)where (Λ Re , Λ Im ) span the complex plane where lie theeigenvalues of the Laplacian operator. The explicit formof the polynomials S , S is described by: S ( x ) = C x + C x + C x + C x + C S ( x ) = C x + C x + C . The above constants are given by: C = σ (1 + σ ) C = (1 + σ ) ( σJ + J ) + 2tr J σ (1 + σ ) C = det J (1 + σ ) + (tr J ) σ + 2tr J (1 + σ ) ( σJ + J ) C = 2tr J (1 + σ ) det J + (tr J ) ( σJ + J ) C = det J (tr J ) C = σ (1 − σ ) C = ( σJ + J ) (1 − σ ) C = J J (1 − σ ) . The above system undergoes a Turing-Hopf instability,induced by the (directed) topology of the interaction net-work. Such mechanism responsible for the emergence ofpattern in directed networks, is known as the topology-driven instability [1]. In Fig. 1 we show the cluster syn-chronisation phenomenon in a directed modular networkwhere the pattern are triggered by the directedness of theunderlying network and the oscillations are as a conse-quence of a complex spectrum of the Laplacian operator.If apart the unstable modular mode, we could chose alsoto destabilise one of the non-modular ones then chimerastates would appear in this case as shown in Fig. 2.
III. ROBUSTNESS OF THESYMMETRY-BREAKING METHOD FORMODULES OF DIFFERENT SIZES
The method that we proposed here has been shown toyield excellent results in the generation and predictionof the cluster synchronization on the one hand and thechimera states on the other. However, in the results pre-sented in the main text and the SM, we have considerednetworks of the same size. This section will discuss theeffectiveness of the symmetry-breaking mechanism whenthe modules’ size is considerably different. We start byconsidering a simple network model of
N W topology of200 nodes divided into two modules with a number ofinter-edges between the modules few on average com-pared to the intra-edges ones. Keeping a small inter vs.intra number of edges is a necessary requisite for the va-lidity of the perturbative approach used to justify ourmethod. Initially, the two modules have 50 nodes each ofas shown in panels a ) − a ) in Fig. 3, and we chose theparameters in such a way to obtain a weak chimera statewhere module I is weakly incoherent due to a mix of un-stable modular and non-modular eigenvectors. Clearly,in this setting, the results are shown throughout this work stand firmly, as confirmed by the clear gap between themodular and non-modular set of eigenvalues, panel a ).Then we tune the network topology by keeping the sametotal number of links but changing the size of the mod-ules in the 150 vs. 50 ratio. We notice that the firsttwo non-zero eigenvalues’ distribution slightly shifts tothe right, panel a ), so that the modular eigenvalue in-creases the spectral gap from the origin. When we raisethe ratio to 190 vs. 10 nodes per module, we can observea neat shift of the modular eigenvalue and a decreaseof the gap with the non-modular one, panel a . Suchbehavior is understandable from the spectral perturba-tive perspective. In fact, when we have smaller modules,i.e., of 10 nodes as in our case, the action of the unitaryweighted interconnection links (that remain invariant innumber in all the three examples) over the small mod-ule can be considered less perturbative compared to theprevious scenarios when such module was of 100 nodes.Nevertheless, the results of Fig. 3, in particular panels a ), b ), and c ) show that the shape of the eigenvectorsremains still sufficiently segregated to yield an explicitcluster synchronization making it an obvious indicator ofthe robustness of our method. [1] M. Asllani, J. D. Challenger, F. S. Pavone, L. Sacconi,D. Fanelli, The theory of pattern formation on directednetworks. Nat. Commun. (1), 1–9 (2014).[2] M. Asllani, T. Carletti, D. Fanelli, P. K. Maini, A universalroute to pattern formation in multicellular systems. Eur.Phys. J. B , 135 (2020).[3] M. Asllani, R. Lambiotte, T. Carletti, Structure and dy-namical behavior of non-normal networks. Sci. Adv. (12), eaau9403 (2018).[4] R. Muolo, M. Asllani, D. Fanelli, P. K. Maini, T. Carletti,Patterns of non-normality in networked systems. J. Theor.Biol. , 81–91 (2019).[5] G. Nicolis and I. Prigogine,
Self-organization in nonequi-librium systems: From dissipative structures to orderthrough fluctuations , (J. Wiley & Sons, 1977).
ContinuousRealImaginary
200 400 600 800 100020406080100120
200 400 600 800 100020406080100120 c ) c ) c )d ) d ) d ) EvectorPattern
ContinuousRealImaginary a ) a ) a )b ) b ) b ) FIG. 1.
Cluster synchronisation in directed modular networks
The 3-species model used to generate the patternsfollows that of Brusselator f ( u, v ) = 1 − ( b + 1) u + cu v and g ( u, v ) = bu − cu v [5] where the parameters . . . are fixed andand only σ z is varying according to the specific case. In the columns from left to right are shown the networks in increasingmodularity, from a ring lattice (left) with σ z = 30 to a strong modular one (right) σ z = 59 passing through a intermediatemodularity (center) σ z = 4 . (a) The dispersion relation for the different setting of clusterisation of the network where asingle unstable mode has been selected. In the inset are shown both the real and the imaginary part of the dispersion relation.Let us notice that for a directed network the (discrete) dispersion relation is always higher then the continuous counterpart. (b)
Comparison between the eigenvectors of the corresponding unstable modes vs. the final pattern. Notice the change in theshape of the eigenvector from a Fourier (discretised) eigenfunction to a cluster segregated one. (c)
Patterns evolution, from atravelling wave in the ring lattice to a cluster synchronised one for the strong modular network. (d)
Graphical representationof snapshots of the oscillatory patterns.
200 400 600 800 100020406080100120 a) b)
FIG. 2.
Chimera states in directed modular networks
We have used the Brusselator model with parameters a = 1 . b = 3 . c = d = 1, D u = 0 . D v = 3 . a) The dispersion relation shows that a modular and a non-modular eigenvalueare simultaneously unstable. In the inset are given the real and the imaginary part of the dispersion relation. b) The chimerastates evolution is shown as a function of time where the last module is incoherent compared to the remaining four ones. a ) a ) a )d ) d ) Mod EigenvectorNon-Mod EigenvectorPattern
III d )
100 200 300 400 50050100150200
III
Mod EigenvectorPattern
III b )
100 200 300 400 50050100150200
III c )
100 200 300 400 50050100150200
III
Mod EigenvectorNon-Mod EigenvectorPattern
I II b )c )b )c ) FIG. 3.
Robustness of the symmetry-breaking mechanism for cluster synchronisation and chimera states inmodular networks
We have used the Zhabotinsky model as defined in the main text, with parameters c = c = 27, c = c = 14 . c = c = 1, c = 24 . c = 2 . σ u = σ v = 0, and σ z = 2 .
6. In the columns from left to right are shown thenetworks in increasing ratio between the modules, respectively from a 100 vs. 100 for left panels to a 150 vs. 50, middle panels,and concluding with 190 vs. 10, right panels. a ) − a ) The dispersion relation for the modules of different sizes. In the insetare shown both the real and the imaginary part of the dispersion relation. b ) − b ) Comparison between the eigenvectors ofthe corresponding unstable modes vs. the final pattern. c ) − c ) Time evolution of patterns, from a weakly chimera patternto a cluster synchronized one. d ) − d3