Analysis of the susceptible-infected-susceptible epidemic dynamics in networks via the non-backtracking matrix

TL;DR

This paper introduces a new, computationally efficient lower bound on the decay rate of infection probabilities in SIS epidemic models on networks, improving upon previous bounds and relating to the network's non-backtracking matrix.

Contribution

The paper presents a novel lower bound based on the non-backtracking matrix that outperforms existing bounds and is more computationally feasible for sparse networks.

Findings

The new bound is tighter than previous bounds.

It is computationally less expensive for sparse networks.

The epidemic threshold is expressed in terms of the line graph and non-backtracking matrices.

Abstract

We study the stochastic susceptible-infected-susceptible model of epidemic processes on finite directed and weighted networks with arbitrary structure. We present a new lower bound on the exponential rate at which the probabilities of nodes being infected decay over time. This bound is directly related to the leading eigenvalue of a matrix that depends on the non-backtracking and incidence matrices of the network. The dimension of this matrix is N+M, where N and M are the number of nodes and edges, respectively. We show that this new lower bound improves on an existing bound corresponding to the so-called quenched mean-field theory. Although the bound obtained from a recently developed second-order moment-closure technique requires the computation of the leading eigenvalue of an N^2 x N^2 matrix, we illustrate in our numerical simulations that the new bound is tighter, while being computationally less expensive for sparse networks. We also present the expression for the corresponding epidemic threshold in terms of the adjacency matrix of the line graph and the non-backtracking matrix of the given network.

Equations

$$\gamma =-\underset{t\to \infty }{\mathrm{limsup}}\frac{\log {\sum }_{i=1}^N{p}_i(t)}{t},$$

$$B_{ij}=\begin{cases}\beta_{ij},&\mbox{if $(v_{i},v_{j})\in\mathcal{E}$,}\\ 0,&\mbox{otherwise,}\end{cases}$$

$$D=\text{diag}({\delta }_1,\dots ,{\delta }_N),$$

$$\gamma \ge {\gamma }_1\equiv -{\lambda }_{\max }(B^{\top }-D),$$

$$C_{i\ell}=\begin{cases}1,&\mbox{if $j_{\ell}=i$},\\ -1,&\mbox{if $i_{\ell}=i$},\\ 0,&\mbox{otherwise}.\end{cases}$$

$$H_{\ell m}=\begin{cases}1,&\mbox{if $j_{\ell}=i_{m}$ and $j_{m}\neq i_{\ell}$},\\ 0,&\mbox{otherwise.}\end{cases}$$

$$\mathcal{A}=\left[\begin{array}{cc}-D& C_{+}{B}^{\prime }\\ D_2^{\prime }{C}_{-}^{\top }& H^{\top }{B}^{\prime }-B^{\prime }-D_1^{\prime }-D_2^{\prime }\end{array}\right],$$

$$B^{\prime }=$$

$$D_1^{\prime }=$$

$$D_2^{\prime }=$$

$$\gamma \ge {\gamma }_2\equiv -{\lambda }_{\max }(\mathcal{A}).$$

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

$$=\left(\sum _{j=1}^N{\beta }_{ji}{q}_{ji}\right)-{\delta }_i{p}_i,$$

$$q_{ji}(t)=E[x_j(t)(1-x_i(t))]$$

$$\sum _{j=1}^N{\beta }_{ji}{q}_{ji}=\sum _{\ell =1;j_{\ell }=i}^M{\beta }_{i_{\ell }{j}_{\ell }}{q}_{i_{\ell }{j}_{\ell }}=\sum _{\ell =1}^M[C_{+}{]}_{i\ell }[B^{\prime }{]}_{\ell \ell }{q}_{i_{\ell }{j}_{\ell }},$$

$$\frac{d{p}_i}{dt}=\left(\sum _{\ell =1}^M[C_{+}{B}^{\prime }{]}_{i\ell }{q}_{i_{\ell }{j}_{\ell }}\right)-{\delta }_i{p}_i(i\in \{1,\dots ,N\}).$$

$$\frac{d\mathbf{p}}{dt}=C_{+}{B}^{\prime }\mathbf{q}-D\mathbf{p},$$

$$\mathbf{q}(t)\equiv {\left[q_{i_1{j}_1}(t),\dots ,q_{i_M{j}_M}(t)\right]}^{\top }.$$

$$\frac{d{q}_{ij}}{dt}$$

$$+\left(\sum _{k=1}^{N}E[(1-x_i)(1-x_j)x_k]{\beta }_{ki}\right)-{\delta }_{i}E[x_i(1-x_j)]$$

$$\le -{\beta }_{ij}{q}_{ij}+{\delta }_j{p}_{ij}+\left(\sum _{k=1;k\ne j}^N{\beta }_{ki}{q}_{ki}\right)-{\delta }_i{q}_{ij}.$$

$$\frac{d{q}_{i_{\ell }{j}_{\ell }}}{dt}\le -({\beta }_{i_{\ell }{j}_{\ell }}+{\delta }_{i_{\ell }}+{\delta }_{j_{\ell }})q_{i_{\ell }{j}_{\ell }}+{\delta }_{j_{\ell }}{p}_{i_{\ell }}+\sum _{k=1;k\ne j_{\ell }}^N{\beta }_{k{i}_{\ell }}{q}_{k{i}_{\ell }}.$$

$$[C_{-}^{\top }\mathbf{p}{]}_{\ell }=\sum _{i=1}^N[C_{-}{]}_{i\ell }{p}_i=p_{i_{\ell }}$$

$$\begin{multlined}{}_{\ell}=\sum_{m=1}^{M}H_{m\ell}B^{\prime}_{mm}q_{m}=\sum^{M}_{\begin{subarray}{c}m=1;j_{m}=i_{\ell},\\ j_{\ell}\neq i_{m}\end{subarray}}\beta_{i_{m}j_{m}}q_{i_{m}j_{m}}\\ =\sum_{k=1;k\neq j_{\ell}}^{N}\beta_{ki_{\ell}}q_{ki_{\ell}},\end{multlined}{}_{\ell}=\sum_{m=1}^{M}H_{m\ell}B^{\prime}_{mm}q_{m}=\sum^{M}_{\begin{subarray}{c}m=1;j_{m}=i_{\ell},\\ j_{\ell}\neq i_{m}\end{subarray}}\beta_{i_{m}j_{m}}q_{i_{m}j_{m}}\\ =\sum_{k=1;k\neq j_{\ell}}^{N}\beta_{ki_{\ell}}q_{ki_{\ell}},$$

$$\frac{d{q}_{i_{\ell }{j}_{\ell }}}{dt}\le$$

$$=$$

$$\frac{d\mathbf{q}}{dt}=D_2^{\prime }{C}_{-}^{\top }\mathbf{p}+(H^{\top }{B}^{\prime }-B^{\prime }-D_1^{\prime }-D_2^{\prime })\mathbf{q}-\mathbf{ϵ}.$$

$$\frac{d}{dt}\left[\begin{array}{c}\mathbf{p}\\ \mathbf{q}\end{array}\right]=\mathcal{A}\left[\begin{array}{c}\mathbf{p}\\ \mathbf{q}\end{array}\right]-\left[\begin{array}{c}0\\ \mathbf{ϵ}\end{array}\right].$$

$$\left[\begin{array}{c}\mathbf{p}(t)\\ \mathbf{q}(t)\end{array}\right]$$

$$\le e^{\mathcal{A}t}\left[\begin{array}{c}\mathbf{p}(0)\\ \mathbf{q}(0)\end{array}\right],$$

$$(B^{\top }-D)\mathbf{v}=-{\gamma }_1\mathbf{v}.$$

$$\mathbf{\xi }=\left[\begin{array}{c}\mathbf{v}\\ \mathbf{w}\end{array}\right],\mathbf{w}=C_{-}^{\top }\mathbf{v}={\left[v_{i_1},\dots ,v_{i_M}\right]}^{\top }.$$

$$\mathbf{\zeta }=\left[\begin{array}{c}{\mathbf{\zeta }}_1\\ {\mathbf{\zeta }}_2\end{array}\right],$$

$$B^{\top }=C_{+}{B}^{\prime }{C}_{-}^{\top }.$$

Full text

bristolbristolfootnotetext: Department of Engineering Mathematics, University of Bristol, Woodland Road, Clifton, Bristol BS8 1UB, UKubmathubmathfootnotetext: Department of Mathematics, University at Buffalo, State University of New York, Buffalo, NY 14260, USA ubcdseubcdsefootnotetext: Computational and Data-Enabled Science and Engineering Program, University at Buffalo, State University of New York, Buffalo, NY 14260, USApennpennfootnotetext: Department of Electrical and Systems Engineering, University of Pennsylvania, Philadelphia, PA 19104, USAApplMath_pennApplMath_pennfootnotetext: Applied Mathematics & Computational Science, University of Pennsylvania, Philadelphia, PA 19104, USA.naistnaistfootnotetext: Division of Information Science, Nara Institute of Science and Technology, Ikoma, Nara 630-0192, Japanosakaosakafootnotetext: Graduate School of Information Science and Technology, Osaka University, Suita, Osaka 565-0871, Japan

1 Introduction

Epidemic processes are probably one of the most extensively studied dynamical processes in complex networks [1, 2, 3, 4, 5]. These processes can be used for modeling the spread of infectious diseases in contact networks, as well as news in (offline or online) social networks, or computer viruses in communication networks, to name a few applications. A fundamental question in the analysis of epidemic processes, in the case of both deterministic and stochastic models, is to quantify the total number of nodes being infected by the spread over time. In most epidemic models, we find two clearly differentiated dynamical phases: one phase in which an initial infection quickly dies out and another phase in which an infection may propagate to a large fraction of the network. The concept of epidemic threshold is used for characterizing the conditions separating these two dynamical phases.

Most of the existing stochastic epidemic models are Markov processes where the disease-free state is a unique absorbing state. This absorbing state is reached with probability one in finite time, regardless of the initial set of infected nodes or the values chosen for the parameters of the model. A critical distinction between the two phases described above is the expected time required to reach the disease-free state. In the first phase mentioned above, the epidemic dynamics converges exponentially fast towards the absorbing state. In contrast, in the second phase, this time can be exponentially long in terms of the number of nodes. It is also worth remarking that this observation is not applicable to stochastic epidemic processes taking place in infinite networks [6, 7] or deterministic models taking place in both finite and infinite networks [8, 2, 4, 5] because in both these cases it is possible for the disease to survive forever. Therefore, for stochastic epidemic processes in finite networks, the exponential rate at which the number of infected individuals decays toward zero, called the decay rate, is a relevant characterization of the dynamics [9, 10, 11]. Intuitively, if the disease-free equilibrium takes a long time to be reached (in expectation), the decay rate would be close to zero. In contrast, if the infection dies out exponentially fast, the decay rate would be a positive value bounded away from zero. The decay rate can e.g. be used to measure the performance of control strategies aiming to eradicate an epidemic exponentially fast [12, 13, 14, 15].

Generally speaking, finding the decay rate of stochastic epidemic processes in a large network is computationally hard. This is because the number of possible states in the Markovian models typically used to model epidemics over networks grows exponentially in terms of the number of nodes in the network. Specifically, the decay rate is given by the leading eigenvalue of the transition-probability matrix of the Markovian model, whose dimension depends exponentially on the number of nodes. For example, in the case of the susceptible-infected-susceptible (SIS) model on $N$ nodes, the exponential rate corresponds to the leading eigenvalue of a $2^{N}\times 2^{N}$ transition-rate matrix [11], which is computationally challenging to calculate for large networks. An alternative approach to the exact computation of the decay rate is to seek computationally feasible bounds. For example, for the SIS model, a lower bound can be obtained using a mean-field approximation. This approximation is based on a first-order moment-closure technique allowing us to compute a bound on the decay rate from the leading eigenvalue of an $N\times N$ matrix [9, 13]. However, this mean-field approximation can result in a loose bound for many networks [16]. To increase the accuracy of the approximation, the authors proposed a tighter bound on the decay rate using a second-order moment-closure techniques [16]. This tighter bound, however, requires the computation of the leading eigenvalue of an $N^{2}\times N^{2}$ matrix, which can be computationally prohibitive when analyzing epidemic processes in large networks.

In the present work, we derive a new lower bound on the decay rate of the stochastic SIS model in an arbitrary finite network. This new bound depends on the leading eigenvalue of an $(N+M)\times(N+M)$ matrix, where $M$ is the number of directed edges; hence, for sparse networks — such as networks with a bounded maximum degree — the proposed lower bound is computationally more tractable than the bound derived in [16]. Our lower bound is based on an alternative second-order moment-closure technique aiming to overcome the computational challenges of existing second-order moment-closure techniques. The new bound depends on the non-backtracking matrix [17, 18]. The non-backtracking matrix has recently gained popularity in the network science community because it is the basis of efficient and theoretically appealing techniques for community detection, network centralities, and others (see references in [22]). We theoretically prove that the new lower bound is tighter than the first-order lower bound. We also show that our new lower bound is numerically more accurate than the bound obtained in [16]. We also present a new epidemic threshold, which corresponds to our lower bound on the decay rate. The new epidemic threshold is given in terms of the adjacency matrix of the line graph and the non-backtracking matrix of the given network.

2 Problem statement

We start with mathematical preliminaries. A directed graph is defined as the pair $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{v_{1},\dotsc,v_{N}\}$ is a finite ordered set of nodes, $N$ is the number of nodes, and $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ is a set of directed edges. By definition, $(v,v^{\prime})\in\mathcal{E}$ indicates that there is an edge from $v$ to $v^{\prime}$. The adjacency matrix of $\mathcal{G}$ is an $N\times N$ matrix of which the $(i,j)$-th entry is equal to $1$ if $(v_{i},v_{j})\in\mathcal{E}$ and [math] otherwise. An in-neighbor of $v$ is a node $v^{\prime}$ such that $(v^{\prime},v)\in\mathcal{E}$.

We denote the identity and the zero matrices by $I$ and $O$, respectively. A real matrix $A$ (or a vector as its special case) is said to be nonnegative, denoted by $A\geq 0$, if all the entries of $A$ are nonnegative. If all the entries of $A$ are positive, then $A$ is said to be positive. We say that $A\leq B$, where $A$ and $B$ are of the same dimension, whenever $B-A\geq 0$. A square matrix $A$ is said to be Metzler if all its off-diagonal entries are nonnegative [19]. If $A$ is Metzler, it holds true that $e^{At}\geq 0$ for all $t\geq 0$ [19]. For a Metzler matrix $A$, the maximum real part of the eigenvalues of $A$ is denoted by $\lambda_{\max}(A)$. For any matrix $A$, the spectral radius is the largest absolute value of its eigenvalues and denoted by $\rho(A)$.

We study the stochastic SIS model on networks, which is also known as the contact process in the probability theory literature [6]. This model is defined as follows: let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a directed graph. At any given continuous time $t\geq 0$, each node is in one of the two possible states, namely, susceptible (i.e. healthy) or infected. An infected node $v_{i}$ stochastically transits to the susceptible state at a constant instantaneous rate of $\delta_{i}>0$, which is called the recovery rate of node $v_{i}$. Whenever $v_{i}$ is susceptible and its infected in-neighbor $v_{j}$ is infected, then $v_{j}$ stochastically and independently infects $v_{i}$ at a constant instantaneous rate of $\beta_{ji}$. We call $\beta_{ji}>0$ the infection rate. Note that the present SIS model effectively accommodates directed and weighted networks because the infection rate $\beta_{ij}$ is allowed to depend on $v_{i}$ and $v_{j}$.

The SIS model is a continuous-time Markov process with $2^{N}$ possible states [11, 4, 5] and has a unique absorbing state in which all the $N$ nodes are susceptible. Because this absorbing state is reachable from any other state, the dynamics of the SIS model reaches the disease-free absorbing equilibrium in finite time with probability one. The aim of the present paper is to study how fast this disease-free equilibrium is reached in expectation. This can be quantified via the following definition:

Definition 1**.**

Let $p_{i}(t)$ be the probability that the $i$th node is infected at time $t$. The decay rate of the SIS model is defined by

[TABLE]

where all nodes are assumed to be infected at $t=0$.

Definition 1 states that $\sum_{i=1}^{N}p_{i}(t)$, which is equal to the expected number of infected nodes at time $t$, roughly decays exponentially in time as $\propto e^{-\gamma t}$. Because the number of infected nodes always becomes zero in finite time, the SIS model always has a positive decay rate (potentially close to zero), even if the infection rate is large.

The decay rate has theoretically been studied in continuous time [11] and discrete time [10] SIS models and is closely related to other quantities of interest, such as the epidemic threshold [11] and the mean time to absorption [9]. However, exact computation of the decay rate is computationally demanding in practice. Even in the homogeneous case, where all nodes share the same infection and recovery rates, the decay rate equals the modulus of the largest real part of the non-zero eigenvalues of a $2^{N}\times 2^{N}$ matrix representing the infinitesimal generator of the Markov chain [11].

Due to the difficulty of its computation, several approaches have been proposed to bound the decay rate. A first-order lower bound, which corresponds to the so-called quenched mean-field approximation [4], is derived as follows [9, 13]: Let $\bm{p}(t)=\left[p_{1}(t),\ldots,p_{N}(t)\right]^{\top}$, where $\top$ represents the matrix transposition. We define the $N\times N$ matrices $B$ and $D$ by

[TABLE]

and

[TABLE]

where $\text{diag}(\alpha_{1},\ldots,\alpha_{N})$ is the $N\times N$ diagonal matrix whose diagonal elements are equal to $\alpha_{1}$, $\ldots$, $\alpha_{N}$. Note that matrix $B$ fully contains the information about the adjacency matrix of $\mathcal{G}$. Then, one can show that $\bm{p}(t)\leq e^{(B^{\top}-D)t}\bm{p}(0)$, which implies that

[TABLE]

where we will call $\gamma_{1}$ the first-order lower bound. Although this lower bound is computationally efficient to find, there can be a large discrepancy between $\gamma_{1}$ and the true decay rate $\gamma$ [16].

A second lower bound on the decay rate was proposed in a recent study [16] and summarized in A. This second bound depends on the leading eigenvalue of an $N^{2}\times N^{2}$ matrix, which is computationally demanding when $N$ is relatively large. In this paper, we propose an alternative lower bound on the decay rate that is computationally efficient, provably more accurate than the first-order bound, and numerically tighter than the second bound described in A.

3 Main results

3.1 A lower bound on the decay rate

To state our mathematical results, we label the directed edges of a given network $\mathcal{G}$ as $\{e_{1},\ldots,e_{M}\}$, where the $\ell$th edge ($1\leq\ell\leq M$) is represented by $e_{\ell}=(i_{\ell},j_{\ell})$, i.e., the edge is directed from node $v_{i_{\ell}}$ to node $v_{j_{\ell}}$. Although we use notations $i$ and $j$ to represent general nodes $v_{i}$ and $v_{j}$ in the following text, $i$ and $j$ with a subscript will exclusively represent the starting and terminating nodes of an edge, thus avoiding confusions. Define the incidence matrix $C\in\mathbb{R}^{N\times M}$ of the network $\mathcal{G}$ by [20, 21]

[TABLE]

Also, define the non-backtracking matrix $H\in\mathbb{R}^{M\times M}$ of $\mathcal{G}$ by [17, 18]

[TABLE]

The main result of the present paper is stated as follows:

Theorem 2**.**

Define the $(N+M)\times(N+M)$ Metzler matrix

[TABLE]

where

[TABLE]

$C_{+}=\max(C,0)$, and $C_{-}=\max(-C,0)$; $C_{+}$ and $C_{-}$ denote the positive and negative parts of the incidence matrix $C$, respectively. Then, we obtain the following lower bound on the decay rate:

[TABLE]

Proof.

Define the binary variable $x_{i}(t)$ such that $x_{i}(t)=0$ or $x_{i}(t)=1$ if node $v_{i}$ is susceptible or infected at time $t$, respectively. The variables $x_{1}(t)$, $\ldots$, $x_{N}(t)$ obey a system of stochastic differential equations with Poisson jumps, and the expectation $p_{i}(t)=E[x_{i}(t)]$ obeys

[TABLE]

where

[TABLE]

is equal to the joint probability that node $v_{j}$ is infected and node $v_{i}$ is susceptible at time $t$.

Using the identities

[TABLE]

one obtains

[TABLE]

Equation (3.11) is equivalent to

[TABLE]

where we remind that $\bm{p}(t)=\left[p_{1}(t),\ldots,p_{N}(t)\right]^{\top}$ and define

[TABLE]

Using the notation $p_{ij}(t)\equiv E[x_{i}(t)x_{j}(t)]$, one obtains

[TABLE]

The first term on the right-hand side of the first line in Eq. (3.14) represents the rate at which node $v_{j}$ is infected when node $v_{i}$ is infected and node $v_{j}$ is susceptible; the second term represents the rate at which $v_{j}$ recovers when both $v_{i}$ and $v_{j}$ are infected; the third term represents the rate at which $v_{i}$ is infected when both $v_{i}$ and $v_{j}$ are susceptible; the fourth term represents the rate at which $i$ recovers when $v_{i}$ is infected and $v_{j}$ is susceptible. To derive the last inequality in Eq. (3.14), for the first term on the right-hand side, we ignored all the $k$ values but $k=i$ in the summation and used $x_{i}^{2}=x_{i}$. For the third term on the right-hand side, we used $E[(1-x_{i})(1-x_{j})x_{k}]\leq E[(1-x_{i})x_{k}]$.

By combining Eq. (3.14) and $p_{ij}=E[x_{i}x_{j}]=E[x_{i}]-E[x_{i}(1-x_{j})]=p_{i}-q_{ij}$, one obtains

[TABLE]

By combining Eq. (3.15) with

[TABLE]

and

[TABLE]

one obtains

[TABLE]

By stacking this inequality with respect to $\ell$, one observes that there exists an $\mathbb{R}^{M}_{+}$-valued function $\bm{\epsilon}(t)$ defined for $t\in[0,\infty)$ such that

[TABLE]

Equations (3.12) and (3.19) imply

[TABLE]

Because $\mathcal{A}$ is Metzler and $\bm{\epsilon}(t)$ is entry-wise nonnegative for every $t\geq 0$, one obtains

[TABLE]

which proves Eq. (3.7). ∎

Next, to prove that the new lower bound is tighter than the first-order lower bound, we start by stating (and proving) a convenient adaptation of the classical Perron–Frobenius theorem [25] for nonnegative matrices to the case of Metzler matrices.

Lemma 3**.**

Let $M$ be an irreducible Metzler matrix.

There exists a positive vector $\bm{v}$ such that $M\bm{v}=\lambda_{\max}(M)\bm{v}$. 2. 2.

Assume that there exist a real number $\mu$ and a nonzero vector $\bm{u}\geq 0$ such that $M\bm{u}\leq\mu\bm{u}$ and $M\bm{u}\neq\mu\bm{u}$. Then, $\lambda_{\max}(M)<\mu$.

Proof.

Let $\nu$ be a real number such that matrix $M^{\prime}=M+\nu I$ is nonnegative. Note that the spectral radius of $M^{\prime}$ satisfies $\rho(M^{\prime})=\lambda_{\max}(M)+\nu$. To prove the first statement, we use the Perron–Frobenius theorem (see Fact 5.b in [25, Chapter 9.2]), which guarantees that $M^{\prime}$ has a positive eigenvector $\bm{v}$ corresponding to the eigenvalue $\rho(M^{\prime})$. The vector $\bm{v}$ satisfies $M\bm{v}=M^{\prime}\bm{v}-\nu\bm{v}=\left[\rho(M^{\prime})-\nu\right]\bm{v}=\lambda_{\max}(M)\bm{v}$.

To prove the second statement, assume that a nonzero vector $\bm{u}\geq 0$ satisfies $M\bm{u}\leq\mu\bm{u}$ and $M\bm{u}\neq\mu\bm{u}$. Then, the nonnegative and irreducible matrix $M^{\prime}$ satisfies $M^{\prime}\bm{u}\leq(\mu+\nu)\bm{u}$ and $M^{\prime}\bm{u}\neq(\mu+\nu)\bm{u}$. Therefore, the Perron–Frobenius theorem (see Fact 7.b in [25, Chapter 9.2]) guarantees that $\rho(M^{\prime})<\mu+\nu$, which yields $\lambda_{\max}(M)=\rho(M^{\prime})-\nu<\mu$. ∎

The following theorem proves that the bound proposed in Eq. (3.7) improves the first-order bound given by Eq. (2.4).

Theorem 4**.**

If the network is strongly connected, then $\gamma_{2}>\gamma_{1}$.

Proof.

Lemma 3.1 implies that the irreducible Metzler matrix $B^{\top}-D$ has a positive eigenvector $\bm{v}$ corresponding to the eigenvalue $-\gamma_{1}$, i.e.,

[TABLE]

Define the positive $(N+M)$-dimensional vector $\bm{\xi}$ as

[TABLE]

Let us define $\bm{\zeta}\equiv\mathcal{A}\bm{\xi}$ and decompose $\bm{\zeta}$ as

[TABLE]

where $\bm{\zeta}_{1}$ and $\bm{\zeta}_{2}$ are $N$- and $M$-dimensional vectors, respectively. Simple algebraic manipulations yield

[TABLE]

Therefore, one obtains $C_{+}B^{\prime}\bm{w}=C_{+}B^{\prime}C_{-}^{\top}\bm{v}=B^{\top}\bm{v}$. This implies that

[TABLE]

One also obtains

[TABLE]

where $A_{L(\mathcal{G})}$ denotes the adjacency matrix of the line graph $L(\mathcal{G})$ defined by

[TABLE]

Matrix $A_{L(\mathcal{G})}$ satisfies

[TABLE]

because

[TABLE]

Because simple algebraic manipulations yield $C_{-}^{\top}D=D_{1}^{\prime}C_{-}^{\top}$, using Eqs. (3.22), (3.25), and (3.29), one obtains

[TABLE]

Using Eqs. (3.27) and (3.31), one obtains

[TABLE]

Because any entry of $B^{\prime}\bm{w}$ is positive, Eqs. (3.26) and (3.32) guarantee that positive vector $\bm{\xi}$ satisfies $\mathcal{A}\bm{\xi}\leq-\gamma_{1}\bm{\xi}$ and $\mathcal{A}\bm{\xi}\neq-\gamma_{1}\bm{\xi}$. Because $\mathcal{A}$ is irreducible, as will be shown later, Lemma 3 guarantees that $\lambda_{\max}(\mathcal{A})<-\gamma_{1}$, which implies that $\gamma_{2}>\gamma_{1}$.

Finally, let us show the irreducibility of matrix $\mathcal{A}$ or, equivalently, the irreducibility of $\mathcal{A}^{\top}$. We regard the matrix $\mathcal{A}^{\top}$ as the adjacency matrix of a directed graph on $N+M$ nodes denoted by $\mathcal{G}^{\prime}$. We label the nodes of $\mathcal{G}^{\prime}$ as $p_{1}$, $\ldots$, $p_{N}$, $q_{i_{1}j_{1}}$, $\ldots$, $q_{i_{M}j_{M}}$. The first term on the right-hand side of Eq. (3.8) implies that $\mathcal{G}^{\prime}$ has an edge $(q_{i_{\ell}j_{\ell}},p_{j_{\ell}})$ for all $\ell$, which corresponds to $C_{+}B^{\prime}$ in Eq. (3.3). The second term on the right-hand side of Eq. (3.15) implies that $\mathcal{G}^{\prime}$ has an edge $(p_{i_{\ell}},q_{i_{\ell}j_{\ell}})$ for all $\ell$, which corresponds to $D_{2}^{\prime}C_{-}^{\top}$ in Eq. (3.3).

To show that $\mathcal{G}^{\prime}$ is strongly connected, we first consider an arbitrary ordered pair of nodes $p_{i}$ and $p_{j}$ in $\mathcal{G}^{\prime}$. Let us take a path $v_{i}=v_{\iota(0)}$, $v_{\iota(1)}$, …, $v_{\iota(s)}=v_{j}$ in the original graph $\mathcal{G}$. Then, from the above observation, we see that the graph $\mathcal{G}^{\prime}$ contains the path $p_{i}=p_{\iota(0)}$, $q_{\iota(0)\iota(1)}$, $q_{\iota(1)}$, $q_{\iota(1)\iota(2)}$, …, $q_{\iota(s-1)\iota(s)}$, $p_{\iota(s)}=p_{j}$. Likewise, for an arbitrary ordered pair of nodes $p_{i}$ and $q_{i_{\ell}j_{\ell}}$ in $\mathcal{G}^{\prime}$, there is a path in $\mathcal{G}^{\prime}$ from $p_{i}$ to $p_{i_{\ell}}$. By appending edge $(p_{i_{\ell}},q_{i_{\ell}j_{\ell}})$ to the end of this path, one obtains a path from $p_{i}$ to $q_{i_{\ell}j_{\ell}}$. A path from arbitrary $q_{i_{\ell}j_{\ell}}$ to $p_{j}$ and one from $q_{i_{\ell}j_{\ell}}$ to $q_{i_{\ell^{\prime}}j_{\ell^{\prime}}}$ can be similarly constructed. Therefore, a path exists between any pair of nodes in $\mathcal{G}^{\prime}$. ∎

3.2 Epidemic threshold

In this section, we assume that $\beta_{ij}=\beta$ and $\delta_{i}=\delta$, where $i,j\in\{1,\ldots,N\}$ and $\beta,\delta>0$, and derive conditions under which the expected number of infected individuals decays exponentially fast. It holds true that having $\gamma_{1}<0$ in Eq. (2.4) is equivalent to the well-known epidemic threshold $\beta/\delta>1/\lambda_{\max}(A)$ [26, 10, 4, 5].

Likewise, Theorem 2 provides a tighter epidemic threshold given by

[TABLE]

where $\gamma_{2}$ is defined in Eq. (3.7). Our following corollary provides an explicit expression of the epidemic threshold in terms of the adjacency matrix of the line graph and the non-backtracking matrix $H$.

Corollary 5**.**

Define the matrix $A_{L(\mathcal{G})}$ by Eq. (3.28). Then,

[TABLE]

Proof.

We decompose $\mathcal{A}$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

Matrix $R$ is a Metzler matrix, all the eigenvalues of $R$ have negative real parts, and matrix $P$ is nonnegative. Therefore, Theorem 2.11 in Ref. [27] implies that $\lambda_{\max}(\mathcal{A})<0$ if and only if $\rho(R^{-1}P)<1$. Because

[TABLE]

one obtains

[TABLE]

where we used Eq. (3.29). Therefore, $\rho(R^{-1}P)<1$ if and only if

[TABLE]

which is equivalent to Eq. (3.34). ∎

Remark: Corollary 5 does not require strong connectedness (i.e., irreducibility of the adjacency matrix) of the network.

4 Numerical results

In this section, we carry out numerical simulations of the stochastic SIS dynamics for several networks to assess the tightness of the different lower bounds on the decay rate. In the following numerical simulations, we set $\beta_{ij}=\beta$ and $\delta_{i}=\delta$, where $i,j\in\{1,\ldots,N\}$, for simplicity. We further assume that $\delta=1$ without loss of generality (because changing $\beta$ and $\delta$ simultaneously by the same factor is equivalent to rescaling the time variable without changing $\beta$ or $\delta$).

We ran the stochastic SIS dynamics $10^{4}$ times starting from the initial condition in which all nodes are infected. For each run of the simulations, we measured the number of infected individuals at every integer time (including time 0) until the infection dies out or the maximum time, which is set to $5\times 10^{4}$, is reached. Then, at each integer time, we summed the number of infected nodes over all the runs and divided it by $N$ and by the number of runs ($=10^{4}$), thus obtaining the average fraction of infected nodes, i.e., $\rho(t)\equiv\sum_{i=1}^{N}p_{i}(t)/N$, where $t=0,1,\ldots$.

We calculated the decay rate from the observed $\{\rho(t):t=0,1,\ldots\}$ as follows: because the fluctuations in $\rho(t)$ are expected to be large when $\rho(t)$ is small, we identified the smallest integer time at which $\rho(t)$ is less than $10^{-4}$ for the first time, and discarded $\rho(t)$ at this and all larger $t$ values. Then, because $\rho(t)$ is expected to decay exponentially in $t$, we calculated a linear regression between $\log\rho(t)$ and $t$ at the remaining integer values of $t$. The sign-flipped slope of this regression provides a numerical estimate of the decay rate. We confirmed that the Pearson correlation coefficient in the linear regression was at least $0.958$ for all networks and all $\beta$ values. The Pearson correlation was typically larger than $0.99$.

We used eight undirected and unweighted networks to compare the numerically obtained decay rate and the rigorous lower bounds. The lower bounds to be compared are $\gamma_{1}$, $\gamma_{2}$, and the one obtained in our previous study [16], which is denoted by $\gamma_{2}^{\prime}$ (see A for a summary).

Four of the eight networks used were created by generative models with $N=100$ nodes. First, we generated a regular random graph in which all nodes had degree six, resulting in $300$ undirected edges (therefore, $M=600$ directed edges). Second, we used the Barabási-Albert (BA) model to generate a power-law degree distribution with an exponent of 3 when $N$ is large [28]. We set the parameters $m_{0}=3$ and $m=3$, where $m_{0}$ is the initial number of nodes forming a clique in the process of growing a network, and $m$ is the number of edges that each new node initially brings into the network. With these parameter values, the mean degree is approximately equal to $2m=6$. The generated network had $294$ undirected edges. Third, we used a cycle graph, where each node had degree two (by definition), and there were $100$ undirected edges. These three models lack community structure that many empirical contact networks have. Therefore, as a fourth network, we used the Lancichinetti–Fortunato–Radicchi (LFR) model that can generate networks with community structure [29]. The LFR model creates networks having a heterogeneous degree distribution and a heterogeneous distribution of community size. A small value of parameter $\mu$ corresponds to a strong community structure. We set $\mu=0.1$. We set the mean degree to six, the largest degree to $N/4=25$, the power-law exponent for the degree distribution to two and the power-law exponent for the distribution of community size to one. The network had $319$ undirected edges.

We also used four real-world networks, for which we ignored the direction and weight of the edge. First, we used the dolphin social network, which has $N=62$ nodes and $159$ undirected edges [30]. A node represents a bottleneck dolphin individual. An edge indicates frequent association between two dolphins. This network is a connected network. Second, we used the largest connected component (LCC) of a coauthorship network of researchers in network science, which has $N=379$ nodes and $914$ undirected edges [31]. A node represents a researcher publishing in fields related to network science. An edge indicates that two researchers have coauthored a paper at least once. Third, we used the LCC of an email network, which has $N=1,133$ nodes and $5,451$ undirected edges [32]. A node represents a member of the University Rovira i Virgili, Tarragona, Spain. An edge is an email exchange relationship between a pair of members. Fourth, we used the LCC of the hamsterter network, which has $N=1,788$ nodes and $12,476$ undirected edges [33]. A node represents a user of the website hamsterster.com. An edge is a friendship relationship between two users.

For a range of values of $\beta$, we compare decay rates obtained numerically with the three lower bounds described in this paper for the eight networks mentioned above. The results are shown in Fig. 1. It should be noted that the decay rate and its bounds are equal to one for $\beta=0$ because we set $\delta=1$ for normalization. The bound $\gamma_{2}$ proposed in the present study is considerably tighter than the first-order bound, $\gamma_{1}$, for some networks, in particular, the cycle (Fig. 1(c)). The improvement tends to be more manifested for smaller networks. We also find that $\gamma_{2}$ is tighter than $\gamma_{2}^{\prime}$ for all the networks and infection rates, despite that $\gamma_{2}$ is easier to calculate than $\gamma_{2}^{\prime}$. For example, for the regular random graph (Fig. 1(a)) and the cycle (Fig. 1(c)), $\gamma_{2}$ is close to the numerically estimated decay rate for small to moderate values of $\beta$, which is not the case for $\gamma_{2}^{\prime}$ as well as for $\gamma_{1}$.

5 Conclusions

We have introduced a lower bound on the decay rate of the SIS model on arbitrary directed and weighted networks. The new bound is based on a new second-order moment-closure technique aiming to improve both the computational cost and the accuracy of existing second-order bound. It is equal to the leading eigenvalue of an $(N+M)\times(N+M)$ Metzler matrix depending on the non-backtracking and incidence matrices of the network (Eq. (3.3)). Therefore, for sparse networks, the dimension of this matrix grows quasi-linearly. Furthermore, we have shown that the new bound, $\gamma_{2}$, is tighter than the first-order lower bound, $\gamma_{1}$, which is equal to the leading eigenvalue of an $N\times N$ matrix depending directly on the adjacency matrix.

Non-backtracking matrices of networks have been employed for analyzing properties of stochastic epidemic processes on networks, such as the epidemic threshold of the SIS model [34, 35] and the susceptible-infected-recovered (SIR) model [36, 37, 38, 39]. The non-backtracking matrix more accurately describes unidirectional state-transition dynamics, such as the SIR dynamics, than the adjacency matrix does because unidirectional dynamics implies that contagions do not backtrack, i.e. if node $v_{i}$ has infected its neighbor $v_{j}$, $v_{j}$ does not re-infect $v_{i}$. For the same reason, the non-backtracking matrix also predicts the percolation threshold for networks better than the adjacency matrix [40, 41]. However, the same logic does not apply to the SIS model, in which re-infection through the same edge can happen indefinitely many times. This is a basis of a recent criticism to the application of the non-backtracking matrix to the SIS model [42]. For some networks, the epidemic threshold of the SIS model that does not take into account backtracking infection paths [34, 35] is not accurate [42]. Although $\gamma_{2}$ and the corresponding epidemic threshold that we have derived use the non-backtracking matrix, they are mathematical bounds and do not suffer from the inaccuracy caused by the neglect of backtracking infection paths. The present study has shown a new and solid usage of the non-backtracking matrix in understanding the SIS model on networks.

By following the derivation of the epidemic threshold via $\gamma_{1}$, we derived the epidemic threshold based on $\gamma_{2}$. The new epidemic threshold is always larger than that based on $\gamma_{1}$, which is the reciprocal of the largest eigenvalue of the adjacency matrix. Because $\gamma_{2}$ improves upon $\gamma_{1}$, we expect that the new epidemic threshold is a better estimate than that based on $\gamma_{1}$. This point warrants future work. Likewise, the eigenvalue statistics for the adjacency matrix of scale-free networks yield intricate relationships between the epidemic threshold based on $\gamma_{1}$ and statistics of the node’s degree in scale-free networks [43]. How such a result translates to the case of the epidemic threshold based on $\gamma_{2}$ also warrants future work.

Appendix A Lower bound on the decay rate derived in Ref. [16]

In the proof of Theorem 2, we have used the following inequality for bounding $q_{ij}=E[x_{i}(1-x_{j})]$ (see Eq. (3.14)):

[TABLE]

in which the inequality $x_{j}\geq 0$ is used; we have presumed that node $j$ is susceptible. In contrast, in our previous study [16], we used

[TABLE]

which was based on  $x_{i}\geq 0$. The use of Eq. (A.2) led to the following lower bound on the decay rate [16]:

Theorem 6**.**

Assume that there exist positive numbers $\beta_{1}$, …, $\beta_{N}$ such that $\beta_{ij}=\beta_{j}$ for all $i,j\in\{1,\ldots,N\}$. Let $A$ be the adjacency matrix of $\mathcal{G}$ and its $(i,j)$th entry be $a_{ij}$. Define the $N^{2}\times N^{2}$ Metzler matrix

[TABLE]

where $\bigoplus_{i=1}^{N}M_{i}$ is the block-diagonal matrix containing matrices $M_{1}$, $\ldots$, $M_{N}$ as the diagonal blocks, $\backslash\{i\}$ denotes all the columns except the $i$th column, $V_{i}\in\mathbb{R}^{(N-1)\times N}$ is the matrix obtained by removing the $i$th row from the $N\times N$ identity matrix, $\Gamma_{i}=\text{diag}(\overline{\gamma}_{i,1},\ldots,\overline{\gamma}_{i,i-1},\overline{\gamma}_{i,i+1},\ldots,\overline{\gamma}_{i,N})$, and $\overline{\gamma}_{i,j}=\delta_{i}+\delta_{j}+a_{ij}\beta_{i}$. Then, the decay rate $\gamma$ satisfies

[TABLE]

Acknowlegdments

We thank Claudio Castellano for valuable comments on the manuscript.

Funding

National Science Foundation (CAREER-ECCS-1651433 to V.M.P.) and Japan Society for the Promotion of Science (JP18K13777 to M.O.).