Analysis and Control of a Continuous-Time Bi-Virus Model

TL;DR

This paper analyzes the stability and equilibrium behaviors of a continuous-time bi-virus model on networks, revealing conditions for virus eradication, dominance, or coexistence, and explores decentralized control limitations.

Contribution

It provides a comprehensive analysis of bi-virus spread dynamics on directed networks, including stability conditions and control limitations, which advances understanding of epidemic control in complex networks.

Findings

Healthy state is globally stable under certain conditions.

Single-virus dominance occurs with almost global stability.

Coexistence of viruses is possible with multiple equilibria.

Abstract

This paper studies a distributed continuous-time bi-virus model in which two competing viruses spread over a network consisting of multiple groups of individuals. Limiting behaviors of the network are characterized by analyzing the equilibria of the system and their stability. Specifically, when the two viruses spread over possibly different directed infection graphs, the system may have (1) a unique equilibrium, the healthy state, which is globally stable, implying that both viruses will eventually be eradicated, (2) two equilibria including the healthy state and a dominant virus state, which is almost globally stable, implying that one virus will pervade the entire network causing a single-virus epidemic while the other virus will be eradicated, or (3) at least three equilibria including the healthy state and two dominant virus states, depending on certain conditions on the healing and infection rates. When the two viruses spread over the same directed infection graph, the system may have zero or infinitely many coexisting epidemic equilibria, which represents the pervasion of the two viruses. Sensitivity properties of some nontrivial equilibria are investigated in the context of a decentralized control technique, and an impossibility result is given for a certain type of distributed feedback controller.

Equations

$$\rho (M)$$

$$s(M)$$

$${\rm sgn}(x)=\left\{\begin{array}[]{lll}-1&\mbox{ if \ $x<0$},\\ 0&\mbox{ if \ $x=0$},\\ 1&\mbox{ if \ $x>0$}.\end{array}\right.$$

$${\dot{S}}_i(t)$$

$$-\sum _{j=1}^n{\alpha }_{ij}^1\frac{S_i(t)}{N_i}{I}_j^1(t)-\sum _{j=1}^n{\alpha }_{ij}^2\frac{S_i(t)}{N_i}{I}_j^2(t)$$

$$=({\mu }_i+{\gamma }_i^1)I_i^1(t)+({\mu }_i+{\gamma }_i^2)I_i^2(t)$$

$$-\sum _{j=1}^n{\alpha }_{ij}^1\frac{S_i(t)}{N_i}{I}_j^1(t)-\sum _{j=1}^n{\alpha }_{ij}^2\frac{S_i(t)}{N_i}{I}_j^2(t),$$

$${\dot{I}}_i^1(t)$$

$$=(-{\gamma }_i^1-{\mu }_i)I_i^1(t)+\sum _{j=1}^n{\alpha }_{ij}^1\frac{S_i(t)}{N_i}{I}_j^1(t),$$

$${\dot{I}}_i^2(t)$$

$$=(-{\gamma }_i^2-{\mu }_i)I_i^2(t)+\sum _{j=1}^n{\alpha }_{ij}^2\frac{S_i(t)}{N_i}{I}_j^2(t),$$

$$x_i^1(t)=\frac{I_i^1(t)}{N_i},x_i^2(t)=\frac{I_i^2(t)}{N_i},$$

$${\beta }_{ij}^1={\alpha }_{ij}^1\frac{N_j}{N_i},{\beta }_{ij}^2={\alpha }_{ij}^2\frac{N_j}{N_i},{\delta }_i^1={\gamma }_i^1+{\mu }_i,{\delta }_i^2={\gamma }_i^2+{\mu }_i.$$

$$\begin{array}{rl}{\dot{x}}_i^1(t)& =-{\delta }_i^1{x}_i^1(t)+(1-x_i^1(t)-x_i^2(t))\sum _{j=1}^n{\beta }_{ij}^1{x}_j^1(t),\\ {\dot{x}}_i^2(t)& =-{\delta }_i^2{x}_i^2(t)+(1-x_i^2(t)-x_i^1(t))\sum _{j=1}^n{\beta }_{ij}^2{x}_j^2(t).\end{array}$$

$$\frac{{\dot{I}}_i^2(t)}{N_i}=(-{\gamma }_i^2-{\mu }_i)\frac{I_i^2(t)}{N_i}+\sum _{j=1}^n{\alpha }_{ij}^2\frac{S_i(t)}{N_i}\frac{N_j}{N_i}\frac{I_j^2(t)}{N_j}.$$

$$\frac{{\dot{I}}_i^2(t)}{N_i}=(-{\gamma }_i^2-{\mu }_i)\frac{I_i^2(t)}{N_i}+\frac{N_i-I_i^1(t)-I_i^2(t)}{N_i}\sum _{j=1}^n{\alpha }_{ij}^2\frac{N_j}{N_i}\frac{I_j^2(t)}{N_j}.$$

$$x_i^1(t)=\frac{I_i^1(t)}{N_i},x_i^2(t)=\frac{I_i^2(t)}{N_i},{\beta }_{ij}^2={\alpha }_{ij}^2\frac{N_j}{N_i},{\delta }_i^2={\gamma }_i^2+{\mu }_i,$$

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

$$\begin{array}{rl}{\dot{x}}^1(t)& =(-D^1+B^1-X^1(t)B^1-X^2(t)B^1)x^1(t),\\ {\dot{x}}^2(t)& =(-D^2+B^2-X^2(t)B^2-X^1(t)B^2)x^2(t),\end{array}$$

$$q_{kl}=\{\begin{array}{ll}{\delta }_i^1,& \text{ if }{s}_i=1,k=l+3^{i-1}\\ {\delta }_i^2,& \text{ if }{s}_i=2,k=l+2(3^{i-1})\\ \sum _{j=1}^n{\beta }_{ij}^1{1}_{s_j=1},& \text{ if }{s}_i=0,k=l-3^{i-1}\\ \sum _{j=1}^n{\beta }_{ij}^2{1}_{s_j=2},& \text{ if }{s}_i=0,k=l-2(3^{i-1})\\ -\sum _{j\ne l}{q}_{jl},& \text{ if }k=l\\ 0,& \text{ otherwise,}\end{array}$$

$$y_k(t)=Pr[Y_k(t)=k],$$

$$\frac{d{y}^{\prime }(t)}{dt}=y^{\prime }(t)Q.$$

$$\begin{array}{r}(v^1{)}^{\prime }(t)=y^{\prime }(t)M^1,\\ (v^2{)}^{\prime }(t)=y^{\prime }(t)M^2,\end{array}$$

$$M^1$$

$$M^2$$

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

$$\dot{x}(t)=\left(-D+B-X(t)B\right)x(t).$$

$${\dot{x}}_i(t)=-{\delta }_i{x}_i(t)+{\beta }_{ii}(1-x_i(t))x_i(t),$$

$$\mathcal{D}=\{(x^1,x^2)|x^1\ge 0,x^2\ge 0,x^1+x^2\le 1\}$$

$$\dot{x}(t)=f(x(t)),$$

$${\dot{x}}_i^1(t)$$

$${\dot{x}}_i^2(t)$$

$$\dot{V}(x^1(t),x^2(t))$$

$${\dot{x}}^1(t)=(-D^1+B^1-X^1(t)B^1)x^1(t)-X^2(t)B^1{x}^1(t).$$

Full text

I Introduction

The spread of epidemic processes over large populations is an important research topic, and is in fact a widely studied one in epidemiology [2]. To model such a process, various epidemic models have been proposed such as the susceptible-infected-recovered (SIR), susceptible-exposed-infected-recovered (SEIR), and susceptible-infected-susceptible (SIS) models [3, 4, 5, 6]. Bernoulli developed one of the first known models inspired by the smallpox virus [7]. The first SIS model was introduced in [8]. In this paper, we focus on the study of distributed SIS epidemic models, where there are two ways to consider such a system: 1) the model consists of $n>1$ interacting individuals and the evolution of the probability of each individual being infected is studied, or 2) the model consists of $n>1$ groups of individuals and the evolution of the percentage of infected members of each group is studied. The first type of SIS model has been studied in both discrete-time [9, 10, 11, 12, 13] and continuous-time [14, 15, 16, 17, 18, 19, 20, 21, 22]. The first multi-agent, probability-based, continuous-time model was proposed by Van Mieghem et al. [14] in which the underlying neighbor graph is assumed to be undirected. The same model on a directed neighbor graph has been recently studied by Khanafer et al. [20] for both strongly and weakly connected neighbor graphs. The second, or group-based type of models, has been studied in [23].

The idea of competing SIS virus models is pursued in [24, 25, 26, 27, 28, 29, 30]. This work is motivated by the competition of different viral strains [24], where there are two competing viruses in a human contact network. These models have a wide range of other applications, including social networks, where the goal is to understand how competing opinions spread on different social networks [26], competing products in a market, and agents’ opinions about politicians from opposing parties [31]. Competing SIS models were first introduced in [24], which is an extension of [8], where the model considers the dynamics of three groups: 1) susceptible, 2) infected with virus one, and 3) infected with virus two. These dynamics are modeled by three differential equations where full connectivity of the agents is assumed (i.e., the infection graph is a complete graph), and it is also assumed that the two viruses are both homogeneous.111 We say that a virus is homogeneous if all agents have the same infection rate and healing rate. Otherwise, the virus is called heterogeneous.

In [25], two competing homogeneous viruses spreading over the same nontrivial (not necessarily fully connected), undirected, connected network is studied. The set of equilibrium points is determined and sufficient conditions for local stability are given for all equilibria except the coexisting equilibrium. In [26], the equilibria of two competing homogeneous virus models over the same as well as different undirected graph structures, are studied. Existence of the coexisting epidemic states, where both viruses are at nontrivial (nonzero) equilibria, is shown, but no stability analysis is provided. Note that all this previous work is conducted for homogeneous viruses over undirected graph structures with limited/local stability analysis. The following are the two exceptions: in [28], a sufficient condition for the global asymptotic survival of a single virus is given for a model of two competing viruses, both homogeneous in the healing rate and propagating over undirected, regular graphs. In [30, 32], a necessary and sufficient condition for local exponential stability of the origin is provided for two competing heterogeneous viruses over strongly connected graphs. In addition, a geometric program is formulated, working toward optimal stabilization and rate control of the virus. However, stability of the epidemic equilibria (i.e., nonzero equilibria) is not explored. Note that none of the existing work considers heterogeneous viruses over directed graph structures and performs global stability analysis, exploring all of the system’s equilibria.

Competing viruses are also explored for an SIR model in [33]. Additionally, recently multiple competing viruses, or multi-virus models, have been explored in [34, 35, 36]. A centralized control technique for multi-virus systems is explored in [35, 36].

Various control strategies have been explored for the single-virus model [20, 37, 38]. In [20], an antidote control technique is proposed. In [37], an optimal vaccination control technique is developed using geometric programming ideas. In [38], a network control scheme is applied to a discretized, linearized version of the single-virus model. To the best of our knowledge, the only control strategies that have been designed for the bi-virus model are the ones in [30] and [32], which use geometric programs and are centralized.

In this paper, we study a distributed continuous-time bi-virus model over directed graphs. The model describes how two competing SIS viruses spread over a network of $n>1$ groups of individuals. By competing we mean that no individual can be infected with the two viruses simultaneously. An individual may be infected with one of the two viruses by individuals in its own group, as well as nearest-neighbor groups. These neighbor relationships among the $n$ groups are described by a directed graph on $n$ vertices with an arc (or a directed edge) from vertex $j$ to vertex $i$ whenever the individuals in group $i$ can be infected by those in group $j$. Thus, the neighbor graph has self-arcs at all $n$ vertices, and the arc direction represents the direction of the contagion. The two viruses may spread over different infection graphs with each infection graph being a spanning subgraph of the neighbor graph. Thus, the neighbor graph is the union of the two infection graphs.222 A spanning subgraph is a subgraph that contains all the vertices of the original graph. The union of two directed graphs with the same vertex set is a directed graph with the same vertex set whose arc set is the union of the arc sets of the two graphs.

For two competing viruses, the SIS model is probably the simplest one, but it has limiting behaviors that are already complicated and challenging to analyze. As we will show, the bi-virus SIS model can predict much richer spreading phenomena compared with the single-virus SIS model. The multi-group bi-virus model can be applied to a number of areas, allowing one to understand dynamics of, for example, two competing products in an economic market [39], two competing memes in a social network [40], and two competing species in an ecological environment [41]. Thus, there is ample motivation to thoroughly understand all possible limiting behaviors of the model, which admit different interpretations in different fields. For example, in an epidemic network, a globally stable healthy state, dominant virus state, and coexisting epidemic state predict the eradication of both of the viruses, the triumph of one virus over the other, and the pervasion of both viruses, respectively. Since different applications may prefer different limiting behaviors, it is of great interest to propose efficient control techniques for promoting and/or precluding the spread process in a competitive environment, preferably in a distributed manner. Apparently, for a bi-virus epidemic network, pervasion of any virus is undesirable, and thus a natural question is how to efficiently attenuate or eliminate epidemic spreading using a distributed controller.

The main contributions of this paper are three-fold. First, we analyze the equilibria of the bi-virus model over directed graphs and their stability for both homogeneous and heterogeneous viruses. Second, we derive a sensitivity result for the nontrivial equilibria with respect to the infection and the healing rates, which demonstrates the effects of the simplest local control (i.e., an individual locally adjusts his/her healing and infection rates). Third, we provide an interesting and surprising impossibility result for a certain type of distributed feedback controller, which reveals why distributed control of (bi-virus) epidemic networks is a challenging problem. All the results are validated by a set of illustrative simulations.

The model and assumptions considered in this paper are more general than those in the existing literature [25, 26, 28, 27, 29, 30, 32] in three aspects. First, the work of [25, 26, 28, 27, 29] only considers undirected graphs, whereas this paper studies directed graphs. Second, the work of [25, 29] only considers homogeneous viruses, whereas this paper studies heterogeneous viruses. Third, the analyses in [25, 26, 28, 27, 29, 30, 32] are limited to local stability of the equilibria, while the majority of the analysis performed in this paper is global. Moreover, this paper explores two possible distributed control techniques for the bi-virus model, an endeavor that has not been widely pursued before.

As we will see shortly, the bi-virus model includes the continuous-time single-virus SIS model as a special case. A byproduct of this paper is thus an analysis of the single-virus system which has been studied earlier in [42, 23, 20]. The analysis herein is performed under weaker assumptions on the infection rates $\beta_{ij}$ and healing rates $\delta_{i}$, therefore generalizing previous results on the single-virus system. Although the system defined in [20] admits the same mathematical expression as the single-virus model considered herein, a key difference between the work of [20] and this paper, other than different physical meaning of the models, is that it is assumed in [20] that $\beta_{ii}=0$, for all $i$, whereas this is not the case here. Another difference is that it is also assumed in [20] that $\beta_{ij}=\beta_{kj}$ for all $i,j,k$ if they are both nonzero. In [42], it is assumed that if $\beta_{ij}>0$, then $\beta_{ji}>0$ for all $i,j$, though they are not necessarily equal. Moreover, in [42, 20, 23], it is assumed that $\delta_{i}>0$ for all $i$. In contrast, this paper does not impose any of the aforementioned assumptions and, thus, considers a more general model than those in [23, 20]. Note that the generalization of allowing $\delta_{i}=0$ includes a susceptible-infected (SI) model.

Some of the material in this paper was presented in preliminary form in [1]; this paper presents a more comprehensive treatment of the work in [1]. Additional contributions of this paper, that are not in [1], include 1) complete proofs of all the results, 2) several extensions, including the establishment of uniqueness of the parallel equilibria in Theorems 6 and 7, 3) viewing the sensitivity analysis from a control perspective, 4) an impossibility result for a certain type of distributed feedback controller, and 5) an in-depth set of simulations, illustrating the results and some unproven phenomena.

The remainder of the paper is organized as follows. The basic properties and assumptions of the system model are given in Section II, and the full probabilistic, $3^{n}$ state model is presented. The system equilibria and their stability are studied in Section III. The sensitivity of the equilibria is investigated in Section IV. An impossibility result for distributed feedback control is provided in Section V. Simulations are given in Section VI. The paper concludes with some remarks in Section VII. The proofs of some assertions in the paper are given in the appendix. In the rest of this section, we introduce some notation and provide a number of preliminary results.

I-A Notation

For any positive integer $n$, we use $[n]$ to denote the set $\{1,2,\ldots,n\}$. We view vectors as column vectors. We use $x^{\prime}$ to denote the transpose of a vector $x$ and, similarly, we use $A^{\prime}$ for the transpose of a matrix $A$. The $i$th entry of a vector $x$ will be denoted by $x_{i}$. The $ij$th entry of a matrix $A$ will be denoted by $a_{ij}$. We use $\mathbf{0}$ and $\mathbf{1}$ to denote the vectors whose entries all equal to [math]’s or $1$’s, respectively, and $I$ to denote the identity matrix, where the dimensions of the vectors and matrices are to be understood from the context. For any vector $x\in\mathbb{R}^{n}$, we use ${\rm diag}(x)$ to denote the $n\times n$ diagonal matrix whose $i$th diagonal entry equals $x_{i}$. For any two sets ${\cal A}$ and ${\cal B}$, we use ${\cal A}\setminus{\cal B}$ to denote the set of elements that are in ${\cal A}$ but not in ${\cal B}$, and ${\cal A}\subset{\cal B}$ to denote equality or a proper subset. The notation $1_{a=b}$ is used as an indicator function which takes value one if $a$ equals $b$ and zero otherwise. For $1_{A=b}$, where $A$ is a matrix and $b$ is a scalar, the result is a binary matrix of the same dimensions as $A$ with entries $1_{a_{ij}=b}$.

For any two real vectors $a,b\in\mathbb{R}^{n}$, we write $a\geq b$ if $a_{i}\geq b_{i}$ for all $i\in[n]$, $a>b$ if $a\geq b$ and $a\neq b$, and $a\gg b$ if $a_{i}>b_{i}$ for all $i\in[n]$. Similarly, for any two real matrices $A,B\in\mathbb{R}^{m\times n}$, we write $A\geq B$ if $a_{ij}\geq b_{ij}$ for all $i\in[m]$ and $j\in[n]$, $A>B$ if $A\geq B$ and $A\neq B$, and $A\gg B$ if $a_{ij}>b_{ij}$ for all $i\in[m]$ and $j\in[n]$.

For a complex number $x$, we use $|x|$ and ${\rm Re}(x)$ to denote its magnitude and real part, respectively. For a real square matrix $M$, we use $\rho(M)$ to denote its spectral radius and $s(M)$ to denote the largest real part among its eigenvalues, i.e.,

[TABLE]

where $\sigma(M)$ denotes the spectrum of $M$.

The sign function of a real number $x$ is defined as follows:

[TABLE]

Note that for any real number $x\neq 0$, $\frac{d|x|}{dx}={\rm sgn}(x)$.

I-B Preliminaries

For any two nonnegative vectors $a$ and $b$ in $\mathbb{R}^{n}$ ($a,b\geq\mathbf{0}$), we say that $a$ and $b$ have the same sign pattern if they have zero entries and positive entries in the same places, i.e., for all $i\in[n]$, $a_{i}=0$ if and only if $b_{i}=0$, and $a_{i}>0$ if and only if $b_{i}>0$. A square matrix is called irreducible if it cannot be permuted to a block upper triangular matrix.

Lemma 1

Suppose that $Mx=y$ where $M\in\mathbb{R}^{n\times n}$ is an irreducible nonnegative matrix and $x,y>\mathbf{0}$ are two vectors in $\mathbb{R}^{n}$. If $x$ has at least one zero entry, then $x$ and $y$ cannot have the same sign pattern. In particular, there exists an index $i\in[n]$ such that $x_{i}=0$ and $y_{i}>0$.

A real square matrix is called Metzler if its off-diagonal entries are all nonnegative. Thus, any nonnegative matrix is Metzler.

Lemma 2

*For any matrix $M$ and any real number $\phi$, if $A:=M-\phi I$, then $\sigma(M)=\sigma(A)+\phi$. *

The proof of this lemma is simple and therefore omitted.

The following results from Chapter 2 of [43] for nonnegative matrices, which also hold for Metzler matrices by Lemma 2, with $\phi=\min\{0,m_{11},...,m_{nn}\}$, will be used in the subsequent analysis.

Lemma 3

(Lemma 2.3 in [43])* Suppose that $M$ is an irreducible Metzler matrix. Then, $s(M)$ is a simple eigenvalue of $M$ and there exists a unique (up to scalar multiple) vector $x\gg\mathbf{0}$ such that $Mx=s(M)x$.*

Lemma 4

(Section 2.1 in [43])* Suppose that $M$ is an irreducible Metzler matrix in $\mathbb{R}^{n\times n}$ and $x>\mathbf{0}$ is a vector in $\mathbb{R}^{n}$. For any $\lambda\in\mathbb{R}$, if $Mx<\lambda x$, then $s(M)<\lambda$. If $Mx=\lambda x$, then $s(M)=\lambda$. If $Mx>\lambda x$, then $s(M)>\lambda$.*

We now introduce the Perron-Frobenius Theorem for irreducible nonnegative matrices.

Lemma 5

(Theorem 2.7 and Lemma 2.4 in [43])* Suppose that $M$ is an irreducible nonnegative matrix. Then,*

$M$* has a positive real eigenvalue equal to its spectral radius, $\rho(M)$.* 2. 2.

$\rho(M)$* is a simple eigenvalue of $M$.* 3. 3.

There is an eigenvector $x\gg\mathbf{0}$ corresponding to $\rho(M)$. 4. 4.

$\rho(M)$* increases when any entry of $M$ increases.* 5. 5.

If $N$ is also an irreducible nonnegative matrix and $M>N$, then $\rho(M)>\rho(N)$.

Proposition 1

Suppose that $\Lambda$ is a negative diagonal matrix in $\mathbb{R}^{n\times n}$ and $N$ is an irreducible nonnegative matrix in $\mathbb{R}^{n\times n}$. Let $M=\Lambda+N$. Then, $s(M)<0$ if and only if $\rho(-\Lambda^{-1}N)<1$, $s(M)=0$ if and only if $\rho(-\Lambda^{-1}N)=1$, and $s(M)>0$ if and only if $\rho(-\Lambda^{-1}N)>1$.

Lemma 6

(Proposition 2 in [44])* Suppose that $M$ is a Metzler matrix such that $s(M)<0$. Then, there exists a positive diagonal matrix $P$ such that $M^{\prime}P+PM$ is negative definite.*

Lemma 7

(Lemma A.1 in [20])* Suppose that $M$ is an irreducible Metzler matrix such that $s(M)=0$. Then, there exists a positive diagonal matrix $P$ such that $M^{\prime}P+PM$ is negative semi-definite.*

II The Bi-Virus Model

We are interested in the following continuous-time distributed model for two competing viruses. Consider a network consisting of $n>1$ groups of individuals, labeled $1$ to $n$. There are two competing viruses spreading over the network. An individual cannot be infected with both viruses simultaneously. An individual may be infected with one of the viruses, only by those in its own and neighboring groups. Neighbor relationships among the $n$ groups are described by a directed graph $\mathbb{G}$ on $n$ vertices with an arc from vertex $j$ to vertex $i$ whenever the individuals in group $i$ can be infected by those in group $j$. Thus, the neighbor graph $\mathbb{G}$ has self-arcs at all $n$ vertices and the directions of arcs in $\mathbb{G}$ represent the directions of contagion. Each virus spreads over a spanning subgraph of $\mathbb{G}$. The two subgraphs can be different. Their union is the neighbor graph $\mathbb{G}$. It will be assumed that the two subgraphs are strongly connected and, thus, so is $\mathbb{G}$.333 A directed graph is strongly connected if for every pair of distinct vertices $i$ and $j$, there is a directed path from $i$ to $j$ in the graph.

Let $S_{i}(t)$ denote the number of susceptible individuals in group $i$ at time $t\geq 0$, and let $I^{1}_{i}(t)$ and $I^{2}_{i}(t)$ respectively denote the number of individuals infected with virus 1 and virus 2 in group $i$ at time $t\geq 0$. Assume that the total number of individuals in each group $i$, denoted by $N_{i}$, does not change over time. In other words, $S_{i}(t)+I^{1}_{i}(t)+I^{2}_{i}(t)=N_{i}$, for all $i\in[n]$ and $t\geq 0$. Several parameters are associated with each group $i$: healing rates $\gamma^{1}_{i}$ and $\gamma^{2}_{i}$ for virus 1 and virus 2 respectively, birth rate $\mu_{i}$, death rate $\bar{\mu}_{i}$, and infection rates $\alpha^{1}_{ij}$ and $\alpha^{2}_{ij}$ for virus 1 and virus 2 respectively, $i,j\in[n]$. Since $N_{i}$ is constant, $\bar{\mu}_{i}=\mu_{i}$. We assume that individuals are susceptible at birth even if their parents are infected. The evolution of the number of infected and susceptible individuals in each group $i$ is as follows:

[TABLE]

[TABLE]

[TABLE]

where the infection of group $i$ is caused by one of its neighboring groups $j$, proportional to the total number of infected individuals in group $j$ and the proportion of susceptible individuals in group $i$, which can be regarded as the probability of contact. It is worth noting that, since group $i$ is defined as a neighboring group of itself, its infected individuals can infect its own susceptible members, which reflects realistic scenarios. It is easy to see that since $S_{i}(t)+I^{1}_{i}(t)+I^{2}_{i}(t)=N_{i}$ is a constant, (1) can be implied by (5) and (7). To simplify the model, define the proportion of infected individuals in group $i$ by

[TABLE]

and let

[TABLE]

From (5) and (7), it follows that

[TABLE]

The progression from (7) to the second equation in (9) is given as follows. Dividing both sides of (7) by $N_{i}$, we have

[TABLE]

Since $S_{i}(t)+I^{1}_{i}(t)+I^{2}_{i}(t)=N_{i}$, we have

[TABLE]

Since

[TABLE]

it follows that

[TABLE]

which is the same as the second equation in (9). The first equation in (9) can be derived from (5) in the same way.

The above derivation generalizes the one in [23] for a single virus to two competing viruses. See [23] for a more detailed explanation of the derivation.

Note that each virus has its own non-symmetric infection rates $\beta^{1}_{ij},\beta^{2}_{ij}$ incorporating the nearest-neighbor graph structures and healing rates $\delta^{1}_{i},\delta^{2}_{i}$. A graphical depiction of this model is given in Figure 1. The model can be written in a matrix form as

[TABLE]

where ${x}^{k}(t)\in[0,1]^{n}$, $B^{k}$ is the matrix of $\beta^{k}_{ij}$’s, $X^{k}(t)={\rm diag}(x^{k}(t))$, and $D^{k}={\rm diag}(\delta^{k})$, with $k=1,2$ indicating virus $1$ or $2$.

The same mathematical model was first proposed in [25] with an alternative interpretation in which the system consists of $n$ agents, and ${x}^{1}_{i}(t)$ and ${x}^{2}_{i}(t)$ are the probabilities that agent $i$ has viruses $1$ and $2$, respectively. The model can be viewed as a simplified model resulting from a mean field approximation on a $3^{n}$ state Markov chain model, similar to what has been done for the single-virus SIS model in [14].

For completeness, we provide here a full description of the $3^{n}$ state Markov model. Each state of the chain, $Y_{k}(t)$, corresponds to a ternary-valued string $s$ of length $n$, where $s_{i}=0$, or $s_{i}=1$, or $s_{i}=2$ indicate that the $i$th agent is either susceptible, or infected with virus 1, or infected with virus 2, respectively. The state transition matrix, $Q$, is defined by

[TABLE]

for $i\in[n]$. Here virus 1 and virus 2 are propagating over a network whose infection rates are given by $\beta^{1}_{ij}$ and $\beta^{2}_{ij}$, respectively (nonnegative with $\beta^{1}_{ii}=\beta^{2}_{ii}=0,\ \forall j$), $\delta_{i}^{1}$ and $\delta_{i}^{1}$ are the respective healing rates of the $i$th agent, and, again, $s_{i}=0$, or $s_{i}=1$, or $s_{i}=2$ indicate that the $i$th agent is either susceptible, or infected with virus 1, or infected with virus 2, respectively. The state vector $y(t)$ is defined as

[TABLE]

with $\sum_{k=1}^{3^{n}}y_{k}(t)=1$. The Markov chain evolves as

[TABLE]

See Figure 2 for an illustration of this chain with $n=2$.

Let $v^{1}_{i}(t)=Pr[X_{i}(t)=1]$ and $v^{2}_{i}(t)=Pr[X_{i}(t)=2]$, where $X_{i}(t)$ is the random variables representing whether the $i$th agent is susceptible or infected with virus 1 or 2. Then

[TABLE]

where the $i$th columns of $M^{1}$ and $M^{2}$ indicate the states in the Markov chain where agent $i$ is infected with virus 1 and 2 (all the ternary strings where $s_{i}=1$ and $s_{i}=2$), respectively, that is,

[TABLE]

where $M\in\mathbb{R}^{3^{n}\times n}$ has rows of lexicographically-ordered ternary numbers, bit reversed.444Matlab code: $M=fliplr(dec2base(0:(3^{n})-1,3)-^{\prime}0^{\prime})$ Therefore, $v^{1}_{i}(t)$ and $v^{2}_{i}(t)$ reflect the summation of all probabilities where $s_{i}=1$ and $s_{i}=2$. Note that the first state of the chain, which corresponds to $s=\mathbf{0}$, the healthy state, for $\delta^{1}_{i},\delta^{2}_{i}>0\ \forall i$, is the absorbing, or sink, state of the chain. This means that the Markov chain will never escape the state once in it, and further, since it is the only absorbing state the system will converge to the healthy state with probability one [45].

In [23], a traditional single-group deterministic SIS model was generalized to a networked multi-group setting, similar to the process in (1)-(9), which results in the same model as the mean-field approximation of the networked Markov chain model proposed in [14], except that $\beta_{ii}$ can be nonzero in the deterministic case, which is not possible in the probabilistic derivation, and therefore is more general. Similar to [14], the model in (10) can derived as a first order approximation of the $3^{n}$ state Markov chain model in (11)-(14), which is the model studied in [32]. Consequently, the states of (10) can be interpreted as the probability of the agents being infected or the proportion of subpopulations that are infected. Since using the deterministic bi-virus model for a networked multi-group setting allows $\beta_{ii}$ to be nonzero and is more general, we focus here on this model. For completeness, to illustrate the effectiveness of the first order approximation, we compare (10) and (11)-(14) via simulations in Section VI.

Note that if ${x}^{2}(t)=\mathbf{0}$ for all $t$, (9) reduces to the single-virus model (we drop the superscript since there is only one virus),

[TABLE]

where $x_{i}(t)$ is the proportion of infected individuals in group $i$ (or the probability that agent $i$ has the virus), $\beta_{ij}$’s are the infection rates, $\delta_{i}$’s are the healing rates, and $x_{i}(0)\in[0,1],i\in[n]$, or in matrix form

[TABLE]

Consider a further special case in which all $n$ groups are isolated, i.e., $\beta_{ij}=0$ for all $i\neq j$. For each group $i$, if $\beta_{ii}>0$, from (15),

[TABLE]

which is an SIS model for a single group of individuals. It is not hard to verify that for this single-group SIS model, if $\delta_{i}\geq\beta_{ii}$, $x_{i}(t)$ will always converge to zero, and if $\delta_{i}<\beta_{ii}$, $x_{i}(t)$ will always converge to a positive value unless $x_{i}(0)=0$. If $\beta_{ii}=0$, then the dynamics reduce to $\dot{x}_{i}(t)=-\delta_{i}x_{i}(t)$ which is not an epidemic model. Therefore, $\beta_{ii}>0$ has an important physical meaning in networked models of groups of individuals (sometimes called metapopulation models in the literature). If we interpret each state of the model to be the infection probability of a single agent, and we further assume that $\beta_{ij}$ can be factored into $\beta_{i}a_{ij}$, where $\beta_{i}$ is the infection rate of agent $i$, and $a_{ij}$ is the connection structure between agents, then (15) becomes the single SIS model proposed in [14].

We impose the following assumptions on the model throughout Sections III and IV.

Assumption 1

For all $i\in[n]$, we have $x^{1}_{i}(0),x^{2}_{i}(0),(1-x^{1}_{i}(0)-x^{2}_{i}(0))\in[0,1]$.

Assumption 2

For all $i\in[n]$, we have $\delta^{1}_{i},\delta^{2}_{i}\geq 0$. The matrices $B^{1}$ and $B^{2}$ are nonnegative and irreducible.

Assumption 1 says that the initial proportions of infected and healthy individuals are in the interval $[0,1]$. The nonnegativity assumption on the matrix $B^{k}$ is equivalent to $\beta^{k}_{ij}\geq 0$ for all $k\in[2]$ and $i,j\in[n]$. Assumption 2 says that all healing and infection rates are nonnegative. The assumption of an irreducible matrix $B^{k}$ is equivalent to a strongly connected spreading graph for virus $k$, $k\in[2]$.

Lemma 8

Under the conditions of Assumptions 1 and 2, $x^{1}_{i}(t),x^{2}_{i}(t),x^{1}_{i}(t)+x^{2}_{i}(t)\in[0,1]$ for all $i\in[n]$ and $t\geq 0$.

Lemma 8 implies that the set

[TABLE]

is positively invariant with respect to the system defined by (10). Since $x^{1}_{i}$ and $x^{2}_{i}$ denote the fractions of group $i$ infected by viruses 1 and 2, respectively, and $1-x^{1}_{i}-x^{2}_{i}$ denotes the fraction of group $i$ that is healthy, it is natural to assume that their initial values are in the interval $[0,1]$, since otherwise the values will lack any physical meaning for the epidemic model considered here. Therefore, in this paper, we focus on the analysis of (10) only on the domain $\mathcal{D}$, as defined in (17).

We are interested in the problem of characterizing limiting behavior of the bi-virus model (10) and its dependence on the network structure (two spreading graphs) and parameters (healing rates $\delta_{i}^{1}$, $\delta_{i}^{2}$ and infection rates $\beta_{ij}^{1}$, $\beta_{ij}^{2}$). The limiting behaviors will be characterized by the equilibria of the system and their stability. The effects of the network structure and parameters on the limiting behavior are important and useful for controlling the epidemic spreading process. Specifically, it will be shown in the following sections that, under certain assumptions, whether the two viruses eventually disappear or not can be determined by checking the network structure and parameters. In the case when at least one virus ultimately spreads over the network, the spreading process can be attenuated or eliminated by modifying the network parameters.

III Equilibria and Their Stability

In this section, we analyze the equilibria of the system (10) and their stability, which characterize limiting behaviors of the bi-virus model.

First, it can be seen that $x^{1}=x^{2}=\mathbf{0}$ is an equilibrium of the system (10), which corresponds to the case when no individual is infected. We call this trivial equilibrium the healthy state. We will show that (10) also admits nonzero equilibria under appropriate assumptions. We call those nonzero equilibria epidemic states. In this section, we study the stability of the healthy state as well as the epidemic states of (10). To state our results, we need the following definition.

Definition 1

Consider an autonomous system

[TABLE]

where $f:{\cal X}\rightarrow\mathbb{R}^{n}$ is a locally Lipschitz map from a domain ${\cal X}\subset\mathbb{R}^{n}$ into $\mathbb{R}^{n}$. Let $z$ be an equilibrium of (18) and ${\cal E}\subset{\cal X}$ be a domain containing $z$. If the equilibrium $z$ is asymptotically stable such that for any $x(0)\in{\cal E}$ we have $\lim_{t\rightarrow\infty}x(t)=z$, then ${\cal E}$ is said to be a domain of attraction for $z$.

The following corollary is a direct consequence of Lyapunov’s stability theorem (see Theorem 4.1 in [46]) and the definition of domain of attraction.

Corollary 1

Let $z$ be an equilibrium of (18) and ${\cal E}\subset{\cal X}$ be a domain containing $z$. Let $V:{\cal E}\rightarrow\mathbb{R}$ be a continuously differentiable function such that $V(z)=0$, $V(x)>0$ in ${\cal E}\setminus\{z\}$, $\dot{V}(z)=0$, and $\dot{V}(x)<0$ in ${\cal E}\setminus\{z\}$. If ${\cal E}$ is a positively invariant set, then the equilibrium $z$ is asymptotically stable with a domain of attraction ${\cal E}$.

The following theorem establishes a sufficient condition for global stability of the healthy state, whereas the work of [30, 32] provides necessary and sufficient conditions for local stability of the healthy state.

Theorem 1

Under Assumption 2, the healthy state is the unique equilibrium of system (10) if, and only if, $s(-D^{1}+B^{1})\leq 0$ and $s(-D^{2}+B^{2})\leq 0$. Furthermore, in this case, the healthy state is asymptotically stable with domain of attraction $\mathcal{D}$, as defined in (17).

To prove the theorem, we need the following result for the single-virus model (16).

Proposition 2

Suppose that $\delta_{i}\geq 0$ for all $i\in[n]$ and that matrix $B$ is nonnegative and irreducible. If $s(-D+B)\leq 0$, then $\mathbf{0}$ is the unique equilibrium of system (16), which is asymptotically stable with domain of attraction $[0,1]^{n}$. If $s(-D+B)>0$, then system (16) has two equilibria, $\mathbf{0}$ and $x^{*}$ which satisfies $x^{*}\gg\mathbf{0}$.

This result has been proved in [42, 23, 20] for the case when $\delta_{i}>0$ for all $i\in[n]$. Additional assumptions on $\beta_{ij}$ are also imposed in [42, 20]. Using a different Lyapunov function from the one used in the aforementioned papers, we extend the result by allowing $\delta_{i}=0$ (see the Appendix), which reflects the situations where certain groups of individuals are unable to heal themselves (due, for example, to the lack of vaccines). Note that this allows for the SI model in which $\delta_{i}=0$ for all $i\in[n]$. In the case when $\delta_{i}>0$ for all $i\in[n]$, the conditions $s(-D+B)\leq 0$ and $\rho(D^{-1}B)\leq 1$ are equivalent, which is a direct consequence of Proposition 1. In the case when $\delta_{i}=0$ for some but not all $i\in[n]$, the two conditions are not equivalent since $\rho(D^{-1}B)$ is not well-posed in this case.

Proof of Theorem 1: We first show that both $x^{1}(t)$ and $x^{2}(t)$ will asymptotically converge to $\mathbf{0}$ as $t\rightarrow\infty$ for any initial condition in $\mathcal{D}$. Since $x^{1}_{i}(t)$ and $x^{2}_{i}(t)$ are always nonnegative by Lemma 8, from (9), we obtain

[TABLE]

which imply that each of the trajectories of $x^{1}_{i}(t)$ and $x^{2}_{i}(t)$ is bounded above by a single-virus model. From Assumption 2 and Proposition 2, both $x^{1}_{i}(t)$ and $x^{2}_{i}(t)$ will asymptotically converge to $\mathbf{0}$ as $t\rightarrow\infty$, and thus the healthy state is the unique equilibrium of (10).

We next show the asymptotic stability of the healthy state. Consider the Lyapunov function candidate $V(x^{1}(t),x^{2}(t))=x^{1}(t)^{\prime}P^{1}x^{1}(t)+x^{2}(t)^{\prime}P^{2}x^{2}(t)$, where $P^{1}$ and $P^{2}$ are positive diagonal matrices chosen in the same way as $P$ is chosen in Proposition 2. Note from (10) that

[TABLE]

Using similar arguments to those in the proof of Proposition 2, it can be shown that $\dot{V}(x^{1}(t),x^{2}(t))$ is a negative definite function over ${\cal D}$ except for the healthy state. By Lemma 8 and Corollary 1, the healthy state is asymptotically stable with domain of attraction ${\cal D}$.

Now we show that if either $s(-D^{1}+B^{1})>0$ or $s(-D^{2}+B^{2})>0$, then system (10) has an epidemic state.

Without loss of generality, suppose that $s(-D^{1}+B^{1})>0$. Set $x^{2}=\mathbf{0}$. Then, the dynamics of $x^{1}$ simplifies to that of the single-virus system, which admits an epidemic state by Proposition 2. Therefore, in the case when $s(-D^{1}+B^{1})>0$, the system (10) always admits an equilibrium of the form $(\tilde{x}^{1},\mathbf{0})$ with $\tilde{x}^{1}\gg\mathbf{0}$.

We have provided a necessary and sufficient condition for the eradication of both of the viruses, i.e., $s(-D^{1}+B^{1})\leq 0$ and $s(-D^{2}+B^{2})\leq 0$. Since larger nonzero entries of $D^{1}$ and $D^{2}$ (i.e., healing rates $\delta_{i}^{1}$, $\delta_{i}^{2}$) will decrease the two quantities, and larger nonzero entries of $B^{1}$ and $B^{2}$ (i.e., infection rates $\beta_{ij}^{1}$, $\beta_{ij}^{2}$) will increase the two quantities, the necessary and sufficient condition can be interpreted as the overall healing capabilities of all the individuals overcoming or balancing out competely the effects of the network infection. This characterization is important for understanding when a system will become completely healthy as illustrated via simulation in Figure 7.

Now we turn to the analysis of epidemic states. We begin with dominant virus states at which one virus is eradicated and the other one pervades in the network.

Theorem 2

Under Assumption 2, if $s(-D^{1}+B^{1})>0$ and $s(-D^{2}+B^{2})\leq 0$, then (10) has two equilibria, the healthy state $(\mathbf{0},\mathbf{0})$, where the system converges to this equilibrium for all initial conditions $(x^{1}(0),x^{2}(0))\in\{(\mathbf{0},x^{2})|x^{2}\in[0,1]^{n}\}$, and a unique epidemic state of the form $(\tilde{x}^{1},\mathbf{0})$ with $\tilde{x}^{1}\gg\mathbf{0}$, which is asymptotically stable with domain of attraction $\mathcal{D}\setminus\{(\mathbf{0},x^{2})|x^{2}\in[0,1]^{n}\}$, with $\mathcal{D}$ defined in (17).

Remark 1

Note that healthy state $(\mathbf{0},\mathbf{0})$ is an unstable equilibrium. A small perturbation of the first virus from the origin will drive the system to the unique epidemic state $(\tilde{x}^{1},\mathbf{0})$. $\Box$

To prove Theorem 2, we need the following result for the single-virus model (16), which builds on Proposition 2.

Proposition 3

Consider the single-virus model (16). Suppose that $\delta_{i}\geq 0$ for all $i\in[n]$, and that the matrix $B$ is nonnegative and irreducible. If $s(-D+B)>0$, then the epidemic state $x^{*}\gg\mathbf{0}$ is asymptotically stable with domain of attraction $[0,1]^{n}\setminus\{\mathbf{0}\}$.

This result has been proved in [23, 20] for the case when $\delta_{i}>0$ for all $i\in[n]$. We extend the result by allowing $\delta_{i}=0$. The analyses in [23, 20] cannot be applied here.

To prove Proposition 3, we need the following lemma.

Lemma 9

Consider the single-virus model (16). Suppose that $\delta_{i}\geq 0$ for all $i\in[n]$, and that the matrix $B$ is nonnegative and irreducible. If $x(0)>\mathbf{0}$, then there exists a $\tau\geq 0$ such that $x(\tau)\gg\mathbf{0}$.

Proofs of these results are provided in the appendix.

We are now in a position to provide the proof for Theorem 2.

Proof of Theorem 2: By Proposition 2, $x^{2}(t)$ will asymptotically converge to $\mathbf{0}$ as $t\rightarrow\infty$ for all initial values $(x^{1}(0),x^{2}(0))\in\{(\mathbf{0},x^{2})|x^{2}\in[0,1]^{n}\}$, since the system reduces to the single virus case for that set of initial conditions.

From (10), we have

[TABLE]

Thus, we can regard the dynamics of $x^{1}(t)$ as an autonomous system

[TABLE]

with a vanishing perturbation $-X^{2}(t)B^{1}x^{1}(t)$, which converges to $\mathbf{0}$ as $t\rightarrow\infty$. From Proposition 3, the autonomous system (19) will asymptotically converge to a unique epidemic state $\tilde{x}^{1}\gg\mathbf{0}$ for any $x^{1}(0)\in[0,1]^{n}\setminus\{\mathbf{0}\}$.

Let $y^{1}_{i}(t)=x^{1}_{i}(t)-\tilde{x}^{1}_{i}$ for all $i\in[n]$, or equivalently, $y^{1}(t)=x^{1}(t)-\tilde{x}^{1}$. Then,

[TABLE]

Let $f(y^{1}(t))=(-D^{1}+(I-\tilde{X}^{1})B^{1}-{\rm diag}(B^{1}x^{1}(t)))y^{1}(t)$ and $g(t,y^{1}(t))=-X^{2}(t)B^{1}(y^{1}(t)+\tilde{x}^{1})$. Consider the Lyapunov function candidate

[TABLE]

Then, $\dot{V}(y^{1}(t))=\frac{\partial V}{\partial y^{1}}f(y^{1}(t))+\frac{\partial V}{\partial y^{1}}g(t,y^{1}(t))$. From the proof of Proposition 3, $\frac{\partial V}{\partial t}+\frac{\partial V}{\partial y^{1}}f(t,y^{1})<0$ unless $y^{1}(t)=\mathbf{0}$ (i.e., $x^{1}(t)=\tilde{x}^{1}$). Since $x^{2}(t)$ asymptotically converges to $\mathbf{0}$, so does $\frac{\partial V}{\partial y^{1}}g(t,y^{1}(t))$. This implies that after a sufficiently long time, $\dot{V}(y^{1}(t))<0$ if $x^{1}(t)$ does not equal $\tilde{x}^{1}$. Using the same argument as in the proof of Proposition 3, $(x^{1}(t),x^{2}(t))$ will asymptotically converge to the unique epidemic state $(\tilde{x}^{1},\mathbf{0})$ for any $(x^{1}(0),x^{2}(0))\in\mathcal{D}\setminus\{(\mathbf{0},x^{2})|x^{2}\in[0,1]^{n}\}$, with $\mathcal{D}$ defined in (17).

Theorem 2 provides conditions under which one virus will pervade the network, and the other one will be driven out. Understanding this condition is useful for characterizing when a designer (marketer, politician, etc.) will consistently dominate a competitor. This behavior is illustrated in Figure 9.

It is clear from the preceding results that as long as one of $s(-D^{k}+B^{k})$, $k\in\{1,2\}$, is less than or equal to zero, at most one virus will ultimately spread over the network. A natural question is thus whether the two viruses can coexist when $s(-D^{k}+B^{k})$, $k\in\{1,2\}$, are both larger than zero. In the following, we will partially answer this question. We begin with a result regarding non-coexisting equilibria.

Let $(\tilde{x}^{1},\tilde{x}^{2})$ be an equilibrium of (10). Here, both $\tilde{x}^{1}$ and $\tilde{x}^{2}$ can be $\mathbf{0}$. Then, the Jacobian matrix of the equilibrium, denoted by $J(\tilde{x}^{1},\tilde{x}^{2})$, with $\tilde{B}^{i}={\rm diag}(B^{i}\tilde{x}^{i})$, $i\in[2]$, is

[TABLE]

Theorem 3

Under Assumption 2, if $s(-D^{1}+B^{1})>0$ and $s(-D^{2}+B^{2})>0$, then (10) has at least three equilibria, the healthy state $(\mathbf{0},\mathbf{0})$, and two epidemic states of the form $(\tilde{x}^{1},\mathbf{0})$ with $\tilde{x}^{1}\gg\mathbf{0}$ and $(\mathbf{0},\tilde{x}^{2})$ with $\tilde{x}^{2}\gg\mathbf{0}$. The healthy state $(\mathbf{0},\mathbf{0})$ is unstable.

Proof: The existence of the two epidemic states is an immediate consequence of Proposition 2. The healthy state $(\mathbf{0},\mathbf{0})$ is always an equilibrium of (10). Since by (20)

[TABLE]

which is unstable as $s(-D^{1}+B^{1})>0$ and $s(-D^{2}+B^{2})>0$, the healthy state $(\mathbf{0},\mathbf{0})$ is unstable.

Analysis of any of the other possible equilibria is challenging for general cases. In the sequel, we will thus consider two special cases, characterized by some additional assumptions on the underlying graphs for the spread of viruses. It turns out that if for both $k\in\{1,2\}$, $s(-D^{k}+B^{k})$ are larger than zero, then the existence of coexisting equilibria is not guaranteed, as shown in the following special case, where two homogeneous viruses spread on the same graph.

Assumption 3

Viruses 1 and 2 spread over the same strongly connected directed graph $\mathbb{G}=([n],{\cal E})$, with $\delta^{1}_{i}=\delta^{1}>0$ and $\delta^{2}_{i}=\delta^{2}>0$ for all $i\in[n]$, and $\beta^{1}_{ij}=\beta^{1}>0$ and $\beta^{2}_{ij}=\beta^{2}>0$ for all $i=j\in[n]$ and $(i,j)\in{\cal E}$.

Under Assumption 3, it should be clear that $D^{1}=\delta^{1}I$, $D^{2}=\delta^{2}I$, $B^{1}=\beta^{1}A$, and $B^{2}=\beta^{2}A$, where $A$ is the adjacency matrix of $\mathbb{G}$, which is an irreducible Metzler matrix.

Theorem 4

Suppose that Assumptions 1 and 3 hold. Then, coexisting equilibria exist only if $\frac{\delta^{1}}{\beta^{1}}=\frac{\delta^{2}}{\beta^{2}}$.

This result has been proved in [25] for the case when $\mathbb{G}$ is an undirected graph. We extend the result by allowing $\mathbb{G}$ to be directed. To prove the theorem, we need the following lemma.

Lemma 10

Suppose that Assumptions 1 and 3 hold. If $(\tilde{x}^{1},\tilde{x}^{2})$ is an equilibrium of (10), then $\tilde{x}^{1}+\tilde{x}^{2}\ll\mathbf{1}$.

Proof of Theorem 4: To prove the theorem, suppose that, to the contrary, there exists an equilibrium $(\tilde{x}^{1},\tilde{x}^{2})$ such that $\tilde{x}^{1},\tilde{x}^{2}>\mathbf{0}$ in the case when $\frac{\delta^{1}}{\beta^{1}}\neq\frac{\delta^{2}}{\beta^{2}}$. From (10) and Assumption 3,

[TABLE]

From Lemma 10, $(I-\tilde{X}^{1}-\tilde{X}^{2})$ is a positive diagonal matrix, and thus $(I-\tilde{X}^{1}-\tilde{X}^{2})A$ is also an irreducible Metzler matrix. Since $\tilde{x}^{1},\tilde{x}^{2}>\mathbf{0}$, from Lemma 4, $s((I-\tilde{X}^{1}-\tilde{X}^{2})A)=\frac{\delta^{1}}{\beta^{1}}=\frac{\delta^{2}}{\beta^{2}}$, which contradicts the hypothesis that $\frac{\delta^{1}}{\beta^{1}}\neq\frac{\delta^{2}}{\beta^{2}}$. Therefore, coexisting equilibria may exist only if $\frac{\delta^{1}}{\beta^{1}}=\frac{\delta^{2}}{\beta^{2}}$.

For the following, without loss of generality, we assume $\frac{\delta^{1}}{\beta^{1}}>\frac{\delta^{2}}{\beta^{2}}$.

Theorem 5

Suppose that Assumptions 1 and 3 hold and that $s(A)>\frac{\delta^{1}}{\beta^{1}}>\frac{\delta^{2}}{\beta^{2}}$. Then, system (10) has three equilibria, the healthy state $(\mathbf{0},\mathbf{0})$ which is unstable, $(\tilde{x}^{1},\mathbf{0})$ with $\tilde{x}^{1}\gg\mathbf{0}$ which is unstable, and $(\mathbf{0},\tilde{x}^{2})$ with $\tilde{x}^{2}\gg\mathbf{0}$ which is locally exponentially stable.

This result has been proved in [25] for the case when $\mathbb{G}$ is an undirected graph. We extend the result by allowing $\mathbb{G}$ to be directed with a proof technique similar to [25]. In [27], a sufficient condition is established for the case in which $\delta^{1}=\delta^{2}=1$ and $\beta_{ij}^{k}$, $k\in\{1,2\}$, are heterogeneous (see Corollary 4 in [27]).

Proof: From Theorem 4, the system (10) cannot have any equilibria of the form $(\tilde{x}^{1},\tilde{x}^{2})$ with $\tilde{x}^{1},\tilde{x}^{2}>\mathbf{0}$. Thus, if $(\tilde{x}^{1},\tilde{x}^{2})$ is an equilibrium of (10), at least one of $\tilde{x}^{1}$ and $\tilde{x}^{2}$ equals $\mathbf{0}$. It is clear that $(\mathbf{0},\mathbf{0})$ is always an equilibrium. Suppose that $\tilde{x}^{1}=\mathbf{0}$ and $\tilde{x}^{2}>\mathbf{0}$. Then, from Proposition 2, $\tilde{x}^{2}\gg\mathbf{0}$ and is unique. Similarly, when $\tilde{x}^{1}>\mathbf{0}$ and $\tilde{x}^{2}=\mathbf{0}$, $\tilde{x}^{1}\gg\mathbf{0}$ and is unique. Thus, the system (10) has exactly three equilibria.

Next we turn to the stability of the three equilibria. Note that from Assumption 3, the hypothesis $s(A)>\frac{\delta^{1}}{\beta^{1}}>\frac{\delta^{2}}{\beta^{2}}$ implies $s(-D^{1}+B^{1}),s(-D^{2}+B^{2})>0$. Then, from Theorem 3, the healthy state $(\mathbf{0},\mathbf{0})$ is unstable.

From (20), the Jacobian at $(\tilde{x}^{1},\mathbf{0})$ equals

[TABLE]

From (10) and Assumption 3, $(I-\tilde{X}^{1})A\tilde{x}^{1}=\frac{\delta^{1}}{\beta^{1}}\tilde{x}^{1}$. It follows from Lemma 10 that $(I-\tilde{X}^{1})A$ is an irreducible Metzler matrix. Then, from Lemma 4, $s((I-\tilde{X}^{1})A)=\frac{\delta^{1}}{\beta^{1}}$. Since $\frac{\delta^{1}}{\beta^{1}}>\frac{\delta^{2}}{\beta^{2}}$, it follows that

[TABLE]

which implies that the Jacobian matrix is unstable. Thus, the equilibrium $(\tilde{x}^{1},\mathbf{0})$ with $\tilde{x}^{1}\gg\mathbf{0}$ is unstable.

From (20), the Jacobian at $(\mathbf{0},\tilde{x}^{2})$ equals

[TABLE]

Using the same arguments as in the previous paragraph, $s(\beta^{1}(I-\tilde{X}^{2})A-\delta^{1}I)<0$. From (10) and Assumption 3,

[TABLE]

Since $\tilde{x}^{2}\gg\mathbf{0}$ and $A$ is irreducible, it must be true that

[TABLE]

It follows from Assumptions 1 and 3 and Lemma 8 that $\beta^{2}(I-\tilde{X}^{2})A-\delta^{2}I-\beta^{2}{\rm diag}(A\tilde{x}^{2})$ is an irreducible Metzler matrix. Then, from Lemma 4, $s(\beta^{2}(I-\tilde{X}^{2})A-\delta^{2}I-\beta^{2}{\rm diag}(A\tilde{x}^{2}))<0$, which implies that the Jacobian matrix is stable. Thus, the equilibrium $(\mathbf{0},\tilde{x}^{2})$ with $\tilde{x}^{2}\gg\mathbf{0}$ is locally exponentially stable.

For the possibility of coexisting equilibria, we have the following interesting result.

Theorem 6

Suppose that Assumptions 1 and 3 hold and that $s(A)>\frac{\delta^{1}}{\beta^{1}}=\frac{\delta^{2}}{\beta^{2}}$. If $(\tilde{x}^{1},\tilde{x}^{2})$ with $\tilde{x}^{1},\tilde{x}^{2}>\mathbf{0}$ is an equilibrium of (10), then $\tilde{x}^{1},\tilde{x}^{2}\gg\mathbf{0}$ and $\tilde{x}^{1}=\alpha\tilde{x}^{2}$ for some constant $\alpha>0$. Furthermore, for each $\alpha>0$ there exists a unique pair $(\tilde{x}^{1},\tilde{x}^{2})$ such that $\tilde{x}^{1}=\alpha\tilde{x}^{2}$.

Proof: From the proof of Theorem 4,

[TABLE]

in which $(I-\tilde{X}^{1}-\tilde{X}^{2})A$ is an irreducible Metzler matrix. From Lemma 3, it must be true that $\tilde{x}^{1},\tilde{x}^{2}\gg\mathbf{0}$ and $\tilde{x}^{1}=\alpha\tilde{x}^{2}$ for some constant $\alpha>0$.

Given some $\alpha>0$, assume $\tilde{x}^{1}=\alpha\tilde{x}^{2}$ and $\hat{x}^{1}=\alpha\hat{x}^{2}$ both satisfy (21). Therefore,

[TABLE]

so that, by Lemma 4,

[TABLE]

Also, without loss of generality, assume there exists $j\in[n]$ such that $\tilde{x}^{1}_{j}>\hat{x}^{1}_{j}$ and $\tilde{x}^{1}_{i}=\hat{x}^{1}_{i}$ for all $i\neq j$. This implies

[TABLE]

which by Lemma 4 implies

[TABLE]

However, this contradicts (22). So for each $\alpha>0$ there exists a unique pair $(\tilde{x}^{1},\tilde{x}^{2})$ such that $\tilde{x}^{1}=\alpha\tilde{x}^{2}$.

Remark 2

Note, from (20) and (21), that when $\frac{\delta^{1}}{\beta^{1}}=\frac{\delta^{2}}{\beta^{2}}$,

[TABLE]

i.e., the Jacobian matrix has a zero eigenvalue. Therefore, linearization says nothing about the local stability of the coexisting equilibria. Simulations indicate that, depending on the initial condition, the system can arrive at different equilibria of the form $\tilde{x}^{1}=\alpha\tilde{x}^{2}$ for different constants $\alpha>0$. $\Box$

Two viruses spreading on the same graph can be thought of as two products spreading in a market or two competing ideas spreading on a social network. Developing an understanding of how the viruses can coexist, and that the equilibrium that is reached is dependent on the initial condition, is vital to deploying initial marketing strategies that result in different market shares. This coexistence behavior is illustrated via simulation in Figures 11 and 12.

A similar result can be established for another special case, where two identical heterogeneous viruses spread on the same graph, as specified by the following assumption.

Assumption 4

Viruses 1 and 2 spread over the same strongly connected directed graph $\mathbb{G}=([n],{\cal E})$, with $\delta^{1}_{i}=\delta^{2}_{i}>0$ for all $i\in[n]$, and $\beta^{1}_{ij}=\beta^{2}_{ij}$ for all $i=j\in[n]$ and $(i,j)\in{\cal E}$.

Under Assumption 4, we have $D^{1}=D^{2}=D$ and $B^{1}=B^{2}=B$, where $D$ is a positive diagonal matrix and $B$ is an irreducible nonnegative matrix.

Theorem 7

Suppose that Assumptions 1 and 4 hold and that $s(-D+B)>0$. If $(\tilde{x}^{1},\tilde{x}^{2})$ with $\tilde{x}^{1},\tilde{x}^{2}>\mathbf{0}$ is an equilibrium of (10), then $\tilde{x}^{1},\tilde{x}^{2}\gg\mathbf{0}$, $\tilde{x}^{1}+\tilde{x}^{2}$ is unique, and $\tilde{x}^{1}=\alpha\tilde{x}^{2}$ for some constant $\alpha>0$. Furthermore, for each $\alpha>0$, there exists a unique pair $(\tilde{x}^{1},\tilde{x}^{2})$ such that $\tilde{x}^{1}=\alpha\tilde{x}^{2}$.

Proof: From (10) and Assumption 4,

[TABLE]

Thus, the dynamics of $x^{1}(t)+x^{2}(t)$ is equivalent to that of the single-virus model (16). By Proposition 2, $x^{1}(t)+x^{2}(t)$ has a unique nonzero equilibrium in $[0,1]^{n}$. Thus, $\tilde{x}^{1}+\tilde{x}^{2}$ is unique. From (10), we have

[TABLE]

Then, $(-D+B-(\tilde{X}^{1}+\tilde{X}^{2})B)(\tilde{x}^{1}-\tilde{x}^{2})=\mathbf{0}$. Using the same arguments as those in the proof of Lemma 10, it can be shown that $\tilde{x}^{1}+\tilde{x}^{2}\ll\mathbf{1}$. Thus, $-D+B-(\tilde{X}^{1}+\tilde{X}^{2})B$ is an irreducible Metzler matrix. By Lemma 4, since

[TABLE]

and $x^{1}(t)+x^{2}(t)\gg\mathbf{0}$, $s(-D+B-(\tilde{X}^{1}+\tilde{X}^{2})B)=0$. From Lemma 3, either $x^{1}(t)=x^{2}(t)$ or $x^{1}(t)-x^{2}(t)=\gamma(x^{1}(t)+x^{2}(t))$ for some constant $\gamma>0$. In both cases, it must be true that $\tilde{x}^{1}=\alpha\tilde{x}^{2}$ for some constant $\alpha>0$, and thus $\tilde{x}^{1},\tilde{x}^{2}\gg\mathbf{0}$.

Since $\tilde{x}^{1}+\tilde{x}^{2}$ is unique and $\tilde{x}^{1}=\alpha\tilde{x}^{2}$ for some constant $\alpha>0$, then $\tilde{x}^{1}+\tilde{x}^{2}=(1+\alpha)\tilde{x}^{2}$ is constant. Therefore, for each $\alpha>0$ there exists a unique $\tilde{x}^{2}$. So for each $\alpha>0$ there exists a unique pair $(\tilde{x}^{1},\tilde{x}^{2})$.

In this section, we have explored the equilibria of the bi-virus model and their local and global stability. It has been shown that unless both $s(-D^{1}+B^{1})$ and $s(-D^{2}+B^{2})$ are less than or equal to zero, at least one virus will pervade the network. Indeed, under different appropriate assumptions, one virus will fully dominate the competition, driving the second virus out, or both viruses will coexist with infinitely many possible equilibria. It is worth noting that the necessary and sufficient condition for the eradication of both viruses (i.e., $s(-D^{1}+B^{1})\leq 0$ and $s(-D^{2}+B^{2})\leq 0$) requires global information on the network, which is inefficient and sometimes impossible. Subsequently, we explore some simple, but local, control techniques, with the expectation of attenuating or eliminating epidemic spreading. We begin with the case when only one virus pervades.

IV Sensitivity

In this section, we regard each healing or infection rate as a local variable that an individual can control. The aim of this section is to understand whether and how local adjustment of healing and infection rates will affect the whole network.

We have already shown that in the case when $s(-D^{1}+B^{1})>0$ and $s(-D^{2}+B^{2})\leq 0$, the system (10) has a unique epidemic state of the form $(\tilde{x}^{1},\mathbf{0})$ with $\tilde{x}^{1}\gg\mathbf{0}$, which is asymptotically stable. It can be seen that the value of $\tilde{x}^{1}$ is independent of the matrices $D^{2}$ and $B^{2}$, but depends on the matrices $D^{1}$ and $B^{1}$, or equivalently, the parameters $\delta^{1}_{i}$ and $\beta^{1}_{ij}$. A natural question is: how does the equilibrium $\tilde{x}^{1}$ change when the values of $\delta^{1}_{i}$ and $\beta^{1}_{ij}$ are perturbed? The aim of this section is to answer this question.

From the proof of Theorem 2, the value of $\tilde{x}^{1}$ equals the unique epidemic state, denoted $x^{*}$, of the single-virus model (16) when $s(-D+B)>0$. Thus, to answer the question just raised, it is equivalent to study how the equilibrium $x^{*}$ changes when the values of $\delta_{i}$ and $\beta_{ij}$ are perturbed.

For our purposes, we assume in this section that $\delta_{i}>0$ for all $i\in[n]$. Then, by Proposition 1, $s(-D+B)>0$ if and only if $\rho(D^{-1}B)>1$.

Suppose that $\rho(D^{-1}B)>1$. By Proposition 2, the epidemic state $x^{*}$ is the unique nonzero equilibrium of (16), which satisfies the equation $(-D+B-X^{*}B)x^{*}=\mathbf{0}$. Define the mapping $\Phi$ as follows:

[TABLE]

Then, the equation $\Phi(x^{*},D,B)=\mathbf{0}$ defines an implicit function $g:\mathbb{R}^{n\times n}\times\mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{n}$ given by $x^{*}=g(D,B)$.

For each pair of matrices $D$ and $B$ for which $\rho(D^{-1}B)>1$, there must exist a small neighborhood ${\cal B}$ such that for any pair of matrices $D+\Delta D$ and $B+\Delta B$ in ${\cal B}$,

[TABLE]

Here $\Delta D$ is the $n\times n$ diagonal matrix whose $i$th diagonal entry equals $\Delta\delta_{i}$, which denotes the perturbation of $\delta_{i}$, and $\Delta B$ is the $n\times n$ matrix whose $ij$th entry equals $\Delta\beta_{ij}$, which denotes the perturbation of $\beta_{ij}$. Let $x^{*}+\Delta x^{*}$ denote the new epidemic state resulting from the perturbations. Then,

[TABLE]

where $\Delta X^{*}={\rm diag}(\Delta x^{*})$. By ignoring the higher-order $\Delta$ terms, it is straightforward to verify that

[TABLE]

where $\Delta\delta$ is the vector in $\mathbb{R}^{n}$ whose $i$th entry equals $\Delta\delta_{i}$ and $\tilde{B}^{*}={\rm diag}(Bx^{*})$. Since we are interested in the local behavior around the equilibrium $x^{*}$, which is equivalent to performing linearization at $x^{*}$, the approximation in (23) is accurate for the following arguments. First note that

[TABLE]

Since $B$ is an irreducible nonnegative matrix and $x^{*}\gg\mathbf{0}$, it follows that $\tilde{B}^{*}$ is a positive diagonal matrix. Let $c>0$ be any positive constant, strictly smaller than the minimal diagonal entry of $\tilde{B}^{*}$. Then, $\tilde{B}^{*}>cI$ and thus $-\tilde{B}^{*}x^{*}<-cx^{*}$. Since $(-D+B-X^{*}B-\tilde{B}^{*})$ is an irreducible Metzler matrix, by Lemma 4, $s(-D+B-X^{*}B-\tilde{B}^{*})<-c<0$, which implies that $(-D+B-X^{*}B-\tilde{B}^{*})$ is nonsingular. Thus, by the Implicit Function Theorem (see, e.g., pages 204-206 in [47]), the function $x^{*}=g(D,B)$ is differentiable in the neighborhood ${\cal B}$. From (23), we have

[TABLE]

To proceed, we need the following lemma.

Lemma 11

(Theorem 2.7 in Chapter 6 of [48])* Suppose that $M$ is a nonsingular, irreducible Hurwitz Metzler matrix. Then, $M^{-1}\ll 0$.*

From this lemma and the preceding discussion, it follows immediately that $(-D+B-X^{*}B-\tilde{B}^{*})^{-1}$ is a strictly negative matrix. Since $\mathbf{0}\ll x^{*}\ll\mathbf{1}$, it follows that all $x_{i}^{*}$’s strictly decrease when any $\delta_{i}$ increases or any $\beta_{ij}$ decreases. We have thus proved the following result.

Proposition 4

Consider the single-virus model (16). Suppose that $\delta_{i}>0$ for all $i\in[n]$, and that the matrix $B$ is nonnegative and irreducible. If $s(-D+B)>0$, then each entry of the epidemic state $x^{*}$ is a strictly decreasing function of $\delta_{i}$, $i\in[n]$, and a strictly increasing function of $\beta_{ij}$, $i,j\in[n]$.

Similarly, we have the following result for the bi-virus model (10).

Theorem 8

Suppose that $\delta^{1}_{i}>0$, $\delta^{2}_{i}\geq 0$ for all $i\in[n]$, and that matrices $B^{1}$ and $B^{2}$ are nonnegative and irreducible. If $s(-D^{1}+B^{1})>0$ and $s(-D^{2}+B^{2})\leq 0$, then each entry of the epidemic state $\tilde{x}^{1}$ is a strictly decreasing function of $\delta^{1}_{i}$, $i\in[n]$, and a strictly increasing function of $\beta^{1}_{ij}$, $i,j\in[n]$.

Theorem 8 characterizes the effects of adjusting each individual’s healing rate and infection rates, which can be regarded as the simplest local control technique. Such adjustments can be achieved, for example, by taking medicine (i.e., increasing the healing rate) or reducing contact with neighbors (i.e., decreasing the infection rates). The result of Theorem 8 shows that any individual’s local adjustment can attenuate every individual’s epidemic state and thus the pervasion of the dominant virus.

It can be seen that if the individuals have sufficiently large healing rates and small infection rates, both $s(-D^{1}+B^{1})$ and $s(-D^{2}+B^{2})$ will be less than zero and thus both viruses will be eradicated. But to decide whether the healing rates are large enough or not, as well as whether the infection rates are small enough or not, requires centralized computation. With this in mind, we are interested in exploring distributed control techniques to eliminate (not only attenuate) epidemic spreading. It turns out that this is a challenging problem as discussed in the next section.

V Distributed Feedback Control

In this section, we regard each healing rate as a local control input of each group (or agent). We begin with the single-virus model (16).

Suppose that the matrix $B$ is fixed. Let $\delta_{i}=\sum_{j=1}^{n}\beta_{ij}$. Then, the row sums of $D^{-1}B$ all equal $1$. By Lemma 3, $\rho(D^{-1}B)=1$, which is equivalent to $s(-D+B)=0$ because of Proposition 1. Thus, by Proposition 2, the healthy state $\mathbf{0}$ is asymptotically stable in this case. This observation implies that in the case when local control inputs $\delta_{i}$’s are constant, there always exist sufficiently large $\delta_{i}$’s which can stabilize the heathy state.

In the following, we will consider local control inputs of the form

[TABLE]

where $k_{i}$ is a feedback gain. Designing the controller as a (linear) function of the infection proportion $x_{i}(t)$ is an intuitive approach since if the virus is eradicated, no control should be necessary. In implementation, this can be regarded as a treatment plan for individuals via administration of antidote or alternate treatment techniques. By (15), the system reduces to

[TABLE]

The resulting $n$ state equations can be combined to yield

[TABLE]

where $K={\rm diag}([k_{1},\dots,k_{n}])$. Similar to the original system (16), we impose the following restrictions on the parameters of the new system (25).

Assumption 5

For all $i\in[n]$, we have $k_{i}>0$ and the matrix $B$ is nonnegative and irreducible.

Using the same arguments as in the proof of Lemma 8, it is straightforward to verify that the set $[0,1]^{n}$ is still positively invariant for the new system (25). Since both $K$ and $X(t)$ are diagonal matrices, they commute. Then, from (25),

[TABLE]

Thus, the system (25) has the same form as the original system (16), with $D$ and $B$ being replaced by $K$ and $K+B$, respectively.

Note that $K^{-1}(K+B)=I+K^{-1}B$. Since by Assumption 5, $K$ is a positive diagonal matrix and $B$ is an irreducible nonnegative matrix, $K^{-1}$ is a positive diagonal matrix and, thus, $K^{-1}B$ is an irreducible nonnegative matrix. By Lemma 3, $\rho(K^{-1}B)>0$ and, thus, $\rho(I+K^{-1}B)>1$. This observation implies, by Proposition 3, that the new system (25) has a unique nonzero equilibrium $x^{*}$ which satisfies $\mathbf{0}\ll x^{*}\ll\mathbf{1}$ and is asymptotically stable with domain of attraction $[0,1]^{n}\setminus\{\mathbf{0}\}$. We are thus led to the following result.

Proposition 5

Let Assumption 5 hold, and let $x(0)>\mathbf{0}$. Then, for any local control inputs of the form (24), the healthy state $\mathbf{0}$ is not a reachable state of the system (25).

Note that instability of the healthy state can also be shown using Jacobian linearization at $\mathbf{0}$. However the above shows that the origin is not only an unstable equilibrium but actually is a repeller, that is, a perturbation in any direction will drive the state away from $\mathbf{0}$ toward equilibrium $x^{*}\gg\mathbf{0}$.

Now we turn to the bi-virus model (10). We consider local control inputs of the form

[TABLE]

where $k^{1}_{i}$ and $k^{2}_{i}$ are feedback gains. By (9), the system reduces to

[TABLE]

The above equations can be written in matrix form:

[TABLE]

where $K^{1}$ and $K^{2}$ are $n\times n$ diagonal matrices with the $i$th diagonal entries equal to $k^{1}_{i}$ and $k^{2}_{i}$, respectively. Similar to the original system (10), we impose the following restrictions on the parameters of the new system (27).

Assumption 6

For all $i\in[n]$, we have $k^{1}_{i},k^{2}_{i}>0$ and the matrices $B^{1}$ and $B^{2}$ are nonnegative and irreducible.

From the preceding discussion, we have the following.

Theorem 9

Let Assumptions 1 and 6 hold. Then, for any local control inputs of the form (26), the healthy state $(\mathbf{0},\mathbf{0})$ is an unstable equilibrium of the system (27).

Theorem 9 implies that the distributed feedback controller (26) can never stabilize the healthy state. This impossibility result is interesting and surprising because proportional controllers, while suboptimal, function fairly well in many applications. The result, as well as Proposition 5, thus partially explains why distributed control of epidemic networks is a challenging open problem.

VI Simulations

In this section we first compare the model in (10) to the full probabilistic $3^{n}$ state model in (11)-(14). We then illustrate some of the results from Section III via simulation. Finally we explore some interesting behavior via simulations.

VI-A Approximation Accuracy

We evaluate the effectiveness of (10) (with $\beta_{ii}=0\ \forall i$) as an approximation of the $3^{n}$ state model in (11)-(14) for line graphs, star (hub–spoke) graphs, and complete graphs; see Figure 3 for examples of each graph structure. All adjacency matrices for these graphs are symmetric and binary-valued, and both viruses spread on the same graph. In the star graph, the central node is the first agent. Each simulation was run for 10,000 time steps (final time $T=10,000$), with three initial conditions: 1) the first node infected with virus 1 and the second node infected with virus 2, $x^{1_{1}}(0)=[1\ 0\cdots\ 0]^{\top},x^{2_{1}}(0)=[0\ 1\ 0\ \cdots\ 0]^{\top}$, 2) the first two nodes infected with virus 1 and the second two nodes infected with virus 2, $x^{1_{2}}(0)=[1\ 1\ 0\cdots\ 0]^{\top},x^{2_{2}}(0)=[0\ 0\ 1\ 1\ 0\ \cdots\ 0]^{\top}$, and 3) the first node infected with virus 1 and the rest of the nodes infected with virus 2, $x^{1_{3}}(0)=[1\ 0\cdots\ 0]^{\top},x^{2_{3}}(0)=[0\ 1\ \cdots\ 1]$. We explore identical homogeneous viruses, $(\beta,\delta)=(\beta^{1},\delta^{1})=(\beta^{2},\delta^{2})$ in these tests. The $(\beta,\delta)$ pairs used are $[(0.1,1),(0.215,1),(0.464,1),(0.5,0.5),(1,0.464),$ $(1,0.215),(1,0.1)]$, and the number of agents, $n=4,6,8$. We limited simulations to these $n$ values since mean field approximations are typically worse for small values of $n$ and there is a computational limitation due to the size of the $3^{n}$ state model.

The results are given in Figures 4-6 in terms of the 2-norm of the difference between the states of (10) at the final time ($[x^{1}(T);x^{2}(T)]$), and the means of the two viruses in the $3^{n}$ state Markov model at the final time ($[v^{1}(T);v^{2}(T)]$ as defined by (14)).

The accuracy of the approximation appears to be very similar to the single virus case [14, 22]. Since (10) is an upper bounding approximation, the results show that the two models converge to the healthy state for the smaller values of $\frac{\beta}{\delta}$, resulting in small errors between the two models. For many of the larger values of $\frac{\beta}{\delta}$, (10) again performs quite well since it is at an epidemic state and the $3^{n}$ state model does not appear to reach the healthy state in the finite time considered in the simulations ($T=10,000$). Therefore for certain values of $\frac{\beta}{\delta}$ and certain time scales, (10) is a sufficient approximation of the $3^{n}$ state Markov model. For values of $\frac{\beta}{\delta}$ that are near one, the models are quite different, similar to the single virus case. The $3^{n}$ state model appears, in most cases, to be at or close to the healthy state while (10) is at an epidemic state, resulting in large errors.

VI-B Illustrative Examples

This section contains several illustrative and insightful simulations. Given the nature of the network-dependent competing viruses and the behavior of the parallel equilibria ($\tilde{x}^{1}=\alpha\tilde{x}^{2}$) shown by Theorems 6 and 7, we employ a unique coloring scheme. In each figure (Figures 4-6), we plot the initial condition on the left and the final state on the right, and we include a link to a video of the full simulation in the caption. The trajectories of the two viruses for each node are then depicted by the color of the corresponding node. Virus 1 is depicted by the color red ($r=[1\ 0\ 0]$) and virus 2 is depicted by the color blue ($b=[0\ 0\ 1]$). For each $i\in[n]$, the color at time $t$ for group (agent) $i$ is

[TABLE]

When $x^{1}_{i}(t)+x^{2}_{i}(t)=0$, the color goes to white, indicating a completely healthy, susceptible group $i$. These colors are used to facilitate the depiction of the parallel equilibria, which will be shown by all nodes converging to the same color. For all $i\in[n]$, the diameter of the node representing group (agent) $i$ is given by

[TABLE]

with $d_{0}$ being the default/smallest diameter and $r_{0}$ being the scaling factor depending on the total sickness of group (agent) $i$. Therefore, the color indicates the proportion of each virus the group (agent) has and the diameter indicates the strength of the viruses or how sick the group (agent) is.

For systems with two different graph structures, the graph over which virus 1 spreads is depicted by gray edges and the graph over which virus 2 spreads is depicted by green edges. If both viruses spread over the same graph, the edges are gray.

First we illustrate Theorem 1, with $s(-D^{1}+B^{1})=-0.1191$ and $s(-D^{2}+B^{2})=-0.0316$. See Figure 7 for the initial and final states. Consistent with the result of the theorem, both viruses are eradicated.

We illustrate Theorem 2, with $s(-D^{1}+B^{1})=0.4145$ and $s(-D^{2}+B^{2})=-0.0802$, with Figure 9 depicting the initial and final states. Consistent with the result of the theorem, one virus reaches an epidemic state, while the other is eradicated.

Theorem 6 is illustrated in Figure 11, with $\frac{\delta^{1}}{\beta^{1}}=\frac{0.4242}{1.9090}=\frac{\delta^{2}}{\beta^{2}}=\frac{0.2121}{0.9545}$ and $s(-\delta^{1}I+\beta^{1})=4.2654=2s(-\delta^{2}I+\beta^{2}A)$. Consistent with the result of the theorem, both viruses reach an epidemic state where one equilibrium is a scaling of the other, depicted by all the colors being the same, by (28).

A question of interest is, when do certain systems converge to different, or the same, parallel equilibrium? In the following simulation we start the system from Figure 11 with three different initial conditions and they all converge to the same equilibrium. This is depicted in Figure 12. However, consistent with Theorem 6, the system has many equilibria which, via simulation, appear to be initial condition dependent.

Finally, there was earlier discussion about coexisting equilibria and it was shown, for particular cases, that convergence to such equilibria can be proven. However simulations show that generic systems can also converge to epidemic equilibria. Consider the system in Figure 14, where $s(-D^{1}+B^{1})=0.2276>0$, $s(-D^{2}+B^{2})=0.3117>0$, $D^{1}\neq D^{2}$, $B^{1}\neq B^{2}$, and the viruses spread over different graphs. Both viruses clearly reach an epidemic equilibrium; see Figure 14(b). Simulations show for the same system, that different initial conditions lead to an epidemic state for virus 1 and a healthy state for virus 2. Therefore, in the generic case the system is initial condition dependent.

VII Conclusion

In this paper, we have thoroughly analyzed the equilibria of a continuous-time bi-virus model and their stability, and in so doing, as a by-product we have improved on the results for the single-virus model. We have provided necessary and sufficient conditions for convergence to the healthy state of (10). We have also provided results on the epidemic states of (10), including several sufficient conditions for stability and instability. We have explored two distributed control techniques for the model. First, we have regarded each healing or infection rate as a local variable that an individual can control, and shown by sensitivity analysis that any small local adjustment of healing or infection rates can attenuate the global pervasion of the dominant virus. Second, we have regarded each healing rate as a local control input, and shown that a distributed proportional controller of the form $\delta_{i}(t)=k_{i}x_{i}(t)$, can never drive the virus model to the healthy state.

For future work, we plan to study bi-virus models with time-varying graph structures, similar to [21, 22] for the single-virus case. We also plan to analyze the multi-virus case, i.e., more than two competing viruses, for the healthy and epidemic states. The sensitivity analysis over the regime $\rho(D_{1}^{-1}B_{1})>0$ and $\rho(D_{2}^{-1}B_{2})>0$, where both strains would be fit to survive, constitutes a future research direction. Distributed control of epidemic networks, including the bi-virus model considered here, is of course another direction for future research.

We note that there are two limitations of the bi-virus model considered in the paper. First, it was assumed that individuals are susceptible at birth, even if their parents are infected with a virus. This assumption is invalid in some realistic situations as individuals can be infected with a virus at birth if their parents are virus carriers. Second, the infection rate within the $i$th group ($\beta_{ii}$) is probably larger than the infection rate between two different groups ($\beta_{ij}\ i\neq j$) in some realistic situations. Modeling and investigating the dynamics of bi-virus spreading in such scenarios are left as an important and interesting area for future research.

VIII Appendix

Proof of Lemma 1: Suppose, to the contrary, that for all $i\in[n]$ such that $x_{i}=0$, we have $y_{i}=0$. Since $x$ has at least one zero entry and $x>\mathbf{0}$, there exists a proper nonempty subset ${\cal E}\subset[n]$ such that for $i\in{\cal E}$, $x_{i}=0$, and for $i\in[n]\setminus{\cal E}$, $x_{i}>0$. By our assumption, it follows that for any $i\in{\cal E}$, $y_{i}=0$. Without loss of generality, let ${\cal E}=\{1,2,\ldots,m\}$ for some $m$, where $1\leq m<n$. Then,

[TABLE]

with $A,D\geq 0$, $B,C>0$, $\bar{x}\gg\mathbf{0}$, and $\bar{y}\geq\mathbf{0}$. This implies the $B=0$, which is a contradiction since $M$ is an irreducible matrix. Therefore, there exists an index $i\in[n]$ such that $x_{i}=0$ and $y_{i}>0$, and thus $x$ and $y$ cannot have the same sign pattern.

Proof of Proposition 1: Suppose that $\Lambda$ is a negative diagonal matrix in $\mathbb{R}^{n\times n}$ and $N$ is an irreducible nonnegative matrix in $\mathbb{R}^{n\times n}$. Let $M=\Lambda+N$. By Theorem 3.29 in [43], $s(M)<0$ if and only if $\rho(-\Lambda^{-1}N)<1$. To prove the proposition, it suffices to show that $s(M)=0$ if and only if $\rho(-\Lambda^{-1}N)=1$.

First suppose that $s(M)=0$. Set $\Lambda_{\varepsilon}=\Lambda-\varepsilon I$ with $\varepsilon>0$. Let $M_{\varepsilon}=\Lambda_{\varepsilon}+N=\Lambda-\varepsilon I+N$. Then, $\lim_{\varepsilon\rightarrow 0^{+}}\rho(-\Lambda_{\varepsilon}^{-1}N)=\rho(-\Lambda^{-1}N)$. Since $\varepsilon>0$, $s(M_{\varepsilon})<0$. Then, $\rho(-\Lambda_{\varepsilon}^{-1}N)<1$ and, therefore, $\lim_{\varepsilon\rightarrow 0^{+}}\rho(-\Lambda_{\varepsilon}^{-1}N)\leq 1$. Thus, $\rho(-\Lambda^{-1}N)\leq 1$. To prove that $\rho(-\Lambda^{-1}N)=1$, suppose that, to the contrary, $\rho(-\Lambda^{-1}N)<1$. Then, $s(M)<0$, which is a contradiction. Therefore, $\rho(-\Lambda^{-1}N)=1$.

Now suppose that $\rho(-\Lambda^{-1}N)=1$. Again set $\Lambda_{\varepsilon}=\Lambda-\varepsilon I$ with $\varepsilon>0$ and $M_{\varepsilon}=\Lambda_{\varepsilon}+N$. Then, $\lim_{\varepsilon\rightarrow 0^{+}}s(M_{\varepsilon})=s(M)$. Since $\varepsilon>0$, $-\Lambda_{\varepsilon}^{-1}N$ is a nonnegative matrix. Since $N$ is irreducible and nonnegative, so is $-\Lambda_{\varepsilon}^{-1}N$. Note that the $i$th diagonal entry of $-\Lambda_{\varepsilon}$ is strictly larger than the $i$th diagonal entry of $-\Lambda$ since $\varepsilon>0$. Thus, $-\Lambda^{-1}N>-\Lambda_{\varepsilon}^{-1}N$. By Lemma 3, $\rho(-\Lambda_{\varepsilon}^{-1}N)<1$. Then, $s(M_{\varepsilon})<0$ and, thus, $\lim_{\varepsilon\rightarrow 0^{+}}s(M_{\varepsilon})\leq 0$. Thus, $s(M)\leq 0$. To prove that $s(M)=0$, suppose that, to the contrary, $s(M)<0$. Then, $\rho(-\Lambda^{-1}N)<1$, which is a contradiction. Therefore, $s(M)=0$.

Proof of Lemma 8: Suppose that at some time $\tau$, $x^{1}_{i}(\tau),x^{2}_{i}(\tau),x^{1}_{i}(\tau)+x^{2}_{i}(\tau)\in[0,1]$ for all $i\in[n]$. Consider an index $i\in[n]$. If $x^{1}_{i}(\tau)=0$, then from (9) and Assumption 2, $\dot{x}^{1}_{i}(\tau)\geq 0$. The same holds for $x^{2}_{i}(\tau)$ and $x^{1}_{i}(\tau)+x^{2}_{i}(\tau)$. If $x^{1}_{i}(\tau)=1$, then from (9) and Assumption 2, $\dot{x}^{1}_{i}(\tau)\leq 0$. The same holds for $x^{2}_{i}(\tau)$ and $x^{1}_{i}(\tau)+x^{2}_{i}(\tau)$. It follows that $x^{1}_{i}(t),x^{2}_{i}(t),x^{1}_{i}(t)+x^{2}_{i}(t)$ will be in $[0,1]$ for all times $t\geq\tau$. Since the above arguments hold for all $i\in[n]$, $x^{1}_{i}(t),x^{2}_{i}(t),x^{1}_{i}(t)+x^{2}_{i}(t)$ will be in $[0,1]$ for all $i\in[n]$ and $t\geq\tau$. Since by Assumption 1 $x^{1}_{i}(0),x^{2}_{i}(0),x^{1}_{i}(0)+x^{2}_{i}(0)\in[0,1]$ for all $i\in[n]$, it follows that $x^{1}_{i}(t),x^{2}_{i}(t),x^{1}_{i}(t)+x^{2}_{i}(t)\in[0,1]$ for all $i\in[n]$ and $t\geq 0$.

Proof of Proposition 2 (case $s(-D+B)\leq 0$): We first consider the case when $s(-D+B)<0$. Since $(-D+B)$ is an irreducible Metzler matrix, by Lemma 6, there exists a positive diagonal matrix $P$ such that $(-D+B)^{\prime}P+P(-D+B)$ is negative definite. Consider the Lyapunov function candidate $V(x(t))=x(t)^{\prime}Px(t)$. From (16), when $x(t)\neq\mathbf{0}$,

[TABLE]

Thus, in this case, $\dot{V}(x(t))<0$ if $x(t)\neq\mathbf{0}$. By Lemma 8 and Corollary 1, $x=\mathbf{0}$ is asymptotically stable with domain of attraction $[0,1]^{n}$.

Next we consider the case when $s(-D+B)=0$. Since $(-D+B)$ is an irreducible Metzler matrix, by Lemma 7, there exists a positive diagonal matrix $P$ such that $(-D+B)^{\prime}P+P(-D+B)$ is negative semi-definite. Consider the Lyapunov function candidate $V(x(t))=x(t)^{\prime}Px(t)$. From (16), we have

[TABLE]

We claim that $\dot{V}(x(t))<0$ if $x(t)\neq\mathbf{0}$. To establish this claim, we first consider the case when $x(t)\gg\mathbf{0}$. Since $B$ is nonnegative and irreducible, $Bx(t)\gg\mathbf{0}$. Since $P$ is a positive diagonal matrix, it follows that $x(t)^{\prime}PX(t)Bx(t)>0$, so $\dot{V}(x(t))<0$. Next we consider the case when $x(t)>\mathbf{0}$ and $x(t)$ has at least one zero entry. Since $(-D+B)$ is an irreducible Metzler matrix and $P$ is a positive diagonal matrix, $(-D+B)^{\prime}P+P(-D+B)$ is a symmetric irreducible Metzler matrix. Since $(-D+B)^{\prime}P+P(-D+B)$ is negative semi-definite, it follows that $s((-D+B)^{\prime}P+P(-D+B))=0$. By Lemma 3, [math] is a simple eigenvalue of $(-D+B)^{\prime}P+P(-D+B)$ and it has a unique (up to scalar multiple) strictly positive eigenvector corresponding to the eigenvalue [math]. That is, there exists $\tilde{x}\gg 0$ such that $A\tilde{x}=0\tilde{x}$. Thus, $x(t)^{\prime}\left((-D+B)^{\prime}P+P(-D+B)\right)x(t)<0$ when $x(t)>\mathbf{0}$ and $x(t)$ has at least one zero entry. Therefore, $\dot{V}(x(t))<0$ if $x(t)\neq\mathbf{0}$. By Lemma 8 and Corollary 1, $x=\mathbf{0}$ is asymptotically stable with domain of attraction $[0,1]^{n}$.

To proceed, we need the following lemma.

Lemma 12

Suppose that $\delta_{i}\geq 0$ for all $i\in[n]$ and that matrix $B$ is nonnegative and irreducible. If $x^{*}$ is a nonzero equilibrium of system (16), then $x^{*}\gg\mathbf{0}$.

Proof of Lemma 12: Suppose that $x^{*}$ is a nonzero equilibrium of (16). By Lemma 8, it must be true that $x^{*}\geq\mathbf{0}$. To prove $x^{*}\gg\mathbf{0}$, suppose that, to the contrary, $x^{*}$ has at least one zero entry. Without loss of generality, set $x^{*}_{1}=0$. Since $x^{*}$ is an equilibrium of (16), from (15),

[TABLE]

It then follows that for any $j\in[n]$ such that $\beta_{1j}>0$, $x^{*}_{j}=0$. By repeating this argument, since $B$ is irreducible, we have $x^{*}_{i}=0$ for all $i\in[n]$. This contradicts the assumption that $x^{*}\geq 0$. Thus, $x^{*}\gg\mathbf{0}$.

Proof of Proposition 2 (case $s(-D+B)>0$): It will suffice to show that if $s(-D+B)>0$, then there exists a unique strictly positive equilibrium. We first show that there exists an $x^{*}\gg\mathbf{0}$ which is an equilibrium of (16).

Let $c>0$ be any positive constant such that

[TABLE]

Such a constant $c$ always exists since $s(-D+B)>0$. Set $\bar{D}=D+cI$. Since $\delta_{i}\geq 0$ for all $i\in[n]$, $D$ is a nonnegative diagonal matrix. Thus, $\bar{D}$ is nonsingular and $\bar{D}^{-1}$ is also a positive diagonal matrix.

Consider the above equation and define a continuous map $f:(0,1]^{n}\rightarrow[0,1]^{n}$ given by

[TABLE]

Since the domain of $f$ is $(0,1]^{n}$, $x$ as the argument of $f$ satisfies $x\gg\mathbf{0}$. Therefore ${\rm diag}(\bar{D}^{-1}Bx)$ is a positive diagonal matrix, and we have that $(I-c\bar{D}^{-1}+{\rm diag}(\bar{D}^{-1}Bx))$ is invertible. Therefore $f$ is well-defined. Note that the $i$th entry of $f(x)$, denoted by $f_{i}(x)$, is given by

[TABLE]

Since $\bar{D}^{-1}$ and $B$ are both nonnegative, for any $y\geq z$ in $(0,1]^{n}$, $f_{i}(y)\geq f_{i}(z)$, so $f(y)\geq f(z)$.

Since $\bar{D}^{-1}B$ is an irreducible nonnegative matrix, by Lemma 3, there exists $v\gg\mathbf{0}$ such that

[TABLE]

where $r=\rho(\bar{D}^{-1}B)$. Since $s(-\bar{D}+B)=s(-D+B)-c$, from (30), $s(-\bar{D}+B)>0$. By Proposition 1, it follows that $r>1$. Then, we can always find an $\varepsilon>0$ such that for each $i\in[n]$,

[TABLE]

From this, it follows that $1\leq\frac{r}{1+\varepsilon rv_{i}}$, and thus, $\varepsilon v_{i}\leq\frac{\varepsilon rv_{i}}{1+\varepsilon rv_{i}}$. From (31), we have

[TABLE]

which implies that $\varepsilon v\leq f(\varepsilon v)$. It follows from (32) that $\varepsilon v\ll\mathbf{1}$. Since we have already shown that for any $y\geq z$ in $(0,1]^{n}$, $f(y)\geq f(z)$, $f$ maps the compact, convex set ${\cal C}=\{x\ |\ \varepsilon v\leq x\leq\mathbf{1}\}$ to itself. By Brouwer’s fixed-point theorem, $f$ has a fixed point in ${\cal C}$, which must be strictly positive. Let $x^{*}\gg\mathbf{0}$ denote this fixed point. Then $f(x^{*})=x^{*}$, i.e.,

[TABLE]

Therefore, we have

[TABLE]

Therefore, $Bx^{*}=\bar{D}x^{*}+X^{*}Bx^{*}-cx^{*}=Dx^{*}+X^{*}Bx^{*}$, by definition of $\bar{D}$. Thus, we have $((-D+B)-X^{*}B)x^{*}=\mathbf{0}$, and therefore $x^{*}\gg\mathbf{0}$ is an equilibrium of (16).

It remains to show that the strictly positive equilibrium is unique. Suppose that $x$ and $y$ are both nonzero equilibria of (16), and let $\varepsilon$ from (32) be sufficiently small such that $x,y\in{\cal C}$. From Lemma 12, it follows that $x,y\gg\mathbf{0}$. Set

[TABLE]

Then, $x\leq\kappa y$, and there exists $j\in[n]$ for which $x_{j}=\kappa y_{j}$. We claim that $\kappa\leq 1$. To establish this claim, suppose that, to the contrary, $\kappa>1$. Since $x$ is a fixed point of $f$ and for any $u\geq v$ in $(0,1]^{n}$, $f_{j}(u)\geq f_{j}(v)$ for all $j\in[n]$, it follows that

[TABLE]

From the assumption that $\kappa>1$, we have

[TABLE]

Since $y$ is a fixed point of $f$,

[TABLE]

Then, it follows that

[TABLE]

which is a contradiction. Therefore, $\kappa\leq 1$, which implies that $x\leq y$. Using the same arguments and by exchanging the roles of $x$ and $y$, it also can be shown that $y\leq x$. Thus, $x=y$, which establishes the uniqueness of the strictly positive equilibrium. This completes the proof.

Proof of Lemma 9: If $x(0)\gg\mathbf{0}$, then the lemma is true with $\tau=0$. Suppose that $x(0)>\mathbf{0}$ with $x_{i}(0)=0$ for at least one $i\in[n]$. Let ${\cal F}(t)$ be the set of all $i\in[n]$ such that $x_{i}(t)=0$. In other words, $x_{i}(t)=0$ for all $i\in{\cal F}(t)$ and $x_{i}(t)>0$ for all $i\in[n]\setminus{\cal F}(t)$. Clearly, the set ${\cal F}(0)$ is nonempty. Moreover, since the matrix $B$ is irreducible, there exists $j\in{\cal F}(0)$ such that $x_{j}(0)=0$, $x_{k}(0)>0$, and $\beta_{jk}>0$. From (15), it follows that $\dot{x}_{j}(0)>0$. Thus, there must exist $\tau_{1}>0$ such that $x_{j}(\tau_{1})>0$ and $x_{i}(\tau_{1})>0$ for all $i\in[n]\setminus{\cal F}(0)$. This implies that ${\cal F}(\tau_{1})$ is a proper subset of ${\cal F}(0)$. Note that ${\cal F}(0)$ is a finite set. By repeating the above arguments, we conclude that there exists $\tau>0$ such that ${\cal F}(\tau)$ is the empty set, which implies that $x_{i}(\tau)>0$ for all $i\in[n]$.

Proof of Proposition 3: Let $y_{i}(t)=x_{i}(t)-x^{*}_{i}$ for all $i\in[n]$, i.e., $y(t)=x(t)-x^{*}$. Let $Y(t)={\rm diag}(y(t))$ and $X^{*}={\rm diag}(x^{*})$. Note that $\left(-D+B-X^{*}B\right)x^{*}=\mathbf{0}$. Then,

[TABLE]

Thus, for all $i\in[n]$,

[TABLE]

Consider the Lyapunov function candidate

[TABLE]

which is well-defined since $x^{*}_{k}>0$ for all $k\in[n]$ by Proposition 2. Note that $V(y(t))\geq 0$, with equality if and only if $y(t)=\mathbf{0}$ (or equivalently, $x(t)=x^{*}$). Moreover, for all $i\in[n]$ and $t\geq 0$, there holds $|y_{i}(t)|\leq V(y(t))x^{*}_{i}$. For any time $t\geq 0$, without loss of generality, let $m\in[n]$ such that

[TABLE]

Then, when $|y_{m}(t)|>0$ (or equivalently, $V(y(t))>0$),

$\displaystyle\dot{V}(y(t))=\frac{1}{x^{*}_{m}}\cdot\frac{d|y_{m}(t)|}{dt}$

[TABLE]

From (VIII) and the definition of $y_{i}(t)$, it is straightforward to verify that $\dot{V}(y(t))=0$ in the case when $x(t)=x^{*}$ (i.e., $y(t)=\mathbf{0}$). Next we consider the case when $x(t)\gg\mathbf{0}$. From (VIII), since the matrix $B$ is irreducible and $x^{*}\gg\mathbf{0}$ by Proposition 2, it can be seen that $\dot{V}(y(t))<0$ if $x(t)\neq x^{*}$.

Recall that by Lemma 9, as long as $x(0)>\mathbf{0}$, there always exists a finite time $t_{0}$ at which $x(t_{0})\gg\mathbf{0}$. From Lemma 8 (with $x^{2}=\mathbf{0}$, this reduces to the single virus model case) and Corollary 1, to prove the proposition, it remains to show that the system (16) is positively invariant on a subset of $(0,1]^{n}$. Without loss of generality, suppose that $x(0)\gg\mathbf{0}$. Let $\varepsilon$ be a nonnegative real number such that

[TABLE]

Consider the set ${\cal B}=\left\{x\ |\ V(x(t))\leq\varepsilon\right\}\cap\ (0,1]^{n}$. Since $\dot{V}(x(t))\leq 0$, the system (16) is positively invariant on $\{x\ |\ V(x(t))\leq\varepsilon\}$. Suppose that there exists some finite time $t_{1}$ at which $x(t_{1})>\mathbf{0}$ and $x(t_{1})\notin(0,1]^{n}$. From Lemma 9, there must exist another finite time $t_{2}>t_{1}$ such that $x(t_{2})\gg\mathbf{0}$. Since we have shown that $\dot{V}(x(t))<0$ if $x(t)\gg\mathbf{0}$ and $x(t)\neq x^{*}$, it follows that $x(t)\in(0,1]^{n}$ for all $t\geq t_{2}$. Note that

[TABLE]

It then follows that for any $z\in(0,1]^{n}$ which is close to $\mathbf{0}$, there holds $V(z)<V(\mathbf{0})$. Since $\dot{V}(x(t))<0$ if $x(t)\gg\mathbf{0}$ and $x(t)\neq x^{*}$, there cannot exist a trajectory from $x(0)$ to $\mathbf{0}$. Thus, the system (16) is positively invariant on $(0,1]^{n}$ for $t\geq t_{2}$, which implies that the system (16) is positively invariant on ${\cal B}$ for $t\geq t_{2}$. Since the above arguments hold for any $x(0)\gg\mathbf{0}$, such a set ${\cal B}$ always exists, and thus after some finite time, the system (16) is positively invariant on ${\cal B}\subset(0,1]^{n}$.

Proof of Lemma 10: Let $(\tilde{x}^{1},\tilde{x}^{2})$ be an equilibrium of (10) and suppose, to the contrary, that there exists some $i\in[n]$ such that $\tilde{x}^{1}_{i}+\tilde{x}^{2}_{i}=1$. Let $(\tilde{x}^{1}(0),\tilde{x}^{2}(0))=(\tilde{x}^{1},\tilde{x}^{2})$. Then, from (9), $\dot{\tilde{x}}^{1}_{i}(0)<0$ and/or $\dot{\tilde{x}}^{2}_{i}(0)<0$, which contradicts the hypothesis that $(\tilde{x}^{1}_{i},\tilde{x}^{2}_{i})$ is an equilibrium. Therefore, $\tilde{x}^{1}(\tau)+\tilde{x}^{2}(\tau)\ll\mathbf{1}$.