Multi-Competitive Viruses over Static and Time--Varying Networks

TL;DR

This paper studies the spread and control of multiple competing viruses over both static and dynamic networks, extending previous work limited to two viruses and static graphs, and introduces stability analysis and an antidote control method.

Contribution

It generalizes the analysis of competing viruses to multiple viruses and dynamic networks, providing new stability results and a control strategy.

Findings

Stability conditions for multi-virus systems are derived.

An antidote control technique is proposed and analyzed.

The model extends previous two-virus static graph studies.

Abstract

Epidemic processes are used commonly for modeling and analysis of biological networks, computer networks, and human contact networks. The idea of competing viruses has been explored recently, motivated by the spread of different ideas along different social networks. Previous studies of competitive viruses have focused only on two viruses and on static graph structures. In this paper, we consider multiple competing viruses over static and dynamic graph structures, and investigate the eradication and propagation of diseases in these systems. Stability analysis for the class of models we consider is performed and an antidote control technique is proposed.

Equations

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

$$\dot{p}(t)=(B-P(t)B-D)p(t),$$

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

$${\dot{p}}_i^k(t)=(1-p_i^1(t)-\cdots -p_i^m(t))\sum _{j=1}^n{\beta }_{ij}^k{p}_j^k(t)-{\delta }_i^k{p}_i^k(t),$$

$${\dot{p}}^k(t)=((I-P^1(t)-\cdots -P^m(t))B^k-D)p^k(t),$$

$$\mathcal{D}=\{(p^1,\dots ,p^m)|p^k\ge 0,k\in [m],\sum _{k=1}^m{p}^k\le 1\}$$

$${\dot{p}}_i^1(t)\le -{\delta }_i^1{p}_i^k(t)+(1-p_i^k(t))\sum _{j=1}^n{\beta }_{ij}^k{p}_j^1(t),$$

$$\begin{array}[]{lcl}\dot{V}(p^{k})&=&(p^{k})^{T}\dot{p}^{k}\\ &=&(p^{k})^{T}(B^{k}-\sum_{l=1}^{m}P^{l}B^{k}-D^{k})p^{k}\\ &\leq&(p^{k})^{T}(B^{k}-D^{k})p^{k}\\ &\leq&\lambda_{1}(B^{k}-D^{k})\|p^{k}\|^{2}<0.\end{array}$$

$${\dot{p}}_i^k(t)=(1-p_i^1(t)-\cdots -p_i^m(t))\sum _{j=1}^n{\beta }_{ij}^k(t)p_j^k(t)-{\delta }_i^k{p}_i^k(t),$$

$$\begin{array}[]{lcl}\dot{V}(p^{k})&=&(p^{k})^{T}\dot{p}^{k}\\ &=&(p^{k})^{T}(B^{k}(t)-\sum_{l=1}^{m}P^{l}B^{k}(t)-D^{k})p^{k}\\ &\leq&(p^{k})^{T}(B^{k}(t)-D^{k})p^{k}\\ &\leq&\lambda_{1}(B^{k}(t)-D^{k})\|p^{k}\|^{2}\\ &\leq&(\sup_{t\geq 0}\lambda_{1}(B^{k}(t)-D^{k}))\|p^{k}\|^{2}<0.\end{array}$$

$$\Vert B^k(t)-D^k\Vert \le c^k(t)e^{-{\lambda }^k(t)t}\forall t\ge 0.$$

$${\gamma }_1^k$$

$${\gamma }_1^k\ge \Vert {\int }_0^{\infty }{e}^{(B^k(t)-D^k{)}^T\tau }{e}^{(B^k(t)-D^k)\tau }d\tau \Vert .$$

$$\begin{array}[]{lcl}\dot{p}^{k}&=&(B^{k}(t)-\sum_{l=1}^{m}P^{l}B^{k}-D^{k})p^{k}\\ &\leq&(B^{k}(t)-D^{k})p^{k}.\end{array}$$

$${\dot{p}}^k=(B^k(t)+{\Delta }^k(t)-P(t)(B^k(t)+{\Delta }^k(t))-D^k)p^k.$$

$$\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_{t_0}^{t_0+T}\Vert B^k(s)-D^k\Vert ds\le a<\infty ,$$

$$\Vert B^k(t+h)-B^k(t)\Vert \le \nu h^{\gamma },$$

$$\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_{t_0}^{t_0+T}{s}_1(B^k(s)-D^k)ds\le \bar{\alpha },$$

$$\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_{t_0}^{t_0+T}\Vert {\Delta }^k(s)\Vert ds\le \eta <\infty ,$$

$$|{\Delta }_{ij}^k(t)|\le {\beta }_{ij}^k(t).$$

$${\dot{p}}^k$$

$$\le (B^k(t)+{\Delta }^k(t))-D^k)p^k.$$

$${\dot{p}}^i(t)=(B^i-D^i-P^i(t)B^i)p^i(t)-\sum _{k\ne i}{P}^k(t)B^k{p}^k(t).$$

$${\dot{p}}^i(t)=(B^i-D^i-P^i(t)B^i)p^i(t),$$

$$\begin{array}{r}(I-{\tilde{P}}^1-\cdots -{\tilde{P}}^m)A{\tilde{p}}^k=\frac{{\delta }^k}{{\beta }^k}{\tilde{p}}^k,\end{array}$$

$${\dot{p}}_i^k(t)=(1-p_i^1(t)-\cdots -p_i^m(t))\sum _{j=1}^n{\beta }_{ij}^k{p}_j^k(t)-({\delta }_i^k+u_i^k(t))p_i^k(t).$$

$$\underset{u_i^k(t)}{\mathrm{minimize}}$$

$$\sum _{i=1}^n{u}_i^k(t)\le c^k,t\ge 0,$$

$$U^k(t)=\text{diag}(u_1^k(t),\dots ,u_n^k(t)),$$

$$u_i^k(t)\ge 0,i=1,\dots ,n,t\ge 0.$$

$$\underset{{\eta }^k,u_i^k(t)}{\mathrm{minimize}}$$

$${\eta }^k\ge \sum _{j=1}^n{\beta }_{ij}^k(t)-({\delta }_i^k+u_i^k(t)),$$

$$\sum _{i=1}^n{u}_i^k(t)\le c^k,$$

Full text

I Introduction

Spread dynamics have been studied for hundreds of years. Bernoulli developed one of the first known models inspired by the smallpox virus [1]. In this paper we focus exclusively on susceptible-infected-susceptible (SIS) models, which have been developed for both continuous [2, 3, 4, 5] and discrete time domains [5, 6, 7]. SIS models consist of a number of agents that are either infected or healthy (susceptible), which may cycle (aperiodically) between these two states. The infection rate combined with the connectivity of the $i$th agent with infected neighbors $j$ (denoted by $\beta_{ij}$) positively affects the probability of being infected, while the healing rate $\delta_{i}$ negatively affects the infection probability. This is depicted in Figure 1a.

The idea of two competing SIS viruses, namely the bi-virus model, has been recently pursued in [8, 9, 10, 11, 12, 13]. The main motivation for such systems is that of competing ideas spreading on different social networks. However these models can have broader applications to political stances, adaptation of competing products, competing practices in farming, etc. and can be generalized to more than two viruses. Consider, for example, the case of three competing viruses; then each state has four possible states: susceptible, infected with virus $1$, $2$, or $3$. The idea of information diffusion on two layered networks has also been explored for a susceptible-infected-recovered (SIR) model in [14]. In [15], a different multi-virus model is considered.

Further, all previous work on competing viruses has focused on viruses over static graph structures. There are recent results for the single virus model over time–varying networks [16, 17, 18, 19, 20]. Some of the ideas from [19, 20] will be employed in this paper and applied to a more general model.

Various control techniques have been applied to SIS virus systems [21, 22, 23, 13]. These techniques assume the healing rate is a control variable. In [13], it is shown that there exists no distributed linear feedback control that can stabilize the system, and in fact, will destabilize the system. Alternative approaches focus on reducing the maximum eigenvalue of the linearized system using the healing rate and/or the infection rate. In [21, 22], distributed control techniques for setting healing rate and quarantine protocols are proposed and implemented on a severe acute respiratory syndrome (SARS) simulation model. In [23], a bound is provided for the cost of fairness of mitigating the spread of disease, that is, the difference between the optimal solution and the fair or homogeneous solution, for several classes of graphs. In [25], geometric programming ideas are used to control single SIS virus systems and the authors present a polynomial time algorithm illustrated on an air transportation network. In [24], similar ideas to [25] are applied to the bi-virus model.

In this paper we present a generalization of the bi-virus model to an arbitrary number, $m$, of competing viruses. We provide conditions for stability of the disease-free equilibrium (DFE) for static as well as time–varying graph structures. We also provide sufficient conditions for stability of the non-disease free equilibrium (NDFE). We provide two control techniques based on minimizing the maximum eigenvalue of the linearized system, appealing to some of the theorems presented herein. These control techniques, which are different from other approaches in the literature, allow every agent to have a base healing rate and an additive control term.

The paper is organized as follows: we first introduce, in Section II, the SIS model and the competing virus model for $m$ viruses. In Sections III and IV we analyze the model, providing conditions for stability of the DFE and the NDFE, and in Section V we provide an antidote control formulation. In Section VI, we present a set of illuminating simulations of various competing virus models over time–varying networks, and we conclude with some discussion in Section VII.

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^{T}$ to denote the transpose of a vector $x$. 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}$ and, also, by $[A]_{ij}$ when convenient. We use $\mathbf{0}$ and $\mathbf{1}$ to denote the vectors whose entries are all equal to [math] and $1$, respectively, and $I$ to denote the identity matrix, while the dimensions of the vectors and matrices are to be understood from the context. For any vector $x\in{\rm I\!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 in ${\cal A}$ but not in ${\cal B}$.

For any two real vectors $a,b\in{\rm I\!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]$. For a real square matrix $M$, we use $s(M)$ to denote the largest real part among its eigenvalues, i.e., $s(M)=\max\left\{{\rm Re}(\lambda)\ :\ \lambda\in\sigma(M)\right\},$ where ${\rm Re}(\cdot)$ is the real part of the argument and $\sigma(M)$ denotes the spectrum of $M$. For a symmetric matrix $M$, we use $\lambda_{1}(M)$ to denote its largest eigenvalue.

II The Model

The generic SIS model, a generalization of the models introduced in [4], is

[TABLE]

where $p_{i}$ is the probability that agent $i$ is infected, the $\beta_{ij}$’s are (possibly asymmetric) infection rates incorporating the nearest-neighbor graph structure, and $\delta_{i}$ is the healing rate. Neighbor relationships among the $n$ agents are described by a directed graph $\mathbb{G}$ on $n$ vertices with an arc from vertex $j$ to vertex $i$ whenever agent $i$ can be infected by agent $j$. The agents can also be thought of as groups of people and $p_{i}$’s as the percentages of the groups that are infected, and therefore the neighbor graph $\mathbb{G}$ can have self-arcs at all $n$ vertices. Hence, $\beta_{ij}$ equals zero if there is not an edge in $\mathbb{G}$ from node $j$ to node $i$. The model in (1) is more general because the underlying graph $\mathbb{G}$ can be directed and the weights given by $\beta_{ij}$ can be any non-negative number. The representation in (1) can be put into matrix form:

[TABLE]

where $p$ is the vector of the $p_{i}$’s, $B$ is the matrix of the $\beta_{ij}$’s, $P=\text{diag}(p)$, and $D=\text{diag}(\delta_{1},\dots,\delta_{n})$. In the analysis that follows, as stated above, $B$ is not assumed to be symmetric unless explicitly stated so.

This model has been extended to have two viruses, providing a generalization of the model introduced in [10],

[TABLE]

where ${p}^{1}_{i}(t)$ and ${p}^{2}_{i}(t)$ are the probabilities that agent $i$ has virus $1$ and $2$ respectively, and each virus has its own infection rates and healing rates. Each virus spreads over a (possibly different) spanning subgraph of $\mathbb{G}$, where their union is the neighbor graph $\mathbb{G}$. It will be assumed that both of the two subgraphs are strongly connected and, thus, so is $\mathbb{G}$.111 A directed graph is strongly connected if for any two distinct vertices $i$ and $j$, there is a directed path from $i$ to $j$.

We need not restrict ourselves to two viruses, however. A direct generalization leads to the following multi-virus model:

[TABLE]

for all $k\in[m]$. This representation can be written in matrix form as:

[TABLE]

where the matrices are the same as in (2), but now they are dependent on which virus they correspond to. Since the subgraph for each virus $k$ is strongly connected, it follows that $B^{k}$ is irreducible, meaning that it cannot be permuted into block triangular matrix form. The assumption that $B^{k}$ is bounded means $\forall i,j,\ \beta^{k}_{ij}<\infty$.

The set

[TABLE]

is invariant with respect to the system defined by (5). If $p^{k}_{i}$ denotes the probability of agent $i$ being infected by virus $k$ and $1-\sum_{k=1}^{m}p_{i}^{k}$ denotes the probability of agent $i$ being healthy, it is natural to assume that their initial values are in $[0,1]$, since otherwise the values will lack any physical meaning for the epidemic model considered herein. Similarly, if the states were representative of the density of infected members of a sub-population, they would also be bounded between zero and one.

Lemma 1**.**

Suppose that for all $i\in[n],k\in[m]$, we have $\delta^{k}_{i}\geq 0$, and the matrices $B^{k}$ are non-negative. If for all $i\in[n],k\in[m]$, we have $p^{k}_{i}(0),(1-p^{1}_{i}(0)-\cdots-p^{m}_{i}(0))\in[0,1]$, then $p^{k}_{i}(t),p^{1}_{i}(t)+\cdots+p^{m}_{i}(t)\in[0,1]$ for all $i\in[n],k\in[m]$ and $t\geq 0$.

Proof.

Suppose that at some time $\tau$, $p^{1}_{i}(\tau)+\cdots+p^{m}_{i}(\tau)\in[0,1]$ and $p^{k}_{i}(\tau)\in[0,1]$ for all $i\in[n],k\in[m]$. Consider an index $i\in[n]$. If $p^{k}_{i}(\tau)=0$, then from (4) and the assumption that the matrices $B^{k}$ are non-negative, $\dot{p}^{k}_{i}(\tau)\geq 0$. The same holds for $p^{1}_{i}(\tau)+\cdots+p^{m}_{i}(\tau)$. If $p^{k}_{i}(\tau)=1$, then from (4) and the assumption that the matrices $B^{k}$ are non-negative, $\dot{p}^{1}_{i}(\tau)\leq 0$. The same holds for $p^{1}_{i}(\tau)+\cdots+p^{m}_{i}(\tau)$. It follows that $p^{k}_{i}(t),p^{1}_{i}(t)+\cdots+p^{m}_{i}(t)\in[0,1]$ for all $i\in[n],k\in[m]$ and $t\geq\tau$.

Since, by assumption, $p^{k}_{i}(0),(1-p^{1}_{i}(0)-\cdots-p^{m}_{i}(0))\in[0,1]$ for all $i\in[n],k\in[m]$, it follows that $p^{k}_{i}(t),p^{1}_{i}(t)+\cdots+p^{m}_{i}(t)\in[0,1]$ for all $i\in[n],k\in[m]$ and $t\geq 0$. ∎

For the rest of the paper we assume $p^{k}_{i}(0),(1-\sum_{j}^{m}p^{j}_{i}(0))\in[0,1]$ for all $i\in[n],k\in[m]$.

It has been shown that there are disease-free equilibrium and non-disease free equilibria for the single virus system [26, 27, 19, 20], as well as for the two-virus system [13]; the same applies to multi-virus systems as well. However, in this case the scenario becomes slightly more complicated because all viruses can reach the DFE, or a NDFE, or there may be some viruses at a DFE and some at a NDFE. We will explore several conditions for convergence to these different equilibria.

III Stability Analysis of the DFE

First, we explore stability of the DFE for both the static and dynamic graph cases.

III-A Static Graph Structure

We first give conditions under which the DFE is asymptotically stable.

Theorem 1**.**

Suppose that for all $i\in[n],k\in[m]$, we have $\delta^{k}_{i}\geq 0$ and the matrices $B^{k}$ are non-negative and irreducible. If $s(B^{k}-D^{k})\leq 0$ for all $k\in[m]$, then the healthy state is the unique equilibrium of (5), which is asymptotically stable with domain of attraction $\mathcal{D}$, as defined in (6).

Proof.

To prove the theorem, it is sufficient to show that for all $k\in[m]$, $p^{k}(t)$ will asymptotically converge to $\mathbf{0}$ as $t\rightarrow\infty$ for any initial condition.

Since for all $k\in[m]$, $p^{k}_{i}(t)$ is always non-negative by Lemma 1, from (4),

[TABLE]

which implies that the trajectories of $p^{k}_{i}(t)$ are bounded above by a single-virus model. Since the $B^{k}$’s are non-negative and irreducible, by Proposition 3 in [13], $p^{k}_{i}(t)$ will asymptotically converge to $\mathbf{0}$ as $t\rightarrow\infty$ for all $k\in[m]$, and thus the healthy state is the unique equilibrium of (5). ∎

We next state a result on global exponential stability for the case when the underlying subgraphs are undirected and the infection rates are symmetric.

Theorem 2**.**

Suppose $B^{k}$ is symmetric, and the maximum eigenvalue of $B^{k}-D^{k}$ is less than zero, that is $\lambda_{1}(B^{k}-D^{k})<0$. Then the DFE is exponentially stable for virus $k$, with domain of attraction $\mathcal{D}$, in (6).

Proof.

Consider the Lyapunov function $V(p^{k})=\frac{1}{2}(p^{k})^{T}p^{k}$. For $p^{k}\neq 0$,

[TABLE]

The first inequality holds because $(P^{l}B^{k})_{ij}\geq 0,\ \forall l,i,j$ by construction since each $p^{l}_{i}(t)$ is a probability. The second inequality holds by the Rayleigh-Ritz Theorem because $B^{k}-D^{k}$ is symmetric. Therefore, the system converges exponentially fast to the origin by Theorem 8.5 in [28]. ∎

Note that this is a generalization of the result in [19, 20].

We can state that the condition in Theorem 1 is necessary and sufficient for eradication of all viruses.

Theorem 3**.**

Suppose $\delta^{k}_{i}\geq 0$, for all $i,k$, and the matrices $B^{k}$ are non-negative and irreducible for all $k$. The DFE (all $k$ viruses eradicated) is the unique equilibrium of (5) if and only if $s(B^{k}-D^{k})\leq 0$ for all $k\in[m]$.

Proof.

Sufficiency has been shown in Theorem 1. Therefore, to prove the theorem, all that needs to be shown is that if for any $j\in[m]$ $s(B^{j}-D^{j})>0$, the system (5) admits a NDFE.

Without loss of generality, suppose that $s(B^{1}-D^{1})>0$. Set $p^{k}=\mathbf{0}$ for all $k=2,\dots,m$. Then, the dynamics of $p^{1}$ simplifies to a single-virus system, which admits a NDFE by Proposition 4 in [13]. Therefore, in the case when $s(B^{1}-D^{1})>0$, the system (5) always admits an equilibrium of the form $(\tilde{p}^{1},\mathbf{0},\dots,\mathbf{0})$ with $\tilde{p}^{1}\gg\mathbf{0}$. ∎

III-B Dynamic Graph Structure

We can generalize the model from (4) to have dynamic graph structure as

[TABLE]

where $\beta^{k}_{ij}(t)$ is a function of time and the equation holds for $k=1,\dots,m$. We now provide a sufficient condition for global exponential stability of the DFE.

Theorem 4**.**

Suppose $B^{k}(t)$ is symmetric, piecewise continuous in $t$, and bounded, and $\sup_{t\geq 0}\lambda_{1}(B^{k}(t)-D^{k})<0$. Then the DFE is exponentially stable for virus $k$, with domain of attraction $\mathcal{D}$, in (6).

Proof.

Consider the Lyapunov function $V(p^{k})=\frac{1}{2}(p^{k})^{T}p^{k}$. For $p^{k}\neq 0$,

[TABLE]

The first inequality holds because $(P^{l}B^{k}(t))_{ij}\geq 0,\ \forall l,i,j,t$ by construction since each $p^{l}_{i}(t)$ is a probability. The second inequality holds by the Rayleigh-Ritz Theorem because $B^{k}(t)-D^{k}$ is symmetric. The last inequality holds by definition of the supremum. Therefore, the system converges exponentially fast to the origin by Theorem 8.5 in [28]. ∎

This result is a generalization of Theorem 1 in [19, 20].

We can also show exponential stability for the case when the infection rates are not symmetric and the underlying subgraphs are undirected, with some added assumptions.

Definition 1**.**

For a given virus $k$, assume that for all $t\geq 0$, there exist $c^{k}(t),\lambda^{k}(t)>0$ such that

[TABLE]

We then define

[TABLE]

Note that

[TABLE]

Theorem 5**.**

Consider the dynamics for virus $k$ in (8) with $B^{k}(t)$ continuously differentiable and $B^{k}(t)-D^{k}$ bounded, that is, there exists an $L>0$ such that $\|B^{k}(t)-D^{k}\|\leq L\ \forall t$. Assume that $\sup_{t\geq 0}s(B^{k}(t)-D^{k})<0$ and $\gamma^{k}_{1}$ in Definition 1 is finite. If $\sup_{t>0}\|\dot{B}^{k}(t)-D^{k}\|<\frac{1}{2(\gamma^{k}_{1})^{2}}$ or $\int_{t}^{t+T}\|\dot{B}^{k}(s)-D^{k}\|ds\leq\mu T+\alpha$ for small enough $\mu>0$, then the DFE is exponentially stable for virus $k$, with domain of attraction $\mathcal{D}$, in (6).

Proof.

Note that since $(P^{l}(t)B^{k}(t))_{ij}\geq 0\ \forall l,i,j$, by construction,

[TABLE]

Therefore, by Grönwall’s Inequality ([28]), the solution of the original system will be bounded above by the solution of the linear system. Thus by Lemma 2 in [20], the DFE is exponentially stable for virus $k$. ∎

Note that this theorem is a generalization of a single virus result provided in [20], where Lemma 2 in [20] is for a less general model; however the same arguments hold by replacing $BA(t)$ with $B^{k}(t)$ and $B\dot{A}(t)$ with $\dot{B}^{k}(t)$.

Theorem 6**.**

Consider the dynamics for virus $k$:

[TABLE]

Assume that

[TABLE]

for all $t_{0}\geq 0$, and for some $\nu>0$ there exists an $h>0$ such that

[TABLE]

for all $t\geq 0$ and some $\gamma,$ $0<\gamma\leq 1$. Assume further that

[TABLE]

for some negative scalar $\bar{\alpha}$ and for all $t_{0}\geq 0$,

[TABLE]

for all $t_{0}\geq 0$, and for all $i,j$ and $t\geq 0$ the perturbation

[TABLE]

Then the origin is exponentially stable for virus $k$.

Proof.

Since $(P(t)(B^{k}(t)+\Delta^{k}(t)))_{ij}\geq 0\ \forall i,j$ by (18) and Lemma 1,

[TABLE]

Therefore, by Grönwall’s Inequality ([28]), the solution of the original system will be bounded above by the solution of the linear system. Thus by Lemma 2 in [20], the origin is exponentially stable for virus $k$. ∎

This result says that if the linearized system is Hurwitz on the average (not strictly Hurwitz for all time, as in the other theorems up to this point), then the system converges to the DFE. This fact is useful in the control design in Section V.

IV NDFE for the Multi-Virus Case

There are a number of different epidemic equilibria. The simplest scenario is when one virus is in an epidemic state and the remaining viruses are eradicated.

Theorem 7**.**

Assume $\delta^{k}_{i}\geq 0$, for all $i,k$, and the matrices $B^{k}$ are non-negative and irreducible for all $k$. If for some $i\in[m]$, $s(B^{i}-D^{i})>0$ and $s(B^{k}-D^{k})\leq 0$ for all $k\neq i$, then (5) has two equilibria, the healthy state $(\mathbf{0},\dots,\mathbf{0})$, which is asymptotically stable with domain of attraction $\{(p^{1},\dots,p^{m})|p^{i}=\mathbf{0}\text{ and }p^{k}\in[0,1]^{n}\ \forall k\neq i\}$, and a unique epidemic state of the form $(\mathbf{0},\dots,\mathbf{0},\tilde{p}^{i},\mathbf{0},\dots,\mathbf{0})$ with $\tilde{p}^{i}\gg\mathbf{0}$, which is asymptotically stable with domain of attraction $\mathcal{D}\setminus\{(p^{1},\dots,p^{m})|p^{i}=\mathbf{0}\text{ and }p^{k}\in[0,1]^{n}\ \forall k\neq i\}$, with $\mathcal{D}$ defined in (6).

Note that this result is an extension of Theorem 3 in [13], and we present a sketch of the proof.

Sketched proof of Theorem 7: From the proof of Theorem 1, $p^{k}(t)$ will asymptotically converge to $\mathbf{0}$ as $t\rightarrow\infty$ for all initial values $(p^{1}(0),\dots p^{m}(0))\in\{(p^{1},\dots,p^{m})|p^{i}=\mathbf{0}\text{ and }p^{k}\in[0,1]^{n}\ \forall k\neq i\}$, for $k\neq i$. From (5),

[TABLE]

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

[TABLE]

with a vanishing perturbation $-\sum_{k\neq i}P^{k}(t)B^{i}p^{i}(t)$, which converges to $\mathbf{0}$ as $t\rightarrow\infty$. From Proposition 5 in [13], the autonomous system (19) will asymptotically converge to a unique epidemic state $(\mathbf{0},\dots,\mathbf{0},\tilde{p}^{i},\mathbf{0},\dots,\mathbf{0})$ for any $(p^{1}(0),\dots,p^{m}(0))\in\mathcal{D}\setminus\{(p^{1},\dots,p^{m})|p^{i}=\mathbf{0}\text{ and }p^{k}\in[0,1]^{n}\ \forall k\neq i\}$, with $\mathcal{D}$ defined in (6). ∎

Another possible NDFE is that of coexisting equilibrium, which is where more than one virus survives. We have the following interesting result similar to Theorem 7 in [13].

Theorem 8**.**

Consider the model in (4) with each virus propagating over the same strongly connected graph $\mathbb{G}$ with the corresponding adjacency matrix $A$, and each virus homogeneous in healing and infection rates, that is, for each $k\in[m]$ $\delta^{k}_{i}=\delta^{k}>0\ \forall i\in[n]$ and $\beta^{k}_{i}=\beta^{k}>0\ \forall i\in[n]$. Suppose that $s(A)>\frac{\delta^{1}}{\beta^{1}}=\cdots=\frac{\delta^{m}}{\beta^{m}}$. If $(\tilde{p}^{1},\dots,\tilde{p}^{m})$ with $\tilde{p}^{k}>\mathbf{0}\ \forall k\in[m]$ is an equilibrium of (4), then $\tilde{p}^{k}\gg\mathbf{0}\ \forall k\in[m]$ and $\tilde{p}^{i}=\alpha^{ik}\tilde{p}^{k}\ \forall i,k\in[m]$, for some constant $\alpha^{ik}>0$.

Proof: The homogeneity assumption on the infection rates allows us to factor $B^{k}=\beta^{k}A$, for each virus $k\in[m]$. To be an equilibrium of (4) the following must hold for all $k\in[m]$

[TABLE]

in which $(I-\tilde{P}^{1}-\cdots-\tilde{P}^{m})A$ is an irreducible Metzler matrix222 A matrix is Metzler if all of the off-diagonal components are non-negative. , since $\tilde{P}^{1}+\cdots+\tilde{P}^{m}$ is diagonal and $[\tilde{P}^{1}+\cdots+\tilde{P}^{m}]_{ii}<1$ for all $i\in[n]$, by Lemma 1. From Lemma 2 in [13], it must be true that $\tilde{p}^{k}\gg\mathbf{0}\ \forall k\in[m]$ and $\tilde{p}^{i}=\alpha^{ik}\tilde{p}^{k}\ \forall i,k\in[m]$, for some constant $\alpha^{ik}>0$. ∎

While the stability of the time–varying case has been explored in [20], the time–varying NDFE or epidemic limit cycle is an open problem, even for the single virus case. Some work has been done to show the existence of a periodic NDFE for a single virus switching system in [18].

V Antidote Control Formulation

Let us assume that for each agent, in addition to the healing rate, there is a control input $u_{i}(t)$ that acts as an additive boost to the healing rate. This implies that the controller can increase the agents’ ability to recover from the virus, which can be thought of as the administration of an antidote or some other type of treatment. This effect is portrayed in the model as

[TABLE]

We define $U(t)=\text{diag}(u)$ with $u=[u_{1}(t),\dots,u_{n}(t)]^{T}$. To simplify the discussion in this section, we assume that $B^{k}(t)$ is symmetric, piecewise continuous in $t$, and bounded $\forall t\geq 0$. Similar to the approaches in [21, 22, 23, 25, 29], we focus on minimizing the maximum eigenvalue of $B^{k}(t)-(D^{k}+U^{k}(t))$. Even though these control techniques are generally effective, we believe the approaches herein are more general and simpler, and therefore more scalable. Also, the assumption that our control input is additive to the base healing rate is novel and more sensible for the main motivating example, that is, every agent should have some inherent healing rate that should not be affected by the controller.

While the solutions to the following posed problems may not meet the conditions of Theorems 2 and 4, that is, they may not result in the maximum eigenvalues being less than zero, they push the system towards those conditions, consistent with the principle of the average being less than zero, presented in Theorem 6. And in practice, illustrated by simulation in the next section, these techniques reduce the spread of the epidemics. Under the aforementioned assumptions we can formulate the following optimization problem for each virus $k$, appealing to Theorems 2 and 4 depending on whether $B^{k}$ is constant or time dependent:

[TABLE]

given that $B^{k}(t)$ is symmetric for all $t\geq 0$.

From the Gershgorin Disc Theorem [30] it is clear that by sufficiently increasing the $u^{k}_{i}$’s, the conditions of Theorems 2 and 4 will be satisfied. Therefore we can relax the above optimization problem to obtain the following:

Problem 1**.**

[TABLE]

This is clearly a linear program and can easily be solved.

To make this a more compelling and realistic problem, we can impose a constraint on the number of agents that can be affected, which is a reasonable assumption because the cost of providing a low-dose treatment to all agents is higher than providing that same treatment dose to a few select members of the population (such as the sickest or most susceptible agents). Define the sparsity metric $\|\cdot\|_{0}$ as the number of the non-zero entries in its argument.

Employing the sparsity metric, we have the following problem, with a capacity constraint and a sparsity constraint:

[TABLE]

where $d^{k}$ is the maximum number of agents that can be treated for virus $k$. At first glance, the second and third constraints may seem redundant; however, the $\ell_{1}$ constraint limits the total amount of antidote that can be used while the sparsity constraint limits the number of agents that can be treated. The inclusion of the $\ell_{1}$ constraint prevents an infinite amount of antidote being administered to the limited number of agents allowed by the sparsity constraint.

It is well known that $\|\cdot\|_{0}$ is highly non-convex [31], making the above problem difficult to solve. Therefore, to solve it we employ another relaxation using the reweighted $\ell_{1}$ norm [32].

Definition 2**.**

The weighted $\ell_{1}$ norm is

[TABLE]

where $w_{i}$’s are positive and can be a constant or depend on time.

In view of this, we can rewrite the above problem as the following:

Problem 2**.**

[TABLE]

where $\kappa$ is a constant weighting factor.

An effective heuristic for the selection of the weights $w_{i}^{k}$’s in (21), proposed in [32], is, for some small $\epsilon>0$,

[TABLE]

For completeness, we include Algorithm 1, which explains the implementation of this heuristic to solve Problem 2. The notation Problem 2$(w_{k-1})$ indicates that $w_{k-1}$ is used for the weighted $\ell_{1}$ norm in the objective function of Problem 2 in the $k$th iteration. Employing this heuristic yields a good solution to Problem 2 but clearly is expensive, since it requires the calculation of multiple solutions. The effectiveness of this approach is illustrated in the following section via simulation.

VI Simulations

In this section we present a set of illuminating simulations of various competing virus models over static and time–varying graph structure networks. Due to limit of dimensions in color and size, for the simulations we will only have three competing viruses. Virus 1 is depicted by the color red ($r$), virus 2 is depicted by the color blue ($b$), and virus 3 is depicted by the color green ($g$). For all $i\in[n]$, the color at each time $t$ for agent $i$ is given by

[TABLE]

When $p^{1}_{i}(t)+p^{2}_{i}(t)+p^{3}_{i}(t)=0$, the color goes to black, indicating completely healthy, susceptible. These are used to facilitate the depiction of the parallel equilibrium ($\tilde{p}^{1}=\alpha^{2}\tilde{p}^{2}=\alpha^{3}\tilde{p}^{3}$), which will be shown by all nodes converging to the same color. For all $i\in[n]$, the diameter of the node representing 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 agent $i$. Therefore the color indicates the sickness each agent has and the diameter indicates how sick each agent is.

For systems that have three different subgraphs, viruses 1, 2, and 3 spread on the graphs depicted by gray, green, and pink edges, respectively. If all viruses spread on the same graph, the edges are gray.

The simulation in Figure 2 has three viruses spreading over the same time–varying graph. Similar to [20], the graph structure is determined by

[TABLE]

where $z_{i}(t)\in\mathbb{R}^{2}$ is the position of agent $i$, with $\hat{r}=10$. The agents have piece-wise constant drifts, that is,

[TABLE]

where $\phi(t)\in\mathbb{R}^{2}$ and is determined, for each dimension $l\in[2]$, by

[TABLE]

where the agents hover around a square, centered at some point $z_{c}$. The initial positions and $\phi$’s are chosen randomly. Each virus is homogeneous in infection rate. The first two viruses meet the assumptions of Theorem 4, while the maximum eigenvalue of the third virus fluctuates between being positive and negative. See Figure 3 for a plot of the maximum eigenvalues of the three-virus dynamics. Consistent with the theorem, the first two viruses are eradicated quite quickly. The third virus is also eliminated, but it takes a little longer. This eradication is illustrated in Figure 2b.

The simulation in Figure 4 meets the assumptions of Theorem 7, where $s(B^{1}-D^{1})>0$, and $s(B^{2}-D^{2})<0$ and $s(B^{3}-D^{3})<0$. Therefore the first virus, depicted in red, reaches an epidemic equilibrium, while the other two viruses are eradicated.

The simulation shown in Figure 5 meets the assumptions of Theorem 8, that is, the three viruses are each homogeneous, with $\frac{\delta^{1}}{\beta^{1}}=\frac{\delta^{2}}{\beta^{2}}=\frac{\delta^{3}}{\beta^{3}}$, and propagate over the same graph structure. There are $15$ agents and the initial conditions are given in Figure 5a. Consistent with the theorem, the system converges to a co-existing parallel equilibria.

We conclude with a simulation that implements the control techniques presented in Section V. Consider the single virus system in Figure 6. This system is homogeneous in infection rate, with $\beta=0.492$. We compare the system with no controller (on the left), a controller using Problem 1 (in the middle), and a controller that uses Algorithm 1 to solve Problem 2 iteratively with $\kappa=.05$ (on the right). The sum of the final probabilities of infection for all agents ($\sum_{i=1}^{n}p_{i}(100)$) for the three plots are 10.7, 4.92, and 3.6, respectively. Therefore, Algorithm 1 performed the best, however both had significant improvements over the uncontrolled simulation. The maximum eigenvalues of the three linearized systems are, from left to right, 1.893, 0.557, and 0.421; so none of the linearized systems are Hurwitz. Therefore, consistent with Theorem 3, the systems are all at NDFE. However, even though the control efforts do not completely eradicate the virus, they do mitigate its effect.

VII Conclusion

We have explored the competing multi-virus SIS model with several theorems exploring stability of the equilibria of the model for the static and time–varying graph cases. We have also proposed several control techniques that appeal to Theorems 2, 4, and 6, providing two efficient centralized antidote distribution/allocation protocols.

In future work we would like to explore more generic cases of co-existing epidemic states. Further, we would like to compare the techniques in Section V to other existing techniques. We would also like to implement the control techniques on large scale systems with at least tens of thousands of nodes.