On the Dynamics of Deterministic Epidemic Propagation over Networks

TL;DR

This paper reviews deterministic nonlinear models for infectious disease spread over contact networks, analyzing equilibria, stability, and thresholds for SI, SIS, and SIR models, with new results on endemic states and transient behavior.

Contribution

It provides novel analytical results and algorithms for the stability, endemic states, and transient dynamics of network-based epidemic models.

Findings

Established equilibria and stability for SI model

Characterized endemic states and thresholds for SIS model

Proposed an iterative algorithm for SIR asymptotic states

Abstract

In this work we review a class of deterministic nonlinear models for the propagation of infectious diseases over contact networks with strongly-connected topologies. We consider network models for susceptible-infected (SI), susceptible-infected-susceptible (SIS), and susceptible-infected-recovered (SIR) settings. In each setting, we provide a comprehensive nonlinear analysis of equilibria, stability properties, convergence, monotonicity, positivity, and threshold conditions. For the network SI setting, specific contributions include establishing its equilibria, stability, and positivity properties. For the network SIS setting, we review a well-known deterministic model, provide novel results on the computation and characterization of the endemic state (when the system is above the epidemic threshold), and present alternative proofs for some of its properties. Finally, for the network SIR setting, we propose novel results for transient behavior, threshold conditions, stability properties, and asymptotic convergence. These results are analogous to those well-known for the scalar case. In addition, we provide a novel iterative algorithm to compute the asymptotic state of the network SIR system.

Figures

Equations

$$x\ll y,$$

$$x\le y,$$

$$x<y,$$

$$v_{\max }^{\top }A={\lambda }_{\text{max}}{v}_{\max }^{\top }\text{and}A{u}_{\max }={\lambda }_{\text{max}}{u}_{\max },$$

$$\dot{x}(t)=\beta s(t)x(t)=\beta\big{(}1-x(t)\big{)}x(t),$$

$$x(t)=\frac{x_0{\mathrm{e}}^{\beta t}}{1-x_0+x_0{\mathrm{e}}^{\beta t}}.$$

$$\dot{x}_{i}(t)=\beta\big{(}1-x_{i}(t)\big{)}\sum_{j=1}^{n}a_{ij}x_{j}(t),$$

$$\dot{x}(t)=\beta\Big{(}I_{n}-\operatorname{diag}\!\big{(}x(t)\big{)}\Big{)}Ax(t),$$

$$\dot{s}(t)=-\beta \mathrm{diag}(s(t))A(\text{mathbbold}{1}_n-s(t)).$$

$$\displaystyle x(t)=\mathrm{e}^{\beta At}x(0)=T\mathrm{e}^{\beta Jt}T^{-1}x(0)=\mathrm{e}^{\beta\lambda_{\textup{max}}t}\big{(}T\mathbbold{e}_{1}\mathbbold{e}_{1}^{\top}T^{-1}x(0)+o(1)\big{)},$$

$$x(t)$$

$$\displaystyle=\mathrm{e}^{\beta\lambda_{\textup{max}}t}\Big{(}\frac{v_{\max}^{\top}x(0)}{v_{\max}^{\top}u_{\max}}u_{\max}+o(1)\Big{)}.$$

$${\dot{s}}_i(t)=-\beta d_i{s}_i(t)⟹s_i(t)={ϵ}_i{\mathrm{e}}^{-\beta d_i(t-T)}.$$

$$\dot{x}(t)=\beta s(t)x(t)-\gamma x(t)=(\beta -\gamma -\beta x(t))x(t).$$

$$x(t)=\frac{(\beta -\gamma )x_0}{\beta x_0-{\mathrm{e}}^{-(\beta -\gamma )t}(\gamma -\beta (1-x_0))};$$

$${\dot{x}}_i(t)=\beta (1-x_i(t))\sum _{j=1}^n{a}_{ij}{x}_j(t)-\gamma x_i(t),$$

$$\dot{x}(t)=\beta\big{(}I_{n}-\operatorname{diag}(x(t))\big{)}Ax(t)-\gamma x(t).$$

$$f_{+}(y)=y/(1+y),\text{and}{f}_{-}(z)=z/(1-z),$$

$$\dot{y}(t)$$

$$a=\frac{v_{\max }^T{u}_{\max }}{v_{\max }^T\mathrm{diag}(u_{\max })u_{\max }},$$

$$y(k+1):=F_{+}\Bigg{(}\frac{\beta}{\gamma}Ay(k)\Bigg{)}.$$

$$\underset{k\to \infty }{\lim }y(k)=x^{\ast }.$$

$$F_{+}(\hat{A}y(0))_{i}=f_{+}\Big{(}\frac{\beta\lambda_{\textup{max}}}{\gamma}y_{i}(0)\Big{)}\geq y_{i}(0).$$

$$\displaystyle\big{(}$$

$$\displaystyle\geq f_{+}\big{(}\alpha(\hat{A}x^{*})_{i}\big{)}-\alpha x^{*}_{i}\qquad\qquad\qquad\text{($\hat{A}\geq\mathbbold{0}_{n\times n}$)}$$

$$\displaystyle>\alpha f_{+}\big{(}(\hat{A}x^{*})_{i}\big{)}-\alpha x^{*}_{i}\qquad\text{($0<\alpha<1$ and $z>0$)}$$

$$\displaystyle=\alpha\big{(}F_{+}(\hat{A}x^{*})-x^{*}\big{)}_{i}=0.\qquad\text{($x^{*}$ is an equilibrium)}$$

$$y(0)=\left(1-\frac{\gamma }{\beta {\lambda }_{\text{max}}}\right)\frac{u_{\max }}{{\max }_i\{{u_{\max }}_{,i}\}}=\frac{\delta }{\delta +1}\frac{u_{\max }}{{\max }_i\{{u_{\max }}_{,i}\}}$$

$$y(0)=\left(1-\frac{\gamma }{\beta {\lambda }_{\text{max}}}\right)\frac{u_{\max }}{{\min }_i\{{u_{\max }}_{,i}\}}=\frac{\delta }{\delta +1}\frac{u_{\max }}{{\min }_i\{{u_{\max }}_{,i}\}}$$

$$\frac{\delta }{\delta +1}\frac{u_{\max }}{{\max }_i\{{u_{\max }}_{,i}\}}\le x^{\ast }\le \frac{\delta }{\delta +1}\frac{u_{\max }}{{\min }_i\{{u_{\max }}_{,i}\}}$$

$$(\delta+1)\big{(}I_{n}-\operatorname{diag}(x^{*})\big{)}Ax^{*}-\lambda_{\textup{max}}x^{*}=0$$

$$A{x}_1-{\lambda }_{\text{max}}{x}_1=0$$

$$A{x}_1+A{x}_2-\mathrm{diag}(x_1)A{x}_1-{\lambda }_{\text{max}}{x}_2=0$$

Full text

Index terms: network propagation model, nonlinear dynamical system, phase transition, mathematical epidemiology

1 Introduction

1.1 Motivation and problem description

Propagation phenomena appear in numerous disciplines. Examples include the spread of infectious diseases in contact networks, the transmission of information in communication networks, the diffusion of innovations in competitive economic networks, cascading failures in power grids, and the spreading of wild-fires in forests. Scalar models of propagation phenomena have been widely studied, e.g., see the survey by Hethcote [12] on scalar epidemic spreading models. These models qualitatively capture some dynamic features, including phase transitions and asymptotic states. However, shortcomings of these scalar models are also prominent: for example, scalar models are typically based on the assumption that individuals in the population have the same chances of interacting with each other. This assumption overlooks the internal structure of the network over which the propagation occurs, as well as the heterogeneity of individuals in the network. Both these aspects play critical roles in shaping the sophisticated dynamical behavior of the propagation processes.

In a general formulation, propagation is a stochastic process on a complex network. Its features are determined by the properties of local node-to-node exchanges as well as of the global contact networks. Stochastic propagation processes can be modeled as Markov chains of exponential dimensions in the size of the network. Due to these properties of propagation processes, many relevant research questions arise naturally. Three key questions are: (1) how to deal with the lack of well-quantified local details of the complex networks; (2) how to identify the models and parameters, and estimate the states of such stochastic large-dimension phenomena, and (3) how to analyze the transient and asymptotic properties of their dynamics.

Three approaches are commonly adopted in the analysis of propagation models. The first is to directly characterize the stochastic process adopting tools from branching processes, percolation, and random graph theory; e.g., see [7]. The second approach proposes a degree-based model based on the approximation that nodes with the same degree exhibit similar behavior; e.g., see [17, Chapter 17]. The first two approaches are both based on the random graphs for which global characteristics can be estimated statistically, such as the degree distribution. The third approach is based on the mean-field approximation of Markov-chain models and algebraic graph theory; this approach results in a deterministic network model. Distinct from the first two approaches, the mean-field models first assume knowledge of the local propagation parameters and then relate the dynamical behavior of the propagation process to some global parameters of the network, which can be estimated without knowing the full local details, e.g., the spectral radius of the network’s adjacency matrix. One of the main advantages of the mean-field approach is that, by representing the network as an adjacency matrix, well-established theorems in matrix analysis and dynamical systems can be applied to the analysis of some sophisticated behavior of the propagation processes.

In this paper we review a class of epidemic propagation models which adopt the third approach. Distinct in the assumptions on the microscopic features of the disease and the individual behavior, the epidemic propagation models we focus on are classified into three types: the Susceptible-Infected (SI) model, the Susceptible-Infected-Susceptible (SIS) model and the Susceptible-Infected-Recovered (SIR) model; basic representations of these models are illustrated in Figure 1. In this work we review models of these three types over networks and characterize their dynamical properties. In short, we study the network SI, network SIS, and network SIR models.

1.2 Literature review

The dynamics of several classic scalar epidemic models, i.e., the population models without network structure, are surveyed in detail by Hethcote [12]. Among the different metrics discussed, identifying the effective reproduction number $R$ is of particular interest to researchers; $R$ is the expected number of individuals that a randomly infected individual can infect during its infection period. In these scalar models, whether an epidemic outbreak occurs or the disease dies down depends upon whether $R>1$ or $R<1$, i.e., upon whether the system is above or below the so-called epidemic threshold. Here by epidemic outbreak we mean an exponential growth of the fraction of the infected population for small time. The basic reproduction number $R_{0}$ is the effective reproduction number in a fully-healthy susceptible population. In what follows we focus our review on deterministic network models.

The earliest work on the (continuous-time heterogeneous) SIS model on networks is [14]. This work proposes an $n$-dimensional model on a contact network and analyzes the system’s asymptotic behavior. This article proposes a rigorous analysis of the threshold for the epidemic outbreak, which depends on both the disease parameters and the spectral radius of the contact network. For the case when the basic reproduction number is above the epidemic threshold, this paper establishes the existence and uniqueness of a nonzero steady-state infection probability, called the endemic state. In what follows we refer to the model proposed by Lajmanovich et al.  [14] as the network SIS model; it is also known as the multi-group or multi-population SIS model. Numerous extensions and variations on these basic results have appeared over the years.

Allen [2] proposes and analyzes a discrete-time network SIS model. This work appears to be the first to revisit and formally reproduce, for the discrete-time case, the earlier results by Lajmanovich et al. [14]; see also the later work by Wang et al. [21]. This work confirms the existence of an epidemic threshold, as a function of the spectral radius of the contact network. Further recent results on the discrete-time model are obtained by Ahn et al. [1] and by Azizan Ruhi et al. [3].

Van Mieghem et al.[20] argue that the (continuous-time) network SIS model is in fact the mean-field approximation of the original Markov-chain SIS model of exponential dimension; this claim is rigorously proven by Sahneh et al. [19]. Van Mieghem et al. [20] refer to this model as the intertwined SIS model and write the endemic state as a continued fraction.

The works by Fall et al. [9] and Khanafer et al. [13] discuss the continuous-time network SIS model in a more modern language. Fall et al. [9] refer to this model as the $n$-group SIS model and apply Lyapunov techniques and Metzler matrix theory to establish existence, uniqueness, and stability of the equilibrium points below and above the epidemic threshold. Khanafer et al. [13] use positive system theory in their analysis and extend the existence, uniqueness, and stability results to the setting of weakly connected digraphs.

An early work by Hethcote [11] proposes a general multi-group SIR model with birth, death, immunization, and de-immunization. The epidemic threshold and the equilibria below/above the threshold are characterized. For the simplified model without birth/death and de-immunization, Hethcote [11] proves that the system converges asymptotically to an all-healthy state. Guo et al. [10] consider a generalized network SIR model with vital dynamics, that is, with birth and death. They characterize the basic reproduction number and, through a careful Lyapunov analysis, show the existence and global asymptotic stability of an endemic state above the threshold. Youssef et al. [22] study a special case of the network SIR model under the name of individual-based SIR model over undirected networks. Through a simulation-based analysis, the epidemic threshold is given as a function of the spectral radius of the network. To the best of our knowledge, no works have comprehensively characterized the properties of the network SI model.

We conclude by mentioning other surveys and textbook treatments. In [15], the stability of equilibria for the SEIR model is reviewed through Lyapunov and graph theory. The additional state $E$ represents the exposed population, i.e., the individuals who are infected but not infectious. The book chapters [17, Chapter 17], [8, Chapter 21], and [4, Chapter 9] review various heterogeneous epidemic models. The recent survey by Nowzari et al. [18] presents various epidemic models and addresses many solved and open problems in the control of epidemic spreading.

1.3 Statement of Contribution

The contributions of this work are as follows: in each section, we start by reviewing the scalar SI, SIS, and SIR models; these are the models in which variables represent an entire “well-mixed” population or nodes of an all-to-all unweighted graph. We then focus our discussion on multi-group network models and provide a tutorial comprehensive treatment with comprehensive statements and proofs for the network SI, SIS and SIR models.

We first introduce and analyze the novel network SI model. We analyze its asymptotic convergence, positivity of infection probabilities, initial growth rate, and the stability of equilibria. We show that in the network SI model, the system does not display a threshold and all the trajectories converge to the full contagion state.

Next we focus on the network SIS model. We review some results in [14, 9, 13] regarding the dynamical behavior of the system below and above the threshold, and present alternative proofs for them. For systems above the epidemic threshold, we present a novel provably-correct iterative algorithm for computing the fraction of infected individuals converging to the endemic state. We present novel Taylor expansions for the endemic state near the epidemic threshold and in the limit of high infection rates. Finally, we show that the spread of infection takes place instantaneously upon infecting at least one node in the network.

Finally, for the network SIR model, we present novel transient behavior and system properties. We propose new threshold conditions above which the epidemic grows initially, and below which it exponentially dies down. We show that, along all system trajectories, the infected population asymptotically vanishes and the epidemic asymptotically dies down. The initial rate of growth above the threshold is given in terms of network characteristics, initial conditions, and infection parameters. We show that our proposed weighted average of the infected population, obtained by the entries of dominant eigenvector of an irreducible quasi-positive matrix, captures information regarding the distribution of infection in the system. We also establish positivity of the infection probabilities and certain monotonicity properties. Moreover, we provide a novel iterative algorithm to compute the asymptotic state of the network SIR model, with any arbitrary initial condition. For the iterative algorithm, the existence and uniqueness of the fixed point, and the convergence of the iteration are rigorously proved. Our results are analogous to the scalar SIR model properties and are valid for any arbitrary network topologies. In comparison with [22], our treatment builds on their numerical results but our result is more general in that it does not depend upon specific initial conditions and graph topologies, and establishes numerous properties, including the novel characterization of epidemic threshold.

Finally, we remark that our deterministic network models are derived from the Markov-chain models through a mean-field approximation. We do not discuss here the Markov-chain model and the approximation process and refer instead to [19] and [6, Chapter 17].

1.4 Organization

Section 2 introduces our model set-up and some preliminary notations. The SI, SIS and SIR models are presented, respectively, in Sections 3, 4, and 5. Section 6 is the conclusion.

2 Model Set-Up and Notations

For the scalar models, we use the notation $x(t)$ ($s(t)$ and $r(t)$ resp.) for the fraction of infected (susceptible and recovered resp.) individuals in the population at time $t$. The rest of this section is about the notations and basic model set-up for the network epidemic model.

a) Contact Network: The epidemics are assumed to propagate over a weighted digraph $G=(V,E)$, where $V=\{1,\dots,n\}$ and $E$ is the set of directed links. Nodes of $G$ can be interpreted as either single individuals in the contact network or as homogeneous populations of individuals at each location/node in the contact network. $A=(a_{ij})_{n\times n}$ denotes the adjacency matrix associated with $G$. For any $i$, $j\in V$, $a_{ij}$ characterizes the contact strength from node $j$ to node $i$. For $(i,j)\in E$, $a_{ij}>0$ and for $(i,j)\notin E$, $a_{ij}=0$. In this paper, $G$ is assumed to be strongly connected.

b) Node States and Probabilities: For different epidemic propagation models, the set of possible node states are distinct. For network SI or SIS models, each node can be in either the “susceptible” or “infected” state, while in the network SIR model, there is an additional possible node state: “recovered.” For a graph in which the nodes are single individuals, let $s_{i}(t)$ ($x_{i}(t)$ and $r_{i}(t)$ resp.) be the probability that individual $i$ is in the susceptible (infected and recovered resp.) state at time $t$. Alternatively, if the nodes are considered to be the populations, then $s_{i}(t)$ ($x_{i}(t)$ and $r_{i}(t)$ resp.) is interpreted as the fraction of susceptible (infected and recovered resp.) individuals in population $i$. In this paper, without loss of generality, we adopt the interpretation of nodes as single individuals.

c) Frequently Used Notations: The symbol $\mathbb{R}$ denotes the set of real numbers, while $\mathbb{R}_{\geq 0}$ denotes the set of non-negative real numbers. The symbol $\phi$ denotes the empty set. For any two vectors $x,y\in\mathbb{R}^{n}$, we write

[TABLE]

We adopt the shorthand notations $\mathbbold{1}_{n}=[1,\dots,1]^{\top}$ and $\mathbbold{0}_{n}=[0,\dots,0]^{\top}$. Given $x=[x_{1},\dots,x_{n}]^{\top}\in\mathbb{R}^{n}$, let $\operatorname{diag}(x)$ denote the diagonal matrix whose diagonal entries are $x_{1},\dots,x_{n}$. For an irreducible nonnegative matrix $A$, let $\lambda_{\textup{max}}(A)$ denote the dominant eigenvalue of $A$ that is equal to the spectral radius $\rho(A)$. Moreover, we let $v_{\max}(A)$ ($u_{\max}(A)$ resp.) denote the corresponding entry-wise strictly positive left (right resp.) eigenvector associated with $\lambda_{\textup{max}}(A)$, normalized to satisfy $\mathbbold{1}_{n}^{\top}v_{\max}(A)=1$ (resp. $\mathbbold{1}_{n}^{\top}u_{\max}(A)=1$). The Perron-Frobenius Theorem for irreducible matrices guarantees that $\lambda_{\textup{max}}(A)$, $v_{\max}(A)$ and $u_{\max}$ are well defined and unique. Where not ambiguous, we will drop the $(A)$ argument and, for example, write

[TABLE]

with $v_{\max}\gg\mathbbold{0}_{n}$ and $\mathbbold{1}_{n}^{\top}v_{\max}=1$; $u_{\max}\gg\mathbbold{0}_{n}$ and $\mathbbold{1}_{n}^{\top}u_{\max}=1$.

3 Susceptible-Infected Model

In this section, we first review the classic scalar susceptible-infected (SI) model, and then present and characterize the network SI model.

3.1 Scalar SI model

The scalar SI model assumes that the growth rate of the fraction of the infected individuals is proportional to the fraction of the susceptible individuals, multiplied by a so-called infection rate $\beta>0$. The model is given by

[TABLE]

and its dynamical behavior is given by the lemma below.

Lemma 1** (Dynamical behavior of the SI model).**

Consider the scalar SI model (1) with $\beta>0$. The solution from initial condition $x(0)=x_{0}\in[0,1]$ is

[TABLE]

All initial conditions $0<x_{0}<1$ result in the solution $x(t)$ being monotonically increasing and converging to the unique equilibrium $1$ as $t\to\infty$.

Solutions to equation (1) with different initial conditions are plotted in Figure 2. The SI model (1) results in an evolution akin to a logistic curve, and is also called the logistic equation for population growth.

3.2 Network SI model

The network SI model on a weighted digraph with the adjacency matrix $A\in\mathbb{R}^{n\times n}_{\geq 0}$ is given by

[TABLE]

or, in equivalent vector form,

[TABLE]

where $\beta>0$ is the infection rate. Alternatively, in terms of the fractions of susceptibile individuals $s(t)=\mathbbold{1}_{n}-x(t)$, the network SI model is

[TABLE]

The following results and their proof are novel.

Theorem 2** (Dynamical behavior of network SI model).**

Consider the network SI model (4) with $\beta>0$. For strongly connected graph with adjacency matrix $A$, the following statements hold:

1. (i.

if $x(0),s(0)\in[0,1]^{n}$, then $x(t),s(t)\in[0,1]^{n}$ for all $t>0$. Moreover, $x(t)$ is monotonically non-decreasing (here by monotonically non-decreasing we mean $x(t_{1})\leq x(t_{2})$ for all $t_{1}\leq t_{2}$). Finally, if $x(0)>\mathbbold{0}_{n}$, then $x(t)\gg\mathbbold{0}_{n}$ for all $t>0$; 2. (ii.

the model (4) has two equilibrium points: $\mathbbold{0}_{n}$ (no epidemic), and $\mathbbold{1}_{n}$ (full contagion);

1. (a)

the linearization of model (4) about the equilibrium point $\mathbbold{0}_{n}$ is $\dot{x}=\beta Ax$ and it is exponentially unstable; 2. (b)

let $D=\operatorname{diag}(A\mathbbold{1}_{n})$ be the degree matrix. The linearization of model (5) about the equilibrium $\mathbbold{0}_{n}$ is $\dot{s}=-\beta Ds$ and it is exponentially stable; 3. (iii.

each trajectory with initial condition $x(0)\neq\mathbbold{0}_{n}$ converges asymptotically to $\mathbbold{1}_{n}$, that is, the epidemic spreads monotonically to the entire network.

Proof.

(i) The fact that, if $x(0),s(0)\in[0,1]^{n}$, then $x(t),s(t)\in[0,1]^{n}$ for all $t>0$ means that $[0,1]^{n}$ is an invariant set for the differential equation (4). This is the consequence of Nagumo’s Theorem (see Theorem 4.7 in [5]), since for any $x$ belonging on the boundary of the set $[0,1]^{n}$, the vector \beta\Big{(}I_{n}-\operatorname{diag}\!\big{(}x\big{)}\Big{)}Ax is either tangent, or points inside the set $[0,1]^{n}$.

Observe that the invariance of the set $[0,1]^{n}$ implies that $\dot{x}(t)\geq\mathbbold{0}_{n}$ and so $x(t_{1})\leq x(t_{2})$ for all $t_{1}\leq t_{2}$.

We want to prove now that,if $x(0)>0_{n}$, then $x(t)\gg 0_{n}$ for all $t>0$. If by contradiction there is $i\in\{1,\ldots,n\}$ and $T>0$ such that $x_{i}(T)=0$, then the monotonicity of $x_{i}(t)=0$ would imply that $x_{i}(t)=0$ for all $t\in[0,T]$, which would yield $\dot{x}_{i}(t)=0$ for all $t\in[0,T]$. By (3) this would imply that $x_{j}(t)=0$ for all $t\in[0,T]$ for all $j$ such that $a_{ij}>0$. We could iterate this argument and using the irreducibility of $A$ we would get the contradiction that $x(t)=0$ for all $t\in[0,T]$ concluding in this way the proof of (i.

(ii) Regarding statement (ii, note that $\mathbbold{0}_{n}$ and $\mathbbold{1}_{n}$ are clearly equilibrium points. Let $\bar{x}\in[0,1]^{n}$ be an equilibrium and assume that $\bar{x}\not=\mathbbold{1}_{n}$. Then there is $i$ such that $\bar{x}_{i}\not=1$. Since \beta\big{(}1-\bar{x}_{i}\big{)}\sum_{j=1}^{n}a_{ij}\bar{x}_{j}=0, then $\sum_{j=1}^{n}a_{ij}\bar{x}_{j}=0$ which implies that $\bar{x}_{j}=0$ for all $j$ such that $a_{ij}>0$. By iterating this argument and using the irreducibility of $A$ we get that $\bar{x}=0$ concluding only $\mathbbold{0}_{n}$ and $\mathbbold{1}_{n}$ are equilibrium points. Statements (iia and (iib are obvious. Exponential stability of the linearization $\dot{s}=-\beta Ds$ is obvious, and the Perron-Frobenius Theorem implies the existence of the unstable positive eigenvalue $\rho(A)>0$ for the linearization $\dot{x}=\beta Ax$.

(iii) Consider the function $V(x)=\mathbbold{1}_{n}^{\top}(\mathbbold{1}_{n}-x)$; this is a smooth function defined over the compact and forward invariant set $[0,1]^{n}$ (see statement (i). Since \dot{V}=-\beta\mathbbold{1}_{n}^{\top}\big{(}I_{n}-\operatorname{diag}(x)\big{)}Ax, we know that $\dot{V}\leq 0$ for all $x$ and $\dot{V}(x)=0$ if and only if $x\in\{\mathbbold{0}_{n},\mathbbold{1}_{n}\}$. The LaSalle Invariance Principle implies that all trajectories with $x(0)$ converge asymptotically to either $\mathbbold{1}_{n}$ or $\mathbbold{0}_{n}$. Additionally, note that $0\leq V(x)\leq n$ for all $x\in[0,1]^{n}$, that $V(x)=0$ if and only if $x=\mathbbold{1}_{n}$ and that $V(x)=n$ if and only if $x=\mathbbold{0}_{n}$. Therefore, all trajectories with $x(0)\neq\mathbbold{0}_{n}$ converge asymptotically to $\mathbbold{1}_{n}$. ∎

For the adjacency matrix $A$, there exists a non-singular matrix $T$ such that $A=TJT^{-1}$, where $J$ is the Jordan normal form of $A$. Since $A$ is non-negative and irreducible, according to Perron-Frobenius theorem, the first Jordan block $J_{1}=(\lambda_{\textup{max}})_{1\times 1}$ and $\lambda_{\textup{max}}>Re(\lambda_{i})$ for any other eigenvalue $\lambda_{i}$ of $A$. Consider now the onset of an epidemic in a large population characterized by a small initial infection $x(0)=x_{0}$ much smaller than $\mathbbold{1}_{n}$. The system evolution is approximated by $\dot{x}=\beta Ax$. This “initial-times” linear evolution satisfies

[TABLE]

where $\mathbbold{e}_{1}$ is the first standard basis vector in $\mathbb{R}^{n}$ and $o(1)$ denotes a time-varying vector that vanishes as $t\to+\infty$. Let $u_{1}$ denote the first column of $T$ and let $v_{1}^{\top}$ denote the first row of $T^{-1}$. Since $AT=TJ$ and $T^{-1}A=JT^{-1}$, one can check that $u_{1}$ ($v_{1}$ resp.) is the right (left resp.) eigenvector of $A$ associated with the eigenvalue $\lambda_{\textup{max}}$. Since $T^{-1}T=I_{n}$, we have $v_{1}^{\top}u_{1}=1$. therefore,

[TABLE]

That is, the epidemic initially experiences exponential growth with rate $\beta\lambda_{\textup{max}}$ and with distribution among the nodes given by the eigenvector $u_{\max}$.

Now suppose that at some time $T$, for all $i$ we have that $x_{i}(T)=1-\epsilon_{i}$, where each $\epsilon_{i}$ is much smaller than $1$. Then, for time $t>T$, the approximated system for $s(t)$ is given by:

[TABLE]

From the discussion above, we conclude that the initial infection rate is proportional to the eigenvector centrality, and the final infection speed is proportional to the degree centrality.

4 Susceptible-Infected-Susceptible model

In this section we review the Susceptible-Infected-Susceptible (SIS) epidemic model. In addition to the existence of an infection process with rate $\beta>0$, this model assumes that the infected individuals recover to the susceptible state at so-called recovery rate $\gamma>0$.

4.1 Scalar SIS model

In the scalar SIS model, the population is divided into two fractions: the infected $x(t)$ and the susceptible $s(t)$, with $x(t)+s(t)=1$, obeying the following dynamics:

[TABLE]

The dynamical behavior of system (7) is given below.

Lemma 3** (Dynamical behavior of the SIS model).**

For the SIS model (7) with $\beta>0$ and $\gamma>0$:

1. (i.

the closed-form solution to equation (7) from initial condition $x(0)=x_{0}\in[0,1]$, for $\beta\neq\gamma$, is

[TABLE] 2. (ii.

if $\beta\leq\gamma$, all trajectories converge to the unique equilibrium $x=0$ (i.e., the epidemic disappears); 3. (iii.

if $\beta>\gamma$, then each trajectory from an initial condition $x(0)>0$ converges to the exponentially stable equilibrium $x^{*}=(\beta-\gamma)/\beta$, which is called the endemic state.

Case (iii corresponds to the case in which epidemic outbreaks take place and a steady-state epidemic contagion persists. The basic reproduction number in this deterministic scalar SIS model is given by $R_{0}=\beta/\gamma$. Simulations regarding to Lemma 3(ii and (iii are shown in Figure 3.

4.2 Network SIS Model

In this section we study the network SIS model which is closely related to the original “multi-group SIS model” proposed by Lajmanovich [14]; see also the intertwined SIS model in [20].

The network SIS model with infection rate $\beta$ and recovery rate $\gamma$ is given by:

[TABLE]

or, in equivalent vector form,

[TABLE]

In the rest of this section we study the dynamical properties of this model. We start by defining the monotonically-increasing functions

[TABLE]

for $y\in\mathbb{R}_{\geq 0}$ and $z\in[0,1[$. Note that $f_{+}(f_{-}(z))=z$ for all $z\in[0,1)$. For vector variables $y\in\mathbb{R}_{\geq 0}^{n}$ and $z\in[0,1)^{n}$, we write $F_{+}(y)=(f_{+}(y_{1}),\dots,f_{+}(y_{n}))$, and $F_{-}(z)=(f_{-}(z_{1}),\dots,f_{-}(z_{n}))$.

Behavior of System Below the Threshold

In this subsection, we characterize the behavior of the network SIS model in a regime we describe as “below the threshold.”

Theorem 4** (Dynamical behavior of the network SIS model: Below the threshold).**

Consider the network SIS model (9), with $\beta>0$ and $\gamma>0$, over a strongly connected digraph with adjacency matrix $A$. Let $\lambda_{\textup{max}}$ and $v_{\max}$ be the dominant eigenvalue of $A$ and the corresponding normalized left eigenvector respectively. If $\beta\lambda_{\textup{max}}/\gamma<1$, then

1. (i.

if $x(0),s(0)\in[0,1]^{n}$, then $x(t),s(t)\in[0,1]^{n}$ for all $t>0$. Moreover, if $x(0)>\mathbbold{0}_{n}$, then $x(t)\gg\mathbbold{0}_{n}$ for all $t>0$; 2. (ii.

there exists a unique equilibrium point $\mathbbold{0}_{n}$, the linearization of (9) about $\mathbbold{0}_{n}$ is $\dot{x}=(\beta A-\gamma I_{n})x$ and it is exponentially stable; 3. (iii.

from any $x(0)\neq\mathbbold{0}_{n}$, the weighted average $t\mapsto v_{\max}^{\top}x(t)$ is monotonically and exponentially decreasing, and all the trajectories converge to $\mathbbold{0}_{n}$.

Historically, it is meaningful to attribute this theorem to [14], even if the language adopted here is more modern.

Proof.

(i As in Theorem 2 the first part is the consequence of Nagumo’s Theorem. Then define $y(t):=e^{\gamma t}x(t)$. Notice that this variable satisfies the differential equation $\dot{y}(t)=\beta\operatorname{diag}(s(t))Ay(t)$. From the same arguments used in the proof of the point (i of Theorem 2 we argue that $y(t)\gg\mathbbold{0}_{n}$ for all $t>0$. From this it follows that also $x(t)\gg\mathbbold{0}_{n}$ for all $t>0$.

(ii Assume that $x^{*}$ is an equilibrium point. It is easy to se that $x^{*}\ll\mathbbold{1}_{n}$. Observe moreover that $x^{*}$ is an equilibrium point if and only if $\hat{A}x^{*}=F_{-}(x^{*})$ or, equivalently, if and only if F_{+}\big{(}\tfrac{\beta}{\gamma}Ax^{*}\big{)}=x^{*}. This means that $x^{*}$ is an equilibrium if and only if it is a fixed point of $\mathcal{F}$, where \mathcal{F}(x):=F_{+}\big{(}\tfrac{\beta}{\gamma}Ax\big{)}. Let $\hat{A}=\beta A/\gamma$. For $x\in[0,1]^{n}$, note $F_{+}(\hat{A}x)\leq\hat{A}x$ because $f_{+}(z)\leq z$. Moreover, $\mathbbold{0}_{n}\leq x\leq y$ implies that $\mathbbold{0}_{n}\leq\mathcal{F}(x)\leq\hat{A}y$. Therefore, if $\mathbbold{0}_{n}\leq x$, then $\mathcal{F}^{k}(x)\leq\hat{A}^{k}x$, for all $k$. Since $\hat{A}$ is Schur stable, then $\lim_{k\to\infty}\mathcal{F}^{k}(x)=0$. This shows that the only fixed point of $\mathcal{F}$ is zero.

Next, the linearization of equation (10) is verified by dropping the second-order terms. The linearized system is exponentially stable at $\mathbbold{0}_{n}$ for $\beta\lambda_{\textup{max}}-\gamma<0$ because $\lambda_{\textup{max}}$ is larger, in real part, than any other eigenvalue of $A$ by the Perron-Frobenius Theorem for irreducible matrices.

(iii Finally, regarding statement (iii, define $y(t)=v_{\max}^{\top}x(t)$ and note that \big{(}I_{n}-\operatorname{diag}(z)\big{)}v_{\max}\leq v_{\max} for any $z\in[0,1]^{n}$. Therefore,

[TABLE]

By the Grönwall-Bellman Comparison Lemma, $y(t)$ is monotonically decreasing and satisfies $y(t)\leq y(0)e^{(\beta\lambda_{\textup{max}}-\gamma)t}$ from all initial conditions $y(0)$. This concludes our proof of statement (iii. ∎

Behavior of System Above the Threshold

We present the dynamical behavior of the network SIS model above the threshold as follows.

Theorem 5** (Dynamical behavior of the network SIS model: Above the threshold).**

Consider the network SIS model (9), with $\beta>0$ and $\gamma>0$, over a strongly connected digraph with adjacency matrix $A$. Let $\lambda_{\textup{max}}$ be the dominant eigenvalue of $A$ and let $v_{\max}$ and $u_{\max}$ be the corresponding normalized left and right eigenvectors respectively. Let $d=A\mathbbold{1}_{n}$. If $\beta\lambda_{\textup{max}}/\gamma>1$, then

1. (i.

if $x(0),s(0)\in[0,1]^{n}$, then $x(t),s(t)\in[0,1]^{n}$ for all $t>0$. Moreover, if $x(0)>\mathbbold{0}_{n}$, then $x(t)\gg\mathbbold{0}_{n}$ for all $t>0$; 2. (ii.

$\mathbbold{0}_{n}$* is an equilibrium point, the linearization of system (10) at $\mathbbold{0}_{n}$ is unstable due to the unstable eigenvalue $\beta\lambda_{\textup{max}}-\gamma$ (i.e., there will be an epidemic outbreak);* 3. (iii.

besides the equilibrium $\mathbbold{0}_{n}$, there exists a unique equilibrium point $x^{*}$, called the endemic state, such that

1. (a)

$x^{*}\gg\mathbbold{0}_{n}$, 2. (b)

$x^{*}=\delta au_{\max}+O(\delta^{2})$* as $\delta\to 0^{+}$, where $\delta:=\beta\lambda_{\textup{max}}/\gamma-1$ and*

[TABLE] 3. (c)

$x^{*}=\mathbbold{1}_{n}-(\gamma/\beta)\operatorname{diag}(d)^{-1}\mathbbold{1}_{n}+O(\gamma^{2}/\beta^{2})$, at fixed $A$, as $\gamma/\beta\to 0^{+}$, 4. (d)

define a sequence $\{y(k)\}_{k\in\mathbb{N}}\subset\mathbb{R}^{n}$ by

[TABLE]

If $y(0)\geq 0$ is a scalar multiple of $u_{\max}$ and satisfies either $0<\max_{i}y_{i}(0)\leq 1-\gamma/(\beta\lambda_{\textup{max}})$ or $\min_{i}y_{i}(0)\geq 1-\gamma/(\beta\lambda_{\textup{max}})$, then

[TABLE]

Moreover, if $\max_{i}y_{i}(0)\leq 1-\gamma/(\beta\lambda_{\textup{max}})$, then $y(k)$ is monotonically non-decreasing; if $\min_{i}y_{i}(0)\geq 1-\gamma/(\beta\lambda_{\textup{max}})$, then $y(k)$ is monotonically non-increasing. 4. (iv.

the endemic state $x^{*}$ is locally exponentially stable and its domain of attraction is $[0,1]^{n}\setminus\mathbbold{0}_{n}$.

Note: statement (ii means that, near the onset of an epidemic outbreak, the exponential growth rate is $\beta\lambda_{\textup{max}}-\gamma$ and the outbreak tends to align with the dominant eigenvector $u_{\max}$; for more details see the discussion leading up to the approximate evolution (6). The basic reproduction number for this deterministic network SIS model is given by $R_{0}=\beta\lambda_{\textup{max}}/\gamma$.

Historically, the existence of a unique endemic state and its global attractivity properties are due to [14]. To the best of our knowledge, the Taylor expansions in parts (iiib and (iiic and the algorithm in part (iiid are novel. The proofs of statements (iii based on the properties of the map $F_{+}$ are novel.

Proof of selected statements in Theorem 5.

(i This point can be proved as done in point (i of Theorem 2.

(ii This follows from the same analysis of the linearized system as in the proof of Theorem 4(ii.

(iii We begin by establishing two properties of the map $x\mapsto F_{+}(\hat{A}x)$, for $\hat{A}=\beta A/\gamma$. First, we claim that, $y\gg z\geq\mathbbold{0}_{n}$ implies $F_{+}(\hat{A}y)\gg F_{+}(\hat{A}z)$. Indeed, note that $G$ being connected implies that the adjacency matrix $A$ has at least one strictly positive entry in each row. Hence, $y-z\gg\mathbbold{0}_{n}$ implies $\hat{A}(y-z)\gg\mathbbold{0}_{n}$ and, since $f_{+}$ is monotonically increasing, $\hat{A}y\gg\hat{A}z$ implies $F_{+}(\hat{A}y)\gg F_{+}(\hat{A}z)$.

Second, we observe that, for any $0<\alpha<1$ and $z>0$, we have $f_{+}(\alpha z)\geq z$ if and only if $z\leq 1-1/\alpha$. Suppose $y(0)$ is a scalar multiple of $u_{\max}$ and $0<\max_{i}y_{i}(0)\leq 1-\gamma/(\beta\lambda_{\textup{max}})$. We have

[TABLE]

Therefore, the sequence $\{y(k)\}_{k\in\mathbb{N}}$ defined by equation (11) satisfies $y(1)\geq y(0)$, which in turn leads to $y(2)=F_{+}(\hat{A}y(1))\geq F_{+}(\hat{A}y(0))=y(1)$, and by induction, $y(k+1)=F_{+}(\hat{A}y(k))\geq y(k)$ for any $k\in\mathbb{N}$. Such sequence $\{y(t)\}$ is monotonically non-decreasing and entry-wise upper bounded by $\mathbbold{1}_{n}$. Therefore, as $k$ diverges, $y(k)$ converges to some $x^{*}\gg\mathbbold{0}_{n}$ such that F_{+}\big{(}\hat{A}x^{*}\big{)}=x^{*}. This proves the existence of an equilibrium $x^{*}=\lim_{k\to\infty}y(k)\gg\mathbbold{0}_{n}$ as claimed in statements (iiia and (iiid.

Similarly, for any $0<\alpha<1$ and $z>0$, $f_{+}(\alpha z)\leq z$ if and only if $z\geq 1-1/\alpha$. Following the same line of argument in the previous paragraph, one can check that the $\{y(k)\}_{k\in\mathbb{N}}$ defined by equation (11) is monotonically non-increasing and converges to some $x^{*}$ , if $y(0)$ is a scalar multiple of $u_{\max}$ and satisfies $\min_{i}y_{i}(0)\geq 1-\gamma/(\beta\lambda_{\textup{max}})$.

Now we establish the uniqueness of the equilibrium $x^{*}\in[0,1]^{n}\setminus\{\mathbbold{0}_{n}\}$. First, we claim that an equilibrium point with an entry equal to [math] must be $\mathbbold{0}_{n}$. Indeed, assume $y^{*}$ is an equilibrium point and assume $y_{i}^{*}=0$ for some $i\in\{1,\dots,n\}$. The equality $y^{*}_{i}=f_{+}(\sum_{j=1}^{n}a_{ij}y^{*}_{j})$ implies that also any node $j$ with $a_{ij}>0$ must satisfy $y^{*}_{j}=0$. Because $G$ is connected, all entries of $y^{*}$ must be zero. Second, by contradiction, we assume there exists another equilibrium point $y^{*}\gg\mathbbold{0}_{n}$ distinct from $x^{*}$. Let $\alpha:=\min_{j}\{y^{*}_{j}/x^{*}_{j}\}$ and let $i$ such that $\alpha=y^{*}_{i}/x^{*}_{i}$. Then $y^{*}\geq\alpha x^{*}\gg\mathbbold{0}_{n}$ and $y_{i}^{*}=\alpha x_{i}^{*}$. Notice that we can assume with no loss of generality that $\alpha<1$ otherwise we exchange $x^{*}$ and $y^{*}$. Observe now that

[TABLE]

Therefore, \big{(}F_{+}(\hat{A}y^{*})-y^{*}\big{)}_{i}>0, which contradicts the fact that $y^{*}$ is an equilibrium.

Now we prove (iiib. Observe first that, since taking

[TABLE]

then $y(k)$ is monotonically non-decreasing and converges to $x^{*}$, and since taking instead

[TABLE]

then $y(k)$ is monotonically non-increasing and converges to $x^{*}$, we can argue that

[TABLE]

This implies that $x^{*}$ is infinitesimal as a function of $\delta$. Consider the expansion $x^{*}(\delta)=x_{1}\delta+x_{2}\delta^{2}+O(\delta^{3})$. Since the equilibrium $x^{*}$ satisfies the equation

[TABLE]

by substituting the expansion and equating to zero the coefficient of the term $\delta$ we obtain the equation

[TABLE]

which proves that $x_{1}$ is a multiple of $u_{\max}$, namely $x_{1}=au_{\max}$ for some constant $a$. By equating to zero the coefficient of the term $\delta^{2}$ we obtain instead the equation

[TABLE]

Using the fact that $x_{1}=au_{\max}$ we argue that

[TABLE]

By multiplying on the left by $v_{\max}^{T}$ we obtain

[TABLE]

which proves that

[TABLE]

Point  (iiic can be proved in a similar way. Indeed, define $\epsilon:=\gamma/\beta$. Since

[TABLE]

we can argue that the expansion $x^{*}(\epsilon)=x_{0}+x_{1}\epsilon+O(\epsilon^{2})$ as $\epsilon$ tends to zero is such that $x_{0}\gg\mathbbold{0}_{n}$. Since the equilibrium $x^{*}$ satisfies the equation

[TABLE]

by substituting the expansion and equating to zero the coefficient of the term $\epsilon^{0}$ we obtain the equation

[TABLE]

which proves that $x_{0}=vectorones[n]$. By equating to zero the coefficient of the term $\epsilon^{1}$ we obtain instead the equation

[TABLE]

Using the fact that $x_{0}=\mathbbold{1}_{n}$ we argue that

[TABLE]

which yieds the thesis.

(iv For this point we refer to [14, 9] or [13, Theorems 1 and 2] in the interest of brevity. ∎

5 Network Susceptible-Infected-Recovered Model

In this section we review the Susceptible-Infected-Susceptible (SIR) epidemic model.

5.1 Scalar SIR model

In this model individuals who recover from infection are assumed not susceptible to the epidemic any more. In this case, the population is divided into three distinct groups: $s(t)$, $x(t)$, and $r(t)$, denoting the fraction of susceptible, infected, and recovered individuals, respectively, with $s(t)+x(t)+r(t)=1$. We write the (Susceptible–Infected–Recovered) SIR model as:

[TABLE]

and present its dynamical behavior in the lemma below.

Lemma 6** (Dynamical behavior of the SIR model).**

Consider the SIR model (12). From each initial condition $s(0)+x(0)+r(0)=1$ with $s(0)>0$, $x(0)>0$ and $r(0)\geq 0$, the resulting trajectory $t\mapsto(s(t),x(t),r(t))$ has the following properties:

1. (i.

$s(t)>0$, $x(t)>0$, $r(t)\geq 0$, and $s(t)+x(t)+r(t)=1$ for all $t\geq 0$; 2. (ii.

$t\mapsto s(t)$* is monotonically decreasing and $t\mapsto r(t)$ is monotonically increasing;* 3. (iii.

$\lim_{t\to\infty}(s(t),x(t),r(t))=(s_{\infty},0,r_{\infty})$, where $r_{\infty}$ is the unique solution to the equality

[TABLE] 4. (iv.

if $\beta s(0)/\gamma<1$, then $t\mapsto x(t)$ monotonically and exponentially decreases to zero as $t\to\infty$; 5. (v.

if $\beta s(0)/\gamma>1$, then $t\mapsto x(t)$ first monotonically increases to a maximum value and then monotonically decreases to [math] as $t\to\infty$; the maximum fraction of infected individuals is given by:

[TABLE]

As mentioned before, we describe the behavior in statement (v as an epidemic outbreak, an exponential growth of $t\mapsto x(t)$ for small times.) The effective reproduction number in the deterministic scalar SIR model is $R=\beta s(t)/\gamma$. Note that the basic reproduction number $R_{0}=\beta/\gamma$ does not have predict power in this model.

5.2 Network SIR model

The network SIR model on a graph with adjacency matrix $A$ is given by

[TABLE]

where $\beta>0$ is the infection rate and $\gamma>0$ is the recovery rate. Note that the third equation is redundant because of the constraint $s_{i}(t)+x_{i}(t)+r_{i}(t)=1$. Therefore, we regard the dynamical system in vector form as:

[TABLE]

We state our main novel results of this section below.

Theorem 7** (Dynamical behavior of the network SIR model).**

Consider the network SIR model (14), with $\beta>0$ and $\gamma>0$, over a strongly connected digraph with adjacency matrix $A$. For $t\geq 0$, let $\lambda_{\textup{max}}(t)$ and $v_{\max}(t)$ be the dominant eigenvalue of the non-negative matrix $\operatorname{diag}(s(t))A$ and the corresponding normalized left eigenvector, respectively. The following statements hold:

1. (i.

if $x(0)>\mathbbold{0}_{n}$, and $s(0)\gg\mathbbold{0}_{n}$, then

1. (a)

$t\mapsto s(t)$* and $t\mapsto x(t)$ are strictly positive for all $t>0$,* 2. (b)

$t\mapsto s(t)$* is monotonically decreasing, and* 3. (c)

$t\mapsto\lambda_{\textup{max}}(t)$* is monotonically decreasing;* 2. (ii.

the set of equilibrium points is the set of pairs $(s^{*},\mathbbold{0}_{n})$, for any $s^{*}\in[0,1]^{n}$, and the linearization of model (14) about $(s^{*},\mathbbold{0}_{n})$ is

[TABLE] 3. (iii.

(behavior below the threshold) let the time $\tau\geq 0$ satisfy $\beta\lambda_{\textup{max}}(\tau)<\gamma$. Then the weighted average $t\mapsto v_{\max}(\tau)^{\top}x(t)$ , for $t\geq\tau$, is monotonically and exponentially decreasing to zero; 4. (iv.

(behavior above the threshold) if $\beta\lambda_{\textup{max}}(0)>\gamma$ and $x(0)>\mathbbold{0}_{n}$, then,

1. (a)

(epidemic outbreak) for small time, the weighted average $t\mapsto v_{\max}(0)^{\top}x(t)$ grows exponentially fast with rate $\beta\lambda_{\textup{max}}(0)-\gamma$, and 2. (b)

there exists $\tau>0$ such that $\beta\lambda_{\textup{max}}(\tau)<\gamma$; 5. (v.

each trajectory converges asymptotically to an equilibrium point, that is, $\lim_{t\mapsto\infty}x(t)=\mathbbold{0}_{n}$ so that the epidemic asymptotically disappears.

The effective reproduction number in the deterministic network SIR model is $R(t)=\beta\lambda_{\textup{max}}(t)/\gamma$. When $R(0)>1$, we have an epidemic outbreak, i.e., an exponential growth of infected individual for short time. In any case, the theorem guarantees that, after at most finite time, $R(t)<1$ and the infected population decreases exponentially fast to zero.

Proof.

Regarding statement (ia, $s(t)\gg\mathbbold{0}_{n}$ is due to the fact that $Ax$ is bounded and $s(t)$ is continuously differentiable to $t$. The statement that $x(t)\gg\mathbbold{0}_{n}$ for all $t>0$ is proved in the same way as Theorem 4 (i. Statement (ib is the immediate consequence of $\dot{s}_{i}(t)$ being strictly negative. From statement (ia we know that each $s_{i}(t)$ is positive, and from $A$ being irreducible and $x(0)\neq\mathbbold{0}_{n}$ we know that $\sum_{j=1}^{n}a_{ij}x_{j}$ is positive. Therefore, $\dot{s}_{i}(t)=-\beta s_{i}(t)\sum_{j=1}^{n}a_{ij}x_{j}(t)<0$ for all $i\in V$ and $t\geq 0$.

For statement (ic, we start by recalling the following property from [16, Example 7.10.2]: for $B$ and $C$ nonnegative square matrices, if $B\leq C$, then $\rho(B)\leq\rho(C)$. Now, pick two time instances $t_{1}$ and $t_{2}$ with $0<t_{1}<t_{2}$. Let $\alpha=\max_{i}s_{i}(t_{2})/s_{i}(t_{1})$ and note $0<\alpha<1$ because $s(t)$ is strictly positive and monotonically decreasing. Now note that,

[TABLE]

so that, using the property above, we know

[TABLE]

This concludes the proof of statement (ic.

Regarding statement (ii, note that a point $(s^{*},x^{*})$ is an equilibrium if and only if:

[TABLE]

Therefore, each point of the form $(s^{*},\mathbbold{0}_{n})$ is an equilibrium. On the other hand, summing the last two equalities we obtain $\mathbbold{0}_{n}=\gamma x^{*}$ and thus $x^{*}$ must be $\mathbbold{0}_{n}$. As a straightforward result, the linearization of model (14) about any equilibrium point $(s^{*},\mathbbold{0}_{n},\mathbbold{1}_{n}-s^{*})$ is given by equation (15).

Regarding statement (iii, multiplying $v_{\max}(\tau)^{\top}$ from the left on both sides of equation (14b) we obtain:

[TABLE]

Therefore, we obtain

[TABLE]

The right-hand side exponentially decays to zero when $\beta\lambda_{\textup{max}}(\tau)<\gamma$. Therefore, $v_{\max}(\tau)^{\top}x(t)$ also decreases monotonically and exponentially to zero for all $t>\tau$.

Regarding statement (iva, note that based on the argument in (ia, we only need to consider the case when $x(0)\gg\mathbbold{0}_{n}$. Left-multiplying $v_{\max}(0)^{\top}$ on both sides of equation (14b), we obtain:

[TABLE]

Since $\beta\lambda_{\textup{max}}(0)-\gamma>0$, the initial time derivative of $v_{\max}(0)^{\top}x(t)$ is positive. Since $t\mapsto v_{\max}(0)^{\top}x(t)$ is a continuously differentiable function, there exists $\tau^{\prime}>0$ such that \frac{d}{dt}\big{(}v_{\max}(0)^{\top}x(t)\big{)}>0 for any $t\in[0,\tau^{\prime}]$.

Regarding statement (ivb, since $\dot{s}(t)\leq\mathbbold{0}_{n}$ and is lower bounded by $\mathbbold{0}_{n}$, we conclude that the limit $\lim\limits_{t\to+\infty}s(t)$ exists. Moreover, since $s(t)$ is monotonically non-increasing, we have $\lim\limits_{t\to+\infty}\dot{s}(t)=0$, which implies either $\lim\limits_{t\to+\infty}s(t)=\mathbbold{0}_{n}$ or $\lim\limits_{t\to+\infty}x(t)=\mathbbold{0}_{n}$. If $s(t)$ converges to $\mathbbold{0}_{n}$, then $\dot{x}(t)$ converges to $-\gamma x(t)$. Therefore, there exists $T>0$ such that $\beta\lambda_{\textup{max}}(T)<\gamma$, which leads to $x(t)\to\mathbbold{0}_{n}$ as $t\to+\infty$; If $s(t)$ converges to some $s^{*}>\mathbbold{0}_{n}$, then $x(t)$ still converges to $\mathbbold{0}_{n}$. Therefore, for any \big{(}s(0),x(0)\big{)}, the trajectory \big{(}s(t),x(t)\big{)} converges to some equilibria with the form $(s^{*},\mathbbold{0}_{n})$, where $s^{*}\geq\mathbbold{0}_{n}$. Let

[TABLE]

We know that $\delta_{s}(t)\geq 0$ and $\delta_{x}(t)\geq 0$ for all $t\geq 0$. Moreover, $\delta_{s}(t)$ is monotonically non-increasing and converges to $\mathbbold{0}_{n}$, and there exists $\tilde{T}>0$ such that, for any $t\geq T$, $\delta_{x}(t)$ is monotonically non-increasing and converges to $\mathbbold{0}_{n}$.

Let $\lambda^{*}$ and $v^{*}$ denote the dominant eigenvalue and the corresponding normalized left eigenvector of matrix $\operatorname{diag}(s^{*})A$, respectively, that is, ${v^{*}}^{\top}\operatorname{diag}(s^{*})A=\lambda^{*}{v^{*}}^{\top}$. First let us suppose $\beta\lambda^{*}-\gamma>0$, then the linearized system of (12) around $(s^{*},\mathbbold{0}_{n})$ is written as

[TABLE]

Since $\beta\lambda^{*}-\gamma>0$, the linearized system is exponentially unstable, which contradicts the fact that \big{(}\delta_{s}(t),\delta_{x}(t)\big{)}\to(\mathbbold{0}_{n},\mathbbold{0}_{n}) as $t\to+\infty$. Alternatively, suppose $\beta\lambda^{*}-\gamma=0$. By left multiplying ${v^{*}}^{\top}$ on both sides of the equation for $\dot{x}(t)$ in (12), we obtain

[TABLE]

which contradicts $\delta_{x}(t)\to\mathbbold{0}_{n}$ as $t\to+\infty$. Therefore, we conclude that $\beta\lambda^{*}-\gamma<0$. Since $\lambda_{\textup{max}}(t)$ is continuous on $t$, we conclude that there exists $\tau<+\infty$ such that $\beta\lambda_{\textup{max}}(t)-\gamma<0$.

∎

In what follows, we present an iterative algorithm that computes the limit state \lim_{t\to\infty}\big{(}s(t),0,r(t)\big{)} of the network SIR model (14) as a function of an arbitrary initial condition \big{(}s(0),x(0),r(0)\big{)}. To our best knowledge, this problem and its solution are novel.

Note that, for the scalar SIR model (12), if we define

[TABLE]

Simple calculations result in dV\big{(}s(t),x(t)\big{)}/dt=0, which implies that the trajectories are on the level sets of $V$ and in the set $\{(s,x)\in\mathbb{R}^{2}\;|\;s\geq 0,x\geq 0,s+x\leq 1\}$. Here, we apply a similar approach to the network SIR system (14). Let

[TABLE]

One can check that, along any trajectory of dynamics (14), $dV_{i}/dt=0$ for any $i\in\{1,\dots,n\}$. Therefore, the trajectories $(s(t),r(t))$ lie on the level curves of the functions $V_{i}(s,r)$ for $i\in\{1,\dots,n\}$.

Let $s(\infty):=\lim_{t\to+\infty}s(t)$, $x(\infty):=\lim_{t\to+\infty}x(t)$, and $r(\infty):=\lim_{t\to+\infty}r(t)$. Notice that $x(\infty)=\mathbbold{0}_{n}$ and so $r(\infty)=\mathbbold{1}_{n}-s(\infty)$. Since $dV_{i}/dt=0$ for any $i\in\{1,\dots,n\}$, we have

[TABLE]

Given any initial condition \big{(}s(0),r(0)\big{)}, the right-hand side of equation (16) defines a map

[TABLE]

and $s(\infty)$ is a fixed point of $H$, that is, s(\infty)=H\big{(}s(\infty)\big{)}.

Theorem 8** (Existence, uniqueness, and algorithm for the asymptotic point).**

Consider the network SIR model (14), with positive rates $\beta$ and $\gamma$ and with initial condition \big{(}s(0),x(0),r(0)\big{)} satisfying $s(0)\geq\mathbbold{0}_{n}$, $x(0)>\mathbbold{0}_{n}$, $r(0)\geq\mathbbold{0}_{n}$ and $s(0)+x(0)+r(0)=\mathbbold{1}_{n}$. Let \big{(}s(\infty),\mathbbold{0}_{n},r(\infty)\big{)} be the asymptotic state of system (14). The map $H:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ has the following properties:

1. (i.

there exists a unique fixed point $s^{*}$ of the map $H$ in the set $\{s\in\mathbb{R}^{n}\;|\;\mathbbold{0}_{n}\leq s\leq\mathbbold{1}_{n}-r(0)\}$. Moreover, $s^{*}=s(\infty)$ and $r(\infty)=\mathbbold{1}_{n}-s^{*}$; and 2. (ii.

any sequence $\{y(k)\}_{k\in\mathbb{N}}$ defined by $y(k+1)=H(y(k))$ and initial condition $\mathbbold{0}_{n}\leq y(0)\leq\mathbbold{1}_{n}-r(0)$ converges to the unique fixed point $s^{*}$.

Proof.

Since $A$ is a non-negative matrix, and $s(0)\leq\mathbbold{1}-r(0)$, one can easily observe that, if $\mathbbold{0}_{n}\leq p\leq q\leq\mathbbold{1}_{n}-r(0)$, then $\mathbbold{0}_{n}\leq H(\mathbbold{0}_{n})\leq H(p)\leq H(q)\leq H(\mathbbold{1}_{n}-r(0))\leq\mathbbold{1}_{n}-r(0)$. According to the Brower Fixed Point Theorem, the map $H$ has at least one fixed point.

Define the sequence $\{p(k)\}_{k\in\mathbb{N}}$ by $p(k+1)=H(p(k))$ and $p(0)=\mathbbold{0}_{n}$. Since

[TABLE]

we have $\mathbbold{1}_{n}-r(0)\geq p(2)=H(p(1))\geq H(p(0))=p(1)$ and, by induction, $\mathbbold{1}_{n}-r(0)\geq p(k+1)\geq p(k)$ for any $k\in\mathbb{N}$. Since $p(k)$ is non-decreasing and upper bounded by $\mathbbold{1}_{n}-r(0)$, we conclude that the limit $p^{*}=\lim_{k\to\infty}p(k)$ exists, and $p^{*}$ is a fixed point of the map $H$.

Similarly, define a sequence $\{q(k)\}_{k\in\mathbb{N}}$ by $q(k+1)=H(q(k))$ and $q(0)=\mathbbold{1}_{n}-r(0)$. One can check that $q(k)$ is non-increasing and that $q^{*}=\lim_{k\to\infty}q(k)$ is a fixed point of map $H$. Moreover, since $p(0)\leq q(0)$, we have $p(k)\leq q(k)$ for any $k\in\mathbb{N}$ and thereby $p^{*}\leq q^{*}$.

If $p^{*}=q^{*}$, then, for any $\mathbbold{0}_{n}\leq y(0)\leq\mathbbold{1}_{n}-r(0)$, the sequence $\{y(k)\}_{k\in\mathbb{N}}$ defined by $y(k+1)=H(y(k))$ satisfies $p(k)\leq y(k)\leq q(k)$ for any $k\in\mathbb{N}$. Therefore, $y^{*}=\lim_{k\to\infty}y(k)$ exists and $y^{*}=p^{*}=q^{*}$, which implies that the fixed point of map $H$ is unique. According to equation (16), $s(\infty)$ is the unique fixed point. This concludes the proof for statement (i) and (ii).

Now we eliminate the case $p^{*}<q^{*}$ by contradiction. First of all we prove that $q^{*}\ll\mathbbold{1}_{n}-r(0)$. Let $N_{i}=\{j\,|\,a_{ij}>0\}$ and \mathcal{I}(k)=\big{\{}i\,\big{|}\,q_{i}(\tau)<1-r_{i}(0)\text{ for any }\tau\geq k\big{\}}. We have $\mathcal{I}(0)=\phi$, and \mathcal{I}(1)=\big{\{}i\,\big{|}\,s_{i}(0)<1-r_{i}(0)\big{\}}, since $q(k+1)=H(q(k))\leq s(0)$ for any $k\in\mathbb{N}$. Moreover, since, for any $i$ such that $N_{i}\cap\mathcal{I}(k)\neq\phi$,

[TABLE]

we have \mathcal{I}(k+1)=\big{\{}i\,\big{|}\,N_{i}\cap\mathcal{I}(k)\neq\phi\big{\}}\cup\mathcal{I}(k). Because the graph associated with $A$ is strongly connected, we can argue that $\mathcal{I}(k)$ contains all the indices when $k$ is large enough. Therefore, $q^{*}\ll\mathbbold{1}_{n}-r(0)$.

Now suppose $p^{*}<q^{*}$. Let

[TABLE]

We have $\alpha>1$, $\mathbbold{0}_{n}\leq w<\mathbbold{1}_{n}-r(0)$, and $w_{i}=1-r_{i}(0)$ for any $i$ such that \alpha_{i}=\big{(}1-r_{i}(0)-p_{i}^{*}\big{)}/(q_{j}^{*}-p_{j}^{*}). Let $\mu=1/\alpha$. Thereby $q^{*}=\mu w+(1-\mu)p^{*}$, where $0<\mu<1$. This means that $q^{*}$ is a convex combination of $p^{*}$ and $w$. Since $H(s)_{i}$ is a strictly convex function of $s$, we obtain that

[TABLE]

In the last inequality, we used the fact that $H(w)_{i}\leq 1-r_{i}(0)$ for any $0\leq w\leq\mathbbold{1}_{n}-r(0)$. The previous inequality yields a contradiction. ∎

In the rest of this section, we present some numerical results for the network SIR model on the undirected unweighted graph illustrated in Figure 5. The adjacency matrix $A$ is binary. Unless otherwise stated, the system parameters are $\beta=0.5$ and $\gamma=0.4$. As initial condition, we select one node fully infected (the dark-gray node in Figure 5, say, with index $1$), 19 fully healthy individuals, and zero recovered fraction — corresponding to $x(0)=\mathbbold{e}_{1}$, $r(0)=\mathbbold{0}_{n}$, and $s(0)=\mathbbold{1}_{n}-x(0)$. These parameters lead to an initial effective reproduction number $R(0)=3.57$.

Figure 6 illustrates the time evolution of $(\beta/\gamma)\lambda_{\textup{max}}(t)$ with varying network parameters. Note that each evolution starts above the threshold, reaches the threshold value $1$ in finite time, and converges to a final value below $1$.

Figure 7 illustrates the behavior of the average susceptible, average infected and average recovered quantities in populations starting from a small initial infection fraction and with an effective reproduction number above $1$ at time [math]. Note that the evolution of the infected fraction of the population displays a unimodal dependence on time, like in the scalar model.

6 Conclusion

This paper provides a comprehensive and consistent treatment of deterministic nonlinear continuous-time SI, SIS, and SIR propagation models over contact networks. We investigated the asymptotic behaviors (vanishing infection, steady-state epidemic, and full contagion). We studied the transient propagation of an epidemic starting from small initial fractions of infected nodes. We presented conditions under which a possible epidemic outbreak occurs or the infection monotonically vanishes for arbitrary fixing topology graphs. We introduced a network SI model and analyzed its behavior. Network SIS model sections includes improved properties over previously proposed works. New transient behavior, threshold condition, and system properties for the network SIR model were proposed. In addition, for the network SIR model, we provide a novel iterative algorithm to compute the asymptotic state of the system. In all cases, we show the results for network models are appropriate generalizations of those for the respective scalar models.