¹¹institutetext: Amity School of Engineering and Technology, Amity University, Bangalore, India. ¹¹email: [email protected] ²²institutetext: Department of Computer Science Engineering, Shiv Nadar Institute of Eminence, Delhi-NCR, India. ²²email: [email protected] ³³institutetext: Ghani Khan Choudhary Institute of Engineering and Technology, Malda, West Bengal, India ³³email: [email protected]

Effect of Perturbation and Topological Structure on Synchronization Dynamics in Multilayer Networks

Rajesh Kumar 11 Suchi Kumari 22 Anubhav Mishra 33

Abstract

The way the topological structure transforms from a decoupled to a coupled state in multiplex networks has been extensively studied through both analytical and numerical approaches, often utilizing models of artificial networks. These studies typically assume uniform interconnections between layers to simplify the analytical treatment of structural properties in multiplex networks. However, this assumption is not applicable for real networks, where the heterogeneity of link weights is an intrinsic characteristic. Therefore, in this paper, link weights are calculated considering the node’s reputation and the impact of the inter-layer link weights are assessed on the overall network’s structural characteristics. These characteristics include synchronization time, stability of synchronization, and the second-smallest eigenvalue of the Laplacian matrix (algebraic connectivity). Our findings reveal that the perturbation in link weights (intra-layer) causes a transition in the algebraic connectivity whereas variation in inter-layer link weights has a significant impact on the synchronization stability and synchronization time in the multiplex networks. This analysis is different from the predictions made under the assumption of equal inter-layer link weights.

Keywords:

Multiplex network Algebraic connectivity Perturbation Synchronization Time scale Stability.

1 Introduction

One fundamental property of connected identical dynamical systems is synchronization. Within multiplex networks, synchronization dynamics span multiple layers and occur within individual network layers. The algebraic connectivity, commonly represented by the second smallest eigenvalue ( $\lambda_{2}$ ) of the multiplex network’s Supra-Laplacian matrix, is closely associated with these synchronization patterns. Research on synchronization dynamics is crucial for various real-world systems in information technology, engineering, biology, physics, and social science [29]. Therefore, in this research, we examined how the algebraic connectivity ( $\lambda_{2}$ ) in multiplex networks is influenced by the topological structure. We then demonstrated the stability of synchronization processes in weighted multiplex networks, followed by an analysis of synchronization time.

The spectral characteristics of the Supra-Laplacian matrix in multiplex networks are crucial for understanding synchronization dynamics. Unlike the adjacency matrix, the eigenvalues of the Laplacian matrix provide clearer insights [22]. Specifically, algebraic connectivity, indicated by the second smallest eigenvalue ( $\lambda_{2}$ ), offers valuable information about the network’s modularity [24] and synchronization capacity [1]. Synchronization dynamics in multiplex networks are influenced by variables such as average node degree, clustering coefficients in each layer, and inter-layer link weights. Erratic shifts in algebraic connectivity often signify faster diffusion due to additional paths between nodes in the multilayer structure [9, 17]. For single-layer networks, synchronization stability is quantified using the eigenratio $R=\lambda_{N}/\lambda_{2}$ , where $\lambda_{2}$ and $\lambda_{N}$ are the second smallest and largest eigenvalues of the Laplacian matrix [6]. Radicchi et al. [25] explored how changes in the second smallest eigenvalue of the supra-Laplacian matrix lead to the emergence of two distinct regimes and a structural transition phase. This observation highlights the significant influence of structural factors on network dynamics.

The effect of inter-layer link patterns on the dynamics of spreading processes in an interconnected network has been studied by researchers in some recent work. Wang et al. [31] demonstrated that inter-layer connections based on node degrees have a relatively smaller effect on the size of infection compared to the density of interconnections. In another study, the authors investigated the correlation between intra-layer and inter-layer degrees [26]. They found that a strong correlation between these degrees can lead to the outbreak state, even if the epidemic threshold is not reached. When dealing with multiplex networks, the influence of their topological properties can become more complex. Simply examining a network layer in isolation, without considering its interactions with other layers, can lead to misleading conclusions.

Aguirre et al. [2] investigated the impact of the connector node degree on the synchronizability of two-star networks with one inter-layer link. They showed that synchronization could be achieved by connecting the high-degree nodes of each network. Xu et al. [32] studied the synchronizability of two-layer multiplex networks for three different coupling patterns. They determined that there exists an optimal value of the inter-layer coupling strength for maximizing complete synchronization. Authors in [27] demonstrated the impact of topological similarity among multiplex layers on synchronization performance. The distance between the pairs of layers is used to measure the similarity between the network layers [4]. They fixed one network layer while rewiring the network layers. The authors discovered that a sizeable inter-layer coupling generally promotes global multiplex synchronization for the fixed intra-layer and inter-layer diffusion coefficients.

In networks like Small World (SW) and Scale-Free (SF) characteristics, larger values of the average clustering coefficient ( $\langle CC\rangle$ ) can hinder global synchronization. This occurs because the network tends to divide into clusters that oscillate at different frequencies [21]. Let us take the example of interconnected Facebook and Twitter social networks. When viewing both networks as interconnected multiplex networks, it has an impact on their structural characteristics. Assume that a node $i$ (part of the connected triangle) has a clustering coefficient value of $1$ in a single independent network layer. However, when node $i$ is connected to another node, say $j$ , in a different network layer, it alters the clustering coefficient of node $i$ . This demonstrates that the topological structures of multiplex networks differ from those of their individual network layers. In the above mentioned research work, in impact of ( $\langle CC\rangle$ ) on the synchronization is studied on the independent networks. However, real-world systems consist multiple subsystems which are inter-connected and inter-dependent e.g., physical system, infra-structure systems, biological system etc. Therefore, in the presented research work, effect of average clustering coefficient ( $\langle CC\rangle$ ) on the stability of the synchronization in the multiplex networks is analyzed. Besides, edge weight in the networks is an abstraction and can be assigned in multiple ways. However, in real world networks, strength of the relationship is function of nodes connecting the edges. For example, number of interaction may be accounted as edge weights. In this research work, a method is proposed to calculated the edge weights which is function of the trust relationship among the pair of nodes. Hence, depending upon the network structure of individual layers, the behavior of $R$ in the multiplex networks changes with the variation in inter-layer link weights [28]. In certain cases, a decrease in the value of eigenratio $R$ = ( $\lambda_{N}/\lambda_{2}$ ) indicates that the two layers of the multiplex network begin to synchronize and behave as if they were a single-layer network [28]. Their investigations showed that for various topologies, for weaker inter-layer coupling weight, the value of parameter $R$ of the multiplex network is approximately the same as the value of $R$ of the individual network layer having the highest value of $\lambda_{2}$ of the Laplacian matrix. For decisive inter-layer coupling weights parameter, $R$ for the multiplex-network is approximated by taking the ratio of $\lambda_{N}$ and $\lambda_{2}$ of the average Laplacian matrix of the layers of the multiplex-network. These findings are a direct consequence of the strong inter-layer coupling strength between two layers[28]. In the previous research works, the behavior of stability has been studied by considering the unweighted and undirected multiplex networks [11, 28]. However, real multiplex networks [18] are often weighted and exhibit distinct topological properties. Hence, in this paper, we address this gap to study the synchronization dynamics in weighted multiplex networks of different topological structures by introducing the perturbation to the inter-layer edge weights. Also, Algebraic connectivity [8] is closely related to synchronization processes in multiplex networks and is used to study diffusion characteristics [11], synchronization ability [1], and modularity [24]. The spectral characteristics of the combinatorial Supra-Laplacian matrix significantly influence the dynamics of these synchronization processes. Therefore, impact of perturbation on the Algebraic connectivity [8] is analyzed by considering the different topological structures of synthetic multiplex (Random and Power law) networks and real world multiplex networks [19].

Major contributions of the presented research work are as follows:

•

A method is proposed to compute the intra-layer and inter-layer link weights.
•

Investigate the effect of perturbation (in the network layers) on the algebraic connectivity ( $\lambda_{2}$ ) for multiplex networks with different topological structures.
•

Analyze the stability of synchronization dynamics on considered multiplex networks by varying the inter-layer link weights.
•

Study the time scale (level of synchronization) for the multiplex networks.
•

Compare the findings of unweighted and weighted synthetic multiplex networks as well as empirical dataset multiplex networks

The remainder of the paper is organized as follows: Section 2 provides information on the proposed methods for computing weights. Section 3 describes the details of the simulation setup used in our study. The analysis and results are presented in Section 4. Finally, Section 5 concludes with the research findings and suggests possible directions for future investigation.

2 Proposed Methodology

A multiplex network is represented as a pair $\mathcal{M}=(\mathcal{G},\mathcal{E})$ , where $\mathcal{G}$ = $\{G_{\alpha};\alpha\in$ $\{1,2,\dots\ ,M\}$ $\}$ is a set of $M$ directed/undirected, weighted/unweighted graphs $G_{\alpha}=(\mathcal{V}_{\alpha},\mathcal{E}_{\alpha})$ (called network layers of $\mathcal{M}$ ) [16]. Each layer of Multiplex networks contains same set of nodes ( $\mathcal{V}_{\alpha}=\mathcal{V}_{\beta}=\mathcal{V}$ ) and each node in network layer $G_{\alpha}$ is connected to its replica node in network layer $G_{\beta}$ i.e., $\mathcal{E}_{\alpha\beta}=\{(i^{\alpha},i^{\beta});i^{\alpha}\in\mathcal{V}_{\alpha},i^{\beta}\in\mathcal{V}_{\beta}$ for every $1\leq\alpha\neq\beta\leq M$ . A good example is social multiplex networks where a group of users are present in each network layer (Facebook, Twitter etc.) shown in Fig. 1, but these users may have different friends in each network layer.

Refer to caption — Figure 1: Schematic representation of multiplex network composed by online social networks [5]

In matrix notation, multiplex networks are represented as supra-adjacency matrix $\mathcal{A}^{\mathcal{M}}$ .

\mathcal{A}^{\mathcal{M}}=\left[\begin{array}[]{c|c}\mathcal{A}^{[\alpha]}&\mathcal{A}^{[\alpha\beta]}\\ \hline\cr\mathcal{A}^{[\alpha\beta]}&\mathcal{A}^{[\beta]}\par\end{array}\right]

where $\mathcal{A}^{[\alpha]}=a_{ij}^{\alpha}\in\mathbb{R}^{N_{\alpha}\times N_{\alpha}}$ is adjacency matrix of network layer $G_{\alpha}$ [7] with element,

a_{ij}^{\alpha}=\left\{\begin{array}[]{ll}1\hskip 5.69054ptif\hskip 2.84526pt(e_{i,j}^{\alpha})\in{E}_{\alpha}\\ 0\hskip 5.69054ptotherwise\end{array}\right\}

and $N_{\alpha}$ is the number of nodes in the network layer $G_{\alpha}$ . The inter-layer adjacency matrix $\mathcal{A}^{[\alpha\beta]}$ with element $(a_{ij}^{\alpha\beta})\in\mathbb{R}^{N_{\alpha}\times N_{\beta}}$ [7] such that,

a_{ij}^{\alpha\beta}=\left\{\begin{array}[]{ll}1\hskip 5.69054ptif\hskip 2.84526pt(e_{i,j}^{\alpha,{\beta}})\in{E}_{\alpha{\beta}}\\ 0\hskip 5.69054ptotherwise\end{array}\right\}

In this section, we first describe the method to compute link weights, and then the stability of synchronization and synchronization time of the entire multiplex network is discussed.

2.1 Intra-layer and inter-layer Edge Weights Calculation

Network edge weights can be conceptually interpreted in various ways depending on the which type pf network we are considering. For example, in social networks, tracking individual interactions over a specified period might provide a representative weight or indicator of link strength. In such networks, the edge weights can be calculated by combining the node’s trust value and the degree of connectivity between nodes in the network. The node’s trust score can be evaluated using the well-known statistical method, Maximum Likelihood Estimation. Th evaluate the overall trust value ( $\phi^{\mathcal{M}}_{i}$ ) associated with a particular node $i$ , we define matrices $\mathcal{S}^{\mathcal{M}}$ and $\mathcal{F}^{\mathcal{M}}$ , of dimensions ( $M\times N$ ), where $M$ is the number of layers in the multi-layer network and $N$ is the total number of nodes. These matrices record the number of “successful” and ”failed” transactions between nodes $i$ and $j$ within the multi-layer network. The formulation is given as follows:

s_{ij}=\begin{cases}T^{\alpha}_{ij}+T^{\alpha\beta}_{ij},&\text{if there is an edge between}\hskip 5.69054pti^{\alpha}\text{and}j^{\alpha},i^{\alpha}\text{and}j^{\beta}:\alpha\neq\beta\\ 0,&\text{otherwise}\end{cases}

f_{ij}=\begin{cases}U^{\alpha}_{ij}+U^{\alpha\beta}_{ij},&\text{if there is an edge between}\hskip 5.69054pti^{\alpha}\text{and}j^{\alpha},i^{\alpha}\text{and}j^{\beta}:\alpha\neq\beta\\ 0,&\text{otherwise}\end{cases}

where ( $s_{ij}\in\mathcal{S}^{\mathcal{M}}$ , $f_{ij}\in\mathcal{F}^{\mathcal{M}}$ ) are the counts of successful and failed transactions between nodes $i$ and $j$ . The terms $\psi^{\mathcal{M}}_{i}=\displaystyle\sum_{i=1,i\neq j}^{M\times N}s_{ij}$ and $\theta^{\mathcal{M}}_{i}=\displaystyle\sum_{i=1,i\neq j}^{M\times N}f_{ij}$ hold the aggregate of successful and failed transaction so that likelihood function of ( $\phi^{\mathcal{M}}_{i}$ ) is given in Eq. (1).

L(\phi^{\mathcal{M}}_{i},y_{i,1},y_{i,2},\dots,y_{i,MN})=\displaystyle\prod_{j=1,j\neq i}^{MN}({\phi^{\mathcal{M}}_{i}})^{\psi^{\mathcal{M}}_{i}}(1-\phi^{\mathcal{M}}_{i})^{\theta^{\mathcal{M}}_{i}}

(1)

Using Eq. (1), the Log likelihood function ( $\zeta$ ) of $\phi^{\mathcal{M}}_{i}$ can be recasted as,

\zeta(\phi^{\mathcal{M}}_{i})=\text{ln}\hskip 2.84526ptL(\phi^{\mathcal{M}}_{i})=\psi^{\mathcal{M}}_{i}\text{ln}\hskip 2.84526pt\phi^{\mathcal{M}}_{i}+\theta^{\mathcal{M}}_{i}\text{ln}\hskip 2.84526pt(1-\phi^{\mathcal{M}}_{i})

(2)

Likelihood expression is obtained by taking derivative of the Eq. (2) and represented in Eq. (3)

\frac{d\zeta(\phi^{\mathcal{M}}_{i})}{d\phi^{\mathcal{M}}_{i}}=\frac{\psi^{\mathcal{M}}_{i}}{\phi^{\mathcal{M}}_{i}}-\frac{\theta^{\mathcal{M}}_{i}}{1-\phi^{\mathcal{M}}_{i}}=0

(3)

The derivative in Eq. (3) provides the the trust score of node $i^{\mathcal{M}}$ using Maximum Likelihood Estimation of $(\phi^{\mathcal{M}}_{i}$ ) to maintain quality of services (QoS). After evaluating the trust score of each node in the multiplex network, the edge weight can be evaluated in inter-layer as well as in intra-layer.

2.1.1 Computation of intra-layer and inter-layer edge weights:

Connection patterns ( $a_{ij}^{\alpha}$ ) and trust values ( $\phi^{\alpha})$ of the nodes are taken into account while calculating intra-layer and inter-layer link weights ( $w_{ij}$ ). The link weights are calculated using $\Gamma$ function in Eq. (5).

w_{ij}^{\alpha}=\mathcal{F}(\Gamma,a_{ij}^{\alpha})=\Gamma+a_{ij}^{\alpha}

(4)

\Gamma(\phi_{i}^{\alpha},\phi_{j}^{\alpha})=\Big{[\Big{(}}\frac{\text{Cos}(\Delta_{ij}^{\alpha})+1}{2}\Big{)}\times(\phi_{i}^{\alpha}\times\phi_{j}^{\alpha})\Big{]}

(5)

The $\Delta$ function is used to identify like-minded nodes in the network. In Eq. (5), $\Delta_{ij}^{\alpha}$ is evaluated by taking the difference ( $|\phi_{i}^{\alpha}-\phi_{j}^{\alpha}|$ ) of the trust scores of nodes $i$ and $j$ . Since the trust score is a probability function, its value will lie in the range $[0,1]$ . The $\cos(\Delta_{ij})$ value, however, will be any real number in the range $[-1,1]$ . To normalize the $\cos(\Delta_{ij})$ score within the range $[0,1]$ , it is transformed to $\frac{\cos(\Delta_{ij}^{\alpha})+1}{2}$ .

Let us examine two different situations:

(i)

$\phi_{i}^{\alpha}=0.1$ and $\phi_{j}^{\alpha}=0.1$ .
(ii)

$\phi_{i}^{\alpha}=0.99$ and $\phi_{j}^{\alpha}=0.99$ .

Since the computed trust difference in both scenarios is $\Delta_{ij}^{\alpha}=0$ , the value of $\frac{\cos(\omega)+1}{2}$ is $1$ . Although both scenarios provide the same score, the trust scores are noticeably greater in the second scenario compared to the first. To rectify this disparity, a scaling factor is implemented, and the trust values are adjusted accordingly. The trust values for nodes $i$ and $j$ are multiplied ( $\phi_{i}^{\alpha}\times\phi_{j}^{\alpha}$ ) by the normalized cosine function $(\frac{\cos(\omega)+1}{2})$ to ensure the trust values are adjusted in a way that appropriately represents their magnitudes.

2.2 Stability of Synchronization

Suppose that each node $i$ at any layer $L^{\alpha}$ has some amount of information $x_{i}^{\alpha}(t)$ at some time $t$ . Assume that this information flows from node $j^{\alpha}$ to node $i^{\alpha}$ at rate $p(x_{j}^{\alpha}-x_{i}^{\alpha})$ in network layer $L^{\alpha}$ and from node $i^{\alpha}\in L^{\alpha}$ to node $i^{\beta}\in L^{\beta}$ with rate $d_{x}(x_{i}^{\beta}-x_{i}^{\alpha})$ , where $p$ and $d_{x}$ are the parameters controlling the intra-layer and inter-layer weights respectively. The dynamical equation [11] characterizing the evolution of the state of node $x_{i}^{\alpha}$ at network layer $L^{\alpha}$ in a multiplex network having $M$ networked layers is given in Eq. (6).

\frac{dx_{i}^{\alpha}}{dt}=p\sum_{j=1}^{N}w_{ij}^{\alpha}(x_{j}^{\alpha}(t)-x_{i}^{\alpha}(t))+d_{x}\sum\limits_{\begin{subarray}{c}\beta=1\\ \alpha\neq\beta\end{subarray}}^{M}w_{ii}^{\alpha\beta}(x_{i}^{\beta}(t)-x_{i}^{\alpha}(t))\hskip 14.22636pt

(6)

where $w_{ij}^{\alpha}$ and $w_{ii}^{\alpha\beta}$ are the edge weights of intra-links and inter-links, respectively. The parameters $p$ and $d_{x}$ are considered as tuning parameters to control $w_{ij}^{\alpha}$ and $w_{ij}^{\alpha\beta}$ , respectively. The rate equation in Eq. (6) can be represented in the matrix form in Eq. (7).

\dot{\textbf{x}}=(p\mathcal{L}^{L}+d_{x}\mathcal{L^{I}})\textbf{x}=\mathcal{L^{M}}\textbf{x}

(7)

where $\mathcal{L}^{L}$ and $\mathcal{L^{I}}$ depicts the supra-Laplacians of the intra-layer and inter-layer networks, respectively and $d_{x}$ controls the inter-layer link weights [11]. The Laplacian matrix of each network layer $\mathcal{L}^{\alpha}$ is formulated as $\mathcal{L}^{\alpha}$ = $\Delta^{\alpha}$ - $W^{\alpha}$ where $W^{\alpha}$ is the weight matrix of network layer $L^{\alpha}$ with elements $w_{ij}^{\alpha}>0$ if there is an edge between $i^{\alpha}$ and $j^{\alpha}$ , otherwise $w_{ij}^{\alpha}=0$ . The matrix ${(\Delta^{\alpha})}$ is diagonal matrix with elements $\delta_{ii}=\sum_{j}w_{ij}^{\alpha}$ . Considering both the parameters ${\Delta^{\alpha}}$ and $W^{\alpha}$ , the matrix $\mathcal{L}^{L}$ is provided in Eq. (8).

\mathcal{L}^{L}=\begin{pmatrix}p^{1}\mathcal{L}^{(1)}&0&\cdots&0\\ 0&p^{2}\mathcal{L}^{(2)}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&p^{M}\mathcal{L}^{(M)}\ \end{pmatrix}\vspace{3mm}

(8)

and inter-layer weight matrix is given in Eq. (9)

\mathcal{W^{I}}=\begin{pmatrix}0&\mathcal{W}_{(1,2)}I&\mathcal{W}_{(1,3)}I&\cdots&\mathcal{W}_{(1,M)}I\\ \mathcal{W}_{(2,1)}I&0&\mathcal{W}_{(2,3)}I&\cdots&\mathcal{W}_{(2,M)}I\\ \mathcal{W}_{(3,1)}I&\mathcal{W}_{(3,2)}I&0&\cdots&\mathcal{W}_{(3,M)}I\\ \vdots&\vdots&&\ddots&\vdots\\ \mathcal{W}_{(M,1)}I&X_{(M,2)}I&X_{(M,3)}I&\cdots&0\ \end{pmatrix}

(9)

where $\mathcal{W}_{(\alpha\beta)}$ is a vector of inter-layer edge weights { $w_{ii}^{\alpha\beta}$ } connecting $i^{\alpha}\in L^{\alpha}$ and $i^{\beta}\in L^{\beta}$ , $I$ is the identity matrix. Inter-layer Laplacian $\mathcal{L^{I}}$ matrix can be determined from $\mathcal{W}^{I}$ . The state of a node $i$ in Eq. 6 can be found using $x_{i}(t)=x(0)e^{-\lambda_{i}t}$ [23, 12], where $\lambda_{i}$ is the eigenvalue of supra-Laplacian matrix $\mathcal{L^{L}}$ . Apart from that, the second smallest eigenvalue ( $\lambda_{2}$ ) of the Laplacian matrix is considered as algebraic connectivity of the graph $G$ because of the following proposition:

Proposition: Let $G$ = $(V,E)$ be a connected graph with positive weights $w_{ij}$ , then the algebraic connectivity of graph $G$ is positive and equal to the minimum of the function [8]

\zeta(x)=N\frac{\sum_{i,j\in E}w_{ij}(x_{i}-x_{j})^{2}}{\sum_{i,j\in E,i<j}w_{ij}(x_{i}-x_{j})^{2}}

(10)

over all non-constant $N$ -tuples $x=(x_{i})$ .

Stability of the synchronization is determined by the parameter $R=\frac{\lambda_{N}}{\lambda_{2}}$ [6], where $\lambda_{N}$ and $\lambda_{2}$ are the largest and smallest non-zero eigenvalues of the matrix $\mathcal{L^{M}}$ . For the lower values of $d_{x}$ , $\lambda_{2}$ of $\mathcal{L^{M}}$ is approximated by $\lambda_{2}(\mathcal{L})\approx\lambda_{2}(\mathcal{L^{I}})$ . Thus the behavior of $R$ can be understood by the approximation [28],

R\approx\frac{max(\lambda_{N}(\mathcal{L}^{\alpha})+s_{i}^{\mathcal{I}})}{\lambda_{2}(\mathcal{L^{I}})d_{x}}

(11)

where $\mathcal{L}^{\alpha}$ is the Laplacian matrix of layer $L^{\alpha}$ , $\lambda_{N}$ is the maximum eigenvalue of $\mathcal{L}^{\alpha}$ and $s_{i}^{\mathcal{I}}=\sum\limits_{\begin{subarray}{c}\beta=1\\ \alpha\neq\beta\end{subarray}}^{M}w_{ii}^{\alpha\beta}$ is the strength of the node $i$ in the inter-layer network having maximum eigenvalue [28].
For the strong inter-layer link weights ( $d_{x}>>1$ ), $R$ is approximated as,

R\approx\frac{d_{x}\lambda_{M}(\mathcal{L^{I}})+\lambda_{N}(\mathcal{L}_{max}^{WA})}{\lambda_{2}(\mathcal{L}^{AV})}

(12)

where $\mathcal{L}^{AV}$ is the average Laplacian of network layers and $\mathcal{L}_{max}^{WA}$ is given by

\mathcal{L}_{max}^{WA}=\frac{1}{\parallel{\textbf{x}}^{{}^{\prime}\mathcal{I}}\parallel^{2}}\sum_{\alpha}(x_{\alpha}^{{}^{\prime}\mathcal{I}})\mathcal{L}^{\alpha}

(13)

where ${\textbf{x}}^{{}^{\prime}\mathcal{I}}$ is the eigenvector corresponding to maximum eigenvalue of inter-layer network matrix[28].

Synchronization Time

Let us assume that each node of the network layer is embedded within an oscillator and the level of synchronization on the $MN$ interaction units can be expressed by quantity $S(\tau)$ for a large timestamps $\tau$ in Eq. (14).

S(\tau)=\frac{1}{MN}\sum_{i=1}^{MN}(1-x_{i}(\tau))

(14)

where $x_{i}(\tau)=x_{i}(0)e^{-\lambda_{i}\tau}$ and $\lambda_{i}$ is the $i^{th}$ eigenvalue of the supra-Laplacian matrix. The value of synchronization score $S\rightarrow 1$ , when the system is completely synchronized and for sufficiently large timestamps $\tau$ . The corresponding time is known as synchronization time $\tau_{s}$ , starting from random phase $x_{i}$ [13] of node $i$ in any network layer.

3 Simulation Setup

For the simulation, two versions of synthetic multiplex networks with two layers are considered. The first is constructed using the Barabasi-Albert model of preferential attachment [3], with network layers $L_{1}$ and $L_{2}$ . Each layer contains 200 nodes, with $m=2$ for layer $L_{1}$ and $m=3$ for layer $L_{2}$ , where $m$ is the number of stubs with which a new node attaches to existing nodes in the network layers. The second network is designed using the Power law network model ( $p(k)\sim k^{-\gamma}$ ), with $\gamma=2.1$ for layer $L_{1}$ and $\gamma=2.2$ for layer $L_{2}$ . Each layer in this network also contains 200 nodes. After constructing the network layers, inter-layer connections are established so that each node in layer $L_{1}$ is connected to its replica node in layer $L_{2}$ , as shown in Fig. 1. For validation, a multiplex network consisting of two layers is designed using the empirical dataset from the CS-Aarhus social multiplex network [19]. The parameters related to the dataset are shown in Table 1. Two network layers, Facebook and Lunch, are considered for analysis due to their single connected component.

Table 1: Network layer parameters of the dataset CS-Aarhus social multiplex network [19]

Name of network layer	Nodes	Edges	Conn. Comp.	Average node degree
Work	61	194	2	6.47
Leisure	61	88	1	3.74
Coauthor	61	21	8	1.68
Lunch	61	193	1	6.43
Facebook	61	124	1	7.75

To analyze the effect of perturbation on the inter-layer links in the behavior of algebraic connectivity ( $\lambda_{2}$ ), stability of the synchronization and synchronization time $(\tau)$ , following scenario is considered for unweighted and weighted multiplex networks.

•

For weighted multiplex networks, intra-layer edge weights $w_{ij}^{\alpha}$ and inter-layer edge weights $w_{ii}^{\alpha\beta}$ are computed according to the method discussed in Section 2.1.1
•

For unweighted case, values of $w_{ij}^{\alpha}$ and $w_{ii}^{\alpha\beta}$ are taken as $1$ .
•

For weighted and unweighted multiplex networks, value of $d_{x}$ and $p$ is taken from $0$ to $2$ in the step-size of $0.001$ .

4 Result and Analysis

In this section, we analyze the effect of perturbations on algebraic connectivity, synchronization stability, and synchronization time in the multiplex network.

4.1 Algebraic Connectivity

We study the behaviour of algebraic connectivity ( $\lambda_{2}$ ) for the considered synthetic as well empirical data set multiplex network (weighted and unweighted). For $d_{x}$ = 0, spectrum of the supra-Laplacian matrix of the multiplex network is given by $\wedge$ ( $\mathcal{{L}^{M}}$ ) = $\{0=\lambda_{1}=\lambda_{2}<\lambda_{3}\leq\lambda_{4}\dots\lambda_{2N}\}$ with $\lambda_{3}$ = min( $\lambda_{2}^{1}$ , $\lambda_{2}^{2}$ ) where $\lambda_{2}^{1}$ and $\lambda_{2}^{2}$ are the second smallest eigen values of Laplacian matrices of first and second network layers of multiplex networks for unweighted and weighted cases [12]. The spectrum of the inter-layer Laplacian matrix is given by $\wedge(\mathcal{L}^{I})={0,2d_{x}}$ . In the absence of inter-layer edges, the two network layers are isolated, resulting in $\lambda_{1}$ and $\lambda_{2}$ both being 0 in the spectrum $\wedge(\mathcal{L}^{M})$ . To study the effect of inter-layer edge weights on the evolution of algebraic connectivity, we consider the following synthetic multiplex networks.

4.1.1 Network Layers of Multiplex Network Designed Using BA Model

Simulation results for $\lambda_{2}(d_{x},p)$ in the unweighted case are plotted in Fig. 2. It is observed that for lower values of $(p=0.2)$ , $\lambda_{2}$ grows monotonically (as a straight line) with increasing $d_{x}$ , and $\lambda_{2}$ saturates at $d_{x}>0.550$ . However, when the network layers are perturbed, the effect of $d_{x}$ becomes significant, and $\lambda_{2}(d_{x}>0.550,p>0.4)$ experiences an abrupt increase and does not stabilize even at $d_{x}=2$ . In the regime where $d_{x}>0.550$ , $\lambda_{2}$ is controlled by a monotonically increasing function that stabilizes at $\lambda_{2}$ of ${(\mathcal{L}^{1}+\mathcal{L}^{2})/2}$ , with ${\mathcal{L}^{1}+\mathcal{L}^{2}}$ being the Laplacian matrix of the aggregated network, as shown in Fig. 2. For the weighted case, not much change is observed in the behavior of $\lambda_{2}(d_{x},p)$ , except at $d_{x}\approx 0.450$ , where $\lambda_{2}(d_{x},p)$ shows an abrupt change, as shown in Fig. 3.

4.1.2 Network Layers of Multiplex Network Designed Using Power Law Network Model

For the unweighted case, it is evident from the simulations results that similar to the BA model case, $\lambda_{2}$ grows monotonically and gets saturated at ( $d_{x}>0.8,p\approx 0.2$ ) as depicted in Fig. 4. Introducing more perturbation to the network layers results in a sudden change in $\lambda_{2}$ at lower values of $d_{x}$ (Fig. 4 shows an abrupt change at $\lambda_{2}(d_{x}\approx 0.1,p\approx 2)$ ). Additionally, variations in the values of ( $d_{x},p$ ) compared to the BA model result in deviations in $\lambda_{2}$ . For the BA model, $\lambda_{2}(d_{x}=2,p=2)\approx 1.8$ , while for the Power law network, $\lambda_{2}(d_{x}=2,p=2)\approx 0.6$ . This difference arises due to the distinct topological structures, specifically the higher average clustering coefficients observed in Table 2. In the weighted case, $\lambda_{2}$ grows and, after a quick change at ( $d_{x}>0.5,p\approx 0.2$ ), stabilizes, as shown in Fig 4. The effect of perturbation (increasing $p$ in the network layers) on $\lambda_{2}$ is similar to the unweighted case, with a slight shift in $\lambda_{2}(d_{x},p)$ . However, unlike the uniform weight case, $\lambda_{2}(d_{x}=2,p=2)\approx 0.38$ .

Indeed, in Small World (SW) and Scale-Free (SF) networks, large values of clustering produce a hindrance to global synchronization since the network partitions into clusters that oscillates at different frequencies[21].

4.1.3 Multiplex Network Constructed From CS-Aarhus Social Multiplex Network [19] Dataset

In the unweighted case, $\lambda_{2}$ shows an abrupt change at lower values of $d_{x}\approx 0.1$ for all considered values of $p$ , and it stabilizes at $d_{x}\approx 0.8$ , achieving a maximum value of $\lambda_{2}\approx 0.24$ , as shown in Fig. 6. The appearance of uniform colors (Yellow, Saffron, Green, Blue) indicates that perturbations in the network layers have little effect, and $\lambda_{2}$ attains nearly the same values for given values of $d_{x}$ and $p$ , as represented in Fig. 6. However, in the weighted case, $\lambda_{2}$ shows an abrupt change at $d_{x}\approx 0.1$ and stabilizes at $\lambda_{2}\approx 1.8$ for $d_{x}>0.1$ at $p=0.2$ . Unlike the unweighted case, variations in $\lambda_{2}$ (different color shades in Fig. 7) are observed for the given values of $p$ and $d_{x}$ , implying that perturbations affect the network layers of the multiplex network, as represented in Fig. 7.

Simulation results obtained from synthetic and empirical multiplex network datasets reveal that the topological structure ( $\langle CC\rangle$ ), intra-layer perturbations in link weights, and inter-layer link weights all affect the algebraic connectivity.

4.2 Stability of Synchronization

In this section, we analyze the effect of inter-layer edge weights (controlled by the parameter $d_{x}$ ) and the average clustering coefficient ( $\langle CC\rangle$ ) on the stability of the synchronization process ( $R$ ) for the considered weighted and unweighted multiplex networks. A lower value of $R$ indicates maximum stability [6]. The optimal value of $R$ for the multiplex network is approximated by $d_{x}$ at which analytical curves plotted using Equations (11) and (12) intersect each other.

When two independent network layers are interconnected to form a multiplex network, the values of $\langle CC\rangle$ differ from those in the individual layers. This occurs due to the addition of new edges between nodes in one layer and their replicas in the other layer. These changes are shown in Table 2 for single-layer networks and multiplex networks. Additionally, there is a variation of $\lambda_{2}$ against $\langle CC\rangle$ . An increase in the value of $\langle CC\rangle$ leads to a decline in the value of $\lambda_{2}$ , as represented in Fig. 8. This occurs because the network splits into dynamical clusters, with nodes within each cluster being more connected compared to the connections among different clusters.

Table 2: Values of average clustering coefficient

(\langle CC\rangle)^{L}

(the individual layers),

(\langle CC\rangle)^{M}

\lambda_{2}^{M}

\lambda_{N}^{M}

and

k^{M}_{max}

of multiplex networks constructed using BA model and Empirical data-set of CS-Aarhus.

Parameters $\Longrightarrow$	${\langle CC\rangle}^{L_{1}}$	${\langle CC\rangle}^{L_{2}}$	${\langle CC\rangle}^{M}$	$\lambda_{2}^{M}$	$\lambda_{N}^{M}$	$k_{max}$
BA model	0.079	0.085	0.05	1.42	48.37	50
Power Law model	0.17	0.13	0.07	0.48	120.8	122
CSA	0.67	0.63	0.41	0.76	28.14	28

4.2.1 Network Layers of Multiplex Network Designed Using BA Model

The values of parameter $R$ obtained from simulation are plotted (Red color curve) against $d_{x}$ for the considered unweighted and weighted multiplex networks, as shown in Figs. 9 (a) and (b), respectively. Theoretical values of $R$ approximated by Equations (11) and (12) are represented by green and blue curves, respectively, in the same figure. As $d_{x}$ increases, the values of $R$ obtained using Equations (11) and (12) decrease and increase continuously, respectively. It is observed from Figs. 9 (a) and (b) that the optimal value of $R\approx 45$ is obtained at $d_{x}=1$ for the unweighted case, and $R\approx 49$ is achieved at $d_{x}=0.7$ for the weighted case.

Interestingly, This behavior is due to the properties of $\lambda_{2}$ as demonstrated in Figs. 2 and 3. For the initial (lower) values of $d_{x}$ , both $\lambda_{2}$ and $\lambda_{2}^{\prime}$ increase, which decreases $R$ . However, at higher values of $d_{x}$ and lower values of $p=0.2$ , $\lambda_{2}$ stabilizes (as shown in Figs. 2 and 3) while $\lambda_{N}$ of the supra-Laplacian matrix continues to increase. Thus, we find that weighted and unweighted multiplex networks differ in terms of $d_{x}$ and $R$ . An interesting point to note that the values of $R$ against $d_{x}$ first decrease and, after reaching a minimum (optimal) value, begin to increase. This behavior is due to the properties of $\lambda_{2}$ as demonstrated in Figs. 2 and 3. For the initial (lower) values of $d_{x}$ , both $\lambda_{2}$ and $\lambda_{2}^{\prime}$ increase, which decreases $R$ . However, at higher values of $d_{x}$ and lower values of $p=0.2$ , $\lambda_{2}$ stabilizes (as shown in Figs. 2 and 3) while $\lambda_{N}$ of the supra-Laplacian matrix continues to increase. Thus, we find that weighted and unweighted multiplex networks differ in terms of $d_{x}$ and $R$ .

Table 3: Represents the optimum values of controlling parameters

d_{x}

for which eigenratio (

R

) is minimum for the BA model and Empirical dataset.

Parameters $\Longrightarrow$	$d_{x}(Unwtd)$	$R(Unwtd)$	$d_{x}(Wtd)$	$R(Wtd)$
BA model	1	44.48	0.7	48.67
Power Law model	0.5	177	0.38	235.61
CSA	0.75	24.98	0.25	66.25

4.2.2 Network Layers of Multiplex Network Designed Using Power Law Network Model

In this case, optimal value of $R\approx 177$ is obtained at $d_{x}\approx 0.5$ for unweighted cases as well as $R\approx 235$ is achieved at $d_{x}\approx 0.38$ for weighted case as shown in the Figs. 10 (a) and (b) respectively. It occurs due to a difference in the topological characteristics as compared to the multiplex network using the BA model. It is observed from Figs. 4 and 5 that for the given values of $d_{x}$ and $p$ , values of $\lambda_{2}$ is lower in the case of the Power Law Network model as compared to the BA model due to the effect of average clustering coefficient ( $\langle CC\rangle$ ). At the same time, from Table 2 we find that $\lambda_{N}$ is higher (approx. three times) in the case of multiplex network constructed using Power Law Network model as compared to BA model. Thus, the difference in $\lambda_{N}$ and $\lambda_{2}$ produces variation in the parameter $R$ .

4.2.3 Multiplex Network Constructed From CS-Aarhus Social Multiplex Network [19] Dataset

It is observed from Figs. 9, 10 and 11 that behaviour of $R$ is approx. same for the synthetic as well as empirical dataset multiplex networks. However, some variation in optimal values $R$ and $d_{x}$ is observed for weighted as well as unweighted cases. In this case, optimal values of $R\approx 25$ is obtained at $d_{x}\approx 0.75$ for unweighted case and $R\approx 66$ is achieved at $d_{x}\approx 0.25$ as shown in Fig. 11. Here, again the variation in $R$ and $d_{x}$ is noticed as compared to a multiplex network constructed using the BA model and Power Law Network model. The reason is the variation ib the topological characteristics of the network layers of the multiplex network.

Simulation results reveal that the stability of the synchronization $R$ is lower in the case of the Power Law Network model as compared to the BA model. However, multiplex networks constructed from the empirical dataset show higher stability $R$ than considered artificial multiplex networks. This variation is due to the topological characteristics (size, sparsity of intra-layer and inter-layer edges, etc.) of the network layers constituting multiplex networks. It is also observed from simulation results that there is variation in $R$ in the cases of the considered unweighted and weighted multiplex network.

4.3 Time of Synchronization

Synchronization time $(\tau)$ is when all the layers of the multiplex networks are fully synchronized $i.e$ when level of synchronizations $S\rightarrow 1$ at sufficiently large $\tau$ . For the considered multiplex networks, values of the parameter $R$ are plotted against time $\tau$ for $d_{x}=0.3$ and $d_{x}=20$ . For the simulation purpose, these values of $d_{x}$ are chosen to quantify lower and higher inter-layer link weights. However, any value of $d_{x}$ can be chosen.

4.3.1 Network Layers of Multiplex Network Designed Using BA Model

The parameter $S$ obtained from Equation (14) is plotted against time-scale $\tau$ as shown in Fig. 12. For the lower values of $d_{x}=0.3$ , the values of $S$ against $\tau$ are nearly same for weighted as well as unweighted cases. At higher $d_{x}=20$ , the parameter $S\rightarrow 1$ at $\tau\simeq 0.2$ whereas for $d_{x}=0.3$ , $S\rightarrow 1$ at $\tau\simeq 0.4$ as demonstrated in Fig. 12. Higher values of $d_{x}$ is an indication of increasing the frequency of inter-layer interaction between the nodes present at the different network layers of the considered multiplex network. Thus, we find that increasing the value of $d_{x}$ lead to establish the synchronization in multiplex networks more quickly.

4.3.2 Network Layers of Multiplex Network Designed Using Power Law Network Model

In this case, for $d_{x}=0.3$ and $d_{x}=20$ the parameter $S$ lags in approaching 1 for weighted case as compared to unweighted cases. It occurs due the fact that in weighted case, non-uniform link weights (intra-layer and inter-layer) cause change in the state (information) of the nodes in the multiplex such that $S\rightarrow 1$ at later time $\tau$ as compared to unweighted case. For $d_{x}=0.3$ , $S\rightarrow 1$ at $\tau\approx 0.4$ for weighted as well unweighted cases. However, for $d_{x}=20$ , $S\rightarrow 1$ at $\tau\approx 0.3$ and $\tau\approx 0.4$ for weighted and unweighted respectively as shown in Fig. 13.

4.3.3 Multiplex Network Constructed From CS-Aarhus Social Multiplex Network [19] dataset

In this case, the behaviour of parameter $S$ against $\tau$ is nearly identical to the BA model. For $d_{x}=0.3$ , $S\approx 0.7$ and for $d_{x}=20$ , $S\approx 0.7$ at $\tau=2$ for weighted and unweighted cases as in Fig. 14. With the increase in time $\tau$ , the gap between the values of $S$ keeps on decreasing and finally $S\rightarrow 1$ at $\tau=6$ for $d_{x}=0.3,20$ . Thus, from Fig. 14 we find that increasing the $d_{x}$ results in synchronization of the considered multiplex network at earlier time $\tau$ .

5 Conclusion and Future Work

In this paper, we examined the impact of perturbations in network layers on the algebraic connectivity ( $\lambda_{2}$ ) and synchronization dynamics in unweighted and weighted multiplex networks, each comprising two network layers. Our simulation results revealed that for both synthetic and real dataset-based multiplex networks, $\lambda_{2}$ varies with the given value of $d_{x}$ and different values of $p$ across the considered individual multiplex networks and different topological multiplex networks. It occurs due the perturbation in the intra-layer link weights (individual multiplex networks) and the differences in the average clustering coefficient ( ${\langle CC\rangle}^{CSA}>{\langle CC\rangle}^{PowerLawNetwork}>{\langle CC\rangle}^{BA}$ ). As far as stability of synchronization is concerned, multiplex networks constructed using the BA model show more stability than the Power Law Network model for weighted and unweighted cases. Maximum stability (minimum value of $R$ ) is obtained at $d_{x}=1,0.5$ for BA and Power Law Network models, respectively, for unweighted cases. Multiplex network constructed from CSA dataset shows the highest stability at $d_{x}=0.75$ for unweighted case. This variation instability occurs due to the characteristics of the network layer. In the case of a synthetic network, although the size (number of nodes) of network layers is the same, these differ in $\langle CC\rangle$ , which produces variation in $\lambda_{2}$ and $\lambda_{N}$ of the supra-Laplacian matrix of the multiplex network thereby showing variation in the parameter $R$ .

Our simulation results show that multiplex networks are slightly less stable in the weighted case compared to the unweighted case due to non-identical intra-layer and inter-layer edge weights. For synchronization time $\tau$ , the parameter $S\rightarrow 1$ is reached earlier when $d_{x}=20$ compared to $d_{x}=0.3$ in both weighted and unweighted cases. The time $\tau$ for $S\rightarrow 1$ is slightly higher in the weighted case for the Power Law Network model and CSA-dataset multiplex networks. This indicates that networks with higher average clustering coefficients ( $\langle CC\rangle$ ) are less stable and take longer to achieve complete synchronization. However, tuning inter-layer edge weights can improve synchronization stability and reduce the time required. Therefore, the characteristics ( $\langle CC\rangle$ ) of individual layers and inter-layer link weights significantly influence the synchronization process. In the future, this work can be extended to analyze graph energy in the synchronization process.

Acknowledgement

We want to thank Dr Anurag Singh, Associate Professor at the Department for Computer Science and Engineering, National Institute of Technology Delhi, India for providing guidance to work in this domain during the doctorate work.

References

[1] Aguirre, J., Sevilla-Escoboza, R., Gutiérrez, R., Papo, D., Buldú, J.: Synchronization of interconnected networks: the role of connector nodes. Physical review letters 112(24), 248701 (2014)
[2] Aguirre, J., Sevilla-Escoboza, R., Gutiérrez, R., Papo, D., Buldú, J.: Synchronization of interconnected networks: the role of connector nodes. Physical review letters 112(24), 248701 (2014)
[3] Albert, R., Barabási, A.L.: Statistical mechanics of complex networks. Reviews of modern physics 74(1), 47 (2002)
[4] Andrade, R.F., Miranda, J.G., Pinho, S.T., Lobão, T.P.: Measuring distances between complex networks. Physics Letters A 372(32), 5265–5269 (2008)
[5] de Arruda, G.F., Rodrigues, F.A., Moreno, Y.: Fundamentals of spreading processes in single and multilayer complex networks. Physics Reports 756, 1–59 (2018)
[6] Barahona, M., Pecora, L.M.: Synchronization in small-world systems. Physical review letters 89(5), 054101 (2002)
[7] Boccaletti, S., Bianconi, G., Criado, R., Del Genio, C.I., Gómez-Gardenes, J., Romance, M., Sendina-Nadal, I., Wang, Z., Zanin, M.: The structure and dynamics of multilayer networks. Physics Reports 544(1), 1–122 (2014)
[8] De Abreu, N.M.M.: Old and new results on algebraic connectivity of graphs. Linear algebra and its applications 423(1), 53–73 (2007)
[9] De Domenico, M., Granell, C., Porter, M.A., Arenas, A.: The physics of spreading processes in multilayer networks. Nature Physics 12(10), 901 (2016)
[10] Gatti, P.L.: Probability theory and mathematical statistics for engineers. CRC Press (2004)
[11] Gomez, S., Diaz-Guilera, A., Gomez-Gardenes, J., Perez-Vicente, C.J., Moreno, Y., Arenas, A.: Diffusion dynamics on multiplex networks. Physical review letters 110(2), 028701 (2013)
[12] Gomez, S., Diaz-Guilera, A., Gomez-Gardenes, J., Perez-Vicente, C.J., Moreno, Y., Arenas, A.: Diffusion dynamics on multiplex networks. Physical review letters 110(2), 028701 (2013)
[13] Guardiola, X., Diaz-Guilera, A., Llas, M., Pérez, C.: Synchronization, diversity, and topology of networks of integrate and fire oscillators. Physical Review E 62(4), 5565 (2000)
[14] Gupta, K., Ambika, G.: Role of time scales and topology on the dynamics of complex networks. Chaos: An Interdisciplinary Journal of Nonlinear Science 29(3), 033119 (2019)
[15] Kivelä, M., Arenas, A., Barthelemy, M., Gleeson, J.P., Moreno, Y., Porter, M.A.: Multilayer networks. Journal of complex networks 2(3), 203–271 (2014)
[16] Kivelä, M., Arenas, A., Barthelemy, M., Gleeson, J.P., Moreno, Y., Porter, M.A.: Multilayer networks. Journal of complex networks 2(3), 203–271 (2014)
[17] Kumar, R., Kumari, S., Bala, M.: Minimizing the effect of cascade failure in multilayer networks with optimal redistribution of link loads. Journal of Complex Networks 9(6), cnab043 (2021)
[18] Kumar, R., Kumari, S., Mishra, A.: Robustness of multilayer networks: A graph energy perspective. Physica A: Statistical Mechanics and its Applications 628, 129160 (2023)
[19] Magnani, M., Micenkova, B., Rossi, L.: Combinatorial analysis of multiple networks. arXiv preprint arXiv:1303.4986 (2013)
[20] Martínez, J.H., Boccaletti, S., Makarov, V.V., Buldú, J.M.: Multiplex networks of musical artists: the effect of heterogeneous inter-layer links. arXiv preprint arXiv:1805.08711 (2018)
[21] McGraw, P.N., Menzinger, M.: Clustering and the synchronization of oscillator networks. Physical Review E 72(1), 015101 (2005)
[22] Mohar, B., Alavi, Y., Chartrand, G., Oellermann, O.: The laplacian spectrum of graphs. Graph theory, combinatorics, and applications 2(871-898), 12 (1991)
[23] Newman, M.: Networks. Oxford university press (2018)
[24] Newman, M.E.: Modularity and community structure in networks. Proceedings of the national academy of sciences 103(23), 8577–8582 (2006)
[25] Radicchi, F., Arenas, A.: Abrupt transition in the structural formation of interconnected networks. Nature Physics 9(11), 717 (2013)
[26] Saumell-Mendiola, A., Serrano, M.Á., Boguná, M.: Epidemic spreading on interconnected networks. Physical Review E 86(2), 026106 (2012)
[27] Serrano, A.B., Gómez-Gardeñes, J., Andrade, R.F.: Optimizing diffusion in multiplexes by maximizing layer dissimilarity. Physical Review E 95(5), 052312 (2017)
[28] Sole-Ribalta, A., De Domenico, M., Kouvaris, N.E., Diaz-Guilera, A., Gomez, S., Arenas, A.: Spectral properties of the laplacian of multiplex networks. Physical Review E 88(3), 032807 (2013)
[29] Strogatz, S.H.: From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D: Nonlinear Phenomena 143(1-4), 1–20 (2000)
[30] Székely, G.J., Alpár, M., Unger, É.: Paradoxes in probability theory and mathematical statistics. Tech. rep., Springer (1986)
[31] Wang, Y., Xiao, G.: Effects of interconnections on epidemics in network of networks. In: Wireless Communications, Networking and Mobile Computing (WiCOM), 2011 7th International Conference on. pp. 1–4. IEEE (2011)
[32] Xu, M., Zhou, J., Lu, J.a., Wu, X.: Synchronizability of two-layer networks. The European Physical Journal B 88(9), 1–6 (2015)