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

Liming Pan  Web Sciences Center, School of Computer Science and Engineering,
University of Electronic Science and Technology of China, Chengdu, 611731,
China  Wei Wang<EMAIL_ADDRESS>Cybersecurity Research Institute, Sichuan
University, Chengdu 610065, China  Web Sciences Center, School of Computer
Science and Engineering, University of Electronic Science and Technology of
China, Chengdu, 611731, China  Shimin Cai  Web Sciences Center, School of
Computer Science and Engineering, University of Electronic Science and
Technology of China, Chengdu, 611731, China  Tao Zhou  Web Sciences Center,
School of Computer Science and Engineering, University of Electronic Science
and Technology of China, Chengdu, 611731, China

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

##  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:

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

We note that $0$ 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

$G=G^{0}+\delta G,$ (2)

where

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

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

$q_{i}(t)=\prod_{j=1}^{N}[1-\lambda G_{ij}^{0}p_{j}(t)].$ (4)

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

$p_{i}(t+1)=(1-\gamma)p_{i}(t)+(1-p_{i}(t))\left(1-q_{i}(t)\right).$ (5)

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

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

and

$q^{*}_{i}=\prod_{j=1}^{N}(1-\lambda G_{ij}^{0}p^{*}_{j}),$ (7)

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

$\mathcal{P}=N^{-1}\mathbf{1}^{\mathrm{T}}p^{*}.$ (8)

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

$\begin{split}p^{*}+\delta p(t+1)&=(1-\gamma)(p^{*}+\delta p(t))\\\
&+(1-p^{*}-\delta p(t))\circ(1-q^{*}-\delta q(t)).\end{split}$ (9)

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

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

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

$q_{i}^{*}+\delta q_{i}(t)=\prod_{j=1}^{N}\left(1-\lambda(G^{0}_{ij}+\delta
G_{ij})(p_{j}^{*}+\delta p_{j}(t))\right).$ (11)

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

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

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

$Z_{ij}=\delta_{ij}\frac{1}{1-\lambda p_{j}^{*}}.$ (13)

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)$:

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

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

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

where

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

and

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

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

$\delta p^{*}=X\delta p^{*}+y,$ (18)

or in the closed form

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

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

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

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

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

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

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

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

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

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

$X_{ij}=G_{ij}\frac{\lambda\left(1-p^{*}_{i}\right)q^{*}_{i}}{1-\lambda
p^{*}_{j}}$ (24)

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

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

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

$y_{i}=-(1-p^{*}_{i})q^{*}_{i}\sum_{j=1}^{N}\delta G_{ij}\log\left(1-\lambda
p^{*}_{j}\right).$ (26)

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

$\left(p^{*}_{a}\right)_{i}=\left(p^{*}\right)_{i}$ (27)

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

$(Z_{a})_{ik}=\delta_{ik}\frac{1}{1-\lambda\left(p^{*}_{a}\right)_{i}}=Z_{ik},$
(28)

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

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

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

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

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

$X^{0}=\left({\begin{array}[]{cc}X^{0}_{a}&0\\\
0&X^{0}_{b}\end{array}}\right),$ (31)

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

$X^{0}_{a}=\lambda\mathrm{diag}(q^{*}_{a}-p^{*}_{a}\circ
q^{*}_{a})G_{a}Z_{a}+\mathrm{diag}(q^{*}_{a}-\gamma).$ (32)

Similarly,

$X^{0}_{b}=\lambda\mathrm{diag}(q^{*}_{b}-p^{*}_{b}\circ
q^{*}_{b})G_{b}Z_{b}+\mathrm{diag}(q^{*}_{b}-\gamma).$ (33)

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

$\delta X=\left({\begin{array}[]{cc}0&\delta X_{ab}\\\ \delta
X_{ba}&0\end{array}}\right),$ (34)

with the off-diagonal blocks given by

$\begin{split}\delta X_{ab}&=\lambda\mathrm{diag}(q^{*}_{a}-p^{*}_{a}\circ
q^{*}_{a})G_{ab}Z_{b}\\\ \delta
X_{ba}&=\lambda\mathrm{diag}(q^{*}_{b}-p^{*}_{b}\circ
q^{*}_{b})G_{ba}Z_{a}.\end{split}$ (35)

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

$\begin{split}\left(\mathbb{I}-X\right)^{-1}=&\left({\begin{array}[]{cc}\mathbb{I}-X^{0}_{a}&-\delta
X_{ab}\\\ -\delta X_{ba}&\mathbb{I}-X^{0}_{b}\end{array}}\right)^{-1}\\\
=&\left({\begin{array}[]{cc}C&C\delta X_{ab}B\\\ D\delta
X_{ba}A&D\end{array}}\right),\end{split}$ (36)

where

$A=\left(\mathbb{I}-X^{0}_{a}\right)^{-1},\
B=\left(\mathbb{I}-X^{0}_{b}\right)^{-1}$ (37)

and

$\begin{split}C&=\left(\mathbb{I}-X^{0}_{a}-\delta X_{ab}B\delta
X_{ba}\right)^{-1},\\\ D&=\left(\mathbb{I}-X^{0}_{b}-\delta X_{ba}A\delta
X_{ab}\right)^{-1}.\end{split}$ (38)

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:

$\begin{split}N\delta\mathcal{P}=&\frac{c_{ij}+c_{ji}x_{ij}B_{jj}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}A\right)_{i}\\\
&+\frac{c_{ji}+c_{ij}x_{ji}A_{ii}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}B\right)_{j},\end{split}$
(39)

where

$\begin{split}x_{ij}&\mathop{:}=\lambda\left[\left(q^{*}_{a}\right)_{i}-\left(p^{*}_{a}\right)_{i}\left(q^{*}_{a}\right)_{i}\right]\left[1-\lambda\left(p^{*}_{b}\right)_{j}\right]^{-1},\\\
x_{ji}&\mathop{:}=\lambda\left[\left(q^{*}_{b}\right)_{j}-\left(p^{*}_{b}\right)_{j}\left(q^{*}_{b}\right)_{j}\right]\left[1-\lambda\left(p^{*}_{a}\right)_{i}\right]^{-1}\end{split}$
(40)

and

$\begin{split}c_{ij}&=-\left[\left(q^{*}_{a}\right)_{i}-\left(p^{*}_{a}\right)_{i}\left(q^{*}_{a}\right)_{i}\right]\log\left[1-\lambda\left(p^{*}_{b}\right)_{j}\right],\\\
c_{ji}&=-\left[\left(q^{*}_{b}\right)_{j}-\left(p^{*}_{b}\right)_{j}\left(q^{*}_{b}\right)_{j}\right]\log\left[1-\lambda\left(p^{*}_{a}\right)_{i}\right].\end{split}$
(41)

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.

Figure 1: (Color online) Incremental spreading prevalence $\delta\mathcal{P}$
by adding each latent interconnection separately. The vertical and horizontal
axes correspond to node IDs in the networks. Thus, each point in the plots
corresponds to an edge connecting $a$ and $b$, and its color represents the
value of $\delta\mathcal{P}$ by adding the edge. (a) Approximate and (b)
numerical predictions of $\delta\mathcal{P}$ with $\lambda=0.3$. (c)
Approximate and (d) numerical predictions of $\delta\mathcal{P}$ with
$\lambda=0.5$. The nodes are arranged in identical order for (a) and (b), as
well as for (c) and (d). Other parameters are set as $N_{a}=N_{b}=100$,
$\alpha_{a}=2.3$, $\alpha_{b}=3.0$ and $\gamma=0.5$.

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

$m_{s}=1-6\frac{\sum_{i=1}^{N_{a}}\sum_{j=1}^{N_{b}}\left(r_{ij}-r^{\prime}_{ij}\right)^{2}}{M_{l}(M_{l}^{2}-1)}.$
(42)

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}}$.

Figure 2: (Color online) Performance of different strategies versus
transmission probability. (a) Spearman’s rank correlation coefficient $m_{s}$
between ranks predicted by the dynamical Katz method and the numerical ranks.
(b) Numerical rank of the optimal edge predicted by the dynamical Katz method.
(c) Normalized numerical rank of the optimal edge predicted by dynamical Katz
(blue dotted line), degree (orange solid line), and eigenvector centrality
(yellow dashed line) methods. Below the spreading threshold, the prevalence
$\mathcal{P}$ is zero, and the rankings are trivial; thus, we consider
$\lambda$ only in the range starting slightly above the threshold. Other
parameters are set as $N_{A}=N_{B}=100$, $\alpha_{a}=3.0$, $\alpha_{b}=2.3$
and $\gamma=0.5$.

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

$\left(\mathbb{I}-X\right)^{-1}\approx\gamma^{-1}\mathbb{I}+\lambda G,$ (43)

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.

Figure 3: (Color online) Incremental spreading prevalence $\delta\mathcal{P}$
versus $\lambda$ in synthetic networks when one interconnecting edge is added.
Networks with power-law degree distributions are considered, where the degree
exponents are (a) $\alpha_{a}=2.3$, $\alpha_{b}=3.0$, (b) $\alpha_{a}=3.0$,
$\alpha_{b}=3.0$, and (c) $\alpha_{a}=4.0$, $\alpha_{b}=3.0$.
$\delta\mathcal{P}$ versus $\lambda$ in semi-log plots are shown in the
insets, where $\delta\mathcal{P}$ is on the logarithm scale. Data points
corresponding to the degree and eigenvector centrality strategies are highly
overlapped. We consider $\lambda$ in the range starting slightly above the
spreading threshold to avoid trivial cases. Other parameters are set as
$N_{A}=N_{B}=5000$ and $\gamma=0.5$.

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.

Figure 4: (Color online) Incremental spreading prevalence $\delta\mathcal{P}$ versus $\lambda$ when one interconnecting edge is added in real-world networks. The network pairs are (a) Advogato, Facebook , (b) OpenFlights, Air traffic control, and (c) Adolescent health, Physicians. Data points for the degree and the eigenvector centrality strategies are highly overlapped for some values of $\lambda$. $\delta\mathcal{P}$ versus $\lambda$ in semi-log plots are shown in the insets, where $\delta\mathcal{P}$ is on the logarithm scale. In (a), for large values of $\lambda$, $\delta\mathcal{P}$ given by Degree and Eigen are so close to zero such that numerically we have $\delta\mathcal{P}=0$ due to the limitation of numerical accuracy. Therefore in the inset of (a), data points for large $\lambda$ cannot be seen on a logarithm scale. Basic statistics of the six real-world networks can be found in TABLE 1. Table 1: Basic statistics of six real-world networks: number of nodes ($N$), number of edges ($M$), maximal degree ($k_{\mathrm{max}}$), first ($\langle k\rangle$) and second ($\langle k^{2}\rangle$) moments of the degree distribution, and the theoretical spreading threshold predicted by DMC $\lambda_{c}^{*}=1/\omega_{1}$. Networks | $N$ | $M$ | $k_{\mathrm{max}}$ | $\langle k\rangle$ | $\langle k^{2}\rangle$ | $\lambda^{*}_{c}$
---|---|---|---|---|---|---
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
Figure 5: (Color online) Normalized degrees and eigenvector centralities of
the optimal edge’s two end-nodes versus $\lambda$ for both synthetic and real-
world networks. (a) $k_{a}^{\mathrm{opt}}/k_{a}^{\mathrm{max}}$, i.e.,
normalized degree of the optimal edge’s end-node in layer $a$, versus
$\lambda$. (b) $k_{b}^{\mathrm{opt}}/k_{b}^{\mathrm{max}}$, i.e., normalized
degree of the optimal edge’s end-node in layer $b$, versus $\lambda$. (c)
$v_{a}^{\mathrm{opt}}/v_{a}^{\mathrm{max}}$, i.e., normalized eigenvector
centrality of the optimal edge’s end-node in layer $a$, versus $\lambda$. (d)
$v_{b}^{\mathrm{opt}}/v_{b}^{\mathrm{max}}$, i.e., normalized eigenvector
centrality of the optimal edge’s end-node in layer $b$, versus $\lambda$.
Different lines correspond to different network pairs, among which the solid
lines correspond to synthetic networks and dashed lines correspond to real-
world networks. The values of $\lambda$ considered start slightly above the
spreading threshold of each pair of 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

$\begin{split}&\left(1-\lambda(G^{0}_{ij}+\delta
G_{ij})\right)\left(p_{j}^{*}+\delta p_{j}(t)\right)\\\ =&\left(1-\lambda
G^{0}_{ij}(p_{j}^{*}+\delta p_{j}(t))\right)\left(1-\lambda\delta
G_{ij}(p_{j}^{*}+\delta p_{j}(t))\right),\end{split}$

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

$\begin{split}1+\frac{\delta
q_{i}(t)}{q_{i}^{*}}&=\prod_{j=1}^{N}\left(1-\frac{\lambda G^{0}_{ij}\delta
p_{j}(t)}{1-\lambda G_{ij}^{0}p_{j}^{*}}\right)\\\
&\times\prod_{j=1}^{N}\left(1-\frac{\lambda\delta G_{ij}\delta
p_{j}(t)}{1-\lambda\delta
G_{ij}p^{*}_{j}}\right)\prod_{j=1}^{N}\left(1-\lambda\delta
G_{ij}p^{*}_{j}\right).\end{split}$ (A.1)

Note that the following relation holds

$\frac{\lambda G^{0}_{ij}\delta p_{j}(t)}{1-\lambda
G_{ij}^{0}p_{j}^{*}}=G^{0}_{ij}\frac{\lambda\delta p_{j}(t)}{1-\lambda
p_{j}^{*}},$ (A.2)

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

$\begin{split}\frac{\delta
q_{i}(t)}{q_{i}^{*}}=&-\sum_{j=1}^{N}G^{0}_{ij}\frac{\lambda\delta
p_{j}(t)}{1-\lambda p_{j}^{*}}-\sum_{j=1}^{N}\delta G_{ij}\frac{\lambda\delta
p_{j}(t)}{1-\lambda p_{j}^{*}}\\\ &+\sum_{j=1}^{N}\log\left(1-\lambda\delta
G_{ij}p^{*}_{j}\right).\end{split}$ (A.3)

Again the terms in the last summation can be checked satisfying

$\log\left(1-\lambda\delta G_{ij}p^{*}_{j}\right)=\delta
G_{ij}\log\left(1-\lambda p^{*}_{j}\right).$ (A.4)

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

$G_{ab}=uv^{\mathrm{T}},$ (B.1)

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

$\delta X_{ab}=x_{ij}uv^{\mathrm{T}},\ \delta X_{ba}=x_{ji}vu^{\mathrm{T}}.$
(B.2)

Thus we have

$\delta X_{ab}B\delta X_{ba}=x_{ij}x_{ji}B_{jj}uu^{\mathrm{T}}$ (B.3)

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

$\begin{split}C&=\left(\mathbb{I}-X^{0}_{a}-x_{ij}x_{ji}B_{jj}uu^{\mathrm{T}}\right)^{-1}\\\
&=A+\frac{x_{ij}x_{ji}B_{jj}Auu^{\mathrm{T}}A}{1-x_{ij}x_{ji}A_{ii}B_{jj}}.\end{split}$
(B.4)

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

$\begin{split}D&=\left(\mathbb{I}-X^{0}_{b}-x_{ij}x_{ji}A_{ii}vv^{\mathrm{T}}\right)^{-1}\\\
&=B+\frac{x_{ij}x_{ji}A_{ii}Bvv^{\mathrm{T}}B}{1-x_{ij}x_{ji}A_{ii}B_{jj}}.\end{split}$
(B.5)

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

$y_{a}=c_{ij}u,\ y_{b}=c_{ji}v.$ (B.6)

Combine the above computations, we arrive at the formula

$\begin{split}N\delta\mathcal{P}=&\mathbf{1}^{\mathrm{T}}Cy_{a}+\mathbf{1}^{\mathrm{T}}C\delta
X_{ab}By_{b}+\\\
&+\mathbf{1}^{\mathrm{T}}Dy_{b}+\mathbf{1}^{\mathrm{T}}D\delta
X_{ba}Ay_{a}.\end{split}$ (B.7)

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

$\begin{split}\mathbf{1}^{\mathrm{T}}Cy_{a}&=\mathbf{1}^{\mathrm{T}}Ay_{a}+\frac{x_{ij}x_{ji}B_{jj}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}A\right)_{i}u^{\mathrm{T}}Ay_{a}\\\
&=\frac{c_{ij}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}A\right)_{i},\end{split}$
(B.8)

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),

$\begin{split}\mathbf{1}^{\mathrm{T}}C\delta
X_{ab}By_{b}=&c_{ji}x_{ij}B_{jj}\mathbf{1}^{\mathrm{T}}Cu\\\
=&\frac{c_{ji}x_{ij}B_{jj}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}A\right)_{i}.\end{split}$
(B.9)

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

$\mathbf{1}^{\mathrm{T}}Dy_{b}=\frac{c_{ji}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}B\right)_{j},$
(B.10)

and

$\mathbf{1}^{\mathrm{T}}D\delta
X_{ba}Ay_{a}=\frac{c_{ij}x_{ji}A_{ii}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}B\right)_{j}.$
(B.11)

Combine the above computations we have

$\begin{split}N\delta\mathcal{P}=&\frac{c_{ij}+c_{ji}x_{ij}B_{jj}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}A\right)_{i}+\\\
&+\frac{c_{ji}+c_{ij}x_{ji}A_{ii}}{1-x_{ij}x_{ji}A_{ii}B_{jj}}\left(\mathbf{1}^{\mathrm{T}}B\right)_{j},\end{split}$
(B.12)

which is Eq. (39).

## References

  * Lü  _et al._ (2016a) L. Lü, D. Chen, X.-L. Ren, Q.-M. Zhang, Y.-C. Zhang, and T. Zhou, Physics Reports 650, 1 (2016a).
  * Kitsak  _et al._ (2010) M. Kitsak, L. K. Gallos, S. Havlin, F. Liljeros, L. Muchnik, H. E. Stanley, and H. A. Makse, Nature Physics 6, 888 (2010).
  * Morone and Makse (2015) F. Morone and H. A. Makse, Nature 524, 65 (2015).
  * Lü  _et al._ (2016b) L. Lü, T. Zhou, Q.-M. Zhang, and H. E. Stanley, Nature Communications 7, 10168 (2016b).
  * Chen  _et al._ (2009) W. Chen, Y. Wang, and S. Yang, in  _Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining_ (ACM, 2009) pp. 199–208.
  * Chen  _et al._ (2010) W. Chen, C. Wang, and Y. Wang, in  _Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining_ (ACM, 2010) pp. 1029–1038.
  * Morone  _et al._ (2016) F. Morone, B. Min, L. Bo, R. Mari, and H. A. Makse, Scientific Reports 6, 30062 (2016).
  * Hu  _et al._ (2018) Y. Hu, S. Ji, Y. Jin, L. Feng, H. E. Stanley, and S. Havlin, Proceedings of the National Academy of Sciences of the United States of America , 115, 7468 (2018).
  * Ji  _et al._ (2017) S. Ji, L. Lü, C. H. Yeung, and Y. Hu, New Journal of Physics 19, 073020 (2017).
  * Ren  _et al._ (2014) Z.-M. Ren, A. Zeng, D.-B. Chen, H. Liao, and J.-G. Liu, Europhysics Letters 106, 48005 (2014).
  * Liao  _et al._ (2017) H. Liao, M. S. Mariani, M. Medo, Y.-C. Zhang, and M.-Y. Zhou, Physics Reports 689, 1 (2017).
  * Pei and Makse (2013) S. Pei and H. A. Makse, Journal of Statistical Mechanics: Theory and Experiment , P12002 (2013).
  * Chen  _et al._ (2012) D. Chen, L. Lü, M.-S. Shang, Y.-C. Zhang, and T. Zhou, Physica A 391, 1777 (2012).
  * Liu  _et al._ (2017) Y. Liu, M. Tang, Y. Do, and P. M. Hui, Physical Review E 96, 022323 (2017).
  * Yang  _et al._ (2008a) R. Yang, T. Zhou, Y.-B. Xie, Y.-C. Lai, and B.-H. Wang, Physical Review E 78, 066109 (2008a).
  * Gao  _et al._ (2016) L. Gao, W. Wang, L. Pan, M. Tang, and H.-F. Zhang, Scientific Reports 6, 38220 (2016).
  * Yang  _et al._ (2008b) R. Yang, L. Huang, and Y.-C. Lai, Physical Review E 78, 026111 (2008b).
  * Roshani and Naimi (2012) F. Roshani and Y. Naimi, Physical Review E 85, 036109 (2012).
  * Gao  _et al._ (2017) L. Gao, W. Wang, P. Shu, H. Gao, and L. A. Braunstein, Europhysics Letters 118, 18001 (2017).
  * Cui  _et al._ (2018) P.-B. Cui, W. Wang, S.-M. Cai, T. Zhou, and Y.-C. Lai, Physical Review E 98, 052311 (2018).
  * Aguirre  _et al._ (2013) J. Aguirre, D. Papo, and J. M. Buldú, Nature Physics 9, 230 (2013).
  * Del Ferraro  _et al._ (2018) G. Del Ferraro, A. Moreno, B. Min, F. Morone, Ú. Pérez-Ramírez, L. Pérez-Cervera, L. C. Parra, A. Holodny, S. Canals, and H. A. Makse, Nature Communications 9, 2274 (2018).
  * Milanese  _et al._ (2010) A. Milanese, J. Sun, and T. Nishikawa, Physical Review E 81, 046112 (2010).
  * Van Mieghem  _et al._ (2010) P. Van Mieghem, H. Wang, X. Ge, S. Tang, and F. Kuipers, The European Physical Journal B 76, 643 (2010).
  * Aguirre  _et al._ (2014) J. Aguirre, R. Sevilla-Escoboza, R. Gutiérrez, D. Papo, and J. Buldú, Physical Review Letters 112, 248701 (2014).
  * Li  _et al._ (2016) Y. Li, X. Wu, J.-A. Lu, and J. Lü, IEEE Transactions on Circuits and Systems II: Express Briefs 63, 206 (2016).
  * Wei  _et al._ (2018a) J. Wei, X. Wu, J.-A. Lu, and X. Wei, Europhysics Letters 120, 20005 (2018a).
  * Wei  _et al._ (2018b) X. Wei, J. Emenheiser, X. Wu, J.-A. Lu, and R. M. D’Souza, Chaos 28, 013110 (2018b).
  * Dai  _et al._ (2019) Y. Dai, J. Zhang, and W. Wang, Europhysics Letters 125, 18003 (2019).
  * Centola (2018) D. Centola,  _How behavior spreads: The science of complex contagions_ , Vol. 3 (Princeton University Press, 2018).
  * Guilbeault  _et al._ (2018) D. Guilbeault, J. Becker, and D. Centola, in  _Complex Spreading Phenomena in Social Systems_ (Springer, 2018) pp. 3–25.
  * Pastor-Satorras and Vespignani (2001) R. Pastor-Satorras and A. Vespignani, Physical Review Letters 86, 3200 (2001).
  * Romero  _et al._ (2011) D. M. Romero, B. Meeder, and J. Kleinberg, in  _Proceedings of the 20th International Conference on World Wide Web_ (ACM, 2011) pp. 695–704.
  * Watts (2002) D. J. Watts, Proceedings of the National Academy of Sciences of the United States of America 99, 5766 (2002).
  * Centola and Macy (2007) D. Centola and M. Macy, American Journal of Sociology 113, 702 (2007).
  * Aral and Nicolaides (2017) S. Aral and C. Nicolaides, Nature Communications 8, 14753 (2017).
  * Lü et al. (2011) L. Lü, D.-B. Chen, and T. Zhou, New Journal of Physics 13, 123005 (2011).
  * Unicomb  _et al._ (2019) S. Unicomb, G. Iñiguez, J. Kertész, and M. Karsai, arXiv: 1901.08306 (2019).
  * Karsai  _et al._ (2014) M. Karsai, G. Iñiguez, K. Kaski, and J. Kertész, Journal of The Royal Society Interface 11, 20140694 (2014).
  * Wang et al. (2019) W. Wang, Q.-H. Liu, J. Liang, Y. Hu, and T. Zhou, arXiv: 1901.02125 (2019).
  * da Silva  _et al._ (2018) P. C. V. da Silva, F. Velásquez-Rojas, C. Connaughton, F. Vazquez, Y. Moreno, and F. A. Rodrigues, arXiv: 1812.01386 (2018).
  * De Domenico  _et al._ (2016) M. De Domenico, C. Granell, M. A. Porter, and A. Arenas, Nature Physics 12, 901 (2016).
  * Granell  _et al._ (2013) C. Granell, S. Gómez, and A. Arenas, Physical Review Letters 111, 128701 (2013).
  * Granell  _et al._ (2014) C. Granell, S. Gómez, and A. Arenas, Physical Review E 90, 012808 (2014).
  * Wang  _et al._ (2014) W. Wang, M. Tang, H. Yang, Y. Do, Y.-C. Lai, and G. Lee, Scientific Reports 4, 5097 (2014).
  * Chen  _et al._ (2018) X. Chen, W. Wang, S. Cai, H. E. Stanley, and L. A. Braunstein, Journal of Statistical Mechanics: Theory and Experiment 053501 (2018).
  * Gao  _et al._ (2012) J. Gao, S. V. Buldyrev, H. E. Stanley, and S. Havlin, Nature Physics 8, 40 (2012).
  * Tejedor  _et al._ (2018) A. Tejedor, A. Longjas, E. Foufoula-Georgiou, T. T. Georgiou, and Y. Moreno, Physical Review X 8, 031071 (2018).
  * Wang  _et al._ (2018) W. Wang, X.-L. Chen, and L.-F. Zhong, Physica A 503, 604 (2018).
  * Radicchi and Arenas (2013) F. Radicchi and A. Arenas, Nature Physics 9, 717 (2013).
  * Reis  _et al._ (2014) S. D. Reis, Y. Hu, A. Babino, J. S. Andrade Jr, S. Canals, M. Sigman, and H. A. Makse, Nature Physics 10, 762 (2014).
  * Van Mieghem (2016) P. Van Mieghem, Physical Review E 93, 042305 (2016).
  * Cozzo  _et al._ (2019) E. Cozzo, G. F. de Arruda, F. A. Rodrigues, and Y. Moreno, arXiv: 1901.04523 (2019).
  * Zhang  _et al._ (2015) X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Physical Review Letters 114, 038701 (2015).
  * Hu  _et al._ (2014) Y. Hu, S. Havlin, and H. A. Makse, Physical Review X 4, 021031 (2014).
  * de Arruda  _et al._ (2017) G. F. de Arruda, E. Cozzo, T. P. Peixoto, F. A. Rodrigues, and Y. Moreno, Physical Review X 7, 011014 (2017).
  * Sanz  _et al._ (2014) J. Sanz, C.-Y. Xia, S. Meloni, and Y. Moreno, Physical Review X 4, 041005 (2014).
  * Saumell-Mendiola  _et al._ (2012) A. Saumell-Mendiola, M. Á. Serrano, and M. Boguná, Physical Review E 86, 026106 (2012).
  * Pan  _et al._ (2019) L. Pan, W. Wang, S. Cai, and T. Zhou, Unpublished (2019).
  * Wang  _et al._ (2016a) Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, and D. Zhao, Physics Reports 664, 1 (2016a).
  * Gómez  _et al._ (2010) S. Gómez, A. Arenas, J. Borge-Holthoefer, S. Meloni, and Y. Moreno, Europhysics Letters 89, 38009 (2010).
  * Wang  _et al._ (2016b) W. Wang, M. Tang, H. E. Stanley, and L. A. Braunstein, Reports on Progress in Physics 80, 036603 (2016b).
  * Newman (2010) M. E. J. Newman,  _Networks: an introduction_ (Oxford University Press, 2010).
  * Lee  _et al._ (2012) K.-M. Lee, J. Y. Kim, W.-k. Cho, K.-I. Goh, and I. Kim, New Journal of Physics 14, 033027 (2012).
  * Goltsev  _et al._ (2012) A. V. Goltsev, S. N. Dorogovtsev, J. Oliveira, and J. F. Mendes, Physical Review Letters 109, 128702 (2012).
  * Massa  _et al._ (2009) P. Massa, M. Salvetti, and D. Tomasoni, in  _Dependable, Autonomic and Secure Computing, 2009. DASC’09. Eighth IEEE International Conference on_ (IEEE, 2009) pp. 658–663.
  * Leskovec and Mcauley (2012) J. Leskovec and J. J. Mcauley, in  _Advances in neural information processing systems_ (2012) pp. 539–547.
  * (68) http://konect.uni-koblenz.de/.
  * Moody (2001) J. Moody, Social Networks 23, 261 (2001).
  * Coleman  _et al._ (1957) J. Coleman, E. Katz, and H. Menzel, Sociometry 20, 253 (1957).