Node differentiation dynamics along the route to synchronization in complex networks
Christophe Letellier, Irene Sendiña-Nadal, Ludovico Minati, I. Leyva
NNode differentiation dynamics along the route to synchronization in complexnetworks
Christophe Letellier, a) Irene Sendi˜na-Nadal,
2, 3
Ludovico Minati,
4, 5 and I. Leyva
2, 3 Rouen Normandie University — CORIA, Campus Universitaire du Madrillet, F-76800 Saint-Etienne du Rouvray,France b)2)
Complex Systems Group & GISC, Universidad Rey Juan Carlos, 28933 M´ostoles, Madrid,Spain c)3)
Center for Biomedical Technology, Universidad Polit´ecnica de Madrid, 28223 Pozuelo de Alarc´on, Madrid,Spain Center for Mind/Brain Sciences (CIMeC), University of Trento, 38123 Trento,Italy Institute of Innovative Research, Tokyo Institute of Technology, Yokohama, 226-8503,Japan d) (Dated: 22 February 2021) Synchronization has been the subject of intense research during decades mainly focused on determining thestructural and dynamical conditions driving a set of interacting units to a coherent state globally stable.However, little attention has been paid to the description of the dynamical development of each individualnetworked unit in the process towards the synchronization of the whole ensemble. In this paper, we show howin a network of identical dynamical systems, nodes belonging to the same degree class differentiate in the samemanner visiting a sequence of states of diverse complexity along the route to synchronization independentlyon the global network structure. In particular, we observe, just after interaction starts pulling orbits fromthe initially uncoupled attractor, a general reduction of the complexity of the dynamics of all units beingmore pronounced in those with higher connectivity. In the weak coupling regime, when synchronization startsto build up, there is an increase in the dynamical complexity whose maximum is achieved, in general, firstin the hubs due to their earlier synchronization with the mean field. For very strong coupling, just beforecomplete synchronization, we found a hierarchical dynamical differentiation with lower degree nodes being theones exhibiting the largest complexity departure. We unveil how this differentiation route holds for severalmodels of nonlinear dynamics including toroidal chaos and how it depends on the coupling function. Thisstudy provides new insights to understand better strategies for network identification and control or to deviseeffective methods for network inference.Keywords: Complex networks — Dynamical complexity — Route to synchronization — Chaos
I. INTRODUCTION
The description of networks can be structural, basedon the characterization and modelling of their topologi-cal properties, or dynamical, usually referring to thecollective behavior emerging from the interaction be-tween node dynamics and the network architecture. Syn-chronization is the most thoroughly investigated ensem-ble dynamics, including a rich variety of related be-haviours such as chimera states, cluster synchroniza-tion or explosive synchronization. The route to the synchronous state has been ap-proached from different points of view. Microscopically,it has been analyzed the different ways local synchro-nization grows as the coupling increases depending onthe connectivity structure.
Globally, the synchro-nizability of a population of identical networked oscilla-tors, that is, the stability of the coherent state, has been a) b) Electronic mail: [email protected] c) Electronic mail: [email protected] d) Electronic mail: [email protected] tackled through the master stability function and fore-sees the emergence of patterns when the stability of thenetwork uniform state is lost.
While these approaches help to understand how syn-chronization clusters grow and eventually merge into amacroscopic coherent state or predict its stability for aparticular network structure and coupling function, fewworks pay attention to the description of the nodal dy-namics along the route to synchronization. Recently, ithas been shown that the degree of complexity of thenode dynamics is strongly dependent on the connectivity;in particular, the local connectedness appears to conferconsiderable node differentiation hallmarks, albeit follow-ing mechanisms that are intricate and delineate a non-monotonic effect.
These ensemble fingerprints in thenode dynamics can be used to infer the network statisti-cal description, or even the detailed connectivity in cer-tain conditions. By node differentiation it is here meantas the dynamical departure of an isolated node due tothe interaction with its neighborhood, being more pro-nounced in networks with nodes having degree hetero-geneity. In particular, a relatively large range of cou-pling strengths over which nodes with a higher degreehave a less complex dynamics has been reported, ex- a r X i v : . [ n li n . AO ] F e b plaining the low complexity observed in hubs of func-tional brain networks. However, the reverse situation,the existence of an even lower range of coupling strengthswhere high degree nodes are more complex than lowerdegree ones, remains under study. The purpose of thepresent work is to elucidate the route of the node differ-entiation dynamics in a complex network before reachinga synchronous state, and, in particular, how it dependson the dynamical system, the coupling function and onthe topological properties of the network structure.In order to characterize and quantify this route tonode differentiation dynamics as the coupling strengthincreases, we used the complexity measure recently intro-duced by Letellier et al . able to distinguish organizedfrom disorganized chaotic behavior together with mapsof the dynamical patterns associated to relevant degreeclasses. Section II is devoted to introducing terminologyand providing a brief description of the dynamical com-plexity measure. We then investigated in Section III thedynamical differences among R¨ossler oscillators in a starnetwork according to their degree, and how those differ-ences dependent upon the coupling function, the nomi-nal dynamics, and the size of the star. In Section IV,we extended the analysis to star networks of dynamicalsystems featuring higher nominal complexity, namely thesymmetric Lorenz system, the high-dimensional Mackey-Glass equation, and the toroidal Saito system, to deter-mine whether a general scenario can be outlined. In Sec-tion V, we considered R¨ossler systems in larger networksto explore the node differentiation induced by embeddingthe oscillators in a much more complex topological envi-ronment. Finally, Section VI offers general conclusions. II. DYNAMICAL COMPLEXITY MEASURE
Let us start by considering the general description of anetwork comprising N diffusively-coupled m -dimensionalidentical dynamical systems, whose state vector x i evolves as ˙ x i = f ( x i ) − d N (cid:88) j =1 L ij h ( x j ) (1)where f : R m → R m and h : R m → R m are the nominaldynamics — here, nominal refers to the vector field of the i -th nodal dynamics when isolated — and the couplingfunction, respectively, and d is the coupling strength. L ij are the elements of the Laplacian matrix encoding thenetwork’s connectivity, with L ii = k i the node degree, L ij = − i and j are connected, and L ij = 0otherwise.To characterize the nodal dynamics in a more reliableway, we computed a Poincar´e section to the trajectory x i ( t ) in the state space to rule out the local linear com-ponent and focus on the nonlinear signatures governingthe dynamics. Thus, for each oscillator, in the follow-ing, we computed a two-dimensional Poincar´e section of the attractor when the dynamics has a toroidal structure,or a first-return map built with one of the coordinates ofthe Poincar´e section in case it is non-toroidal. Then, eachoscillator’s map is used to obtain a complexity coefficient C D introduced in Ref. and defined as C D = S p + ∆ , (2)being S p a permutation entropy and ∆ a structuralitymarker. We chose this complexity definition because ofits ability to distinguish organized unpredictable behav-ior (dissipative chaos) from disorganized unpredictablephenomena (noise or conservative chaos). While the en-tropy quantifies the unpredictability of the dynamics, thestructurality accounts for its undescribability, two verydifferent aspects contributing to the complexity of a dy-namics.The permutation entropy S p is defined as in Ref. : S p = − N s !) N s ! (cid:88) π p π log p π ∈ [0 , . (3)It is based on the probability distribution p π of the N s !possible ordinal patterns constructed from the order rela-tions of N s successive data-points in the Poincar´e section.The structurality is computed by dividing the first-returnmap (or the two-dimensional projection of the Poincar´esection) into a N q × N q boxes and∆ = N q (cid:88) i,j =1 q ij N q ∈ [0 , , (4)where q ij = 1 if the box ( i, j ) is visited at least once, and q ij = 0 otherwise. The structurality ∆ is typically low(∆ < .
2) for organized dynamics, and large (∆ > . C D ≈ S p ≈ ≈ . < C D < . C D > . N s = 6 as the lengthof the ordinal patterns. The number N d of data pointsin the Poincar´e section to compute S p and ∆ was N d =10 , N d > N s ! rec-ommended by Riedl et al. The number of boxes was N q = 5 log( N d ) chosen such that N d (cid:29) N q , as pre-scribed in Ref. . The box width δ p is determined bythe largest range visited in our simulations for a givensystem as δ p = ( x max − x min ) /N q , where x max and x min are the maximal and minimal values recorded along oneaxis of the first-return map. To avoid excessively pro-moting noise contamination or small fluctuations aroundperiod-1 limit cycle, we introduced a “noise filter” to in-terpret the ordered patterns of the N s data points toallow permutation only if | x i − x j | > δ p .Let us consider the paradigmatic R¨ossler system toexemplify the ability of Eq.(2) to discriminate betweendifferent dynamics. The vector flow in Eq. (1) is f ( x ) = [ − y − z, x + ay, b + z ( x − c )] (5) FIG. 1. Dynamical characterization of a single R¨ossler systemupon variation of the a parameter and fixing b = 0 . c = 5 .
7. Bifurcation diagram using the Poincar´e section P x (a) and dynamical complexity C D (b). The dashed verticallines in (a) and (b) are located at a = 0 . a = 0 .
31, and a = 0 .
38. The corresponding first-return maps from y n areshown in (c)-(e) and their complexity values at the bottom ofeach panel. taking b = 0 . c = 5 . a ∈ [0 . , .
38] as the bifurcation parameter. Equa-tion (1) was integrated using a fourth-order Runge-Kuttascheme with a time step δt = 0 .
01. Initial conditionsare randomly selected within a small neighborhood ofradius 0 . a for a single R¨ossler system ( d = 0) com-puted from the y n coordinate of the Poincar´e section P x ≡ (cid:8) ( y n , z n ) ∈ R | x n = x − , ˙ x n > (cid:9) , (6)where x − = ( c − √ c − ab ) / The classical period-doubling cas-cade as a route to chaos is quantitatively characterizedby the dynamical complexity C D in Fig. 1(b). It per-fectly captures the periodic windows ( C D < .
2) inter-mingled with chaotic behavior. Panels (c)-(e) show thefirst-return maps for the three a values marked in (a)-(b) with vertical lines. The first value, a = 0 .
26 is lo-cated between a period-3 and a period-2 window and itsfirst-return map is bimodal [Fig. 1(c)] but with a slightlydeveloped third branch appearing just after the period-3 window. For a = 0 .
31, Fig. 1(d), the bimodal maphas the three branches. A more developed chaos, char-acterized with three monotone branches, is observed for a = 0 .
38 in Fig. 1(e). In particular, the right end of thethird branch approaches the vicinity of the bisecting line,indicating that a new period-1 orbit is about to be cre-
FIG. 2. Colormaps of the complexity values for a star of N = 16 coupled R¨ossler oscillators in the d - a parameter plane.Actual complexity change of the hub C hub − C (a)-(b) andthe difference C hub −(cid:104) C leaves (cid:105) (c)-(d) when oscillators are cou-pled through the y (a,c), and x (b,d) variables, respectively.Bottom panels (e)-(f) show the synchronization error E . Thecolor scales are provided on the right sides. In (a)-(f) panels,white curves are the corresponding null isolines, and horizon-tal dashed lines in (a)-(b) are plotted at a = 0 . , .
31 and a = 0 .
38. White domains correspond to the ejection of thetrajectory to infinity. Other parameter values are set as inFig. 1. ated in the population of unstable periodic orbits. Thisis the most developed chaos that can be observed alongthis line of the parameter space, as indicated by the dy-namical complexity ( C D = 0 . III. A STAR NETWORK OF R ¨OSSLER SYSTEMS
Let us now consider a star network of N R¨ossler sys-tems coupled according to Eq. (1), with N − d − a . Top panels (a) and (b) of Fig. 2 showthe actual change in complexity C hub − C of the hub ina N = 16 star with respect to the reference complexityvalue C (which is a function of the parameter a ) of anuncoupled node, for coupling schemes through variables y and x , respectively. Middle panels (c) and (d) providethe complexity difference between the hub and the leaves C hub − (cid:104) C leaves (cid:105) .From the colormaps, the first clear observation is thatthe coupling function has a role in the node differentia-tion dynamics, as reflected by the different distributionof complexity change on the a - d plane, which stronglydepends on the network synchronizability. Panels (e,f)in Fig.2 show the time-averaged synchronization errorcomputed as E = 2 N ( N − (cid:88) i (cid:54) = j (cid:107) x i − x j (cid:107) . (7)which is in agreement with the prediction given by themaster stability function (MSF) for R¨ossler systems cou-pled through the y variable, type I synchronizability class-the MSF becomes negative above a critical value, and x variable, type II class -the MSF is negative in a boundedinterval. For the particular network structure of astar of size N = 16, complete synchronization is indeedreached for a larger number of pairs ( a , d ) when R¨osslersystems are y -coupled than when x -coupled. In the lattercase, full synchronization cannot be reached for a (cid:38) . dλ i , being λ i the i -th eigenvalue of theLaplacian matrix, is outside the range where the MSF isnegative for any value of d , and the synchronous solutionis always unstable.Roughly, in the regions where the synchronization er-ror is not null, for both x and y coupling, there is alarge region of the parameter space with yelowish areasdelimited by the white null isolines in panels in Figs.2(a) and 2(b) where the effect of coupling is to renderthe dynamics of the hub more complex with respect toits uncoupled regime ( C hub > C ). Outside this region,there are smaller islands (in blue) where the dynamicsturns out to be less complex. However, when comparingthis departure from the uncoupled dynamics between thehub and the leaves, Figs. 2(c) and 2(d), there are regionswhere the hub dynamics appears to be more complexthan that of the leaves — for intermediate values of thecoupling and of the parameter a — and regions wherethe opposite is realized. For example, for the y coupling, C hub < (cid:104) C leaves (cid:105) at the left and right side of the yellowregion.To analyse with further detail the dynamical differen-tiation of hub and leaves along the route to synchroniza-tion, we show in Fig. 3 cuts of the complexity difference C j = C j − C , with j = { hub , leaf } (black for the huband red for one leaf) along the three horizontal dashedlines shown in Fig. 2(a). These curves clearly show thathub and leaf experience unalike differentiation routes be-fore either the star synchronizes at d c to the same dy-namical state as in d = 0 for y -coupling [Fig. 3(a)-(d)]or it reaches a boundary crisis at d b for the x coupling[Fig. 3(e)-(h)]. In all cases, the maximal differentiationoccurs always before for the hub than for the leaf, that is,the maximum of the C j curves is located at a couplingthat is lower for the hub than for the leaf. As it is shownat the inset of each panel, this shift coincides with theearlier synchronization of the hub to the mean field, heremeasured as S j = (cid:104) Re( e i( θ j − Φ) ) (cid:105) t , (8) -1-0.500.5(a) hubleaf -1-0.500.5 ABC D E F (b) -1-0.500.5(c) 0 0.1 0.2 0.3 0.4-1-0.500.5(d) -1-0.500.5(e) hubleaf -1-0.500.5(f) -1-0.500.5(g) 0 0.1 0.2 0.3 0.4-1-0.500.5(h)
FIG. 3. Evolution of the relative complexity C j as a functionof the coupling strength d for the hub (in black) and a leaf(in red) of a star network of y coupled (a)-(d) and x coupled(e)-(h) R¨ossler oscillators. (a,e) a = 0 .
26, (b,f) a = 0 .
31, and(c,g) a = 0 .
38 and N = 16. (d,h) a = 0 .
31 and N = 32.Insets in all panels show the phase synchronization S j of thehub (black) and the leaf (red) with the mean field. Otherparameter values: b = 0 .
2, and c = 5 .
7. One run per d -value. being (cid:104) . . . (cid:105) t the time average, Re stands for the real part, θ j = arctan ( y j /x j ) the phase of the j = hub / leaf, andΦ the mean-field global phase. By definition, S j ∈ [0 , j -th oscillator islocked to the phase of the mean field.Another observed systematic behavior is that, for in-creasing values of a , as the uncoupled dynamics is moredeveloped (see Fig. 1), the relative complexity of the hubdiminishes while increases for the leaf. Therefore, the rel-ative position between the hub and leaf curves changessuch that for a = 0 .
26 the hub’s curve is almost alwaysabove the leaf’s one, C hub , > C leaf , , while for a = 0 . C hub , < C leaf , . Anothersalient feature is that in the very weak coupling regime( d < . C forboth hub and leaf. Finally, the impact of the size of thestar is analysed in Fig. 3(d,h) for N = 32 and a = 0 . N = 16. Precisely, the hub com-plexity is slightly affected and only in the region of lowcoupling. As predicted by the MSF, for the y coupling,the threshold for synchronization d c does only depend onthe smallest non-zero eigenvalue λ , which equals λ = 1in both configurations. On the other hand, for the x cou-pling, the boundary crisis is anticipated to smaller cou-pling values as the largest eigenvalue λ N = N increaseswith the size.To picture the different dynamical states the hub andthe leaf are visiting, we monitored their first-return mapsat a specific points along the route to synchronization asshown in Fig. 3(b) with vertical dashed lines. Having inmind the map characterizing the uncoupled dynamics for a = 0 .
31 and shown in Fig. 1(d), for very low coupling[Fig. 4(A)], the leaf dynamics is nearly unaffected whilethe hub is exhibiting the presence of incoherent smallperturbations and a shorter third branch. A slight in-crease in the coupling strength [Fig. 4(B)] amplifies thiseffect turning the hub’s attractor into a very thick uni-modal map and leaving the leaf’s map almost unaltered.The next scenario [Fig. 4(C)] corresponds to the largestdeviation from the uncoupled dynamics with both typesof nodes displaying a much less complex dynamics. Thehub is constrained to a small neighborhood of the innersingular point as revealed by the feeble cloud of points inthe bisecting line. The leaf is locked on a period-2 limitcycle. Beyond this drop of complexity, the hub dynam-ics develops into a single huge cloud of points, displacedalong the bisecting line towards the period-1 limit cy-cle [Fig. 4(D)]. The large complexity value attained bythe hub C hub = 1 . > d > d c in Fig. 4(F) and all nodes share the same dynamics asthey had when uncoupled. IV. STAR NETWORKS MADE OF OTHERDYNAMICAL SYSTEMSA. Lorenz systems
Let us now consider the Lorenz 63 system whose vec-tor flow in Eq. (1) is f ( x ) = [ σ ( y − x ) , Rx − y − xz, − bz + xy ] (9)which is equivariant under a rotation symmetry aroundthe z -axis . This global property is the main differencewith respect to the R¨ossler dynamics since, when thesymmetry is modded out, the Lorenz attractor is topolog-ically equivalent to the R¨ossler one. Lorenz systemscoupled through variable y present a type- I synchroniz-ability class and, therefore, complete synchronization isstable above a critical coupling threshold.To compute the first-return map we proceeded as fol-lows. The common Lorenz attractor for R = 28, re-sembling the two wings of a butterfly, is bounded by agenus-3 torus and the Poincar´e section is the union of thetwo components P ± ≡ (cid:8) ( y n , z n ) ∈ R | x n = ± x ± , ˙ x n ≶ (cid:9) (10)where x ± = ± (cid:112) b ( R − The interval visited byeach one of these components is normalized to the unitinterval, and the component P − is shifted by −
1, leadingto variable ρ . A quite close example of the first-returnmap for R = 28 is shown at the bottom of Fig. 5A featur-ing four branches paired due to the rotation symmetry. The two increasing branches correspond to the reinjec-tion of the trajectory into the wing from which it is is-sued, while the two decreasing branches are associatedwith the transition from one wing to the other. The twoleft (right) branches correspond to the n th intersectionin the left (right) wing and to the ( n + 1)th intersectionin the left (right) or right (left) wing depending on thesign of the slope.Figure 5 shows the relative complexity C j of the huband a leaf of a N = 32 star network y -coupled Lorenzsystems for two values of the parameter R . For R = 28,full synchronization is reached for d c ≈ .
22. As with theR¨ossler systems, the complexity of the hub and leavesdrop below C in the weak coupling regime ( d (cid:46)
1) andbeyond that point the hub becomes more complex thanthe nominal dynamics before full synchronization is even-tually reached. For this parameter setting, the Lorenzhub experiences a sudden drop in complexity as illus-trated in the first-return map [Fig. 5(A)]: while the hubdynamics is characterized by small fluctuations aroundthe singular point in the centre of one of the wings, theleaves are nearly unperturbed with their four branches. A: d = 0.004 B: d = 0.014 C: d = 0.024 D: d = 0.100 E: d = 0.240 F: d = 0.300 H ub L ea f C D = 0.65 C D = 0.83 C D = 0.72 C D = 0.89 C D = 0.02 C D = 0.12 C D = 1.18 C D = 0.98 C D = 1.07 C D = 1.27 C D = 0.82 C D = 0.82 FIG. 4. First-return maps to the Poincar´e section P for the hub (top row) and one of the leaves (bottom row) in a star networkof N = 16 y -coupled R¨ossler oscillators for a = 0 .
31 and for different d values. The corresponding complexities C D are reported.The nominal complexity is C = 0 . C j HubLeaf0 1 2 3 4 5 6
Coupling strength d -0,500,5 C j R = R = A B d c d c C (a)(b) ρ n ρ n + A : d = 0.20 ρ n ρ n + B : d = 1.06 ρ n ρ n + C : d = 1.20 ρ n ρ n + ρ n ρ n + ρ n ρ n + H ub L ea f C D = 0.03 C D = 0.02 C D = 0.65 C D = 0.40 C D = 0.02 C D = 0.50 FIG. 5. (a)-(b) Relative complexity for the hub (black) andone of the leaves (red) as a function of the coupling strength d for a star network with N = 32 y -coupled Lorenz systemsfor two R -values. First-return maps of the hub (top) and ofone leaf (bottom) for the cases A, B and C indicated by theblue dashed lines in the top panel. Other parameter values: σ = 10, b = . A magnified view of the maps is plotted forthe case B: the range used is not the same between the huband leaf. When d is further increased [Fig. 5(B)], the leaf dynam-ics also collapses, with all nodes surprisingly exhibiting aquasi-periodic dynamics, a regime not observed in an iso-lated Lorenz system. Notice that the size of the two toriare different. When the nodes’ dynamics have topologi-cally equivalent attractors (here, tori) but with a scalingfactor, the phenomenon is known as amplitude enveloppesynchronization . Finally, beyond this point of reduceddynamics and just before full synchronization [Fig. 5(C)],the interaction between hub and leaves moves the dy-namics of the leaves again to the four paired branches al-though slightly more noisy than the uncoupled one, whilethe hub is strongly perturbed with a very disorganizeddynamics. This scenario resembles the one observed forthe R¨ossler system depicted in panel D of Fig. 4. Fi-nally, we explored a second regime with an even moredeveloped uncoupled dynamics for R = 175 whose first-return map has six branches (not shown). As in Fig. 3(c)for the R¨ossler system, the switch between the relativecomplexities of hub and leaf is no longer observed in Fig.5 and the leaf curve is always above C leaf , > C hub , . B. Mackey-Glass delay differential equation
It is possible to produce a chaotic attractor character-ized by a smooth unimodal map via one-dimensional de-lay differential equations, for instance, the Mackey-Glass(MG) equation ˙ x = µ x τ x pτ − x t (11)where x is the population of blood cells, x t = x ( t ) and x τ = x ( t − τ ). This equation was initially proposed forthe control of hematopoiesis (the production of bloodcells). Typically, the delay τ is the time-scale for pro-liferation and maturation of these blood cells. The di- FIG. 6. (a) Bifurcation diagram for the Mackey-Glass delaydifferential equation (11) versus the time delay τ . (b-c) Bifur-cation diagrams computed versus the coupling strength d forthe hub and one of the leaves of a star network of N = 32 cou-pled Mackey-Glass equations (11) for τ = 4 . µ = 1 . p = 18 . mension of the effective state space associated with a de-lay differential equation is dependent on the delay. Parameter values for parameters µ and p are such thata period-doubling cascade as a route to chaos is ob-served [Fig. 6(a)]. We chose a delay τ = 4 . τ further leads toa much more complex behavior and computing a reliablePoincar´e section is rather tricky.Our goal here is to investigate the route to synchroniza-tion for unimodal dynamics produced by a potentiallyhigh-dimensional system. A network of MG systems canbe fully synchronized using bidirectional coupling. The attractor is bounded by a genus-1 torus and thesingle-component Poincar´e section can be defined as P MG ≡ { x n , ∈ R | ˙ x n = 0 . , ¨ x n > , x n < . } . (12)The first-return map is one-dimensional and slightly fo-liated and quite similar to the one shown for the leaf inFig.7(A) for N = 16 coupled MG systems. Figure 7(a)shows the relative complexity of the hub and leaf as afunction of the coupling strength. As observed in theR¨ossler and Lorenz systems, for low d -values, the hubexperiences a sudden complexity drop while the leaveskeep almost their nominal dynamics. As the coupling isincreased, the hub starts synchronizing with the meanfield and its relative complexity rises above 0. Beyondpoint B, all nodes reduce their complexity down to a Coupling strength d -1-0.500.5 C j HubLeaf d c A B C (a) x n x n + A : d = 0.011 x n x n + B : d = 0.038 x n x n + C : d = 0.150 x n x n + x n x n + x n x n + H ub L ea f C D = 0.47 C D = 0.86 C D = 0.22 C D = 0.23 C D = 0.73 C D = 0.70 FIG. 7. (a) Relative complexity as a function of the couplingstrength d for a star network of N = 16 coupled Mackey-Glassdelay differential equations (11) for τ = 4 .
8. (A)-(C) First-return maps of the hub (top) and of one leaf (bottom) forthree different values of d . Other parameter values: µ = 1 . p = 18 . minimal value (point C) before finally increasing up tofull synchronization is reached. The noisy curves in thatregion is due to the extreme sensitivity to initial condi-tions found in this system. The maps for the hub andleaf [Figs. 7(A)-7(C)] corresponding to the setting pointsA,B, and C marked in Fig. 7(a) sketch the dynamicaldifferentiation route.The bifurcation diagrams for the hub and leaf respec-tively are computed as a function of the coupling d [Figs.6(b)-6(c)]. In both cases, an inverse period-doubling cas-cade is observed up to a period-2 limit cycle is settled at d = 0 .
18. Curiously, the size of both limit cycles is differ-ent which, again, is an example of an amplitude envelopesynchronization. The crisis leading to a larger chaotic at-tractor is strongly dependent on the initial conditions. Insummary, the prevalent lines of the route to synchroniza-tion previously sketched are visible observed, albeit withsome differences. Further investigations are necessary tofully understand their origin.
C. Four-dimensional Saito model
In this section, to confirm the generality of the nodaldynamical differentiation route, we will investigate acompletely different model of nonlinear dynamics, theSaito model, which is able to generate a very rich dy-namical behavior including toroidal chaos. The vector Coupling strength d -1-0.500.5 C j -1-0.500.5 C j HubLeaf d c d c δ = . δ = . AB (a)(b) x n z n x n z n x n z n x n z n C Hub = 0.54 C Leaf = 0.84 C Leaf = 0.96 C Hub = 0.74 A : δ = . , d = . B : δ = . , d = . Hub Leaf
FIG. 8. Relative complexity C k as a function of the couplingstrength d of a star network with N = 16 y -coupled Saitomodels (13) for (a) δ = 0 .
58 and (b) δ = 0 .
82. Poincar´esections for the hub and a leaf for (A) δ = 0 .
58 and for (B) δ = 0 .
82, both with d = 0 .
07. Other parameter values: ρ =14, η = 1, and (cid:15) = 0 . flow in Eq. (1) reads f ( x ) = [ ρ ( − y + z ) , x + 2 δy, − x − w, ( z − Φ( w )) /(cid:15) ] (13)where x = ( x, y, z, w ) andΦ( w ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) w − (1 + η ) w ≥ η − wη if | w | < ηw + (1 + η ) w ≤ − η . (14)This is a four-dimensional system involving a linear piece-wise function as a switch mechanism. The dynamicsproduced by this model can be chaotic, quasi-periodic,toroidal chaotic, or even hyperchaotic, also structuredaround a torus. Parameters ρ , η , and (cid:15) are fixed and δ isused as a bifurcation parameter. Here we chose δ = 0 .
58 for toroidal chaos, and δ = 82 to produce hyperchaotictoroidal chaos. Close examples of the Poincar´e sectionsof these attractors can be grasped, respectively, in Fig.8(A) and 8(B) for a leaf of a star of N = 16 y -coupledSaito models whose dynamics are almost identical to theuncoupled scenarios. The main difference between thesetwo dynamics is the “thickness” of the Poincar´e sectionbeing much thiner for δ = 0 .
58 than for δ = 0 .
82. Thehyperchaotic nature is revealed by the overlapping struc-tures of the Poincar´e section as observed in the folded-towel map introduced by R¨ossler. This is partly con-firmed with the Lyapunov exponents which for δ = 0 . λ = 0 . > λ = 0 . ≈ λ = − . > λ = − . , with two null exponents as expected for toroidal chaosstructured around a torus T , and for δ = 0 .
82 are λ = 0 . > λ = 0 . ≈ λ = − . > λ = − . , with two positive and one null exponents as needed fortoroidal hyperchaos. Note that in the latter case it is stillunclear whether the second positive exponent is mergedor not with the third null exponent and this is an issuecurrently under study. Investigations, which are out ofthe scope of the present work, would allow determiningwhether the fourth dimension is required for embeddingthe toroidal chaos ( δ = 0 . δ = 0 .
82, the dynamics are necessarily four-dimensional.When these models are coupled through the variable y , they synchronize as shown in Figs. 8(a) and 8(b) withthe convergence of their relative complexity to zero, be-ing the critical coupling d c larger for the hyperchaoticdynamics. Nevertheless, in both cases, what is preservedalong the route to synchronization with these toroidalchaotic and hyperchaotic dynamics is the node differen-tiation mainly of the hub whose dynamics turns less de-veloped than the uncoupled one while the leaves sustainit over the whole range of coupling strengths [comparethe first-return maps in Fig. 8(A), and in Fig. 8(B)]. Theroute to synchronization appears very simple, most likelydue to the constrained toroidal structure of the nominaldynamics. V. NETWORKS OF R ¨OSSLER OSCILLATORS
Having characterized the relationship between the de-gree centrality and dynamical complexity in star net-works, we move forward to generalize our results to net-works with a broader degree distribution. This issue waspartially tackled in Refs. , where a strong correlationbetween degree k and complexity C D allowed establishinga node hierarchy. However, given the present novel re-sults revealing the important role played by the nominalnodal dynamics, we extend our present study to largernetworks.We maximize the degree heterogeneity by usingBarabasi-Albert scale-free (SF) networks of N identical S k a=0.26 (a) k=2k=5k=10k=15k=24 a=0.31 (b)0 0.2 0.4 0.6 d -0.6-0.4-0.200.20.4 C k (c) 0 0.2 0.4 0.6 d (d)0.7 0.8 0.9 1 S k -0.500.5 C k (e) 0.7 0.8 0.9 1 S k (f) FIG. 9. (a,b) Average k -class S k phase synchronization pa-rameter for several values of degree k for N = 300, (cid:104) k (cid:105) = 4SF networks of y -coupled R¨ossler systems when a = 0 .
26 (leftpanels) and a = 0 .
31 (right panels), (c,d) averaged relativecomplexity C k . The vertical dotted lines mark the cou-plings analyzed in Fig. 10 the lower panels: d = 0 .
04 (black), d = 0 .
08 (blue) and d = 0 .
22 (magenta). (e,f) Relative com-plexity for different k -classes as a function of the phase syn-chronization S k . y -coupled R¨ossler oscillators retaining the parameter set-tings provided in Section III. Since we expect that nodeshaving the same degree k play equivalent roles in the net-work, we calculate the evolution of C k within a degreeclass k by averaging over the N k nodes having degree k ,that is, C k = 1 N k (cid:88) [ j | k j = k ] C j , (15)where C j is the dynamical complexity of the j th node.Here, we use the relative complexity C k = C k − C whichhelps to better assess the effects of both the coupling andthe topology in the complexity of the k -class nodes. Inaddition, to evaluate the impact of the nodal environ-ment, we perform our calculations for networks made of N = 300 nodes, wherein (cid:104) k (cid:105) = 4. All the results areaveraged over 10 different networks realizations.First, in Fig. 9(a-b) we plot the time averaged phasesynchronization S k of k -class nodes with respect to thephase of the mean field for N = 300 SF networks as afunction of the coupling d . We define S k as S k = 1 N k N (cid:88) [ j | k j = k ] S j . (16)where S j is the time averaged phase synchronizationof node j with the mean-field global phase, defined in Eq. (8). The result is then ensemble averaged over 10network realizations. Along the route to synchronization,the k -classes synchronize hierarchically to the mean fieldbefore all of them lock at the critical coupling ( d c = 0 . a = 0 .
26 and d c = 0 .
52 for a = 0 . therefore existing node differentiation also in heteroge-neous networks. While the phase synchronization S k ismonotonously increasing for most of the coupling rangeand classes, the weakly coupled regime presents anoma-lous synchronization over which most of the S k arebelow the basal, uncoupled level. This anomalous range,more prominent for a = 0 .
31, is associated with the max-imal node differentiation [compare Fig. 9(b) and 9(d)].As already observed in Section III, the less developedthe dynamics, the smaller the critical value d c [Fig. 9(c)-(d)]. In the weakly coupled regime, the relative com-plexity C k shows a strong node differentiation: whilethe hubs present a markedly reduced complexity withrespect to the nominal value (with a clear negative min-imal), for the smaller degrees the complexity increases.This increment in the less connected nodes ( k = 2 and k = 5 in the example) is non-monotonous with d and k ,as observed in the stars in Sections III and IV.For stronger coupling all the nodes increase their com-plexity well above the nominal value, ordered followingthe reverse degree ranking [Fig. 9(c)-(d)], recovering thescenario observed for the hub in Section III. All thesefeatures are qualitatively shared for the two different a -values, but the deviations from the uncoupled value arelarger for the more developed chaos, a = 0 . C k as a function ofthe phase synchronization S k reveals a dependency whichis stronger for nodes with a larger degree [Fig. 9(e)-(f)].Typically, low-degree nodes exhibit a dynamics which isnearly independent of the synchronization while it is theopposite for large-degree nodes. It also clearly showsthat the relative complexity converges to zero for larger d -values when k increases. This delineates another sig-nature of node differentiation.Therefore, we conclude that in most of the regimes itis possible to correlate the node degree with the relativecomplexity. Furthermore, the degree centrality is the sin-gle structural parameter that affects the node behaviour,while the rest of environmental topological features hasno impact. This is shown in Fig. 10, where we plot thevalue of the relative complexity C k as a function of k for three representative values of d , both for SF networksand ER networks with (cid:104) k (cid:105) = 4. The ER and SF curvesoverlap and, therefore, the dependence of C k on a and d are the same regardless of the topology. This is quiteremarkable since the ER and SF networks have a differ-ent critical coupling d c : for a given d -value, their globaldynamics are different, but the nodes of degree k haveequivalent dynamics independently of the environment.The same result is obtained for different sizes of both ERand SF networks (not shown).To better illustrate the node differentiation in largernetworks, we plotted the first-return maps for three dif-0 k -0.500.5 C k (a)(a) SF d=0.04ER d=0.04SF d=0.08ER d=0.08SF d=0.22ER d=0.22 k (b)(b) FIG. 10. Relative complexity C k as a function of the de-gree k for (a) a = 0 .
26, and (b) a = 0 .
31. In each panel,curves correspond to SF (void symbols) and ER (full sym-bols) with the coupling strengths marked in Fig.10(c,d) withvertical lines: d = 0 .
04 (circles), d = 0 .
08 (squares), and d = 0 .
22 (triangles). Results are averaged over 10 differentnetwork realizations of N = 300. Other parameter values asin Fig. 1. ferent k -classes of nodes along the route to synchroniza-tion (Fig. 11). Low-degree nodes ( k = 2) in SF networksproduce first-return maps whose thickness increases withthe coupling strength, up to d < d c (top row in Fig.11): for these nodes, the relative complexity is alwayspositive, and the first maps resemble the map of the un-coupled dynamics but thicker (compare with the map for a = 0 .
31 in Fig. 1). For large degree nodes (third rowin Fig. 11), once the minimum complexity is reached inthe weakly coupled regime, the complexity increases withthe coupling strength d . Around the minimum, the mapscomprise a small cloud of points in the neighborhood ofthe inner singular point (bottom left of the first-returnmap); this is progressively transformed into a “noisy”period-1 limit cycle, which is characterized by a cloudof points elongated perpendicularly to the bisecting lineand located around the centre of the map ( d = 0 .
08 with C = 0 . d = 0 .
22 with C = 0 . k = 8 in the ex-ample) produce a map which has features intermediatebetween the two extreme cases previously discussed: thenode differentiation is, thus, evidently correlated with thenode degree. As previously discussed, in large networks,node differentiation is mostly a monotonous function ofthe coupling strength and of the degree. When SF net-works are replaced by ER ones, the degree range is nar-rower, but the maps for k = 2 and k = 8 are very similarto their counterparts in SF networks. We conclude thatnode degree is clearly the most important factor for thenode differentiation along the route to synchronization. VI. CONCLUSION
Routes to synchronization need to be elucidated to at-tain a better understanding and knowledge of the possi-ble scenarios that may be encountered depending on thecoupling, nominal node dynamics, topology, and network k = d = 0 : C =0.94 d = 0 : C =1.00 d = 0 : C =1.11 k = C =0.67 C =0.92 C =1.00 k = C =0.17 C =0.48 C =0.88 k = C =0.98 C =1.06 C =1.19 k = C =0.91 C =0.94 C =1.04 S F n e t w o r k s E R n e t w o r k s FIG. 11. Different first-return maps to a Poincar´e section forthree k -classes of nodes along a route to synchronization ina N = 300 ( (cid:104) k (cid:105) = 4) SF network of y -coupled R¨ossler nodes( a = 0 .
31 and other parameters as in Fig. 1). The complexity C D is reported in each case. size. Here, we confirmed and extended previous work de-picting a non-trivial effect of connectedness on node dy-namics, particularly the existence of a non-monotonic re-lationship between the complexity of node dynamics andcoupling strength. There is, indeed, a node differentia-tion that evolves with the coupling strength. Typically,when the coupling function provides type- i synchroniz-ability, low coupling strengths induce a significant reduc-tion in the complexity of the large-degree nodes, whilethose having a small degree are left nearly unaffected.Consequently, increasing the d -value, all the nodes insmall networks or those with a large degree in large net-works present a minimal dynamical complexity. In everynetwork, all nodes have a complexity that increases, of-ten reaching a greater level than the nominal one, beforethe onset of full synchronization. This sketch for theroute to synchronization is clearly observed for the two1three-dimensional systems (R¨ossler and Lorenz). Withmore complex node dynamics (Mackey-Glass and Saito),the decrease towards a minimal complexity is only ob-served in the hub of star networks; further studies withother network types are still needed for attaining a moregeneral view of the latter system. The node differenti-ation is not a monotonic function of the coupling norof the synchrony. When the coupling provides a type- ii synchronizability, the route to synchronization is anabridged version of the route observed with a type- i syn-chronizability. ACKNOWLEDGMENTS
ISN and IL acknowledge financial support from theMinisterio de Econom´ıa, Industria y Competitividad ofSpain under project FIS2017-84151-P. M. E. J. Newman, “The structure and function of complex net-works,” SIAM Review , 167–256 (2003). S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U.Hwang, “Complex networks: Structure and dynamics,” PhysicsReports , 175–308 (2006). S. Boccaletti, G. Bianconi, R. Criado, C. del Genio, J. G´omez-Gard˜nes, M. Romance, I. Sendi˜na-Nadal, Z. Wang, andM. Zanin, “The structure and dynamics of multilayer networks,”Physics Reports , 1–122 (2014). A. Arenas, A. D´ıaz-Guilera, and C. J. P´erez-Vicente, “Synchro-nization processes in complex networks,” Physica D , 27–34(2006). F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, “TheKuramoto model in complex networks,” Physics Reports ,1–98 (2016). Y. Kuramoto and D. Battogtokh, “Coexistence of coherence andincoherence in nonlocally coupled phase oscillators,” NonlinearPhenomena in Complex Systems , 380–385 (2002). D. M. Abrams and S. H. Strogatz, “Chimera states for coupledoscillators,” Physical Review Letters , 174102 (2004). A. M. Hagerstrom, T. E. Murphy, R. Roy, P. H¨ovel,I. Omelchenko, and E. Sch¨oll, “Experimental observation ofchimeras in coupled-map lattices,” Nature Physics , 658–661(2012). L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy,and R. Roy, “Cluster synchronization and isolated desynchro-nization in complex networks with symmetries,” Nature Com-munications , 4079 (2014). J. G´omez-Gardenes, S. G´omez, A. Arenas, and Y. Moreno, “Ex-plosive synchronization transitions in scale-free networks,” Phys-ical Review Letters , 128701 (2011). I. Leyva, I. Sendi˜na-Nadal, J. Almendral, A. Navas, S. Olmi, andS. Boccaletti, “Explosive synchronization in weighted complexnetworks,” Physical Review E , 042808 (2013). J. G. Restrepo, E. Ott, and B. R. Hunt, “Synchronization inlarge directed networks of coupled phase oscillators,” Chaos ,015107 (2006). J. G´omez-Garde˜nes, Y. Moreno, and A. Arenas, “Paths to syn-chronization on complex networks,” Physical Review Letters ,034101 (2007). T. Pereira, S. van Strien, and M. Tanzi, “Heterogeneously cou-pled maps: hub dynamics and emergence across connectivity lay-ers,” Journal of the European Mathematical Society , 2183–2252 (2020). L. M. Pecora and T. L. Carroll, “Master stability functionsfor synchronized coupled systems,” Physical Review Letters ,2109–2112 (1998). M. Barahona and L. M. Pecora, “Synchronization in small-worldsystems,” Physical Review Letters , 054101 (2002). L. M. Pecora, “Synchronization of oscillators in complex net-works,” Pramana , 1175–1198 (2008). A. Tlaie, I. Leyva, R. Sevilla-Escoboza, V. P. Vera-Avila, andI. Sendi˜na Nadal, “Dynamical complexity as a proxy for the net-work degree distribution,” Physical Review E , 012310 (2019). A. Tlaie, L. M. Ballesteros-Esteban, I. Leyva, and I. Sendi˜na-Nadal, “Statistical complexity and connectivity relationship incultured neural networks,” Chaos, Solitons & Fractals , 284–290 (2019). L. Minati, H. Ito, A. Perinelli, L. Ricci, L. Faes, N. Yoshimura,Y. Koike, and M. Frasca, “Connectivity influences on nonlin-ear dynamics in weakly-synchronized networks: Insights fromR¨ossler systems, electronic chaotic oscillators, model and biolog-ical neurons,” IEEE Access , 174793–174821 (2019). D. Eroglu, M. Tanzi, S. van Strien, and T. Pereira, “Reveal-ing dynamics, communities, and criticality from data,” PhysicalReview X , 021047 (2020). J. H. Mart´ınez, M. E. L´opez, P. Ariza, M. Chavez, J. A. Pineda-Pardo, D. L´opez-Sanz, P. Gil, F. Maest´u, and J. M. Buld´u,“Functional brain networks reveal the existence of cognitive re-serve and the interplay between network topology and dynam-ics,” Scientific reports , 10525 (2018). C. Letellier, I. Leyva, and I. Sendi˜na Nadal, “Dynamical com-plexity measure to distinguish organized from disorganized dy-namics,” Physical Review E , 022204 (2020). C. Letellier, “Estimating the Shannon entropy: Recurrence plotsversus symbolic dynamics,” Physical Review Letters , 254102(2006). C. Bandt and B. Pompe, “Permutation entropy: A natural com-plexity measure for time series,” Physical Review Letters ,174102 (2002). M. Riedl, A. M¨uller, and N. Wessel, “Practical considerations ofpermutation entropy,” European Physical Journal Specical Top-ics , 249–262 (2013). O. E. R¨ossler, “An equation for continuous chaos,” Physics Let-ters A , 397–398 (1976). W. H. Press, S. A. Teukolsky, and W. T. V. B. P. Flannery,
Numerical Recipes in C. The Art of Scientific Computing , 2nded. (Cambdridge University Press, Cambridge, New York, PortChester, Melbourne, Sydney, 1992). C. Letellier, P. Dutertre, and B. Maheu, “Unstable periodicorbits and templates of the R¨ossler system: Toward a systematictopological characterization,” Chaos , 271–282 (1995). L. Huang, Q. Chen, Y.-C. Lai, and L. M. Pecora, “Genericbehavior of master-stability functions in coupled nonlinear dy-namical systems,” Physical Review E , 036204 (2009). I. Sendi˜na-Nadal, S. Boccaletti, and C. Letellier, “Observabilitycoefficients for predicting the class of synchronizability from thealgebraic structure of the local oscillators,” Physical Review E , 042205 (2016). E. N. Lorenz, “Deterministic nonperiodic flow,” Journal of theAtmospheric Sciences , 130–141 (1963). C. Letellier and R. Gilmore, “Covering dynamical systems: Two-fold covers,” Physical Review E , 016206 (2001). C. Letellier, P. Dutertre, and G. Gouesbet, “Characterization ofthe Lorenz system, taking into account the equivariance of thevector field,” Physical Review E , 3492–3495 (1994). C. Letellier,
Caract´erisation topologique et reconstruction des at-tracteurs ´etranges , Ph.D. thesis, University of Paris VII, Paris,France (1994). T. D. Tsankov and R. Gilmore, “Topological aspects of the struc-ture of chaotic attractors in R ,” Physical Review E , 056206(2004). C. Letellier, L. A. Aguirre, and J. Maquet, “Relation betweenobservability and differential embeddings for nonlinear dynam-ics,” Physical Review E , 066213 (2005). G. Byrne, R. Gilmore, and C. Letellier, “Distinguishing betweenfolding and tearing mechanisms in strange attractors,” Physical Review E , 056214 (2004). J. M. Gonzalez-Miranda, “Amplitude envelope synchronizationin coupled chaotic oscillators,” Physical Review E , 036232(2002). M. C. Mackey and L. Glass, “Oscillation and chaos in physiolog-ical control systems,” Science , 287–289 (1977). L. Glass and M. C. Mackey, “Pathological conditions resultingfrom instabilities in physiological control system,” Annals of theNew York Academy of Sciences , 214–235 (1979). I. Gumowski, “Sensitivity of certain dynamic systems with re-spect to a small delay,” Automatica , 659–674 (1974). J. D. Farmer, “Chaotic attractors of an infinite-dimensional dy-namical system,” Physica D , 366–393 (1982). K. Pyragas, “Synchronization of coupled time-delay systems:Analytical estimations,” Physical Review E , 3067–3071(1998). T. Saito, “An approach toward higher dimensional hysteresischaos generators,” IEEE Transactions on Circuits and Systems , 399–409 (1990). O. E. R¨ossler, “Chaos,” in
Structural Stability in Physics , editedby G. G¨uttinger and H. Eikemeier (Springer, Berlin Heidelberg,1979) pp. 290–309, proceedings of Two International Symposiaon Applications of Catastrophe Theory and Topological Conceptsin Physics T¨ubingen, May 2-6 and December 11-14, 1978. C. Letellier and O. E. R¨ossler, “Chaos: The world of nonperiodicoscillations,” (Springer, Cham, Switzerland, 2020) Chap. Anupdated hierarchy of chaos, pp. 181–203. C. Zhou and J. Kurths, “Hierarchical synchronization in complexnetworks with heterogeneous degrees,” Chaos , 015104 (2006). T. Pereira, “Hub synchronization in scale-free networks,” Physi-cal Review E , 036201 (2010). B. Blasius, E. Montbri´o, and J. Kurths, “Anomalous phase syn-chronization in populations of nonidentical oscillators,” PhysicalReview E , 035204 (2003). B. Boaretto, R. Budzinski, T. Prado, J. Kurths, and S. Lopes,“Neuron dynamics variability and anomalous phase synchroniza-tion of neural networks,” Chaos28