## Abstract

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

## INTRODUCTION

Network theory (*1*–*9*) 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 (*10*–*12*). 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 (*13*–*19*).

Because networks often support dynamical processes, the interplay between the structure and the unfolding of collective phenomena has been the subject of numerous studies (*20*–*22*). 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 (*27*–*29*). 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 (*30*–*42*).

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 (*43*–*46*). 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 (*53*–*57*), epidemics spreading (*58*–*61*), controllability (*62*), evolutionary games (*63*–*66*), 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 (*70*–*73*).

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.

## RESULTS

### 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 sumof the weights of all the interactions of node *i* in layer α is the strength of the node in that layer.

Regarding the dynamics, each node represents a *d*-dimensional dynamical system. Thus, the state of node *i* is described by a vector **x**_{i} with *d* components. The local dynamics of the nodes is captured by a set of differential equations of the formwhere the dot indicates time derivative and **F** is an arbitrary *C*^{1}-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 **x**_{j} and **x**_{i}. Then, the dynamics of the whole system is described by the following set of equations(1)where **L**^{(α)} is the graph Laplacian of layer α, whose elements are(2)

Note that our treatment of this setting is valid for all possible choices of **F** and **H**_{α}, so long as they are *C*^{1}, 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** ≡ {**x**_{i}(*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 δ**x**_{i} ≡ **x**_{i} − **s** and δ**X** ≡ δ**x**_{1}, δ**x**_{2},…, δ**x**_{N})^{T}(3)where 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 as(4)for *j* = 2,…, *N*, where **η**_{j} is the vector coefficient of the eigendecomposition of δ**X**, and is the *r*th eigenvalue of the Laplacian of layer α, sorted in a nondecreasing order; we have put , 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 Г^{(α)} = for all α, and Eq. 4 becomesrecovering 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 , where we have put *x* ≡ *x*_{1}, *y* ≡ *x*_{2}, and z ≡ *x*_{3}. 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 ;

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

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, **H**_{1}(**x**) = (0, 0, *z*) and **H**_{2}(**x**) = (0, *y*, 0);

**Case 2:** layer 1 in class I and layer 2 in class III, that is, **H**_{1}(**x**) = (0, 0, *z*) and **H**_{2}(**x**) = (*x*, 0, 0);

**Case 3:** layer 1 in class II and layer 2 in class III, that is, **H**_{1}(**x**) = (0, *y*, 0) and **H**_{2}(**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 here(5)(6)(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).

Case 2.. For case 2, the system of equations (4) reads(8)(9)(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) become(11)(12)(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.

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.

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.

## DISCUSSION

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 *C*^{1} 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 (*77*–*80*). 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*, *81*–*84*), 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.

## MATERIALS AND METHODS

### Linearization around the synchronized solution

To linearize the system in Eq. 4 around the synchronization manifold, use the fact that for any *C*^{1}-vector field **f** we can write

Using this relation, we can expand **F** and **H** around **s** in the system of equations 4 to obtain(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 _{d} to be the *d*-dimensional identity matrix and multiply Eq. 3 on the left by

Now, we use the relation(15)to obtainwhere **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 in the right-hand side by , and right-multiply **F** and **H**_{1} by . Then, using again Eq. 15, it becomes

Factor out to get

The relationimplies that is the *mN*-dimensional identity matrix. Then, we left-multiply the last δ**X** by this expression, obtaining

Now, we define the vector-of-vectors

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 to obtain

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 , then consider the nonvariational part. Its contribution to the *j*th component of is given by the product of the *j*th 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 *j*th row of multiplied by *J***H**_{α} (**s**)

Summing over all the components **η**_{k} yieldswhich 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 Г^{(α)} = . One can generalize the definition of Г^{(α)} to consider any two layers, introducing the matrices that can even be used to define a measure ℓ_{D} of “dynamical distance” between two layers α and β

### 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 , where **K**_{i} ≡ *J***F**(**s**) – σλ_{i}*J***H**(**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 , 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 . Because is entirely determined by the network topology and 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 *m*_{0} nodes and a set χ containing *N* – *m*_{0} 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 *N* – *m*_{0} 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 *t*_{trans}, 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 *t*_{trans} = 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 , making **Ω** a unit vector. At the same time, we continue the integration of the single Rössler unit to provide for *s*_{1} and *s*_{3}, 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 as

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/11/e1601679/DC1

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2016, The Authors