Partial Synchronization and Partial Amplitude Death in Mesoscale Network Motifs
PPartial Synchronization and Partial Amplitude Death in Mesoscale Network Motifs
Winnie Poel, Anna Zakharova, and Eckehard Sch¨oll ∗ Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, Hardenbergstr. 36, 10623 Berlin, Germany
We study the interplay between network topology and complex space-time patterns and introduce a conceptto analytically predict complex patterns in networks of Stuart–Landau oscillators with linear symmetric andinstantaneous coupling based solely on the network topology. These patterns consist of partial amplitude deathand partial synchronization and are found to exist in large variety for all undirected networks of up to 5 nodes.The underlying concept is proved to be robust with respect to frequency mismatch and can also be extended tolarger networks. In addition it directly links the stability of complete in-phase synchronization to only a smallsubset of topological eigenvalues of a network.
I. INTRODUCTION
In recent years the study of networks has increasingly im-proved the understanding of complex systems across scien-tific fields. Interpreting a complex system as a network whosenodes are nonlinear oscillators interacting via the edges of thenetwork graph, one can discover an abundance of collectivephenomena that have also been widely observed in experi-ments. Among others, global synchronization and oscillationquenching mechanism, i.e., amplitude death (AD) and oscil-lation death (OD), have been intensely studied both theoreti-cally and experimentally [1–3]. Amplitude death is associatedwith the stabilization of an already existing trivial steady statewhile oscillation death is characterized by a newly born in-homogeneous steady state. Applications of synchronizationrange from neuronal and genetic networks to coupled electri-cal and laser networks, coupled mechanical systems or evennetworks in social sciences [1, 4–9] while oscillation quench-ing has been observed across many man-made and naturalsystems ranging from lasers and electronic circuits [10, 11],chemical and biological networks including neurons [12–14]to climate systems [15]. Applications of AD are mainlyin controlling physical and chemical systems (e.g., coupledlasers [16]) and suppressing neuronal oscillations [17, 18],while OD has been suggested as a mechanism to generate het-erogeneity in homogeneous systems (e.g., stem cell differen-tiation [19] in morphogenesis). Synchronized behavior is insome cases considered essential for the proper functioning ofthe network (e.g., power distribution networks [20] or coher-ent circadian output in mammals [21, 22]) while in others itis undesired and harmful (e.g., in several pathological neu-ronal states such as Parkinson’s disease and essential tremor[23, 24]). Thus the control of synchronization has also widelybeen studied [25–28].The extensions of the studies of these global phenomena to-wards cluster and group synchronization [29–35], partial am-plitude and oscillation death [3, 36, 37] and chimera states[38–41] are first steps towards a better understanding of theemergence of local mesoscale structures on networks and theirinfluence on the global dynamics [42].The identical temporal behavior of parts of a network iscommon to all of these phenomena and might be a universal ∗ [email protected] concept in coupled systems that one can try to understand in-dependently of the dynamical phenomenon itself on the basisof the network topology. Efforts have recently been made tounderstand the spectral properties and symmetries of complexnetworks [43–45], and from some of the results dynamical im-plications have been inferred [32, 42]. Ground breaking workconnecting the topological properties of the network with dy-namical features was the introduction of the master stabilityfunction for the completely synchronous state [46] and its ex-tension to time-delayed coupling [25, 47–49] and group andcluster synchronization [29, 30, 32].The interplay of local dynamics and the network topologyand their role in the formation of dynamic patterns is nev-ertheless only poorly understood in general, and a systematicand large-scale comparison of the different dynamical patternswhich a certain network topology can display has not beenmade. Two of the issues complicating this task are the com-plex local dynamics each oscillator already exhibits on its ownand the complexity and sheer size of real world networks.This paper addresses this question of a possible connectionbetween the network topology and the dynamical patterns itcan support. The model system used for this purpose shouldbe as simple as possible so that general basic mechanisms maybecome apparent. Therefore a symmetric instantaneous cou-pling and identical elements with simple generic local dynam-ics will be considered. The Stuart-Landau oscillator is such ageneric model of a nonlinear oscillator close to a Hopf bi-furcation. It has been widely studied, and shows the abovephenomena, i.e., synchronization, amplitude and oscillationdeath, and chimera states [50–55]. It is used for the local dy-namics of the model in this paper.The size of the network is the other factor making it difficultto understand general mechanisms clearly. The study of net-work motifs has suggested that these small networks may havedynamical implications or functions on their own [42, 56–58].Therefore, this paper uses the whole variety of undirected mo-tifs with up to five nodes as a reference system for the occur-rence of dynamical patterns. Having such an abundance ofexamples helps to distinguish properties that belong to onespecific topology only from general mechanisms.The paper is structured as follows. Sect. II introduces themodel we use. Sect. III presents the analytical approach de-veloped in our work that can be used to deduce dynamicalpatterns from a network’s topology. The stability of these pat-terns is analyzed in Sect. IV before their robustness with re- a r X i v : . [ n li n . AO ] N ov spect to mismatch in the oscillators’ frequencies is studied inSect. V. Finally, Sect. VI points out how to extend our con-cept to larger networks. II. MODEL
We use the Stuart-Landau oscillator, a generic model of asystem close to a Hopf bifurcation. The dynamics of the i -thoscillator, z i ∈ C , i ∈ { , . . . , N } is given by ˙ z i = f ( z i ) + σe iβ N (cid:88) j =1 A ij z j (1) f ( z ) = ( λ + iω − | z | ) z (2)where ω ∈ R is the oscillator frequency, and λ ∈ R is thebifurcation parameter. For λ > there exists a limit cycle ofradius √ λ that is born in a supercritical Hopf bifurcation at λ = 0 for a single Stuart–Landau oscillator. The real param-eters σ and β are the coupling strength and coupling phase,respectively.The oscillators are connected via a linear bidirectional in-stantaneous coupling using an adjacency matrix, A , with unityrow sum, N (cid:88) j =1 A ij = 1 , (3)to ensure that complete in-phase synchronization can bereached. This matrix A has only real eigenvalues − ≤ η ≤ (proof see appendix A). An example of this normalized adja-cency matrix can be found in Fig. 1b) for a five-node motif.In this paper we focus on the motif depicted in Fig. 1a) as anillustrative example, but similar results have been obtained forall undirected network motifs with up to N = 5 nodes in ourwork [59]. a) b) A =
13 13
13 13
14 14
14 14
13 13
13 13 FIG. 1. a) Motif of 5 coupled oscillators given by the normalized ad-jacency matrix, A, shown in b) . It is used in this paper for illustratinggeneral results that can be applied to all undirected motifs. III. EIGENSOLUTIONS – AN ANALYTIC CONCEPT
This section introduces an analytic solution for the tempo-ral behavior of the oscillators z i ( t ) that are described by thecoupled ordinary differential equations (ODEs) in Eq. (1). Itis obtained by using the following ansatz. An eigenvector v of the normalized adjacency matrix A with eigenvalue ηA v = η v . (4) and eigenvector components v i ∈ {− , , +1 } (5)is used to make the ansatz z i ( t ) = v i z η ( t ) . (6)The function z η ( t ) is to be determined. The following proper-ties of the Stuart-Landau dynamics f (0) = 0 f ( − z i ) = − f ( z i ) (7)are combined in f ( z η v i ) = v i f ( z η ) (8)and used to decouple the system of ODEs in Eq. (1) yielding ˙ z η v i = (cid:0) f ( z η ) + σe iβ z η η (cid:1) v i . (9)If v i = 0 then Eq. (9) is trivially true independently of z η . For v i (cid:54) = 0 dividing it by v i yields ˙ z η = f ( z η ) + ησe iβ z η . (10)We introduce new parameters ˜ λ = λ + ησ cos β ˜ ω = ω + ησ sin β (11)and using that η ∈ R one can write Eq. (10) as ˙ z η = (˜ λ + i ˜ ω − | z η | ) z η . (12)This has the form of the local dynamics of a single Stuart-Landau oscillator given in Eq. (2) and is solved by z η ≡ z ∗ η =0 or the limit cycle z η ≡ z LCη = (cid:112) ˜ λe i ˜ ωt . (13)In summary the system has, in addition to the trivial steadystate, z ∗ η = 0 , a new periodic attractor given by z i ( t ) = v i z LCη ( t ) . (14)if the motif’s topology is such that the adjacency matrix hasan eigenvector consisting only of components v i ∈ {− , , } with eigenvalue η . All considered topologies (undirected mo-tifs up to N = 5 ) have at least one eigenvector of the specialtype needed for the eigensolution. Most motifs exhibit morethan one such eigenvector, and therefore more than one eigen-solution can be found for them.In what follows, we explain how the interplay betweentopology and complex space-time patterns becomes accessi-ble by the concept of an eigensolution.According to the three different values of v i , the oscilla-tors z i can be partitioned into three groups with different tem-poral behavior as sketched in Fig. 2. All oscillators of onegroup behave identically, i.e., they are fully synchronized inamplitude and phase. The oscillators for which v i = 0 holds(green/light gray group), are in a steady state z i = 0 . Oscil-lators with v i = ± (red/dark gray and blue/black group) areon the limit cycle z η given by Eq. (13) with a fixed phase dif-ference of φ − φ = π between the two groups. Thus eigen-solutions are spatial patterns of amplitude death (green/lightgray group) and partial in-phase and anti-phase synchroniza-tion (blue/black and red/dark gray). We call two oscillators, z i = r i e iφ i and z j = r j e iφ j with r i = r j , in-phase syn-chronized if | φ i − φ j | = 0 and anti-phase synchronized if | φ i − φ j | = π . r = (cid:112) ˜ λ, φ = φ r = 0 r = (cid:112) ˜ λ, φ = φ = φ + π FIG. 2. (color online) Sketch of the different oscillator groups oc-curring in an eigensolution. The red (dark grey) and blue (black)group are on the limit cycle given by Eq. (13) with phase difference | φ − φ | = π . Oscillators of the green (light gray) group are in thesteady state at z i = 0 . The spatial distribution of amplitude death and partial syn-chronization is governed solely by the eigenvector defining theeigensolution and is therefore purely determined by the topol-ogy. The motif shown in Fig. 1 has the following eigenvectors v k and corresponding topological eigenvalues η k v = (1 , , , , T η = 1 v = (1 , − , , , − T η = − / v = (0 , , , , − T η = 0 v = (1 , , , − , T η = 0 v = (1 , , − , , T η = − / . (15)Of these eigenvectors all but v can be used to obtain aneigensolution. They yield the patterns shown in Fig. 3 a) to d).The pattern shown in c) is a linear combination of the eigen-vector v and v , which are equivalent if one renumbers theoscillators, and observes that they have the same eigenvalue η = 0 . In addition, the system can be in the trivial steadystate, shown in part e) of Fig. 3.A time series of the pattern PAD+PS 1 (Fig. 3c)) for eachoscillator is shown in Fig. 4. The separation of the oscilla-tors into three groups is clearly visible. Oscillators of the blue(black) and the red (dark gray) group oscillate in anti-phase onthe limit cycle while the green (light gray) one is in the steadystate at z i = 0 . As can be seen in Eqs. (11) and (13), the limitcycle of the coupled system has a direct dependence on thetopological eigenvalue η via its radius r LC = (cid:112) ˜ λ and phase φ LC = ˜ ωt . In comparison to the single Stuart-Landau oscilla-tor’s limit cycle, z SL = √ λe iωt , the coupled system’s radiusis altered by the additional term ησ cos β and its frequency ischanged by ησ sin β . a) b) c) d) e)CIS PAD+PS 1 PAD+PS 2 PAD+PS 3 CADFIG. 3. (color online) All eigensolutions for one motif as definedby the eigenvectors in Eq. (15): a) v , b) v , c) v + v , d) v .The eigensolutions are patterns of partial amplitude death (green orgray) and partial synchronization (red or light gray, and blue or darkgray). For more details on the color code refer to legend in Fig. 2.Additionally, the steady state at z i = 0 ∀ i ∈ { , . . . , N } is shownin e). These patterns will further be referred to by the acronyms CIS(complete in-phase synchronization), PAD (partial amplitude death),PS (partial synchronization), CAD (complete amplitude death). Note that the example in Fig. 4 is for negative λ where thesingle oscillator is in a stable steady state. The limit cycleof the coupled system is shifted by the additional term in ˜ λ ,and therefore the coupled system can exhibit oscillations eventhough the uncoupled system is in a stable steady state. Ingeneral the existence of an eigensolution is determined by ˜ λ ≥ and thus λ + ησ cos β ≥ . (16)This inequality defines a region in the parameter space of λ , σ and β in which the eigensolution corresponding to η exists.Since one motif may have several eigensolutions with differ-ent topological eigenvalues, it crucially depends on the pointin parameter space which spatio-temporal patterns a givenmotif can exhibit.
990 992 994 996 998 1000time1.00.50.00.51.0 analyticalOscillator 1Oscillator 2Oscillator 3Oscillator 4Oscillator 5 R e ( z i ) FIG. 4. (color online) Time series for the PAD+PS 1 pattern shownin Fig. 3b). The black line gives the analytical eigensolution whilethe colored dots are the result of numerically integrating the coupledODEs from Eq. (1) from random initial conditions up to t = 1000 time steps. Oscillator 3 is in the steady state at z = 0 while theothers are on the limit cycle z LCη given by Eq. (13). Oscillator 1and 4 have a fixed phase difference of ∆ φ = π to oscillator 2 and 5.Parameters: β = 0 , σ = − , λ = − , ω = 2 . In conclusion, this section has shown that eigensolutionsare an analytical method to infer complex space-time patternsfrom topology. They are built from three groups of oscilla-tors, one being amplitude death and the other two oscillatingin anti-phase. The spatial ordering of these groups depends ona topological eigenvector while temporal dynamics and exis-tence in parameter space depend on the corresponding topo-logical eigenvalue.
IV. LINEAR STABILITY ANALYSIS
In this section, we study the stability of the patterns con-structed in Sect. III. We start with patterns without partialamplitude death (PAD) which can be treated in polar coordi-nates in contrast to the case including PAD where the deadoscillator’s phase is undefined. These will be treated sepa-rately in Sect. IV B. Finally Sect. IV C analyzes the stabilityof the trivial steady state, and our results are illustrated for anexample in Sect. IV D.
A. Patterns Without Partial Amplitude Death
Introducing polar coordinates z j = r j e iφ j in the system ofcoupled oscillators given by Eqs. (1) and (2) yields ˙ r i = ( λ − r i ) r i + σ (cid:88) j A ij r j cos( β + φ j − φ i ) ,r i ˙ φ i = ωr i + σ (cid:88) j A ij r j sin( β + φ j − φ i ) . (17)We introduce small perturbations, δr i and δφ i , of the limitcycle solution z i = v i z LCη in radius and phase r i = r LCi + δr i = r η + δr i ,φ i = φ LCi + δφ i = φ η + π v i + δφ i . (18)Inserting Eq. (18) into Eq. (17) and expanding all right handsides around δr i = 0 , δφ i = 0 up to first order yields ˙ δr i =( λ − r η ) δr i + σ N (cid:88) j =1 A ij (cid:104) cos (cid:16) β + π v j − v i ) (cid:17) δr j − r η sin (cid:16) β + π v j − v i ) (cid:17) ( δφ j − δφ i ) (cid:105) , ˙ δφ i r η = − ησ sin βδr i + σ N (cid:88) j =1 A ij (cid:104) sin (cid:16) β + π v j − v i ) (cid:17) δr j + r η cos (cid:16) β + π v j − v i ) (cid:17) ( δφ j − δφ i ) (cid:105) . (19)Using v i , v j ∈ {− , +1 } and therefore π ( v j − v i ) ∈{− π, , π } one can write cos (cid:16) β + π v j − v i ) (cid:17) = v i v j cos β , sin (cid:16) β + π v j − v i ) (cid:17) = v i v j sin β . Making use of (cid:80) j A ij v j = ηv i , as well as v i = 1 and ad-ditionally introducing ξ i = ( δr i , δφ i ) T , the Eqs. (19) can be rewritten as ˙ ξ i = (cid:18) λ − r η r η ση sin β − r η ση sin β − ση cos β (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) J ξ i + σ N (cid:88) j =1 v i A ij v j (cid:18) cos β − r η sin β r η sin β cos β (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) R ξ j . (20)After introduction of the diagonal matrix V = diag ( v ) (withcomponents of the eigenvector v in its main diagonal) and ξ = ( ξ , ξ , . . . , ξ N ) T , the system of equations takes the form ˙ ξ = [ I N ⊗ J + σ ( V AV ) ⊗ R ] (cid:124) (cid:123)(cid:122) (cid:125) C ξ . (21)The dynamical eigenvalues of the N × N matrix C (de-fined in Eq. (21)) determine the stability of the eigensolution.Fortunately, one can show that the matrix V AV in the secondterm of C is diagonalizable and has the same eigenvalues as A , see appendix B.We can therefore transform Eq. (20) in a way that diago-nalizes V AV . We denote the transformed perturbations as χ i with corresponding dynamics ˙ χ i = ( J + σµ i R ) χ i . (22)The i th topological eigenvalue of A is now denoted as µ i (in-stead of η i ) in order not to confuse it with the particular eigen-value of A that was used to obtain the eigensolution and whichwe will keep referring to as η . So there is (at least) one index i ∈ { , . . . , N } for which µ i = η is true. Without loss ofgenerality one may assume the topological eigenvalues to besorted such that µ ≤ µ ≤ · · · ≤ µ N . (23)Now in general the real parts of the eigenvalues of the × matrices M i = J + σµ i R, i ∈ { , · · · , N } (24)need to be negative to secure stability of the eigensolution.The characteristic equation for the dynamical eigenvalues α i of this 2 × α i − Tr [ M i ] α i + Det [ M i ] = 0 . (25)According to the Hurwitz criterion this equation has solutionswith negative real part if and only ifTr [ M i ] = − (cid:2) λ + σ cos( β )(2 η − µ i ) (cid:3) < , (26)Det [ M i ] = ( η − µ i ) σ (cid:2) λ cos( β ) + σ { η cos(2 β ) . +2 η − µ i } (cid:3) > . (27)We demand that the above inequalities Eq. (26) and Eq. (27)are fulfilled for alla) µ i (cid:54) = η, b) σ (cid:54) = 0 , c) λ > − ησ cos β . (28) σ cos β < σ cos β > µ i > η λ > h ( µ i , η, β ) σ λ < h ( µ i , η, β ) σµ i < η λ < h ( µ i , η, β ) σ λ > h ( µ i , η, β ) σ TABLE I. Inequality Eq. (27) expressed in terms of the function h defined in Eq. (30) (white cells). Four different cases of parametervalues σ and β have to be distinguished as indicated by the top rowand left column (gray cells). These conditions restrict the stability analysis to a) perturba-tions perpendicular to the synchronization manifold, b) thecase of coupled oscillators and c) the region of parameterspace in which the eigensolution exists. If η is degenerate, µ i = η has to be considered since it no longer solely belongsto a perturbation inside the synchronization manifold. In thiscase there will always be a dynamical eigenvalue which is zeroand the linear stability analysis remains inconclusive.There are thus N − inequalities given by Eqs. (26) and(27) determining the stability of the eigensolution belongingto the topological eigenvalue η . They define regimes of sta-bility in the parameter space of λ, σ and β . These inequalitieswill now be used to analyze the boundaries of stability in the( σ , λ ) plane (thus for constant β ) using the definitions g ( µ i , η, β ) = ( µ i − η ) cos β, (29) h ( µ i , η, β ) = 12 cos β (cid:2) µ i − η (cos 2 β + 2) (cid:3) . (30)The condition from Eq. (26) then reads λ > g ( µ i , η, β ) σ (31)whereas inequality Eq. (27) yields the four cases given inTable I.For constant β the two columns in Table I distinguish be-tween σ < and σ > . Hence in the ( σ , λ ) plane thetwo cases that are in the same column define a set of lower( λ > h ( µ i ) ) and a set of upper ( λ < h ( µ i ) ) boundaries for λ that are given by the different µ i > η (top row) or µ i < η (bottom row).We thus obtain a series of lower and upper bounds of λ forpositive σ and a possibly different series of upper and lowerbounds of λ for negative σ . These need to be combined withyet another set of lower bounds for λ defined by the inequali-ties from Eq. (31) for all the µ i (cid:54) = η and the criterion for theexistence of the solution given in Eq. (28).These boundaries are a collection of lines in the ( σ, λ ) planethat are each going through the origin with a slope determinedby either g ( µ i , η, β ) , h ( µ i , η, β ) , or − η cos β . One can findtwo of these lines that yield the maximum constraint for λ inthe regime of σ > or σ < such that in each of theseregimes just one upper and one lower boundary for λ can begiven, each only depending on one selected topological eigen-value µ i and fixed β and η . The inequalities limiting the re-gion of stability can be found in the white cells of Table IIfor different cases of parameters β and σ which are given in the gray cells. Dark gray cells apply to all white cells in theircolumn or row while parameter conditions in light gray cellsonly apply to the white cell right below them. Crosses in thewhite cells indicate that for this parameter combination theeigensolution is always unstable.As can be seen by the different columns of Table II, thestability region of the eigensolution belonging to the topolog-ical eigenvalue η differs if η is the smallest ( η = µ ), thelargest ( η = µ N ), or some intermediate eigenvalue ( η = µ j , j / ∈ { , N } ) . Note that there are critical values of the couplingphase β for each of these cases defined by the inequalities inthe light gray cells.To illustrate Table II we consider the motif from Fig. 1. Theonly pattern without partial amplitude death is the completelyin-phase synchronized state (CIS), Fig. 3a). It belongs to η = µ N = 1 , the largest topological eigenvalue. According to thesecond column of Table II critical values for β are determinedby cos(2 β ) = µ − µ N − µ N − µ and cos(2 β ) = 0 . (32)The eigenvalues µ i of this motif are µ = − , µ = − , µ = µ = 0 , µ = 1 . (33)The critical values of β ∈ { , π } are thus ( β ≈ . ∨ β ≈ . ⇔ cos(2 β ) = − (34) ( β = π ∨ β = 3 π ⇔ cos(2 β ) = 0 . (35)Increasing the coupling phase from β = 0 when cross-ing β = 0 . , a region of stability between λ = g ( µ ) σ and λ = h ( µ N − ) σ is created and later destroyed again at β = 2 . as given by the upper part of the second column inTable II where η = µ N and σ cos β < .When the second set of critical values, β = π/ and β = 3 π/ , is crossed, the lower boundary of the region of sta-bility changes between λ = g ( µ N − ) σ and λ = h ( µ N − ) σ as indicated by the lower part of the second column in TableII where η = µ N and σ cos β < . Note that for fixed β therows of Table II differentiate between positive and negative σ .Depending on the value of β we thus have either just one sta-ble region limited by the axis at σ = 0 and a lower boundary( β < . or β > . ) or we get an additional region onthe other side of the σ = 0 axis that has an upper and a lowerboundary ( . < β < . ).These regions in which CIS is stable can be seen for dif-ferent values of β in Fig. 5. One clearly sees the appearanceof an additional regime of stability when the critical value of β = 0 . is crossed as discussed above. B. Patterns Including Partial Amplitude Death
Here we consider patterns including partial amplitudedeath. We need to keep cartesian coordinates z j = x j + i y j to TABLE II. Conditions for the regimes of stability of an eigensolution belonging to η in the parameter space of λ , σ and β . This table onlyapplies to eigensolutions that do not include partial amplitude death. The inequalities limiting the stable region are given in the white cells.Which of these white cells needs to be considered for a fixed set of β , σ and η is determined by the conditions in the gray cells. Dark graycells apply to all white cells in their row and column while light gray cells only apply to the white cell directly below them. η = µ η = µ N η = µ j cos(2 β ) < β ) < µ − µ N − µ N − µ cos(2 β ) < µ − µ j +1 µ j − µ λ > g ( µ ) σ g ( µ ) σ < λ < h ( µ N − ) σ h ( µ j +1 ) σ < λ < h ( µ j − ) σσ cos β < β ) > β ) > µ − µ N − µ N − µ µ − µ j +1 µ j − µ < cos(2 β ) < µ − µ j − µ j − µ g ( µ ) σ < λ < h ( µ j − ) σλ > h ( µ ) σ cos(2 β ) > µ − µ j − µ j − µ cos(2 β ) < µ N − µ µ − µ N cos(2 β ) < β ) < µ N − µ j − µ j − µ N g ( µ N ) σ < λ < h ( µ ) σ λ > g ( µ N − ) σ h ( µ j − ) σ < λ < h ( µ j +1 ) σσ cos β > β ) > µ N − µ µ − µ N cos(2 β ) > µ N − µ j − µ j − µ N < cos(2 β ) < µ N − µ j +1 µ j − µ N g ( µ N ) σ < λ < h ( µ j +1 ) σλ > h ( µ N − ) σ cos(2 β ) > µ N − µ j +1 µ j − µ N σ a) β =0 . σ b) β =1 . σ c) β =1 . σ d) β =1 . FIG. 5. (color online) Stability diagrams for the state of complete in-phase synchronization (CIS, Fig. 3a)) for the motif shown in Fig. 1.The plots depict the region in the ( σ , λ ) plane in which CIS is stableas colored in red (gray) for a) β = 0 , b) β = 1 . , c) β = 1 . ,d) β = 1 . . The regions are given by the inequalities in the whitecells in Table II. They are whitened wherever the eigensolution doesnot exist according to Eq. (16). One clearly sees the appearance of asecond stable regime when the critical value of β = 0 . is crossed. analyze the stability of solutions containing dead oscillatorsvia a linear stability analysis. Using these coordinates equa-tions Eq. (1) read ˙ x j = ( λ − x j − y j ) x j − ωy j + σ N (cid:88) k =1 A jk (cos βx k − sin βy k ) , ˙ y j = ( λ − x j − y j ) y j + ωx j + σ N (cid:88) k =1 A jk (sin βx k + cos βy k ) . (36)We write the eigensolution as z j ( t ) = v j z η ( t ) = v j x η ( t ) + iv j y η ( t ) (37) where x η = r η cos φ η y η = r η sin φ η with r η = (cid:112) ˜ λφ η = ˜ ωt (38)and introduce small pertubations δx j and δy j x j = x η v j + δx j y j = y η v j + δy j . (39)Linearizing around δx j = 0 , δy j = 0 leads to ˙ δx j ˙ δy j = (cid:20) ( λ − v j ˜ λ ) + ω −
11 0 (cid:124) (cid:123)(cid:122) (cid:125) J (cid:21) δx j δy j − v j ˜ λ cos (2˜ ωt ) sin (2˜ ωt )sin (2˜ ωt ) − cos (2˜ ωt ) (cid:124) (cid:123)(cid:122) (cid:125) J ( t ) δx j δy j + σ (cid:88) k A jk cos β − sin β sin β cos β (cid:124) (cid:123)(cid:122) (cid:125) R δx k δy k . (40)In contrast to the calculations in polar coordinates the Jaco-bian here depends on time as it contains the matrix J ( t ) . Toeliminate this dependence we need to change to a co-rotatingframe using the rotation matrix S ( t ) = cos(˜ ωt ) sin(˜ ωt ) − sin(˜ ωt ) cos(˜ ωt ) . (41)In the transformed coordinates δ ¯ x i δ ¯ y i = S δx i δy i (42)and after introducing ξ j = ( δ ¯ x j , δ ¯ y j ) T and ξ =( ξ , ξ , . . . , ξ n ) T Eq. (40) reads ˙ ξ = (cid:104) I N ⊗ J A − ˜ λV ⊗ J B + σA ⊗ R (cid:105)(cid:124) (cid:123)(cid:122) (cid:125) C ξ , (43)with V = diag ( v i ) and J A = λ ησ sin β − ησ sin β λ , J B = . (44)Unfortunately, in general V and A do not commute so thatwe can not reduce the eigenvalue problem as we did in the pre-vious subsection. The eigenvalues of C whose real part, theLyapunov exponent, determines the stability of the eigensolu-tion therefore need to be calculated numerically for the set ofparameters one is interested in. The range of λ ∈ ( − , , σ ∈ ( − , , β ∈ (0 , π ) , ω = 2 was systematically scannedfor all motifs. The only significantly large stability region ofa PAD+PS state found in this scan is the pattern shown in Fig.3b) named PAD+PS1. This stable region is colored in green(light gray) in Fig. 6 for two different values of β . C. The Steady State
The stability of the trivial steady state, z i = 0 ∀ i ∈{ , . . . , N } , can be understood completely analytically. If weconsider it as a special case of a partially dead state, namelythe one with v i = 0 ∀ i , we can use some results of the previ-ous section.We start from Eq. (40) where we insert v j = 0 to obtain ˙ δx j ˙ δy j = λ − ωω λ (cid:124) (cid:123)(cid:122) (cid:125) J δx j δy j + σ N (cid:88) k =1 A jk R δx k δy k (45)or with ξ j = ( δx j , δy j ) T and ξ j = ( ξ , ξ , . . . , ξ N ) T ˙ ξ = [ I N ⊗ J + σA ⊗ R ] (cid:124) (cid:123)(cid:122) (cid:125) C ξ . (46) We switch to new coordinates that diagonalize A , and the tran-formed perturbations χ i then obey ˙ χ i = λ + µ i σ cos β − ( ω + µ i σ sin β ) ω + µ i σ sin β λ + µ i σ cos β (cid:124) (cid:123)(cid:122) (cid:125) C i χ i (47)where µ i are the eigenvalues of the adjacency matrix A . Wecan apply the criteria from Eq. (26) and Eq. (27)Tr [ C i ] = 2( λ + µ i σ cos β ) < , Det [ C i ] = ( λ + µ i σ cos β ) + ( ω + µ i σ sin β ) > . The first condition can be rewritten as λ < − µ i σ cos β . (48)If this condition is met, then ( λ + µ i σ cos β ) (cid:54) = 0 , and thusthe second condition is automatically fulfilled as well. Whichof the eigenvalues µ ≤ µ ≤ · · · ≤ µ n imposes the strongestrestriction on λ coming from Eq. (48) is determined by thesign of σ cos β . It is either the largest, µ N = 1 , or the smallest, µ . TABLE III. Analytic stability conditions of the trivial steady state fortwo different conditions on the coupling parameters σ and β .parameter condition σ cos β < σ cos β > stability condition λ < − µ cos βσ λ < − cos βσ Thus the stable regime of the trivial steady state in the ( σ , λ )plane is bounded by two lines through the origin, one withslope − cos β and one with − µ cos β as given in Table III.Now remember that the eigensolution to µ i exists if λ ≥− µ i σ cos β (compare Eq. (16)). Thus the above conditionsare fulfilled and the trivial steady state is stable if and only ifnone of the eigensolutions exists. This stable region is coloredin blue (black) in Fig. 6. D. Example
The results of the above stability analysis are combined forthe motif of Fig. 1 in the stability diagrams in Fig. 6 for β = 0 and β = 1 . . The two values of β were chosen to representthe case without and with the additional region of stability forCIS (see Fig. 5a) and c)). The stability of CAD and CIS areanalytical results given in Table II and III respectively, whilethe stable region of PAD+PS 1 was obtained by numericallycalculating the dynamical eigenvalues from Eq. (43) for thedifferent sets of parameters λ , σ and β for η = − / . V. ROBUSTNESS UNDER FREQUENCY MISMATCH
Real world systems are typically subject to heterogeneitiesin the parameters while in our model the oscillators are homo- a)
10 5 0 5 10 σ β =0 . λ b)
10 5 0 5 10 σ β =1 . λ CAD CIS PAD + PS 1 other
FIG. 6. (color online) Semi-analytic stability diagram for the homo-geneous motif of Stuart-Landau oscillators as shown in Fig. 1 for a) β = 0 and b) β = 1 . . The regions are colored according to the pat-tern from Fig. 3 that is stable there. The dark grey region denoted as‘other’ indicates that none of the known patterns is stable. The stabil-ity of CAD and CIS is known purely analytically from tables II andIII, while the stability of the partially dead solutions PAD+PS 1 to 3was calculated numerically from Eq. (43). Of these only PAD+PS 1is found to be stable. geneous. To address this problem we numerically test our pre-vious observations for heterogeneous oscillator frequencies.The system of coupled differential equations obtained fromEq. (1) by introducing individual oscillator frequencies ω i is ˙ z i = f i ( z i ) + σe iβ N (cid:88) j =1 A ij z j f i ( z ) = ( λ + iω i − | z | ) z . (49)The previously introduced analytical method fails for het-erogeneous frequencies ω i . For the numerical investigationof coupled Stuart-Landau oscillators the ordinary differentialequations given in Eq. (49) that describe the system’s dynam-ics are integrated numerically until transients have vanishedand the system assumes an asymptotic state. The last few pe-riods of the obtained time series are then used to classify theasymptotic state reached for the chosen parameter values of λ , σ , β .Ideally, one would systematically scan the entire parame-ter space, but here we will focus on the cases of β = 0 and β = 1 . in the coupling phase to keep the numerical expensereasonable for this study. These cases refer to a real coupling( β = 0 ) and an example of complex coupling ( β = 1 . )for which the state of complete in-phase synchronization isknown to have an additional region of stability as can be seenin Fig. 6 b) for the homogeneous system. For σ and λ theparameter space is scanned systematically within the range ( − , . This range may even be sufficient since the stabil-ity analysis in Sect. IV shows that knowledge of the systemclose to the origin explains stability for all other values of λ and σ .Because there might be multistability in our system we haveto choose initial conditions carefully. We are looking for theeigensolutions sketched in Fig. 3 and thus use states closeto them as initial conditions. To do so, one snapshot in time of the analytic solution for z i ( t ) is taken and a random sig-nal with entries of the order of 5% of the solution’s radiusis added. Additionally random initial conditions are tested,where the real and imaginary part of each oscillator are cho-sen randomly from the interval (-1, 1).Once the numerical integration is finished, the very last partof the time series is used to classify the behavior. The differentcategories are: • Complete Amplitude Death (CAD): All oscillatorshave a radius equal to zero. • Complete In-phase Synchronization (CIS): All oscil-lators have the same constant finite radius and are in-phase synchronized. • Partial Amplitude Death and Partial Synchroniza-tion (PAD + PS): Some of the oscillators have a ra-dius equal to zero and the others all have the same con-stant finite radius and are anti-phase or in-phase syn-chronized. The spatial distribution of these groups isaccording to one of the patterns in Fig. 3. There is thusone category for each of the patterns that include partialamplitude death. They are labeled as shown in Fig. 3. • Partial Amplitude Death (PAD): Some of the oscilla-tors have a radius equal to zero. The others are not syn-chronized as in one of the previously defined patterns.The algorithm checks the last 50 out of 8000 time steps of thesimulation for these properties as follows. A critical thresh-old r crit is defined as 1% of the maximum radius occurring inall the oscillators in this time interval. If this is smaller than10 times an internal numerical threshold, the machine epsilon[60], the state is considered to be CAD. If not, each oscil-lator is classified by its radius. A radius smaller than r crit is considered zero and a signal with r min ≥ . r max and r max > r crit is identified as a constant finite radius. If themean radii of all oscillators that are not dead do not differamong each other for more than r crit , these oscillators are re-garded to have the same radius.As a measure for synchronization among the active ( r (cid:54) = 0 )oscillators a slightly modified version of the Kuramoto or-der parameter [61] is used. The Kuramoto order parameteris equal to unity for complete synchronization and zero forcomplete desynchronization. It is commonly defined as R = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) j e iφ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (50)where N is the number of oscillators whose synchronizationbehavior one is interested in. We introduce a slight changethat takes into account the fixed phase differences of the os-cillators in the investigated patterns and can define a similarorder parameter for the synchronization pattern belonging tothe eigensolution given by v k R k = 1 N a (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) j =1 v kj e iφ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (51)where N a is the number of active oscillators in the pattern, N a = (cid:80) j | v j | , and v kj is the j -th component of the k -th eigen-vector of A . The parameter R k becomes unity if the oscilla-tors show the pattern of in-phase and anti-phase synchroniza-tion defined by v k . We consider the Stuart-Landau oscillatorssynchronized if R ≥ . over the last 50 time steps.The plots obtained for random initial conditions are dis-played in Fig. 7 for β = 0 and β = 1 . together with theplot with initial conditions close to the CIS state for β = 1 . ,which shows the additional stable CIS region already seen inFig. 6 for homogeneous oscillators. All other initial condi-tions yield results similar to panels a) and b) in Fig. 7. Thus,except for the additional CIS region for β = 1 . , the eigenso-lutions do not seem to have any coexisting stable attractors.The numerical stability diagrams for the heterogeneous sys-tem exhibit a structure that is very similar to the analytic dia-grams for the homogeneous case depicted in Fig. 6. All cate-gories are found in both diagrams and in roughly the same re-gion of parameter space. The first main differences lies in thePAD+PS 1 region. It is slightly smaller for the heterogeneousthan for the homogeneous system and is surrounded by a re-gion in which PAD without any special type of synchronizedbehavior is present. This suggests that even though PAD maybe robust to frequency mismatch, the partial synchronizationmay not. On the other hand it is also possible that in the het-erogeneous system very long transients are encountered andthat thus the PAD+PS 1 state may still be reached after longerintegration times in the PAD region.The other main difference, the change in the CIS region, isonly seen for β = 1 . . For β = 0 the border of the CIS andCAD state shows some states classified as ‘other’ but these aremost likely due to very long transients at the border of stabil-ity where the Lyapunov exponent is close to zero. For β = 1 . though the shape of the CIS region changes significantly andin a similar way for all initial conditions and randomly drawnsets of frequencies. Thus this effect might be a general resultof the heterogeneity and not be due to merely numerical is-sues as the measure of synchronization or long transients. Thenew shape suggests that with heterogeneity a critical couplingstrength (depending on λ ) has to be reached before the stateof complete in-phase synchronization is stable. This seemsplausible since the uncoupled oscillators all have different fre-quencies and would never synchronize unless forced by theirinteraction. VI. TOWARDS COMPLEX NETWORKS
Here we present first steps in extending the concept ofeigensolutions to larger networks.For this purpose two approaches are chosen. First, we willextend the network by connecting motifs or adding a deadnode, and second, we will identify special eigenvectors thatare present independently of system size for the coupling ma-trix A [62]. a)
10 5 0 5 10 σ β =0 . λ b)
10 5 0 5 10 σ β =1 . λ otherPADPAD + PS 1CISCAD c)
10 5 0 5 10 σ β =1 . λ FIG. 7. (color online) Numerical stability diagrams for the heteroge-neous motif of Stuart-Landau oscillators as shown in Fig. 1. Panelsa) and b) have random initial conditions while panel c) was initiallyclose to the CIS state. Each panel has its own set of random fre-quencies ω i = ω · rand (0 . , . , ω = 2 . The color of each pointindicates the state that has been reached after numerically integratinguntil t = 8000 . These panels differ from the analytic stability dia-gram of the homogeneous system, Fig. 6, in the size of the PAD+PS1 region and the shape of the CIS region. One can also see multista-bility in the additional CIS region for β = 1 . . A. Connecting Motifs Via a Dead Node
If two motifs each have an eigensolution with a dead nodethese nodes can be combined to form a larger network as il-lustrated in Fig. 8. The combined network then has a spatio-temporal pattern given by the two motifs’ eigensolutions. Informulas this reads as follows. + =
13 2 45
FIG. 8. (color online) Two motifs with partially dead eigensolutionscan be connected at their dead nodes without changing their patterns.This may create patterns with more than one frequency.
Let the two single motifs have eigensolutions z = z η v η ∈ R M , z = z µ v µ ∈ R N containing each a dead node. Without loss of generality wecan fix them as v ηM = 0 and v µ = 0 . Then the resulting largermotif with M + N − nodes has the solution z ∈ R M + N − z i = z η v ηi if i ∈ { , . . . , M − } if i = Mz µ v µi − ( M − if i ∈ { M + 1 , . . . , M + N − } . (52)This is also true if one or both of the eigenvectors are replacedby the zero-vector of the same dimension. Thus one can attachan arbitrarily shaped network to an eigensolution’s dead node,and still the eigensolution remains valid on the small motifwhile the attached network is in the trivial steady state.If the two eigensolutions belong to different eigenvalues, η (cid:54) = µ , the combined solution will have more than one fre-quency and amplitude because of z η = (cid:112) λ + ησ cos βe i ( ω + ησ sin β ) t z µ = (cid:112) λ + µσ cos βe i ( ω + µσ sin β ) t . (53)For example, in the case shown in Fig. 8, the left (trian-gular) motif’s eigensolution belongs to η = − / while theright (linear) motif’s eigensolution belongs to µ = 0 . Thusin the combined motif the oscillators 1 and 2 coming origi-nally from the triangular motif have a different frequency andamplitude than oscillators 4 and 5 coming originally from thelinear motif z = z = (cid:114) λ − σ cos β e i ( ω − σ sin β ) t ,z = 0 ,z = z √ λ e iωt . This combined pattern is no longer based on a single eigen-vector of the combined networks adjacency matrix, but ratheron a linear combination of eigenvectors.
B. Attaching a Dead Node to a Set of Active Nodes
If an eigensolution corresponds to η = 0 , one can attacha dead oscillator by connecting it such that its neighbors sumup to zero. If we call the newly created network’s adjacencymatrix (containing the additional dead node) B this translatesto (cid:88) j B ij z j = 0 . For other values of η this does not work. Some examples canbe seen in Fig. 9 where a dead node is attached in two differentmanners to a 4-node ring.Together these two mechanism make it possible to createlarge networks containing synchronization clusters with dif-ferent frequencies and partial amplitude death. Given an arbi-trary network one can always place as many motif eigensolu-tions on it as possible under the constraint of connecting themvia dead nodes. All the other nodes can then be set to zeroand thus one has constructed a solution of the large arbitrarynetwork. a) b) c)FIG. 9. (color online) A dead node (green/light gray) can be attachedto a network without changing an eigensolution belonging to η = 0 if all neighbors of the dead node sum up to zero. In this example thedead node is attached to the 4-node motif with partially synchronizedeigensolution shown in a). There are two possible ways to do thiswhich are depicted in b) and c). C. Eigensolutions Existing For All Network Sizes
As pointed out in Sect. II the CIS state exists for all net-work topologies independently of size because of our choiceof the coupling matrix A . Also, one can show that for bi-partite graphs the matrix A always has an eigenvector withinterchanging entries v i = 1 and v j = − such that oscilla-tors with v i = 1 are only connected to those with v j = − .This is another partially synchronized eigensolution that is in-dependent of the network size and it always corresponds tothe topological eigenvalue η = − .Our analytic stability analysis from Sect. IV holds for allsystem sizes. According to Table II the stability of the CISstate ( η = µ N = 1 ) only depends on the second largest, µ N − and the smallest topological eigenvalue, µ . The stability ofthe PS state on a bipartite network belonging to η = µ = − is governed by the second smallest, µ and the largest, µ N ,topological eigenvalue. Since by construction of A we alwaysget µ N = 1 the stability really only depends on µ in thiscase.Thus our work analytically describes two general typesof behavior for an arbitrarily large number of homogeneousStuart-Landau oscillators coupled via a normalized adjacencymatrix. One could possibly think of an application where theCIS state can be (de)stabilized by adding or removing certainnodes or edges as observed, e.g., for power grids [9]. In orderto do this one would not need to compute the entire spectrumof topological eigenvalues but only the second largest and thesmallest. An example of how one may change the topology a) b) c) d)FIG. 10. (color online) The topology of the motif from Fig. 1 ischanged step by step by detaching the middle node. This processchanges the topological eigenvalues (Table IV) and thus the stabilityof the CIS state as can be seen in Fig. 11. to influence the stability of CIS is shown in Fig. 10. The pre-viously discussed five-node motif is altered step by step byslowly detaching the middle node. The change in the small-est and second largest eigenvalue introduced this way can beseen in Table IV. The different regions of stability defined bythese topological eigenvalues according to the upper part of1the second column in Table II are presented in Fig. 11. TABLE IV. Smallest, µ , and second largest, µ N − , topologicaleigenvalue of the motifs shown in Fig. 10. These eigenvalues de-termine the stability of the CIS state that can be seen in Fig. 11.a) b) c) d) µ N − ( − √
13 1 √ µ − − (1 + √ − - 1
10 8 6 4 2 0 σ λ β =1 . β =1 . β =1 . β =1 . a)b)c)d) FIG. 11. (color online) Analytic stability diagram for the motifsshown in Fig. 10 a) to d). The CIS state is stable in the regionshaded according to the legend on the right. There is just one con-nected region for each color with dark colors overlaying the lighterones. Changing the network topology clearly changes the stabilityof the CIS state. The change is governed solely by the change of thesmallest and second largest topological eigenvalue given in Table IV.
VII. CONCLUSION
The main achievement of this paper is the development ofeigensolutions as a concept to analytically predict the coexis-tence of partial in-phase and anti-phase synchronization withpartial amplitude death. The connection between topologyand dynamics becomes clearly visible and analytically acces-sible in this concept. One could use eigensolutions to derivedynamical patterns of partial amplitude death and partial syn-chronization solely from the network topology. This, how-ever, requires a thorough understanding of the special eigen-vectors needed for its construction, especially of the condi-tions under which they are present in a network. Thereforeinvestigating the occurrence of these eigenvectors is the nextimportant step that should be taken in future work.The stability of eigensolutions without partial amplitudedeath was determined to be dependent directly on a small sub-set of the network’s topological eigenvalues, Sect. IV. Af-ter the introduction of eigensolutions this is the second strongconclusion about the direct interplay of dynamics and topol-ogy made in this work. For the general case of eigensolutionscontaining partial amplitude death and partial synchronizationthe stability for a set of parameters can be calculated numer-ically from the analytically obtained time independent Jaco- bian. This Jacobian should be further analyzed and may yieldanalytic stability conditions in the form of inequalities as well.The prediction of dynamical patterns via eigensolutionshas proved not to be restricted to networks of homogeneousStuart-Landau systems. The same patterns of partial ampli-tude death and partial synchronization have been observedfor heterogeneous frequencies. This broad appearance of thedynamical patterns predicted by the eigensolutions raises thequestion if there is a general concept present that one may beable to use in more complex systems.Finally, the stability and occurrence has been fully linkedto the network topology for complete in-phase synchroniza-tion as well as for patterns of interchanging in-phase and anti-phase synchronized oscillators on bipartite graphs. Their sta-bility is shown to depend solely on two of the network’s topo-logical eigenvalues while their existence is independent of thenetwork’s size and shape. This implies that the topology canbe used to (de)stabilize these states, an aspect linked to thecontrol on complex systems that should be further investi-gated.
VIII. ACKNOWLEDGMENTS
This work was supported by DFG in the framework of SFB910. Helpful discussions with Wolfram Just, Kenneth Showal-ter and Mark Tinsley are gratefully acknowledged.
Appendix A: Eigenvalues of the Coupling Matrix
The adjacency matrix describes the topology of a network.It is usually defined as ˜ A ij = (cid:40) if node i and j are connected otherwise. (A1)where i, j ∈ { , . . . , N } . For the purpose of this study we usea definition where A is normalized to unity row sum A = D − ˜ A , (A2)with the degree matrix DD = diag ( d , d , · · · , d n ) where d i is the degree of node i . A is in general not sym-metric, but D − and (for undirected networks) ˜ A are. Notethat ˜ A is symmetric with only real entries and that D is obvi-ously symmetric, real and positive definite since all eigenval-ues d i > . It therefore can be written in Cholesky decompo-sition as D = LL T where L is real and L T is the transpose of L . Now we write A = D − ˜ A = LL T ˜ A B = L T ˜ AL . (A3)It is similar to A , i.e., there exists a matrix P such that B = P − AP . (A4)It is easy to see that P = L satisfies this condition P − AP = L − (cid:16) LL T ˜ A (cid:17) L = L T ˜ AL = B .
Now we show that two similar matrices have the same eigen-values. If ν is an eigenvector of A with eigenvalue λAν = λν then using A = P BP − obtained from Eq. (A4) we can write P BP − ν = λν ⇒ BP − ν = λP − ν . Thus, if ν is an eigenvector of A with eigenvalue λ , then P − ν is an eigenvector of B with the same eigenvalue λ . Similarly,if Bν = λν we can use B = P − AP and obtain P − AP ν = λν ⇒ AP ν = λP ν . Hence it is true that A and B have the same eigenvalues. Nownote that according to Eq. (A3) B is real and as can be seenfrom B T = (cid:16) L T ˜ AL (cid:17) T = L T A T ( L T ) T = L T AL = B it is also symmetric. This means B and therefore the normal-ized adjacency matrix A has only real eigenvalues.In addition, the Gershgorin circle theorem yields that since (cid:80) j A ij = 1 ∀ i ∈ { , . . . , N } the eigenvalues of A lie insidethe unit circle in the complex plane. In summary: − ≤ η ≤ . (A5) Appendix B: Eigenvalues of the matrix
V AV