Optimizing spreading dynamics in interconnected networks

TL;DR

This paper develops a theoretical strategy for adding multiple interconnecting edges between networks to optimize spreading dynamics, outperforming simpler methods near the critical point.

Contribution

It introduces a novel optimal edge-adding strategy for interconnected networks based on SIS model analysis, improving spreading prevalence and outbreak threshold.

Findings

The strategy maximizes spreading prevalence near the critical point.

It minimizes the outbreak threshold for interconnected networks.

Numerical simulations confirm superior performance over baseline methods.

Abstract

Adding edges between layers of interconnected networks is an important way to optimize the spreading dynamics. While previous studies mostly focus on the case of adding a single edge, the theoretical optimal strategy for adding multiple edges still need to be studied. In this study, based on the susceptible-infected-susceptible (SIS) model, we investigate the problem of maximizing the stationary spreading prevalence in interconnected networks. For two isolated networks, we maximize the spreading prevalence near the critical point by choosing multiple interconnecting edges. We present a theoretical analysis based on the discrete-time Markov chain approach to derive the approximate optimal strategy. The optimal inter-layer structure predicted by the strategy maximizes the spreading prevalence, meanwhile minimizes the spreading outbreak threshold for the interconnected network simultaneously. Numerical simulations on synthetic and real-world networks show that near the critical point, the proposed strategy gives better performance than connecting large degree nodes and randomly connecting.

Tables

Table 1. Table 1: Some basic statistics of four real-world networks. The statistics includes the number of nodes ( N 𝑁 N ), the number of edges ( M 𝑀 M ), the maximal degree ( k max subscript 𝑘 max k_{\mathrm{max}} ), the first moment of the degree distribution ( ⟨ k ⟩ delimited-⟨⟩ 𝑘 \langle k\rangle ), the clustering coefficient ( c 𝑐 c ), the theoretical epidemic threshold predicted by λ c ∗ = 1 / ω 1 superscript subscript 𝜆 𝑐 1 subscript 𝜔 1 \lambda_{c}^{*}=1/\omega_{1} and the modularity ( Q 𝑄 Q )  Newman ( 2006 ) .

	$N$	$M$	$k_{\max }$	$⟨k⟩$	$c$	${\lambda }_c^{\ast }$	$Q$
Advogato	5042	39227	803	15.56	0.092	0.014	0.337
Facebook	2888	2981	769	2.0644	0.0004	0.036	0.809
HepPh	34401	420784	846	24.463	0.280	0.013	0.408
HepTh	27400	352021	2468	25.695	0.269	0.009	0.328

Equations

$$G^{0}=\left({\begin{array}[]{cc}G_{a}&0\\ 0&G_{b}\end{array}}\right).$$

$$\delta G=\left({\begin{array}[]{cc}0&C\\ C^{\mathrm{T}}&0\end{array}}\right).$$

$$⟨1_{N_a},C{1}_{N_b}⟩=\delta M,$$

$${\rho }_i(t+1)=(1-\mu ){\rho }_i(t)+\left(1-{\rho }_i(t)\right)\left(1-q_i(t)\right)$$

$$q_i(t)=\prod _{j=1}^N\left(1-\lambda G_{ij}{\rho }_j(t)\right).$$

$${\rho }_i=\underset{t\to \infty }{\lim }{\rho }_i(t).$$

$$G\rho =(\mu /\lambda )\rho .$$

$$⟨\rho ⟩\approx \left({\lambda }^{\ast }{\omega }_1-1\right)\frac{{\sum }_{i=1}^N{u}_i}{N{\sum }_{i=1}^N{u}_i^3}$$

$$\{{\omega }_k^a:k=1,\cdots ,N_a\}\cup \{{\omega }_l^b:l=1,\cdots ,N_b\}.$$

$${\hat{v}}_k^a={\left(v_k^a,0\right)}^{\mathrm{T}},$$

$${\hat{v}}_l^b={\left(0,v_l^b\right)}^{\mathrm{T}}$$

$$u^0={\beta }^a{\hat{v}}_1^a+{\beta }^b{\hat{v}}_1^b$$

$$u=c\underset{n\to \infty }{\lim }{\left({\omega }_1^a\right)}^{-n}{G}^n{u}^0$$

$$u\approx \frac{c}{{\omega }_1^a}\left({\beta }^a{\omega }_1^a{\hat{v}}_1^a+{\beta }^b{\omega }_1^b{\hat{v}}_1^b+{\beta }^a\delta G{\hat{v}}_1^a+{\beta }^b\delta G{\hat{v}}_1^b\right).$$

$$c{\beta }^a\to {\beta }^a,\frac{c}{{\omega }_1^a}{\beta }^b{\omega }_1^b\to {\beta }^b,$$

$$u\approx u^0+\delta u,$$

$$\delta u=\frac{{\beta }^a}{{\omega }_1^a}\delta G{\hat{v}}_1^a+\frac{{\beta }^b}{{\omega }_1^b}\delta G{\hat{v}}_1^b.$$

$$\left(G^0+\delta G\right)u^0=\left({\omega }_1^a+\delta {\omega }_1\right)u^0.$$

$$G^0{u}^0={\beta }^a{\omega }_1^a{\hat{v}}_1^a+{\beta }^b{\omega }_1^b{\hat{v}}_1^b,$$

$$\left(\delta G-\delta {\omega }_1\mathbb{I}\right)u^0={\beta }^{b}g{\hat{v}}_1^b,$$

$$\begin{array}{rl}& {\beta }^b⟨v_1^a,C{v}_1^b⟩=\delta {\omega }_1{\beta }^a,\\ & {\beta }^a⟨v_1^a,C{v}_1^b⟩-g{\beta }^b=\delta {\omega }_1{\beta }^b.\end{array}$$

$$\left({\begin{array}[]{cc}0&E(C)\\ E(C)&-g\end{array}}\right)\left({\begin{array}[]{c}\beta^{a}\\ \beta^{b}\end{array}}\right)=\delta\omega_{1}\left({\begin{array}[]{c}\beta^{a}\\ \beta^{b}\end{array}}\right)$$

$$\delta {\omega }_1=\frac{1}{2}\left(\sqrt{4E^2(C)+g^2}-g\right).$$

$$\begin{array}{rl}& {\beta }^a=\sqrt{1+{\left(\frac{g}{2E(C)}\right)}^2}+\frac{g}{2E(C)},\\ & {\beta }^b=1.\end{array}$$

$$\begin{array}{rl}& {\omega }_1={\omega }_1^a+\delta {\omega }_1,\\ & u={\beta }^a{\hat{v}}_1^a+{\beta }^b{\hat{v}}_1^b,\end{array}$$

$$\begin{array}{rl}⟨\rho ⟩=& \left({\lambda }^{\ast }{\omega }_1^a+{\lambda }^{\ast }\delta {\omega }_1-1\right)\\ & \times \frac{\left({\left({\beta }^a\right)}^3+{\beta }^a\right){\theta }_1^a+\left({\left({\beta }^a\right)}^2+1\right){\theta }_1^b}{N\left({\left({\beta }^a\right)}^3{\theta }_3^a+{\theta }_3^b\right)}\end{array}$$

$${\theta }_1^a=⟨1_{N_a},v_1^a⟩,{\theta }_3^a=⟨1_{N_a},v_1^a\cdot v_1^a\cdot v_1^a⟩$$

$$\Delta =\sqrt{\overline{⟨\rho {⟩}^2}/{\overline{⟨\rho ⟩}}^2-1},$$

$$\cos \theta =⟨u^{\ast },u⟩.$$

Full text

Spreading dynamics in interconnected networks relay on the inter-layer structure apart from the structure within each layer. For two given networks, how to design the interconnecting structure to optimize the spreading dynamics is a very appealing topic. Previous studies obtained the optimal strategy when considering adding a single edge in two-layer interconnected networks, while the optimal strategy of adding multiple edges lacks theoretical studies. Therefore in this study, a novel strategy is proposed to promote the spreading dynamics by adding multiple interconnecting edges for two isolated networks. Near the critical point, the spreading prevalence can be written in terms of the leading eigenvalue and the corresponding eigenvector of the adjacency matrix. Basing on an approximation scheme for the leading eigenvalue and eigenvector of the interconnected network, we optimize the spreading prevalence among all candidate edges. The optimal inter-layer structure is achieved by selecting edge candidates that are top-ranked by the product of eigenvector centrality of nodes in its two ends. Meanwhile, the optimal strategy also minimizes the outbreak threshold of the interconnected network. Numerical simulations on three pairs of synthetic networks and two pairs of real-world networks show that the strategy gives better performance than the heuristic strategies of connecting large degree nodes and of randomly connecting, especially near the critical point.

I Introduction

In real-world social systems, individuals might communicate with others via multiple possible channels (such as Twitter, Facebook, LinkedIn). The communication relations with a specific channel can be represented by a network, where the nodes correspond to individuals, and the edges correspond to the communication relations. Therefore all the communication relations combined can be described by a multilayer network Boccaletti et al. (2014); Gao et al. (2012); Kivelä et al. (2014); De Domenico et al. (2016), where each layer corresponds to one of the communications channels. Multilayer networks and the dynamics on them have attracted attentions from diverse areas. It has been observed that multilayer networks display distinct collective behaviors compared to that of single-layer networks Cai et al. (2018); Solé-Ribalta et al. (2016); Buono and Braunstein (2015); Buldyrev et al. (2010); Cozzo et al. (2013); Kim and Goh (2013); Lee et al. (2012); Wang et al. (2015); Zhang et al. (2015). As an example, percolation processes on multilayer networks display first-order phase transitions Buldyrev et al. (2010); Valdez et al. (2013); Di Muro et al. (2016a, b); Valdez et al. (2016), which is intrinsically different from the second-order phase transitions on single-layer networks Cohen et al. (2002); Serrano and Boguná (2006); Goltsev et al. (2008). For evolutionary games, multilayer networks promote cooperation better than single-layer networks Perc et al. (2017); Wang et al. (2013a, 2015). For synchronization processes, explosive synchronization and hysteresis loop are observed on multilayer networks Zhang et al. (2015).

Spreading dynamics on multilayer networks have attracted considerable attention in recent studies De Domenico et al. (2016); Wang et al. (2019); Brummitt et al. (2012); Lee et al. (2014); Yağan and Gligor (2012); Wang et al. (2018a, b); Shu et al. (2018); Chen et al. (2018). Saumell-Mendiola, Serrano and Boguná Saumell-Mendiola et al. (2012) have studied susceptible-infected-susceptible (SIS) model Pastor-Satorras and Vespignani (2001) on multilayer networks. They found that adding a small fraction of edges between layers can lead to the outbreak of epidemics while without these edges the epidemic extinct. The epidemic threshold depends on the structure of multilayer networks Cozzo et al. (2013); Wang et al. (2013b), and it is possible to observe localization phenomenons de Arruda et al. (2017). For susceptible-infected-recovered (SIR) model, Dickison, Havlin and Stanley Dickison et al. (2012) found that the system might exhibit a mixed phase, i.e., the epidemic outbreaks in one layer but not others. Refs. Wang et al. (2014, 2016a); Liu et al. (2016) studied a model with the spreading of information and epidemics simultaneously. The studies found that the diffusion of information can inhibit the spreading of epidemics.

The inter-layer structure in multilayer networks has significant impacts on the dynamics. Parshani et al. Parshani et al. (2011) found that a positive inter-layer degree correlation will inhibit large scale cascading failures. Ref. Saumell-Mendiola et al. (2012) found that those degree correlations will make the epidemic outbreak more easily, while Ref. Wang et al. (2014) found that the positive inter-layer degree correlations can also inhibit the outbreak of the epidemic. Understanding what kind of interconnecting structure will lead to better performance for specific dynamics is an essential task for understanding dynamics on multilayer network and for designing better network structures. For spreading dynamics, when considering adding a single edge between layers, the optimal solution was given in Ref. Aguirre et al. (2013) analytically. Based on matrix perturbation theory, Ref. Aguirre et al. (2013) derive that connecting the two nodes with the largest eigenvector centralities in each layer minimizes the epidemic threshold while maximizes the spreading prevalence. For better synchronizability, Aguirre et al. Aguirre et al. (2014) studied the optimal strategy when adding one single edge between layers analytically. Based on matrix perturbation theory, they found that connecting large degree nodes will give better synchronizability. Li et al. Li et al. (2016) further generalized the optimal strategy when adding two edges. When consider adding multiple interconnecting edges, current results are mostly based on numerical methods. Wei et al. Wei et al. (2018a) studied the interconnecting strategy numerically when adding a small number of edges for two-layer networks with random regular networks in each layer. Their studies suggest that adding inter-layer connections gives a more significant contribution to synchronizability compared to inner-layer connections. Wei et al. Wei et al. (2018b) also did numerical simulations for the optimal strategy for general multilayer networks.

In this study, we investigate the problem of optimizing the spreading prevalence in two-layer networks by adding multiple inter-layer edges. For two isolated networks, we try to understand how to add a small number of edges to maximize the stationary spreading prevalence in the interconnected network. We mainly focus on the SIS model near the critical point. With a known formula given in Goltsev et al. (2012), the epidemic prevalence near the critical point can be written in terms of the leading eigenvalue and eigenvector. We first develop a scheme for approximating the new leading eigenvalue and eigenvector for the interconnected network after adding those interconnecting edges. With this approximation scheme, we obtain a formula that predicts the stationary epidemic prevalence in the two-layer interconnected network. Then this approximated prevalence can be optimized among all possible inter-layer structures. The optimal inter-layer structure that maximizes the spreading prevalence will found to minimizes the spreading outbreak threshold for the interconnected network simultaneously. Numerical simulations are performed to compare the strategy with some other heuristic strategies. The proposed strategy gives a better performance at least near the critical point when adding a small number of interconnecting edges.

The rest of the paper is organized as follows. In Sec. II, we introduce the basic setups of the model and some notations. In Sec. III, the theoretical derivations of the strategy is given. Then in Sec. IV, the strategy is tested and compared to some other heuristic strategies and finally in Sec. V we give some conclusions and discussions.

II Model

Starting with two isolated networks $a$ and $b$, we add a fixed number of edges to interconnect the two networks. The way of adding the edges will affect the dynamical behaviors on the interconnected network. We focus on the way of adding these edges that maximizing the spreading prevalence.

Let the adjacency matrices of the networks $a,b$ be $G_{a}$ and $G_{b}$ respectively. The number of nodes in $a(b)$ is $N_{a}(N_{b})$ and the number of edges is $M_{a}(M_{b})$. The total number of nodes is denoted by $N=N_{a}+N_{b}$ and the total number of edges by $M=M_{a}+M_{b}$. The adjacency matrix of the combined network then is

[TABLE]

A set of edges with fixed cardinality $\delta M$ will be added between the two networks. After the operation the new network combining $a,b$ and interconnecting edges will have adjacency matrix $G=G^{0}+\delta G$, where

[TABLE]

Here $C$ is an $N_{a}\times N_{b}$ matrix which indicates how the inter layer connections are added. Its elements take values $C_{ij}\in\{0,1\}$, where $C_{ij}=1$ if and edge is added between the $i$-th node of network $a$ and $j$-th node of network $b$ and $C_{ij}=0$ otherwise. The matrix satisfies the constraint on total number of added edges as

[TABLE]

where $\mathbf{1}_{N_{a}}=\left[1,\cdots,1\right]^{T}$ is all-one vector of length $N_{a}$ and $\langle\cdot\rangle$ is the inner product of two vectors.

By choosing among all possible assignment of $C$, the epidemic prevalence can be maximized. In this paper, we consider the SIS model. For the SIS model, each node is in either the susceptible or infected state. In a discrete-time setting, at time step $t$, infected nodes have a probability $\lambda$ to infect their susceptible neighbors independently. Then the infected nodes (not including nodes get infected at the current time step) become susceptible again with probability $\mu$. In large time limit, the density of infected nodes will converge to its stationary value. The target of the paper is to maximize the stationary spreading prevalence.

III Theory

We employ the discrete time Markov chain approach Gómez et al. (2010) to describe SIS model on general networks. For the discrete time Markov chain approach, the status of a node $i$ is characterized by $\rho_{i}(t)$, which is the probability that $i$ is infected at time step $t$. Then $\rho_{i}(t)$ evolves according to the following discrete time equations

[TABLE]

for $i\in\{1,\cdots,N\}$ and

[TABLE]

The first term on the r.h.s. of Eq. (4) corresponds to the probability that node $i$ is infected at $t$ and not recovered in $t+1$, and second term is the probability that $i$ is susceptible at $t$ and get infected by at least one infected neighbor. $q_{i}(t)$ in Eq. (5) is the probability that $i$ is not get infected by any of its neighbors at $t$.

The stationary solution is given by the limit in $t$

[TABLE]

Let $\rho$ and $q$ be the vector with elements $\rho_{i}$ and $q_{i}$ for $i\in\{1,\cdots,N\}$. Near the critical point $\rho$ is expected to be small and the equation can be linearized as $q\approx 1-\lambda G\rho$, the stationary equation reads

[TABLE]

The spreading outbreaks only when $\lambda/\mu>1/\omega_{1}$ where $\omega_{1}$ is the leading eigenvalue of $G$. It has been shown in Goltsev et al. (2012) the stationary prevalence of epidemic is approximately

[TABLE]

where $\left\langle\rho\right\rangle=(1/N)\sum_{i=1}^{N}\rho_{i}$, $\lambda^{*}=\lambda/\mu$ is the effective infection probability and $u$ is the eigenvector corresponding to the eigenvalue $\omega_{1}$. It can be seen that $\left\langle\rho\right\rangle$ is determined by the leading eigenvalue and eigenvector approximately near the critical point. Since $G$ is obtained by adding a small number of edges to $G^{0}$, its spectra should be closely related to that of $G^{0}$. Next we develop a scheme to approximate the spectra of $G$ from that of $G^{0}$.

First consider the spectra of matrix $G^{0}$. Let the $\omega^{a}_{k}$ for $k\in\{1,\cdots,N_{a}\}$ be the $k$-th eigenvalue of $G_{a}$ and $v^{a}_{k}$ the corresponding eigenvector. Similarly $\omega^{b}_{l}$, $v^{b}_{l}$ with $l\in\{1,\cdots,N_{b}\}$ are eigenvalues and eigenvectors for $G_{b}$. The adjacency matrix of networks $a$ and $b$ combined $G^{0}$ is a diagonal block matrix by putting $G_{a}$ and $G_{b}$ in the diagonal, thus with eigenvalues

[TABLE]

Clearly the leading eigenvalue of $G^{0}$ is $\max\{\omega^{a}_{1},\omega^{b}_{1}\}$. Without losing of generality we assume $\omega^{a}_{1}\geq\omega^{b}_{1}$. The corresponding eigenvector for an eigenvalue $\omega^{a}_{k}$ is

[TABLE]

which is by combining $v^{a}_{1}$ and all-zero vector of length $N_{b}$. Similarly, for eigenvalue $\omega^{b}_{l}$, the corresponding eigenvector of $G^{0}$ is

[TABLE]

with zero vector of length $N_{a}$.

Now we consider adding a small number of interconnecting edges. By assuming these edges won’t shift the spectra of the two isolated networks too much, a first approximation of $u$ would be proportional to $v_{1}^{a}$ in the first $N_{a}$ elements, and to $v_{1}^{b}$ in the last $N_{b}$ elements, thus written as

[TABLE]

where $\beta^{a},\beta^{b}\in\mathbb{R}$ are coefficients to be determined. The eigenvector $u$ is given by the following limit

[TABLE]

for some constant $c$ and $G^{n}$ denotes self matrix multiplication of $G$ for $n$ times. Since the number of interconnecting edges is small, it can be assumed that $u^{0}$ is already close to $u$ by choosing $\beta^{a}$ and $\beta^{b}$ properly. The limit thus approximated by setting $n=1$ and this gives

[TABLE]

Rescale the parameters by

[TABLE]

the eigenvector $u$ is approximated by the form

[TABLE]

where

[TABLE]

For scale-free networks which we mainly consider in the paper, it has been shown that the leading eigenvalue diverges in the thermodynamic limit Goh et al. (2001). Thus $\delta u$ can be ignored for large enough networks and $u$ is approximated as $u\approx u^{0}$.

Denote the gap by $g=\omega^{a}_{1}-\omega^{b}_{1}$. The leading eigenvalue $\omega_{1}$ of $G$ can be written as $\omega_{1}=\omega^{a}_{1}+\delta\omega_{1}$, which is $\omega^{a}_{1}$ plus a correction term $\delta\omega_{1}$. With these approximations we arrive at the following eigenvalue equation

[TABLE]

By definition

[TABLE]

and after some algebra Eq. (18) becomes

[TABLE]

where $\mathbb{I}$ is the $N$ by $N$ identity matrix. Multiplying $\hat{v}^{a}_{1}$ and $\hat{v}^{b}_{1}$ from the left separately gives the following equations

[TABLE]

The equations can be written in the form of an eigenvalue equation

[TABLE]

where we have denoted by $E(C)=\left\langle v_{1}^{a},Cv_{1}^{b}\right\rangle$. The equation gives two eigenvalues and pick the larger one which is

[TABLE]

The corresponding eigenvectors gives

[TABLE]

In summary, the leading eigenvalue and eigenvector of the interconnected network are approximated as

[TABLE]

where $\delta\omega_{1}$, $\beta^{a}$ and $\beta^{b}$ are given in Eq. (23) and Eq. (24). Clearly when $g/E(C)\to 0$, $\beta^{a}\to 1$ and the two networks play equal role. Meanwhile if $g/E(C)\gg 0$, network $a$ dominates $b$. Normalizing $u$ to unity and substituting back into (8) gives

[TABLE]

where

[TABLE]

and similarly for $\theta^{b}_{1}$, $\theta^{b}_{3}$.

With this approximation scheme, we can optimize $\left\langle\rho\right\rangle$ over $C$. Since $\left\langle\rho\right\rangle$ depends on the interconnecting matrix $C$ only in the form of $E(C)$, it is sufficient optimize $\left\langle\rho\right\rangle$ as a function of $E(C)$. After determining the optimal $E(C)$, we choose the matrix $C$ such that $E(C)$ is closest to its optimal value. Empirically, for all the networks considered in this study (see Sec. IV), $\left\langle\rho\right\rangle$ is an increasing function of $E(C)$ for small $E(C)$. Thus it is sufficient to choose $C$ such that $E(C)$ is maximized, or in other words, to choose the top $\delta M$ edges ranked by $v^{a}_{1}(i)\times v^{b}_{1}(j)$ for $i\in\{1,\cdots,N_{a}\},\ j\in\{1,\cdots,N_{b}\}$. Here $v^{a}_{1}(i)$ denotes the $i$th element of $v^{a}_{1}$ and similarly for $v^{b}_{1}(j)$. Note that a node can be connected to multiple nodes in the other layer. After determine $C$ in this way, we perform the check to ensure that $\left\langle\rho\right\rangle$ is non-decreasing in the region $(0,E(C)]$. From Eq. (23), $\delta\omega_{1}$ is also an increasing function of $E(C)$. Therefore, the optimal strategy also maximizes $\delta\omega_{1}$ as well as $\left\langle\rho\right\rangle$. Since the outbreak threshold for the interconnected network is given by $1/\left(\omega^{a}_{1}+\delta\omega_{1}\right)$, as a consequence, the optimal inter-layer structure also minimizes the outbreak threshold. In the rest of the paper, the proposed strategy is called large eigenvector connecting (LEC), since it is by connecting inter-layer node pairs with larger product of the two nodes’ eigenvector centrality.

IV Results

In this section we test the strategy proposed in Sec. III both for synthetic and real-world networks. For synthetic networks, we consider scale-free networks generated by uncorrelated configuration model (UCM). Denote by ‘SF $\alpha$’ for scale-free network with degree distribution $p(k)\sim k^{-\alpha}$ and degree exponent $\alpha$, we consider the following three pairs of networks: (i) SF $3.0$-SF $3.0$, (ii) SF $3.0$-SF $2.3$, (iii) SF $3.0$-SF $4.0$ and two pairs of real world networks: (iv) Advogato-Facebook, (v) HepPh-HepTh. Here Advogato Massa et al. (2009) and Facebook Leskovec and Mcauley (2012) are two online social networks, and HepPh and HepTh Leskovec et al. (2007) are citation networks of papers from the high energy physics-theory and high energy physics-phenomenology sections of the e-print arXiv. The networks are downloaded from Kon . TABLE 1 shows some basic statistics of the four real-world networks. When performing Monte-Carlo simulations on these networks, we first choose $\lceil 0.02N\rceil$ nodes as infected seeds. Then the spreading runs at discrete time steps. At each time step, all the infected nodes infect each of their susceptible neighbors with probability $\lambda$ independently and then returns to susceptible state with probability $\mu$. The steps are then repeated until the density of infected nodes in the network reaches the stationary value.

First, we verify the accuracy of epidemic prevalence predicted by Eqs. (4). Define the ratio $p=\delta M/M$. For the five network pairs considered, $\lceil pM\rceil$ of edges are added uniformly random between the two networks. Starting with random initial conditions, we iterate Eqs. (4) until reaching a tolerated error. We perform Monte-Carlo simulations on the networks generated to compare with the theoretical predictions. Fig. 1 shows how $\langle\rho\rangle$ changes as a function of $\lambda$. The recovery rate is fixed to $\mu=0.5$ and as well in the rest of the paper. The density of interconnecting edges is $p=0.01$. In the insets of Fig. 1, the theoretical epidemic threshold $\lambda_{c}^{*}=1/\omega_{1}$ is compared to the simulation results. The variability measure $\Delta$ Shu et al. (2015, 2016) is applied to determine the spreading threshold in simulations. In particular

[TABLE]

where the overline $\overline{\bullet}$ denotes the average over independent runs of the simulations. As a convention $\Delta=0$ when $\overline{\langle\rho\rangle}=0$. The simulated outbreak threshold is given by where the variability measure reaches its maximal. We mark the two thresholds by vertical lines both in Fig. 1 and the insets. From Fig. 1, for all the network pairs considered expect Advogato-Facebook (Fig. 1(d)), $\langle\rho\rangle$ and the epidemic threshold predicted by Eqs. (4) agree with simulations in good accuracy. Meanwhile, for the network pair Advogato-Facebook, predictions of Eqs. (4) are less accurate. Since mean-field theory ignores dynamical correlations Wang et al. (2017), i.e., ignores dependencies among the states of each node and its neighbors, it has been shown that mean-field theory could be less accurate for real-world networks with high clustering coefficient Radicchi and Castellano (2016), high modularity Wang et al. (2016b) and low average degree Gleeson et al. (2012). As shown in TABLE 1, the network Facebook has high modularity and low average degree, therefore, this could be the possible origins that Eqs. (4) are less accurate for Advogato-Facebook. In general, quenched mean-field theory still predicts well the simulations, therefore can provide a necessary guarantee for the accuracy of the theory developed in Sec. III.

Next we test the accuracy of predictions for the leading eigenvalue $\omega_{1}$ and eigenvector $u$ of the interconnected network. We consecutively add interconnecting edges between the two networks according to the LEC strategy and compare the eigenvalue and eigenvector predicted by Eqs. (25) to the true values $\omega_{1}^{*}$ and $u^{*}$. Start with two isolated networks, we consecutively add the edges and do the comparison after each adding. To quantify the accuracy, for the leading eigenvalue we consider the relative error $\left(\omega_{1}^{*}-\omega_{1}\right)/\omega_{1}^{*}$, where $\omega_{1}^{*}$ is the true eigenvalue and $\omega_{1}$ is that predicted by Eqs. (25). For leading eigenvector the accuracy is measured by the cosine similarity between the true one $u^{*}$ and predicted one $u$. Let $\theta$ be the angle between $u^{*}$ and $u$,

[TABLE]

Note that we have assumed both $u$ and $u^{*}$ are normalized to unity.

Still denote the fraction of edges added by $p=\delta M/M$. $\left(\omega_{1}^{*}-\omega_{1}\right)/\omega_{1}^{*}$ and $\cos\theta$ versus $p$ are shown in Figs. 2(a) and (c) respectively, while $\omega_{1}^{*}$ and $\omega_{1}$ versus $p$ are shown in Fig. 2(b). From Figs. 2, the theoretical predictions gives close approximations to it’s true value. As can observed in Figs. 2(a) and (b), errors in the eigenvalue $\omega_{1}^{*}-\omega_{1}$ have a small magnitude compared to the value of $\omega_{1}^{*}$, at least for small values of $p$. Also $\cos\theta$ is close to $1$ which indicates $u$ and $u^{*}$ are closely aligned. After tested the accuracy of Eqs. (4) and Eqs. (25), the interconnecting strategy proposed in Sec. III is expected to be reliable.

Next we test the performance of the strategy of optimizing spreading. For comparison, the two following heuristic strategies are considered. (i) Large degree connecting (LDC). This is by choosing top candidate edges that ranked by the total degree of its two ends $k_{i}+k_{j}$. (ii) Random connecting (RC). The edges are chosen uniformly random among all node pairs. Remind that the strategy proposed in Sec. III is called largest eigenvector connecting (LEC).

As the model depends on the parameter $p$, we first consider comparing the three strategies with fixed $p=0.005$. The results for different small values of $p$ are similar as will be shown in the next. After choosing $C$ according to the strategies, Eqs. (4) are iterated to stationary. For the convenience of visualization, the performance is evaluated by $\delta\langle\rho\rangle$, which is the difference of stationary $\langle\rho\rangle$ after and before adding the interconnecting edges in $C$. $\delta\langle\rho\rangle$ versus the infect rate $\lambda$ are shown in Fig. 3 for the five pairs of networks. Still, the recovery rate is fixed to $\mu=0.5$. The threshold value of $\lambda_{c}=\mu/\omega^{a}_{1}$ for the isolated network $G^{0}$ is denoted by a black vertical line. For $\lambda$ near the critical point, LEC strategy gives highest $\delta\langle\rho\rangle$ for all the five pairs of networks. For the three model network pairs in Figs. 3(a)-(c), LDC gives a very close performance to LEC. This is because that LEC and LDC predict similar $C$ for these three network pairs, as the degree centrality and eigenvalue centrality are strongly correlated in ranking for the UCM Lü et al. (2016). For real-world networks, the structure becomes complex; therefore, degree and eigenvector centrality are less correlated. Thus a better performance of LEC is observed near the critical point for the two real-world network pairs.

As $\lambda$ becomes large and deviates the critical value, those nodes with large eigenvalue centrality have a high probability of to be infected. Thus interconnecting these nodes becomes unnecessary. Meanwhile, high degree nodes still play a central role in maintaining the epidemic. Thus gradually, better performance of the LDC strategy is observed. Eventually, when $\lambda$ become very large, the infected state is prevalent; therefore, all central nodes have a high probability to be infected. Meanwhile, the RC strategy more likely picks low degree nodes due to the power-law degree distribution. Thus in this region of $\lambda$, RC is expected to give the best performance among the three.

To get a more comprehensive picture, we study how the performance of the three strategies depend on $p$ and $\lambda$. For each point in the parameter space of $p$ and $\lambda$, we compute $\langle\rho\rangle$ for the three strategies and find the differences of LEC minus LDC and LEC minus RC. In other words, a positive difference indicates that LEC performs better than the other one compared. The results are shown in Fig. 4. The first row in Fig. 4, i.e. Fig. 4(a1)-(e1), corresponds to LEC minus LDC for the five network pairs, and the second row (Figs. 4(a2)-(e2)) correspond to LEC minus RC. The contours of $\epsilon=1\mathrm{e}{-6}$ are shown in the plots to indicate the region when LEC performs better. Here $\epsilon$ is not set to be $\epsilon=0$ to eliminate possible numeric errors. In Fig. 4, the position of $\lambda_{c}=\mu/\omega^{a}_{1}$ for the unconnected network is marked by a vertical black dashed line. For all the cases in Fig. 4, there is a vertical band in the parameter space where LEC gives better performance. The vertical band locates in where $\lambda$ is close to the critical value. Thus we can conclude that LEC is better than the other two near the critical point when the number of added edges is small.

V Discussion

In the study, we investigated the optimal strategy for interconnecting two isolated networks in order to maximize the spreading prevalence of SIS model. We develop a scheme to approximate the leading eigenvalue and eigenvector for the interconnected network. This approximation gives a formula that predicts the epidemic prevalence for the interconnected network. By maximizing the spreading prevalence over the interconnecting matrix, we obtain the optimal inter-layer structure. By numerically iterating the discrete-time Markov equations, we find the strategy gives better performance than some other heuristic strategies.

By maximizing the spreading prevalence $\langle\rho\rangle$ among all interconnecting matrices, it turns out that the strategy is equivalent to select edge candidates that are top-ranked by the product of eigenvector centrality of nodes in its two ends. Actually, for the networks considered in this study, $\langle\rho\rangle$ is not always an increasing function of $E(C)$ for large $E(C)$ and large $\lambda$. Nevertheless, in this region, the number of edges added must be very large, and the infection probability is far from critical; as a consequence, the eigenvalue approximation scheme becomes inaccurate. What is a better strategy when adding a large number of edges cannot be answered by the current approach.

Besides, the current method mainly focuses on the region near the epidemic threshold. In summary, we mainly work in regions where linear approximations are reliable. For larger effective spreading rates and larger numbers of added edges, the problem might become nonlinear and more complex. If it is possible to get a more comprehensive theory in all parameter regions requires further studies.

Acknowledgements.

This work was partially supported by National Natural Science Foundation of China (Nos. 61673150, 11622538 and 61673086), China Postdoctoral Science Foundation (No. 2018M631073) and the Fundamental Research Funds for the Central Universities.