Synchronization in networks with multiple interaction layers

See allHide authors and affiliations

Science Advances  16 Nov 2016:
Vol. 2, no. 11, e1601679
DOI: 10.1126/sciadv.1601679


The structure of many real-world systems is best captured by networks consisting of several interaction layers. Understanding how a multilayered structure of connections affects the synchronization properties of dynamical systems evolving on top of it is a highly relevant endeavor in mathematics and physics and has potential applications in several socially relevant topics, such as power grid engineering and neural dynamics. We propose a general framework to assess the stability of the synchronized state in networks with multiple interaction layers, deriving a necessary condition that generalizes the master stability function approach. We validate our method by applying it to a network of Rössler oscillators with a double layer of interactions and show that highly rich phenomenology emerges from this. This includes cases where the stability of synchronization can be induced even if both layers would have individually induced unstable synchrony, an effect genuinely arising from the true multilayer structure of the interactions among the units in the network.

  • Multi-layer networks
  • synchronization
  • stability
  • master stability function


Network theory (19) has proved a fertile ground for the modeling of a multitude of complex systems. One of the main appeals of this approach lies in its power to identify universal properties in the structure of connections among the elementary units of a system (1012). In turn, this enables researchers to make quantitative predictions about the collective organization of a system at different length scales, ranging from the microscopic to the global scale (1319).

Because networks often support dynamical processes, the interplay between the structure and the unfolding of collective phenomena has been the subject of numerous studies (2022). Many relevant processes and their associated emergent phenomena, such as social dynamics (23), epidemic spreading (24), synchronization (25), and controllability (26), have been proven to significantly depend on the complexity of the underlying interaction backbone. Synchronization of systems of dynamical units is a particularly noteworthy topic because synchronized states are at the core of the development of many coordinated tasks in natural and engineered systems (2729). Thus, in the past two decades, considerable attention has been paid to shedding light on the role that network structure plays in the onset and stability of synchronized states (3042).

However, in the past years, the limitations of the simple network paradigm have become increasingly evident, as the unprecedented availability of large data sets with ever-higher resolution levels has revealed that real-world systems can seldom be described by an isolated network. Several works have proved that mutual interactions between different complex systems cause the emergence of networks composed of multiple layers (4346). This way, nodes can be coupled according to different kinds of ties so that each of these interaction types defines an interaction layer. Examples of multilayer systems include social networks, in which individual people are linked and affiliated by different types of relations (47), mobility networks, in which individual nodes may be served by different means of transport (48, 49), and neural networks, in which the constituent neurons interact over chemical and ionic channels (50). Multilayer networks have thus become the natural framework to investigate new collective properties arising from the interconnection of different systems (51, 52). The multilayer studies of processes, such as percolation (5357), epidemics spreading (5861), controllability (62), evolutionary games (6366), and diffusion (67), have all evidenced a very different phenomenology from the one found on monolayer structures. For example, whereas isolated scale-free networks are robust against random failures of nodes or edges (68), interdependent ones are very fragile (69). Nonetheless, the interplay between multilayer structure and dynamics remains, under several aspects, still unexplored, and in particular, the study of synchronization is still in its infancy (7073).

Here, we present a general theory that fills this gap and generalizes the celebrated master stability function (MSF) approach in complex networks (30) to the realm of multilayer complex systems. Our aim is to provide a full mathematical framework that allows one to evaluate the stability of a globally synchronized state for nonlinear dynamical systems evolving in networks with multiple layers of interactions. To do this, we perform a linear stability analysis of the fully synchronized state of the interacting systems and exploit the spectral properties of the graph Laplacians of each layer. The final result is a system of coupled linear ordinary differential equations for the evolution of the displacements of the network from its synchronized state. Our setting does not require (nor assume) special conditions concerning the structure of each single layer, except that the network is undirected and that the local and interaction dynamics are described by continuous and differentiable functions. Because of this, the evolutionary differential equations are nonvariational. We validate our predictions in a network of chaotic Rössler oscillators with two layers of interactions featuring different topologies. We show that, even in this simple case, there is the possibility of inducing the overall stability of the complete synchronization manifold in regions of the phase diagram where each layer, taken individually, is unstable.


The model

From the structural point of view, we consider a network composed of N nodes, which interact through M different layers of connections, each layer generally having different links and representing a different kind of interactions among the units (see Fig. 1 for a schematic illustration of the case of M = 2 layers and N = 7 nodes). Note that in our setting, the nodes interacting in each layer are literally the same elements. Node i in layer 1 is precisely the same node as node i in layer 2, 3, or M. This contrasts with other works in which there is a one-to-one correspondence between nodes in different layers, but these represent potentially different states. The weights of the connections between nodes in layer α (α = 1,…, M) are given by the elements of the matrix W(α), which is, therefore, the adjacency matrix of a weighted graph. The sumEmbedded Imageof the weights of all the interactions of node i in layer α is the strength of the node in that layer.

Fig. 1 Schematic representation of a network with two layers of interaction.

The two layers (corresponding here to solid violet and dashed orange links) are made of links of different types for the same nodes, such as different means of transport between two cities or chemical and electric connections between neurons. Note that the layers are fully independent, in that they are described by two different Laplacians L(1) and L(2), so that the presence of a connection between two nodes in one layer does not affect their connection status in the other.

Regarding the dynamics, each node represents a d-dimensional dynamical system. Thus, the state of node i is described by a vector xi with d components. The local dynamics of the nodes is captured by a set of differential equations of the formEmbedded Imagewhere the dot indicates time derivative and F is an arbitrary C1-vector field. Similarly, the interaction in layer α is described by a continuous and differentiable vector field Hα (generally different from layer to layer), possibly weighted by a layer-dependent coupling constant σα. We assume that the interactions between node i and node j are diffusive, that is, for each layer in which they are connected, their coupling depends on the difference between Hα evaluated on xj and xi. Then, the dynamics of the whole system is described by the following set of equationsEmbedded Image(1)where L(α) is the graph Laplacian of layer α, whose elements areEmbedded Image(2)

Note that our treatment of this setting is valid for all possible choices of F and Hα, so long as they are C1, and for any particular undirected structure of the layers. This stands in contrast to other approaches to the study of the same equation set (1) proposed in previous works (and termed as dynamical hypernetworks), which, although based on ingenious techniques such as simultaneous block diagonalization, can be applied only to special cases [such as commuting Laplacians, unweighted and fully connected layers, and nondiffusive coupling (74)] or cannot guarantee to always provide a satisfactory solution (75).

Stability of complete synchronization in networks with multiple layers of interactions

We are interested in assessing the stability of synchronized states, which means determining whether a system eventually returns to the synchronized solution after a perturbation. For further details of the following derivations, we refer to Materials and Methods.

First, note that because Laplacians are zero-row-sum matrices, they all have a null eigenvalue, with corresponding eigenvector N−1/2 (1, 1,…, 1)T, where T indicates transposition. This means that the general system of equations (1) always admits an invariant solution S ≡ {xi(t) = s(t), ∀i = 1, 2,…, N}, which defines the complete synchronization manifold in ℝdN.

Because one does not need a very strong forcing to destroy synchronization in an unstable state, we aim at predicting the behavior of the system when the perturbation is small. First, we linearize Eq. 1 around the synchronized manifold S by obtaining the equations ruling the evolution of the local and global synchronization errors δxixis and δX ≡ δx1, δx2,…, δxN)TEmbedded Image(3)where Embedded Image is the N-dimensional identity matrix, ⨂ denotes the Kronecker product, and J is the Jacobian operator.

Second, we spectrally decompose δX in the equation above and project it onto the basis defined by the eigenvectors of one of the layers. The particular choice of layer is completely arbitrary because the eigenvectors of the Laplacians of each layer form M equivalent bases of ℝN. In the following, to fix the ideas, we operate this projection onto the eigenvectors of L(1). With some algebra, the system of equations 3 can then be expressed asEmbedded Image(4)for j = 2,…, N, where ηj is the vector coefficient of the eigendecomposition of δX, and Embedded Image is the rth eigenvalue of the Laplacian of layer α, sorted in a nondecreasing order; we have put Embedded Image, in which V(α) indicates the matrix of eigenvectors of the Laplacian of layer α. Note that to obtain this result, one must ensure that the Laplacian eigenvectors of each layer are orthonormal, a choice that is always possible because all the Laplacians are real symmetric matrices. Thus, the sums run from 2 rather than 1 because the first eigenvalue of the Laplacian, corresponding to r = 1, is always 0 for all layers, and the first eigenvector, to which all others are orthogonal, is common to all layers. Equation 4 is notable in that it includes previous results about systems with commuting Laplacians as a special case. If the Laplacians commute, they can be simultaneously diagonalized by a common basis of eigenvectors. Thus, in this case, V(α) = V(1)V for all α. In turn, this implies that Г(α) = Embedded Image for all α, and Eq. 4 becomesEmbedded Imagerecovering an M-parameter variational form as in (74).

Note that the stability of the synchronized state is completely specified by the maximum conditional Lyapunov exponent Λ, corresponding to the variation of the norm of Ω ≡ (η2,…, ηN). Because Ω will evolve on average as |Ω|(t)~exp(Λt), the fully synchronized state will be stable against small perturbations only if Λ < 0.

Case study: Networks of Rössler oscillators

To illustrate the predictive power of the framework described above, we apply it to a network of identical Rössler oscillators, with two layers of connections. Note that our method is fully general, and it can be applied to systems composed by any number of layers and containing oscillators of any dimensionality d. The particular choice of M = 2 and d = 3 for our example allows us to study a complex phenomenology while retaining ease of illustration. The dynamics of the Rössler oscillators is described by Embedded Image, where we have put xx1, yx2, and z ≡ x3. The parameters are fixed to the values a = 0.2, b = 0.2, and c = 9, which ensure that the local dynamics of each node is chaotic.

Considering each layer of connections individually, it is known that the choice of the function H allows (for an ensemble of networked Rössler oscillators) the selection of one of the three classes of stability (see Materials and Methods for more details), which are

(i) H(x) = (0, 0, z), for which synchronization is always unstable;

(ii) H(x) = (0, y, 0), for which synchronization is stable only for Embedded Image;

(iii) H(x) = (x, 0, 0) for which synchronization is stable only for Embedded Image.

Because of the double layer structure, one can now combine together different classes of stability in the two layers, studying how one affects the other and identifying new stability conditions arising from the different choices. In the following, we consider three combinations, namely

Case 1: layer 1 in class I and layer 2 in class II, that is, H1(x) = (0, 0, z) and H2(x) = (0, y, 0);

Case 2: layer 1 in class I and layer 2 in class III, that is, H1(x) = (0, 0, z) and H2(x) = (x, 0, 0);

Case 3: layer 1 in class II and layer 2 in class III, that is, H1(x) = (0, y, 0) and H2(x) = (x, 0, 0).

As for the choices of the Laplacians L(1,2), we consider three possible combinations: (i) both layers as Erdős-Rényi networks of equal mean degree (ER-ER); (ii) both layers as scale-free networks with power-law exponent 3 (SF-SF); and (iii) layer 1 as an Erdős-Rényi network and layer 2 as a scale-free network (ER-SF). In all cases, the graphs are generated using the algorithm of Gómez-Gardeñes and Moreno (76), which allows a continuous interpolation between scale-free and Erdős-Rényi structures (see Materials and Methods for details). Therefore, in the following, we will consider nine possible scenarios, that is, the three combinations of stability classes for each of the three combinations of layer structures.

Case 1.. Rewriting the system of equations (4) explicitly for each component of ηj, we obtain hereEmbedded Image(5)Embedded Image(6)Embedded Image(7)from which the maximum Lyapunov exponent can be numerically calculated. As shown in the top panel of Fig. 2, we observe that for ER-ER topologies, the first layer is dominated by the second because the stability region of the whole system appears to be almost independent of σ1, disregarding a slight increase of the critical value of σ2 as σ1 increases. This demonstrates the ability of class II systems to control the instabilities inherent to systems in class I. This result appears to be robust with respect to the choice of underlying structures because qualitatively similar results are obtained for SF-SF, ER-SF, and SF-ER topologies (fig. S1).

Fig. 2 Maximum Lyapunov exponent for ER-ER topologies in case 1 (top panel) and case 2 (bottom panel).

The darker blue lines mark the points in the (σ1, σ2) space where Λ vanishes, whereas the striped lines indicate the critical values of σ2 if layer 2 is considered in isolation (or, equivalently, if σ1 = 0).

Case 2.. For case 2, the system of equations (4) readsEmbedded Image(8)Embedded Image(9)Embedded Image(10)

As shown in the bottom panel of Fig. 2, also in this case, the second layer strongly dominates the whole system because the overall stability window is almost independent from the value of σ1. This result, together with that obtained for case 1, suggests that class I systems, although intrinsically preventing synchronization, are easily controllable by both class II and class III systems, even though, analogous to case 1, we observe a slight widening of the stability window for increasing values of σ1. Again, the results are almost independent from the choice of the underlying topologies (fig. S2).

Case 3.. Finally, for case 3, equations (4) becomeEmbedded Image(11)Embedded Image(12)Embedded Image(13)

Here, the system reveals its most striking features. In particular, for ER-ER topologies (see top panel of Fig. 3), we observe six different regions, identified in the figure by Roman numerals. In region I, synchronization is stable in both layers taken individually (or, equivalently, for either σ1 = 0 or σ2 = 0), and the full bilayered network is also stable. Regions II, III, and IV correspond to scenarios qualitatively similar to the ones seen previously, that is, where stability properties of one layer dominate over those of the other. Finally, regions V and VI are the most important because within them, one finds effects that are genuinely arising from the multilayered nature of the interactions. In these regions, both layers are individually unstable, and synchronization would not be observed at all for either σ1 = 0 or σ2 = 0. However, the emergence of a collective synchronous motion is remarkably obtained with a suitable tuning of the parameters. In these regions, it is therefore the simultaneous action of the two layers that induces stability.

Fig. 3 Maximum Lyapunov exponent in case 3 for ER-ER and SF-SF topologies (top and bottom panel, respectively).

The darker blue lines mark the points in the (σ1, σ2) plane where the maximum Lyapunov exponent is 0, whereas the striped lines indicate the stability limits for the σ1 = 0 and σ2 = 0. The points marked in the top panel indicate the choices of coupling strengths used for the numerical validation of the model. Note that for SF networks in class III, the stability window disappears.

Together, the results we obtained for the three cases indicate that a multilayer interaction topology enhances the stability of the synchronized state, even allowing for the possibility of stabilizing systems that are unstable when considered isolated.

Numerical validation

We validate the stability predictions derived from equations (4) by simulating the full nonlinear system of equations (1) for an ER-ER topology in case 3, with three different choices of coupling constants σ1 and σ2. The three specific sets of coupling values (shown in the top panel of Fig. 3) correspond to situations in which either one or both layers are unstable when isolated but yield a stable synchronized state when coupled. More specifically, we have chosen (σ1 = 0.04, σ2 = 0.3), (σ1 = 0.15, σ2 = 0.5), and (σ1 = 0.04, σ2 = 0.5) corresponding to regions II, IV, and VI, respectively.

For all three cases, we run the simulations initially with the presence of only the unstable layer, by setting either σ1 = 0 or σ2 = 0 depending on the set of couplings considered. Note that for the third set of couplings (region VI), either layer can be the initially active one because both are unstable when isolated. After 100 integration steps, we then activate the other layer by setting its interaction strength to the (nonzero) value corresponding to the region for which we predicted a stable synchronized state. As the systems evolve, we monitor the evolution of the norm |Ω|(t) to evaluate the deviation from the synchronized solution with time.

The results in Fig. 4 show that when only the unstable interaction layer is active, |Ω|(t) never vanishes. However, as soon as the other layer is switched on, the norm of Ω undergoes a sudden change of behavior, starting an exponential decay toward 0. This confirms the prediction that the unstable behavior induced by each layer is compensated by the mutual presence of two interaction layers.

Fig. 4 Numerical validation of the stability analysis.

The error of synchronization increases as long as the only active layer is the one predicted to be unstable. When the other layer is switched on, at time 100, the error of synchronization decays exponentially toward 0, as predicted by the model. With respect to Fig. 3, the top left panel corresponds to region II, where layer 1 is unstable and layer 2 is stable, and the interaction strengths used were σ1 = 0.04 and σ2 = 0.3. The bottom left panel corresponds to region IV, where layer 1 is stable and layer 2 is unstable, and the interaction strengths were σ1 = 0.15 and σ2 = 0.5. The top right and bottom right panels correspond to region VI, where both layers are unstable. The layer active from the beginning was layer 1 for the top right panel and layer 2 for the bottom right panel. In both cases, the interaction strengths were σ1 = 0.04 and σ2 = 0.5.

Qualitatively similar scenarios are observed in case 3 for SF-SF topologies, as well as for ER-SF and SF-ER structures (fig. S3). Again, they confirm the correctness of the predictions, showing that in region I, layer 1 dominates over layer 2, and that in region II, overall stability can be induced even when both layers are unstable in isolation.

To provide an even stronger demonstration of the predictive power of our method, we simulate the full system for the ER-ER topology in case 3, fixing the value of σ2 to 1 and varying the value of σ1 from 0 to 0.2. Starting from an initial perturbed synchronized state, after a transient of 100 time units, we measure the average of |Ω| over the next 20 integration steps. The results in Fig. 5 show agreement between the simulations and the theoretical prediction (see Fig. 4). For values of σ1 less than a critical value of approximately 0.04, the system never synchronizes. Conversely, when σ1 crosses the critical value, the system can reach a synchronized state. Repeating the simulation with σ2 = 0, one recovers the monoplex case. Also in this instance, we find agreement between theoretical prediction and simulation, with a critical coupling value of approximately 0.08.

Fig. 5 Identification of the critical points.

For a system with ER-ER topology in case 3 and fixed σ2 = 1, the synchronization error never vanishes if σ1 < σc ≈ 0.04. Conversely, as soon as σ1 > σc, the system is again able to synchronize (green line). One recovers the monolayer case by imposing σ2 = 0, for which similar results are found, with a critical coupling strength of approximately 0.08 (red line). Both results are in perfect agreement with the theoretical predictions (see Fig. 4).


The results shown above clearly illustrate the rich dynamical phenomenology that emerges when the multilayer structure of real networked systems is taken into account. In an explicit example, we have observed that synchronization stability can be induced in unstable networked layers by coupling them with stable ones. In addition, we have shown that stability can be achieved even when all the layers of a complex system are unstable, if considered in isolation. This latter result constitutes a clear instance of an effect that is intrinsic to the true multilayer nature of the interactions among the dynamical units. Similarly, we expect that the opposite could also be observed, namely, that the synchronizability of a system decreases or even disappears when two individually synchronizable layers are combined.

On more general grounds, the theory developed here allows one to assess the stability of the synchronized state of coupled nonlinear dynamical systems with multilayer interactions in a fully general setting. The system can have any arbitrary number of layers, and perhaps more importantly, the network structures of each layer can be fully independent because we do not exploit any special structural or dynamical property to develop our theory. This way, our approach generalizes the celebrated MSF (30) to multilayer structures, retaining the general applicability of the original method. The complexity in the extra layers is reflected in the fact that the formalism yields a set of coupled linear differential equations (Eq. 4), rather than a single parametric variational equation, which is recovered in the case of commuting Laplacians. This system of equations describes the evolution of a set of d-dimensional vectors that encode the displacement of each dynamical system from the synchronized state. The solution of the system gives a necessary condition for stability: if the norm of these vectors vanishes in time, then the system gets progressively closer to synchronization, which is therefore stable; if, instead, the length of the vectors always remains greater than 0, then the synchronized state is unstable.

The generality of the method presented, which is applicable to any undirected structure, and its straightforward implementation for any choice of C1 dynamical setup pave the way for the exploration of synchronization properties on multilayer networks of arbitrary size and structure. Thus, we are confident that our work can be used in the design of optimal multilayered synchronizable systems, a problem that has attracted much attention in monolayer complex networks (7780). The straightforward nature of our formalism makes it suitable for efficient use together with successful techniques, such as the rewiring of links or the search for an optimal distribution of link weights, in the context of multilayer networks. In turn, these techniques may help in addressing the already mentioned question of the suppression of synchronization due to the interaction between layers, unveiling possible combinations of stable layers that, when interacting, suppress the dynamical coherence that they show in isolation. Also, we believe that the reliability of our method will provide aid to the highly current field of multiplex network controllability (26, 8184), enabling researchers to engineer control layers to drive the system dynamics toward a desired state.

In addition, several extensions of our work toward more general systems are possible. A particularly relevant one is the study of multilayer networks of heterogeneous oscillators, which have a rich phenomenology and whose synchronizability has been shown to depend on all the Laplacian eigenvalues (85) in a way similar to the results presented here. Relaxing the requirement of an undirected structure, our approach can also be used to study directed networks. The graph Laplacians in this case are not necessarily diagonalizable, but a considerable amount of information can be still extracted from them using singular value decomposition. For example, it is already known that directed networks can be rewired to obtain an optimal distribution of in-degrees for synchronization (86). Further areas that we intend to explore in future work are those of almost identical oscillators and almost identical layers, which can be approached using perturbative methods and constitute more research directions with even wider applicability.

Finally, as our method allows one to study the rich synchronization phenomenology of general multilayer networks, we believe that it will find application in technological, biological, and social systems where synchronization processes and multilayered interactions are at work. Some examples are coupled power grid and communication systems, some brain neuropathologies such as epilepsy, and the onset of coordinated social behavior when multiple interaction channels coexist. As mentioned above, these applications will demand further advances to include specific features such as the nonhomogeneity of interacting units or the possibility of directional interactions.


Linearization around the synchronized solution

To linearize the system in Eq. 4 around the synchronization manifold, use the fact that for any C1-vector field f we can writeEmbedded Image

Using this relation, we can expand F and H around s in the system of equations 4 to obtainEmbedded Image(14)

Now, we use the Kronecker matrix product to decompose the equation above into self-mixing and interaction terms and introduce the vector δX, to get the final system of equations 3. The system 3 can be rewritten by projecting δX onto the Laplacian eigenvectors of a layer. The choice of layer to carry out this projection is entirely arbitrary because the Laplacian eigenvectors are always a basis of ℝN. Without the loss of generality, we choose layer 1 and ensure that the eigenvectors are orthonormal. We then define Embedded Imaged to be the d-dimensional identity matrix and multiply Eq. 3 on the left by Embedded ImageEmbedded Image

Now, we use the relationEmbedded Image(15)to obtainEmbedded Imagewhere D(α) is the diagonal matrix of the eigenvalues of layer α, and we have split the sum into the first term and the remaining M − 1 terms. Left-multiply the first occurrence of Embedded Image in the right-hand side by Embedded Image, and right-multiply F and H1 by Embedded Image. Then, using again Eq. 15, it becomesEmbedded Image

Factor out Embedded Image to getEmbedded Image

The relationEmbedded Imageimplies that Embedded Image is the mN-dimensional identity matrix. Then, we left-multiply the last δX by this expression, obtainingEmbedded Image

Now, we define the vector-of-vectorsEmbedded Image

Each component of η is the projection of the global synchronization error vector δX onto the space spanned by the corresponding Laplacian eigenvector of layer 1. The first eigenvector, which defines the synchronization manifold, is common to all layers, and all other eigenvectors are orthogonal to it. Thus, the norm of the projection of η over the space spanned by the last N – 1 eigenvectors is a measure of the synchronization error in the directions transverse to the synchronization manifold. Because of how η is built, this projection is just the vector Ω, consisting of the last N – 1 components of η. With this definition of η, we left-multiply L(α) by the identity expressed as Embedded Image to obtainEmbedded Image

In this vector equation, the first part is purely variational because it consists of a block-diagonal matrix that multiplies the vector-of-vectors η. However, the second part mixes different components of η. This can be seen more easily expressing the vector equation as a system of equations, one for each component j of η.

To write this system, it is convenient to first define Embedded Image, then consider the nonvariational part. Its contribution to the jth component of Embedded Image is given by the product of the jth row of blocks of the block-matrix multiplied by η. In turn, each element of this row of blocks consists of the corresponding element of the jth row of Embedded Image multiplied by JHα (s)Embedded Image

Summing over all the components ηk yieldsEmbedded Imagewhich is Eq. 4. Note that the sums over r and k start from 2 because the first eigenvalue is always 0, and the orthonormality of the eigenvectors guarantees that all the elements of the first column of Г(α), except the first, are 0. Each matrix Г(α) effectively captures the alignment of the Laplacian eigenvectors of layer α with those of layer 1. If the eigenvectors for layer α are identical to those of layer 1, as it happens when the two Laplacians commute, then Г(α) = Embedded Image. One can generalize the definition of Г(α) to consider any two layers, introducing the matrices Embedded Image that can even be used to define a measure ℓD of “dynamical distance” between two layers α and βEmbedded Image

MSF and stability classes

A particular case of the treatment we considered above happens when M = 1. In this case, the second term on the right-hand side of Eq. 4 disappears, and the system takes the variational form Embedded Image, where KiJF(s) – σλiJH(s) is an evolution kernel evaluated on the synchronization manifold. Because λ1 = 0, this equation separates the contribution parallel to the manifold, which reduces to Embedded Image, from the other N – 1, which describe perturbations in the directions transverse to the manifold and which have to be damped for the synchronized state to be stable. Because the Jacobians of F and H are evaluated on the synchronized state, the variational equations differ only in the eigenvalues λi. Thus, one can extract from each of them a set of d-conditional Lyapunov exponents, evaluated along the eigenmodes associated to λi. Putting ν ≡ σλi, the parametrical behavior of the largest of these exponents Λ(ν) defines the so-called MSF (30). If the network is undirected, then the spectrum of the Laplacian is real, and the MSF is a real function of ν. Crucially, for all possible choices of F and H, the MSF of a network falls into one of three possible behavior classes, defined as follows (6)

class I: Λ(ν) never intercepts the x axis;

class II: Λ(ν) intercepts the x axis at a single point at some νc ≥ 0;

class III: Λ(ν) is a convex function with negative values within some window νc1 < ν < νc2; in general, νc1 ≥ 0, with the equality holding when F supports a periodic motion.

The elegance of the MSF formalism manifests itself at its finest for systems in class III, for which synchronization is stable only if σλ2 > νc1 and σλN < νc2 hold simultaneously. This condition implies Embedded Image. Because Embedded Image is entirely determined by the network topology and Embedded Image depends only on the dynamical functions F and H, one has a simple stability criterion in which structure and dynamics are decoupled.

Network generation

To generate the networks for our simulations, we used the algorithm described by Gómez-Gardeñes and Moreno (76), which creates a one-parameter family of complex networks with a tunable degree of heterogeneity. The algorithm works as follows: start from a fully connected network with m0 nodes and a set χ containing Nm0 isolated nodes. At each time step, select a new node from χ and link it to other m nodes, selected among all other nodes. The choice of the target nodes happens uniformly at random with probability α and following a preferential attachment rule with probability 1 – α. Repeating these steps Nm0 times, one obtains networks with the same number of nodes and links, whose structure interpolates between ER, for α = 1, and SF, for α = 0.

Numerical calculations

To compute the maximum Lyapunov exponent for a given pair of coupling strengths σ1 and σ2, we first integrate a single Rössler oscillator from an initial state (0, 0, 0) for a transient time ttrans, sufficient to reach the chaotic attractor. The integration is carried out using a fourth-order Runge-Kutta integrator with a time step of 5 × 10–3, for which we choose a transient time ttrans = 300. Then, we integrate the systems for the perturbations (Eqs. 5 to 7, 8 to 10, and 11 to 13) using Euler’s method, again with a same time step of 5 × 10–3. The initial conditions are such that all components of all the ηj are Embedded Image, making Ω a unit vector. At the same time, we continue the integration of the single Rössler unit to provide for s1 and s3, which appear in the perturbation equations. This process is repeated for 500 time windows, each with the duration of 1 unit (200 steps). After each window n, we compute the norm of the overall perturbation |Ω|(n) and rescale the components of the ηj so that at the start of the next time window, the norm of Ω is again set to 1. Finally, when the integration is completed, we estimate the maximum Lyapunov exponent asEmbedded Image


Supplementary material for this article is available at

fig. S1. Maximum Lyapunov exponent Λ for systems falling into case 1 (layer 1 in stability class I and layer 2 in stability class II) for SF-SF, ER-SF, and SF-ER topologies (left, center, and right panels, respectively).

fig. S2. Maximum Lyapunov exponent Λ for systems falling into case 2 (layer 1 in stability class I and layer 2 in stability class III) for SF-SF, ER-SF, and SF-ER topologies (left, center, and right panels, respectively).

fig. S3. Maximum Lyapunov exponent Λ for systems falling into case 3 (layer 1 in stability class II and layer 2 in stability class III) for ER-SF and SF-ER topologies (left and right panels, respectively).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank A. Arenas and J. Buldú for interesting and fruitful discussions. Funding: The work of J.G.-G. was supported by the Spanish MINECO through grants FIS2012-38266-C02-01 and FIS2011-25167 and by the European AQ38 Union through FET Proactive Project MULTIPLEX (Multilevel Complex Networks and Systems) contract no. 317532. Author contributions: C.I.d.G. and S.B. developed the theory. S.B. designed the simulations. J.G.-G. implemented and carried out the simulations and analyzed the results. All authors wrote the manuscript. Competing interests: The authors declare they have no competing interests. Data and material availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article