Optimal interlayer structure for promoting spreading of SIS model in two-layer networks

TL;DR

This paper introduces a strategy to enhance SIS epidemic spreading in multilayer networks by adding a single interlayer edge, using a derived index based on a modified Katz centrality, which outperforms traditional methods.

Contribution

The study develops a novel index for optimal interlayer edge placement to promote spreading in multilayer networks, validated through both small and large network experiments.

Findings

The proposed index effectively predicts the impact of interlayer edges on spreading.

Adding a single optimal interlayer edge significantly increases spreading prevalence.

The method outperforms traditional static strategies like degree and eigenvector centrality.

Abstract

Real-world systems, ranging from social and biological to infrastructural, can be modeled by multilayer networks. Promoting spreading dynamics in multilayer networks may significantly facilitate electronic advertising and predicting popular scientific publications. In this study, we propose a strategy for promoting the spreading dynamics of the susceptible-infected-susceptible model by adding one interconnecting edge between two isolated networks. By applying a perturbation method to the discrete Markovian chain approach, we derive an index that estimates the spreading prevalence in the interconnected network. The index can be interpreted as a variant of Katz centrality, where the adjacency matrix is replaced by a weighted matrix with weights depending on the dynamical information of the spreading process. Edges that are less infected at one end and its neighborhood but highly infected at the other will have larger weights. We verify the effectiveness of the proposed strategy on small networks by exhaustively examining all latent edges and demonstrate that performance is optimal or near-optimal. For large synthetic and real-world networks, the proposed method always outperforms other static strategies such as connecting nodes with the highest degree or eigenvector centrality.

Tables

Table 1. Table 1: Basic statistics of six real-world networks: number of nodes ( N 𝑁 N ), number of edges ( M 𝑀 M ), maximal degree ( k max subscript 𝑘 max k_{\mathrm{max}} ), first ( ⟨ k ⟩ delimited-⟨⟩ 𝑘 \langle k\rangle ) and second ( ⟨ k 2 ⟩ delimited-⟨⟩ superscript 𝑘 2 \langle k^{2}\rangle ) moments of the degree distribution, and the theoretical spreading threshold predicted by DMC λ c ∗ = 1 / ω 1 superscript subscript 𝜆 𝑐 1 subscript 𝜔 1 \lambda_{c}^{*}=1/\omega_{1} .

Networks	$N$	$M$	$k_{\max }$	$⟨k⟩$	$⟨k^2⟩$	${\lambda }_c^{\ast }$
Advogato	5042	39227	803	15.56	1284.00	0.014
Facebook	2888	2981	769	2.06	528.13	0.036
OpenFlights	2905	15645	242	10.77	601.45	0.016
Air traffic control	1226	2408	34	3.928	28.90	0.109
Adolescent health	2539	10455	27	8.24	86.41	0.076
Physicians	117	465	26	7.95	79.16	0.099

Equations

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

$$G=G^0+\delta G,$$

$$\delta G=\left({\begin{array}[]{cc}0&G_{ab}\\ G_{ba}&0\end{array}}\right)$$

$$q_i(t)=\prod _{j=1}^N[1-\lambda G_{ij}^0{p}_j(t)].$$

$$p_i(t+1)=(1-\gamma )p_i(t)+(1-p_i(t))\left(1-q_i(t)\right).$$

$$p^{\ast }=(1-\gamma )p^{\ast }+(1-p^{\ast })\circ (1-q^{\ast })$$

$$q_i^{\ast }=\prod _{j=1}^N(1-\lambda G_{ij}^0{p}_j^{\ast }),$$

$$\mathcal{P}=N^{-1}{1}^{\mathrm{T}}{p}^{\ast }.$$

$$\begin{array}{rl}{p}^{\ast }+\delta p(t+1)& =(1-\gamma )(p^{\ast }+\delta p(t))\\ & +(1-p^{\ast }-\delta p(t))\circ (1-q^{\ast }-\delta q(t)).\end{array}$$

$$\delta p(t+1)=(q^{\ast }-\gamma )\delta p(t)-(1-p^{\ast })\cdot \delta q(t).$$

$$q_i^{\ast }+\delta q_i(t)=\prod _{j=1}^N\left(1-\lambda (G_{ij}^0+\delta G_{ij})(p_j^{\ast }+\delta p_j(t))\right).$$

$$\begin{array}{rl}\delta q(t)=& -\lambda q^{\ast }\circ \left(G^0+\delta G\right)Z\delta p(t)\\ & +q^{\ast }\circ \delta G\log (1-\lambda p^{\ast }),\end{array}$$

$$Z_{ij}={\delta }_{ij}\frac{1}{1-\lambda p_j^{\ast }}.$$

$$\begin{array}{rl}\delta p(t+1)=& (q^{\ast }-\gamma )\delta p(t)\\ & +(1-p^{\ast })\circ \lambda q^{\ast }\circ \left(G^0+\delta G\right)Z\delta p(t)\\ & -(1-p^{\ast })\circ q^{\ast }\circ \delta G\log (1-\lambda p^{\ast }).\end{array}$$

$$\delta p(t+1)=X\delta p(t)+y,$$

$$X=\lambda \mathrm{diag}(q^{\ast }-p^{\ast }\circ q^{\ast })(G^0+\delta G)Z+\mathrm{diag}(q^{\ast }-\gamma )$$

$$y=-(1-p^{\ast })\circ q^{\ast }\circ \delta G\log (1-\lambda p^{\ast }).$$

$$\delta p^{\ast }=X\delta p^{\ast }+y,$$

$$\delta p^{\ast }={\left(\mathbb{I}-X\right)}^{-1}y.$$

$$\delta \mathcal{P}:=N^{-1}{1}^{\mathrm{T}}\delta p=N^{-1}{1}^{\mathrm{T}}{\left(\mathbb{I}-X\right)}^{-1}y$$

$$S_{\mathrm{Katz}}=1^{\mathrm{T}}{\left(\mathbb{I}-\beta G\right)}^{-1},$$

$${\left(\mathbb{I}-\beta G\right)}^{-1}=\mathbb{I}+\beta G+{\beta }^2{G}^2+\cdots ,$$

$$S_{\mathrm{Dyn}}=1^{\mathrm{T}}{\left(\mathbb{I}-X\right)}^{-1}.$$

$$X_{ij}=G_{ij}\frac{\lambda \left(1-p_i^{\ast }\right)q_i^{\ast }}{1-\lambda p_j^{\ast }}$$

$${\left(\mathbb{I}-X\right)}^{-1}\approx {\gamma }^{-1}{\left(\mathbb{I}-\frac{\lambda }{\gamma }G\right)}^{-1}.$$

$$y_i=-(1-p_i^{\ast })q_i^{\ast }\sum _{j=1}^N\delta G_{ij}\log \left(1-\lambda p_j^{\ast }\right).$$

$${\left(p_a^{\ast }\right)}_i={\left(p^{\ast }\right)}_i$$

$$(Z_a{)}_{ik}={\delta }_{ik}\frac{1}{1-\lambda {\left(p_a^{\ast }\right)}_i}=Z_{ik},$$

$$X^0=\lambda \mathrm{diag}(q^{\ast }-p^{\ast }\circ q^{\ast })G^{0}Z+\mathrm{diag}(q^{\ast }-\gamma )$$

$$\delta X=\lambda \mathrm{diag}(q^{\ast }-p^{\ast }\circ q^{\ast })\delta GZ$$

Full text

I Introduction

Promoting spreading dynamics in networked systems is attracting considerable attention in network science, statistical physics, and computer science Lü et al. (2016a). Maximizing spreading prevalence is of both theoretical and practical importance for achieving better information spreading and providing vaccination guidance. Strategies for maximizing spreading prevalence can be roughly divided into three categories: identifying vital nodes Kitsak et al. (2010); Morone and Makse (2015); Lü et al. (2016b); Chen et al. (2009, 2010); Morone et al. (2016); Hu et al. (2018); Ji et al. (2017); Ren et al. (2014); Liao et al. (2017); Pei and Makse (2013); Chen et al. (2012); Liu et al. (2017), designing effective transmission strategies Yang et al. (2008a); Gao et al. (2016); Yang et al. (2008b); Roshani and Naimi (2012); Gao et al. (2017); Cui et al. (2018), and performing network structural perturbations Aguirre et al. (2013); Del Ferraro et al. (2018); Milanese et al. (2010); Van Mieghem et al. (2010). For vital node identification, centrality measures, such as K-core, H-index, betweenness, and degree centrality, are assigned to network nodes. Nodes with high centrality are then chosen to be initial seeds for spreading. For effective transmission, spreading protocols have been designed to avoid invalid contacts (i.e., contacts among infected nodes). For performing structural perturbations, networks are modified slightly to promote spreading Aguirre et al. (2013). Structural perturbations are also widely applied to enhance network synchronizability  Aguirre et al. (2014); Li et al. (2016); Wei et al. (2018a, b); Dai et al. (2019).

The effectiveness of strategies for promoting spreading relies on the underlying spreading models. Various models, such as the susceptible-infected-susceptible (SIS), susceptible-infected-recovered (SIR), and threshold models, have been employed to test the effectiveness of such strategies Lü et al. (2016a). These spreading models can be divided into two classes, namely, simple and complex contagions Centola (2018); Guilbeault et al. (2018). In simple contagions, a susceptible individual could be infected by a single contact with an infected individual. Simple contagions are usually applied to model disease spreading Pastor-Satorras and Vespignani (2001) and simple information spreading (e.g., hashtags Romero et al. (2011)). In complex contagions, individuals evaluate the legitimacy of the information and make a risk assessment; the probability of infection increases with the cumulated number of contacts with other infected social peers. This mechanism is called social reinforcement  Watts (2002); Centola and Macy (2007); Aral and Nicolaides (2017); Lü et al. (2011). Complex contagions are usually applied to model complex information spreading (e.g., political information Romero et al. (2011)) and behavior adoption Centola and Macy (2007); Aral and Nicolaides (2017); Unicomb et al. (2019); Karsai et al. (2014). More complex spreading mechanisms, such as the coevolution of multiple diseases and/or information, are discussed in the recent review Wang et al. (2019).

The spreading dynamics in multilayer networks can be fundamentally different from that in single-layer networks da Silva et al. (2018); De Domenico et al. (2016); Granell et al. (2013, 2014); Wang et al. (2014); Chen et al. (2018); Gao et al. (2012); Tejedor et al. (2018); Wang et al. (2018). For instance, Granell et al. Granell et al. (2013) demonstrated that epidemic spreading has a metacritical point defined by the awareness dynamics and the topology of multilayer networks. The structure of the interconnections between two networks significantly affects robustness Radicchi and Arenas (2013); Reis et al. (2014); Van Mieghem (2016); Cozzo et al. (2019), synchronization Aguirre et al. (2014); Zhang et al. (2015) and spreading dynamics Hu et al. (2014); de Arruda et al. (2017); Sanz et al. (2014). Saumell-Mendiola et al. Saumell-Mendiola et al. (2012) demonstrated that interlayer degree correlations might trigger epidemic outbreaks. Wang et al. Wang et al. (2014) considered the coevolution of epidemics and information spreading in multilayer networks, and it was demonstrated that the interlayer degree correlations can also suppress epidemic outbreaks without altering the outbreak threshold.

A natural question is to determine the optimal interlayer structure for spreading in multilayer networks. To address this, Aguirre et al. Aguirre et al. (2013) applied a matrix perturbation approach and demonstrated that adding a connection between two nodes with large eigenvector centrality is more likely to promote the spreading dynamics for two competing networks. Recently, Pan et al. Pan et al. (2019) suggested applying perturbation theory to the adjacency matrix to obtain the optimal interconnections between two networks. This method is effective near the spreading threshold when a small number of edges are added.

In this study, we consider the problem of choosing a single interlayer edge that maximizes the spreading prevalence of the SIS model in two-layer networks. The SIS model can be applied to simple information or disease spreading. Therefore, understanding the maximization of the spreading prevalence may facilitate the promotion of information spreading or provide vaccination guidance Wang et al. (2016a). We develop a theoretical framework that provides the optimal or near-optimal interconnecting edge for all parameter regions. Starting with the discrete Markovian chain approach for the SIS model in two isolated networks, we propose a perturbation method so that the spreading prevalence in the interconnected network may be accurately approximated. The edge with the largest incremental spreading prevalence is then chosen as the interconnecting edge. The incremental spreading prevalence incorporates information regarding both network structure and spreading dynamics. Moreover, it has a simple physical interpretation as a variant of Katz centrality Newman (2010), where the adjacency matrix is replaced by a matrix with weights depending on the dynamical information of the spreading process.

The paper is organized as follows. We present the model in Sec. II and then develop a theory for obtaining the optimal interconnecting strategy in Sec. III. In Sec. IV, we perform extensive numerical simulations to verify the effectiveness of the proposed strategy. Sec. V concludes the paper.

II Model description

We consider the SIS model in two-layer networks. Let $a$ and $b$ be the two layers respectively. The number of nodes in $a$ and $b$ is denoted by $N_{a}$ and $N_{b}$, respectively, and the number of edges by $M_{a}$ and $M_{b}$, respectively. The adjacency matrices of the two layers are $G_{a}$ and $G_{b}$, and we assume that there are no interconnecting edges between them. Let $N={N_{a}+N_{b}}$. Then, the adjacency matrix $G^{0}$ of the two isolated layers combined is the following $N\times N$ matrix:

[TABLE]

We note that [math] in the off-diagonal part denotes zero matrices. There are multiple ways to interconnect the two isolated networks, and the dynamics of the interconnected network rely on the interlayer structure. Our aim is to determine the optimal interconnecting edge such that the spreading prevalence is maximized.

By adding the interconnecting edge, the adjacency matrix becomes

[TABLE]

where

[TABLE]

is the adjacency matrix for the interconnection between the two isolated networks. When $(G_{ab})_{ij}=(G_{ba})_{ji}=1$ for $i\in\{1,\cdots,N_{a}\}$ and $j\in\{1,\cdots,N_{b}\}$, an undirected edge is added between nodes $i$ and $j$.

We adopt the classical SIS model as the spreading model. Thus, each node can be in either the susceptible or infected state. Initially, a small fraction of nodes are selected as infected seeds, and the remaining nodes are susceptible. At each time step, every infected node in $a$ ($b$) tries to infect the susceptible neighbors in the same network with probability $\lambda_{a}$ ($\lambda_{b}$) and infect susceptible neighbors in $b$ ($a$) with probability $\lambda_{ab}$ ($\lambda_{ba}$). Then, all the infected nodes return to the susceptible state with probability $\gamma_{a}\ (\gamma_{b})$. We assume $\lambda_{a}=\lambda_{b}=\lambda_{ab}=\lambda_{ba}=\lambda$ and $\gamma_{a}=\gamma_{b}$. In the limit, the system reaches the steady state, and the fraction of infected nodes fluctuates around a stable value. Our aim is to choose an interconnecting edge such that the infected density of the new steady state in the interconnected network is maximized.

III Theoretical analysis

To study the SIS model in networks, we adopt the discrete Markovian chain (DMC) approach Gómez et al. (2010), which assumes that there are no dynamical correlations among the states of neighbors Wang et al. (2016b). In this section, we first present the DMC approach for the SIS model in the network $G^{0}$ when there are no interconnections between the networks $a$ and $b$. Then, using a perturbation method for DMC, we derive a formula that approximately provides the spreading prevalence in the interconnected network. Subsequently, we discuss physical interpretations of this formula, and finally, we study the problem of determining the optimal interconnecting edge based on the obtained formula.

III.1 Perturbation method for the discrete Markovian chain

Let $p_{i}(t)$ be the probability that node $i$ is infected at time $t$. Then, the node is susceptible with probability $1-p_{i}(t)$. If $i$ is in infected state at $t+1$, then either it was infected at $t$ and has not recovered, or it was susceptible at $t$ and has been infected by at least one infected neighbor. The former case occurs with probability $(1-\gamma)p_{i}(t)$ and the latter with probability $(1-p_{i}(t))\left(1-q_{i}(t)\right)$. Here, $1-q_{i}(t)$ is the probability that node $i$ is infected by at least one infected neighbor at time $t$, which is given by

[TABLE]

Combining the two cases, the evolution equation of $p_{i}(t)$ can be written as

[TABLE]

In the steady state, we have $p_{i}(t)=p_{i}(t+1)=p_{i}^{*}$ and $q_{i}(t)=q_{i}(t+1)=q_{i}^{*}$. Writing Eqs. (4) and (5) in terms of vectors in the steady state yields

[TABLE]

and

[TABLE]

where $p^{*}$, $q^{*}$ are vectors of length $N$ with entries $p^{*}=(p_{1}^{*},\cdots,p_{N}^{*})^{\mathrm{T}}$, $q^{*}=(q_{1}^{*},\cdots,q_{N}^{*})^{\mathrm{T}}$, and $\circ$ denotes component-wise vector product. The expected number of infected nodes in the steady state is

[TABLE]

Previous studies Gómez et al. (2010); de Arruda et al. (2017) demonstrated that a globally spreading outbreak occurs when the effective transmission probability $\lambda^{*}=\lambda/\gamma$ is larger than $1/\omega_{1}$, where $\omega_{1}$ is the leading eigenvalue of adjacency matrix $G^{0}$. That is, the spreading outbreak threshold is $\lambda^{*}_{c}=1/\omega_{1}$, whereas if $\lambda^{*}_{c}\leq 1/\omega_{1}$, then no outbreaks will be observed.

We now add one interconnecting edge between the two networks. Clearly, the spreading prevalence will increase after the edge is added. Subsequently, we develop a perturbation method to obtain an approximate estimate of the incremental spreading prevalence in the interconnected network.

When an interconnection is added between the two isolated networks, the adjacency matrix becomes $G=G^{0}+\delta G$. The fixed point of $p(t)$ in the interconnected network deviates from $p^{*}$ and the magnitude of deviation depends on where the interconnection is added. Nevertheless, as long as the two isolated networks are large enough, the modification in network structure can be regarded as small. As a consequence, the fixed point of $p(t)$ in the interconnected network should stay close to $p^{*}$. Since we focus on the case of adding one interconnection, this assumption should be valid for moderate network size. The actual magnitude of incremental spreading prevalence by interconnecting the networks can be seen from numerical results in Sec. IV, for example, in Fig. 1.

We now iterate the DMC equations in the interconnected network with initial condition $p(0)=p^{*}$, and then we use the decompositions $p(t)=p^{*}+\delta p(t)$ and $q(t)=q^{*}+\delta q(t)$ for some small $\delta p(t)$ and $\delta q(t)$. More explicitly, Eq. (5) in the interconnected network becomes

[TABLE]

Expanding Eq. (9) and substituting Eq. (6) yields

[TABLE]

We note that as $\delta p(t)$ and $\delta q(t)$ are assumed small, the second-order term $\delta p(t)\circ\delta q(t)$ is ignored. Similarly, Eq. (4) (the iteration equation for $q(t)$) in the interconnected network becomes

[TABLE]

As before, by expanding this equation up to first-order terms in $\delta p(t)$, we obtain

[TABLE]

where $\log(1-\lambda p^{*})$ is the vector obtained by taking the logarithm in each entry of $1-\lambda p^{*}$, and $Z$ is the $N\times N$ diagonal matrix with entries

[TABLE]

The detailed derivation of Eq. (12) is provided in Appendix A.

Substituting Eq. (12) back into Eq. (10) yields the following iteration formula for $\delta p(t)$:

[TABLE]

This equation can be written in terms of matrix multiplication as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

Here, $\mathrm{diag}\left(\cdot\right)$ denotes the diagonal matrix with the elements of the input vector as diagonal entries. The stationary solution $\delta p^{*}$ of the perturbed system satisfies

[TABLE]

or in the closed form

[TABLE]

This provides an explicit relation between the interconnection edge and the stationary infected density increment. Therefore, it remains to choose $\delta G$ such that the incremental spreading prevalence

[TABLE]

is maximized. We note that Eq. (20) holds even when we add multiple interconnecting edges.

III.2 Physical interpretations

Before analytically studying the optimization of Eq. (20), we should intuitively understand which interconnecting edge will give larger incremental spreading prevalence. We recall the Katz centrality Newman (2010) $S_{\mathrm{Katz}}$, which is defined by

[TABLE]

where $G$ is the adjacency matrix and $\beta$ a tunable parameter. Then, $S_{\mathrm{Katz}}$ is a vector with entries representing the centrality of the corresponding nodes. Katz centrality is defined by considering the number of weighted walks between nodes, where $\beta$ is the attenuation factor of walk length Newman (2010). The matrix inverse in $S_{\mathrm{Katz}}$ has expansion

[TABLE]

where $G^{t}$ is the matrix multiplication of $G$ by itself $t$ times. The entry $G^{t}_{ij}$ of $G^{t}$ then counts the number of walks of length $t$ between nodes $i$ and $j$.

We define a row vector

[TABLE]

Then, Eq. (20) can be written as $\delta\mathcal{P}=N^{-1}S_{\mathrm{Dyn}}y$. The vector $S_{\mathrm{Dyn}}$ has the same form as $S_{\mathrm{Katz}}$, with $\beta G$ in $S_{\mathrm{Katz}}$ replaced by $X$. We now further explore this connection and interpret $S_{\mathrm{Dyn}}$ as a weighted version of $S_{\mathrm{Katz}}$.

By the definition in Eq. (16), the entries of $X$ are given by

[TABLE]

for $i\neq j$. We note that $X_{ij}$ is nonzero only when $G_{ij}$ is nonzero; thus, $X$ can be understood as a weighted network with edge weights defined in terms of the dynamical information provided by $p^{*}$ and $q^{*}$. By Eq. (24), it is straightforward that the edge weight $X_{ij}$ is a decreasing function of $p^{*}_{i}$ and an increasing function of $q^{*}_{i},p^{*}_{j}$. That is, the edge connecting $i$ and $j$ will have a larger weight if it is less infected at $i$ (small $q^{*}_{i}$) and a neighborhood of $i$ (large $q^{*}_{i}$) but highly infected at $j$ (large $p^{*}_{j}$). More briefly, edges connecting less infected and highly infected regions will have larger weights.

Accordingly, $S_{\mathrm{Dyn}}$ is simply the Katz index defined on the weighted graph. $S_{\mathrm{Katz}}$ does not distinguish walks with the same length, as can be seen from Eq. (22), whereas $S_{\mathrm{Dyn}}$ further weights walks using dynamical information. When the transmission probability is below the spreading threshold, we have $p^{*}\approx 0$ and $q^{*}\approx 1$. Then, $X\approx 1-\gamma+\lambda G$ and

[TABLE]

In this case, $S_{\mathrm{Dyn}}$ reduces to $S_{\mathrm{Katz}}$ with $\beta=\lambda/\gamma$ (up to a constant factor $\gamma^{-1}$).

The incremental spreading prevalence $\delta\mathcal{P}=N^{-1}S_{\mathrm{Dyn}}y$ is then the weighted average of $S_{\mathrm{Dyn}}$, with nodes again weighted by the vector $y$. By the definition in Eq. (17), the entries of $y$ are

[TABLE]

Similarly, $y_{i}$ is a decreasing function of $p^{*}_{i}$ and an increasing function of $q^{*}_{i}$ and $p^{*}_{j}$. Thus, $y_{i}$ takes lager values if $i$ and its neighborhood are less infected and is connected to a highly infected node by an interconnection in $\delta G$.

By combining the discussions on $X$ and $y$, the optimal strategy can be understood as selecting an edge such that the infection is more easily transmitted from highly infected to less infected regions, which is consistent with intuition. We will refer to the method proposed in this section as dynamical Katz method.

III.3 Choosing the optimal edge

We now discuss the optimization of Eq. (20). We will consider only one connecting edge between the two isolated networks, that is, the optimal edge. We first introduce some notations. For the vector $p^{*}$, let $p^{*}_{a}$ be its part corresponding to network $a$. Specifically, $p^{*}_{a}$ is a vector of length $N_{a}$ with elements

[TABLE]

for $1\leq i\leq N_{a}$. $p^{*}_{b}$, $q^{*}_{a}$, $q^{*}_{b}$, $y_{a}$, and $y_{b}$ are defined analogously. We define the $N_{a}\times N_{a}$ diagonal matrix $Z_{a}$ with entries

[TABLE]

for $1\leq i,k\leq N_{a}$. $Z_{b}$ is defined analogously and corresponds to $b$.

We decompose $X$ as $X=X^{0}+\delta X$, where

[TABLE]

depends only on $G^{0}$, and

[TABLE]

depends only on $\delta G$. We note that $X^{0}$ is a diagonal block matrix and can be further written as

[TABLE]

where $X^{0}_{a}$ is the block diagonal part of $X^{0}$ that depends only on $G_{a}$, with

[TABLE]

Similarly,

[TABLE]

$\delta X$ is an off-diagonal block matrix

[TABLE]

with the off-diagonal blocks given by

[TABLE]

Using the properties of block matrices, the matrix inverse in Eq. (19) can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

We now assume that we add an interconnecting edge between node $i$ of network $a$ and node $j$ of $b$. By the Sherman–Morrison formula, the resulting increment in the spreading prevalence can be written in the following explicit form:

[TABLE]

where

[TABLE]

and

[TABLE]

The detailed derivation of Eq. (39) is given in Appendix B.

This provides a simple formula for the spreading prevalence in the interconnected network. The optimal strategy is simply to select the edge with the highest corresponding $\delta\mathcal{P}$. This strategy relies not only on the network topology (i.e., the adjacency matrices $G_{a}$ and $G_{b}$) but also on the dynamical information of the spreading process when the two networks are isolated (i.e., $\lambda$, $\gamma$, $q^{*}$, and $p^{*}$).

IV Numerical simulations

In this section, we perform extensive numerical simulations on both synthetic and real-world networks to verify the performance of the strategy. We note that we do not compare the DMC predictions with Monte Carlo simulations because the DMC approach can accurately predict the simulations Gómez et al. (2010). In the following, the numerical value of $\delta\mathcal{P}$ obtained by iterating the DMC is denoted by $\delta\mathcal{P}^{\mathrm{num}}$. The optimal edge predicted by the DMC equations is called the numerical optimal edge. The approximate $\delta\mathcal{P}$ predicted by Eq. (39) is denoted by $\delta\mathcal{P}^{\mathrm{approx}}$.

For two networks with number of nodes $N_{a}$ and $N_{b}$, there are in total $M_{l}=N_{a}\times N_{b}$ latent interconnections. For small networks, it is possible to check all the latent connections exhaustively so that the optimal may be determined. However, for large $N$, an exhaustive search is slow and gradually becomes impossible. We first use small networks to verify the accuracy of $\delta\mathcal{P}^{\mathrm{approx}}$ predicted by Eq. (39) and compare the optimal edge by this strategy with the numerical optimal edge (the optimal edge predicted by the DMC equations).

To construct synthetic networks, we adopt the uncorrelated configuration model with power-law degree distributions. Specifically, we set the degree distributions of networks $a$ and $b$ to $P(k)\sim k^{-\alpha_{a}}$ and $P(k)\sim k^{-\alpha_{b}}$ respectively, where $\alpha_{a}$ and $\alpha_{b}$ are the degree exponents. The network sizes are set to $N_{a}=N_{b}=100$. Without loss of generality, we set the recovery probability of the SIS model to $\gamma=0.5$ and make the infection probability $\lambda$ a tunable parameter.

We first compare $\delta\mathcal{P}^{\mathrm{approx}}$ predicted by Eq. (39) with the predictions by the DMC approach. For each latent edge connecting node $i\in\{1,\cdots N_{a}\}$ and node $j\in\{1,\cdots,N_{b}\}$, we compute $\delta\mathcal{P}^{\mathrm{approx}}$ using Eq. (39) for $\lambda=0.3$ (Fig. 1(a)) and $\lambda=0.5$ (Fig. 1(c)). Then, we add the edge to the network and iterate the DMC to obtain $\delta\mathcal{P}^{\mathrm{num}}$, which is shown in Figs. 1(b) and  1(d) for $\lambda=0.3$ and $\lambda=0.5$, respectively. Nodes are arranged in identical order in Figs. 1(a), 1(b), also in identical order in Figs. 1(c), 1(d). The approximate values usually appear higher than the numerical values, but intuitively, they are strongly correlated. The maximum relative error $\left(\delta\mathcal{P}^{\mathrm{approx}}-\delta\mathcal{P}^{\mathrm{num}}\right)/\delta\mathcal{P}^{\mathrm{num}}$ for all edges is $0.315$ in Figs. 1(a) and (b), and $0.396$ in Figs. 1(c) and (d). However, we will demonstrate that they are almost linearly correlated in order, which suggests the approximate value is sufficient to obtain the optimal edge.

To see the correlations, we compute the Spearman’s rank correlation coefficient Lee et al. (2012); Wang et al. (2014) between the approximate and numerical $\delta\mathcal{P}$. We score all the latent interconnecting edges by their $\delta\mathcal{P}^{\mathrm{num}}$ and $\delta\mathcal{P}^{\mathrm{approx}}$; then, two rankings can be obtained. Let $r_{ij}$ and $r^{\prime}_{ij}$ be the rank of the edge connecting node $i$ in network $a$ and node $j$ in network $b$ scored by $\delta\mathcal{P}^{\mathrm{num}}$ and $\delta\mathcal{P}^{\mathrm{approx}}$, respectively. Spearman’s rank correlation coefficient is defined as

[TABLE]

We plot $m_{s}$ as a function of $\lambda$ in Fig. 2(a). It can be observed that Spearman’s rank correlation coefficients are close to $1$ for all $\lambda$. The minimum value of $m_{s}$ for all $\lambda$ in Fig. 2(a) is $0.9968$. This suggests that the proposed strategy accurately predicts the overall order of $\delta P^{\mathrm{num}}$.

In addition to the strong overall correlations for the approximate and numerical values of $\delta\mathcal{P}$, we are particularly concerned with the top-ranked edge. We further verify the performance of the strategy by comparing the predicted optimal edge with the numerical optimal edge. For each $\lambda$, we select the edge with the highest $\delta\mathcal{P}^{\mathrm{approx}}$ predicted by the dynamical Katz method and compute its numerical rank in all the latent edges. The edge rank versus $\lambda$ is shown in Fig. 2(b). It can be seen that the rank is $1$ or near $1$ for all values of $\lambda$. When the rank is exactly $1$, the optimal edge predicted by the dynamical Katz method coincides with the numerical optimal edge, and this is the case for most values of $\lambda$.

As the dynamical Katz strategy incorporates information regarding both the network structure and spreading dynamics, it is useful to compare it with simple strategies that consider only the static network structure to understand the role of dynamical information. Specifically, we consider the strategy of connecting the two nodes with the highest degree or eigenvector centrality. The normalized ranks (ranks divided by $M_{l}=N_{a}\times N_{b}$) by the dynamical Katz and the two static strategies are shown in Fig. 2(c). All three strategies are optimal or near-optimal when the transmission probability $\lambda$ is slightly above the critical value, but the two static strategies fail quickly when $\lambda$ becomes large, whereas the dynamical Katz method still performs well.

As discussed in Sec. III.2, when $p^{*}\approx 0$, the dynamical Katz matrix reduces to the Katz matrix. When $\lambda/\gamma$ is small, we have

[TABLE]

and this reduces to degree centrality. For uncorrelated configuration models, degree and eigenvector centrality are strongly correlated. When $\lambda$ is small, nodes with high centrality values (i.e., degree and eigenvector centrality) have a larger probability to be infected. If we connect them with an edge, then high-centrality nodes together with their neighbors could form an infected cluster Goltsev et al. (2012) and further transmit the infection to other nodes. Thus, for small $\lambda$, the degree and eigenvector strategies perform well. For large values of $\lambda$, globally spreading outbreaks occur, and nodes with small centrality have a higher probability to be susceptible. In this case, additional connections to these nodes are required for promoting the spreading dynamics. Therefore, both the degree and eigenvector strategies fail, and the dynamical information should be considered.

For large networks, exhaustive searching becomes impossible. In this case, we compare the performance of the dynamical Katz method with that of the two static methods based on degree and eigenvector centrality. For the three methods, we add the predicted optimal edge separately and compare the resulting $\delta\mathcal{P}$. We first consider synthetic networks. We construct three pairs of networks with power-law degree distributions, with degree exponents (i) $\alpha_{a}=2.3$, $\alpha_{b}=3.0$, (ii) $\alpha_{a}=3.0$, $\alpha_{b}=3.0$, and (iii) $\alpha_{a}=4.0$, $\alpha_{b}=3.0$. The graphs of $\delta\mathcal{P}$ versus $\lambda$ for the three network pairs are shown in Fig. 3. We further add the semi-log plot in the insets for the first two network pairs for better visualization, with $\delta\mathcal{P}$ on a logarithm scale. When $\lambda$ is close to the critical point, all three strategies exhibit highly similar performance. As in small networks, this could also be near the maximal possible value of $\delta\mathcal{P}$. When $\lambda$ becomes large, dynamical Katz outperforms the other two static methods for all network pairs. In this case, connecting nodes with large degree or eigenvector centrality yields almost zero $\delta\mathcal{P}$, which decays with $\lambda$, as can be seen in, e.g., the insets in Figs. 3(a) and (b). Moreover, it is worth noticing that $\delta\mathcal{P}$ is always maximized slightly above the spreading threshold, which suggests that the marginal improvement is maximized near the critical point.

We now test the dynamical Katz method on real-world networks. Three pairs of networks are considered: (i) Advogato Massa et al. (2009), Facebook Leskovec and Mcauley (2012), (ii) OpenFlights Kon , Air traffic control Kon , and (iii) Adolescent health Moody (2001), Physicians Coleman et al. (1957). The first pair (Advogato and Facebook) are two online social networks, the second pair (OpenFlights and Air traffic control) are infrastructure networks of airports and flights, and the third pair (Adolescent health and Physicians) are two offline social networks. The networks were downloaded from Kon , and details can be found therein. Some basic statistics are shown in TABLE 1.

$\delta\mathcal{P}$ versus $\lambda$ for the three network pairs are shown in Fig. 4. As in the case of the synthetic networks, it can be seen that the dynamical Katz method performs best for all values of $\lambda$. However, for small values of $\lambda$, the three methods are quite close. For larger $\lambda$, the two static method yield $\delta\mathcal{P}$ very close to zero, whereas the dynamical Katz exhibits significant improvement. The results further confirmed the effectiveness of the dynamical Katz method on real-world networks.

For both synthetic and real-world networks, the Degree and Eigen methods work well near the critical $\lambda$ and fail when $\lambda$ becomes large, while the dynamical Katz method performs well in all the parameter region. To better understand the structural properties of the optimal edge predicted by the dynamical Katz method, we study how the degrees and eigenvector centralities of the optimal edge’s two end-nodes change with $\lambda$. Let $k_{a}^{\mathrm{opt}}$ and $k_{b}^{\mathrm{opt}}$ be the degrees of the optimal edge’s two end-nodes in layer $a$ and layer $b$ respectively. Similarly, we define $v_{a}^{\mathrm{opt}}$ and $v_{b}^{\mathrm{opt}}$ as the two nodes’ eigenvector centralities. Let $k_{a}^{\mathrm{max}}$ and $k_{b}^{\mathrm{max}}$ be the maximum degree of $a$ and $b$ respectively, while $v_{a}^{\mathrm{max}}$ and $v_{b}^{\mathrm{max}}$ be the maximum eigenvector centrality of $a$ and $b$ respectively. $k_{a}^{\mathrm{opt}}/k_{a}^{\mathrm{max}}$, $k_{b}^{\mathrm{opt}}/k_{b}^{\mathrm{max}}$, $v_{a}^{\mathrm{opt}}/v_{a}^{\mathrm{max}}$ and $v_{b}^{\mathrm{opt}}/v_{b}^{\mathrm{max}}$ versus $\lambda$ are shown in Figs. 5(a)-(d) respectively. When $\lambda$ is near the critical point, nodes with high degree and eigenvector centrality are chosen to be connected. Near the critical point, the spreading prevalence is small, connecting nodes with high centrality will help to maintain the infected cluster and further transmit the infection to other nodes. When $\lambda$ becomes large, nodes with high centrality have a very high probability to be infected, therefore connecting these nodes becomes unnecessary. As shown in Fig. 5, the degrees and eigenvector centralities of the optimal edge become small when $\lambda$ becomes large. The numerical results have further verified the discussions about the relations between dynamical Katz and other two static methods (below Eq. (43)).

V Discussion

We studied the problem of determining the optimal interconnecting edge for promoting spreading dynamics. By applying a perturbation method to the DMC equations, we obtained a Katz-like index for predicting the spreading prevalence in the interconnected networks. This index accurately predicts the optimal interconnecting edge for promoting spreading over all parameter regions, as demonstrated in small networks. For large synthetic and real-world networks, the method outperformed certain static strategies, namely, connecting nodes with highest degree or eigenvector centrality. For small $\lambda$, the three strategies had similar performance. For large $\lambda$, the two heuristic strategies yielded almost zero incremental spreading prevalence, whereas the dynamical Katz method performed well. In addition to accurately predicting the optimal edge, the dynamical Katz method provides a clear physical interpretation of how the optimal edge is chosen.

We considered the addition of only one interconnecting edge, but real-world multilayer networks usually have multiple interconnecting edges. We note that Eq. (20), which estimates the incremental spreading prevalence in terms of interconnections, is valid for general interconnecting structures. This could provide the foundation for further study in the case of multiple edges. For the single-edge case, the interconnection matrix $C$ can be written as the outer product of two vectors. By applying the Sherman–Morrison formula, $\delta\mathcal{P}$ can take a simple form that is easy to optimize. When multiple edges are added, the outer product decomposition of $C$ cannot be used in general.

A simple heuristic method for adding multiple edges is to add edges one by one using the proposed method. Specifically, at each step, one edge is added using the dynamical Katz method, and then the DMC equations are iterated to converge in this new network. The procedure is repeated until all edges are added. By adding one edge, the dynamical Katz method is likely to be optimal or near-optimal; therefore, this heuristic algorithm can be considered a greedy algorithm. The performance of such a greedy algorithm can be further analyzed, and more sophisticated algorithms could be designed. We leave this as an open issue for future exploration. Moreover, the perturbation method developed in this study could also be extended to other types of networks (e.g., temporal networks) and spreading models (e.g., social contagions and cascading failures).

Acknowledgements.

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

Appendix A: Derivation of the perturbuted equation for $\delta q(t)$

In this section we detail the derivation of the perturbed equation for $\delta q(t)$ in Eq. (12) starting from Eq. (11). Since $G^{0}$ is a diagonal block matrix, while $\delta G$ is an off-diagonal block matrix, then $G^{0}_{ij}=1$ and $\delta G_{ij}=1$ can not be observed simultaneously. Thus the following equation holds

[TABLE]

which can be checked by substituting all possible combinations of $G^{0}_{ij}$ and $\delta G_{ij}$. Divide by $q^{*}_{i}$ for both sides and substitute Eq. (7) gives

[TABLE]

Note that the following relation holds

[TABLE]

since $G^{0}_{ij}\in\{0,1\}$ and similarly when replacing $G^{0}_{ij}$ by $\delta G_{ij}\in\{0,1\}$. Take the logarithm on both sides of Eq. (A.1), expand to the first orders of $\delta p_{i}(t)$, $\delta q_{i}(t)$, and apply the above relation gives

[TABLE]

Again the terms in the last summation can be checked satisfying

[TABLE]

With the above calculations, Eq. (A.3) can be written in matrix form as Eq. (12). This completes the derivation of Eq. (12).

Appendix B: Derivation of the incremental spreading prevalence by adding an edge

In this section we give the detailed derivation of Eq. (39), i.e., the explicit formula for incremental spreading prevalence when only adding one interconnecting edge.

Since we only add one edge, then the matrix $G_{ab}=G_{ba}^{\mathrm{T}}$ can be written as an outer product

[TABLE]

where $u$ is a vector of length $N_{a}$ with $u_{k}=\delta_{k,i}$ for $1\leq k\leq N_{a}$, and $v$ a length $N_{b}$ vector with $v_{k}=\delta_{k,j}$ for $1\leq k\leq N_{b}$. Recall $x_{ij}$ and $x_{ji}$ defined in Eq. (40), then it’s easy to check that

[TABLE]

Thus we have

[TABLE]

In other words, $\delta X_{ab}B\delta X_{ba}$ is an zero matrix expect in the $j$th element in the diagonal. The Sherman-Morrison formula says that

[TABLE]

With this formula we can construct $\left(\mathbb{I}-X\right)^{-1}$ easily from $C$. Similarly,

[TABLE]

Again recall the definitions of $c_{ij}$ and $c_{ij}$ in Eq. (41), then $y_{a},y_{b}$ can written as

[TABLE]

Combine the above computations, we arrive at the formula

[TABLE]

The first term on the r.h.s. of Eq. (B.7) can be written as

[TABLE]

where the first line is by substituting Eq. (B.4), and second line is by using definition of $y_{a}$ and $u$. For the second term of r.h.s. in Eq. (B.7),

[TABLE]

With similar computations, we can can obtain the expressions for the rest two terms in the r.h.s. of Eq. (B.7), which are

[TABLE]

and

[TABLE]

Combine the above computations we have

[TABLE]

which is Eq. (39).