Synergistic Effects in Networked Epidemic Spreading Dynamics

TL;DR

This paper introduces a rigorous method to analyze how synergistic interactions influence epidemic spreading in complex networks, providing bounds on growth rates and validating results with real social network data.

Contribution

It develops a novel eigenvalue-based approach to quantify synergy effects in epidemic models on arbitrary networks, advancing beyond previous approximation methods.

Findings

Derived an upper bound on epidemic growth rate considering synergy.

Validated the theoretical bounds with simulations on real social networks.

Showed that synergy significantly accelerates epidemic spreading.

Abstract

In this brief, we study epidemic spreading dynamics taking place in complex networks. We specifically investigate the effect of synergy, where multiple interactions between nodes result in a combined effect larger than the simple sum of their separate effects. Although synergistic effects play key roles in various biological and social phenomena, their analyses have been often performed by means of approximation techniques and for limited types of networks. In order to address this limitation, this paper proposes a rigorous approach to quantitatively understand the effect of synergy in the Susceptible-Infected-Susceptible model taking place in an arbitrary complex network. We derive an upper bound on the growth rate of the synergistic Susceptible-Infected-Susceptible model in terms of the eigenvalues of a matrix whose size grows quadratically with the number of the nodes in the network. We confirm the effectiveness of our result by numerical simulations on empirically observed human and animal social networks.

Equations

$${\beta }_i+{\gamma }_j{m}_j(t),$$

$$\frac{d{p}_i}{dt}$$

$$\frac{d{p}_{i\ell }}{dt}$$

$$\displaystyle\quad\quad\quad-\delta_{i}p_{i\ell}+\beta_{i}a_{i\ell}(p_{\ell}-p_{i\ell})+\sum_{j<\ell}a_{ji}\beta_{i}p_{j\ell}+\sum_{j>\ell}a_{ji}\beta_{i}p_{\ell j}+\sum_{j=1}^{N}\Biggl{(}\sum_{\begin{subarray}{c}k<j\\ k\neq i\end{subarray}}a_{ji}a_{kj}\gamma_{j}p_{kj}+\sum_{\begin{subarray}{c}k>j\\ k\neq i\end{subarray}}a_{ji}a_{kj}\gamma_{j}p_{jk}\Biggr{)}.$$

$$\begin{multlined}dx_{i}=-x_{i}d\Pi_{\delta_{i}}+(1-x_{i})\sum_{j\in\mathscr{N}_{i}}x_{j}d\Pi_{\beta_{i}}\\ +(1-x_{i})\sum_{j\in\mathscr{N}_{i}}\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}x_{j}x_{k}d\Pi_{\gamma_{j}},\end{multlined}dx_{i}=-x_{i}d\Pi_{\delta_{i}}+(1-x_{i})\sum_{j\in\mathscr{N}_{i}}x_{j}d\Pi_{\beta_{i}}\\ +(1-x_{i})\sum_{j\in\mathscr{N}_{i}}\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}x_{j}x_{k}d\Pi_{\gamma_{j}},$$

$$\begin{multlined}\frac{d}{dt}\mathbb{E}[x_{i}]=-\delta_{i}\mathbb{E}[x_{i}]+\sum_{j\in\mathscr{N}_{i}}\beta_{i}\mathbb{E}[(1-x_{i})x_{j}]\\ +\sum_{j\in\mathscr{N}_{i}}\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}\gamma_{j}\mathbb{E}[(1-x_{i})x_{j}x_{k}].\end{multlined}\frac{d}{dt}\mathbb{E}[x_{i}]=-\delta_{i}\mathbb{E}[x_{i}]+\sum_{j\in\mathscr{N}_{i}}\beta_{i}\mathbb{E}[(1-x_{i})x_{j}]\\ +\sum_{j\in\mathscr{N}_{i}}\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}\gamma_{j}\mathbb{E}[(1-x_{i})x_{j}x_{k}].$$

$$\begin{multlined}\frac{dp_{i}}{dt}=-\delta_{i}p_{i}+\sum_{j\in\mathscr{N}_{i}}\beta_{i}(p_{j}-p_{ij})\\ +\sum_{j\in\mathscr{N}_{i}}\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}\gamma_{j}(p_{jk}-p_{ijk}).\end{multlined}\frac{dp_{i}}{dt}=-\delta_{i}p_{i}+\sum_{j\in\mathscr{N}_{i}}\beta_{i}(p_{j}-p_{ij})\\ +\sum_{j\in\mathscr{N}_{i}}\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}\gamma_{j}(p_{jk}-p_{ijk}).$$

$$\frac{d{p}_i}{dt}\le -{\delta }_i{p}_i+{\beta }_i\sum _{j=1}^N{a}_{ij}{p}_j+\sum _{j=1}^N\sum _{k\ne i}{a}_{ij}{a}_{kj}{\gamma }_j{p}_{jk}.$$

$$\begin{multlined}\frac{d}{dt}\mathbb{E}[x_{i}x_{\ell}]=\!\!-\delta_{\ell}\mathbb{E}[x_{i}x_{\ell}]+\sum_{m=1}^{N}a_{m\ell}\beta_{\ell}\mathbb{E}[x_{i}(1-x_{\ell})x_{m}]\\ +\sum_{m=1}^{N}\sum_{n\neq\ell}a_{m\ell}a_{nm}\gamma_{m}\mathbb{E}[x_{i}(1-x_{\ell})x_{m}x_{n}]\\ -\delta_{i}\mathbb{E}[x_{i}x_{\ell}]+\sum_{j=1}^{N}a_{ji}\beta_{i}\mathbb{E}[(1-x_{i})x_{j}x_{\ell}]\\ +\sum_{j=1}^{N}\sum_{k\neq i}a_{ji}a_{kj}\gamma_{j}\mathbb{E}[(1-x_{i})x_{j}x_{k}x_{\ell}].\end{multlined}\frac{d}{dt}\mathbb{E}[x_{i}x_{\ell}]=\!\!-\delta_{\ell}\mathbb{E}[x_{i}x_{\ell}]+\sum_{m=1}^{N}a_{m\ell}\beta_{\ell}\mathbb{E}[x_{i}(1-x_{\ell})x_{m}]\\ +\sum_{m=1}^{N}\sum_{n\neq\ell}a_{m\ell}a_{nm}\gamma_{m}\mathbb{E}[x_{i}(1-x_{\ell})x_{m}x_{n}]\\ -\delta_{i}\mathbb{E}[x_{i}x_{\ell}]+\sum_{j=1}^{N}a_{ji}\beta_{i}\mathbb{E}[(1-x_{i})x_{j}x_{\ell}]\\ +\sum_{j=1}^{N}\sum_{k\neq i}a_{ji}a_{kj}\gamma_{j}\mathbb{E}[(1-x_{i})x_{j}x_{k}x_{\ell}].$$

$$\mathbb{E}[x_i(1-x_{\ell })x_m{x}_n]\le p_{mn},\mathbb{E}[(1-x_i)x_j{x}_k{x}_{\ell }]\le p_{jk}.$$

$$\begin{multlined}\frac{dp_{i\ell}}{dt}\leq-\delta_{\ell}p_{i\ell}+a_{i\ell}\beta_{\ell}(p_{i}-p_{i\ell})\\ +\sum_{m\neq i}a_{m\ell}\beta_{\ell}p_{im}+\sum_{m=1}^{N}\sum_{n\neq\ell}a_{m\ell}a_{nm}\gamma_{m}p_{mn}\\ -\delta_{i}p_{i\ell}+a_{\ell i}\beta_{i}(p_{\ell}-p_{i\ell})\\ +\sum_{j\neq\ell}a_{ji}\beta_{i}p_{j\ell}+\sum_{j=1}^{N}\sum_{k\neq i}a_{ji}a_{kj}\gamma_{j}p_{kj}.\end{multlined}\frac{dp_{i\ell}}{dt}\leq-\delta_{\ell}p_{i\ell}+a_{i\ell}\beta_{\ell}(p_{i}-p_{i\ell})\\ +\sum_{m\neq i}a_{m\ell}\beta_{\ell}p_{im}+\sum_{m=1}^{N}\sum_{n\neq\ell}a_{m\ell}a_{nm}\gamma_{m}p_{mn}\\ -\delta_{i}p_{i\ell}+a_{\ell i}\beta_{i}(p_{\ell}-p_{i\ell})\\ +\sum_{j\neq\ell}a_{ji}\beta_{i}p_{j\ell}+\sum_{j=1}^{N}\sum_{k\neq i}a_{ji}a_{kj}\gamma_{j}p_{kj}.\vspace{3mm}$$

$$\frac{dp}{dt}\le \mathcal{M}p$$

$${\lambda }_{\max }(\mathcal{M})<0.$$

$$y^{\ast }=y(t)-1,$$

$$\mathcal{E}=\{(\delta ,\beta ):y^{\star }<1\}$$

$$\underaccent{\bar}{\mathcal E}=\{(\delta,\beta)\colon\lambda_{\max}(\mathscr{M})<0\}.$$

$$\underaccent{\bar}{\mathcal E}_{\operatorname*{SIS}}=\{(\delta,\beta)\colon\bar{\rho}_{\operatorname*{SIS}}<0\}$$

Full text

I Introduction

Understanding epidemic spreading dynamics taking place in complex networks is of utmost importance with applications in epidemiology [1], opinion formation, rumor propagation, and cyber security [2, 3]. During the last two decades, several important advances have been made toward the analysis, modeling, and control of networked spreading processes [4]. For example, fundamental connections between the eigenvalues of a network and epidemic thresholds were identified in [5]. A unified framework for the analysis of spreading processes evolving on multilayer networks was presented in [6]. The dynamical change of network structures dependent on the evolution of spreading processes has been actively investigated in the context of adaptive networks [7]. Recently, control and optimization tools have been proven to be effective for containing spreading processes in static [8], multi-layer [9], and temporal [10] networks. A prominent example of the achievements of network epidemiology can be found in, for example, effective predictions of the pandemic of H1N1 influenza on 2009 [11, 12].

All the aforementioned results assume nearest-neighbor interactions, in which individual nodes are directly influenced by only their neighbors. However, in realistic epidemic spreading processes, we frequently observe multi-node interactions, in which a node can affect not only its neighbors but also other nodes such as the neighbors of its neighbors. A typical example of such phenomena can be found in the colonization by fungal and bacterial pathogens [13] and tumor growth [14]. In the context of human social networks, it has been experimentally confirmed that the spreads of opinions [15] and behaviors [16] are effectively accelerated by synergistic effects, which make the combined effect from multiple interactions greater than the simple sum of the effect of individual interactions.

We find in the literature several works providing theoretical insight for understanding the effect of synergistic interactions on spreading processes. The authors in [17] considered an extended Susceptible-Infected-Recovered (SIR) model with synergistic effects, and numerically confirmed on lattices that synergy can significantly alter asymptotic behaviors of the spread. Similar models were further investigated analytically on lattices [18] and small-world networks [19], respectively. The authors in [20] studied the effect of synergy in the Susceptible-Infected-Susceptible (SIS) model under the assumption that the nodes having the same degree behave in the same manner. The effect of synergy on deterministic spreading processes was investigated in [21] with a mean-field approximation technique.

This brief proposes a rigorous but tractable approach to quantitatively analyze the effect of synergy in the SIS model taking place in arbitrary networks, without resorting to mean-field approximations. A primary challenge in this context is to obtain computationally tractable results, because the state space of the spreading model grows exponentially with the number of nodes in the network. To overcome this difficulty, in this brief we show that the growth rate of the size of the infected population in the network can be bounded from above by the maximum real part of the eigenvalues of a matrix whose size grows only quadratically with the number of the nodes. In the derivation of the upper bound, we use Poisson-type stochastic differential equations describing the spreading model. The effectiveness of our theoretical results is numerically illustrated by simulations on empirical human and animal social networks.

This brief is organized as follows. After preparing mathematical notations, in Section II we introduce the synergistic SIS model over arbitrary complex networks. In Section III, we derive our upper-bound on the growth rate of the model. Numerical simulations are presented in Section IV.

Notations

The maximum real part of the eigenvalues of a square matrix $A$ is denoted by $\lambda_{\max}(A)$. An undirected graph is defined as the pair $G=(V,E)$, where $V=\{{1},\dotsc,{N}\}$ is a set of nodes and $E$ is a set of edges, i.e., unordered pairs of nodes. We say that nodes $i$ and $j$ are adjacent if $\{i,j\}\in\mathscr{E}$. The neighborhood of node $i$, denoted by $\mathscr{N}_{i}$, is defined as the set of nodes adjacent to node $i$. The adjacency matrix $A=[a_{ij}]_{i,j}$ of $G$ is defined as the $N\times N$ matrix such that $a_{ij}=1$ if nodes $i$ and $j$ are adjacent, and $a_{ij}=0$ otherwise.

II Synergistic SIS model

We introduce the synergistic SIS model over complex networks in this section. Let us start by reviewing the standard SIS model without synergy (see, e.g., [22]). Let $G$ be an undirected graph having $N$ nodes labeled as $i=1,\dotsc,N$, with edges representing interactions between nodes. At time $t\geq 0$, each node can be either susceptible or infected. When a node $i$ is infected, the node transitions into the susceptible state with an instantaneous rate $\delta_{i}$. Also, if node $i$ is susceptible and has an infected neighbor, then the neighbor can infect node $i$ with an instantaneous rate $\beta_{i}$. The rates $\delta_{i}$ and $\beta_{i}$ are called the recovery and transmission rate of node $i$, respectively. We suppose that the above model evolves in the continuous-time; therefore, the rates represent probabilities per unit time [1] and, therefore, can be greater than one.

The SIS model uses a pair-wise description for infection events and, therefore, does not allow synergistic effects. In order to study the effect of synergy, we study the synergistic SIS model [23], which we describe below. In the model, while an infected node recovers with the rate $\delta_{i}$ as in the SIS model, the rate at which a susceptible node $i$ is infected by its infected neighbor, say, node $j$, is equal to

[TABLE]

where $\gamma_{j}$ denotes the strength of synergy and is dependent on individual nodes, and $m_{j}(t)$ denotes the number of the infected neighbors of node $j$ at time $t$. The second term in (1) is introduced to describe synergistic effects: a susceptible node is directly affected not only by its infected neighbors but also their neighbors (see Fig. 1 for a schematic picture). We remark that, according to the classification in [17], the synergistic SIS model incorporates only the d-synergies, which are caused by adjacent pairs of infected nodes. Due to the limitations of the space, this brief does not deal with the other category of synergies called r-synergies arising in the interactions between a susceptible node and its infected neighbors.

In this brief, we are interested in quantifying the rate at which the spreading process described by the synergistic SIS model grows. Therefore, we introduce the following definition:

Definition 1**.**

For all node $i$ and time $t\geq 0$, let $p_{i}(t)$ denote the probability that node $i$ is infected at time $t$. Define the growth rate of the synergistic SIS model by $\rho=\sup_{V_{0}\subset\{1,\dotsc,n\}}\limsup_{t\to\infty}t^{-1}{\log\sum_{i=1}^{N}p_{i}(t)}$, where $V_{0}$ denotes the set of initially infected nodes.

In Definition 1, the sum $\sum_{i=1}^{N}p_{i}(t)$ represents the expected number of infected nodes. Therefore, the expression $\limsup_{t\to\infty}t^{-1}\log\sum_{i=1}^{N}p_{i}(t)$ measures the asymptotic growth of the expected number of infected nodes. Hence, the growth rate in Definition 1 quantifies the worst-case growth of the average size of the infected population.

The growth rate is closely related to fundamental notions in the network epidemiology such as mean-time-to-absorption [5] and epidemic thresholds [22]. However, it is hard in practice to exactly compute the growth rate. Since the synergistic SIS model, as a Markov process, has a total of $2^{N}$ states, the computation of the growth rate requires finding the eigenvalues of a $2^{N}\times 2^{N}$ matrix representing the infinitesimal generator of the Markov process [22]. Therefore, the computational complexity of exactly calculating the growth rate increases exponentially with the number of the nodes [24]. To overcome this difficulty, in the next section we introduce an upper-bound on the growth rate in terms of the eigenvalues of a matrix whose size grows only quadratically with $N$.

III Growth rate

In this section, we derive an upper bound on the growth rate of the synergistic SIS model. For each node $i$ and time $t$, let us define the variable $x_{i}(t)$ as $x_{i}(t)=1$ if node $i$ is infected at time $t$, and $x_{i}(t)=0$ otherwise. Then, since $m_{j}(t)=\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}x_{k}(t)$, the synergistic rate of infection (1) shows that the transition probability of the variable $x_{i}$ is described as $\Pr(x_{i}(t+h)=1\mid x_{i}(t)=0)=\sum_{j\in\mathscr{N}_{i}}x_{j}(t)(\beta_{i}+\gamma_{j}\sum_{k\in\mathscr{N}_{j}\backslash\{i\}}x_{k}(t))h+o(h)$ and $\Pr(x_{i}(t+h)=0\mid x_{i}(t)=1)=\delta_{i}h+o(h)$ for $h>0$. This implies that the variable $x_{i}$ obeys the following stochastic differential equation with Poisson jumps (see, e.g., [25]):

[TABLE]

where $\Pi_{\delta_{i}}$, $\Pi_{\beta_{i}}$, and $\Pi_{\gamma_{j}}$ are independent Poisson counters with rates $\delta_{i}$, $\beta_{i}$, and $\gamma_{j}$, respectively. In the stochastic differential equation (2), the first, second, and third terms represent the recovery, direct infection from neighbors, and synergy effects.

By (2), the expectation $\mathbb{E}[x_{i}]$ satisfies

[TABLE]

Notice that we have $p_{i}(t)=\mathbb{E}[x_{i}(t)]$. Also, let us introduce the notations $p_{ij}(t)=\mathbb{E}[x_{i}(t)x_{j}(t)]$ and $p_{ijk}(t)=\mathbb{E}[x_{i}(t)x_{j}(t)x_{k}(t)]$ for all $i,j,k\in\{1,\dotsc,N\}$ and $t\geq 0$. Then, the differential equation (3) is rewritten as

[TABLE]

This differential equation rigorously describes the evolution of the infection probabilities $p_{i}$, which we require to compute the growth rate $\rho$. However, the differential equation contains the higher-order terms $p_{ij}$ and $p_{ijk}$ and, therefore, is not closed in the first-order term $p_{i}$. For this reason, we cannot readily use the differential equation for calculating the growth rate.

We avoid this difficulty by ignoring the third order terms $p_{ijk}$ and deriving a set of differential inequalities for upper-bounding the growth rate of the synergistic SIS model. Intuitively, we can expect that ignoring the higher-order terms brings relatively small errors in the regime of extinction (i.e., when the epidemics dies out), in which lower-order terms become dominant. From differential equation (4), we obtain

[TABLE]

It should be emphasized that, in deriving inequality (5), we chose to ignore the negative second order terms $-\beta_{i}p_{ij}$ for $j\in\mathscr{N}_{i}$ as well as the third-order terms. We discuss the reason for and consequence from this manipulation later in Remark 2.

Then, let us derive differential equations for the evolution of the second-order terms. If we apply Itô’s formula [25] to the mapping $(x_{i},x_{\ell})\mapsto x_{i}x_{\ell}$ and, then, take mathematical expectations in the resulting stochastic differential equation, we obtain

[TABLE]

Since $x_{i}$ is $\{0,1\}$-valued, the second term on the right-hand side of the differential equation (6) is bounded from above as $\sum_{m=1}^{N}a_{m\ell}\beta_{\ell}\mathbb{E}[x_{i}(1-x_{\ell})x_{m}]\leq a_{i\ell}\beta_{\ell}(p_{i}-p_{i\ell})+\sum_{m\neq i}a_{m\ell}\beta_{\ell}p_{im}$. Likewise, the fifth term on the right-hand side of (6) can be bounded as $\sum_{j=1}^{N}a_{ji}\beta_{i}\mathbb{E}[(1-x_{i})x_{j}x_{\ell}]\leq a_{\ell i}\beta_{i}(p_{\ell}-p_{i\ell})+\sum_{j\neq\ell}a_{ji}\beta_{i}p_{j\ell}$. Furthermore, we have

[TABLE]

These observations yield that

[TABLE]

Remark 1*.*

One can derive different bounds on the derivative $dp_{i\ell}/dt$ by replacing the reductions in (7) with others such as $\mathbb{E}[x_{i}(1-x_{\ell})x_{m}x_{n}]\leq p_{im}$ and $\mathbb{E}[(1-x_{i})x_{j}x_{k}x_{\ell}]\leq p_{j\ell}$. It is not theoretically easy to identify the best reduction giving the tightest bound on the derivative. However, we have numerically observed that using different reductions does not significantly alter the results of numerical simulations performed in Section IV.

By replacing second-order moments $p_{ij}$ having indices $i>j$ with their equivalent counterparts $p_{ji}=E[x_{j}x_{j}]=E[x_{i}x_{j}]=p_{ij}$, we rewrite the differential inequality (5) as (9) for all $i$. In the same manner, we rewrite the differential inequality (8) as (10), which contains only the variables belonging to the set $\mathscr{P}=\{p_{i},p_{jk}\colon i,j,k\in\{1,\dotsc,N\}\mbox{ and }j<k\}$. Now, let $p$ denote the $\frac{N^{2}+N}{2}$-dimensional vectorial variable that can be obtained by stacking the variables in the set $\mathscr{P}$. Then, we rewrite the system of differential inequalities (9) and (10) as

[TABLE]

for a uniquely determined $\frac{N^{2}+N}{2}\times\frac{N^{2}+N}{2}$ matrix $\mathscr{M}$.

Remark 2*.*

Since we ignored second order terms $-\beta_{i}p_{ij}$ for deriving (5), one should expect a certain conservativeness in the bound (11) and, therefore, the bound of the decay rate to be later derived in Theorem 1. However, it is necessary to ignore the negative second order terms because, if not ignored, we cannot use the fact that the matrix $\mathscr{M}$ is Metzler for the proof of the theorem.

Inequality (11) shows that, for all $t\geq 0$, there exists a nonnegative vector $\epsilon(t)$ such that ${dp}/{dt}=\mathscr{M}p-\epsilon$. From differential inequalities (9) and (10), we see that the off-diagonals of $\mathscr{M}$ are nonnegative, i.e., $\mathscr{M}$ is a Metzler matrix. This implies that the exponential $\exp(\mathscr{M}t)$ has only nonnegative entries for all $t\geq 0$ (see, e.g., [26]). Furthermore, by definition, each entry of the vector $\epsilon(t)$ is nonnegative for all $t\geq 0$. Therefore, we obtain $p(t)=\exp(\mathscr{M}t)p(0)-\int_{0}^{t}\exp(\mathscr{M}(t-\tau))\epsilon(\tau)\,d\tau\leq\exp(\mathscr{M}t)p(0)$. This inequality bounds the vector $p$, which contains the infection probabilities $p_{1}$, …, $p_{N}$, by the matrix exponential $\exp(\mathscr{M}t)$. We have thus proved the following theorem, which is the main theoretical result of this brief:

Theorem 1**.**

The growth rate $\rho$ of the synergistic SIS model satisfies $\rho\leq\lambda_{\max}(\mathscr{M})$.

IV Numerical simulations

In this section, we illustrate the effectiveness of Theorem 1 by numerical simulations. From Theorem 1, we can obtain a sufficient condition for the spreading process to become extinct (i.e., the spread is not invasive) as

[TABLE]

This section is devoted to numerically validating the sufficiency of the extinction condition (12). To determine if the spreading process becomes extinct, we use the notion of the meta-stable number of infected nodes that is commonly used in the Network Science. To identify the meta-stable number of infected nodes, we begin by computing the following long-time average of the number of infected nodes $y(t)=t^{-1}\int_{0}^{t}\sum_{i=1}^{N}x_{i}(\tau)\,d\tau$ for a sufficiently large $t$. In simulations, the spreading process can die out due to random fluctuations. To prevent this from happening, we use the re-infection procedure [27] in which a randomly chosen node is infected immediately after the spreading process dies. After we finish the simulation, we determine the meta-stable number of infected nodes, denoted by $y^{\star}$, by

[TABLE]

where we subtract one to compensate for the effect of the re-infection procedure. Finally, when $y^{\star}<1$, it is determined that the spreading process becomes extinct.

In our numerical simulations, we use the following three empirical human and animal social networks: (1) a social network of Zachary’s Karate club [28] (Karate, $N=34$), (2) a social network of bottlenose dolphins [29] (Dolphin, $N=62$), and (3) a friendship network in Facebook (Facebook, $N=247$). For simplicity of illustration, we assume that the network consists of homogeneous nodes; i.e., we assume $\beta_{i}=\beta$, $\delta_{i}=\delta$, and $\gamma_{j}=\gamma$ for constants $\beta$, $\delta$, and $\gamma$. In this numerical simulation, we fix $\gamma=10^{-2}$. For the above three networks and various values of the transmission rate $\beta$ and the recovery rate $\delta$, we compute the meta-stable number $y^{\star}$ of infected nodes using the procedure described above, as well as our upper-bound $\lambda_{\max}(\mathscr{M})$ on the growth rate. We then obtain the true extinction region

[TABLE]

and its inner estimate from the sufficient condition (12):

[TABLE]

In Fig. 3, we illustrate the regions $\mathscr{E}$ and \underaccent{\bar}{\mathcal E} for the three empirical networks. We confirm that, for any of the networks, our inner-estimate \underaccent{\bar}{\mathcal E} is contained in the true extinction region $\mathscr{E}$ (i.e., the boundary \partial\underaccent{\bar}{\mathcal E} indicated by the white line passes through the colored region). For comparison, in Fig. 3, we also illustrate the estimate

[TABLE]

of the extinction region based on the conventional upper-bound $\bar{\rho}_{\operatorname*{SIS}}=\beta\lambda_{\max}(A)-\delta$ on the growth rate of the SIS model without synergy [22]. Since the boundary \partial\underaccent{\bar}{\mathcal E}_{\operatorname*{SIS}} passes outside the true extinction region for any of the three graphs, we observe that the use of the conventional method, which does not take into account synergistic effects, can lead to the overestimate of the extinction regions. In Fig. 3, we show sample paths of the time-evolution of the number of infected nodes for $\beta=0.02$ and $\delta=0.1$, $0.3$, $1$, and $3$. From Fig. 3.b, we confirm that the conventional method does overestimate the region of extinction.

V Conclusion

In this brief, we have studied the synergistic SIS model over arbitrary complex networks. Since the state space of the model grows exponentially fast with the number of nodes, it is almost impossible to directly analyze the model except when the network size is very small. To overcome this difficulty, we have derived an upper-bound on the growth rate the spreading model in terms of the eigenvalues of a matrix whose size grows quadratically with the number of the nodes. We have illustrated the effectiveness of our theoretical results by numerical simulations on empirical human and animal social networks. A possible direction for future research is eliminating the computational complexity for calculating our upper-bound. Since the size of $\mathscr{M}$ grows quadratically, it is still not easy to apply the current bound for large-scale networks. It is of theoretical interest to examine if we can exploit the sparsity structure of $\mathscr{M}$ to efficiently calculate the upper bound.