# Convergence of Bi-Virus Epidemic Models with Non-Linear Rates on Networks -
A Monotone Dynamical Systems Approach

Vishwaraj Doshi, Shailaja Mallick, and Do Young Eun
Accepted for publication at IEEE/ACM Transactions on Networking, in September
2022. A subset of the material in this paper appears in [1]. Vishwaraj Doshi
is with the Operations Research Graduate Program, Shailaja Mallick is with the
Department of Computer Science, and Do Young Eun is with the Department of
Electrical and Computer Engineering, North Carolina State University, Raleigh,
NC. Email: {vdoshi, smallic<EMAIL_ADDRESS>This work was supported in part
by National Science Foundation under Grant Nos. CNS-2007423 and CNS-1824518.

###### Abstract

We study convergence properties of competing epidemic models of the
Susceptible-Infected-Susceptible ($SIS$) type. The SIS epidemic model has seen
widespread popularity in modelling the spreading dynamics of contagions such
as viruses, infectious diseases, or even rumors/opinions over contact networks
(graphs). We analyze the case of two such viruses spreading on overlaid
graphs, with non-linear rates of infection spread and recovery. We call this
the non-linear bi-virus model and, building upon recent results, obtain
precise conditions for global convergence of the solutions to a trichotomy of
possible outcomes: a virus-free state, a single-virus state, and to a
coexistence state. Our techniques are based on the theory of monotone
dynamical systems (MDS), in contrast to Lyapunov based techniques that have
only seen partial success in determining convergence properties in the setting
of competing epidemics. We demonstrate how the existing works have been
unsuccessful in characterizing a large subset of the model parameter space for
bi-virus epidemics, including all scenarios leading to coexistence of the
epidemics. To the best of our knowledge, our results are the first in
providing complete convergence analysis for the bi-virus system with non-
linear infection and recovery rates on general graphs.

###### Index Terms:

Epidemics on networks, bi-virus models, multi-layer graphs, monotone dynamical
systems.

##  I Introduction and overview

Graph-based epidemic models are widely employed to analyze the spread of real
world phenomena such as communicable diseases [2, 3], computer viruses,
malware [4, 5, 6], product adoption [7, 8, 9], opinions, and rumors [10, 11,
12, 13]. The propagation of such phenomenon (which we cumulatively refer to as
epidemics or viruses) usually takes place via processes such as human contact,
word-of-mouth, exchange of emails or even in social media platforms. Graph
based techniques, with edge based mechanisms to model information spread, have
therefore proven to be effective in capturing such epidemic dynamics, and have
been a research focus over the past few decades [14, 15, 16, 17]. In recent
years, the development of models which capture the competition of two or more
of such epidemics has seen a surge of interest. In particular, models
capturing the behavior of two competing epidemics of the Susceptible-Infected-
Susceptible (SIS) types, also known as the bi-virus or bi-SIS models, have
garnered significant attention over the years [8, 18, 19, 20, 21].

Epidemic models take the form of ordinary differential equations (ODEs) and
their analysis involves the identification of fixed points of the system,
their uniqueness properties, and ultimately showing the convergence of the
solution trajectories to those fixed points. The technique via Lyapunov
functions has historically been a popular method to prove convergence to fixed
points and was also used in epidemiology literature to derive the convergence
properties of the SIS epidemic model. The SIS model was originally introduced
in [2] to capture the spread of Gonorrhea due to contact between individuals
in a population, and was further developed in [22, 23, 24, 25, 26, 27, 28,
29]. The central result for SIS epidemics, originally proved using Lyapunov
functions in [2], is a dichotomy arising from the relation between model
parameter ($\tau\\!>\\!0$) representing the effective infection rate or
strength of the virus,111$\tau=\beta/\delta$, where $\beta>0$ stands for the
infection rate of the virus and $\delta>0$ the recovery rate from the virus.
Section II provides a detailed explanation. and a threshold value
($\tau^{*}\\!>\\!0$). When $\tau\\!\leq\\!\tau^{*}$, the virus spread is not
strong enough and the system converges to a ‘virus-free’ state. When
$\tau\\!>\\!\tau^{*}$, it converges to a state where the virus infects a non-
zero portion of the population. Attempts have also been made to perform
similar convergence analysis for the bi-virus epidemic model [8, 19, 20, 21].
The key questions posed in such literature are: Can both competing epidemics
coexist over the network? If not, which one prevails? Or do both die out? This
trichotomy of possible results is what the recent literature has been trying
to characterize.

When the propagation of the two epidemics occurs over the same network [8,
30], it has been established that coexistence of two viruses is impossible
except in the rare cases where their effective strengths
($\tau_{1},\tau_{2}\\!>\\!0$ for viruses 1, 2, respectively) are equal [21, 8,
20, 19, 18]; the virus with the larger effective strength otherwise wiping out
the other, a phenomenon sometimes referred to as winner takes all [8]. The
situation is much more complicated when the two viruses spread over two
distinct networks overlaid on the same set of nodes. This modeling approach is
more representative of the real world, where competing rumors/products/memes
may not use the same platforms to propagate, though they target the same
individuals. Recent works [18, 21, 19, 20, 31, 32, 33, 34] therefore consider
this more general setting, but unfortunately, a complete characterization of
the trichotomy of outcomes has still proven to be elusive and remains open as
of now.

While the original SIS model introduced in [2] had the aggregate infection and
recovery rates of a node as linear functions of the number of infected
neighbors, there has been a push towards studying more generalized models
where these rates are made heterogeneous (across nodes) and _non-linear_ [35,
36, 37, 38, 39]. Realistic assumptions such as infection rates tending to
saturation with continual increase in neighborhood infection [40, 41, 42, 43]
have become more commonplace, implying that the models employing strictly
linear spreading dynamics often provide overestimates to the real world
infection rates [24, 20]. This paper does not concern itself with answering
which non-linear infection rate best captures the exact dynamics, but we
direct the readers to [20] which provides simulation results comparing non-
linear rate functions to the exact Markovian dynamics for some special
randomly generated graph topologies. In some special cases, non-linear
recovery rates also have an interpretation linking them to reliability theory
in the form infection duration with increasing failure rates (failure here
being the recovery of an infected node). Allowing for non-linear infection and
recovery rates leads to a more general version of the bi-virus model on
overlaid graphs, albeit much more complicated, and the complete convergence
criterion is yet to be fully established [19, 20]. It should be noted that
while we extensively refer to the infection and recovery rates being either
linear or non-linear in this paper, the bi-virus epidemic model itself will
always be a system of non-linear ODEs.

#### Limitations of existing works

Of all the recent works concerning the spread of SIS type bi-virus epidemics
on overlaid networks, [20] and [19] provide conditions under which the system
globally converges to the state where one virus survives while the other dies
out. [20] approaches the problem of showing global convergence by employing
the classic technique via Lyapunov functions. However, finding appropriate
Lyapunov functions is a highly non-trivial task, and as mentioned in [19], is
even more difficult due to the coupled nature of the bi-virus ODE system. This
can be seen in the condition they derive in [20] for the case where, say,
Virus 1 dies out and Virus 2 survives. When $\tau_{1}$ and $\tau_{2}$
represent the effective strengths of Virus 1 and Virus 2, respectively, their
condition translates to $\tau_{1}\\!\leq\\!\tau_{1}^{*}$ where $\tau_{1}^{*}$
is the threshold corresponding to the single-virus case, meaning that Virus 1
would not have survived even if it was the only epidemic present on the
network. More importantly, [20] is unable to characterize convergence
properties for $\tau_{1}\\!>\\!\tau_{1}^{*}$ and
$\tau_{2}\\!>\\!\tau_{2}^{*}$.

The authors in [19] take a different approach and tackle this problem by
applying their ‘qualitative analysis’ technique, which uses results from other
dynamical systems that bound the solutions of the bi-virus ODE; and provide
conditions under which the system globally converges to single-virus
equilibria. As we show later in Section V-B, however, their conditions not
only characterize just a subset of the actual space of parameters that lead to
global convergence to the single-virus equilibria (which they themselves
pointed out), but the size of this subset is highly sensitive to the graph
topology, often much smaller than what it should be in general. In other
words, a complete characterization of the _entire_ space of model parameters,
on which the system globally converges to one of the trichotomic states, has
still been recognized as an open problem in the bi-virus literature [20, 19,
21].

#### Our contributions

In this paper, we analyze the bi-virus model with _non-linear_ infection and
recovery rates (or the _non-linear bi-virus model_ in short) and provide the
complete characterization of the trichotomy of the outcomes with necessary and
sufficient conditions under which the system globally converges to one of the
three possible points: (i) a ‘virus-free’ state, (ii) a ‘single-virus’
equilibrium, or (iii) an equilibrium where both viruses coexist over the
network. While the result for convergence to the virus-free state of the bi-
SIS model is not new for non-linear infection and linear recovery rates, our
proof for the same is the most general form known to date, covering the case
with both infection _and_ recovery rates being non-linear. The proof of
convergence to the virus-free state of the bi-virus model is straightforward,
and directly follows from the convergence criterion for the single-virus SIS
model with non-linear rates. However, the convergence results for fixed points
where only one of the two viruses survives, or to the equilibrium where both
viruses coexist, are not as straightforward to establish, rendering the
typical Lyapunov based approach largely inapplicable.

In proving these results, we first show, using a specially constructed cone
based partial ordering, that the bi-virus epidemic model possesses some
inherent monotonicity properties. We then use novel techniques from the theory
of _monotone dynamical systems_ (MDS) [44] to prove our main results. In
recent control systems literature [45, 46, 47, 48, 49], techniques based on
the construction of cone based partial orderings that leverage the
monotonicity properties of dynamical systems have indeed been studied.
Dynamical systems exhibiting such monotonicity properties are also sometimes
called deferentially positive systems [50] and cooperative systems [51] in the
ODE setting, with interesting applications in consensus problems for
distributed systems [52] and even neural networks [53]. In this paper, we
utilize these MDS techniques in the setting of competing epidemics, and as a
result demonstrate an alternative to Lyapunov based approaches to analyze
convergence properties of epidemic models. The novelty of using the MDS
approach for analysis also lies with [54], which uses similar techniques to
analyze the bi-virus system for the special case of linear infection and
recovery rates, and was developed concurrently and independently with the
initial version of this work [1]. This further highlights the utility of MDS
techniques for the analysis of epidemic models on graphs.

This paper is an extension of our previous work [1], which gives necessary and
sufficient conditions for convergence to the three types of equilibria only
for the special case of the bi-virus model with _linear_ infection and
recovery rates (or the _linear bi-virus model_ in short). Our conditions
therein take a more precise form in terms of the model parameters $\tau_{1}$
and $\tau_{2}$ and one can visualize an exact partition of the model parameter
space into regions corresponding to various convergence outcomes. We note that
this partition of the model parameter space coincides with that in [18],
wherein they employed only _local_ stability results via bifurcation analysis
– concerning only solution trajectories that originate from a small
neighborhood of those fixed points. In contrast, our results in this paper
concern global stability of the system with any combination of linear as well
as more general, non-linear infection and recovery rates.

#### Structure of the paper

In Section II, we first introduce the basic notation used throughout the
paper, along with the classical (single-virus) SIS model and the bi-virus
model. We then provide the generalization to non-linear infection and recovery
rates in Section III with some key assumptions on the infection and recovery
rate functions, complimented by a discussion in Appendix C regarding a special
class of recovery rates. In Section IV, we provide a primer to the MDS theory,
and establish monotonicity results for the single-virus SIS model, proving the
convergence result for the single-virus model with non-linear infection and
recovery rates whose proofs are deferred to Appendix E. We then go on to show
in Section V-A that the non-linear bi-virus model is also a monotone dynamical
system with respect to a specially constructed cone-based partial ordering,
and include the main convergence results in Section V-B. In Section VI we take
the opportunity to provide a more intuitive version of our results by
considering the special case of linear infection and recovery rates, along
with brief comparisons with the existing literature. In Section VII, we
provide numerical results which confirm our theoretical findings. We then
conclude in Section VIII.

For better readability of the paper, all technical proofs of the main results
are deferred to Appendix F. The appendices also include some selected
definitions and results from matrix theory (Appendix A), ODE theory (Appendix
B), and from MDS theory (Appendix D), which we use as part of our proofs of
the Theorems in Section V-B.

##  II Preliminaries

###  II-A Basic Notations

We standardize the notations of vectors and matrices by using lower case,
bold-faced letters to denote vectors ($\mathbf{v}\\!\in\\!\mathbb{R}^{N}$),
and upper case, bold-faced letters to denote matrices
($\mathbf{M}\\!\in\\!\mathbb{R}^{N\times N}$). We denote by
$\lambda(\mathbf{M})$ the largest real part222We use the $\lambda$ notation
instead of something like $\lambda_{Re}$, since it will mostly be used in
cases where the largest eigenvalue is real, for which $\lambda$ itself is the
largest real eigenvalue. For example, $\lambda(\mathbf{A})$ becomes the
spectral radius for any non-negative matrix $\mathbf{A}$ [55, 56]. of all
eigenvalues of a square matrix $\mathbf{M}$. We use $\text{diag}(\mathbf{v})$
or $\mathbf{D}_{\mathbf{v}}$ to denote the $N\\!\\!\times\\!\\!N$ diagonal
matrix with entries of vector $\mathbf{v}\in\mathbb{R}^{N}$ on the diagonal.
Also, we denote $\mathbf{1}\\!\triangleq\\![1,\\!\cdots\\!,1]^{T}$ and
$\mathbf{0}\\!\triangleq\\![0,\\!\cdots\\!,0]^{T}$, the $N$-dimensional vector
of all ones and zeros, respectively. For vectors, we write
$\mathbf{x}\\!\leq\\!\mathbf{y}$ to indicate that $x_{i}\\!\leq\\!y_{i}$ for
all $i$; $\mathbf{x}\\!<\\!\mathbf{y}$ if $\mathbf{x}\\!\leq\\!\mathbf{y}$ and
$\mathbf{x}\\!\neq\\!\mathbf{y}$; $\mathbf{x}\\!\ll\\!\mathbf{y}$ when all
entries satisfy $x_{i}\\!<\\!y_{i}$. We use
$\mathcal{G}(\mathcal{N},\mathcal{E})$ to represent a general, undirected,
connected graph with $\mathcal{N}\triangleq\\{1,2,\cdots,N\\}$ being the set
of nodes and $\mathcal{E}$ being the set of edges. When we refer to a matrix
$\mathbf{A}\\!=\\![a_{ij}]$ as the adjacency matrix of some graph
$\mathcal{G}(\mathcal{N},\mathcal{E})$, it satisfies
$a_{ij}\triangleq\mathds{1}_{\\{(i,j)\in\mathcal{E}\\}}$ for any
$i,j\in\mathcal{N}$; we use $d_{min}(\mathbf{A})$ and $d_{max}(\mathbf{A})$ to
denote the minimum and maximum degrees of the nodes of the corresponding
graph. Since we only consider connected graphs, all the adjacency matrices in
this paper are automatically considered to be irreducible (see Definition A.1
in Appendix A).

###  II-B $SIS$ Model with Linear rates

Consider the graph $\mathcal{G}(\mathcal{N},\mathcal{E})$, and assume that at
any given time $t\geq 0$, each node $i\in\mathcal{N}$ of the graph is either
in an _infected (I)_ , or in a _susceptible (S)_ state. An infected node can
infect each of its susceptible neighbors with rate $\beta>0$.333We say an
event occurs with some _rate_ $\alpha>0$ if it occurs after a random amount of
time, exponentially distributed with parameter $\alpha>0$. It can also, with
rate $\delta>0$, be cured from its infection and revert to being susceptible
again. We write $\mathbf{x}(t)=[x_{i}(t)]\in\mathbb{R}^{N}$, where $x_{i}(t)$
represents the probability that node $i\in\mathcal{N}$ is infected at any
given time $t\geq 0$. Then, the dynamics of the $SIS$ model can be captured
via the system of ODEs given by

$\frac{dx_{i}(t)}{dt}\triangleq\beta(1-x_{i}(t))\sum_{j\in\mathcal{N}}a_{ij}x_{j}(t)-\delta
x_{i}(t)$ (1)

for all $i\in\mathcal{N}$ and $t\geq 0$. In a matrix-vector form, this can be
written as

$\frac{d\mathbf{x}}{dt}\triangleq\beta\text{diag}(\mathbf{1}-\mathbf{x})\mathbf{A}\mathbf{x}-\delta\mathbf{x}$
(2)

where we suppress the $(t)$ notation for brevity. The system (2) is positively
invariant in the set $[0,1]^{N}$, and has $\mathbf{0}$ as a fixed point (the
virus-free equilibrium). The following result is well known from [2], which we
will generalize in Section IV-B.

######  Theorem II.1 (Theorem 3.1 in [2])

Let $\tau\\!\triangleq\\!\beta/\delta$. Then,

  1. (i)

either $\tau\leq 1/\lambda(\mathbf{A})$, and $\mathbf{x}^{*}=\mathbf{0}$ is a
globally asymptotically stable fixed point of (2);

  2. (ii)

or $\tau>1/\lambda(\mathbf{A})$, and there exists a unique, strictly positive
fixed point $\mathbf{x}^{*}\in(0,1)^{N}$ such that $\mathbf{x}^{*}$ is
globally asymptotically stable in
$[0,1]^{N}\setminus\\{\mathbf{0}\\}$.$\hfill\square$

###  II-C Bi-Virus Model with Linear rates

Consider two graphs $\mathcal{G}_{1}(\mathcal{N},\mathcal{E}_{1})$ and
$\mathcal{G}_{2}(\mathcal{N},\mathcal{E}_{2})$, on the same set of nodes
$\mathcal{N}$ but with different edge sets $\mathcal{E}_{1}$ and
$\mathcal{E}_{2}$. At any given time $t\geq 0$, a node $i\in\mathcal{N}$ is
either infected by Virus 1, infected by Virus 2, or is susceptible. A node
infected by Virus 1 infects each of its susceptible neighbors with rate
$\beta_{1}>0$, just like in the $SIS$ model, but does so only to nodes which
are its neighbors with respect to the graph
$\mathcal{G}_{1}(\mathcal{N},\mathcal{E}_{1})$. Nodes infected by Virus 1 also
recover with rate $\delta_{1}>0$, after which they enter the susceptible
state. Similarly, nodes infected by Virus 2 infect their susceptible
neighbors, this time with respect to the graph
$\mathcal{G}_{2}(\mathcal{N},\mathcal{E}_{2})$, with rate $\beta_{2}>0$, while
recovering with rate $\delta_{2}>0$. This competing bi-virus model of epidemic
spread, also referred to as the $SI_{1}I_{2}S$ model, can be represented by
the following ODE system:

$\begin{split}\frac{dx_{i}}{dt}&\triangleq\beta_{1}\left(1-x_{i}-y_{i}\right)\sum_{j\in\mathcal{N}}a_{ij}x_{j}-\delta_{1}x_{i}\\\
\frac{dy_{i}}{dt}&\triangleq\beta_{2}\left(1-x_{i}-y_{i}\right)\sum_{j\in\mathcal{N}}b_{ij}y_{j}-\delta_{2}y_{i}\end{split}$
(3)

for all $i\in\mathcal{N}$ and $t\geq 0$. In matrix-vector form, (3) becomes:

$\begin{split}\frac{d\mathbf{x}}{dt}&\triangleq\beta_{1}\text{diag}\left(\mathbf{1}-\mathbf{x}-\mathbf{y}\right)\mathbf{A}\mathbf{x}-\delta_{1}\mathbf{x}\\\
\frac{d\mathbf{y}}{dt}&\triangleq\beta_{2}\text{diag}\left(\mathbf{1}-\mathbf{x}-\mathbf{y}\right)\mathbf{B}\mathbf{y}-\delta_{2}\mathbf{y},\end{split}$
(4)

where $\mathbf{A}=[a_{ij}]$ and $\mathbf{B}=[b_{ij}]$ are the adjacency
matrices of graphs $\mathcal{G}_{1}(\mathcal{N},\mathcal{E}_{1})$ and
$\mathcal{G}_{2}(\mathcal{N},\mathcal{E}_{2})$, respectively.

##  III Epidemic Models with Non-linear Infection and Recovery rates

In this section, we introduce the single-virus and bi-virus SIS models with
non-linear infection and recovery rates. Non-linearities can be attributed to
the spread and recovery from the virus being related to the susceptibility of
the disease (or its prevalence in the population) in a more complicated
manner. This is more general than simply exponential random variables with
constant rates used to model the spreading and recovery processes, which in
aggregate scale linearly with the infection probabilities.444‘Aggregate’ here
refers to the mean field approximation which is one way to derive SIS-type
ODEs. Another way is the large population mean field limit of a stochastic
process, where the connection to the corresponding ODE system is formed via
the Kurtz’s theorem [16]. In this case, linearity is induced by the uniform or
homogeneous mixing assumption which is also a subject of criticism in
epidemiology literature [35, 36, 37, 38]. This is shown to be limiting in
accurately modelling the trajectories of an infection spread; the linear
scaling of the infection and recovery rates shown to being an overestimate to
what is observed in reality [20, 37]. Many works thus argue for the modelling
of these spreading processes with non-linear functions [38, 35, 36, 40]. We
first present the more general single-virus SIS model with a set of intuitive
assumptions (A1)–(A5) for the non-linear infection and recovery rates.

###  III-A $SIS$ Model with Non-linear rates

In (1) the term $\sum_{j\in\mathcal{N}}a_{ij}x_{j}(t)$ denotes the overall
rate at which a susceptible node $i\in\mathcal{N}$ gets infected by its
neighbors. In what follows, we replace this by a generic function
$f_{i}(\mathbf{x}(t))$, thereby allowing the overall infection rate for each
node to be any non-linear function of $x_{j}(t)$ for all neighbors $j$ of $i$.
Similarly, we replace the term $\delta x_{i}(t)$, denoting the overall
recovery rate for any node $i\in\mathcal{N}$, by a non-linear function
$q_{i}(\mathbf{x}(t))$. This generic version of the SIS model, allowing for
non-linear infection and recovery rates, is given by the ODE

$\frac{dx_{i}(t)}{dt}=\bar{f}_{i}(\mathbf{x}(t))\triangleq(1-x_{i}(t))f_{i}(\mathbf{x}(t))-q_{i}(\mathbf{x}(t))$
(5)

for all $i\in\mathcal{N}$ and $t\geq 0$. In a matrix-vector form, this can be
written as

$\frac{d\mathbf{x}}{dt}=\bar{F}(\mathbf{x})\triangleq\text{diag}(\mathbf{1}-\mathbf{x})F(\mathbf{x})-Q(\mathbf{x})$
(6)

where $F(\mathbf{x})=[f_{i}(\mathbf{x})]\in\mathbb{R}^{N}$, and
$Q(\mathbf{x})=[q_{i}(\mathbf{x})]\in\mathbb{R}^{N}$ are the vectors of non-
linear infection and recovery rate functions, respectively. We assume that
they are continuous and twice differentiable in $[0,1]^{N}$, with
$\mathbf{J}_{F}(\mathbf{x})$ and $\mathbf{J}_{Q}(\mathbf{x})$ denoting the
Jacobians of $F$ and $Q$ respectively, evaluated at any point
$\mathbf{x}\in[0,1]^{N}$. We now make the following key assumptions:

  * (A1)

$F(\mathbf{0})=\mathbf{0}$ and $Q(\mathbf{0})=\mathbf{0}$;

  * (A2)

$\left[\mathbf{J}_{F}(\mathbf{x})\right]_{ij}=\frac{\partial
f_{i}(\mathbf{x})}{\partial x_{j}}>0~{}~{}\forall i\neq j$ with $a_{ij}>0$,
otherwise $\left[\mathbf{J}_{F}(\mathbf{x})\right]_{ij}=0$;

  * (A3)

$\left[\mathbf{J}_{Q}(\mathbf{x})\right]_{ii}\\!=\\!\frac{\partial
q_{i}(\mathbf{x})}{\partial x_{i}}\\!>\\!0$, and
$\left[\mathbf{J}_{Q}(\mathbf{x})\right]_{ij}\\!=\\!\frac{\partial
q_{i}(\mathbf{x})}{\partial x_{j}}\\!\leq\\!0$ for all $i\\!\neq\\!j$,
$\mathbf{x}\in[0,1]^{N}$. Moreover, $\sum\limits_{j\neq
i}\left[\mathbf{J}_{Q}(\mathbf{x})\right]_{ij}\\!<\\!\left[\mathbf{J}_{Q}(\mathbf{x})\right]_{ii}$;

  * (A4)

$f_{i}(\mathbf{x})$ is concave in $[0,1]^{N}\\!\\!$, that is,
$\frac{\partial^{2}f_{i}}{\partial x_{j}\partial x_{k}}\\!\leq\\!0$ for all
$i,\\!j,\\!k\\!\in\\!\mathcal{N}$;

  * (A5)

$q_{i}(\mathbf{x})$ is convex function of $x_{i}\in[0,1]^{N}$, and a concave
function of $x_{j}$ for all $j\neq i$. That is,
$\frac{\partial^{2}q_{i}}{\partial^{2}x_{i}}\geq 0$ and
$\frac{\partial^{2}q_{i}}{\partial x_{j}\partial x_{k}}\leq 0$ for all
$i\\!\in\\!\mathcal{N}$, and $j,k\\!\in\\!\mathcal{N}\\!\setminus\\!\\{i\\}$.

Assumption (A1) ensures that the virus-free state is a fixed point of (6),
while (A2) is a proximity assumption that models infection spread only through
edges of the underlying graph. Assumption (A3) concerns with the recovery
rate, allowing it to be reduced by infected neighbors while still being no-
negative. (A4) and (A5) assume concavity properties of the functions
$f_{i}(\mathbf{x})$ and $q_{i}(\mathbf{x})$ in $x_{j}$ for any neighbor $j$ of
$i$. This allows the effect of neighborhood infection $x_{j}$ to saturate555As
$x_{j}$ increases for any neighbor $j$ of node $i$, the magnitude of the
resulting change in both infection rate $f_{i}(\mathbf{x})$ and recovery rate
$q_{i}(\mathbf{x})$ decreases. This is similar to the case of diminishing
returns. as $x_{j}$ increases. Assumption (A5) also assumes convexity of
$q_{i}(\mathbf{x})$ in local infection $x_{i}$, which means that increase in
recovery rate caused by $x_{i}$ can be larger as $x_{i}$ increases.

Examples for non-linear infection rates satisfying (A1)–(A5) include
logarithmic functions $f_{i}(\mathbf{x})=\sum_{j}a_{ij}\ln{(1+x_{j})}$,
similar to those in [20]. Examples of non-linear recovery rates include
polynomial functions such as $q_{i}(\mathbf{x})=(1+x_{i})^{k}-1$ for any
$k\geq 1$. A special class of the permissible non-linear recovery rates, where
the infection duration is dependent solely on local infection $x_{i}$, is
related to processes that have decreasing failure rates (DFR)666Failure rate
for a non-negative random variable is defined as the ratio between its
probability density function (PDF) and its complimentary cumulative
distribution function (CCDF). In the context of infection duration, decreasing
failure rate means that nodes recover at a decreased rate the longer they stay
continuously infected. A more detailed discussion regarding the connection to
SIS recovery rates can be found in Appendix C.. This special class of recovery
processes that are DFR also includes the case of linear recovery rates. Note
that our assumptions allow $f_{i}(\mathbf{x})$ and $q_{i}(\mathbf{x})$ to be
heterogeneous across all nodes $i\in\mathcal{N}$, and the case with linear
rates in (2) readily satisfies (A1)–(A5). This also extends to the linear bi-
virus model (4) being a special case of the non-linear bi-virus model
introduced in the next subsection, with infection and recovery rate functions
therein satisfying the same assumptions (A1)–(A5).

Figure 1: Bi-Virus epidemic spread across overlaid graphs sharing the same set
of nodes. Red and Blue arrows denote the spread of Virus 1 and 2, respectively
from infected nodes $j$ and $k$ (coloured Red and Blue) to the susceptible
node $i$ (uncoloured) with the instantaneous rates as shown. The infected Red
and Blue nodes also recover with a total rate of $r_{i}(\mathbf{x})$ and
$s_{i}(\mathbf{y})$ for any node $i\in\mathcal{N}$, respectively.

###  III-B Bi-Virus Model with Non-linear rates

The Bi-Virus model with non-linear infection and recovery rates is given by
the following coupled system of ODEs:

$\begin{split}\frac{dx_{i}}{dt}=\bar{g}_{i}(\mathbf{x},\mathbf{y})&\triangleq\left(1-x_{i}-y_{i}\right)g_{i}(\mathbf{x}(t))-r_{i}(\mathbf{x})\\\
\frac{dy_{i}}{dt}=\bar{h}_{i}(\mathbf{x},\mathbf{y})&\triangleq\left(1-x_{i}-y_{i}\right)h_{i}(\mathbf{y}(t))-s_{i}(\mathbf{y})\end{split}$
(7)

for all $i\in\mathcal{N}$ and $t\geq 0$. In a matrix-vector form, (7) becomes:

$\begin{split}\frac{d\mathbf{x}}{dt}=\bar{G}(\mathbf{x},\mathbf{y})&\triangleq\text{diag}\left(\mathbf{1}-\mathbf{x}-\mathbf{y}\right)G(\mathbf{x})-R(\mathbf{x})\\\
\frac{d\mathbf{y}}{dt}=\bar{H}(\mathbf{x},\mathbf{y})&\triangleq\text{diag}\left(\mathbf{1}-\mathbf{x}-\mathbf{y}\right)H(\mathbf{y})-S(\mathbf{y}),\end{split}$
(8)

Where $G(\mathbf{x})=[g_{i}(\mathbf{x})]$,
$R(\mathbf{x})=[r_{i}(\mathbf{x})]$, and $H(\mathbf{y})=[h_{i}(\mathbf{y})]$,
$S(\mathbf{y})=[s_{i}(\mathbf{y})]$ are the non-linear infection and recovery
rate functions for viruses 1 and 2, respectively. The pairs $(G,R)$ and
$(H,S)$ each satisfy the assumptions (A1)–(A5); where $G$ and $H$ specifically
satisfy (A2) with respect to their corresponding graphs with adjacency
matrices $\mathbf{A}$ and $\mathbf{B}$, respectively. Figure 1 illustrates of
how these competing epidemics spread over the corresponding overlaid graphs.

Assumptions (A1)–(A5) are also more general (weaker) than those assumed in
[19, 20], where the recovery rates are restricted to being linear functions
and are thus a special case of our model. We emphasize that while the set off
assumptions for non-linear rates are mostly similar to (slightly more general
than) those in literature, the characterization of all convergence scenarios
for their respective bi-virus models is incomplete, as we shall discuss later
in Section VI.

##  IV Monotone Dynamical Systems and the Single Virus Epidemic

In this section, we provide a succinct introduction to monotone dynamical
systems (MDS) and some important definitions therein. We go on to show that
the $SIS$ model (6) is a monotone dynamical system (specifically a cooperative
system) and briefly apply these MDS techniques to epidemic models by deriving
the exact convergence result of the non-linear $SIS$ model. We also observe
that Theorem II.1 is a special case for when the infection and recovery rates
are linear.

###  IV-A Monotone Dynamical Systems - A Primer

A well known result from real analysis is that monotone sequences in compact
(closed and bounded) subsets of $\mathbb{R}^{n}$ converge in $\mathbb{R}^{n}$
[57]. This simple, yet powerful result has been fully integrated with the
theory of dynamical systems in a series of works [51, 58, 59, 60, 61, 62, 63,
64, 65, 66], which cumulatively form the theory of monotone dynamical systems
(MDS). The foundations of MDS were laid down in [51, 58, 59, 60, 61] which
study ordinary differential equations, specifically _cooperative_ ODE systems.
We here provide a brief, informal introduction to such ODE systems, with more
details in Appendix D.

A central tool in the theory of MDS is the notion of generalized cone-
orderings, which extends the concept of monotonicity in vector spaces.

###### Definition IV.1

Given a convex cone $K\subset X$ for any vector space $X$, the cone-ordering
$\leq_{K}$ ($<_{K}$, $\ll_{K}$) generated by $K$ is an order relation that
satisfies

  1. (i)

$~{}\mathbf{x}\\!\leq_{K}\\!\mathbf{y}\\!\iff\\!(\mathbf{y}\\!-\\!\mathbf{x})\in
K$;

  2. (ii)

$~{}\mathbf{x}\\!<_{K}\\!\mathbf{y}\\!\iff\\!\mathbf{x}\\!\leq_{K}\\!\mathbf{y}$
and $\mathbf{x}\\!\neq\\!\mathbf{y}$; and

  3. (iii)

$~{}\mathbf{x}\\!\ll_{K}\\!\mathbf{y}\\!\iff\\!(\mathbf{y}\\!-\\!\mathbf{x})\in\text{int}(K)$,
for any $\mathbf{x},\mathbf{y}\in X$.

Note that, ‘$\ll_{K}$’ implies ‘$<_{K}$’ and is a stronger relation. Cone-
orderings generated by the positive orthant $K\\!=\\!\mathbb{R}^{n}_{+}$ are
simply denoted by $\leq$ ($<,\ll$), that is, without the ‘$K$’ notation.

Let $\phi_{t}(\mathbf{x})$ denote the solution of a dynamical system at some
time $t\\!>\\!0$ starting from an initial point
$\phi_{0}(\mathbf{x})\\!=\\!\mathbf{x}\\!\in\\!\mathbb{R}^{n}$.

###### Definition IV.2

Given a cone-ordering $\leq_{K}$ ($<_{K}$, $\ll_{K}$), the dynamical system is
said to be _monotone_ if for every
$\mathbf{x},\mathbf{y}\\!\in\\!\mathbb{R}^{n}$ such that
$\mathbf{x}\\!\leq_{K}\\!\mathbf{y}$, we have
$\phi_{t}(\mathbf{x})\\!\leq_{K}\\!\phi_{t}(\mathbf{y})$ for all $t\\!>\\!0$.
The system is called _strongly monotone_ if for all
$\mathbf{x},\mathbf{y}\\!\in\\!\mathbb{R}^{n}$ such that
$\mathbf{x}\\!<_{K}\\!\mathbf{y}$, we have
$\phi_{t}(\mathbf{x})\\!\ll_{K}\\!\phi_{t}(\mathbf{y})$ for all $t\\!>\\!0$.

The main result from MDS theory says that (almost) every solution trajectory
of a strongly monotone system always converges to some equilibrium point of
the system [58, 64, 65, 44]. If the system has only one stable fixed point,
then this in itself is enough to prove global convergence. Monotonicity
properties of a dynamical system can therefore be leveraged as an alternative
to constructing Lyapunov functions, which is often intractable.

Consider the following autonomous ODE system

$\dot{\mathbf{x}}=\bar{F}(\mathbf{x}),$ (9)

where $\bar{F}(\mathbf{x})=[\bar{f}_{i}(\mathbf{x})]\in\mathbb{R}^{n}$ is the
vector field. If $\phi_{t}(\mathbf{x})$ is the solution of this ODE system, we
say the system is _co-operative_ if it is monotone. There are ways to find out
whether an ODE system is co-operative or not. In particular, one can answer
this by observing the Jacobian of the vector field [67]. The so-called Kamke
condition [66] says that (9) is co-operative with respect to the cone-ordering
generated by the positive orthant $K=\mathbb{R}^{n}_{+}$ if and only if

$\frac{\partial\bar{f}_{i}}{\partial x_{i}}\geq
0,~{}~{}~{}~{}~{}~{}~{}\text{for all }i\neq j.$ (10)

While it is not straightforward to obtain such a clean condition for any
general convex cone $K$, one can still deduce the co-operative property of the
ODE with respect to any one of the other orthants of $\mathbb{R}^{n}$ by
observing the signed entries of the Jacobian. We will show how this is done
for the bi-virus system (4) later in Section V-A.

If the Jacobian of an ODE system is an irreducible matrix in a subset $D$ of
the state space, we say that the ODE system is irreducible in $D$ (Definition
D.2 in Appendix D). If the ODE system is co-operative in $D$ as well as
irreducible in $D$, then it is strongly monotone in $D$ (Theorem D.4 in
Appendix D). To prove convergence properties, we should ideally be able to
show that our system is strongly monotone in the entirety of the state space
it is contained in, for which we can directly apply the main MDS convergence
result. However, this is often not the case, and one needs additional results
from MDS literature to prove convergence. These details are deferred to
Appendix D.

###  IV-B Monotonicity and convergence of SIS epidemic models

The following proposition establishes the monotonicity of the single-virus SIS
model with non-linear infection and recovery rates with respect to the regular
ordering relationship (cone-ordering generated by $R^{N}_{+}$).

###### Proposition IV.3

The ODE system (6) is cooperative in $[0,1]^{N}$ and irreducible in
$(0,1)^{N}$ with respect to the cone-ordering generated by the positive
orthant $\mathbb{R}^{N}_{+}$.$\hfill\square$

We now state the convergence criterion for the non-linear single-virus $SIS$
model.

###### Theorem IV.4

Let $\mathbf{J}_{F}(\mathbf{x})$ and $\mathbf{J}_{Q}(\mathbf{x})$ denote the
Jacobian matrices of the vector valued infection and recovery rate functions
$F(\mathbf{x})$ and $Q(\mathbf{x})$ from (6), respectively. Then,

  1. (i)

either $\lambda(\mathbf{J}_{F}(\mathbf{0})-\mathbf{J}_{Q}(\mathbf{0}))\leq 0$,
and $\mathbf{x}^{*}=0$ is the globally asymptotically stable fixed point of
(6);

  2. (ii)

or $\lambda(\mathbf{J}_{F}(\mathbf{0})-\mathbf{J}_{Q}(\mathbf{0}))>0$, and
there exists a unique, strictly positive fixed point $\mathbf{x}^{*}\gg 0$
such that $\mathbf{x}^{*}$ is globally asymptotically stable in
$[0,1]^{N}\setminus\\{\mathbf{0}\\}$.$\hfill\square$

The proof for Theorem IV.4 utilizes a result from the monotone dynamical
systems literature, provided as Theorem E.1 in Appendix E. It was originally
proved and applied to linear SIS epidemics in [68] as an alternate proof of
the convergence properties of the model for Gonorrhea spread in [2], which is
a special case of our non-linear model (6). We can also see this in the
following remark.

###### Remark IV.5

For the single-virus SIS model with linear infection and recovery rates (2),
the conditions derived in Theorem IV.4 reduce to those in Theorem II.1.

###### Proof:

By substituting $F(\mathbf{x})=\beta\mathbf{A}\mathbf{x}$ and
$Q(\mathbf{x})=\delta\mathbf{x}$ in (21) (Jacobian of the single-virus system
(6), mentioned in the proof of Theorem IV.4) and evaluating at
$\mathbf{x}=\mathbf{0}$, we get
$\mathbf{J}_{\bar{F}}(\mathbf{0})=\mathbf{J}_{F}(\mathbf{0})\\!-\\!\mathbf{J}_{Q}(\mathbf{0})=\beta\mathbf{A}\\!-\\!\delta\mathbf{I}$.
The condition
$\lambda(\mathbf{J}_{F}(\mathbf{0})\\!-\\!\mathbf{J}_{Q}(\mathbf{0}))=\lambda(\beta\mathbf{A}\\!-\\!\delta\mathbf{I})>0$
$(\leq 0)$ can be rewritten as $\tau>1/\lambda(\mathbf{A})$ $\left(\leq
1/\lambda(\mathbf{A})\right)$ where $\tau=\beta/\delta$, which as the same as
in Theorem II.1. ∎

While Theorem IV.4 could be proved using the steps in [2], which were
recreated again in [20], it requires first the application of two different
Lyapunov functions and also requires proving the uniqueness of the positive
fixed point. Alternatively, one could apply Theorem 1 in [69] to establish the
uniqueness of the positive fixed point by first showing that the Jacobian of
$\bar{F}(\mathbf{x})$ evaluated at any point $\mathbf{x}\gg\mathbf{0}$
satisfying $\bar{F}(\mathbf{x})=\mathbf{0}$, is Hurwitz. This, combined with
Proposition IV.3, could then provide the necessary convergence criterion.
However, we maintain that using Theorem E.1 would be a simpler way to derive
the same results, whose proof is deferred to Appendix E.

##  V Main results for the non-linear Bi-Virus model

We provide the necessary and sufficient results on the non-linear infection
and recovery rates of the bi-virus system (8) for convergence to each of the
three different kinds of equilibria: the virus-free, the single-virus
equilibrium, and the co-existence equilibrium. However, before stating the
main convergence results (proofs deferred to Appendix F in [70]), we establish
the monotonicity of the non-linear bi-virus model.

###  V-A Monotonicity of the Bi-Virus epidemic models

We first revisit the Kamke condition from Section IV-A, in this instance given
for a the _southeast cone-ordering_ as stated below.

#### Southeast cone-ordering and the Kamke condition

Consider the cone-ordering generated by the convex cone
$K=\\{\mathbb{R}^{N}_{+}\times\mathbb{R}^{N}_{-}\\}\subset\mathbb{R}^{2N}$.
This cone is one of the orthants of $\mathbb{R}^{2N}$, and for $N=1$, it would
correspond to the southeast orthant of $\mathbb{R}^{2}$
$\left(K=\\{\mathbb{R}_{+}\times\mathbb{R}_{-}\\}\subset\mathbb{R}^{2}\right)$.
For any two points $(\mathbf{x},\mathbf{y})$,
$(\bar{\mathbf{x}},\bar{\mathbf{y}})\in\mathbb{R}^{2N}$, it satisfies the
following:

  1. (i)

$(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\iff\\!x_{i}\\!\leq\\!\bar{x}_{i}$
and $y_{i}\\!\geq\\!\bar{y}_{i}$ for all $i\\!\in\\!\mathcal{N}$;

  2. (ii)

$(\mathbf{x},\mathbf{y})\\!\\!<_{K}\\!\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\\!\iff\\!\\!(\mathbf{x},\mathbf{y})\\!\\!\leq_{K}\\!\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})$
and $(\mathbf{x},\mathbf{y})\\!\neq\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})$;

  3. (iii)

$(\mathbf{x},\mathbf{y})\\!\ll_{K}\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\iff\\!x_{i}\\!<\\!\bar{x}_{i}$
and $y_{i}\\!>\\!\bar{y}_{i}$ for all $i\\!\in\\!\mathcal{N}$.

This type of cone-ordering is often referred to as the southeast cone-
ordering, and the corresponding cone $K$ is the southeast orthant of
$\mathbb{R}^{2N}$. As shown in [67], the Kamke condition for determining
whether an ODE system is cooperative or not with respect to the positive
orthant $\mathbb{R}^{2N}_{+}$ can be generalised for cone-orderings generated
by any orthant of $\mathbb{R}^{2N}$, including the southeast orthant. Once
again, this is done by observing the Jacobian of the respective ODE system.
Consider the $2N$ dimensional system given by

$\dot{\mathbf{x}}=\bar{G}(\mathbf{x},\mathbf{y})~{}~{}\text{and}~{}~{}\dot{\mathbf{y}}=\bar{H}(\mathbf{x},\mathbf{y}),$

where $\bar{G}(\mathbf{x},\mathbf{y})=[\bar{g}_{i}(\mathbf{x},\mathbf{y})]$
and $\bar{H}(\mathbf{x},\mathbf{y})=[\bar{h}_{i}(\mathbf{x},\mathbf{y})]$ are
vector-valued functions in $\mathbb{R}^{N}$. The Kamke condition for this
system with respect to the southeast cone-ordering [67] is

$\frac{\partial\bar{g}_{i}}{\partial x_{j}}\geq
0,~{}\frac{\partial\bar{h}_{i}}{\partial y_{j}}\geq 0,~{}\forall i\neq
j,~{}~{}~{}\text{and}~{}~{}~{}\frac{\partial\bar{g}_{i}}{\partial y_{j}}\leq
0,~{}\frac{\partial\bar{h}_{i}}{\partial x_{j}}\leq 0,~{}\forall i,j.$

Roughly speaking, the Jacobian $\mathbf{J}_{GH}(\mathbf{x},\mathbf{y})$ of the
system, evaluated at all points in the state space, should be in the following
block matrix form (where the signs are not strict):

$\mathbf{J}_{\bar{G}\bar{H}}=\begin{bmatrix}*&+&+&-&-&-\\\ +&*&+&-&-&-\\\
+&+&*&-&-&-\\\ -&-&-&*&+&+\\\ -&-&-&+&*&+\\\ -&-&-&+&+&*\end{bmatrix}$ (11)

Note that the state space of the ODE system (4) is given by
$D\triangleq\left\\{(\mathbf{x},\mathbf{y})\in[0,1]^{2N}~{}|~{}\mathbf{x}+\mathbf{y}\leq\mathbf{1}\right\\}$.

###### Proposition V.1

The ODE system (8) (the non-linear bi-virus model) is cooperative in $D$ with
respect to the southeast cone-ordering. It is also irreducible in
$\text{Int}(D)$.

###### Proof:

For all $(\mathbf{x},\mathbf{y})\in D$ and $i\neq j\in\mathcal{N}$, we have

$\frac{\partial\bar{g}_{i}(\mathbf{x},\mathbf{y})}{\partial
x_{j}}=(1-x_{i}-y_{i})\frac{\partial g_{i}(\mathbf{x})}{\partial
x_{j}}-\frac{\partial r_{i}(\mathbf{x})}{\partial x_{j}}\geq 0,$
$\frac{\partial\bar{h}_{i}(\mathbf{x},\mathbf{y})}{\partial
y_{j}}=(1-x_{i}-y_{i})\frac{\partial h_{i}(\mathbf{y})}{\partial
y_{j}}-\frac{\partial s_{i}(\mathbf{x})}{\partial y_{j}}\geq 0$

since $\frac{\partial g_{i}(\mathbf{x})}{\partial x_{j}}\geq 0$,
$\frac{\partial r_{i}(\mathbf{x})}{\partial x_{j}}\leq 0$ and $\frac{\partial
h_{i}(\mathbf{y})}{\partial y_{j}}\geq 0$, $\frac{\partial
s_{i}(\mathbf{y})}{\partial y_{j}}\leq 0$ from assumptions (A2) and (A3), and
$(1-x_{i}-y_{i})\geq 0$. Moreover for all $i\in\mathcal{N}$,

$\frac{\partial\bar{g}_{i}}{\partial y_{i}}=-g_{i}(\mathbf{x})\leq
0~{}~{}\text{and}~{}~{}\frac{\partial\bar{h}_{i}}{\partial
x_{i}}=-h_{i}(\mathbf{y})\leq 0,$

with ${\partial\bar{g}_{i}}/{\partial y_{j}}={\partial\bar{h}_{i}}/{\partial
x_{j}}=0$. Thus, the Kamke conditions are satisfied and the system is
cooperative in $D$.

The Jacobian $\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{x},\mathbf{y})$ of system
(4) is written as

$\begin{split}&\\!\\!\\!\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{x},\mathbf{y})=\\\
&\\!\\!\\!\\!\\!\begin{bmatrix}\mathbf{S}_{\mathbf{x}\mathbf{y}}\mathbf{J}_{G}(\mathbf{x})\\!\\!-\\!\\!\mathbf{D}_{G(\mathbf{x})}\\!\\!-\\!\\!\mathbf{J}_{R}(\mathbf{x})&\\!\\!\\!\\!-\\!\mathbf{D}_{G(\mathbf{x})}\\\
-\\!\mathbf{D}_{H(\mathbf{y})}&\\!\\!\\!\\!\mathbf{S}_{\mathbf{x}\mathbf{y}}\mathbf{J}_{H}(\mathbf{y})\\!\\!-\\!\\!\mathbf{D}_{H(\mathbf{y})}\\!\\!-\\!\\!\mathbf{J}_{S}(\mathbf{y})\end{bmatrix}\\!\\!,\end{split}$
(12)

where
$\mathbf{S}_{\mathbf{x},\mathbf{y}}\triangleq\text{diag}(\mathbf{1}-\mathbf{x}-\mathbf{y})$,
$\mathbf{D}_{G(\mathbf{x})}\triangleq\text{diag}(G(\mathbf{x}))$ and
$\mathbf{D}_{H(\mathbf{y})}\triangleq\text{diag}(H(\mathbf{y}))$. Since the
infection rate functions satisfy assumption (A2) for their corresponding
underlying graphs, $\mathbf{J}_{G}(\mathbf{x})$ and
$\mathbf{J}_{H}(\mathbf{y})$ follow the sign structure of $\mathbf{A}$ and
$\mathbf{B}$ respectively and are irreducible. The off-diagonal blocks of
$\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{x},\mathbf{y})$ are diagonal matrices
with non-zero diagonal entries for $(\mathbf{x},\mathbf{y})\in\text{Int}(D)$,
and there does not exist a permutation matrix that would transform this into a
block upper triangular matrix. Hence, by Definition D.2, the system is
irreducible in $\text{Int}(D)$, and this completes the proof. ∎

From Proposition V.1, we deduce that the non-linear bi-virus system of ODEs
(8) is co-operative in $D$, and thus strongly monotone in $\text{Int}(D)$ in
view of Theorem D.4 in Appendix D. This property also extends to the linear
bi-virus system (4) which is a special case of (8).

###  V-B Convergence and Coexistence properties of the Bi-Virus model

We are now ready to establish results on convergence properties of the bi-
virus model and provide conditions for coexistence of two viruses in the non-
linear bi-virus model as in (8).

Let $\mathbf{x}^{*}$ and $\mathbf{y}^{*}$ be the globally attractive fixed
points of the single-virus SIS models that system (8) would reduce to when
Virus 2 and 1, respectively, are not present over the network. These systems
are given by

$\dot{\mathbf{x}}=F^{x}(\mathbf{x})\triangleq\bar{G}(\mathbf{x},\mathbf{0})=\text{diag}(\mathbf{1}-\mathbf{x})G(\mathbf{x})-R(\mathbf{x}),$
(13)
$\dot{\mathbf{y}}=F^{y}(\mathbf{y})\triangleq\bar{H}(\mathbf{0},\mathbf{y})=\text{diag}(\mathbf{1}-\mathbf{y})H(\mathbf{y})-S(\mathbf{y});$
(14)

and by Theorem IV.4, $\mathbf{x}^{*}\\!=\\!\mathbf{0}$
($\mathbf{y}^{*}\\!=\\!\mathbf{0}$) if
$\lambda\left(\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!\leq\\!0$
(if
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!\leq\\!0$),
and $\mathbf{x}^{*}\\!\gg\\!\mathbf{0}$ ($\mathbf{y}^{*}\\!\gg\\!\mathbf{0}$)
otherwise.

We first state the result when the virus-free equilibrium is globally
attractive. We prove this by presenting simple arguments which require only
Theorem IV.4 for SIS model along with the monotonicity properties derived in
the previous section, eliminating the need of a Lyapunov based approach.

######  Theorem V.2 (Convergence to virus-free equilibria)

If
$\lambda\left(\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!\leq\\!0$
and
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!\leq\\!0$,
trajectories of (8) starting from any point in $D$ converge to
$(\mathbf{0},\mathbf{0})$.$\hfill\square$

We next characterize the conditions when the system globally converges to
equilibria when only one of the viruses survives over the network. Let
$\mathbf{S}_{\mathbf{x}}\\!\triangleq\\!\text{diag}(\mathbf{1}\\!-\\!\mathbf{x})$
and
$\mathbf{S}_{\mathbf{y}}\\!\triangleq\\!\text{diag}(\mathbf{1}\\!-\\!\mathbf{y})$
for any $\mathbf{x},\mathbf{y}\in\mathbb{R}^{N}$. Also denote by
$B_{x}\\!\triangleq\\!\left\\{(\mathbf{x},\mathbf{y})\in
D~{}|~{}\mathbf{x}\\!>\\!\mathbf{0}\right\\}$ the set of all points
$(\mathbf{x},\mathbf{y})\\!\in\\!D$ for which $x_{i}\\!>\\!0$ for some
$i\in\mathbb{N}$, and let
$B_{y}\\!\triangleq\\!\left\\{(\mathbf{x},\mathbf{y})\in
D~{}|~{}\mathbf{y}\\!>\\!\mathbf{0}\right\\}$ be a similar set for the $y_{i}$
entries.

######  Theorem V.3 (Convergence to single-virus equilibria)

When
$\lambda\\!\left(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!>\\!0$
and
$\lambda\\!\left(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!\leq\\!0$,
$(\mathbf{x}^{*},\mathbf{0})$ is globally attractive in $B_{x}$;777We consider
$B_{x}$ as the global domain of attraction instead of $D$ because
$\mathbf{x}=0$ for all points in the set $D\setminus B_{x}$. Starting from
such points the system is no longer a bi-virus epidemic, but a single-virus
SIS system for Virus 2. that is, every trajectory of system (8) starting from
points in $B_{x}$ converges to $(\mathbf{x}^{*},\mathbf{0})$.

Similarly, when
$\lambda\left(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(\mathbf{0})-\mathbf{J}_{R}(\mathbf{0})\right)\leq
0$ and
$\lambda\left(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{J}_{H}(\mathbf{0})-\mathbf{J}_{S}(\mathbf{0})\right)>0$
is globally attractive in $B_{y}$. $\hfill\square$

###### Proof:

The idea behind the proof is illustrated in Figure 2. For every
$(\mathbf{x},\mathbf{y})\\!\in\\!B_{x}$ (for example $p_{1}$ and $p_{2}$ in
Figure 2), we construct a point $(\mathbf{x}_{r},\mathbf{y}_{s})$ which
eventually bounds the trajectory starting from $(\mathbf{x},\mathbf{y})$; that
is, we have
$(\mathbf{x}_{r},\mathbf{y}_{s})\\!\ll_{K}\\!\phi_{t_{1}}(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\mathbf{x}^{*},\mathbf{0})$888$\phi_{t}(\mathbf{x},\mathbf{y})$
denotes the solution of (4) at $t\\!\geq\\!0$, with initial point
$(\mathbf{x},\mathbf{y})$. for some $t_{1}\\!\geq\\!0$. From the monotonicity
shown in Proposition V.1, we have
$\phi_{t}(\mathbf{x}_{r},\mathbf{y}_{s})\\!\ll_{K}\\!\phi_{t+t_{1}}(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\mathbf{x}^{*},\mathbf{0})$
for all time $t\\!\geq\\!0$. We prove that the trajectory starting from
$(\mathbf{x}_{r},\mathbf{y}_{s})$ converges to $(\mathbf{x}^{*},0)$
monotonically, with respect to the southeast cone-ordering (Figure 2(a)).
Using this, we show the convergence of trajectories starting from
$(\mathbf{x},\mathbf{y})$ via a sandwich argument (Figure 2(b)). See Appendix
F in [70] for detailed proof. ∎

(a) For every point $p_{k}$, there is a point
$(\mathbf{x}_{rk},\mathbf{y}_{sk})$ starting from which, trajectories converge
monotonically $(\leq_{K})$ to $(\mathbf{x}^{*},0)$. (b) Trajectories starting
from $p_{k}$ eventually bounded by $(\mathbf{x}_{rk},\mathbf{y}_{sk})$;
monotonicity of the system gives convergence to $(\mathbf{x}^{*},0)$.

Figure 2: Illustration of the convergence to $(\mathbf{x}^{*},0)$

(a) Limitations of the literature. (b) Complete characterization of the
convergence trichotomy.

Figure 3: Characterization of the parameter space

Finally, we give the necessary and sufficient conditions that guarantee the
co-existence of the two viruses in the long run. Let $E$ denote the set of all
fixed points of the system in (8).

######  Theorem V.4 (Convergence to coexistence equilibria)

If
$\lambda\left(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!>\\!0$
and
$\lambda\left(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!>\\!0$,
there exist fixed points of system (8)
$(\hat{\mathbf{x}},\hat{\mathbf{y}})\\!\gg\\!(\mathbf{0},\mathbf{0})$ and
$(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\gg\\!(\mathbf{0},\mathbf{0})$ such
that

$(\mathbf{0},\mathbf{y}^{*})\ll_{K}(\hat{\mathbf{x}},\hat{\mathbf{y}})\leq_{K}(\bar{\mathbf{x}},\bar{\mathbf{y}})\ll_{K}(\mathbf{x}^{*},\mathbf{0}),$

with the possibility that
$(\hat{\mathbf{x}},\hat{\mathbf{y}})=(\bar{\mathbf{x}},\bar{\mathbf{y}})$. All
trajectories of system (8) starting from $B_{x}\cap B_{y}$ converge to the set
of coexistence fixed points
$S\triangleq\left\\{(\mathbf{x}_{e},\mathbf{y}_{e})\\!\in\\!E~{}|~{}(\hat{\mathbf{x}},\hat{\mathbf{y}})\\!\leq_{K}\\!(\mathbf{x}_{e},\mathbf{y}_{e})\\!\leq_{K}\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})\right\\}$.$\hfill\square$

The proof of Theorem V.4 follows similar arguments to that of the previous
theorem, and is the first convergence result for coexistence fixed points in
the competing SIS literature. Note that while we have convergence to ‘a’
coexistence equilibrium, it may or may not be unique in the state space. The
global convergence is therefore to the set of possible coexistence equilibria,
and not necessarily a singular point. Thus, via Theorems V.2, V.3 and V.4 we
cover all possible convergence scenarios of the bi-virus SIS system (8), and
successfully establish the complete theoretical characterization for the
trichotomy of possible outcomes.

##  VI Linear Infection and Recovery rates - Discussion and Comparison to
Literature

We now take a look at the special case of the bi-virus epidemic model where
infection and recovery rates scale linearly with the local infection
probability. This is the most commonly analysed setting in literature [21, 54,
31, 32, 33, 34], and allows us to provide a comprehensive discussion on the
related works. With the exception of [54], a line of work seemingly developed
concurrently to ours, we observe that most existing works only provide limited
results regarding convergence to coexistence equilibria. In what follows, we
provide corollaries of Theorems V.2, V.3 and V.4 which characterize
convergence to the trichotomy of possible outcomes for the special case of
linear infection and recovery rates. These results, along with Figure 3, are
reproduced here as they originally were in our previous work [1] which focused
only on characterizing the convergence properties in the case of linear
infection and recovery rates.

The model considered in this section is the bi-virus system (4) with
homogeneous infection and recovery rates999every infected node
$i\in\mathcal{N}$ infects its susceptible neighbor with the same rate
$\beta_{1}>0$ or $\beta_{2}>0$, and in turn recovers with the same rate
$\delta_{1}>0$ or $\delta_{2}>0$, depending on whether it is infected by Virus
1 or 2 respectively.. While at first this may seem too simplistic compared to
the case of linear, heterogeneous rates101010The adjacency matrices
$\mathbf{A}$ and $\mathbf{B}$ in (4) can be symmetric, irreducible, weighted;
with $a_{ij},b_{ij}\geq 0$ (not necessarily $0/1$ valued) multiplied by
$\beta_{1}$ and $\beta_{2}$ respectively, being the infection rates from node
$j\to i$ for Viruses 1 and 2. Recovery rates can similarly be heterogenized as
$\boldsymbol{\delta}_{1}=[\delta_{1}^{i}]$ and
$\boldsymbol{\delta}_{2}=[\delta_{2}^{i}]$ for Viruses 1 and 2; written as
recovery rate matrices $\text{diag}(\boldsymbol{\delta}_{1})$ and
$\text{diag}(\boldsymbol{\delta}_{1})$, respectively., and even generic, non-
linear rates analyzed in literature [19, 20, 21, 54, 31, 32, 33, 34], the
discussions in the ‘Comparison to existing ilterature’ subsection will still
hold for these more general cases. We only stick to the bi-virus system with
homogeneous rates as in (4) to be able to illustrate our results in the form
of Figure 3; the axes capturing the parameters of the system. This enables us
to better explain our contribution, using visual aids in the form of Figure 3,
helping us compare our work with some of the existing literature more
effectively, as opposed to presenting any other special case of the bi-virus
model.

Consider the linear bi-virus system (4). By setting
$G(\mathbf{x})=\beta_{1}\mathbf{A}\mathbf{x}$,
$R(\mathbf{x})=\delta_{1}\mathbf{x}$ and
$H(\mathbf{y})=\beta_{2}\mathbf{B}\mathbf{y}$,
$S(\mathbf{y})=\delta_{2}\mathbf{y}$, we get

$\mathbf{J}_{G}(\mathbf{0})\\!=\\!\beta_{1}\mathbf{A},~{}~{}\mathbf{J}_{R}(\mathbf{0})\\!=\\!\delta_{1}\mathbf{I},$

and

$\mathbf{J}_{H}(\mathbf{0})\\!=\\!\beta_{2}\mathbf{B},~{}~{}\mathbf{J}_{S}(\mathbf{0})\\!=\\!\delta_{2}\mathbf{I}.$

Defining $\tau_{1}\triangleq\beta_{1}/\delta_{1}$,
$\tau_{2}=\triangle\beta_{2}/\delta_{2}$, and plugging in the above
expressions for the Jacobians in Theorems V.2 and V.3, we have the following
Corollaries.

###### Corollary VI.1

If $\tau_{1}\lambda(\mathbf{A})\\!\leq\\!1$ and
$\tau_{2}\lambda(\mathbf{B})\\!\leq\\!1$, trajectories of (4) starting from
any point in $D$ converge to $(\mathbf{0},\mathbf{0})$.$\hfill\square$

###### Corollary VI.2

When $\tau_{1}\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})\\!>\\!1$ and
$\tau_{2}\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})\\!\leq\\!1$,
$(\mathbf{x}^{*},\mathbf{0})$ is globally attractive in $B_{x}$;111111We
consider $B_{x}$ as the global domain of attraction instead of $D$ because
$\mathbf{x}=0$ for all points in the set $D\setminus B_{x}$. Starting from
such points the system is no longer a bi-virus epidemic, but a single-virus
SIS system for Virus 2. that is, every trajectory of system (4) starting from
points in $B_{x}$ converges to $(\mathbf{x}^{*},\mathbf{0})$.

Similarly, when
$\tau_{1}\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})\\!\leq\\!1$ and
$\tau_{2}\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})\\!>\\!1$,
$(\mathbf{0},\mathbf{y}^{*})$ is globally attractive in $B_{y}$.
$\hfill\square$

From Corollary VI.2, we can deduce that the threshold values for $\tau_{1}$
and $\tau_{2}$ below which each of the viruses will die out are given by the
equations $\tau_{1}\\!=\\!1/\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})$
and $\tau_{2}\\!=\\!1/\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})$,
respectively. Figure 3(b) plots these threshold values for Virus 1 (in blue)
and Virus 2 (in red) for varying values of $\tau_{1}$ and $\tau_{2}$, and
partitions the entire parameter space into regions R1 – R6 as shown. When
$\tau_{1}\\!>\\!1/\lambda(\mathbf{A})$ and
$\tau_{2}\\!>\\!1/\lambda(\mathbf{B})$, for which values of
$\tau_{1},\tau_{2}$ do not lie in regions R1, R2 or R3, the blue curve lies
above the red curve as in Figure 3(b). This was originally shown in [18] by
deducing that the ratio of slopes of the red and blue curves at point
$(\tau_{1},\tau_{2})=\left(1/\lambda(\mathbf{A}),1/\lambda(\mathbf{B})\right)$
is less than one. This means there exist combinations of $\tau_{1},\tau_{2}$
for which $\tau_{1}$ lies to the right of the blue curve
($\tau_{1}\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})\\!>\\!1$), and
$\tau_{2}$ lies above the red curve
($\tau_{2}\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})\\!>\\!1$).121212Note
that $\tau_{1}\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})\\!\leq\\!1$ and
$\tau_{2}\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})\\!\leq\\!1$ is only
possible in region R1, since it is the only region where $\tau_{1}$ can lie to
the left of the blue curve, and $\tau_{2}$ can lie below the red curve. This
effectively reduces the expressions to
$\tau_{1}\lambda(\mathbf{A})\\!\leq\\!1$ and
$\tau_{2}\lambda(\mathbf{B})\\!\leq\\!1$, the conditions for convergence to
the virus-free equilibrium as in Corollary VI.1. This corresponds to region R6
in Figure 3(b), and our final corollary (derived from Theorem V.4) shows that
for values of $\tau_{1},\tau_{2}$ which lie in R6, we observe convergence to
coexistence equilibria.

######  Corollary VI.3 (Convergence to coexistence equilibria)

If $\tau_{1}\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})\\!>\\!1$ and
$\tau_{2}\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})\\!>\\!1$, there exist
fixed points of system (4)
$(\hat{\mathbf{x}},\hat{\mathbf{y}})\\!\gg\\!(\mathbf{0},\mathbf{0})$ and
$(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\gg\\!(\mathbf{0},\mathbf{0})$ such
that

$(\mathbf{0},\mathbf{y}^{*})\ll_{K}(\hat{\mathbf{x}},\hat{\mathbf{y}})\leq_{K}(\bar{\mathbf{x}},\bar{\mathbf{y}})\ll_{K}(\mathbf{x}^{*},\mathbf{0}),$

with the possibility that
$(\hat{\mathbf{x}},\hat{\mathbf{y}})=(\bar{\mathbf{x}},\bar{\mathbf{y}})$. All
trajectories of system (4) starting from $B_{x}\cap B_{y}$ converge to the set
of coexistence fixed points
$S\triangleq\left\\{(\mathbf{x}_{e},\mathbf{y}_{e})\\!\in\\!E~{}|~{}(\hat{\mathbf{x}},\hat{\mathbf{y}})\\!\leq_{K}\\!(\mathbf{x}_{e},\mathbf{y}_{e})\\!\leq_{K}\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})\right\\}$.$\hfill\square$

| $g_{i}(\mathbf{x})$ | $h_{i}(\mathbf{y})$ | $r_{i}(\mathbf{x})$ | $s_{i}(\mathbf{y})$
---|---|---|---|---
CASE 1 | $\sum_{j}a_{ij}x_{j}$ | $\sum_{j}b_{ij}y_{j}$ | $\delta_{1}x_{i}$ | $\delta_{2}y_{i}$
CASE 2 | $\sum_{j}a_{ij}\ln(1+\alpha_{1}x_{j})$ | $\sum_{j}b_{ij}\ln(1+\alpha_{2}y_{j})$ | $\delta_{1}x_{i}$ | $\delta_{2}y_{i}$
CASE 3 | $\sum_{j}a_{ij}\ln(1+\alpha_{1}x_{j})$ | $\sum_{j}b_{ij}\ln(1+\alpha_{2}y_{j})$ | $(1+x_{i})^{2}-1$ | $(1+y_{i})^{2}-1$

Table I: Summary of infection and recovery rate functions chosen.

#### Comparison to existing literature

Now that we have established all our results, we briefly compare our work with
results from [20, 19], which also talk about global convergence to single-
virus equilibria. To this end, we first illustrate the limitations of the
existing conditions for global convergence in [20, 19] in Figure 3(a); and use
Figure 3(b), where we provide complete characterization of the parameter
space, to draw comparisons with our results. We then discuss the works [31,
34, 32, 33] which consider more general models where there can be more than
two viruses, but present sharper results in the bi-virus setting. Finally, we
will briefly comment on the finiteness of the coexistence equilibria, citing
results from [54].

When translated to the setting of linear infection and recovery rates as in 4,
the result from [19] says that when
$\tau_{1}d_{min}(\mathbf{A})\\!>\\!\tau_{2}d_{max}(\mathbf{B})$, the Virus 2
is sure to die out (Virus 1 could persist or die out), and similarly when
$\tau_{1}d_{max}(\mathbf{A})\\!<\\!\tau_{2}d_{min}(\mathbf{B})$, the Virus 1
is sure to die out. We illustrate these conditions in Figure 3(a), where Virus
1 (Virus 2) is sure to die out if parameters ($\tau_{1},\tau_{2}$) lie above
(below) the blue (red) line. Therefore, the entire yellow-shaded region in
Figure 3(a), between the blue and red lines, is left uncharacterized in [19].
When $\mathbf{A}$ and $\mathbf{B}$ are regular graphs with the same degree
($d_{min}\\!=\\!d_{max}\\!=\\!d$), the blue and red lines coincide, making
coexistence infeasible. This is also mentioned in [18] where they show that
for regular graphs with same degree, the system behaves as if the two graphs
were the same - rendering coexistence impossible (which is also in line with
results in [8]). In contrast, the maximum degree of graphs can also be much
larger than the minimum degree (e.g., power law graphs), causing the yellow-
shaded space to become very large, possibly spanning almost the entire
parameter space.

The main result in [20], when similarly translated to our setting as above,
says that when $\tau_{1}\lambda(\mathbf{A})\\!>\\!1$ and
$\tau_{2}\lambda(\mathbf{B})\\!\leq\\!1$, Virus 1 survives and Virus 2 dies
out. Similarly, when $\tau_{2}\lambda(\mathbf{B})\\!>\\!1$ and
$\tau_{1}\lambda(\mathbf{A})\\!\leq\\!1$, Virus 2 survives and Virus 1 dies
out. These correspond to regions R2 and R3 in Figure 3(b). However, their
results do not cover the convergence properties for $\tau_{1},\tau_{2}$ which
lie in regions R4 – R6. Our Theorems V.3 and V.4, through their corresponding
corollaries, do account for these values of $\tau_{1},\tau_{2}$, and show
convergence to $(\mathbf{0},\mathbf{y}^{*})$, $(\mathbf{x}^{*},\mathbf{0})$ or
to a coexistence fixed point whenever they lie in regions R4, R5, or R6,
respectively.

The works [32, 33] consider the bi-virus epidemic model with heterogeneous
linear infection and recovery rates as a special case of their respective
multi-virus models. Corollary 2 in [33], a more general version of Theorem 5
in [32] which considers the case where $N=2$, establishes existence conditions
for the coexistence equilibria. These conditions are identical to the ones
emerging out of Theorem V.4 when applied to the bi-virus model considered
therein (also identical to the conditions in Corollary VI.3 for the special
case of homogeneous, linear infection and recovery rates), and our result can
therefore be considered as an extension of those in [32, 33]; providing
_convergence_ results in addition to their existence results. Theorem 6 in
[34] (Theorem 8 in [31]) is another interesting result concerning coexistence
equilibria, where they show for the special case of viruses spreading over the
same (possibly weighted) graph that the survival probability vectors of both
the viruses are the same up to a constant multiple; that is, they are
parallel.

The finiteness of the number of single-virus equilibria is evident from
Theorem IV.4, which proves its uniqueness. However, Theorem V.4 and Corollary
VI.3 do not explicitly show that coexistence equilibria are finitely many, let
alone uniqueness131313In Section VII, we show with the aid of simulation
results that the coexistence equilibria are indeed not unique in general.. For
linear, heterogeneous infection and recovery rates, Theorem 3.6 in [54] uses
novel techniques from algebraic geometry to prove that the coexistence
equilibria are finitely many for all possible values of infection and recovery
rates that do not lie in an algebraic set of measure zero. However, this
remains an open problem for general, non-linear infection and recovery rate
functions satisfying (A1)–(A5).

In summary, without our Theorems V.3 and V.4, convergence results from
literature fail to characterize a sizeable portion of the parameter space as
shown in Figure 3(a) by the ‘?’ region (part of the shaded region surrounded
by the arrows). The parameters leading to coexistence are entirely contained
in this region as well - explaining the dearth of convergence results for such
equilibria in the existing literature.

##  VII Numerical Results

In this section, we present simulation results to support our theoretical
findings for the bi-virus SIS model for combinations of non-linear as well as
linear infection and recovery rates. To this end, we consider an undirected,
connected graph (103 nodes, 239 edges), called Autonomous System (AS-733),
from the SNAP repository [71]. For both the linear and non-linear bi-virus
model, we generate an additional graph, overlaid on the same set of nodes, by
modifying the original graph (AS-733-A with
$\lambda(\mathbf{A})\\!=\\!12.16$), removing and adding edges while ensuring
connectivity between the nodes. The new additional graph, AS-733-B, has 741
edges with $\lambda(\mathbf{B})\\!=\\!15.53$. Note that since our theoretical
results hold for any general graphs, we only use this set as example graphs to
numerically demonstrate the convergence properties. Similar numerical results
can indeed be obtained for any other networks (such as social networks).

We test the convergence dynamics of the bi-virus model over a range of
combinations of linear and non-linear infection and recovery rates. To this
end, we consider three different bi-virus models, and Table I summarizes the
three cases with the corresponding infection and recovery rate functions as
shown. Note that for non-linear infection and recovery rates, we consider the
logarithmic and polynomial functions briefly mentioned in Section III, to
ensure that our three cases satisfy assumptions (A1)–(A5).

For each of the three cases, we construct combinations of parameters
($\tau_{1}$ or $\tau_{2}$ for linear rates, and $\alpha_{1}$ or $\alpha_{2}$
for non-linear rates), to develop three convergence scenarios, that satisfy
the assumptions of Theorems V.3 and V.4. These three scenarios correspond to
global convergence of the bi-virus system to fixed points where (a) Virus 1 is
the surviving epidemic (which spreads on graph AS-733-A), (b) Virus 2 is the
surviving epidemic (which spreads on graph AS-733-B), (c) both viruses
coexist, (where Virus 1 spreads on graph AS-733-A and Virus 2 on AS-733-B).
Parameters corresponding to these three scenarios are provided in the table
inset in Figures 4–6(a)–(c) corresponding to the three cases.

To visualize our system in two dimensions, we use
$avgX\\!\triangleq\\!(1/N)\sum_{i\in\mathcal{N}}x_{i}$ on the x-axis, and
$avgY\\!\triangleq\\!(1/N)\sum_{i\in\mathcal{N}}y_{i}$ on the y-axis. We plot
trajectories of the bi-virus system starting from different initial points in
the state space $D$ to observe their convergence, with red arrows representing
the trajectories’ direction of movement at various time intervals. Here, the
state space $D$ is the region that lies below the dotted-line (for example, in
Figure 4), ensuring $x_{i}+y_{i}\\!\leq\\!1$ for all $i\in\mathcal{N}$, for
every initial point. To ensure that the convergences observed in our phase
plots match the conditions laid out in Theorems V.3 and V.4, we track the
eigenvalues
$\lambda(\mathbf{U})\triangleq\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(0)-\mathbf{J}_{R}(0))$
and
$\lambda(\mathbf{V})\triangleq\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{J}_{H}(0)-\mathbf{J}_{S}(0))$.
$\lambda(\mathbf{U})$ ($\lambda(\mathbf{V})$) being positive or negative
corresponds to Virus 1 (Virus 2) surviving or dying out, respectively.

(a) $\lambda(U)>0$, $\lambda(V)<0$; Virus 1 survives (b) $\lambda(U)<0$,
$\lambda(V)>0$; Virus 2 survives

(c) $\lambda(U)>0$, $\lambda(V)>0$; Both coexist

Figure 4: Phase plots for a system with linear infection and recovery rates
(CASE 1) on the AS-733 graph.

(a) $\lambda(U)>0$, $\lambda(V)<0$; Virus 1 survives (b) $\lambda(U)<0$,
$\lambda(V)>0$; Virus 2 survives

(c) $\lambda(U)>0$, $\lambda(V)>0$; Both coexist

Figure 5: Phase plots for a system with non-linear infection and linear
recovery rates (CASE 2) on the AS-733 graph.

(a) $\lambda(U)>0$, $\lambda(V)<0$; Virus 1 survives (b) $\lambda(U)<0$,
$\lambda(V)>0$; Virus 2 survives

(c) $\lambda(U)>0$, $\lambda(V)>0$; Both coexist

Figure 6: Phase plots for a system with non-linear infection and recovery
rates (CASE 3) on the AS-733 graph. Figure 7: Coexistence condition with
Multiple equilibrium points

In Figures 4–6(a)–(c), we show numerical results for the three cases,
respectively. Figures 4–6(a) and 4–6(b) show convergence to the two different
single-virus equilibria, where the parameters therein satisfy the two set of
conditions as in Theorem V.3. Figures 4–6(c) show convergence to the
coexistence equilibria, which also satisfies the coexistence conditions as
outlined in Theorem V.4. We observe a unique coexistence equilibrium when the
viruses are competing over graphs AS-733-A and AS-733-B, for which the
eigenvalues $\lambda(\mathbf{A})$ and $\lambda(\mathbf{B})$ are significantly
different. Interestingly, we also observe multiple coexistence equilibria as
shown in Figure 7. We obtain this result by creating another additional graph
by modifying the original graph AS-733-A such that the eigenvalue of this new
graph is as close to the original one where this new graph AS-733-C has 259
edges with $\lambda(\mathbf{C})\\!\\!=\\!\\!12.26$. The ‘upper left’ and
‘lower right’ coexistence fixed points characterize the set $S$ of all such
equilibria, as in Theorem V.4. This can be seen more closely in the inset in
Figure 7, where the number beside each fixed point (in red) corresponds to the
different initial starting points (in blue) of the trajectories. Thus,
convergence to set $S$ occurs globally over the state space, but exactly which
coexistence fixed point the system converges to is dependent on the initial
point. We are thus able to observe all possible convergence scenarios from
Section V-B, including multiple coexistence equilibria.

##  VIII Concluding Remarks

By utilizing the techniques from Monotone Dynamical Systems (MDS), in this
paper, we show that a generic bi-virus epidemic model with non-linear
infection and recovery rates is monotone with respect to a specially
constructed partial ordering. This monotonicity allows us to give necessary
and sufficient conditions on the non-linear infection and recovery rates, and
thus completely characterize the entire parameter space of the bi-virus
system, a contrast to the usual Lyapunov based approach. We bridge the gap
between linear stability properties and global convergence results (or lack
thereof) for the bi-virus model with non-linear rates (including the special
case with linear rates) in the literature, and succeed in providing a complete
characterization of the trichotomy of possible outcomes for such competing
epidemics - a well known open problem. Our results demonstrate how powerful
these alternative proving techniques can be, compared to classical Lyapunov
approaches; and we note that it may be worth exploring such monotonicity
properties in other dynamics on graphs as well, where competition is a general
theme. Additionally, establishing a rigorous relationship between the SIS ODE
models with non-linear rates as studied in this paper, and the correct
probabilistic dynamics describing these non-linear rates, is of interest in
order to complete the theoretical pictures for SIS models with non-linear
rates.

## References

  * [1] V. Doshi, S. Mallick, and D. Y. Eun, “Competing Epidemics on Graphs - Global Convergence and Coexistence,” in _IEEE INFOCOM_ , 2021.
  * [2] A. Lajmanovich and J. A. Yorke, “A deterministic model for gonorrhea in a nonhomogeneous population,” _Mathematical Biosciences_ , vol. 28, no. 3, pp. 221 – 236, 1976.
  * [3] H. W. Hethcote, “The mathematics of infectious diseases,” _SIAM Review_ , vol. 42, no. 4, pp. 599–653, 2000.
  * [4] M. Garetto, W. Gong, and D. Towsley, “Modeling malware spreading dynamics,” in _IEEE INFOCOM_ , San Francisco, CA, 2003.
  * [5] L.-X. Yang, X. Yang, J. Liu, Q. Zhu, and C. Gan, “Epidemics of computer viruses: a complex-network approach,” _Applied Mathematics and Computation_ , vol. 219, no. 16, pp. 8705–8717, 2013.
  * [6] S. Hosseini and M. A. Azgomi, “A model for malware propagation in scale-free networks based on rumor spreading process,” _Computer Networks_ , vol. 108, pp. 97–107, 2016.
  * [7] K. R. Apt and E. Markakis, “Diffusion in social networks with competing products,” in _International Symposium on Algorithmic Game Theory_ , 2011\.
  * [8] B. A. Prakash, A. Beutel, R. Rosenfeld, and C. Faloutsos, “Winner takes all: competing viruses or ideas on fair-play networks,” in _ACM World Wide Web_ , 2012.
  * [9] S. F. Ruf, K. Paarporn, P. E. Pare, and M. Egerstedt, “Dynamics of opinion-dependent product spread,” in _IEEE Conference on Decision and Control_ , Melbourne, Australia, 2017.
  * [10] D. Trpevski, W. K. Tang, and L. Kocarev, “Model for rumor spreading over networks,” _Physical Review E_ , vol. 81, no. 5, p. 056102, 2010.
  * [11] L. Zhao, H. Cui, X. Qiu, X. Wang, and J. Wang, “SIR rumor spreading model in the new media age,” _Physica A: Statistical Mechanics and its Applications_ , vol. 392, no. 4, pp. 995–1003, 2013.
  * [12] X. Lin, Q. Jiao, and L. Wang, “Opinion propagation over signed networks: Models and convergence analysis,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 8, pp. 3431–3438, 2018.
  * [13] I. Koprulu, Y. Kim, and N. B. Shroff, “Battle of opinions over evolving social networks,” _IEEE/ACM Transactions on Networking_ , vol. 27, no. 2, pp. 532–545, 2019.
  * [14] S. Banerjee, A. Chatterjee, and S. Shakkottai, “Epidemic thresholds with external agents,” in _IEEE INFOCOM_ , Toronto, ON, 2014.
  * [15] A. Ganesh, L. Massoulie, and D. Towsley, “The effect of network topology on the spread of epidemics,” in _IEEE INFOCOM_ , Miami, FL, 2005.
  * [16] M. Draief and L. Massoulié, _Epidemics and Rumours in Complex Networks_ , 1st ed. Cambridge University Press, 2010\.
  * [17] F. Darabi Sahneh, C. Scoglio, and P. Van Mieghem, “Generalized epidemic mean-field model for spreading processes over multilayer complex networks,” _IEEE/ACM Transactions on Networking_ , vol. 21, no. 5, pp. 1609–1620, 2013\.
  * [18] F. D. Sahneh and C. Scoglio, “Competitive epidemic spreading over arbitrary multilayer networks,” _Physical Review E_ , vol. 89, no. 6, p. 062817, 2014\.
  * [19] A. Santos, J. M. F. Moura, and J. M. F. Xavier, “Bi-virus SIS epidemics over networks: Qualitative analysis,” _IEEE Transactions on Network Science and Engineering_ , vol. 2, no. 1, pp. 17–29, 2015.
  * [20] L.-X. Yang, X. Yang, and Y. Y. Tang, “A bi-virus competing spreading model with generic infection rates,” _IEEE Transactions on Network Science and Engineering_ , vol. 5, no. 1, pp. 2–13, 2017.
  * [21] J. Liu, P. E. Paré, A. Nedich, C. Y. Tang, C. L. Beck, and T. Basar, “Analysis and control of a continuous-time bi-virus model,” _IEEE Transactions on Automatic Control_ , 2019.
  * [22] P. Van Mieghem, “The n-intertwined SIS epidemic network model,” _Computing_ , vol. 93, no. 2–4, p. 147–169, 2011.
  * [23] J. Omic and P. Van Mieghem, “Epidemic spreading in networks—variance of the number of infected nodes,” _Delft University of Technology, Report_ , 2009\.
  * [24] P. Van Mieghem, J. Omic, and R. Kooij, “Virus spread in networks,” _IEEE/ACM Transactions on Networking_ , vol. 17, no. 1, pp. 1–14, 2009.
  * [25] A. Gray, D. Greenhalgh, L. Hu, X. Mao, and J. Pan, “A stochastic differential equation SIS epidemic model,” _SIAM Journal on Applied Mathematics_ , vol. 71, no. 3, pp. 876–902, 2011.
  * [26] C. Li, R. van de Bovenkamp, and P. Van Mieghem, “Susceptible-infected-susceptible model: A comparison of n-intertwined and heterogeneous mean-field approximations,” _Physical Review E_ , vol. 86, no. 2, p. 026116, 2012.
  * [27] Y. Wang, Z. Jin, Z. Yang, Z.-K. Zhang, T. Zhou, and G.-Q. Sun, “Global analysis of an SIS model with an infective vector on complex networks,” _Nonlinear Analysis: Real World Applications_ , vol. 13, no. 2, pp. 543–557, 2012.
  * [28] D. Guo, S. Trajanovski, R. van de Bovenkamp, H. Wang, and P. Van Mieghem, “Epidemic threshold and topological structure of susceptible-infectious-susceptible epidemics in adaptive networks,” _Physical Review E_ , vol. 88, no. 4, p. 042802, 2013.
  * [29] M. Benaïm and M. W. Hirsch, “Differential and stochastic epidemic models,” _Fields Institute communications_ , vol. 21, pp. 31–44, 1999.
  * [30] Y. Wang, G. Xiao, and J. Liu, “Dynamics of competing ideas in complex social systems,” _New Journal of Physics_ , vol. 14, no. 1, p. 013015, 2012.
  * [31] P. E. Paré, J. Liu, C. L. Beck, A. Nedić, and T. Başar, “Multi-competitive viruses over static and time-varying networks,” in _IEEE American Control Conference_ , Seattle, WA, 2017.
  * [32] A. Janson, S. Gracy, P. E. Paré, H. Sandberg, and K. H. Johansson, “Analysis of a Networked SIS Multi-Virus Model with a Shared Resource,” _IFAC-PapersOnLine_ , vol. 53, no. 5, pp. 797–802, 2020, 3rd IFAC Workshop on Cyber-Physical and Human Systems CPHS 2020.
  * [33] A. Janson, S. Gracy, P. E. Paré, H. Sandberg, and K. H. Johansson, “Networked Multi-Virus Spread with a Shared Resource: Analysis and Mitigation Strategies,” _ArXiv_ , vol. abs/2011.07569, 2020.
  * [34] P. E. Paré, J. Liu, C. L. Beck, A. Nedić, and T. Başar, “Multi-competitive viruses over time-varying networks with mutations and human awareness,” _Autom._ , vol. 123, p. 109330, 2021.
  * [35] S. Bansal, B. Grenfell, and L. Meyers, “When individual behaviour mattersl homogeneous and network models in epidemiology,” _Journal Royal Society, Interface_ , vol. 4, no. 16, pp. 879–891, 2007.
  * [36] M. E. Hochberg, “Non-linear transmission rates and the dynamics of infectious disease,” _Journal of Theoretical Biology_ , vol. 153, no. 3, pp. 301–321, 1991.
  * [37] H. Hu, K. Nigmatulina, and P. Eckhoff, “The scaling of contact rates with population density for the infectious disease models,” _Mathematical Biosciences_ , vol. 244, no. 2, pp. 125–134, 2013.
  * [38] N. D. Barlow, “Non-linear transmission and simple models for bovine tuberculosis,” _Journal of Animal Ecology_ , vol. 69, no. 4, pp. 703–713, 2000.
  * [39] C. Gan, X. Yang, W. Liu, Q. Zhu, and X. Zhang, “An epidemic model of computer viruses with vaccination and generalized nonlinear incidence rate,” _Applied Mathematics and Computation_ , vol. 222, pp. 265–274, 2013.
  * [40] W. Liu, H. Hetchote, and S. Levin, “Dynamical behavior of epidemiological models with nonlinear incidence rates.” _Journal of Mathematical Biology_ , vol. 25, no. 4, pp. 359–380, 1987.
  * [41] L.-X. Yang and X. Yang, “The impact of nonlinear infection rate on the spread of computer virus,” _Nonlinear Dynamics_ , vol. 82, 05 2015.
  * [42] H. Yuan, G. Liu, and G. Chen, “On modeling the crowding and psychological effects in network-virus prevalence with nonlinear epidemic model,” _Applied Mathematics and Computation_ , vol. 219, p. 2387–2397, 11 2012\.
  * [43] S. Ruan and W. Wang, “Dynamical behavior of an epidemic model with a nonlinear incidence rate,” _Journal of Differential Equations_ , vol. 188, no. 1, pp. 135–163, 2003.
  * [44] H. L. Smith, “Monotone dynamical systems: Reflections on new advances and applications,” _Discrete and Continuous Dynamical Systems - A_ , vol. 37, p. 485, 2017.
  * [45] P. De Leenheer and D. Aeyels, “Stability properties of equilibria of classes of cooperative systems,” _IEEE Transactions on Automatic Control_ , vol. 46, no. 12, pp. 1996–2001, 2001.
  * [46] D. Angeli and E. D. Sontag, “Monotone control systems,” _IEEE Transactions on Automatic Control_ , vol. 48, no. 10, pp. 1684–1698, 2003.
  * [47] V. S. Bokharaie, O. Mason, and M. Verwoerd, “D-stability and delay-independent stability of homogeneous cooperative systems,” _IEEE Transactions on Automatic Control_ , vol. 55, no. 12, pp. 2882–2885, 2010.
  * [48] L. Van Hien and H. Trinh, “Exponential stability of two-dimensional homogeneous monotone systems with bounded directional delays,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 8, pp. 2694–2700, 2018.
  * [49] D. Efimov, T. Raissi, and A. Zolghadri, “Control of nonlinear and lpv systems: Interval observer-based framework,” _IEEE Transactions on Automatic Control_ , vol. 58, no. 3, pp. 773–778, 2013.
  * [50] F. Forni and R. Sepulchre, “Differentially positive systems,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 2, pp. 346–359, 2016.
  * [51] M. W. Hirsch, “Systems of differential equations which are competitive or cooperative: I. limit sets,” _SIAM Journal on Mathematical Analysis_ , vol. 13, no. 2, pp. 167–179, 1982.
  * [52] C. Altafini, “Consensus problems on networks with antagonistic interactions,” _IEEE Transactions on Automatic Control_ , vol. 58, no. 4, pp. 935–946, 2013.
  * [53] M. D. Marco, M. Forti, M. Grazzini, and L. Pancioni, “Limit set dichotomy and multistability for a class of cooperative neural networks with delays,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 23, no. 9, pp. 1473–1485, 2012.
  * [54] M. Ye, B. Anderson, and J. Liu, “Convergence and equilibria analysis of a networked bivirus epidemic model,” _arXiv preprint arXiv:2111.07507_ , 2021\.
  * [55] C. D. Meyer, _Matrix analysis and applied linear algebra_. SIAM, 2000, vol. 71.
  * [56] A. Berman and R. J. Plemmons, _Nonnegative Matrices in the Mathematical Sciences_. SIAM, 1994.
  * [57] J. Yeh, _Real Analysis_ , 2nd ed. WORLD SCIENTIFIC, 2006.
  * [58] M. W. Hirsch, “Systems of differential equations that are competitive or cooperative: II. convergence almost everywhere,” _SIAM Journal on Mathematical Analysis_ , vol. 16, no. 3, pp. 423–439, 1985.
  * [59] ——, “Systems of differential equations which are competitive or cooperative: III. competing species,” _Nonlinearity_ , vol. 1, no. 1, pp. 51–71, 1988.
  * [60] ——, “System of differential equations that are competitive or cooperative: IV. structural stability in three-dimensional systems,” _SIAM Journal on Mathematical Analysis_ , vol. 21, no. 5, p. 1225–1234, 1990.
  * [61] ——, “Systems of differential equations that are competitive or cooperative: V. convergence in 3-dimensional systems,” _Journal of Differential Equations_ , vol. 80, no. 1, pp. 94 – 106, 1989.
  * [62] H. L. Smith, “Systems of ordinary differential equations which generate an order preserving flow. a survey of results,” _SIAM Review_ , vol. 30, no. 1, pp. 87–113, 1988.
  * [63] H. L. Smith and H. R. Thieme, “Quasi convergence and stability for strongly order-preserving semiflows,” _SIAM Journal on Mathematical Analysis_ , vol. 21, no. 3, pp. 673–692, 1990.
  * [64] ——, “Convergence for strongly order-preserving semiflows,” _SIAM Journal on Mathematical Analysis_ , vol. 22, no. 4, pp. 1081–1101, 1991.
  * [65] M. W. Hirsch and H. L. Smith, “Generic Quasi-convergence for Strongly Order Preserving Semiflows: A New Approach,” _Journal of Dynamics and Differential Equations_ , vol. 16, pp. 433–439, 2004.
  * [66] H. L. Smith, _Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems_. American Mathematical Society, 2014.
  * [67] ——, “Is my system of ODEs cooperative?” 2012. [Online]. Available: https://math.la.asu.edu/ halsmith/identifyMDS.pdf
  * [68] U. Krause and P. Ranft, “A limit set trichotomy for monotone nonlinear dynamical systems,” _Nonlinear Analysis: Theory, Methods & Applications_, vol. 19, no. 4, pp. 375 – 392, 1992.
  * [69] M. Ye, J. Liu, B. D. Anderson, and M. Cao, “Applications of the Poincare-Hopf Theorem: Epidemic Models and Lotka-Volterra Systems,” _IEEE Transactions on Automatic Control_ , 2021.
  * [70] V. Doshi, S. Mallick, and D. Y. Eun, “Convergence of bi-virus epidemic models with non-linear rates on networks - a monotone dynamical systems approach: Supplementary material.”
  * [71] J. Leskovec and A. Krevl, “SNAP Datasets: Stanford large network dataset collection,” http://snap.stanford.edu/data, jun 2014.
  * [72] L. Perko, _Differential Equations and Dynamical Systems_ , 3rd ed. Springer Science & Business Media, 2001.
  * [73] S. Ross, _Stochastic Processes_. Wiley, 1996.

##  Appendix A Basic Definitions and Results from Matrix Theory

We first provide some well known results surrounding irreducible square
matrices.

###### Definition A.1

[55] A square matrix $\mathbf{A}$ is reducible if there exists a permutation
matrix $\mathbf{P}$ such that $\mathbf{P}^{T}\mathbf{A}\mathbf{P}$ is a block
diagonal matrix. If no such permutation matrix exists, we say that
$\mathbf{A}$ is irreducible.

One way to check if a matrix is irreducible is by observing the underlying
directed graph, where there is an edge between two nodes only if $a_{ij}\neq
0$. The matrix $A$ is irreducible if and only if this underlying directed
graph is strongly connected.

###### Definition A.2

[56] A M-matrix is a matrix with non-positive off-diagonal elements with
eigenvalues whose real parts are non-negative.

We use the following well known result for non-negative, irreducible matrices
heavily throughout the paper.

###### Theorem A.3

(Perron-Frobenius)[55] Let $\mathbf{A}$ be a non-negative, irreducible matrix.
Then, $\lambda(\mathbf{A})$ is a strictly positive real number, and the
corresponding eigenvector $\mathbf{v}$ where
$\mathbf{A}\mathbf{v}=\lambda(\mathbf{A})\mathbf{v}$ is also strictly
positive. We call $\lambda(\mathbf{A})>0$ and $\mathbf{v}\gg\mathbf{0}$ the PF
eigenvalue and PF eigenvector of the matrix respectively.$\hfill\square$

The following result is on irreducible M-matrices.

###### Lemma A.4

[56] Given an irreducible and non-singular M-matrix $\mathbf{M}$, its inverse
$\mathbf{M}^{-1}$ has strictly positive entries.$\hfill\square$

##  Appendix B Definitions and results from ODE literature

We use the following definitions and results from the ODE literature
throughout the paper.

###### Definition B.1

The ‘flow’ of a dynamical system in a metric space $X$ is a map
$\phi:X\\!\times\\!\mathbb{R}\\!\to\\!X$ such that for any $x_{0}\\!\in\\!X$
and all $s,t\in\mathbb{R}$, we have $\phi_{0}(x_{0})\\!=\\!x_{0}$ and
$\phi_{s}\left(\phi_{t}(x_{0})\right)\\!=\\!\phi_{t+s}(x_{0})$.

###### Definition B.2

A flow $\phi:X\times\mathbb{R}\to X$ is positively invariant in set $P\subset
X$ if for every $x_{0}\in P$, $\phi_{t}(x_{0})\in P$ for all $t>0$.

###### Definition B.3

Given a flow $\phi$, an ‘equilibrium’ or a ‘fixed point’ of the system is a
point $x^{*}\in X$ such that $\\{x^{*}\\}$ is a positively invariant set. For
the ODE system $\dot{x}=F(x)$, we have $F(x^{*})=0$ at the equilibrium.

For an equilibrium point $x^{*}\in X$ we say that the trajectory starting at
$x_{0}\in X$ converges to $x^{*}$ if $\lim_{t\to\infty}\phi_{t}(x_{0})=x^{*}$.
The following result is true for stable fixed points of the ODE system from
Definition B.1.

###### Proposition B.4

[72] Let $\mathbf{J}F(x_{0})$ be the Jacobian of the ODE system evaluated at a
fixed point $x_{0}$ and assume it to be an irreducible matrix. Let
$\lambda\left(\mathbf{J}F(x_{0})\right)<0$ and suppose the corresponding
eigenvalue $\mathbf{v}$ is strictly positive ($\mathbf{v}\gg\mathbf{0}$).
Then, there exists an $\epsilon>0$ such that $F(x_{0}+r\mathbf{v})\ll 0$ for
all $r\in(0,\epsilon]$ and $F(x_{0}+r\mathbf{v})\gg 0$ for all
$r\in(0,-\epsilon]$141414In other words eigenvector $\mathbf{v}$ is tangent to
the stable manifold of the ODE system at the stable fixed points
$x_{0}$..$\hfill\square$

##  Appendix C DFR Processes as Non-Linear Recovery Rates

In this appendix, we form the connection between failure rates from
reliability theory [73], and the infection duration at any node in SIS type
epidemics. To this end, we start by formally defining the term failure rate.

###### Definition C.1

[73] Let $T>0$ be any continuous random variable with distribution
$F_{T}(s)=\mathbb{P}(T\leq s)$, and density function $f_{T}(s)$ for all $s>0$,
with $\bar{F}_{T}(s)=1-F_{T}(s)=\mathbb{P}(T>s)$. Then, the failure rate at
any given time $s>0$ is defined as

$r_{T}(s)\triangleq\frac{f_{T}(s)}{\bar{F}_{T}(s)}.$ (15)

We say $T$ has a decreasing/increasing failure rate (DFR/IFR) if $r_{T}(s)$ is
a decreasing/increasing function of $s>0$.

When $T$ is the lifetime of a system, the DFR case corresponds to the system
aging negatively. This means that as time elapses, the residual time (time
till the system fails) is more likely to increase rather than decrease. $T$
could also have an interpretation in the context of node recovery. For the
linear SIS epidemic model as in (1), consider an infected node
$i\in\mathcal{N}$ and define $T\triangleq$ time taken for node $i$ to recover
(random), with $f_{T}(s)$ and $\bar{F}_{T}(s)$ as in Definition C.1. Loosely
speaking, we can ignore the infection rate terms in (1) to take a closer look
at the recovery process via the ODE

$\dot{x}_{i}(s)=-\delta x_{i}(s),$ (16)

with the initial condition $x_{i}(0)=1$ (implying that node $i$ is last
infected at time $s=0$). The ODE (16) has an exact solution for all $s>0$,
given by $x_{i}(s)=e^{-\delta s}.$ This solution allows us to interpret
$x_{i}$ as the cumulative distribution function (CCDF) of an exponential
random variable151515When $T\sim\exp(\delta)$, we have
$\bar{F}_{T}(s)=P(T>s)=e^{-\delta s}$. with rate $\delta>0$. Using this
interpretation, we have $x_{i}(s)=P(T>s)=\bar{F}_{T}(s)$, and
$-\dot{x}_{i}(s)=f_{T}(s)$. (16) can then be rewritten as

$r_{T}(s)=\frac{-\dot{x}_{i}(s)}{x_{i}(s)}=\delta,$

for any $s>0$. $T$ is thus exponentially distributed, and has a constant
failure rate (it is both DFR and IFR).

We now consider the case where the random variable $T$ is defined for the more
general SIS epidemic model with non-linear recovery rate $q_{i}(x_{i})$ for
node $i$.161616Note that this is the special case where $q_{i}$ is only a
function of $x_{i}$, not of $x_{j}$ for neighbors $j$ of node $i$. Ignoring
the infection rate terms in (5) like before, we obtain

$\dot{x}_{i}(s)=-q_{i}\left(x_{i}(s)\right),$ (17)

retaining the previous interpretation of $x_{i}$ as the CCDF of $T$. This can
be further rearranged to obtain an expression for the failure rate as

$r_{T}(s)=\frac{-\dot{x}_{i}(s)}{x_{i}(s)}=\frac{q_{i}\left(x_{i}(s)\right)}{x_{i}(s)}$

for any $s>0$. From Definition C.1 we know $T$ is DFR if $r_{T}(s)$ is
decreasing in $s>0$. Supposing $q_{i}$ is such that $T$ is indeed DFR,
$\log(r_{T}(s))$ is also decreases in $s$, and we get

$\frac{d}{ds}\log\left(r_{T}(s)\right)=\frac{q_{i}^{\prime}(x_{i}(s))\dot{x}_{i}(s)}{q(x_{i}(s))}-\frac{\dot{x}_{i}(s)}{x_{i}(s)}\leq
0,$

where $q_{i}^{\prime}(x_{i}(s))$ denotes the derivative with respect to
$x_{i}$. Since $\dot{x}_{i}(s)=-q_{i}(x_{i}(s))$ from (17) and
$q_{i}^{\prime}(x(s))\geq 0$ from (A3), rearranging the previous equation
gives us following the condition for $T$ to be DFR

$x_{i}q_{i}^{\prime}(x_{i})-q_{i}(x_{i})\geq 0.$ (18)

In (18), the $(s)$ notation has been suppressed for clarity. Since
$q_{i}(0)=0$, the convexity of $q_{i}$ with respect to $x_{i}$ implies (18).

Roughly speaking, the DFR case (which also includes linear recovery rates as
in (1)) is a subclass of recovery rate functions $q_{i}(\mathbf{x})$
satisfying assumptions (A1)–(A5). Even though the above steps may not be
exact, they provide intuition on how infections which fester and grow worse
with time form part of our modelling assumptions in Section III.

##  Appendix D Results from MDS and Cooperative Systems

###### Definition D.1

[51, 62, 44] A flow $\phi$ is said to be monotone if for all
$\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}$ such that
$\mathbf{x}\leq_{K}\mathbf{y}$ and any $t\geq 0$, we have
$\phi_{t}(\mathbf{x})\leq_{K}\phi_{t}(\mathbf{y}).$

If the flow represents the solution of an ODE system, we say that the ODE
system is co-operative.

###### Definition D.2

Consider the system (9) and let
$\mathbf{J}F(\mathbf{x})\\!\triangleq\\!\left[{df_{i}(\mathbf{x})}/{dx_{j}}\right]$
be the Jacobian of the right hand side evaluated at any point
$\mathbf{x}\\!\in\\!\mathbb{R}^{n}$. We say that (9) is an irreducible ODE in
set $D\in\mathbb{R}^{n}$ if for all $\mathbf{x}\in D$,
$\mathbf{J}F(\mathbf{x})$ is an irreducible matrix.

###### Definition D.3

[62, 66, 44] The flow $\phi$ is said to be strongly monotone if it is
monotone, and for all $\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}$ such that
$\mathbf{x}<_{K}\mathbf{y}$, and time $t\geq 0$, we have
$\phi_{t}(\mathbf{x})\ll_{k}\phi_{t}(\mathbf{y}).$

###### Theorem D.4

[62, 66, 44] Let (9) be irreducible and co-operative in some set
$D\subset\mathbb{R}^{n}$. Then the solution $\phi$ (restricted to $t\geq 0$)
is strongly monotone.$\hfill\square$

As part of the main result of monotone dynamical systems, trajectories of
strongly monotone systems, starting from almost anywhere (in the measure
theoretic sense) in the state space, converge to the set of equilibrium points
[58, 64, 65, 44]. However, often the systems are strongly monotone only in the
interior of the state spaces instead of the entirety of the state space. In
such cases, the following results are useful.

###### Proposition D.5

(Proposition 3.2.1 in [66]) Consider the ODE system (9) which is cooperative
in a compact set $D\subset\mathbb{R}^{n}$ with respect to some cone-ordering,
and let $<_{r}$ stand for any of the order relations $\leq_{K},<_{K},\ll_{K}$.
Then,
$P_{+}\\!\triangleq\\!\left\\{\mathbf{x}\\!\in\\!D~{}|~{}\mathbf{0}\\!<_{r}\\!F(\mathbf{x})\right\\}$
and
$P_{-}\\!\triangleq\\!\left\\{\mathbf{x}\\!\in\\!D~{}|~{}F(\mathbf{x})\\!<_{r}\\!\mathbf{0}\right\\}$
are positively invariant, and the trajectory
$\left\\{\phi_{t}(\mathbf{x})\right\\}_{t\geq 0}$ for any point
$\mathbf{x}\\!\in\\!P_{+}$ or $\mathbf{x}\\!\in\\!P_{-}$ converges to an
equilibrium.$\hfill\square$

###### Theorem D.6

(Theorem 4.3.3 in [66]) Let (9) be cooperative (with respect to some cone-
ordering $\leq_{K}$) in a compact set $D\subset\mathbb{R}^{n}$ and let
$\mathbf{x}_{0}\in D$ be an equilibrium point. Suppose that
$s\triangleq\lambda(\mathbf{J}F(\mathbf{x}_{0}))>0$ (i.e. $\mathbf{x}_{0}$ is
an unstable fixed point) and there is an eigenvector
$\mathbf{v}\gg_{K}\mathbf{0}$ such that
$\mathbf{J}F(\mathbf{x}_{0})\mathbf{v}=s\mathbf{v}$. Then, there exists
$\epsilon_{0}\in(0,\epsilon]$ and another equilibrium point $\mathbf{x}_{e}$
such that for each $r\in(0,\epsilon_{0}]$, the solution
$\phi_{t}(\mathbf{x}_{r})$ has the following properties:

  * (1)

$\mathbf{x}_{r}\\!\ll_{K}\\!\phi_{t_{1}}(\mathbf{x}_{r})\\!\ll_{K}\\!\phi_{t_{2}}(\mathbf{x}_{r})\\!\ll_{K}\\!\mathbf{x}_{e}$,
for any $0\\!<\\!t_{1}\\!<\\!t_{2}$.

  * (2)

${d\phi_{t}(\mathbf{x}_{r})}/{dt}\gg_{K}\mathbf{0}$, for any $t>0$.

  * (3)

$\phi_{t}(\mathbf{x}_{r})\rightarrow\mathbf{x}_{e}$, as
$t\rightarrow\infty$.$\hfill\square$

##  Appendix E Proofs of the results in Section IV

###### Proof:

To prove that system (6) is co-operative with respect to the positive orthant,
we show that it satisfies Kamke’s condition in (10). Differentiating the right
hand side of (5) with respect to $x_{j}$, we get

$\displaystyle\frac{\partial\bar{f}_{i}(\mathbf{x})}{\partial
x_{j}}=(1-x_{i})\frac{\partial f_{i}(x)}{\partial x_{j}}=\frac{\partial
q_{i}(\mathbf{x})}{\partial x_{j}}.$

This corresponds to the $(ij)$’th off-diagonal entry of the Jacobian
$\mathbf{J}_{\bar{F}}(\mathbf{x})$ evaluated at $\mathbf{x}\in[0,1]^{N}$. It
is non-negative for any $i\neq j\in\mathcal{N}$ since $(1-x_{i})\geq 0$ and
due to assumption (A3), and the ODE (6) is therefore co-operative in
$[0,1]^{N}$ with respect to the regular cone ordering.

From assumption (A3), $\mathbf{J}_{\bar{F}}(\mathbf{x})_{ij}$ is also strictly
positive for any $\mathbf{x}\in(0,1)^{N}$ whenever $a_{ij}>0$. This means that
$\mathbf{J}_{\bar{F}}(\mathbf{x})$, and as a consequence the ODE system, is
irreducible for any $\mathbf{x}\in(0,1)^{N}$. ∎

To derive the convergence properties of the non-linear $SIS$ model, we make
use of a result form [68], rewritten below in a simpler form suitable for our
setting.

###### Theorem E.1

(Theorem 4 in [68]) Consider a generic ODE system (9) invariant to some subset
$S\subset\mathbb{R}^{N}_{+}$, and let $\mathbf{J}_{\bar{F}}$ stand for its
Jacobian matrix. Suppose that:

  * (C1)

$f_{i}(\mathbf{x})\geq 0$ for all $\mathbf{x}\geq 0$ with $x_{i}=0$;

  * (C2)

for all $\mathbf{x}\gg\mathbf{0}$ in $S$, $\alpha\in(0,1)$, it satisfies
$\mathbf{J}_{\bar{F}}(\mathbf{x})_{ij}\leq\mathbf{J}_{\bar{F}}(\alpha\mathbf{x})_{ij}$
for all $i,j\in\mathcal{N}$, with strict inequality for at least one pair of
$i,j$;

  * (C3)

for all $\mathbf{u}\ll\mathbf{w}$ in $S$, it satisfies
$\mathbf{J}_{\bar{F}}(\mathbf{w})\leq\mathbf{J}_{\bar{F}}(\mathbf{u})$;

  * (C4)

it is co-operative in $S$ with respect to the regular ordering relation, and
irreducible in $\text{Int}(S)$.

Then, exactly one of the following outcomes occurs:

  1. (i)

$\phi_{t}(\mathbf{x})$ is unbounded for all $\mathbf{x}\in
S\setminus\\{\mathbf{0}\\}$;

  2. (ii)

$\phi_{t}(\mathbf{x})\rightarrow\mathbf{0}$ as $t\rightarrow\infty$, for all
$\mathbf{x}\in S\setminus\\{\mathbf{0}\\}$;

  3. (iii)

There exists a unique, strictly positive fixed point
$\mathbf{x}^{*}\gg\mathbf{0}$ such that
$\phi_{t}(\mathbf{x})\rightarrow\mathbf{x}^{*}$ as $t\rightarrow\infty$, for
all $\mathbf{x}\in S\setminus\\{\mathbf{0}\\}$.$\hfill\square$

We now use the above to prove Theorem IV.4.

###### Proof:

We prove Theorem IV.4 by showing that it satisfies conditions (C1)-(C4) of
Theorem E.1, and then performing stability analysis to evaluate conditions for
each of the three possible outcomes therein.

From Proposition (IV.3), we know that (6) already satisfies (C4). The right
hand side of (5) satisfies (C1) because $q_{i}(x_{i})=0$ when $x_{i}=0$, and
because $(1-x_{i})$ and $f_{i}(\mathbf{x})$ are all non-negative for any
$\mathbf{x}\in[0,1]^{N}$. To check whether (C2) and (C3) is satisfied, observe
that from assumptions (A2)–(A5), we have

$\displaystyle\mathbf{J}_{F}(\mathbf{u})>\mathbf{J}_{F}(\mathbf{w})$ (19)
$\displaystyle\mathbf{J}_{Q}(\mathbf{u})<\mathbf{J}_{Q}(\mathbf{w})$ (20)

for all $\mathbf{u}<\mathbf{w}$.171717Here, the ordering between matrices
$\mathbf{M}^{a}<\mathbf{M}^{b}$ means
$\mathbf{M}^{a}_{ij}\leq\mathbf{M}^{b}_{ij}$ with the inequality being strict
for at least one pair of $i,j$. Here, $\mathbf{J}_{Q}$ is a diagonal matrix
since $\partial q_{i}/\partial x_{j}=0$ for all $i\neq j\in\mathcal{N}$.

Denote by $\mathbf{J}_{\bar{F}}$ the Jacobian matrix of system (6). Note that
for any point $\mathbf{x}\in[0,1]^{N}$, we have

$\mathbf{J}_{\bar{F}}(\mathbf{x})=\text{diag}(\mathbf{1}-\mathbf{x})\mathbf{J}_{F}(\mathbf{x})-\text{diag}\left(F(\mathbf{x})\right)-\mathbf{J}_{Q}(\mathbf{x})$
(21)

Combining the above with (19) and (20), we have for any points
$\mathbf{u}<\mathbf{w}$ that

$\displaystyle\mathbf{J}_{\bar{F}}(\mathbf{u})$
$\displaystyle=\text{diag}(\mathbf{1}-\mathbf{u})\mathbf{J}_{F}(\mathbf{u})-\text{diag}\left(F(\mathbf{u})\right)-\mathbf{J}_{Q}(\mathbf{u})$
$\displaystyle>\text{diag}(\mathbf{1}-\mathbf{w})\mathbf{J}_{F}(\mathbf{w})-\text{diag}\left(F(\mathbf{u})\right)-\mathbf{J}_{Q}(\mathbf{w})$
$\displaystyle\geq\text{diag}(\mathbf{1}-\mathbf{w})\mathbf{J}_{F}(\mathbf{w})-\text{diag}\left(F(\mathbf{w})\right)-\mathbf{J}_{Q}(\mathbf{w})$
$\displaystyle=\mathbf{J}_{\bar{F}}(\mathbf{w}),$

where the first inequality is due to
$(\mathbf{1}-\mathbf{u})>(\mathbf{1}-\mathbf{w})$ and (19) and (20). The
second inequality is from the non-negativity and monotonicity assumptions (A2)
and (A3) implying $F(\mathbf{u})\leq F(\mathbf{w})$. Since
$\mathbf{J}_{\bar{F}}(\mathbf{u})>\mathbf{J}_{\bar{F}}(\mathbf{w})$ for any
$\mathbf{u}<\mathbf{w}$, this is enough to satisfy both conditions (C2) and
(C3).

Since system (6) satisfies (C1)–(C4), Theorem E.1 applies. Since the system is
invariant in $[0,1]^{N}$, which is a bounded subset of $\mathbb{R}^{N}$,
outcome (i) of Theorem E.1 never occurs. From assumption (A1), the vector
$\mathbf{0}=[0,\cdots,0]^{T}$ (the virus-free equilibrium) is always a fixed
point of the system. We now find conditions under which trajectories of (6)
starting from anywhere in $[0,1]^{N}\setminus\\{\mathbf{0}\\}$ converge to
either zero, or to a unique strictly positive fixed point (outcomes (ii) and
(iii) in Theorem E.1 respectively), by check the stability properties of the
system.

The virus-free fixed point zero is unstable [72] when
$\lambda(\mathbf{J}_{\bar{F}}(\mathbf{0}))=\lambda(\mathbf{J}_{F}(\mathbf{0})-\mathbf{J}_{Q}(\mathbf{0}))\leq
0$. Under this condition, outcome (ii) in Theorem E.1 is not possible, and
there exists a unique, strictly positive fixed point
$\mathbf{x}^{*}\gg\mathbf{0}$ which is globally asymptotically stable in
$[0,1]^{N}\setminus\\{\mathbf{0}\\}$. Conversely when zero is a stable fixed
point, that is when
$\lambda(\mathbf{J}_{\bar{F}}(\mathbf{0}))=\lambda(\mathbf{J}_{F}(\mathbf{0})-\mathbf{J}_{Q}(\mathbf{0}))>0$,
it is globally attractive. ∎

##  Appendix F Proofs of the Main Results

Throughout this Section, we use $\phi_{t}(\mathbf{x}_{0},\mathbf{y}_{0})$ to
represent the solution of (8) at time $t\geq 0$, starting from
$(\mathbf{x}_{0},\mathbf{y}_{0})\in D$. We will need the following results to
prove the theorems from Section V-B.

###### Proposition F.1

Starting from any point $D\setminus\left\\{(\mathbf{0},\mathbf{0})\right\\}$,
trajectories of (8) converge to the set

$Z\triangleq\left\\{(\mathbf{u},\mathbf{w})\in
D~{}|~{}(\mathbf{0},\mathbf{y}^{*})\leq_{K}(\mathbf{u},\mathbf{w})\leq_{K}(\mathbf{x}^{*},\mathbf{0})\right\\}.$

###### Proof:

For any $(\mathbf{r},\mathbf{s})\in D\setminus\\{(\mathbf{0},\mathbf{0})\\}$,
there exists points $\mathbf{x},\mathbf{y}\in[0,1]^{N}$ such that
$(\mathbf{0},\mathbf{y})\leq_{K}(\mathbf{r},\mathbf{s})\leq_{K}(\mathbf{x},\mathbf{0})$.
Then, from Definition D.1 of a monotone system, we have
$\phi_{t}(\mathbf{0},\mathbf{y})\leq_{K}\phi_{t}(\mathbf{r},\mathbf{s})\leq_{K}\phi_{t}(\mathbf{x},\mathbf{0})$
for any $t>0$. Since
$\phi_{t}(\mathbf{x},\mathbf{0})\rightarrow(\mathbf{x}^{*},\mathbf{0})$ and
$\phi_{t}(\mathbf{0},\mathbf{y})\rightarrow(\mathbf{0},\mathbf{y}^{*})$, we
get
$(\mathbf{0},\mathbf{y}^{*})\leq_{K}\lim_{t\rightarrow\infty}\phi_{t}(\mathbf{r},\mathbf{s})\leq_{K}(\mathbf{x}^{*},\mathbf{0})$.
Thus the trajectory $\left\\{\phi_{t}(\mathbf{r},\mathbf{s})\right\\}_{t\geq
0}$ converges to $Z$, completing the proof. ∎

Since the set $Z$ depends on $\mathbf{x}^{*}$ and $\mathbf{y}^{*}$, the fixed
points of systems (13) and (14), and we can determine when these fixed points
are positive or zero, Proposition F.1 helps us to quickly point out a subset
of the state space to which trajectories starting from any point in
$D\\!\setminus\\!\left\\{(\mathbf{0},\mathbf{0})\right\\}$ converge.

###### Proof:

When
$\lambda\left(\mathbf{J}_{G}(\mathbf{0})-\mathbf{J}_{R}(\mathbf{0})\right)\\!\leq\\!0$
and
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})-\mathbf{J}_{S}(\mathbf{0})\right)\\!\leq\\!0$,
we know from Theorem IV.4 that $\mathbf{x}^{*}=\mathbf{y}^{*}=0$. Therefore,
trajectories of (8) starting from any point in
$D\setminus\left\\{(\mathbf{0},\mathbf{0})\right\\}$ converge to the set
$Z\triangleq\left\\{(\mathbf{u},\mathbf{w})\in
D~{}|~{}(\mathbf{0},\mathbf{0})\leq_{K}(\mathbf{u},\mathbf{w})\leq_{K}(\mathbf{0},\mathbf{0})\right\\}=\left\\{(\mathbf{0},\mathbf{0})\right\\}$.
Hence, the virus-free equilibrium is globally asymptotically stable in $D$,
which completes the proof. ∎

Proposition F.1 can also be applied to show that $(\mathbf{x}^{*},\mathbf{0})$
where $\mathbf{x}^{*}\\!\gg\\!\mathbf{0}$ is globally attractive when
$\lambda\left(\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!>\\!0$
and
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!\leq\\!0$.
This is because from Theorem IV.4, we know that
$\mathbf{x}^{*}\\!\gg\\!\mathbf{0}$ and $\mathbf{y}^{*}\\!=\\!\mathbf{0}$. We
then have $Z\triangleq\left\\{(\mathbf{u},\mathbf{w})\in
D~{}|~{}(\mathbf{0},\mathbf{0})\leq_{K}(\mathbf{u},\mathbf{w})\leq_{K}(\mathbf{x}^{*},\mathbf{0})\right\\}$,
implying that the system (8) ultimately reduces to the single $SIS$ system
(13), which we know globally converges to $\mathbf{x}^{*}$. By a symmetric
argument, we also have that $(\mathbf{0},\mathbf{y}^{*})$ where
$\mathbf{y}^{*}\\!\gg\\!\mathbf{0}$ is globally attractive when
$\lambda\left(\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!\leq\\!0$
and
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!>\\!0$.
Therefore these cases are easily analyzed by applying Proposition F.1 in
conjunction with Theorem IV.4. In terms of the linear bi-virus model whose
parameters are easier to visualize, values of $\tau_{1}$ and $\tau_{2}$ which
satisfy these conditions, lie in regions R2 and R3 of Figure 3(b) and we
henceforth exclude them from our analysis, considering only those values of
$\tau_{1}$ and $\tau_{2}$ for which $\tau_{1}\lambda(\mathbf{A})\\!>\\!1$ and
$\tau_{2}\lambda(\mathbf{B})\\!>\\!1$ always holds; equivalently considering
only the cases where
$\lambda\left(\mathbf{J}_{G}(\mathbf{0})-\mathbf{J}_{R}(\mathbf{0})\right)>0$
and
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})-\mathbf{J}_{S}(\mathbf{0})\right)>0$
always hold for nonlinear infection and recovery rates. Thus, $\mathbf{x}^{*}$
and $\mathbf{y}^{*}$ are henceforth implied to be strictly positive vectors.

Before formally proving Theorems V.3 and V.4, we provide some additional
constructions and notations which will help simplify the proofs. As in the
proof of Theorem IV.4, the Jacobians $\mathbf{J}F^{x}(\mathbf{x})$ and
$\mathbf{J}F^{y}(\mathbf{y})$ of systems (13) and (14), respectively, are

$\displaystyle\mathbf{J}F^{x}(\mathbf{x})$
$\displaystyle=\text{diag}(\mathbf{1}\\!-\\!\mathbf{x})\mathbf{J}_{G}(\mathbf{x})-\text{diag}(G(\mathbf{x}))-\mathbf{J}_{R}(\mathbf{x}),$
$\displaystyle\mathbf{J}F^{y}(\mathbf{y})$
$\displaystyle=\text{diag}(\mathbf{1}\\!-\\!\mathbf{y})\mathbf{J}_{H}(\mathbf{y})-\text{diag}(H(\mathbf{y}))-\mathbf{J}_{S}(\mathbf{y}),$

for all $\mathbf{x},\mathbf{y}\in[0,1]^{N}$. Now recall the Jacobian
$\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{x},\mathbf{y})$ of the bi-virus ODE (8)
from (12). When evaluated at $(\mathbf{x}^{*},\mathbf{0})$ and at
$(\mathbf{0},\mathbf{y}^{*})$, we get

$\begin{split}\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{x}^{*},\mathbf{0})=\begin{bmatrix}\mathbf{J}F^{x}(\mathbf{x}^{*})&\mathbf{K}\\\
\mathbf{0}&\mathbf{J}_{y}\end{bmatrix}\end{split}$ (22)

where $\mathbf{K}\\!=\\!-\text{diag}(G(vx^{*}))$,
$\mathbf{J}_{y}\\!=\\!\text{diag}(\mathbf{1}\\!-\\!\mathbf{x}^{*})\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})$,
and

$\begin{split}\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{0},\mathbf{y}^{*})=\begin{bmatrix}\mathbf{J}_{x}&\mathbf{0}\\\
\mathbf{L}&\mathbf{J}F^{y}(\mathbf{y}^{*})\end{bmatrix}\end{split}$ (23)

where $\mathbf{L}\\!=\\!-\text{diag}(H(\mathbf{y}^{*}))$,
$\mathbf{J}_{x}\\!=\\!\text{diag}(\mathbf{1}\\!-\\!\mathbf{y}^{*})\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})$.
This leads us to the following proposition, where the ordering $\leq_{K}$
($<_{K},\ll_{K}$) stands for the south east cone-ordering.

###### Proposition F.2

When
$\lambda\\!\left(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!>\\!0$,
we have
$\lambda\left(\mathbf{J}_{\bar{G}\bar{H}}(0,\mathbf{y}^{*})\right)\\!=\\!\lambda(\mathbf{J}_{x})\\!>\\!0$,
and the corresponding eigenvector
$(\mathbf{u},\mathbf{v})\\!\in\\!\mathbb{R}^{2N}$ of
$\mathbf{J}_{\bar{G}\bar{H}}(0,\mathbf{y}^{*})$ satisfies
$(\mathbf{u},\mathbf{v})\\!\gg_{K}\\!(\mathbf{0},\mathbf{0})$.

###### Proof:

First, recall that $\mathbf{y}^{*}\gg\mathbf{0}$ is the asymptotically stable
fixed point of (14). This implies that the real parts of all eigenvalues of
the Jacobian $\mathbf{J}F^{y}(\mathbf{y}^{*})$ of (14) evaluated at
$\mathbf{y}^{*}$ are negative. Since $\mathbf{J}F^{y}(\mathbf{y}^{*})$ is an
irreducible matrix as discussed in Section V-A, with non-negative off-diagonal
elements, its PF eigenvalue (obtained by perturbing with a large multiple of
the identity matrix) is real and negative, that is
$\lambda\left(\mathbf{J}F^{y}(\mathbf{y}^{*})\right)<0$.

From the assumption, we have
$\lambda\left(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(\mathbf{0})-\mathbf{J}_{R}(\mathbf{0})\right)=\lambda(\mathbf{J}_{x})>0$.
Since $\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{0},\mathbf{y}^{*})$ is a block
triangle matrix, we have
$\lambda\left(\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{0},\mathbf{y}^{*})\right)\\!=\\!\max\\!\left\\{\lambda(J_{x}),\lambda\left(\mathbf{J}F^{y}(\mathbf{y}^{*})\right)\right\\}$,
and since $\lambda\left(\mathbf{J}F^{y}(\mathbf{y}^{*})\right)<0$, we obtain
$\lambda\left(\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{0},\mathbf{y}^{*})\right)=\lambda(\mathbf{J}_{x})>0$.
Then, the corresponding eigenvector $(\mathbf{u},\mathbf{v})$ satisfies

$\mathbf{J}_{x}\mathbf{u}\\!=\\!\lambda(\mathbf{J}_{x})\mathbf{u}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}\mathbf{L}\mathbf{u}\\!+\\!\mathbf{J}F^{y}(\mathbf{y}^{*})\mathbf{v}\\!=\\!\lambda(\mathbf{J}_{x})\mathbf{v}.$

From the first equation, we can tell that $\mathbf{u}$ is the eigenvector of
$\mathbf{J}_{x}$ corresponding to its PF eigenvalue, and thus satisfies
$\mathbf{u}\\!\gg\\!\mathbf{0}$. Now recall that
$\mathbf{J}F^{y}(\mathbf{y}^{*})$ had eigenvalues with strictly negative real
parts.
$\lambda(\mathbf{J}_{x})\mathbf{I}\\!-\\!\mathbf{J}F^{y}(\mathbf{y}^{*})$ is
then a matrix with eigenvalues having strictly positive real parts (since
$\lambda(\mathbf{J}_{x})\\!>\\!0$). The matrix
$\mathbf{M}\triangleq\lambda(\mathbf{J}_{x})\mathbf{I}\\!-\\!\mathbf{J}F^{y}(\mathbf{y}^{*})$
is then, by Definition A.2, an M-matrix. By construction, it is also
irreducible and invertible and from Lemma A.4, we obtain that
$\mathbf{M}^{-1}$ is a (strictly) positive matrix. The second equation in the
above can then be rewritten as
$\mathbf{v}=\mathbf{M}^{-1}\mathbf{L}\mathbf{u}\ll\mathbf{0}$, where the
inequality is because $\mathbf{L}\\!=\\!-\text{diag}(H(\mathbf{y}^{*}))$ has
strictly negative diagonal elements ($H(\mathbf{y}^{*})$ being positive from
assumptions (A2) and (A3)). Therefore, since $\mathbf{u}\gg\mathbf{0}$ and
$\mathbf{v}\ll\mathbf{0}$, we have $(\mathbf{u},\mathbf{v})\gg_{K}\mathbf{0}$,
completing the proof. ∎

The intention behind introducing Proposition F.2 was to satisfy the
assumptions of Theorem D.6. In particular, when
$\lambda\left(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(\mathbf{0})-\mathbf{J}_{R}(\mathbf{0})\right)\\!>\\!0$,
$(0,\mathbf{y}^{*})$ is an unstable fixed point; by Proposition F.2 and
Theorem D.6, there exists an $\epsilon_{1}>0$ and another fixed point
$(\mathbf{x}_{e},\mathbf{y}_{e})$ such that for any point
$(\mathbf{x}_{r},\mathbf{y}_{r})\triangleq(\mathbf{0},\mathbf{y}^{*})+r(\mathbf{u},\mathbf{v})$
where $r\in(0,\epsilon_{1}]$, we have

$\begin{split}(0,\mathbf{y}^{*})\\!\ll\\!(\mathbf{x}_{r},\mathbf{y}_{r})\\!\ll_{K}\\!\phi_{t}(\mathbf{x}_{r},\mathbf{y}_{r})\\!\ll_{K}\\!\phi_{s}(\mathbf{x}_{r},\mathbf{y}_{r})\\!\leq_{K}\\!(\mathbf{x}^{*},\mathbf{0})\end{split}$

for all $s\\!>\\!t\\!>\\!0$. Moreover, for all $(\mathbf{x},\mathbf{y})$ such
that
$(\mathbf{0},\mathbf{y}^{*})\\!\ll_{K}\\!(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\mathbf{x}_{e},\mathbf{y}_{e})$,
there exists an $r\\!\in\\!(0,\epsilon]$ sufficiently small such that
$(\mathbf{x}_{r},\mathbf{y}_{r})\\!\leq_{K}\\!(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\mathbf{x}_{e},\mathbf{y}_{e})$.
Since
$\phi_{t}(\mathbf{x}_{r},\mathbf{y}_{r})\\!\rightarrow\\!(\mathbf{x}_{e},\mathbf{y}_{e})$,
monotonicity implies
$\phi_{t}(\mathbf{x},\mathbf{y})\\!\rightarrow\\!(\mathbf{x}_{e},\mathbf{y}_{e})$
as $t\\!\to\\!\infty$.

Now, we can either have
$(\mathbf{x}_{e},\mathbf{y}_{e})\\!=\\!(\mathbf{x}^{*},\mathbf{0})$, which
occurs when $(\mathbf{x}^{*},\mathbf{0})$ is the other stable fixed point of
(8), or
$(\mathbf{x}_{e},\mathbf{y}_{e})\\!=\\!(\hat{\mathbf{x}},\hat{\mathbf{y}})\\!\gg\\!\mathbf{0}$
which occurs when $(\mathbf{x}^{*},\mathbf{0})$ is an unstable fixed point.
Note that $(\mathbf{x}^{*},\mathbf{0})$ is stable (unstable) if and only if
$\lambda\left(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{J}_{G}(\mathbf{0})-\mathbf{J}_{R}(\mathbf{0})\right)\\!\leq\\!0$
($>\\!0$). We will talk about both these possibilities one by one and
exploring these will eventually lead to Theorems V.3 and V.4. But before we do
that, we first prove the following proposition about convergence to the fixed
point $(\mathbf{x}_{e},\mathbf{y}_{e})$ (whichever of the two it may be).

###### Proposition F.3

Trajectories of the system (4) starting from any point
$(\mathbf{x},\mathbf{y})$ such that
$(\mathbf{0},\mathbf{y}^{*})<_{K}(\mathbf{x},\mathbf{y})\leq_{K}(\mathbf{x}_{e},\mathbf{y}_{e})$
converge to $(\mathbf{x}_{e},\mathbf{y}_{e})$.$\hfill\square$

###### Proof:

Recall that we already know that for all
$(\mathbf{0},\mathbf{y}^{*})\\!\ll_{K}\\!(\mathbf{x},\mathbf{y})\\!\leq\\!(\mathbf{x}^{*},\mathbf{y}^{*})$,
$\phi_{t}(\mathbf{x},\mathbf{y})\\!\rightarrow\\!(\mathbf{x}^{*},\mathbf{0})$.
We would however like to show this for all $(\mathbf{x},\mathbf{y})\in
Z\setminus(\mathbf{0},\mathbf{y}^{*})$, that is even when
$(\mathbf{x},\mathbf{y})$ satisfies
$(\mathbf{0},\mathbf{y}^{*})<_{K}(\mathbf{x},\mathbf{y})\leq(\mathbf{x}^{*},\mathbf{0})$.
To do this, we create a set of points which converge to
$(\mathbf{x}^{*},\mathbf{0})$, just like we created
$(\mathbf{x}_{r},\mathbf{y}_{r})$ before, and then use a monotonicity argument
to show convergence to $(\mathbf{0},\mathbf{y}^{*})$ of trajectories starting
for all points $(\mathbf{x},\mathbf{y})$ satisfying
$(\mathbf{0},\mathbf{y}^{*})<_{K}(\mathbf{x},\mathbf{y})\leq(\mathbf{x}^{*},\mathbf{0})$.

Recall that, $\mathbf{y}^{*}$ is an asymptotically stable fixed point of (14),
and from the proof of Proposition F.2 we know that
$\lambda\left(\mathbf{J}F^{y}(\mathbf{y}^{*})\right)<0$. Let
$\mathbf{w}\gg\mathbf{0}$ be the corresponding PF eigenvector. Then by
Proposition B.4, there exists an $\epsilon_{2}>0$ such that for all
$s\in(0,\epsilon_{2}]$, $F^{y}(\mathbf{y}^{*}+s\mathbf{w})\ll\mathbf{0}$. We
can then define points
$(\mathbf{x}_{r},\mathbf{y}_{s})\triangleq(r\mathbf{u},\mathbf{y}^{*}+s\mathbf{w})$
for any $r\in(0,\epsilon_{1}]$ and $s\in(0,\epsilon_{2}]$, where
$\mathbf{u}\gg\mathbf{0}$ is the eigenvector of $\mathbf{J}_{x}$ from
Proposition F.2. We will first show that trajectories starting from these
points converge to $(\mathbf{x}_{e},\mathbf{y}_{e})$. By rearranging the terms
of (8), we can rewrite it as

$\displaystyle\dot{\mathbf{x}}=$
$\displaystyle~{}\text{diag}(\mathbf{1}-\mathbf{y}^{*})G(\mathbf{x})-R(\mathbf{x})+\text{diag}(\mathbf{y}^{*}-\mathbf{x}-\mathbf{y})G(\mathbf{x})$
$\displaystyle=$
$\displaystyle~{}\text{diag}(\mathbf{1}-\mathbf{y}^{*})\mathbf{J}_{G}(\mathbf{0})\mathbf{x}-\mathbf{J}_{R}(\mathbf{0})\mathbf{x}$
$\displaystyle+\text{diag}(\mathbf{y}^{*}-\mathbf{x}-\mathbf{y})G(\mathbf{x})+O\left(\|\mathbf{x}\|^{2}\right)$
$\displaystyle=$
$\displaystyle~{}\mathbf{J}_{x}\mathbf{x}+O\left(\|\mathbf{x}\|\left[\|\mathbf{y}-\mathbf{y}^{*}\|+\|\mathbf{x}\|\right]\right),$
$\displaystyle\dot{\mathbf{y}}=$
$\displaystyle~{}\text{diag}(\mathbf{1}-\mathbf{y})H(\mathbf{y})-S(\mathbf{y})-\text{diag}(\mathbf{x})H(\mathbf{y})$
$\displaystyle=$
$\displaystyle~{}F^{y}(\mathbf{y})+O\left(\|\mathbf{y}\|\right),$

for all $(\mathbf{x},\mathbf{y})\in D,$181818Here, $O(x)$ is used to represent
terms which satisfy $O(x)\to 0$ as $x\to 0$.where the first equality is from a
Taylor series expansion of $G$ and $R$ around $\mathbf{0}$. For any point
$(\mathbf{x}_{r},\mathbf{y}_{s})=(r\mathbf{u},\mathbf{y}^{*}+s\mathbf{w})$,
the above equations can be written as

$\displaystyle\dot{\mathbf{x}}$
$\displaystyle=r\lambda(\mathbf{J}_{x})\mathbf{u}+rO\left(\|\mathbf{u}\|\left[s\|\mathbf{w}\|+r\|\mathbf{u}\|\right]\right)$
$\displaystyle=r\left[\lambda(\mathbf{J}_{x})\mathbf{u}+O(r+s)\right]\ $
$\displaystyle\dot{\mathbf{y}}$
$\displaystyle=F^{y}(\mathbf{y}^{*}+s\mathbf{w})+O\left(\|s\mathbf{y}\|\right)$
$\displaystyle=F^{y}(\mathbf{y}^{*}+s\mathbf{w})+O\left(s\right).$

For sufficiently small $r$ and $s$, we have $\dot{\mathbf{x}}\gg\mathbf{0}$
(since $\lambda(\mathbf{J}_{x})>0$ and $\mathbf{u}\gg\mathbf{0}$) and
$\dot{\mathbf{y}}\ll\mathbf{0}$ (since
$F^{y}(\mathbf{y}^{*}+s\mathbf{w})\ll\mathbf{0}$ for all
$s\in(0,\epsilon_{2}]$). This satisfies the conditions for Proposition D.5,
and trajectories starting from such points will be monotonically increasing
(according to the south-east cone ordering), eventually converging to the
fixed point $(\mathbf{x}_{e},\mathbf{y}_{e})$.

Now see that for any point $(\mathbf{x},\mathbf{y})$ such that
$(\mathbf{0},\mathbf{y}^{*})\\!<_{K}\\!(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\mathbf{x}_{e},\mathbf{y}_{e})$,
where $\mathbf{x}\\!>\\!\mathbf{0}$ and $\mathbf{y}\\!\leq\\!\mathbf{y}^{*}$,
by the nature of the ODE system (4) all zero entries of the $\mathbf{x}$ term
will eventually become positive (if it isn’t already). Therefore, there exists
a time $t_{1}>0$ such that $\mathbf{x}(t_{1})\\!\gg\\!\mathbf{0}$, and there
exist $r,s$ small enough such that
$(\mathbf{x}_{r},\mathbf{y}_{s})\\!\ll_{K}\\!\phi_{t_{1}}(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\mathbf{x}_{e},\mathbf{y}_{e})$.
Again by monotonicity, since
$\phi_{t}(\mathbf{x}_{r},\mathbf{y}_{s})\rightarrow(\mathbf{x}_{e},\mathbf{y}_{e})$,
we have
$\phi_{t+t_{1}}(\mathbf{x},\mathbf{y})\rightarrow(\mathbf{x}_{e},\mathbf{y}_{e})$
as $t\rightarrow\infty$, completing the proof. ∎

We now consider the case where
$(\mathbf{x}_{e},\mathbf{y}_{e})\\!=\\!(\mathbf{x}^{*},\mathbf{0})$ and give
the proof for Theorem V.3. We prove it only for when
$\tau_{1}\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})\\!>\\!1$ and
$\tau_{2}\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})\\!\leq\\!1$, since the
other case follows by a symmetric argument.

###### Proof:

When
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!\leq\\!0$,
$(\mathbf{x}^{*},\mathbf{0})$ is a stable fixed point of system (8), since all
eigenvalues of $\mathbf{J}_{\bar{G}\bar{H}}(\mathbf{x}^{*},\mathbf{0})$ have
non-positive real parts, and we have
$(\mathbf{x}_{e},\mathbf{y}_{e})=(\mathbf{x}^{*},\mathbf{0})$. Proposition F.3
then implies that trajectories starting from all points in
$Z\setminus\left\\{(\mathbf{0},\mathbf{y}^{*})\right\\}$ converge to
$(\mathbf{x}^{*},\mathbf{0})$. According to Proposition F.1, trajectories
starting from all points $(\mathbf{x},\mathbf{y})\in B_{x}$ in the system
eventually enter the set $Z$, thereby eventually converging to
$(\mathbf{x}^{*},\mathbf{0})$, giving us global convergence in $B_{x}$. ∎

Similarly, we use Propositon F.3 to prove Theorem V.4.

###### Proof:

When
$\lambda\left(\mathbf{J}_{G}(\mathbf{0})\\!-\\!\mathbf{J}_{R}(\mathbf{0})\right)\\!>\\!0$
and
$\lambda\left(\mathbf{J}_{H}(\mathbf{0})\\!-\\!\mathbf{J}_{S}(\mathbf{0})\right)\\!>\\!0$,
both $(\mathbf{0},\mathbf{y}^{*})$ and $(\mathbf{x}^{*},\mathbf{0})$ are
unstable fixed points, and $(\mathbf{x}_{e},\mathbf{y}_{e})$ takes the form of
a positive fixed point $(\hat{\mathbf{x}},\hat{\mathbf{y}})\gg\mathbf{0}$ (it
cannot be $(\mathbf{x}^{*},\mathbf{0})$, which is unstable). Then from
Proposition F.3, it attracts trajectories beginning from all points
$(\mathbf{x},\mathbf{y})$ satisfying
$(\mathbf{0},\mathbf{y}^{*})\\!<_{K}\\!(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\hat{\mathbf{x}},\hat{\mathbf{y}})$.

Similarly, we have a symmetric result beginning from
$\tau_{2}\lambda(\mathbf{S}_{\mathbf{x}^{*}}\mathbf{B})\\!>\\!1$ (symmetric to
Proposition F.2 which assumes
$\tau_{1}\lambda(\mathbf{S}_{\mathbf{y}^{*}}\mathbf{A})\\!>\\!1$ instead), and
we can say that there exists another fixed point
$(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\gg\\!\mathbf{0}$ which attracts all
points $(\mathbf{x},\mathbf{y})$ satisfying
$(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\leq_{K}\\!(\mathbf{x},\mathbf{y})\\!<_{K}\\!(\mathbf{x}^{*},\mathbf{0})$.
By construction, we then have
$(\hat{\mathbf{x}},\hat{\mathbf{y}})\\!\leq_{K}\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})$,
with the possibility of being equal.

To prove global convergence of the system to the set
$S\\!=\\!\left\\{(\mathbf{x}_{e},\mathbf{y}_{e})\\!\in\\!E~{}|~{}(\hat{\mathbf{x}},\hat{\mathbf{x}})\\!\leq_{K}\\!(\mathbf{x}_{e},\mathbf{y}_{e})\\!\leq_{K}\\!(\bar{\mathbf{x}},\bar{\mathbf{y}})\right\\}$,
observe first that as part of the proof of Proposition F.3 we showed that for
trajectories starting from any point $(\mathbf{x},\mathbf{y})$ in the state
space, there exists $r\\!>\\!0$ and $s\\!>\\!0$ small enough, and
$t_{1}\\!>\\!0$ such that
$(\mathbf{x}_{r},\mathbf{y}_{s})\\!\ll_{K}\\!\phi_{t_{1}}(\mathbf{x},\mathbf{y})\\!\leq_{K}\\!(\hat{\mathbf{x}},\hat{\mathbf{y}})$
where $(\mathbf{x}_{r},\mathbf{y}_{s})$ is a point very close to
$(\mathbf{x}^{*},\mathbf{0})$. By a parallel argument, we can find a similar
point $(\mathbf{x}_{p},\mathbf{y}_{q})$ very close to
$(\mathbf{0},\mathbf{y}^{*})$ and a time $t_{2}$ such that
$(\bar{\mathbf{x}},\bar{\mathbf{y}})\\!\leq_{K}\\!\phi_{t_{2}}(\mathbf{x},\mathbf{y})\\!\ll_{K}\\!(\mathbf{x}_{p},\mathbf{y}_{q})$.
Then, we have
$(\mathbf{x}_{r},\mathbf{y}_{s})\\!\ll_{K}\\!\phi_{\max\\{t_{1},t_{2}\\}}(\mathbf{x},\mathbf{y})\\!\ll_{K}\\!(\mathbf{x}_{p},\mathbf{y}_{q})$.
Since
$\phi_{t}(\mathbf{x}_{r},\mathbf{y}_{s})\\!\rightarrow\\!(\hat{\mathbf{x}},\hat{\mathbf{x}})\in
S$, and
$\phi_{t}(\mathbf{x}_{p},\mathbf{y}_{q})\rightarrow(\bar{\mathbf{x}},\bar{\mathbf{x}})\in
S$, we can once again, due to monotonicity of the system and by invoking a
sandwich argument, say that
$\phi_{t+\max\\{t_{1},t_{2}\\}}(\mathbf{x},\mathbf{y})$ converges to an
equilibrium point in $S$ as $t\\!\rightarrow\\!\infty$. This completes the
proof. ∎

|  Vishwaraj Doshi received his B.E. degree in mechanical engineering from the
University of Mumbai, Mumbai, MH, India, and Masters degree in Operations
Research from North Carolina State University, Raleigh, NC, USA, in 2015 and
2017 respectively. He completed his Ph.D. degree with the Operations Research
Graduate Program at North Carolina State University in 2022, and is now a part
of the Data Science and Advanced Analytics team at IQVIA. His primary research
interests include design of randomized algorithms on graphs, and epidemic
models on networks.
---|---
|  Shailaja Mallick is a Ph.D. student in the Computer Science Department at
North Carolina State University. She received her B.Tech in Computer Science
from UCE, Burla, India and Masters in Computer Systems and Networks from
Chalmers University of Technology, Sweden. Her current research interests are
in the area of social network analysis, network and performance modeling using
techniques from mathematical biology, graph theory, stochastic modeling and
simulation.
---|---
|  Do Young Eun (Senior Member, IEEE) received his B.S. and M.S. degree in
Electrical Engineering from Korea Advanced Institute of Science and Technology
(KAIST), Taejon, Korea, in 1995 and 1997, respectively, and Ph.D. degree from
Purdue University, West Lafayette, IN, in 2003. Since August 2003, he has been
with the Department of Electrical and Computer Engineering at North Carolina
State University, Raleigh, NC, where he is currently a professor. His research
interests include distributed optimization for machine learning, machine
learning algorithms for networks, distributed and randomized algorithms for
large social networks and wireless networks, epidemic modeling and analysis,
graph analytics and mining techniques with network applications. He has been a
member of Technical Program Committee of various conferences including IEEE
INFOCOM, ICC, Globecom, ACM MobiHoc, and ACM Sigmetrics. He is serving on the
editorial board of IEEE Transactions on Network Science and Engineering, and
previously served for IEEE/ACM Transactions on Networking and Computer
Communications Journal, and was TPC co-chair of WASA’11. He received the Best
Paper Awards in the IEEE ICCCN 2005, IEEE IPCCC 2006, and IEEE NetSciCom 2015,
and the National Science Foundation CAREER Award 2006. He supervised and co-
authored a paper that received the Best Student Paper Award in ACM MobiCom
2007.
---|---