\section{Introduction}
In February $1918$ a deadly influenza pandemic (popularly known as the Spanish flu) swept across the globe. It lasted until $1920$, and caused 
approximately
$50$ million deaths \cite{johnson2002updating}. Influenza viruses have continued to spread across the globe in recurring epidemics \cite{potter2001history}.
Given that the spread of infectious diseases has an enormous impact on society,
the study of spread
has been an active area of research since Bernoulli's seminal paper \cite{bernoulli1760essai}. The overarching goal
of these research directions is to
find conditions that would cause an epidemic to become eradicated, and leverage the knowledge of these conditions to design spread control strategies. To this end, various infection models have been proposed and studied in the literature;
susceptible-infected-susceptible (SIS), susceptible-infected-removed
(SIR), susceptible-exposed-infected-removed
(SEIR),
etc. In this paper,
we focus on SIS models.  

More specifically, we consider networked SIS models, in which a population of individuals is divided into subpopulations,
and the infection can spread both within and between these subpopulations. 
Networked SIS models have been studied extensively using discrete-time \cite{ahn2013global,wang2003epidemic,chakrabarti2008epidemic,BokharaieMTNS10} and continuous-time dynamics
\cite{van2008virus,khanafer2016stability,van2009virus}.
In the present paper, we will focus on continuous-time dynamics. Although both time-invariant and time-varying continuous-time models have been studied in the literature (see for instance \cite{pare2018epidemic,ogura2016stability} and references therein), we will restrict ourselves to the time-invariant case.

The spread represented by the SIS model has typically been understood as a consequence of
human contact. However, 
the spread of infectious diseases
can significantly worsen 
due to the presence of a shared resource. For example, waterborne pathogens spread via water distribution systems \cite{kough2015modelling} and droplet-transmitted pathogens spread via surfaces in public transit vehicles \cite{hertzberg2018behaviors}. Such observations have motivated the development of the susceptible-infected-water-susceptible (SIWS) model \cite{liu2019networked,pare2019multi} -- essentially, a networked continuous-time \emph{single-virus} SIS model that incorporates shared resources. For the SIWS model with a \emph{single} shared resource,
sufficient conditions for asymptotic convergence to the healthy state (where each subpopulation is infection-free, and the shared resource is contamination-free) have been provided  \cite[Theorem~1]{liu2019networked}. This result has been generalized to also account for \emph{multiple} shared resources
in \cite[Theorem~1]{pare2019multi}, where
virus transmission could be due to individual-to-individual contact, due to individual-to-resource contac
, and, in contrast to the single resource case, also due to resource-to-resource contact. However,
no theoretical guarantees for the endemic behavior (i.e., where the virus persists) of the SIWS model have been provided

\par The analysis of epidemic spread using SIWS models has been restricted to the single-virus cas
. While such analysis provides insights into how to battle
an epidemic,
they do not account for settings where multiple strains of a virus could be simultaneously active within a population.
In particular, it is possible that viral strains \emph{compete} with each other to infect the population: 
each individual
can be infected by one, and only one, of the multiple viral strains prevalent \cite{castillo1989epidemiological}, i.e., the strains impose \emph{cross-immunity}. Such a phenomenon is not only restricted to epidemics, but may also appear
in the context of incompatible  information transmission on social networks \cite{sahneh2014competitive}, or of two competing products within the same market \cite{li2015marketing}. It is well-known that competitive multi-virus propagation exhibits
rich behaviour in comparison to single-virus propagation \cite{newman2005threshold}. One possible outcome of competitive multi-virus propagation is \emph{coexistence} (i.e., multiple strains coexist in a population by infecting
separate fractions of each subpopulation), while another is \emph{competitive exclusion} (i.e., the spread parameters of one strain dominate those of the other strains, thereby causing those strains to become eradicated). In particular,
\cite{bremermann1989competitive} established that cross-immunity
between strains typically leads to competitive exclusion in single-population SIS models, whereas \cite{li2004coexist,liu2019analysis,pareautomatica,sahneh2014competitive}
provided conditions for coexistence in networked SIS models.
In particular, analysis of the various equilibria of a competing continuous-time time-invariant bi-virus model has been provided in \cite{liu2019analysis},
whereas a necessary and sufficient condition for
a coexisting equilibrium has been established \cite[Theorem~6]{pareautomatica}.
However, the results obtained in \cite{liu2019analysis,pareautomatica} are restrictive in the following sense:
i)
\cite[Theorems~6 and 7]{liu2019analysis} rely on the assumption that the spread parameters with respect to each virus is the \emph{same} for every subpopulation; and ii)
\cite[Theorem~6]{pareautomatica}
is reliant on the assumption that the set of spread parameters for each virus
is a scaled version of that of other viruses.
Moreover, none of these works account for the presence of a shared resource in the network. Thus, to the best of our knowledge, for the multi-virus SIWS model, a detailed analysis of the various equilibria and their stability properties remains missing in the literature. The present paper aims to fill this gap and in the process, for bi-virus SIS models, establish conditions for the existence of a coexisting equilibrium under less stringent assumptions.

The overarching goal is to develop 
strategies for the eradication of epidemics.
Several control strategies have been formulated in the context of networked SIS models:
optimal antidote distribution
\cite{mai2018distributed},
contact reduction \cite{theodorakopoulos2012selfish, tomovski2011simple}, etc. In particular, for a directed network with heterogeneous spread, \cite{enyioha2015distributed} provides a fully distributed Alternating Direction Method of Multipliers (ADMM) algorithm that allows for local computation of optimal investment required to boost the healing rate at each node. Note that this ADMM algorithm requires every node to communicate with its neighbors its local estimate of the full network, which
can be undesirable to do in practice.
A decentralized algorithm that involves disconnecting nodes and increasing the healing rates, subject to resource constraints, has been proposed in \cite{torres2016sparse}.
More recently, a distributed algorithm that, given resource limitations, guarantees the eradication of an epidemic with a \emph{specified} rate has been proposed in~\cite{mai2018distributed}. In the case of multiple cross-immune strains, eradication of all viruses via optimal antidote distribution is featured in \cite{pare2017multi}. In certain multi-strain epidemics, one strain can impose more severe symptoms than other strains \cite{lancet2020coronastrains}. If resource limitations do not permit us to eradicate all strains, we might prioritize eradicating a malignant strain first. The idea of eradicating one strain while sustaining another has been explored in \cite{watkins2016optimal}.
Note that the
algorithms provided in
\cite{liu2019analysis,enyioha2015distributed,torres2016sparse,mai2018distributed}, and the results in \cite{nowzari2016analysis,pareautomatica,watkins2016optimal,gracy2020asymptotic} account only for
SIS models.
















\subsection*{Paper Contributions}
The paper studies the multi-virus SIWS model. We establish existence, uniqueness and stability of certain equilibria, and develop mitigation strategies based on these results, as follows:
\begin{enumerate}[label=(\roman*)]
   
   
   
   
    \item
   
    Since the persistence of a virus has an enormous impact on society, it is important to know under what conditions a virus persists and, assuming at least one individual in a subpopulation is infected initially, if the virus can become endemic in the population. To this end, we provide a sufficient condition for the existence, uniqueness, and asymptotic stability of a single-virus endemic equilibrium; see Theorem~\ref{thm:equi}.
   
   
    \item \label{cont:suff} In the bi-virus case with (resp. without) a shared resource, as discussed previously, a natural question is whether both viruses can persist simultaneously. We establish a sufficient condition for the existence of a coexisting equilibrium;
    see Theorem~\ref{thm:joint_eq_exist_shared} (resp. Corollary~\ref{cor:joint_eq_exist}). 
    \item Similarly, in the bi-virus case with (resp. without) a shared resource, 
   
    it is also important to determine conditions that lead to competitive exclusion. Thus, we establish a necessary condition
    for the existence of a coexisting equilibrium; see Theorem~\ref{thm:noequi_single-virus} (resp. Corollary~\ref{cor:noequi_single-virus_noshared}).
   
     
     
     
     
   
   
    \item Utilizing our analytical results, we have developed two mitigation strategies: one
   
    for the stabilization (exponential or asymptotic) of the healthy state, and the other for the eradication of a malignant virus in the bi-virus case; see Theorems~\ref{thm:all_deltas_chosen_to_heal} and~\ref{thm:virus-as-vaccine}, respectively.
\end{enumerate}


Additionally, we present results on the stability and uniqueness of the healthy state, see Theorems~\ref{thm:expo},~\ref{thm:asymp} and~\ref{thm:uniqueness:healthy:state}. A preliminary version of some of the results in this paper
is scheduled for presentation
at the $3^{\rm{rd}}$ IFAC Workshop on Cyber-Physical \& Human Systems; see \cite{axel2020multi}. 

The aforementioned findings have several consequences. In an epidemiological context, 
the results on the healthy state 
shed light on the conditions that ensure the eradication of a virus in the population and the shared resource, and, therefore, aid in developing a mitigation strategy (see Section~\ref{sect:dist:control}). The results on the single-virus endemic equilibrium are useful
for learning the spread parameters (i.e., healing rate and infection rate), as discussed for the special case of SIS models in~\cite{pare2018analysis}.
One of the results concerning the coexisting equilibrium inspires a novel virus eradication strategy, where the goal is to push the dynamics towards the single-virus endemic equilibrium of the benign virus, as opposed to the healthy state (see Section~\ref{sect:dist:control}).
The results on the existence of a coexisting equilibrium
find applications 
not just in epidemiology, but also in opinion dynamics in social networks,
where it provides conditions under which
competing ideas (for instance, pro-Covid-19 restrictions versus anti-Covid-19 restrictions) coexist in different groups of the same population.

   
\subsection*{Paper Outline}
The paper is organized as follows:
We introduce the model
and formally present the main problems being investigated
in Section~\ref{sect:Model}. We collect the necessary technical background needed for establishing our main results in Section~\ref{sect:prelims}. The main results
are split up across Sections~\ref{sect:stab:analysis:healthystate
--\ref{sect:dist:control}. More specifically, Sections~\ref{sect:stab:analysis:healthystate
--\ref{sect:coexist} concern the analysis of equilibria, whereas mitigation strategies for eradicating or mitigating the spread of viruses are provided in Section~\ref{sect:dist:control}. We illustrate our
findings via simulations in Section~\ref{sect:simulations}. Finally, we summarize the main findings of our paper, and highlight relevant problems of future interest in Section~\ref{sect:conclusions}.

\subsection*{Notations}
We denote the set of real numbers by $\mathbb{R}$, and the set of nonnegative real numbers by $\mathbb{R}_+$. For any positive integer $n$, we use $[n]$ to denote the set $\{1,2,...,n\}$. The $i^{\rm{th}}$ entry of a vector $x$ is denoted by $x_i$. The element in the $i^{\rm{th}}$ row and $j^{\rm{th}}$ column of a matrix $M$ is denoted by $M_{ij}$. We use $\textbf{0}$ and $\textbf{1}$ to denote the vectors whose entries all equal $0$ and $1$, respectively, and use $I$ to denote the identity matrix, while the sizes of the vectors and matrices are to be understood from the context. For a vector $x$ we denote the square matrix with $x$ along the diagonal by $\diag(x)$. 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]$. Likewise, for any two real matrices $A, B \in \mathbb{R}^{n \times m}$, we write $A \geq B$ if $A_{ij} \geq B_{ij}$ for all $i \in [n]$, $j \in [m]$, and $A>B$ if $A \geq B$ and $A \neq B$.
For a square matrix $M$, we use $\sigma(M)$ to denote the spectrum of $M$, $\rho(M)$ to denote the spectral radius of $M$, and $s(M)$ to denote the largest real part among the eigenvalues of $M$, i.e., $s(M) = \max\{\rm{Re}(\lambda) : \lambda \in \sigma(M)\}$. Given a matrix $A$, $A \prec 0$ (resp. $ A\preccurlyeq 0 $) indicates that $A$ is negative definite (resp. negative semidefinite), whereas $A \succ 0$ (resp. $ A\succcurlyeq 0 $) indicates that $A$ is positive definite (resp. positive semidefinite).
We denote a subset by $P \subseteq Q$, a proper subset by $P \subset Q$, and set difference by $P \setminus Q$

\begin{figure}
\centering
    \scalebox{0.6}[0.6]{\includegraphics[width=\columnwidth]{multivirus_shared_model.png}}
    \caption{Visualization of the model for the case when $m=2$. An individual is either susceptible (S), infected with virus~$1$ (I$^1$), or infected with virus~$2$ (I$^2$).
    The shared resource (W) is contaminated by individuals infected with either virus, and in turn augments the corresponding infection rate.}
    \label{fig:model_visualized}
\end{figure}

\section{Problem Formulation}\label{sect:Model}
In this section, we detail a model of multi-viral spread across a population network with a shared resource.
We then formally state the problems being investigated. Finally, we introduce pertinent assumptions and definitions for use in later sections.

\subsection{The multi-virus SIWS model} \label{subsec:model}


Consider a population of individuals, subdivided into $n$
\emph{population nodes} in a network, with a resource $W$ being shared among some or all of the population nodes. Suppose that $m$ viruses are active in the population. An individual can become infected by a virus, either by coming into contact with an infected individual, or due to interaction with the (possibly) contaminated shared resource. We make the assumption that viral infection causes cross-immunity, meaning that an individual can be infected by no more than one virus at the same time. An infected individual can then recover, returning
to the susceptible state. The model is visualized from an individual's perspective in Figure~\ref{fig:model_visualized}. The spread of the $m$ viruses across the population can be represented by a multi-layer network $G$ with $m$ layers, where the vertices correspond to population nodes and the shared resource, and each layer contains a set of directed edges, $E^k$, specific to each virus~$k$.
There exists a directed edge from node $j$ to node $i$ in $E^k$ if an individual, infected by virus~$k$ in node $j$,
can directly infect individuals in node $i$.
Furthermore, the existence of a directed edge in $E^k$ from node $j$ (resp. from the shared resource $W$) to the shared resource $W$ (resp.  to the node $j$) signifies that the shared resource $W$ (resp. node $j$) can be contaminated with virus~$k$ by infected individuals in node $j$ (resp. by the shared resource $W$).
We say that the $k^{\rm{th}}$ layer of $G$ is strongly connected if there is a path via the directed edges in $E^k$ from each node, and from the shared resource, to every node, and to the shared resource. In real-world scenarios, it is often the case that viral epidemics can spread from each subpopulation to every other subpopulation,
in which case we assume that each layer is strongly connected.
\par Each population node $i$ contains $N_i$ individuals, with a birth rate $\mu_i$ equal to its death rate $\Bar{\mu}_i$. At any time $t \geq 0$, $S_i(t)$ is the number of susceptible individuals in node $i$, while $I_i^k(t)$ is the number of individuals infected by virus~$k$ in node $i$, while $S_i(t) + \textstyle \sum_{k=1}^m I_i^k(t) = N_i$. The rate at which individuals infected by virus~$k$ in node $j$ infect susceptible individuals in node $i$ is denoted by $\alpha_{ij}^k$, where $\alpha_{ij}^k = 0$ corresponds to the absence of a directed edge from node $j$ to node $i$ in $E^k$. In node $i$, individuals infected by virus~$k$ recover to the susceptible state at a rate $\gamma_i^k$. The shared resource contains a viral mass with respect to each virus~$k$, denoted by $W^k(t)$, representing the level of contamination at time $t \geq 0$. The viral mass of virus~$k$ grows at a rate proportional to all $I_i^k(t)$ scaled by their corresponding rates $\zeta_i^k$, and decays at a rate $\delta_w^k$. The resource-to-node infection rate to node $i$ with respect to virus~$k$ is denoted by $\alpha_{iw}^k$.
The time evolution of the number of susceptible and infected individuals (with respect to each virus~$k \in [m]$) in population node $i \in [n]$ is given by
\begin{equation} \label{eq:orig}
 \begin{split}
  \dot{S}_i(t) = & \, \mu_i N_i - \Bar{\mu}_i S_i(t) + \textstyle \sum_{k=1}^{m} \gamma_i^k I_i^k(t) \\
  &
  - \textstyle\sum_{k=1}^{m} \big{(} \alpha_{iw}^k W^k(t) - \textstyle \sum_{j=1}^{n} \alpha_{ij}^k \frac{I_j^k(t)}{N_i} \big{)} S_i(t), \\
  \dot{I}_i^k(t) = &- \, (\Bar{\mu}_i + \gamma_i^k) I_i^k(t) \\
  &+ \big{(} \alpha_{iw}^k W^k(t) + \textstyle \sum_{j=1}^{n} \alpha_{ij}^k \frac{I_j^k(t)}{N_i} \big{)} S_i(t), \\
  \dot{W}^k(t) = &- \delta_w^k W^k(t) +\textstyle \sum_{j=1}^{n} \zeta_j^k I_j^k(t).
 \end{split}
\end{equation}
We define new variables to simplify the system. Let
\begin{align*}
&p_i^k(t) \coloneqq \frac{I_i^k(t)}{N_i}, \, z^k(t) \coloneqq \frac{\delta_w^k W^k(t)}{\sum_{j=1}^{n} \zeta_j^k N_j}, \, \delta_i^k \coloneqq \gamma_i^k + \mu_i, \\
&\beta_{ij}^k \coloneqq \alpha_{ij}^k \frac{N_j}{N_i}, \, \beta_{iw}^k \coloneqq \frac{\alpha_{iw}^k}{\delta_w^k} \textstyle \sum_{j=1}^{n} \zeta_j^k N_j, c_i^k \coloneqq \frac{\zeta_i N_i}{\sum_{j=1}^{n} \zeta_j^k N_j},
\end{align*}
where the variables can be interpreted as follows: With respect to virus~$k$, $p_i^k(t)$ is the fraction of currently infected individuals in node $i$, $z^k(t)$ is a scaled contamination level in the shared resource, $\delta_i^k$ is the healing rate in node $i$, $\beta_{ij}^k$ is the node-to-node infection rate from node $j$ to $i$, scaled with respect to population ratios, $\beta_{iw}^k$ is a scaled resource-to-node infection rate to node $i$, and $c_i^k$ is a scaled node-to-resource contamination rate from node $i$.
Then, assuming that the birth rates and death rates are equal for each node,~\eqref{eq:orig} can be rewritten as
\begin{align} 
  \dot{p}_i^k(t) = - \delta_i^k p_i^k(t) + \big{(} &1 - \textstyle \sum_{l=1}^m p_i^l(t) \big{)} \nonumber \\ 
  &\times \big{(} \beta_{iw}^k z^k(t) + \textstyle \sum_{j=1}^{n} \beta_{ij}^k p_j^k(t) \big{)}, \nonumber\\
  \dot{z}^k(t) = \delta_w^k \big{(} - z^k(t) +& \textstyle \sum_{i=1}^{n} c_i^k p_i^k(t) \big{)}.  \label{eq:split}
\end{align}
Using vector notation,~\eqref{eq:split} can be rewritten as
\begin{equation} \label{eq:vec}
 \begin{split}
  \dot{p}^k(t) = & \, \Big{(} \big{(} I - \textstyle \sum_{l=1}^m \diag(p^l(t)) \big{)} B^k - D^k \Big{)} p^k(t) \\
  &\,\,+ \big{(} I - \textstyle \sum_{l=1}^m \diag(p^l(t)) \big{)} b^k z^k(t), \\
  \dot{z}^k(t) = & \, \delta_w^k \big{(} - z^k(t) + c^{k} p^k(t) \big{)}, \\
 \end{split}
\end{equation}
where 
$B^k$ is an $n \times n$-matrix with 
$\beta_{ij}^k$ as the $(i,j)^{\rm{th}}$ element,
$D^k$ is a diagonal $n \times n$-matrix with $D_{ii}^k = \delta_i^k$ for all $i \in [n]$, $b^k$ is a column vector with $\beta_{iw}^k$ as the $i^{\rm{th}}$ element, and $c^k$ is a row vector with $c_i^k$ as the $i^{\rm{th}}$ element.
To simplify notation further, we define
\begin{gather*}
 y^k(t)
 \coloneqq
 \begin{bmatrix}
 p^k(t)\\
 z^k(t)
 \end{bmatrix}
 ,
 y(t)
 \coloneqq
 \begin{bmatrix}
 y^1(t) \\
 \vdots \\
 y^m(t) \\
 \end{bmatrix}
 ,
 \normalsize B_w^k
 \coloneqq
 \begin{bmatrix}
 B^k & b^k \\
 \delta_w^k c^k & 0 
 \end{bmatrix}
 ,
\end{gather*}

\begin{gather*}
 X(y^k(t))
 \coloneqq
 \begin{bmatrix}
 \diag(p^k(t))  & 0 \\
 0 & 0 
 \end{bmatrix}
 ,
 \,\,
 D_w^k
 \coloneqq
 \begin{bmatrix}
 D^k & 0 \\
 0 & \delta_w^k 
 \end{bmatrix}
 .
\end{gather*}
With these variables in place, 
we can rewrite~\eqref{eq:vec} as
\begin{equation}\label{eq:yk}

 \dot{y}^k(t) =\big{(} - D_w^k + (I - \textstyle \sum_{l=1}^m X(y^l(t)))B_w^k \big{)}y^k(t).
\end{equation}
Defining $A_w^k(y(t)) \coloneqq \big{(} - D_w^k + (I - \textstyle \sum_{l=1}^m X(y^l(t)))B_w^k \big{)}$, the dynamics of the system of all $m$ viruses are given by
\begin{gather} \label{eq:full}
 \dot{y}(t)
 =
 \begin{bmatrix}
 A_w^1 \big{(} y(t) \big{)} & 0 & \dots & 0 \\
 0 & A_w^2 \big{(} y(t) \big{)} & \dots & 0 \\
 \vdots & \vdots & \ddots & \vdots \\
 0 & 0 & \dots & A_w^m \big{(} y(t) \big{)}
 \end{bmatrix}
 y(t).
\end{gather}
We say that virus~$k \in [m]$ is \emph{eradicated}, or in its \emph{eradicated state}, if $y^k = \textbf{0}$, which is clearly an equilibrium of~\eqref{eq:yk}.
When considering the system of $m$ viruses~\eqref{eq:full}, we say that the system is in the \emph{healthy state} if
all viruses are eradicated, i.e., $y = \textbf{0}$.
If~\eqref{eq:full} has an endemic (non-zero) equilibrium, it can belong to one of two types: \emph{single-virus endemic equilibrium}, where $y^k>\textbf{0}$ for some $k \in [m]$ and all other viruses in $[m]$ are eradicated; or \emph{coexisting equilibrium}, where $y^k>\textbf{0}$ for multiple $k \in [m]$.
\subsection{Problem Statements}\label{subsec:pb:statement}
For the model~\eqref{eq:full}, we formally state the problems being investigated in this paper.
\begin{enumerate}[label=(\roman*)]
    \item \label{q1} Under what conditions does $y^k(t)$ converge exponentially, or asymptotically, to its eradicated state, i.e., $y^k=\textbf{0}$, for some $k \in [m]$?
    \item \label{q2} For a single-virus setup, i.e., $m=1$, under what conditions does the system have a unique single-virus endemic equilibrium, $y^* > \textbf{0}$, and under such conditions, does the system converge asymptotically to $y^*$ from any non-zero initial condition?
   
    \item \label{q3} What is a necessary and sufficient condition for the healthy state, i.e., $y=\textbf{0}$, to be the unique equilibrium?
    \item \label{q4} For a bi-virus setup, i.e., $m=2$, under what conditions does the system have a coexisting equilibrium, i.e., $(\hat{y}^1, \hat{y}^2)$ such that $\hat{y}^1>\textbf{0}$ and $\hat{y}^2>\textbf{0}$?
    
    \item \label{q5} For a bi-virus setup, under what conditions does the system \emph{not} have a coexisting equilibrium, i.e., $(\hat{y}^1, \hat{y}^2)$ such that $\hat{y}^1>\textbf{0}$ and $\hat{y}^2>\textbf{0}$?
    
    \item \label{q6} How can the healing rates, i.e., $\delta_i^k$, be chosen to ensure that the system converges exponentially, or asymptotically, to the healthy state, i.e., $y=\textbf{0}$?
    \item \label{q7} For a bi-virus setup, how can the healing rates of virus~$2$, i.e., $\delta_i^2$, be chosen to ensure that the system converges to the single-virus endemic equilibrium of virus~$1$?
\end{enumerate}

   
    
    
   
   
     
     
      
      
   
   

\noindent Before we address these questions, we point out two connections between the considered setup and the existing literature.
\begin{rem} \label{rem:m=1}
If $m=1$, 
system~\eqref{eq:yk} coincides with the single-virus model proposed in \cite{liu2019networked
~$\blacksquare$
\end{rem}
Hereafter, when a single-virus system is considered, we drop the superscripts identifying the virus from all variables:

\noindent $B_w\coloneqq B_w^1$, $D_w\coloneqq D_w^1$, etc.





\begin{rem} \label{rem:no_W}
Note that if $b^k = \textbf{0}$ for all $k \in [m],$ the influence of the shared resource on the population is nullified. Then, the multi-virus dynamics of the population nodes, $p^k(t)$, in~\eqref{eq:full} are equivalent to the time-invariant multi-virus setup of \cite{pare2017multi}.~$\blacksquare$
\end{rem}
\subsection{Positivity assumptions}

In order for~\eqref{eq:full} to be well-defined and realistic, we make
the following assumption.
\begin{assum} \label{assum:base}
Suppose that $\delta_i^k>0, \delta_w^k > 0, \beta_{ij}^k \geq 0, \beta_{iw}^k \geq 0$ and $c_i^k \geq 0$ for all $i, j \in [n]$ and $k \in [m]$, with $c_l^k > 0$ for at least one $l \in [n]$.~$\blacksquare$ 
\end{assum}
Note that if Assumption~\ref{assum:base} holds, then for all $k \in [m]$, $B_w^k$ is a nonnegative matrix and $D_w^k$ is a positive diagonal matrix
Moreover, recall that a square matrix $M$ is said to be irreducible if, replacing the non-zero elements of $M$ with ones and interpreting it as an adjacency matrix, the corresponding graph is strongly connected. Then, noting that non-zero elements in $B_w^k$ represent directed edges in the set $E^k$, we see that $B_w^k$ is irreducible whenever the $k^{\rm{th}}$ layer of the multi-layer network $G$ is strongly connected.

Thanks to
Assumption~\ref{assum:base}, we can
restrict our analysis to the sets $\mathcal{D} \coloneqq \{y(t): p^k(t) \in [0,1]^n, z^k(t) \in [0, \infty) \, \forall k \in [m]\}$ and $\mathcal{D}^k \coloneqq \{y^k(t): p^k(t) \in [0,1]^n, \, z^k(t) \in [0, \infty)\}$.
Since $p_i^k(t)$ is to be interpreted as a fraction of a population, and $z^k(t)$ is a nonnegative quantity, these sets represent the
sensible domain of the
system. 
That is,
if $y^k(t)$ takes values outside of $\mathcal{D}^k$, then those values would lack physical meaning. The following lemma shows that once $y(t)$ enters $\mathcal{D}$,
it never leaves $\mathcal{D}$.
\begin{lem}\label{lem:pos}
Let Assumption~\ref{assum:base} hold. Then $\mathcal{D}$ is positively invariant with respect to~\eqref{eq:full
~$\blacksquare$
\end{lem}
\textit{Proof:}
Consider $y(t) \in \mathcal{D}$. If $p_i^k(t) = 1$, then $\dot{p}_i^k(t) < 0$, so if $p_i^k(0) \leq 1$ then $p_i^k(t) \leq 1$, for all $t \geq 0$, $k \in [m], i \in [n]$. Further, if $p_i^k(t) = 0$ and $y(t) \geq \textbf{0}$, then $\dot{p}_i^k(t) \geq 0$, for all $t \geq 0$. Similarly, if $z^k(t) = 0$ and $y(t) \geq \textbf{0}$, then $\dot{z}^k(t) \geq 0$. Thus, $y(t) \geq \textbf{0}$ if $y(0) \geq \textbf{0}$, for all $t \geq 0$.~$\square$
\section{Preliminaries} \label{sect:prelims}
In this section, we recall some preliminary results, pertinent to the analysis of system~\eqref{eq:full}. A real square matrix is said to be Metzler if all elements outside the diagonal are nonnegative. We require the following results for Metzler matrices.
\begin{lem} \label{lem:metz_irreduc}
\cite[Lemma A.1]{khanafer2016stability} Suppose that $M$ is an irreducible Metzler matrix such that $s(M) = 0$. Then there exists a positive diagonal matrix $P$ such that $M^T P + P M \preccurlyeq 0$.~$\blacksquare$ \end{lem}
\begin{lem} \label{lem:eigspec}
\cite[Proposition~1]{liu2019analysis} Suppose that $\Lambda$ is a negative diagonal matrix and $N$ is an irreducible nonnegative matrix. 
Let $M$ be the irreducible Metzler matrix $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$.~$\blacksquare$
\end{lem}
We will also be making use of the following variants of the Perron-Frobenius theorem for irreducible matrices.
\begin{lem} \label{lem:perron_frob}
\cite[Chapter 8.3]{meyer2000matrix} \cite[Theorem~2.7]{varga1999matrix}
Suppose that $N$ is an irreducible nonnegative matrix. Then,
\begin{enumerate}[label=(\roman*)]
    \item $r = \rho(N)$ is a simple eigenvalue of $N$. \label{item:perfrob_simpleeig}
    \item There is an eigenvector $\zeta \gg \textbf{0}$ corresponding to the eigenvalue $r$. \label{item:perfrob_pos_exists}
    \item $x > \textbf{0}$ is an eigenvector of $N$ only if $Nx = rx$ and $x \gg \textbf{0}$.
    \label{item:perfrob_pos_necess}
    \item If $A$ is a nonnegative matrix such that $A < N$, then $\rho(A) < \rho(N)$. \label{item:perfrob_matrix_ineq}~$\blacksquare$
\end{enumerate}
\end{lem}
\begin{lem} \label{lem:perron_frob_metz}
\cite[Lemma~2.3]{varga1999matrix} Suppose that $M$ is an irreducible Metzler matrix. Then $r = s(M)$ is a simple eigenvalue of $M$, and a corresponding eigenvector is $\zeta \gg \textbf{0}$.~$\blacksquare$
\end{lem}
\noindent Let $\mathbb{G} \subseteq \mathbb{R}^n$. Consider 
\vspace{-1ex}
\begin{equation} \label{eq:genericsys}
    \dot{x}(t) = f(x),
\end{equation}

\vspace{-1ex}

\noindent where $f : \mathbb{G} \rightarrow \mathbb{R}^n$ is a locally Lipschitz map.
\begin{prop} \label{prop:asymp_domatt}
Let $\textbf{0}$ be an equilibrium of~\eqref{eq:genericsys} and $\mathbb{E} \subseteq \mathbb{G}$ be a positively invariant and connected set with respect to~\eqref{eq:genericsys}, containing $\textbf{0}$. Let $V: \mathbb{E} \rightarrow \mathbb{R}$ be a continuously differentiable and positive definite function, such that $\dot{V}(x(t))$ is negative definite.
Further, let it hold that
\begin{equation} \label{eq:bounded_contours}
\{ x : V(x) \leq c\} \cap \mathbb{E}
\end{equation}
is a bounded set for any constant $c > 0$, and is equal to $\mathbb{E}$ as $c \rightarrow \infty$. Then the equilibrium $\textbf{0}$ is asymptotically stable, with domain of attraction containing $\mathbb{E}$.~$\blacksquare$
\end{prop}
This proposition can be proven using \cite[Theorem~4.1]{khalil2002nonlinear}.
The following lemma pertains to system~\eqref{eq:full}, providing a constraint on any endemic equilibrium.
\begin{lem} \label{lem:equi_non-zero_nonone}
Consider system~\eqref{eq:full} under Assumption~\ref{assum:base}. Suppose, for all $k \in [m]$, that $B_w^k$ is irreducible. If $y = (y^1, \dots, y^m) \in \mathcal{D}$ is an equilibrium of~\eqref{eq:full}, then, for each $k \in [m]$, either $y^k = \textbf{0}$, or $\textbf{0} \ll y^k \ll \textbf{1}$. Moreover,
$\textstyle \sum_{k=1}^m y^k \ll \textbf{1}$.~$\blacksquare$
\end{lem}
 \textit{Proof:} See Appendix.~$\square$

Lemma~\ref{lem:equi_non-zero_nonone} states that when the underlying network is strongly connected, any endemic equilibrium involves each active virus infecting a separate fraction of each population node, and contaminating the shared resource to some degree. 






\section{Analysis of the eradicated state of a virus}\label{sect:stab:analysis:healthystate}
In this section, we present
sufficient conditions for the exponential (resp. asymptotic)
stability of the eradicated state of a virus.
The key condition is found through eigenvalue analysis of $(B_w^k - D_w^k)$, as seen in the following theorem.
\begin{thm} \label{thm:expo}
Consider the SIWS model~\eqref{eq:yk} under Assumption~\ref{assum:base} 
with
$y(0) \in \mathcal{D}$. Suppose that for some virus~$k \in [m]$ we have $s(B_w^k - D_w^k) < 0$. Then the eradicated state of virus~$k$ is exponentially stable, with domain of attraction containing~$\mathcal{D}^k$.~$\blacksquare$
\end{thm}
\textit{Proof:} See Appendix.~$\square$ 


Theorem~\ref{thm:expo} states that if
the linearized state matrix of virus~$k$
is Hurwitz, then, for all initial conditions in the sensible domain, virus~$k$ is eradicated exponentially fast.
Theorem~\ref{thm:expo} answers the first part of question~\ref{q1} in Section~\ref{subsec:pb:statement}. 
With respect to \cite[Proposition~2]{liu2019analysis}, Theorem~\ref{thm:expo} is an improvement in the sense that it
holds globally (on the sensible domain), and accounts for the multi-virus case, whereas
\cite[Proposition~2]{liu2019analysis} established \emph{local} exponential stability for the \emph{single-virus} case

Theorem~\ref{thm:expo} indeed guarantees exponential eradication of virus~$k$, however, the condition is quite strict. For certain viruses it suffices to know whether or not the virus will be eradicated, but the \emph{speed} with which this eradication takes place is of less importance.
Indeed, it turns out that a relaxation of the strict inequality of the eigenvalue condition in Theorem~\ref{thm:expo} guarantees asymptotic eradication of a virus, as stated in the following theorem.


\begin{thm} \label{thm:asymp}
Consider the SIWS model~\eqref{eq:yk} under Assumption~\ref{assum:base} with $y(0) \in \mathcal{D}$. Suppose that for some virus~$k \in [m]$ we have $s(B_w^k - D_w^k) \leq 0$, and that the matrix $B_w^k$ is irreducible. Then the eradicated state of virus~$k$ is asymptotically stable, with domain of attraction containing $\mathcal{D}^k$.~$\blacksquare$
\end{thm}
\textit{Proof:} See Appendix.~$\square$\\
Assuming that the $k^{\rm{th}}$ layer of the graph is strongly connected, Theorem~\ref{thm:asymp} states that if the largest real part of any eigenvalue of the linearized state matrix of virus~$k$
is non-positive, then, for all initial conditions in the sensible domain, virus~$k$ is eradicated asymptotically. Theorem~\ref{thm:asymp}
answers the second part of question~\ref{q1} in Section~\ref{subsec:pb:statement}. 
Observe that  
for the single-virus case ($m=1$), Theorem~\ref{thm:asymp} 
improves \cite[Theorem~1]{liu2019networked} by relaxing the requirements $\beta_{iw} > 0$ and $c_i > 0$ for all $i \in [n]$.


\begin{rem}[Epidemiological Interpretation]
\label{rem:asymp:stab}
Observe that, due to Lemma~\ref{lem:eigspec}, the conditions in Theorem~\ref{thm:expo} (Theorem~\ref{thm:asymp})
are equivalent to $\rho((D_w^k)^{-1} B_w^k) < 1$ ($\rho((D_w^k)^{-1} B_w^k) \leq 1$).
The fact that Theorem~\ref{thm:expo} (resp. Theorem~\ref{thm:asymp}) guarantees exponential (asymptotic) eradication of a virus~$k$ whenever $\rho((D_w^k)^{-1} B_w^k) < 1$ ($\rho((D_w^k)^{-1} B_w^k) \leq 1$) is consistent with an interpretation of $\rho((D_w^k)^{-1} B_w^k)$ as the basic reproduction number of the virus in the network, typically denoted by $\mathcal{R}_0$. In epidemiology, this is the average number of secondary infections caused by an infected individual, before recovering.
Thus, Theorem~\ref{thm:expo} (Theorem~\ref{thm:asymp}) states that whenever the basic reproduction number of a virus
is strictly less than (less than or equal to) one, the virus will exponentially (asymptotically) converge to its eradicated state.~$\blacksquare$
\end{rem
\section{Persistence of a virus}\label{sec:persist}
In this section, we study the possibility of viruses persisting in the network, corresponding to non-zero equilibria of~\eqref{eq:full}. Naturally, the persistence of a virus must follow from the violation of
the conditions of Theorem~\ref{thm:asymp}.
The resulting behavior is detailed in the rest of this section.
Before we proceed, we state the following result:

\begin{lem} \label{lem:pos_never_zero}
Let Assumption~\ref{assum:base} hold. Then $\mathcal{D} \setminus \{ \textbf{0} \}$ is positively invariant with respect to system~\eqref{eq:full}.~$\blacksquare$
\end{lem}
 \textit{Proof:} See Appendix.~$\square$


With Lemma~\ref{lem:pos_never_zero} in place, we have the following theorem for a single-virus system, guaranteeing existence of a unique, asymptotically stable, single-virus endemic equilibrium when the eigenvalue condition in Theorem~\ref{thm:asymp} is violated
\begin{thm} \label{thm:equi}
Consider the SIWS model~\eqref{eq:yk} under Assumption~\ref{assum:base} with $m=1$. Suppose that $B_w$ is irreducible and $s(B_w - D_w) > 0$. Then there exists a unique single-virus endemic equilibrium $\Tilde{y} \in \mathcal{D}$, with $\textbf{0} \ll \Tilde{y} \ll \textbf{1}$, and it is asymptotically stable with domain of attraction containing $\mathcal{D} \setminus \{ \textbf{0} \}$.~$\blacksquare$
\end{thm}
\textit{Proof:} See Appendix.~$\square$

For a single-virus system, Theorem~\ref{thm:equi} states that when the eigenvalue condition in Theorem~\ref{thm:asymp} is violated, then as long as some viral infection is present initially, the viral spreading process will converge to a unique infection ratio in each population node and a unique contamination level in the shared resource. This answers question~\ref{q2} in Section~\ref{subsec:pb:statement}.
\begin{rem}[Epidemiological Interpretation] \label{rem:endemic:equi}
Applying Lemma~\ref{lem:eigspec}, we see that the conditions of Theorem~\ref{thm:equi} are equivalent to $\rho((D_w)^{-1}B_w) > 1$. This is again consistent with the interpretation of $\rho((D_w)^{-1}B_w)$ as the basic reproduction number of the virus, since a persisting virus should have a basic reproduction number greater than one.~$\blacksquare$
\end{rem}
Theorem~\ref{thm:equi} establishes the existence, uniqueness and asymptotic stability of a single-virus endemic equilibrium, extending \cite[Theorem~2.4.]{fall2007epidemiological}
to the setting with a shared resource, whereas \cite{liu2019networked} 
illustrated this extension in simulations, without providing theoretical guarantees.
We can partially extend Theorem~\ref{thm:equi} to the multi-virus case, specifically the existence and uniqueness of a single-virus endemic equilibrium for a virus, resulting in the following proposition.
\begin{prop} \label{prop:necessity}
Consider the SIWS model~\eqref{eq:full} under Assumption~\ref{assum:base}. For each $k \in [m]$ such that $B_w^k$ is irreducible and $s(B_w^k - D_w^k) > 0$, there is a unique single-virus endemic equilibrium $(\textbf{0}, \dots, \Tilde{y}^k, \dots, \textbf{0})$ in $\mathcal{D}$, with $\textbf{0} \ll \Tilde{y}^k \ll \textbf{1}$.~$\blacksquare$
\end{prop}
\textit{Proof:} See Appendix.~$\square$

Proposition~\ref{prop:necessity} states that each virus violating the eigenvalue condition in Theorem~\ref{thm:asymp} has a unique single-virus endemic equilibrium.
This result is unsurprising, since nullifying the other viruses in the system reduces it to a single-virus system. However, note that Proposition~\ref{prop:necessity} does not say anything about the stability of these equilibria. Nor does the proposition say anything about the case of more than one virus persisting in the population simultaneously.
Due to Theorem~\ref{thm:asymp} and Proposition~\ref{prop:necessity}, we have the following necessary and sufficient condition for the healthy state to be the unique equilibrium
\begin{thm}\label{thm:uniqueness:healthy:state}
Consider the SIWS model~\eqref{eq:full} under Assumption~\ref{assum:base}. Suppose that, for all $k \in [m]$, $B_w^k$ is irreducible. Then the healthy state is the unique equilibrium in $\mathcal{D}$ if, and only if, for all $k \in [m]$, $s(B_w^k - D_w^k) \leq 0$.~$\blacksquare$
\end{thm}
\noindent
Theorem~\ref{thm:uniqueness:healthy:state} states that as long as the largest real part of the eigenvalue of the linearized state matrix of \emph{each} virus is non-positive, the healthy state is the only equilibrium of~\eqref{eq:full}.
Theorem~\ref{thm:uniqueness:healthy:state} answers question~\ref{q3} in Section~\ref{subsec:pb:statement}. 
Note that Theorem~\ref{thm:uniqueness:healthy:state} extends \cite[Theorem~1]{liu2019analysis} to the setting with more than two viruses and a shared resource, 
albeit under the assumption that the healing rate of each agent with respect to each virus is strictly positive.




\section{Coexistence of viruses}\label{sect:coexist}

Beyond the single-virus endemic equilibria from 
Section~\ref{sec:persist}, we would like to know when multiple viruses can persist in a population simultaneously, corresponding to a coexisting equilibrium.
The following theorem makes use of a particular
eigenvalue condition to show the existence of a coexisting
equilibrium in a bi-virus system. Before proceeding to the statement of the theorem, recall from Proposition~\ref{prop:necessity} that, if i) $B_w^1$, $B_w^2$ are irreducible, ii) $s(B_w^1-D_w^1) > 0$, and iii) $s(B_w^2-D_w^2) > 0$, then there are exactly two single-virus endemic equilibria, namely
\begin{equation}\label{eq:endemicm2}
    (\Tilde{y}^1, \textbf{0}) \text{ and } (\textbf{0}, \Tilde{y}^2). 
\end{equation}
Moreover, we know that $\textbf{0} \ll \Tilde{y}^1 \ll \textbf{1}$ and $\textbf{0} \ll \Tilde{y}^2 \ll \textbf{1}$.



\begin{thm} \label{thm:joint_eq_exist_shared}
Consider the SIWS model~\eqref{eq:full} under Assumption~\ref{assum:base} with $m=2$. Suppose that $B_w^1$ and $B_w^2$ are irreducible matrices,
and that $s(B_w^1 - D_w^1) > 0$ and $s(B_w^2 - D_w^2) > 0$. 
If $s(- D_w^1 + (I - X(\Tilde{y}^2))B_w^1) > 0$ and $s(- D_w^2 + (I - X(\Tilde{y}^1))B_w^2) > 0$, with $\Tilde{y}^1,\Tilde{y}^2$ defined in~\eqref{eq:endemicm2},
then there exists at least one coexisting equilibrium $(\hat{y}^1, \hat{y}^2) \gg \textbf{0}$ in $\mathcal{D}$ such that $\hat{y}^1 + \hat{y}^2 \leq \textbf{1}$.~$\blacksquare$
\end{thm}























\textit{Proof:} See Appendix.~$\square$


\noindent
With each virus satisfying the condition for the existence of its single-virus endemic equilibrium,
Theorem~\ref{thm:joint_eq_exist_shared} states that if, for each virus, the largest real part of any eigenvalue of the matrix of the dynamics linearized around the single-virus endemic equilibrium of the other virus
is positive, then both the viruses can \emph{simultaneously} infect separate fractions of each population node.
Theorem~\ref{thm:joint_eq_exist_shared} answers question~\ref{q4} in Section~\ref{subsec:pb:statement}

\begin{rem} [Epidemiological Interpretation] \label{rem:coexist:equi}
Applying Lemma~\ref{lem:eigspec} to the conditions in Theorem~\ref{thm:joint_eq_exist_shared}, we see that they are equivalent to having $\rho((I-X(\Tilde{y}^2))(D_w^1)^{-1}B_w^1)>1$ and $\rho((I-X(\Tilde{y}^1))(D_w^2)^{-1}B_w^2)>1$. This is consistent with an interpretation of $\rho((I-X(\Tilde{y}^2))(D_w^1)^{-1}B_w^1)$ and $\rho((I-X(\Tilde{y}^1))(D_w^2)^{-1}B_w^2)$ as the invasion reproduction numbers of virus~1 invading virus~2 and virus~2 invading virus~1, respectively. The invasion reproduction number is defined for an invading pathogen, introduced into a setting with another, endemic pathogen at equilibrium. It is defined as the average number of secondary 
infections caused by an individual infected by the invading pathogen, at the time of introduction \cite{porco1998designing}. In line with this interpretation, Theorem~\ref{thm:joint_eq_exist_shared} shows that coexistence is possible whenever both invasion reproduction numbers are greater than one.~$\blacksquare$
\end{rem}
While conditions that guarantee existence of coexisting equilibria may be found in \cite{li2004coexist}
and \cite{pareautomatica}, these references do not account for the presence of a shared resource.
In order to compare
the results in \cite{li2004coexist,pareautomatica} with Theorem~\ref{thm:joint_eq_exist_shared}, we particularize our model to the setting without a shared resource.
Specifically,
\eqref{eq:yk} reduces to
\vspace{-1ex}
\begin{equation} \label{eq:nosharedsys}
\dot{p}^k(t) =
\Big{(} - D^k + \big{(} I - \textstyle \sum_{l=1}^m \diag(p^l(t)) \big{)} B^k \Big{)} p^k(t).
\end{equation} 

\vspace{-1ex}

\noindent Furthermore, we can let $p(t) \coloneqq [p^1(t),  \dots , p^m(t)]^T$. Then, with $A^k(p(t)) \coloneqq \big{(}- D^k + (I - \textstyle \sum_{l=1}^m \diag(p^l(t)))B^k \big{)}$, the dynamics of $p(t)$ are given by
\begin{gather} \label{eq:fullnosharedsys}
 \dot{p}(t)
 = 
 \begin{bmatrix}
 A^1 \big{(} p(t) \big{)} & 0 & \dots & 0 \\
 0 & A^2 \big{(} p(t) \big{)} & \dots & 0 \\
 \vdots & \vdots & \ddots & \vdots \\
 0 & 0 & \dots & A^m \big{(} p(t) \big{)}
 \end{bmatrix}
 p(t).
\end{gather}
Assumption~\ref{assum:base}, particularized for the setting without a shared resource, is given as follows:
\begin{assum} \label{assum:noshared}
Suppose that $\delta_i^k>0$ and $\beta_{ij}^k \geq 0$ for all $i, j \in [n]$ and $k \in [m]$.~$\blacksquare$
\end{assum}
Applying Proposition~\ref{prop:necessity}
to the bi-virus setting without a shared resource, we see that if i) $B^1$ and $B^2$ are irreducible, ii) $s(B^1-D^1) > 0$, and iii) $s(B^2-D^2) > 0$, then there exist exactly two single-virus endemic equilibria, namely
\begin{equation}\label{eq:endemicm2_noshared}
    (\Tilde{p}^1, \textbf{0}) \text{ and } (\textbf{0}, \Tilde{p}^2), 
\end{equation}
such that $\textbf{0} \ll \Tilde{p}^1 \ll \textbf{1}$ and $\textbf{0} \ll \Tilde{p}^2 \ll \textbf{1}$. Hence, we have the following corollary to Theorem~\ref{thm:joint_eq_exist_shared}, for the setting without a shared resource.
\begin{cor} \label{cor:joint_eq_exist}
Consider the SIS model~\eqref{eq:fullnosharedsys} under Assumption~\ref{assum:noshared} with $m=2$. 
Suppose that $B^1$ and $B^2$ are irreducible matrices, and that $s(B^1 - D^1) > 0$ and $s(B^2 - D^2) > 0$. 
If $s(- D^1 + (I - diag(\Tilde{p}^2))B^1) > 0$ and $s(- D^2 + (I - diag(\Tilde{p}^1))B^2) > 0$, with $\Tilde{p}^1,\Tilde{p}^2$ defined in~\eqref{eq:endemicm2_noshared}, then there exists at least one coexisting equilibrium $(\hat{p}^1, \hat{p}^2) \gg \textbf{0}$ in $[0,1]^{2n}$ such that $\hat{p}^1 + \hat{p}^2 \leq \textbf{1}$.~$\blacksquare$
\end{cor}
Corollary~\ref{cor:joint_eq_exist} guarantees existence of a coexisting equilibrium in the bi-virus setup. It is an improvement of a similar result in \cite[Theorem~5.2]{li2004coexist}, wherein the same is established for $n=2$. However, \cite[Theorem~5.2]{li2004coexist}
establishes uniqueness of the coexisting equilibrium, whereas Corollary~\ref{cor:joint_eq_exist} provides no such guarantees.
Note that \cite[Theorem~6]{pareautomatica}, particularized for $m=2$, establishes the existence of not just one but infinitely many coexisting equilibria, thus implying that a coexisting equilibrium in a bi-virus setup
is not necessarily unique.  However, it turns out that the conditions in Corollary~\ref{cor:joint_eq_exist} do \emph{not} coincide with the conditions in \cite[Theorem~6]{pareautomatica}, as discussed in the following remark.


\begin{rem} \label{rem:not_compatible_conditions}
Suppose that, for $m=2$, the conditions in \cite[Theorem~6]{pareautomatica} are satisfied. That is, for some $v > 0$, $B^1$ is an irreducible nonnegative matrix with $B^1 = v B^2$, and $D^1$ is a positive diagonal matrix with $D^1 = v D^2$, $s(B^1 - D^1) > 0$ and $s(B^2 - D^2) > 0$. Let $\Tilde{p}^1$ and $\Tilde{p}^2$ be the unique single-virus endemic equilibria defined in~\eqref{eq:endemicm2_noshared}.
It follows from~\eqref{eq:nosharedsys}, and the fact that $D^1, D^2$ are invertible, 
that
$\Tilde{p}^1, \Tilde{p}^2$ fulfill
\begin{equation} \label{eq:rem_not_compat_equiequa}
\begin{split}
    (I - \diag(\Tilde{p}^1))(D^1)^{-1} B^1 \Tilde{p}^1 &= \Tilde{p}^1, \\
    (I - \diag(\Tilde{p}^2))(D^2)^{-1} B^2 \Tilde{p}^2 &= \Tilde{p}^2. 
\end{split}
\end{equation}
Since $B^1 = \nu B^2$, $D^1 = \nu D^2$ gives $(D^1)^{-1} B^1 = (D^2)^{-1} B^2$, it 
follows from~\eqref{eq:rem_not_compat_equiequa} that $\Tilde{p}^1 = \Tilde{p}^2 = \Tilde{p}$, and therefore
\begin{equation} \label{eq:rem_not_compat_equiequa_equal}
\begin{split}
    (I - \diag(\Tilde{p}))(D^1)^{-1} B^1 \Tilde{p} &= \Tilde{p}, \\
    (I - \diag(\Tilde{p}))(D^2)^{-1} B^2 \Tilde{p} &= \Tilde{p}.  
\end{split}
\end{equation}
By Proposition~\ref{prop:necessity} we have $\Tilde{p} \ll \textbf{1}$, implying that $(I-\diag(\Tilde{p}))$ is a positive diagonal matrix. Then, for $k \in [2]$, it follows that $(I - \diag(\Tilde{p}))(D^k)^{-1} B^k$ is an irreducible nonnegative matrix. Therefore, item~\ref{item:perfrob_pos_necess} in Lemma~\ref{lem:perron_frob} can be applied to~\eqref{eq:rem_not_compat_equiequa_equal}, from which it follows that $\rho ((I - \diag(\Tilde{p}))(D^1)^{-1} B^1) = 1$ and $\rho ((I - \diag(\Tilde{p}))(D^2)^{-1} B^2) = 1$.
Applying Lemma~\ref{lem:eigspec} we obtain
\begin{equation} \label{eq:rem_cond_zero}  
\begin{split}
    s(-D^1 + (I - \diag(\Tilde{p}^2)) B^1) &= 0, \\
    s(-D^2 + (I - \diag(\Tilde{p}^1)) B^2) &= 0. 
\end{split}
\end{equation}
Observe ~\eqref{eq:rem_cond_zero} is incompatible with the following  condition in Corollary~\ref{cor:joint_eq_exist},
\vspace{-3ex}
\begin{equation} \label{eq:rem_cond_pos} 
\begin{split}
    s(-D^1 + (I - \diag(\Tilde{p}^2)) B^1) &> 0, \\
    s(-D^2 + (I - \diag(\Tilde{p}^1)) B^2) &> 0.
\end{split}
\end{equation}
Hence, it follows that the conditions in Corollary~\ref{cor:joint_eq_exist} and \cite[Theorem~6]{pareautomatica} are mutually exclusive.~$\blacksquare$
\end{rem}














Since the conditions in Corollary~\ref{cor:joint_eq_exist} are not compatible
with the conditions in \cite[Theorem~6]{pareautomatica}, the
question: \enquote{do the conditions in Corollary~\ref{cor:joint_eq_exist} guarantee uniqueness of the coexisting equilibrium?} is worth investigating. Our simulations show that such a coexisting equilibrium may indeed be unique, and further, asymptotically stable; see Section~\ref{sect:simulations} and Figure~\ref{fig:coexistence}.

While Theorem~\ref{thm:joint_eq_exist_shared} provides conditions for the existence of coexisting equilibria, a related problem is finding conditions under which no coexisting equilibria can exist. The following theorem makes use of a nontrivial condition to eliminate the possibility of coexisting equilibria in a bi-virus setting, and establishes one virus as being dominant.
\begin{thm} \label{thm:noequi_single-virus}
Consider the SIWS model~\eqref{eq:full} under Assumption~\ref{assum:base} with $m=2$. Suppose that $B_w^1$ and $B_w^2$ are irreducible matrices, and that $s(B_w^1 - D_w^1) > 0$ and $s(B_w^2 - D_w^2) > 0$. 
If $(D_w^1)^{-1}B_w^1 > (D_w^2)^{-1}B_w^2$
then there are exactly three equilibria in $\mathcal{D}$, namely the healthy state, which is unstable, $(\textbf{0}, \Tilde{y}^2)$ with $\textbf{0} \ll \Tilde{y}^2 \ll \textbf{1}$, which is unstable, and $(\Tilde{y}^1, \textbf{0})$ with $\textbf{0} \ll \Tilde{y}^1 \ll \textbf{1}$, which is locally exponentially stable.~$\blacksquare$
\end{thm}
\textit{Proof:} See Appendix.~$\square$




\noindent
Theorem~\ref{thm:noequi_single-virus} states that if one virus has a stronger set of spread and healing parameters than the other virus, then these two viruses cannot coexist in the population.
Theorem~\ref{thm:noequi_single-virus} answers question~\ref{q5} in Section~\ref{subsec:pb:statement}. The following remark aids in understanding the result in Theorem~\ref{thm:noequi_single-virus}.
\begin{rem} \label{rem:no-coexist}
Underlying Theorem~\ref{thm:noequi_single-virus} is the so-called competitive exclusion principle, which states that complete competitors cannot coexist \cite{hardin1960competitive}. For instance, if two strains of viruses (say virus~$1$ and virus~$2$) compete with each other to infect the same population, and if virus~$1$
has a slight advantage over virus~$2$
(for example higher infection rates or lower healing rates), then virus~$1$
will eventually displace virus~$2$,
which will get eradicated
~$\blacksquare$ 
\end{rem}
Similar results to Theorem~\ref{thm:noequi_single-virus} can be found in \cite[Theorem~5]{liu2019analysis}, for the setting without a shared resource. In order to make a comparison, we employ the SIS model~\eqref{eq:fullnosharedsys} to state the following corollary of Theorem~\ref{thm:noequi_single-virus}.
\begin{cor} \label{cor:noequi_single-virus_noshared}
Consider the SIS model~\eqref{eq:fullnosharedsys} under Assumption~\ref{assum:noshared} with $m=2$. Suppose that $B^1$ and $B^2$ are irreducible matrices, and that $s(B^1 - D^1) > 0$ and $s(B^2 - D^2) > 0$. 
If $(D^1)^{-1}B^1 > (D^2)^{-1}B^2$ 
then there are exactly three equilibria in $[0,1]^{2n}$, namely the healthy state, which is unstable, $(\textbf{0}, \Tilde{p}^2)$ with $\textbf{0} \ll \Tilde{p}^2 \ll \textbf{1}$, which is unstable, and $(\Tilde{p}^1, \textbf{0})$ with $\textbf{0} \ll \Tilde{p}^1 \ll \textbf{1}$, which is locally exponentially stable.~$\blacksquare$
\end{cor}
Corollary~\ref{cor:noequi_single-virus_noshared} is directly comparable to \cite[Theorem~5]{liu2019analysis}. The main difference is that \cite[Theorem~5]{liu2019analysis} requires that $\beta_{ij}^k = \beta^k$ or $\beta_{ij}^k = 0$, and $\delta_i^k = \delta^k$, for all $i,j \in [n]$, $k \in [2]$, and some $\beta^k > 0$, $\delta^k > 0$, whereas Corollary~\ref{cor:noequi_single-virus_noshared} has no such restrictions. Hence, Corollary~\ref{cor:noequi_single-virus_noshared} subsumes and improves upon \cite[Theorem~5]{liu2019analysis}, while Theorem~\ref{thm:noequi_single-virus} extends the result to the setting with a shared resource.


\section{Mitigation strategies} \label{sect:dist:control}
In this section, we present two strategies for ensuring that some (or all) viruses are eradicated. First, we
establish that it is possible to cause convergence to the eradicated state of a virus~$k$ by boosting the corresponding healing rates~$\delta_i^k$. Applying this technique to all $m$ viruses, the system converges to the healthy state. Second, we show that it is possible to leverage a benign virus in order to eradicate a malignant virus, in a bi-virus setting.
\subsection{Boosting the healing rates} \label{subsec:boost_healing}
When managing the epidemic spread of a virus, a natural strategy is to boost the healing rates in the population nodes. The following result shows that the healing rates can always be chosen to ensure asymptotic or exponential eradication of a virus, using an algorithm inspired by \cite{liu2019analysis,gracy2020asymptotic}.
\begin{prop} \label{prop:deltas_chosen_to_heal}
Consider the SIWS model~\eqref{eq:yk} under Assumption~\ref{assum:base} with
$y(0) \in~\mathcal{D}$. Suppose, for some~$k \in [m]$, that $B_w^k$ is irreducible,
and that the healing rates are of the form
\begin{equation} \label{eq:deltas_chosen_to_heal}
\delta_i^k = \beta_{iw}^k + \textstyle \sum_{j=1}^n \beta_{ij}^k + \epsilon_i^k,
\end{equation}
for each $i \in [n]$, where $\epsilon_i^k \geq 0$. Then, if $\epsilon_i^k > 0$ for some $i \in [n]$, the eradicated state of virus~$k$ is exponentially stable, with domain of attraction containing $\mathcal{D}^k$. Otherwise,
the eradicated state of virus~$k$ is asymptotically stable, with domain of attraction containing $\mathcal{D}^k$.~$\blacksquare$
\end{prop}
The proof is straightforward, following from Theorem~\ref{thm:asymp} and similar arguments as in \cite[Section~V]{liu2019analysis}, and \cite{gracy2020asymptotic}.\\
\textit{Proof:} See Appendix.~$\square$


Proposition~\ref{prop:deltas_chosen_to_heal} represents one strategy to ensure eradication of a virus. By applying this strategy to all viruses,
we obtain the following result.
\begin{thm} \label{thm:all_deltas_chosen_to_heal}
Consider the SIWS model~\eqref{eq:full} under Assumption~\ref{assum:base}. Suppose, for all~$k \in m$, that $B_w^k$ is irreducible and \eqref{eq:deltas_chosen_to_heal} is fulfilled.
Then, if $\epsilon_i^k > 0$
for some $i \in [n]$ and all $k \in [m]$, the healthy state is exponentially stable, with domain of attraction containing $\mathcal{D}$. If not, the healthy state is asymptotically stable, with domain of attraction containing~$\mathcal{D}$.~$\blacksquare$
\end{thm}
\noindent
Theorem~\ref{thm:all_deltas_chosen_to_heal} represents a strategy to eradicate all viruses in a system, which equivalently ensures exponential (resp. asymptotic) convergence to the healthy state. Thus, Theorem~\ref{thm:all_deltas_chosen_to_heal} addresses question~\ref{q6}
in Section~\ref{subsec:pb:statement}.
 
\begin{rem}[Epidemiological interpretation]\label{rem:boost:healing:rates}
The mitigation strategy outlined in Theorem~\ref{thm:all_deltas_chosen_to_heal} could be understood as follows: with respect to each virus, if the healing rate of each subpopulation is sufficiently increased, which could be accomplished
by prescribing high dosages of drugs, by administering vaccines, etc.,  then each of the viruses get eradicated exponentially (resp. asymptotically) fast. Observe that this 
strategy is extreme in the sense that it does not factor in limitations on the availability of resources, and essentially encourages health administration officials to amass (possibly) excessive amounts of resources in order to prevent epidemic outbreaks.~$\blacksquare$
\end{rem}
Note that in the absence of sufficient resources, implementing the aforementioned strategy in practice is impossible. Hence, we are motivated to seek different strategies.
\subsection{Leveraging one virus to eradicate another}\label{subsect:competive:exclusion}

It turns out that in a bi-virus setting, where one virus is malignant and the other virus is benign, we can leverage the benign virus in order to help eradicate the malignant virus, as stated in the following theorem. 
\begin{thm} \label{thm:virus-as-vaccine}
Consider the SIWS model~\eqref{eq:full} under Assumption~\ref{assum:base} with $m=2$. Suppose that $B_w^1$ and $B_w^2$ are irreducible matrices; 
$s(B_w^1 - D_w^1) > 0$; $s(B_w^{2} - D_w^{2}) > 0$; 
$c^1=c^2$; and $E^2 \subseteq E^1$, where $E^1$ (resp. $E^2$) is the set of directed edges between the population nodes and the shared resource, with respect to virus~$1$ (resp. virus~$2$).
If the healing rates for virus~$2$ fulfill
\vspace{-1ex}
\begin{equation} \label{eq:deltas_chosen_to_leverage}
\delta_i^2 > \max_{j\in[n+1]} \left\{\frac{(B_w^2)_{ij}}{((D_w^1)^{-1}B_w^1)_{ij}} \right\}_{(B_w^1)_{ij}>0},
\end{equation}

\vspace{-1ex}

\noindent for all $i \in [n]$, then the only locally asymptotically stable equilibrium in $\mathcal{D}$ is $(\Tilde{y}^1, \textbf{0})$ with $\textbf{0} \ll \Tilde{y}^1 \ll \textbf{1}$.~$\blacksquare$
\end{thm}
\textit{Proof:} See Appendix.~$\square$

\noindent
Theorem~\ref{thm:virus-as-vaccine} represents a strategy to eradicate one of the viruses in a bi-virus system, made possible by leveraging the fact that one virus has a stronger set of spread parameters than the other. Thus, Theorem~\ref{thm:virus-as-vaccine} addresses question~\ref{q7}
in Section~\ref{subsec:pb:statement}. 

\begin{rem}[Virus as vaccine]\label{rem:virus:vaccine}
Since the strategy given in Theorem~\ref{thm:virus-as-vaccine}
ensures local asymptotic convergence to the single-virus endemic equilibrium of the benign virus,
it 
could also be interpreted in the following sense: the benign virus effectively acts as a \emph{vaccine} against the malignant virus. In the context of 
battling
epidemic outbreaks, where the goal is to minimize the mortality 
rate, this strategy could potentially provide health administration officials with an effective tool.
~$\blacksquare$
\end{rem}

The mitigation strategies detailed in this section can be compared as follows.
On the one hand, assuming that the
objective of public health  officials 
is solely to eradicate one virus 
in a bi-virus system 
while considering
resource constraints (e.g., availability of vaccines, drugs, ventilators, etc.),
it may be more feasible
to implement the strategy given in Theorem~\ref{thm:virus-as-vaccine} instead of that in Proposition~\ref{prop:deltas_chosen_to_heal}. 
On the other hand, 
the strategy in Theorem~\ref{thm:virus-as-vaccine}, 
opposed to
that in Theorem~\ref{thm:all_deltas_chosen_to_heal} particularized for a bi-virus setting, requires the persistence of one virus, 
which may be undesirable. 
Moreover, with respect to a given virus, if the healing rate of at least one node is sufficiently boosted, then Proposition~\ref{prop:deltas_chosen_to_heal} guarantees eradication of said virus exponentially fast, whereas Theorem~\ref{thm:virus-as-vaccine} guarantees only asymptotic eradication of a virus. We 
explore the advantages and disadvantages of the strategies outlined in Proposition~\ref{prop:deltas_chosen_to_heal} and Theorem~\ref{thm:virus-as-vaccine} via simulations in Section~\ref{sect:simulations}

 







\section{Simulations}
\label{sect:simulations}
In this section, we present a number of simulations to illustrate our theoretical findings, using the city of Stockholm as an example setting.
In particular,~15 districts in and around Stockholm are taken to be the population nodes of the network, thus, $n = 15$. All of these districts are connected to the Stockholm metro, which we view as the shared resource of the network; see Figure~\ref{fig:stockholm}.


\begin{figure}[h!] 
\centering
    \frame{{\includegraphics[width=0.6\columnwidth]{Stockholm_with_T.png}}}
    \caption{A map of~15 city districts in Stockholm. All districts are connected to the metro system.}
    \label{fig:stockholm}
\end{figure} 

In all simulated scenarios we consider two competing viruses, namely virus~$1$ and virus~$2$. We denote the average infection ratio of virus~$k$, i.e., $\textstyle \tiny{\frac{1}{n}} \sum_i^n p_i^k(t)$, by $\Bar{p}^k(t)$. The contact network of each virus is the same, i.e., $E^1 = E^2$, and is represented in Figure~\ref{fig:stockholm}. We set $p_i^k(0) = 0.5$ and $z^k(0) = 0.5$ for all $i \in [15]$ and $k \in [2]$, and use the same
in all simulated scenarios throughout this section. 
Observe that the aforementioned choice of initial states represents the case where half of the population in each node is infected by virus~$1$, while the other half is infected by virus~$2$, and the shared resource is contaminated by both viruses. The spread parameters $\beta_{ij}^k$ are taken to be~$1$ if district $i$ is adjacent to
district $j$, or if $i = j$, and $0$ otherwise, for $k \in [2]$.
In all scenarios, we also assume that each district is bi-directionally connected to the shared resource, i.e., the Stockholm metro, with $\beta_{iw}^k = 1$ and $c_i^k = 1/15$, for all $i \in [15]$ and $k \in [2]$. As a consequence, $B_w^k$ is irreducible for $k \in [2]$. The following scenarios differ only in terms of the choice of $\delta_i^k$ and $\delta_w^k$. 



\begin{figure}[h!] 
\centering
    \scalebox{0.8}[0.8]{\includegraphics[width=\columnwidth]{strong_vs_weak.png}}
    \caption{Simulation with two viruses, one (blue) converging to eradication, while the other (red) converges to its single-virus endemic equilibrium.
    }
    \label{fig:strong_vs_weak}
\end{figure}

In the simulation depicted in Figure~\ref{fig:strong_vs_weak}, we chose $\delta_i^1 = 4.6$, $\delta_w^1 = 4$ and $\delta_i^2 = 10$, $\delta_w^2 = 10$, for all $i \in [15]$. As a consequence, $s(B_w^2 - D_w^2) = -4.2$, so Theorem~\ref{thm:expo} applies to virus~$2$. Consistent with the result in Theorem~\ref{thm:expo}, virus~$2$ becomes eradicated exponentially fas
. However, since $s(B_w^1-D_w^1) = 0.3$, virus~$1$ fulfills the conditions for Proposition~\ref{prop:necessity}, and furthermore, once virus~$1$ is eradicated, the system is essentially a single-virus system. Treated as such, virus~$1$ satisfies the conditions in Theorem~\ref{thm:equi}. In line with the result in Theorem~\ref{thm:equi},
virus~$1$ converges to a single-virus endemic equilibriu
. Insofar as it can be determined by varying $y^1(0)$ in $\mathcal{D}^1$, this single-virus endemic equilibrium appears to be unique and asymptotically stable. 

\begin{figure}[h!] 
\centering
    \scalebox{0.8}[0.8]{\includegraphics[width=\columnwidth]{coexistence.png}}
    \caption{Simulation with two viruses (red and blue), converging to a coexisting equilibrium.
    }
    \label{fig:coexistence}
\end{figure}

In the simulation depicted in Figure~\ref{fig:coexistence}, we chose $\delta_i^1 = 1.5$ for $i \in [8]$, and $\delta_i^1 = 2$ for 
$i \in [15]\setminus [8]$, with $\delta_w^1 = 1$. Relating these choices of parameters for virus~$1$ to Figure~\ref{fig:stockholm}, this corresponds to a low healing rate north of the river Mälaren and a high healing rate south of the river.
Mirroring this pattern, we chose $\delta_i^2 = 2$ for $i \in [8]$, and $\delta_i^2 = 1.5$ for 
$i \in [15]\setminus [8]$, with $\delta_w^2 = 1$. With these parameters, it follows that $s(B_w^1-D_w^1) = 2.8$, and $s(B_w^2 - D_w^2) = 2.8$. Hence, both viruses fulfill the conditions for Proposition~\ref{prop:necessity}, providing the existence of exactly two single-virus endemic equilibria, namely $(\Tilde{y}^1, \textbf{0})$ and $(\textbf{0}, \Tilde{y}^2)$. $\Tilde{y}^1$ can be approximated by setting $y^1(0) > \textbf{0}$ and $y^2(0) = \textbf{0}$, and running the simulation for a sufficiently long period of time $T$. Then, assuming that $\Tilde{y}^1 \approx y^1(T)$, and, with an analogous approximation for virus~$2$, $\Tilde{y}^2 \approx y^2(T)$,
we obtain $s((I-X(\Tilde{y}^2))B_w^1-D_w^1) = 0.2$, and $s((I-X(\Tilde{y}^1))B_w^2 - D_w^2) = 0.2$. As a consequence, this pair of viruses fulfills the conditions for Theorem~\ref{thm:joint_eq_exist_shared}. In line with the result in Theorem~\ref{thm:joint_eq_exist_shared}, we see that there exists a coexisting equilibrium. Moreover, our simulations show that the viral infection levels appear to converge to this coexisting equilibriu
. Additionally, regardless of how the initial condition is varied within $\mathcal{D}$, barring $y^1(0) = \textbf{0}$ or $y^2(0) = \textbf{0}$, we observe that all simulations converge to the same coexisting equilibrium. This observation suggests that the coexisting equilibrium might be unique, as well as asymptotically stable.

\begin{figure}[h!] 
\centering
    \scalebox{0.8}[0.8]{\includegraphics[width=\columnwidth]{strong_vs_stronger.png}}
    \caption{Simulation with two viruses, with one (blue) converging to eradication, while the other (red) eventually converges to its single-virus endemic equilibrium.
    }
    \label{fig:strong_vs_stronger}
\end{figure}

In the simulation depicted in Figure~\ref{fig:strong_vs_stronger}, we chose $\delta_i^1 = 3$, $\delta_w^1 = 3$ and $\delta_i^2 = 4$, $\delta_w^2 = 4$, for all $i \in [15]$. It follows that $s(B_w^1-D_w^1) = 1.7$, and $s(B_w^2 - D_w^2) = 0.9$. Moreover, we have $(D_w^1)^{-1}B_w^1 > (D_w^2)^{-1}B_w^2$. Therefore, this pair of viruses fulfills the conditions for Theorem~\ref{thm:noequi_single-virus}. In line with the result in Theorem~\ref{thm:noequi_single-virus}, we observe that
virus~$1$ has a competitive edge, allowing it to persist and converge to its single-virus endemic equilibrium. Meanwhile, virus~$2$ is eradicated, despite $s(B_w^2 - D_w^2)$ being positive.


\begin{figure}[h!] 
\centering
    \scalebox{0.8}[0.8]{\includegraphics[width=\columnwidth]{coexistence_heal.png}}
    \caption{Simulation with two viruses (red and blue), initialized at a coexisting equilibrium. At $t = 0.5$ the healing rates of virus~$2$ are changed to fulfill Proposition~\ref{prop:deltas_chosen_to_heal}.
    }
    \label{fig:control_heal}
\end{figure}

\begin{figure}[h!] 
\centering
    \scalebox{0.8}[0.8]{\includegraphics[width=\columnwidth]{coexistence_exclusion.png}}
    \caption{Simulation with two viruses (red and blue), initialized at a coexisting equilibrium. At $t = 0.5$ the healing rates of virus~$2$ are changed to fulfill Theorem~\ref{thm:virus-as-vaccine}.
    }
    \label{fig:control_exclusion}
\end{figure}

The simulations depicted in Figures~\ref{fig:control_heal} and~\ref{fig:control_exclusion} were both initialized at the coexisting equilibrium from Figure~\ref{fig:coexistence}, with all parameters as in that simulation.
At $t = 0.5$ in the simulation in Figure~\ref{fig:control_heal},
the healing rates of virus~2, i.e., $\delta_i^2$,
are chosen as in~\eqref{eq:deltas_chosen_to_heal}, 
with $\epsilon_i^2 = 0$ for all $i \in [15]$. Similarly, at $t = 0.5$ in the simulation in Figure~\ref{fig:control_exclusion}, the healing rates of virus~2, i.e., $\delta_i^2$, are chosen as in~\eqref{eq:deltas_chosen_to_leverage}.  
More specifically, $\delta_i^2 = 2$ for $i \in [8]$ as before, but $\delta_i^2 = 2.1$ for 
$i \in [15]\setminus [8]$.
We assume that the cost of a strategy $\ell$ (denoted as $\cost^{\ell}$) is given by the sum of all healing rates (i.e., $\cost^{\ell} = \sum_{i=1}^{15} \delta_i^2$, where $\delta_i^2$ is chosen as in~\eqref{eq:deltas_chosen_to_heal} 
when $\ell =1$, and as in~\eqref{eq:deltas_chosen_to_leverage} when $\ell =2$).
Our simulations show that $\cost^{1} = 61$,
whereas $\cost^{2} = 33$.
However, as seen in Figures~\ref{fig:control_heal} and~\ref{fig:control_exclusion}, the end result of both strategies is the same. That is, virus~1 persists and reaches its single-virus endemic equilibrium, whereas virus~2 is eradicated. Hence, the strategy in Theorem~\ref{thm:virus-as-vaccine} allows us to eradicate virus~2 at a lower cost than the strategy in Proposition~\ref{prop:deltas_chosen_to_heal}, in this case.
However, the convergence to the single-virus endemic equilibrium of virus~1 is faster with Proposition~\ref{prop:deltas_chosen_to_heal}, as expected from the excessive healing rates.


\section{Conclusions}
\label{sect:conclusions}

In this paper, we introduced and analyzed a novel SIWS model of competitive, multi-viral spread across a network of population nodes with a shared resource. We established conditions under which a virus is eradicated exponentially fast (resp. asymptotically). We also provided a condition that ensures a virus can reach and sustain a unique, single-virus endemic equilibrium of infection levels across the network. Taken together, these results allowed us to state a necessary and sufficient condition for convergence to the healthy state.
Moreover, we presented sufficient conditions under which two viruses can sustain a coexisting equilibrium, 
that is, where
neither virus is eradicated.
Conversely, we also provided a necessary condition for the existence of such a coexisting equilibrium. To mitigate the spreading process, we proposed two strategies.
The first strategy involved choosing the healing rate of each subpopulation, with respect to each virus, in a suitable manner so as to ensure that all subpopulations converge to the healthy state.
The second strategy exploited the notion of competitive exclusion to ensure that in a bi-virus setup, with one virus being malignant and the other benign, the malignant virus becomes eradicated.



The present paper concerns time-invariant SIWS models, whereas in order to obtain a better understanding of real-world scenarios (viz. mobile agents, mutating viruses), it is natural to consider time-varying SIWS models.
Secondly, in order to generalize the result on the existence of a coexisting equilibrium (cf. Theorem~\ref{thm:joint_eq_exist_shared}) to more than two viruses, we would need to establish uniqueness of the coexisting equilibrium for the bi-virus case, which is non-trivial.

\vspace{-1ex}

\bibliographystyle{IEEEtran}
