Mean-field nature of synchronization stability in networks with multiple interaction layers

The interactions between the components of many real-world systems are best modelled by networks with multiple layers. Different theories have been proposed to explain how multilayered connections affect the linear stability of synchronization in dynamical systems. However, the resulting equations are computationally expensive, and therefore difficult, if not impossible, to solve for large systems. To bridge this gap, we develop a mean-field theory of synchronization for networks with multiple interaction layers. By assuming quasi-identical layers, we obtain accurate assessments of synchronization stability that are comparable with the exact results. In fact, the accuracy of our theory remains high even for networks with very dissimilar layers, thus posing a general question about the mean-field nature of synchronization stability in multilayer networks. Moreover, the computational complexity of our approach is only quadratic in the number of nodes, thereby allowing the study of systems whose investigation was thus far precluded. Multilayer networks can achieve synchronization, both for homogeneous and heterogeneous layers, whose dynamics is described by a system of equations often computationally complex and expensive. Here, the authors propose a mean-field approach for estimating the stability of the synchronized state of multilayer networks and show this applies to both homogeneous and heterogeneous layers, lowering computational complexity.

T he use of a network structure, consisting of a list of pairwise connections called edges between discrete elements called nodes or vertices, has long been a powerful abstraction to model and investigate the behavior of complex systems [1][2][3] . One area in which the network paradigm has proved especially useful is the study of distributed coupled dynamical systems 4 , where the nodes correspond to dynamical systems that interact across the edges. Collective and organized dynamics is widely studied in such networks, with a particular regard to the phenomenon of synchronization, which holds fundamental importance in numerous natural instances [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] .
Within this topic, considerable attention has been paid to the stability of synchronized states. A powerful method to estimate it is the so-called Master Stability Function (MSF) 6 , which allows one to compute the value of the largest Lyapunov exponent Λ max of the system upon perturbation from the synchronized state. Then, one can determine whether synchronization is stable simply from its sign: the perturbed system will asymptotically resynchronize only if Λ max is negative. The MSF has maintained its status as method of choice even when applied to models that extend the network approach to the multilayer case [21][22][23] and to networks with higher-order interactions 24 , which better capture the many levels of complexity of real-world systems.
In the dynamical multilayer networks we consider here, the nodes are allowed to interact over multiple layers, each representing a different type of interaction. Thus, they can be considered a special case of the traditional multiplex networks, in which each node is replicated in every layer, and interlayer links determine how the state in each layer affects those of the neighbouring ones. Networks with multiple interaction layers are often encountered in natural and engineered complex systems of wide relevance. The paradigmatic example is that of the nervous system, in which the same neurons are connected to each other via two different types of synapses 25 . Another instance is provided by the transport infrastructure of large cities where the same station is served by multiple means of mass transport. In both cases, the value of a variable of interest, such as the potential accumulated in a neuron or the state of congestion of a station, affects equally all the processes involving the node 26 .
The increase in complexity of the system with respect to singlelayer networks is reflected by the fact that the use of the MSF no longer results in a single equation, but rather produces a set of coupled linear differential equations 27 . Apart from their generality, these master stability equations allowed researchers to identify synchronized states whose stability properties are inherently due to the multiplex architecture, such as stable coherent dynamics in networked layers that are unstable when studied in isolation. Also, the versatility of the MSF made it a method of choice to study this type of networks under different conditions, such as time-varying structures 28 , as well as with the help of other approaches, such as fast-switching techniques 29 . Unfortunately, despite the promising preliminary results, the master stability equations also have a drawback, as the computational complexity of solving the system makes them unwieldy for networks larger than a couple of hundred nodes. This difficulty is due to the mathematical structure that results from the MSF approach, and therefore it does not depend on the specific network model studied or on the particular formalism used upon it. In fact, even when the MSF takes a fairly workable form, such as in traditional multiplex networks with dynamics formalized via supra-Laplacians 23 , the explicit numerical calculations reach a high computational complexity if the supra-Laplacians do not commute.
In this article, we introduce a mean-field theory of synchronization for networks with multiple interaction layers. We derive our theory in the assumption that the layers are quasi-identical, but show that its range of applicability encompasses the case of very dissimilar layers. In particular, the estimate of the stability of synchronization obtained with our theory virtually never changes with respect to the exact result, suggesting that global synchronization stability in multilayer networks is inherently a mean-field phenomenon. In addition, the numerical complexity of our method is lower than that of the exact formulation, making it applicable to large systems whose study would otherwise be prevented by computational costs.

Results
The model. The derivation of a MSF on a networked system is effectively a decomposition of the dynamics into eigenmodes. In a network with N nodes and M interaction layers, the connection weights of each layer α are the elements of the weighted adjacency matrix W (α) . Let x i be an m-dimensional vector describing the state of node i, and F : R m ! R m and H α : R m ! R m be continuous and differentiable vector fields describing the local dynamics and the interactions in layer α, respectively. Assuming diffusive coupling between the nodes, the global dynamics is determined by the system In Eq. (1), σ α is the interaction strength within layer α, and L α ð Þ is the graph Laplacian of layer α, whose diagonal elements are L ðαÞ i;j , and whose off-diagonal elements are L ðαÞ i;j ¼ ÀW ðαÞ i;j . Following the approach by del Genio et al. 27 , one can linearize Eq. (1), obtaining expressions for the evolution of the global synchronization error vector, which can be projected onto the eigenvectors of the Laplacian of one of the layers. Choosing layer α = 1 (without loss of generality) results in the system here, η is a vector of vectors whose jth component is the projection of the global synchronization error vector onto the space spanned by the jth Laplacian eigenvector of layer 1; J is the Jacobian operator; s is the state vector corresponding to the synchronized state; λ ðαÞ j is the jth Laplacian eigenvalue of layer α; and Γ (α) is the spectral overlap matrix between layer 1 and layer α,

is the matrix of Laplacian eigenvectors of layer α and T indicates transposition.
Mean-field theory. To develop our theory, we start from the assumption that the interaction layers are quasi-identical. Then, due to their adjacency matrices and their Laplacians being very similar, one can expect their Laplacian eigenvectors to be equal up to some small perturbation. Our goal is to spread the effect of this perturbation in an equal way over all the directions transverse to the synchronization manifold, which is identified by the Laplacian eigenvector corresponding to the null eigenvalue. To do so, compute the dynamical distance between layer α and layer 1 27 : Note that, while here we consider dynamical distances between any layer and the reference layer 1, in principle this quantity is an indicator of the similarity between the dynamics of any two layers. In fact, its general definition is the sum of the squares of the off-diagonal terms in the spectral overlap matrix between the two layers considered, so that if the dynamics of the layers are identical (i.e., their Laplacians commute), their dynamical distance vanishes. In the equation above, the matrices Ξ (α, β) are defined as Ξ ðα; Also, the sums in Eq. (3) start from 2 because the first Laplacian eigenvector is always N À 1 = 2 1; ; 1 ð Þ T . Therefore, Γ ðαÞ 1;1 ¼ 1 and Γ ðαÞ 1;k ¼ Γ ðαÞ k;1 ¼ 0 for all k > 1. Now note that each Γ (α) is an orthogonal matrix, since it is the product of two orthogonal matrices. Moreover, from the definition, it follows that V ðαÞ ¼ V ð1Þ Γ ðαÞ T . In other words, Γ (α) is the transformation matrix from the Laplacian eigenvectors of the first layer to those of layer α. Then, our aim is replacing Γ (α) with a Γ (α),MF whose action is to change each eigenvector of the first layer in the same fashion.
As Γ (α) is a rotation, a natural choice is to make Γ (α),MF rotate the eigenvectors of the first layer by the same angle in every direction. More precisely, we consider all possible 2-dimensional subspaces of R NÀ1 determined by choosing any two Laplacian eigenvectors of layer 1 except the first, and then construct a matrix that rotates all the Laplacian eigenvectors of layer 1 except the first by the same amount in each of these subspaces. Note that, in principle, this problem is underspecified, as rotation matrices do not commute in 3 or more dimensions, and the specific form of Γ (α),MF depends on the order in which the rotations in the individual subspaces are performed. However, since the layers are quasi-identical, the rotation angle needed is very small, as we will justify quantitatively later on. Thus, the rotations that compose Γ (α),MF are infinitesimal. In turn, this means that each one of them, and indeed Γ (α),MF itself, can be written as the sum of the identity matrix and an element of the Lie algebra of O N À 1 ð Þ, or, more precisely, of SO N À 1 ð Þ, since the rotations are proper. Thus, for all 1 < r < s, the elements of the matrix R (r,s) that operates the rotation in the subspace spanned by the rth and sth eigenvectors are where ε (α) is the rotation angle, which depends on the layer α. Importantly, all the R (r,s) commute, which removes the problem of considering the order of the constitutive rotations of Γ (α),MF . Taking into account all subspaces, and neglecting terms of order higher than ε (α) , yields the following form for the matrix: . . .
Notice that Γ (α),MF is an orthogonal matrix, as it should be, but only to first order, as expected from the approximations used. Also note that the final results one obtains using this formulation are the same that would be found using the exact expressions on layers with a different structure. This mean-field-equivalent structures can be found by first computing the mean-fieldequivalent Laplacian eigenvectors using the same equation as above, namely V ðαÞ;MF ¼ V ð1Þ Γ ðαÞ;MF À Á T , and then using them to recover the mean-field-equivalent Laplacian itself. To find the value of ε (α) to be used in Eq. (4), one can notice that ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi . But from Eq. (3) and the fact that Ξ (α, 1) = Γ (α) it follows that this quantity has to be equal to ffiffiffiffiffiffiffi ' α;1 p . Thus, Note that Eq. (5) effectively provides a mean of the dynamical distance over all the directions transverse to the synchronization manifold, thereby justifying our definition of the method as a mean-field theory. Finally, to find the mean-field form of the MSF, one can notice that Eq. (2) contains the product Γ ðαÞ r;k Γ ðαÞ r;j within its innermost sum. Then, replacing the matrices Γ (α) with their mean-field version, one obtains if j < i and j > 1 0 otherwise: Since, in Eq. (2), j, r and k are always greater than 1, one can write if r ¼ k and j > r; or if k > r and r ¼ j ε ðαÞ if r ¼ k and j < r; or if k < r and r ¼ j ε ðαÞ 2 % 0 if k > r and j > r; or if k < r and j < r Àε ðαÞ 2 % 0 if k > r and j < r; or if k < r and j > r: Using the expression above, the second term on the right-hand side of Eq. (2) becomes where we have split the sum over k into k = j, k < j and k > j. The expression (6) greatly simplifies the calculations with respect to the original formulation. In fact, the computational complexity of calculating each component of η according to Eq. (2) is O N 2 M À Á , whereas using Eq. (6), this reduces to O NM ð Þ. Also note that in Eq. (6), the first term inside the sum over α corresponds to the case of commuting Laplacians. Thus, while the transverse modes are still not completely decoupled, which would only be possible in the commuting case, our theory effectively consists in a first-order correction, obtained by a mean-field perturbative approximation of the dynamics. This becomes even more evident when rewriting the whole Eq. (2) as for which we assume that the Laplacian eigenvalues are sorted in a non-decreasing way.
The derivation of Eq. (7), which constitutes the mean-field approximation of Eq. (2), not only offers a decrease in the complexity of evaluating the linear stability of this type of multiplexes, but also, more importantly, paves the way to compute the stability diagram of systems whose size is too large to permit the use of Eq. (2). For example, a scaling test reveals that simulation of a single duplex random network with 10000 nodes would take approximately 103 years on a workstation with a 16-core Intel® Xeon® Gold 6130 CPU if using Eq. (2). This time would be reduced to just above 5 days if using Eq. (6).
Numerical validation. To demonstrate the validity and applicability of our theory, we carried out extensive numerical simulations on 2-layer random networks of chaotic Rössler oscillators. Note that the dynamical distance between layers 2 and 1, ℓ 2,1 , is bound between 0 and N − 1. Thus, to allow for a direct comparison between systems of different size, we henceforth plot our results as a function of the normalized dynamical distancẽ ' 2;1 ¼ ' 2;1 = NÀ1 ð Þ . In Fig. 1 we show the average error in the determination of the largest Lyapunov exponent Λ max as a function of' . Our results indicate that the accuracy of the theory is generally high, but with a marked dependence on the region of the layer-stability phase diagram. More specifically, the accuracy increases with the ratio σ 1 = σ 2 , which is on average different between the regions of the phase diagram (schematically illustrated in the inset of Fig. 1). This is consistent with our theory being a perturbative correction to the commuting case, since the relative contribution of such a correction is indeed inversely proportional to the ratio of the interaction strengths within the layers. At the same time, and as we stressed above, the choice of the reference layer 1 is entirely arbitrary. As such, one can always choose the layer with the largest interaction strength to be layer 1, which implies that the actual worst case scenario is the one in which all layers have the same interaction strength.
At the light of these considerations, we performed additional simulations, now imposing σ 1 = σ 2 , to evaluate the performance of our theory for increasing system sizes. The results, illustrated in Fig. 2, show that the theory provides accurate results for a wide range of normalized dynamical distances between layers.
To further explore the scope of applicability of our method, we carried out extra simulations on networks consisting of three random layers, as well as on preferential-attachment networks 30 . The results are presented in three figures ( Supplementary  Figs. 1-3) of the Supplementary Information. The case of three-layer networks shows a general improvement with respect to the two-layer case, which is particularly evident for larger systems sizes (see Supplementary Fig. 1). Also, the errors in the assessment of synchronization stability completely vanish, even for smaller networks. The situation is slightly different for preferential-attachment networks, where the errors in assessing stability are not identically null, but they remain fairly negligible, reaching a maximum of 0.8% at normalized dynamical distances of 0.3 or more (see Supplementary Fig. 2). The average relative errors in estimating the maximum Lyapunov exponent remain instead roughly in line with those for random networks, but without a clear trend with network size (see Supplementary  Fig. 3).

Discussion
The simulation results illustrated above show that the error in the estimation of Λ max is never larger than 30%, even for normalized dynamical distances as high as' ¼ 0:6 . Importantly, wheñ ' < 0:1, the accuracy of the estimate increases with the system size. This is most likely a direct consequence of Eq. (5), which shows that, for the same value of the normalized dynamical distance, ε (α) decreases proportionally to ffiffiffiffi N p . The applicability of our theory increases even more when considering the typical use of the MSF, namely the assessment of the stability of synchronization. In this case, our method virtually never fails in correctly identifying stable and unstable synchronized states. In our simulations, we only found occasional errors for at most 0.2% of the cases, and exclusively in networks of only 10 nodes. Given the very small size of such systems, we believe these errors to have been caused by finite-size effects, indicating that our theory always allows a correct determination of synchronization stability, and it holds true far beyond infinitesimal dynamical distances.
Note that, as this method is an approximation, we do expect that it will present limitations. At the current stage, we can only speculate that, in a manner similar to other mean-field theories, it may stop providing correct results in the vicinity of critical transitions, or for networks with pathological structures. Such Fig. 1 The accuracy of our theory increases with the ratio of the interaction strengths of the first layer σ 1 to that of the second layer σ 2 . The relative error in the estimate of the largest Lyapunov exponent Λ max of the perturbed system decreases from a maximum of approximately 40% in region 6 ( σ 1 = σ 2 % 0:0917) to a minimum of approximately 0.7% in region 3 ( σ 1 = σ 2 % 3:5). Within each region, the error increases sublinearly with the normalized dynamical distance. Each point is averaged over 1000 realizations; error bars are smaller than the symbol size. Inset: schematic illustration of the six regions (adapted from del Genio et al. 27 ). Layer 1 is individually stable only when σ 1 is greater than a critical value (red striped regions); layer 2 is individually stable only when σ 2 lies between two critical values (blue striped regions); region 1 is the only zone of the phase diagram where both layers are already individually stable. Fig. 2 Synchronization stability of networks with multiple interaction layers is virtually always mean-field. In networks with N nodes per layer, for normalized dynamical distances ℓ 2,1 /(N − 1) smaller than 0.1, the accuracy of the estimate of the maximum Lyapunov exponent Λ max improves with system size. Also, the error in the estimate never increases beyond 0.3, even for very large dynamical distances. Each point is averaged over 1,000 realizations; when not drawn, error bars are smaller than the symbol size. Inset: our mean-field theory assesses synchronization stability always correctly, with possible occasional exceptions only occurring in very small systems.
cases are however beyond the scope of the present work, and they will be explored in future studies.
In conclusion, we developed a mean-field theory of synchronization stability for networks with multiple interaction layers. In principle, the same approach could be applied to different types of networks, most likely resulting in similar but different equations. While the theory has been derived under the assumption of quasi-identical layers, we have shown that its range of validity and applicability includes the case of very different layers. In fact, our theory provides an accurate assessment of synchronization stability in networks whose layers are actually substantially different from each other. Moreover, the accuracy of the predictions increases with system size, raising the question of whether, in the thermodynamic limit, the linear stability of the globally synchronized state becomes a pure mean-field effect. In addition, the numerical complexity of our approach is lower than that of the exact solution of the problem. These considerations make the application of our theory particularly attractive in the case of large natural systems, whose studies have been so far frustrated, if not completely inhibited, and which may now become tractable both analytically and computationally.

Methods
Network structures. The simulations performed on Erdős-Rényi used a probability of occurrence of edges p = 0.4. This choice was made to guarantee a high probability that the layers were connected, which was tested for in each case, and which is a necessary condition for the existence of a globally synchronized state. Simulations on layers with heterogeneous topology were carried out on preferential-attachment layers with the same density. After creating the reference layer, the second layer and third layer were obtained by perturbing the edges of the first using a doubling-bisecting scheme to target the desired normalized dynamical distance with a tolerance of 10 −4 .
Layer dynamics. The local dynamics of Rössler oscillators is described by x ¼ ðx 1 ; x 2 ; x 3 Þ; F x ð Þ ¼ ðÀx 2 À x 3 ; x 1 þ ax 2 ; b þ ðx 1 À cÞx 3 Þ T . Here, we chose a = b = 0.2 and c = 9 to ensure chaotic local dynamics. We also let the interaction functions for the layers be H 1 x ð Þ ¼ ð0; x 2 ; 0Þ and H 2 x ð Þ ¼ ðx 1 ; 0; 0Þ, because these choices are known to create a rich phase diagram with 6 distinct regions of behaviour determined by the combinations of synchronization stability of the individual layers 27 . The third layer in the 3-layer simulations was given interaction function H 3 x ð Þ ¼ ðx 1 ; 0; 0Þ.
ODE integration. The differential equation systems were integrated using a Runge-Kutta-Fehlberg 4(5) method, with step size 0.01 and tolerance 10 −6 . The value of each initial component of the projected global synchronization error vector transverse to the synchronization manifold was chosen independently from a Gaussian distribution. This guarantees that, upon normalization, the initial projected transverse synchronization error vector was a uniform random unit vector, due to the spherical symmetry of multivariate normal distributions.
Estimate of maximum Lyapunov exponents. To compute the maximum Lyapunov exponents, after a transient time of 50, we evolved the systems for 500 windows of 100 integration steps each. After each window, we computed the logarithm of the norm of the components of the projected global synchronization errors transverse to the synchronization manifold, and normalized them back to a unit vector. Their averages provided estimates for the maximum Lyapunov exponents sought.