In February 1918 a deadly influenza pandemic (popularly known as the Spanish flu) swept across the globe. It lasted until 1920, and caused approximately 50 million deaths [1] . Influenza viruses have continued to spread across the globe in recurring epidemics [2] . Given that the spread of infectious diseases has an enormous impact on society, the study of spread has been an active area of research since Bernoulli's seminal paper [3] . The overarching goal of these research directions is to find conditions that would cause an epidemic to become eradicated, and leverage the knowledge of these conditions to design spread control strategies. To this end, various infection models have been proposed and studied in the literature; susceptible-infected-susceptible (SIS), susceptible-infectedremoved (SIR), susceptible-exposed-infected-removed (SEIR), etc. In this paper, we focus on SIS models. Axel More specifically, we consider networked SIS models, in which a population of individuals is divided into subpopulations, and the infection can spread both within and between these subpopulations. Networked SIS models have been studied extensively using discrete-time [4] - [7] and continuoustime dynamics [8] - [10] . In the present paper, we will focus on continuous-time dynamics. Although both time-invariant and time-varying continuous-time models have been studied in the literature (see for instance [11] , [12] and references therein), we will restrict ourselves to the time-invariant case.
The spread represented by the SIS model has typically been understood as a consequence of human contact. However, the spread of infectious diseases can significantly worsen due to the presence of a shared resource. For example, waterborne pathogens spread via water distribution systems [13] and droplet-transmitted pathogens spread via surfaces in public transit vehicles [14] . Such observations have motivated the development of the susceptible-infected-watersusceptible (SIWS) model [15] , [16] -essentially, a networked continuous-time single-virus SIS model that incorporates shared resources. For the SIWS model with a single shared resource, sufficient conditions for asymptotic convergence to the healthy state (where each subpopulation is infection-free, and the shared resource is contamination-free) have been provided [15, Theorem 1] . This result has been generalized to also account for multiple shared resources in [16, Theorem 1] , where virus transmission could be due to individual-to-individual contact, due to individual-to-resource contact, and, in contrast to the single resource case, also due to resource-to-resource contact. However, no theoretical guarantees for the endemic behavior (i.e., where the virus persists) of the SIWS model have been provided.
The analysis of epidemic spread using SIWS models has been restricted to the single-virus case. While such analysis provides insights into how to battle an epidemic, they do not account for settings where multiple strains of a virus could be simultaneously active within a population. In particular, it is possible that viral strains compete with each other to infect the population: each individual can be infected by one, and only one, of the multiple viral strains prevalent [17] , i.e., the strains impose cross-immunity. Such a phenomenon is not only restricted to epidemics, but may also appear in the context of incompatible information transmission on social networks [18] , or of two competing products within the same market [19] . It is well-known that competitive multi-virus propagation exhibits rich behaviour in comparison to singlevirus propagation [20] . One possible outcome of competitive multi-virus propagation is coexistence (i.e., multiple strains coexist in a population by infecting separate fractions of each subpopulation), while another is competitive exclusion (i.e., the spread parameters of one strain dominate those of the other strains, thereby causing those strains to become eradicated). In particular, [21] established that cross-immunity between strains typically leads to competitive exclusion in singlepopulation SIS models, whereas [18] , [22] - [24] provided conditions for coexistence in networked SIS models. In particular, analysis of the various equilibria of a competing continuoustime time-invariant bi-virus model has been provided in [23] , whereas a necessary and sufficient condition for a coexisting equilibrium has been established [24, Theorem 6] . However, the results obtained in [23] , [24] are restrictive in the following sense: i) [23, Theorems 6 and 7] rely on the assumption that the spread parameters with respect to each virus is the same for every subpopulation; and ii) [24, Theorem 6] is reliant on the assumption that the set of spread parameters for each virus is a scaled version of that of other viruses. Moreover, none of these works account for the presence of a shared resource in the network. Thus, to the best of our knowledge, for the multivirus SIWS model, a detailed analysis of the various equilibria and their stability properties remains missing in the literature. The present paper aims to fill this gap and in the process, for bi-virus SIS models, establish conditions for the existence of a coexisting equilibrium under less stringent assumptions.
The overarching goal is to develop strategies for the eradication of epidemics. Several control strategies have been formulated in the context of networked SIS models: optimal antidote distribution [25] , contact reduction [26] , [27] , etc. In particular, for a directed network with heterogeneous spread, [28] provides a fully distributed Alternating Direction Method of Multipliers (ADMM) algorithm that allows for local computation of optimal investment required to boost the healing rate at each node. Note that this ADMM algorithm requires every node to communicate with its neighbors its local estimate of the full network, which can be undesirable to do in practice. A decentralized algorithm that involves disconnecting nodes and increasing the healing rates, subject to resource constraints, has been proposed in [29] . More recently, a distributed algorithm that, given resource limitations, guarantees the eradication of an epidemic with a specified rate has been proposed in [25] . In the case of multiple cross-immune strains, eradication of all viruses via optimal antidote distribution is featured in [30] . In certain multi-strain epidemics, one strain can impose more severe symptoms than other strains [31] . If resource limitations do not permit us to eradicate all strains, we might prioritize eradicating a malignant strain first. The idea of eradicating one strain while sustaining another has been explored in [32] . Note that the algorithms provided in [23] , [25] , [28] , [29] , and the results in [24] , [32] - [34] account only for SIS models.
The paper studies the multi-virus SIWS model. We establish existence, uniqueness and stability of certain equilibria, and develop mitigation strategies based on these results, as follows:
(i) Since the persistence of a virus has an enormous impact on society, it is important to know under what conditions a virus persists and, assuming at least one individual in a subpopulation is infected initially, if the virus can become endemic in the population. To this end, we provide a sufficient condition for the existence, uniqueness, and asymptotic stability of a single-virus endemic equilibrium; see Theorem 3. (ii) In the bi-virus case with (resp. without) a shared resource, as discussed previously, a natural question is whether both viruses can persist simultaneously. We establish a sufficient condition for the existence of a coexisting equilibrium; see Theorem 5 (resp. Corollary 1). (iii) Similarly, in the bi-virus case with (resp. without) a shared resource, it is also important to determine conditions that lead to competitive exclusion. Thus, we establish a necessary condition for the existence of a coexisting equilibrium; see Theorem 6 (resp. Corollary 2). (iv) Utilizing our analytical results, we have developed two mitigation strategies: one for the stabilization (exponential or asymptotic) of the healthy state, and the other for the eradication of a malignant virus in the bi-virus case; see Theorems 7 and 8, respectively. Additionally, we present results on the stability and uniqueness of the healthy state, see Theorems 1, 2 and 4. A preliminary version of some of the results in this paper is scheduled for presentation at the 3 rd IFAC Workshop on Cyber-Physical & Human Systems; see [35] .
The aforementioned findings have several consequences. In an epidemiological context, the results on the healthy state shed light on the conditions that ensure the eradication of a virus in the population and the shared resource, and, therefore, aid in developing a mitigation strategy (see Section VII). The results on the single-virus endemic equilibrium are useful for learning the spread parameters (i.e., healing rate and infection rate), as discussed for the special case of SIS models in [36] . One of the results concerning the coexisting equilibrium inspires a novel virus eradication strategy, where the goal is to push the dynamics towards the single-virus endemic equilibrium of the benign virus, as opposed to the healthy state (see Section VII). The results on the existence of a coexisting equilibrium find applications not just in epidemiology, but also in opinion dynamics in social networks, where it provides conditions under which competing ideas (for instance, pro-Covid-19 restrictions versus anti-Covid-19 restrictions) coexist in different groups of the same population.
The paper is organized as follows: We introduce the model and formally present the main problems being investigated in Section II. We collect the necessary technical background needed for establishing our main results in Section III. The main results are split up across Sections IV-VII. More specifically, Sections IV-VI concern the analysis of equilibria, whereas mitigation strategies for eradicating or mitigating the spread of viruses are provided in Section VII. We illustrate our findings via simulations in Section VIII. Finally, we summarize the main findings of our paper, and highlight relevant problems of future interest in Section IX. Fig. 1 . Visualization of the model for the case when m = 2. An individual is either susceptible (S), infected with virus 1 (I 1 ), or infected with virus 2 (I 2 ). The shared resource (W) is contaminated by individuals infected with either virus, and in turn augments the corresponding infection rate.
We denote the set of real numbers by R, and the set of nonnegative real numbers by R + . For any positive integer n, we use [n] to denote the set {1, 2, ..., n}. The i th entry of a vector x is denoted by x i . The element in the i th row and j th column of a matrix M is denoted by M ij . We use 0 and 1 to denote the vectors whose entries all equal 0 and 1, respectively, and use I to denote the identity matrix, while the sizes of the vectors and matrices are to be understood from the context. For a vector x we denote the square matrix with x along the diagonal by diag(x). For any two real vectors
For a square matrix M , we use σ(M ) to denote the spectrum of M , ρ(M ) to denote the spectral radius of M , and s(M ) to denote the largest real part among the eigenvalues of M , i.e., s(M ) = max{Re(λ) : λ ∈ σ(M)}. Given a matrix A, A ≺ 0 (resp. A 0) indicates that A is negative definite (resp. negative semidefinite), whereas A 0 (resp. A 0) indicates that A is positive definite (resp. positive semidefinite). We denote a subset by P ⊆ Q, a proper subset by P ⊂ Q, and set difference by P \ Q.
In this section, we detail a model of multi-viral spread across a population network with a shared resource. We then formally state the problems being investigated. Finally, we introduce pertinent assumptions and definitions for use in later sections.
Consider a population of individuals, subdivided into n population nodes in a network, with a resource W being shared among some or all of the population nodes. Suppose that m viruses are active in the population. An individual can become infected by a virus, either by coming into contact with an infected individual, or due to interaction with the (possibly) contaminated shared resource. We make the assumption that viral infection causes cross-immunity, meaning that an individual can be infected by no more than one virus at the same time. An infected individual can then recover, returning to the susceptible state. The model is visualized from an individual's perspective in Figure 1 . The spread of the m viruses across the population can be represented by a multi-layer network G with m layers, where the vertices correspond to population nodes and the shared resource, and each layer contains a set of directed edges, E k , specific to each virus k. There exists a directed edge from node j to node i in E k if an individual, infected by virus k in node j, can directly infect individuals in node i. Furthermore, the existence of a directed edge in E k from node j (resp. from the shared resource W ) to the shared resource W (resp. to the node j) signifies that the shared resource W (resp. node j) can be contaminated with virus k by infected individuals in node j (resp. by the shared resource W ). We say that the k th layer of G is strongly connected if there is a path via the directed edges in E k from each node, and from the shared resource, to every node, and to the shared resource. In real-world scenarios, it is often the case that viral epidemics can spread from each subpopulation to every other subpopulation, in which case we assume that each layer is strongly connected.
Each population node i contains N i individuals, with a birth rate µ i equal to its death rateμ i . At any time t ≥ 0, S i (t) is the number of susceptible individuals in node i, while I k i (t) is the number of individuals infected by virus k in node i,
The rate at which individuals infected by virus k in node j infect susceptible individuals in node i is denoted by α k ij , where α k ij = 0 corresponds to the absence of a directed edge from node j to node i in E k . In node i, individuals infected by virus k recover to the susceptible state at a rate γ k i . The shared resource contains a viral mass with respect to each virus k, denoted by W k (t), representing the level of contamination at time t ≥ 0. The viral mass of virus k grows at a rate proportional to all I k i (t) scaled by their corresponding rates ζ k i , and decays at a rate δ k w . The resource-to-node infection rate to node i with respect to virus k is denoted by α k iw . The time evolution of the number of susceptible and infected individuals (with respect to each virus k ∈ [m]) in population node i ∈ [n] is given bẏ
We define new variables to simplify the system. Let
where the variables can be interpreted as follows: With respect to virus k, p k i (t) is the fraction of currently infected individuals in node i, z k (t) is a scaled contamination level in the shared resource, δ k i is the healing rate in node i, β k ij is the node-tonode infection rate from node j to i, scaled with respect to population ratios, β k iw is a scaled resource-to-node infection rate to node i, and c k i is a scaled node-to-resource contamination rate from node i. Then, assuming that the birth rates and death rates are equal for each node, (1) can be rewritten aṡ
Using vector notation, (2) can be rewritten aṡ
where B k is an n × n-matrix with β k ij as the (i, j) th element, D k is a diagonal n × n-matrix with D k ii = δ k i for all i ∈ [n], b k is a column vector with β k iw as the i th element, and c k is a row vector with c k i as the i th element. To simplify notation further, we define
With these variables in place, we can rewrite (3) aṡ
Defining A k w (y(t)) := − D k w + (I − m l=1 X(y l (t)))B k w , the dynamics of the system of all m viruses are given bẏ
We say that virus k ∈ [m] is eradicated, or in its eradicated state, if y k = 0, which is clearly an equilibrium of (4). When considering the system of m viruses (5), we say that the system is in the healthy state if all viruses are eradicated, i.e., y = 0. If (5) has an endemic (non-zero) equilibrium, it can belong to one of two types: single-virus endemic equilibrium, where y k > 0 for some k ∈ 
For the model (5), we formally state the problems being investigated in this paper.
(i) Under what conditions does y k (t) converge exponentially, or asymptotically, to its eradicated state, i.e., y k = 0, for some k ∈ [m]? (ii) For a single-virus setup, i.e., m = 1, under what conditions does the system have a unique single-virus endemic equilibrium, y * > 0, and under such conditions, does the system converge asymptotically to y * from any non-zero initial condition? (iii) What is a necessary and sufficient condition for the healthy state, i.e., y = 0, to be the unique equilibrium? (iv) For a bi-virus setup, i.e., m = 2, under what conditions does the system have a coexisting equilibrium, i.e., (ŷ 1 ,ŷ 2 ) such thatŷ 1 > 0 andŷ 2 > 0? (v) For a bi-virus setup, under what conditions does the system not have a coexisting equilibrium, i.e., (ŷ 1 ,ŷ 2 ) such thatŷ 1 > 0 andŷ 2 > 0? (vi) How can the healing rates, i.e., δ k i , be chosen to ensure that the system converges exponentially, or asymptotically, to the healthy state, i.e., y = 0? (vii) For a bi-virus setup, how can the healing rates of virus 2,
i.e., δ 2 i , be chosen to ensure that the system converges to the single-virus endemic equilibrium of virus 1? Before we address these questions, we point out two connections between the considered setup and the existing literature.
If m = 1, system (4) coincides with the singlevirus model proposed in [15] Hereafter, when a single-virus system is considered, we drop the superscripts identifying the virus from all variables:
, the influence of the shared resource on the population is nullified. Then, the multi-virus dynamics of the population nodes, p k (t), in (5) are equivalent to the time-invariant multi-virus setup of [30] .
In order for (5) to be well-defined and realistic, we make the following assumption.
and k ∈ [m], with c k l > 0 for at least one l ∈ [n].
Note that if Assumption 1 holds, then for all k ∈ [m], B k w is a nonnegative matrix and D k w is a positive diagonal matrix Moreover, recall that a square matrix M is said to be irreducible if, replacing the non-zero elements of M with ones and interpreting it as an adjacency matrix, the corresponding graph is strongly connected. Then, noting that non-zero elements in B k w represent directed edges in the set E k , we see that B k w is irreducible whenever the k th layer of the multi-layer network G is strongly connected.
Thanks to Assumption 1, we can restrict our analysis to the sets D :
is to be interpreted as a fraction of a population, and z k (t) is a nonnegative quantity, these sets represent the sensible domain of the system. That is, if y k (t) takes values outside of D k , then those values would lack physical meaning. The following lemma shows that once y(t) enters D, it never leaves D.
In this section, we recall some preliminary results, pertinent to the analysis of system (5) . A real square matrix is said to be Metzler if all elements outside the diagonal are nonnegative. We require the following results for Metzler matrices. 
We will also be making use of the following variants of the Perron-Frobenius theorem for irreducible matrices. 
where f : G → R n is a locally Lipschitz map. Proposition 1. Let 0 be an equilibrium of (6) and E ⊆ G be a positively invariant and connected set with respect to (6), containing 0. Let V : E → R be a continuously differentiable and positive definite function, such thatV (x(t)) is negative definite. Further, let it hold that
is a bounded set for any constant c > 0, and is equal to E as c → ∞. Then the equilibrium 0 is asymptotically stable, with domain of attraction containing E.
This proposition can be proven using [39, Theorem 4.1]. The following lemma pertains to system (5), providing a constraint on any endemic equilibrium.
D is an equilibrium of (5), then, for each k ∈ [m], either y k = 0, or 0 y k 1. Moreover, m k=1 y k 1.
Proof: See Appendix. Lemma 6 states that when the underlying network is strongly connected, any endemic equilibrium involves each active virus infecting a separate fraction of each population node, and contaminating the shared resource to some degree.
In this section, we present sufficient conditions for the exponential (resp. asymptotic) stability of the eradicated state of a virus. The key condition is found through eigenvalue analysis of (B k w − D k w ), as seen in the following theorem.
Then the eradicated state of virus k is exponentially stable, with domain of attraction containing D k .
Proof: See Appendix. Theorem 1 states that if the linearized state matrix of virus k is Hurwitz, then, for all initial conditions in the sensible domain, virus k is eradicated exponentially fast. Theorem 1 answers the first part of question (i) in Section II-B. With respect to [23, Proposition 2], Theorem 1 is an improvement in the sense that it holds globally (on the sensible domain), and accounts for the multi-virus case, whereas [23, Proposition 2] established local exponential stability for the single-virus case.
Theorem 1 indeed guarantees exponential eradication of virus k, however, the condition is quite strict. For certain viruses it suffices to know whether or not the virus will be eradicated, but the speed with which this eradication takes place is of less importance. Indeed, it turns out that a relaxation of the strict inequality of the eigenvalue condition in Theorem 1 guarantees asymptotic eradication of a virus, as stated in the following theorem. Remark 3 (Epidemiological Interpretation). Observe that, due to Lemma 3, the conditions in Theorem 1 (Theorem 2) are equivalent to ρ(
. The fact that Theorem 1 (resp. Theorem 2) guarantees exponential (asymptotic) eradication of a virus k whenever
as the basic reproduction number of the virus in the network, typically denoted by R 0 . In epidemiology, this is the average number of secondary infections caused by an infected individual, before recovering. Thus, Theorem 1 (Theorem 2) states that whenever the basic reproduction number of a virus is strictly less than (less than or equal to) one, the virus will exponentially (asymptotically) converge to its eradicated state.
In this section, we study the possibility of viruses persisting in the network, corresponding to non-zero equilibria of (5). Naturally, the persistence of a virus must follow from the violation of the conditions of Theorem 2. The resulting behavior is detailed in the rest of this section. Before we proceed, we state the following result: Proof: See Appendix. With Lemma 7 in place, we have the following theorem for a single-virus system, guaranteeing existence of a unique, asymptotically stable, single-virus endemic equilibrium when the eigenvalue condition in Theorem 2 is violated. Proof: See Appendix. For a single-virus system, Theorem 3 states that when the eigenvalue condition in Theorem 2 is violated, then as long as some viral infection is present initially, the viral spreading process will converge to a unique infection ratio in each population node and a unique contamination level in the shared resource. This answers question (ii) in Section II-B.
Remark 4 (Epidemiological Interpretation). Applying Lemma 3, we see that the conditions of Theorem 3 are equivalent to ρ((D w ) −1 B w ) > 1. This is again consistent with the interpretation of ρ((D w ) −1 B w ) as the basic reproduction number of the virus, since a persisting virus should have a basic reproduction number greater than one.
Theorem 3 establishes the existence, uniqueness and asymptotic stability of a single-virus endemic equilibrium, extending [40, Theorem 2.4.] to the setting with a shared resource, whereas [15] illustrated this extension in simulations, without providing theoretical guarantees. We can partially extend Theorem 3 to the multi-virus case, specifically the existence and uniqueness of a single-virus endemic equilibrium for a virus, resulting in the following proposition.
Proof: See Appendix. Proposition 2 states that each virus violating the eigenvalue condition in Theorem 2 has a unique single-virus endemic equilibrium. This result is unsurprising, since nullifying the other viruses in the system reduces it to a single-virus system. However, note that Proposition 2 does not say anything about the stability of these equilibria. Nor does the proposition say anything about the case of more than one virus persisting in the population simultaneously. Due to Theorem 2 and Proposition 2, we have the following necessary and sufficient condition for the healthy state to be the unique equilibrium. 
. Theorem 4 states that as long as the largest real part of the eigenvalue of the linearized state matrix of each virus is non-positive, the healthy state is the only equilibrium of (5). Theorem 4 answers question (iii) in Section II-B. Note that Theorem 4 extends [23, Theorem 1] to the setting with more than two viruses and a shared resource, albeit under the assumption that the healing rate of each agent with respect to each virus is strictly positive.
Beyond the single-virus endemic equilibria from Section V, we would like to know when multiple viruses can persist in a population simultaneously, corresponding to a coexisting equilibrium. The following theorem makes use of a particular eigenvalue condition to show the existence of a coexisting equilibrium in a bi-virus system. Before proceeding to the statement of the theorem, recall from Proposition 2 that, if 
Moreover, we know that 0 ỹ 1 1 and 0 ỹ 2 1.
Theorem 5. Consider the SIWS model (5) under Assumption 1 with m = 2. Suppose that B 1 w and B 2 w are irreducible matrices, and that s( (8), then there exists at least one coexisting equilibrium (ŷ 1 ,ŷ 2 ) 0 in D such thatŷ 1 +ŷ 2 ≤ 1.
Proof: See Appendix. With each virus satisfying the condition for the existence of its single-virus endemic equilibrium, Theorem 5 states that if, for each virus, the largest real part of any eigenvalue of the matrix of the dynamics linearized around the single-virus endemic equilibrium of the other virus is positive, then both the viruses can simultaneously infect separate fractions of each population node. Theorem 5 answers question (iv) in Section II-B.
Remark 5 (Epidemiological Interpretation). Applying Lemma 3 to the conditions in Theorem 5, we see that they are equivalent to having ρ((I − X(ỹ 2 ))(D 1
as the invasion reproduction numbers of virus 1 invading virus 2 and virus 2 invading virus 1, respectively. The invasion reproduction number is defined for an invading pathogen, introduced into a setting with another, endemic pathogen at equilibrium. It is defined as the average number of secondary infections caused by an individual infected by the invading pathogen, at the time of introduction [41] . In line with this interpretation, Theorem 5 shows that coexistence is possible whenever both invasion reproduction numbers are greater than one.
While conditions that guarantee existence of coexisting equilibria may be found in [22] and [24] , these references do not account for the presence of a shared resource. In order to compare the results in [22] , [24] with Theorem 5, we particularize our model to the setting without a shared resource. Specifically, (4) reduces tȯ
Assumption 1, particularized for the setting without a shared resource, is given as follows:
Applying Proposition 2 to the bi-virus setting without a shared resource, we see that if i) B 1 and B 2 are irreducible, ii) s(B 1 − D 1 ) > 0, and iii) s(B 2 − D 2 ) > 0, then there exist exactly two single-virus endemic equilibria, namely
such that 0 p 1 1 and 0 p 2 1. Hence, we have the following corollary to Theorem 5, for the setting without a shared resource. , particularized for m = 2, establishes the existence of not just one but infinitely many coexisting equilibria, thus implying that a coexisting equilibrium in a bi-virus setup is not necessarily unique. However, it turns out that the conditions in Corollary 1 do not coincide with the conditions in [24, Theorem 6] , as discussed in the following remark.
Remark 6. Suppose that, for m = 2, the conditions in [24, Theorem 6] are satisfied. That is, for some v > 0, B 1 is an irreducible nonnegative matrix with B 1 = vB 2 , and D 1 is a positive diagonal matrix with D 1 = vD 2 , s(B 1 − D 1 ) > 0 and s(B 2 −D 2 ) > 0. Letp 1 andp 2 be the unique single-virus endemic equilibria defined in (11) . It follows from (9) , and the fact that D 1 , D 2 are invertible, thatp 1 ,p 2 fulfill
Since
it follows from (12) thatp 1 =p 2 =p, and therefore
By Proposition 2 we havep 1, implying that (I − diag(p)) is a positive diagonal matrix. Then, for k ∈ [2] , it follows that (I − diag(p))(D k ) −1 B k is an irreducible nonnegative matrix. Therefore, item (iii) in Lemma 4 can be applied to (13) , from which it follows that ρ((I − diag(p))(D 1 ) −1 B 1 ) = 1 and ρ((I −diag(p))(D 2 ) −1 B 2 ) = 1. Applying Lemma 3 we obtain
Observe (14) is incompatible with the following condition in Corollary 1, s(−D 1 + (I − diag(p 2 ))B 1 ) > 0,
Hence, it follows that the conditions in Corollary 1 and [24, Theorem 6] are mutually exclusive.
Since the conditions in Corollary 1 are not compatible with the conditions in [24, Theorem 6] , the question: "do the conditions in Corollary 1 guarantee uniqueness of the coexisting equilibrium?" is worth investigating. Our simulations show that such a coexisting equilibrium may indeed be unique, and further, asymptotically stable; see Section VIII and Figure 4 .
While Theorem 5 provides conditions for the existence of coexisting equilibria, a related problem is finding conditions under which no coexisting equilibria can exist. The following theorem makes use of a nontrivial condition to eliminate the possibility of coexisting equilibria in a bi-virus setting, and establishes one virus as being dominant. 
w then there are exactly three equilibria in D, namely the healthy state, which is unstable, (0,ỹ 2 ) with 0 ỹ 2 1, which is unstable, and (ỹ 1 , 0) with 0 ỹ 1 1, which is locally exponentially stable.
Proof: See Appendix.
Theorem 6 states that if one virus has a stronger set of spread and healing parameters than the other virus, then these two viruses cannot coexist in the population. Theorem 6 answers question (v) in Section II-B. The following remark aids in understanding the result in Theorem 6.
Remark 7. Underlying Theorem 6 is the so-called competitive exclusion principle, which states that complete competitors cannot coexist [42] . For instance, if two strains of viruses (say virus 1 and virus 2) compete with each other to infect the same population, and if virus 1 has a slight advantage over virus 2 (for example higher infection rates or lower healing rates), then virus 1 will eventually displace virus 2, which will get eradicated.
Similar results to Theorem 6 can be found in [23, Theorem 5], for the setting without a shared resource. In order to make a comparison, we employ the SIS model (10) to state the following corollary of Theorem 6. 1, which is unstable, and (p 1 , 0) with 0 p 1 1, which is locally exponentially stable.
Corollary 2 is directly comparable to [23, Theorem 5] . The main difference is that [23, Theorem 5] requires that β k ij = β k or β k ij = 0, and δ k i = δ k , for all i, j ∈ [n], k ∈ [2] , and some β k > 0, δ k > 0, whereas Corollary 2 has no such restrictions. Hence, Corollary 2 subsumes and improves upon [23, Theorem 5] , while Theorem 6 extends the result to the setting with a shared resource.
In this section, we present two strategies for ensuring that some (or all) viruses are eradicated. First, we establish that it is possible to cause convergence to the eradicated state of a virus k by boosting the corresponding healing rates δ k i . Applying this technique to all m viruses, the system converges to the healthy state. Second, we show that it is possible to leverage a benign virus in order to eradicate a malignant virus, in a bi-virus setting.
When managing the epidemic spread of a virus, a natural strategy is to boost the healing rates in the population nodes. The following result shows that the healing rates can always be chosen to ensure asymptotic or exponential eradication of a virus, using an algorithm inspired by [23] , [34] . 
, the eradicated state of virus k is exponentially stable, with domain of attraction containing D k . Otherwise, the eradicated state of virus k is asymptotically stable, with domain of attraction containing D k .
The proof is straightforward, following from Theorem 2 and similar arguments as in [23, Section V] , and [34] . Proof: See Appendix.
Proposition 3 represents one strategy to ensure eradication of a virus. By applying this strategy to all viruses, we obtain the following result. Theorem 7 represents a strategy to eradicate all viruses in a system, which equivalently ensures exponential (resp. asymptotic) convergence to the healthy state. Thus, Theorem 7 addresses question (vi) in Section II-B.
The mitigation strategy outlined in Theorem 7 could be understood as follows: with respect to each virus, if the healing rate of each subpopulation is sufficiently increased, which could be accomplished by prescribing high dosages of drugs, by administering vaccines, etc., then each of the viruses get eradicated exponentially (resp. asymptotically) fast. Observe that this strategy is extreme in the sense that it does not factor in limitations on the availability of resources, and essentially encourages health administration officials to amass (possibly) excessive amounts of resources in order to prevent epidemic outbreaks.
Note that in the absence of sufficient resources, implementing the aforementioned strategy in practice is impossible. Hence, we are motivated to seek different strategies.
It turns out that in a bi-virus setting, where one virus is malignant and the other virus is benign, we can leverage the benign virus in order to help eradicate the malignant virus, as stated in the following theorem. 
is the set of directed edges between the population nodes and the shared resource, with respect to virus 1 (resp. virus 2). If the healing rates for virus 2 fulfill
for all i ∈ [n], then the only locally asymptotically stable equilibrium in D is (ỹ 1 , 0) with 0 ỹ 1 1.
Proof: See Appendix. Theorem 8 represents a strategy to eradicate one of the viruses in a bi-virus system, made possible by leveraging the fact that one virus has a stronger set of spread parameters than the other. Thus, Theorem 8 addresses question (vii) in Section II-B.
Remark 9 (Virus as vaccine). Since the strategy given in Theorem 8 ensures local asymptotic convergence to the singlevirus endemic equilibrium of the benign virus, it could also be interpreted in the following sense: the benign virus effectively acts as a vaccine against the malignant virus. In the context of battling epidemic outbreaks, where the goal is to minimize the mortality rate, this strategy could potentially provide health administration officials with an effective tool.
The mitigation strategies detailed in this section can be compared as follows. On the one hand, assuming that the objective of public health officials is solely to eradicate one virus in a bi-virus system while considering resource constraints (e.g., availability of vaccines, drugs, ventilators, etc.), it may be more feasible to implement the strategy given in Theorem 8 instead of that in Proposition 3. On the other hand, the strategy in Theorem 8, opposed to that in Theorem 7 particularized for a bi-virus setting, requires the persistence of one virus, which may be undesirable. Moreover, with respect to a given virus, if the healing rate of at least one node is sufficiently boosted, then Proposition 3 guarantees eradication of said virus exponentially fast, whereas Theorem 8 guarantees only asymptotic eradication of a virus. We explore the advantages and disadvantages of the strategies outlined in Proposition 3 and Theorem 8 via simulations in Section VIII.
In this section, we present a number of simulations to illustrate our theoretical findings, using the city of Stockholm as an example setting. In particular, 15 districts in and around Stockholm are taken to be the population nodes of the network, thus, n = 15. All of these districts are connected to the Stockholm metro, which we view as the shared resource of the network; see Figure 2 . In all simulated scenarios we consider two competing viruses, namely virus 1 and virus 2. We denote the average infection ratio of virus k, i.e., 1 n n i p k i (t), byp k (t). The contact network of each virus is the same, i.e., E 1 = E 2 , and is represented in Figure 2 . We set p k i (0) = 0.5 and z k (0) = 0.5 for all i ∈ [15] and k ∈ [2] , and use the same in all simulated scenarios throughout this section. Observe that the aforementioned choice of initial states represents the case where half of the population in each node is infected by virus 1, while the other half is infected by virus 2, and the shared resource is contaminated by both viruses. The spread parameters β k ij are taken to be 1 if district i is adjacent to district j, or if i = j, and 0 otherwise, for k ∈ [2] . In all scenarios, we also assume that each district is bi-directionally connected to the shared resource, i.e., the Stockholm metro, with β k iw = 1 and c k i = 1/15, for all i ∈ [15] and k ∈ [2] . As a consequence, B k w is irreducible for k ∈ [2] . The following scenarios differ only in terms of the choice of δ k i and δ k w . In the simulation depicted in Figure 3 , we chose δ 1 i = 4.6, δ 1 w = 4 and δ 2 i = 10, δ 2 w = 10, for all i ∈ [15] . As a consequence, s(B 2 w − D 2 w ) = −4.2, so Theorem 1 applies to virus 2. Consistent with the result in Theorem 1, virus 2 becomes eradicated exponentially fast. However, since s(B 1 w − D 1 w ) = 0.3, virus 1 fulfills the conditions for Proposition 2, and furthermore, once virus 1 is eradicated, the system is essentially a single-virus system. Treated as such, virus 1 satisfies the conditions in Theorem 3. In line with the result in Theorem 3, virus 1 converges to a single-virus endemic equilibrium. Insofar as it can be determined by varying y 1 (0) in D 1 , this single-virus endemic equilibrium appears to be unique and asymptotically stable. In the simulation depicted in Figure 4 , we chose δ 1 i = 1.5 for i ∈ [8] , and δ 1 i = 2 for i ∈ [15] \ [8] , with δ 1 w = 1. Relating these choices of parameters for virus 1 to Figure 2 , this corresponds to a low healing rate north of the river Mälaren and a high healing rate south of the river. Mirroring this pattern, we chose δ 2 i = 2 for i ∈ [8] , and δ 2 i = 1.5 for i ∈ [15] \ [8] , with δ 2 w = 1. With these parameters, it follows that s(B 1 w − D 1 w ) = 2.8, and s(B 2 w − D 2 w ) = 2.8. Hence, both viruses fulfill the conditions for Proposition 2, providing the existence of exactly two single-virus endemic equilibria, namely (ỹ 1 , 0) and (0,ỹ 2 ).ỹ 1 can be approximated by setting y 1 (0) > 0 and y 2 (0) = 0, and running the simulation for a sufficiently long period of time T . Then, assuming that y 1 ≈ y 1 (T ), and, with an analogous approximation for virus 2, y 2 ≈ y 2 (T ), we obtain s((I − X(ỹ 2 ))B 1 w − D 1 w ) = 0.2, and s((I − X(ỹ 1 ))B 2 w − D 2 w ) = 0.2. As a consequence, this pair of viruses fulfills the conditions for Theorem 5.
In line with the result in Theorem 5, we see that there exists a coexisting equilibrium. Moreover, our simulations show that the viral infection levels appear to converge to this coexisting equilibrium. Additionally, regardless of how the initial condition is varied within D, barring y 1 (0) = 0 or y 2 (0) = 0, we observe that all simulations converge to the same coexisting equilibrium. This observation suggests that the coexisting equilibrium might be unique, as well as asymptotically stable. In the simulation depicted in Figure 5 , we chose δ 1 i = 3, δ 1 w = 3 and δ 2 i = 4, δ 2 w = 4, for all i ∈ [15] . It follows that s(
Therefore, this pair of viruses fulfills the conditions for Theorem 6. In line with the result in Theorem 6, we observe that virus 1 has a competitive edge, allowing it to persist and converge to its single-virus endemic equilibrium. Meanwhile, virus 2 is eradicated, despite s(B 2 w − D 2 w ) being positive. The simulations depicted in Figures 6 and 7 were both initialized at the coexisting equilibrium from Figure 4 , with all parameters as in that simulation. At t = 0.5 in the simulation in Figure 6 , the healing rates of virus 2, i.e., δ 2 i , are chosen as in (16), with 2 i = 0 for all i ∈ [15] . Similarly, at t = 0.5 in the simulation in Figure 7 , the healing rates of virus 2, i.e., δ 2 i , are chosen as in (17) . More specifically, δ 2 i = 2 for i ∈ [8] as before, but δ 2 i = 2.1 for i ∈ [15] \ [8] . We assume that the cost of a strategy (denoted as cost ) is given by the sum of all healing rates (i.e., cost =
i is chosen as in (16) when = 1, and as in (17) when = 2). Our simulations show that cost 1 = 61, whereas cost 2 = 33. However, as seen in Figures 6 and 7 , the end result of both strategies is the same. That is, virus 1 persists and reaches its single-virus endemic equilibrium, whereas virus 2 is eradicated. Hence, the strategy in Theorem 8 allows us to eradicate virus 2 at a lower cost than the strategy in Proposition 3, in this case. However, the convergence to the single-virus endemic equilibrium of virus 1 is faster with Proposition 3, as expected from the excessive healing rates.
In this paper, we introduced and analyzed a novel SIWS model of competitive, multi-viral spread across a network of population nodes with a shared resource. We established conditions under which a virus is eradicated exponentially fast (resp. asymptotically). We also provided a condition that ensures a virus can reach and sustain a unique, single-virus endemic equilibrium of infection levels across the network. Taken together, these results allowed us to state a necessary and sufficient condition for convergence to the healthy state. Moreover, we presented sufficient conditions under which two viruses can sustain a coexisting equilibrium, that is, where neither virus is eradicated. Conversely, we also provided a necessary condition for the existence of such a coexisting equilibrium. To mitigate the spreading process, we proposed two strategies. The first strategy involved choosing the healing rate of each subpopulation, with respect to each virus, in a suitable manner so as to ensure that all subpopulations converge to the healthy state. The second strategy exploited the notion of competitive exclusion to ensure that in a bi-virus setup, with one virus being malignant and the other benign, the malignant virus becomes eradicated.
The present paper concerns time-invariant SIWS models, whereas in order to obtain a better understanding of real-world scenarios (viz. mobile agents, mutating viruses), it is natural to consider time-varying SIWS models. Secondly, in order to generalize the result on the existence of a coexisting equilibrium (cf. Theorem 5) to more than two viruses, we would need to establish uniqueness of the coexisting equilibrium for the bi-virus case, which is non-trivial.
Consider a virus k such that s(B k w −D k w ) < 0. By Lemma 1 we know that D is positively invariant with respect to (5), so because y(0) ∈ D we have y(t) ∈ D for all t ≥ 0. Then, (4) with respect to virus k can be bounded by:
Since y k (t) ≥ 0 for all t ≥ 0, it follows from Grönwall-Bellman's Inequality [39, pg 651 ] that the solution of (4) will be bounded above by the solution of the linear system corresponding to (18) with equality. Since s(B k w − D k w ) < 0, we know that 0 is globally exponentially stable in the linear system. Therefore, the eradicated state of virus k is exponentially stable, with domain of attraction containing D k .
Consider an equilibrium y ∈ D of system (5) . Assume, by way of contradiction, that m k=1 p k i ≥ 1 for some i ∈ [n]. Plugged into (2) under Assumption 1, we obtain m k=1ṗ
where (19) follows from i) Assumption 1, ii) m k=1 p k i ≥ 1 and iii) that y ∈ D. Note that (19) is a contradiction of the fact that y is an equilibrium, following from the assumption m k=1 p k i ≥ 1 for some i ∈ [n]. Therefore, m k=1 p k 1. Now, assume by way of contradiction that m k=1 z k ≥ 1. Plugged into (2) under Assumption 1, we obtain m k=1ż 
The inequality (21) follows from Assumption 1, and the assumption that m k=1 z k ≥ 1. Note that (21) is a contradiction of the fact that y is an equilibrium of system (5) . Since this contradiction follows from the assumption m k=1 z k ≥ 1, we must have m k=1 y k 1, and therefore y k 1, for all k ∈ [m], for any equilibrium y ∈ D. Now, for all k ∈ [m], y k is a equilibrium of (4), so we have
, assume by way of contradiction that y k > 0, with y k i = 0 for all i ∈ W , where W ⊂ [n + 1] is nonempty. Then, by the properties of irreducible nonnegative matrices, ((I − m l=1 X(y l ))(D k w ) −1 B k w y k ) j > 0 for some j ∈ W . Since y k j = 0, this contradicts (22) , and therefore we must either have y k 0, or y k = 0, for each k ∈ [m].
Consider a virus k such that s(B k w −D k w ) ≤ 0. Since y(0) ∈ D, Lemma 1 states that we have y(t) ∈ D for all t ≥ 0, and further that D k is positively invariant with respect to (4) . Now, note that if s(B k w − D k w ) < 0, Theorem 1 implies that the eradicated state is exponentially stable for virus k, with domain of attraction containing D k . Hence, the rest of the proof considers the case when
is an irreducible Metzler matrix, by Lemma 2 there exists a positive diagonal matrix Q k such that
Define the Lyapunov function candidate V (y k (t)) = y k (t) T Q k y k (t) with D k as the domain. Note that V (y k (t)) 0, and that V (y k (t)) fulfills (7) .
We will now show thatV (y k ) < 0 if y k > 0. First, consider the case where y k i = 0 for some i ∈ [n + 1]. Then, noting that y kT Q k X(y k )B k w y k ≥ 0, we see that (23) is bounded byV (y k ) ≤ y kT Λ k y k . (24) Given that Λ k 0, we have s(Λ k ) ≤ 0. Then, since Λ k is an irreducible Metzler matrix, it follows from Lemma 5 that r := s(Λ k ) is a simple eigenvalue of Λ k , with a corresponding eigenvector ζ 0. Since we consider the case where y k i = 0 for some i ∈ [n + 1], y k can not be parallel to ζ. By the Rayleigh-Ritz Theorem [43, Theorem 4.2.2] , y kT Λ k y k = ry kT y k only if y k is parallel to ζ, and y kT Λ k y k < ry kT y k otherwise. Since r ≤ 0, it follows from (24) thatV (y k ) < 0 when y k i = 0 for some i ∈ [n + 1]. Now, consider the case when y k 0. Recall that Λ k 0, hence, (23) is bounded bẏ
w is an irreducible nonnegative matrix, Q k is a positive diagonal matrix and y k (t) 0, it follows that y kT Q k X(y k )B k w y k > 0. Hence, (25) gives usV (y k ) < 0 when y k 0. Then we haveV (y k ) < 0 for all y k > 0, and it is clear thatV (0) = 0. Therefore,V (y k ) ≺ 0. Finally, since D k is positively invariant with respect to (4), we see that V (y k (t)) meets the conditions for Proposition 1. This shows that the eradicated state of virus k is asymptotically stable, with domain of attraction containing D k .
Lemma 1 states that D is positively invariant with respect to system (5) . It remains to be shown that if y(0) ∈ D, and y(0) > 0, it follows that y(t) > 0 for all t > 0. With y(0) > 0 we have y k (0) > 0 for some k ∈ [m], and by Lemma 1 we know that y k (t) ≥ 0 for all t > 0. Then, (4) is bounded bẏ
due to i) y k (t) ≥ 0, and ii) that (I − m l=1 X(y l (t)))B k w is a nonnegative matrix. Integrating (26) yields
for all t > 0. The final inequality in (27) follows from the fact that D k w is a positive diagonal matrix. Since y k (t) > 0, it follows that y(t) > 0 for all t > 0. Therefore, D \ {0} is positively invariant with respect to system (5).
Part 1 -Proof of existence: Note that if y > 0, diag(D −1 w B w y) is a nonnegative diagonal matrix, and therefore the inverse of (I + diag(D −1 w B w y)) exists. Define a map T (y) : R n+1
Observe that the components of T (y) are
Note that the scalar function s/(1 + s) is increasing in s, and that D −1 w B w is a nonnegative matrix. Therefore, v ≥ z implies T (v) ≥ T (z). Now, observe that a fixed point of T (y) fulfills
For a given y ∈ R n+1 + , define X(y) to be its diagonalization with the final element y n+1 set to zero. As such, subtracting
Since X(y) and D −1 w are diagonal matrices, they commute. Furthermore, by pre-multiplying (31) with D w , and suitably rearranging terms, we obtain (32) is clearly an equilibrium of (4) with m = 1. As such, it suffices to show that T (y) has a nonzero fixed pointỹ in D. We will now show that at least one such fixed point exists. Since s(B w − D w ) > 0, by Lemma 3,
Furthermore, by item (ii) in Lemma 4, we know that the eigenspace of λ * is spanned by a vector y * 0. Then, since λ * > 1, there exists some constant > 0 such that, for all i ∈ [n + 1], we have
Noting that (D −1 w B w y * ) n+1 y * n+1 > 0, we also have
Due to the inequalities (33) and (34), we have T ( y * ) ≥ y * .
Since y ≥ z implies T (y) ≥ T (z), it follows that for any y ≥ y * we have T (y) ≥ y * . Consider T (1) for i ∈ [n],
For i = n + 1, we have
Due to (35) and (36), we have T (1) ≤ 1.
Since v ≥ w implies T (v) ≥ T (w), it follows that T (y) ≤ 1 if y ≤ 1. By Brouwer's fixed-point theorem [44, Theorem 9.3] , there is at least one fixed point of T (y) in the domain {y : y * ≤ y ≤ 1}. Since a fixed point of T (y) is equivalent to an equilibrium of (5), by Lemma 6, any fixed point must fulfill y 1. In conclusion, the map T (y) has at least one fixed point in the domain {y : y * ≤ y 1}, and therefore system (5) has at least one equilibriumỹ ∈ D such that 0 ỹ 1.
We will now prove that the single-virus endemic equilibrium is unique. Suppose that there are two single-virus endemic equilibria,ỹ andỹ. By Lemma 6 we haveỹ 0 andỹ 0. Let κ = max i∈[n+1]ỹi /ỹ i . First we show that κ is given by
To do this, assume by way of contradiction that κ = y n+1 /ỹ n+1 , and κ >ỹ i /ỹ i , for all i ∈ [n]. Note that since bothỹ andỹ are equilibria of system (4), we haveỹ
Since we assume that κ >ỹ i /ỹ i , we have κỹ i >ỹ i , for all i ∈ [n]. Then, (38) gives us y n+1 = n i c iỹi < κ n i c iỹi = κỹ n+1 , Hence, κ >ỹ n+1 /ỹ n+1 , which contradicts the assumption that κ =ỹ n+1 /ỹ n+1 . Therefore, κ must be given by equation (37) . Now, by (37) we know thatỹ ≤ κỹ. For some j ∈ [n] we haveỹ j = κỹ j . Assume, by way of contradiction, that κ > 1. Then, using the fact that an equilibrium of (4) also constitutes a fixed point of T (y) (see part 1 of this proof), we havẽ
where (39) follows fromỹ ≤ κỹ and that T (v) ≥ T (w) whenever v ≥ w, (40) follows from the assumption κ > 1, and (41) follows from the fact thatỹ is an equilibrium of (4). Note that (42) is a contradiction, following from our assumption that κ > 1. Hence, κ ≤ 1, meaning thatỹ ≤ỹ.
Switching the roles ofỹ andỹ, we see thatỹ ≤ỹ. Therefore, y =ỹ, and thus the equilibrium is unique.
Recall thatỹ is the unique single-virus endemic equilibrium of system (4) 
where (44) follows from (43) , and (45) follows from X(u)v = X(v)u. Now, since D w is invertible, (43) is equivalent to
Since B w is an irreducible nonnegative matrix, and D w is a positive diagonal matrix,ỹ 1 ensures that (I − X(ỹ))D −1 w B w is an irreducible nonnegative matrix, and in turn that (−D w + (I − X(ỹ))B w ) is an irreducible Metzler matrix. Therefore, sinceỹ 0, item (iii) in Lemma 4 applied to (46) gives us ρ((I − X(ỹ))D −1 w B w ) = 1, which by Lemma 3 is equivalent to s(−D w + (I − X(ỹ))B w ) = 0. Then, Lemma 2 guarantees the existence of a positive diagonal matrix Q such that Ψ := (−D w + (I − X(ỹ))B w ) T Q + Q(−D w + (I − X(ỹ))B w ) 0. Define the Lyapunov function candidate V (∆y(t)) = ∆y(t) T Q∆y(t), with y(t) ∈ D \ {0} as the domain. Note that V (∆y(t)) 0 and that V (∆y(t)) fulfills (7) . Differentiating V (∆y(t)) with respect to t yieldṡ V (∆y(t)) = 2∆y T Q(−D w + (I − X(ỹ))B w )∆y − 2∆y T QX(B w y)∆y (47)
= ∆y T Ψ∆y − 2∆y T QX(B w y)∆y (48) where (47) makes use of (45). We want to show thaṫ V (∆y(t)) < 0 for all y(t) ∈ D \ {0} such that ∆y(t) = y(t) −ỹ = 0. First, consider all y(t) ∈ D \ {0} such that y(t) 0. Combining (48) and the fact that Ψ 0 gives uṡ V (∆y(t)) ≤ −2∆y(t) T QX(B w y(t))∆y(t).
(49) Note that, given y(t) 0, QX(B w y(t)) is a positive diagonal matrix, and thus QX(B w y(t)) 0. Therefore, (49) gives uṡ V (∆y(t)) < 0 for all y(t) ∈ D \ {0} such that y(t) 0, ∆y(t) = y(t) −ỹ = 0. Now, consider all y(t) ∈ D \ {0} such that y(t) > 0, y i (t) = 0 for some i ∈ [n + 1]. Note that, given y(t) > 0, QX(B w y(t)) is a nonnegative diagonal matrix. Then (48) can be bounded bẏ V (∆y(t)) ≤ ∆y(t) T Ψ∆y(t).
(50) Given that (−D w + (I − X(ỹ))B w ) is an irreducible Metzler matrix and Q is a positive diagonal matrix, Ψ is an irreducible Metzler matrix. Employing (43) , we see that y T Ψỹ = 0.
(51) Item (i) in Lemma 5 stipulates that r := s(Ψ) is a simple eigenvalue of Ψ. Due to Ψ 0 and the Rayleigh-Ritz Theorem [43, Theorem 4.2.2] , it follows from (51) that r = 0, and that y spans the eigenspace of r. As such, due toỹ 0, and y(t) > 0 with y i (t) = 0 for some i ∈ [n + 1], y(t) can not be parallel toỹ. As a consequence, ∆y(t) can not be parallel toỹ. By the Rayleigh-Ritz Theorem [43, Theorem 4.2.2 ], x T Ψx = rx T x only if x is parallel toỹ, and x T Ψx < rx T x otherwise. Therefore, r = 0 together with (50) gives usV (∆y(t)) < 0 for all y(t) ∈ D \ {0} such that y(t) 0 and y(t) −ỹ = 0. Thus, we haveV (∆y(t)) < 0 for all y(t) ∈ D \ {0} such that ∆y(t) = y(t) −ỹ = 0, and it is clear thatV (0) = 0. Therefore,V (∆y(t)) ≺ 0 for y(t) ∈ D \ {0}.
Finally, from Lemma 7 we have that D \ {0} is a positively invariant set with respect to (5) . Thus, we see that V (∆y(t)) meets the conditions for Proposition 1 with respect to the shifted coordinates ∆y(t) = y(t) −ỹ, for all y(t) ∈ D \ {0}. This shows that the unique single-virus endemic equilibriumỹ is asymptotically stable, with domain of attraction containing D \ {0}. Thus, with parts 1, 2 and 3 in place, the proof of Theorem 3 is concluded.
Suppose that, for some k ∈ [m], B k w is irreducible and s(B k w − D k w ) > 0, and that y l = 0 for all l ∈ [m], l = k. Then the dynamics of virus k can be written aṡ
Note that (52) corresponds to the dynamics of the single-virus case. Therefore, since B k w is irreducible and s(B k w −D k w ) > 0, it follows from the first and second parts of the proof of Theorem 3 that there exists a unique single-virus endemic equilibrium of the form (0, . . . ,ỹ k , . . . , 0), with 0 ỹ k 1 in D. This holds for each k ∈ [m] such that B k w is irreducible and s(B k w − D k w ) > 0, by repeating the arguments above.
Recall that for k ∈ [2] , Assumption 1 implies that D k w is a positive diagonal matrix, and therefore invertible. Furthermore, note that (I + diag((D 1 w ) −1 B 1 w y 1 )) and (I + diag((D 2 w ) −1 B 2 w y 2 )) are positive diagonal matrices whenever y 1 ≥ 0 and y 2 ≥ 0, and are then also invertible. For y ∈ R n+1 + , define X(y) to be diag(y) with y n+1 set to zero. Define the maps T 
For i ∈ [n], the i th components of the maps are
Furthermore, the (n + 1) th components of the maps are
Note that the scalar function s/(1 + s) is increasing in s, and for k ∈ [2] , the matrix (D k w ) −1 B k w is nonnegative, therefore, T k i is an increasing function in y k j for all i, j ∈ [n + 1]. Moreover, T 1 i is a decreasing function in y 2 i and T 2 i is a decreasing function in y 1 i , for all i ∈ [n]. Hence, for any
The inequalities in (53) state that T k (y 1 , y 2 ) is increasing in its k th argument and decreasing in its other argument, for k ∈ [2] . Let y = (y 1 , y 2 ), and let T (y) : [0, 1] 2(n+1) → [0, 1] 2(n+1) be the map T (y) = (T 1 (y), T 2 (y)). A fixed point of T (y) fulfills
Pre-multiplying the first line (resp. second line) of (54) by (I + diag((D 1 w ) −1 B 1 w y 1 )) (resp. (I + diag((D 2 w ) −1 B 2 w y 2 ))) gives us
from the first (resp. second line) yields (I − X(y 1 ) − X(y 2 ))(D 1 w ) −1 B 1 w y 1 = y 1 ,
Making use of the fact that diagonal matrices commute, premultiplying the first line (resp. second line) of (56) by D 1 w (resp. D 2 w ), and rearranging terms gives us (−D 1 w + (I − X(y 1 ) − X(y 2 ))B 1 w )y 1 = 0, (−D 2 w + (I − X(y 1 ) − X(y 2 ))B 2 w )y 2 = 0.
Comparing (57) with (4), it follows that a fixed point of T (y) constitutes an equilibrium of system (5) and vice versa. It suffices to show that T (y) has a fixed pointŷ = (ŷ 1 ,ŷ 2 ) 0, such thatŷ 1 +ŷ 2 ≤ 1.
Recall that (ỹ 1 , 0) and (0,ỹ 2 ) are single-virus endemic equilibria of system (5) . Consider T 1 (ỹ 1 , y 2 ). By assumption, (ỹ 1 , 0) is an equilibrium of (5), therefore T 1 (ỹ 1 , 0) =ỹ 1 . By the inequalities in (53) we have T 1 (ỹ 1 , y 2 ) ≤ỹ 1 , and thus T 1 (y 1 , y 2 ) ≤ỹ 1 , for all y 1 ≤ỹ 1 . Analogously, it can be shown that we have T 2 (y 1 , y 2 ) ≤ỹ 2 , for all y 2 ≤ỹ 2 . Thus,
whenever (y 1 , y 2 ) ≤ (ỹ 1 ,ỹ 2 ). Now, by assumption, s(−D 1 w + (I − X(ỹ 2 ))B 1 w ) > 0, and, since D 1 w and (I − X(ỹ 2 )) are positive diagonal matrices and B 1 w is an irreducible nonnegative matrix, (−D 1 w + (I − X(ỹ 2 ))B 1 w ) is an irreducible Metzler matrix. Then, by Lemma 3 and the fact that diagonal matrices commute, we have ρ((I − X(ỹ 2 ))(D 1
is an irreducible nonnegative matrix, by item (i) in Lemma 4 we know that
is a simple eigenvalue of this matrix. Furthermore, by item (ii) in Lemma 4, we know that the eigenspace of λ 1 is spanned by a vectorȳ 1 0. Analogously we get λ 2 = ρ((I − X(ỹ 1 ))(D 2 w ) −1 B 2 w ) > 1, and the corresponding eigenvectorȳ 2 0. With the eigenvectorsȳ 1 ,ȳ 2 in place, we see that since
. Further, given thatȳ 1 0,ȳ 2 0,ỹ 1 0, andỹ 2 0, we haveỹ 1 i /ȳ 1 i > 0, andỹ 2 i /ȳ 2 i > 0, for all i ∈ [n + 1]. Moreover, note that λ 1 − 1 > 0 and λ 2 − 1 > 0. Hence, there exist 1 > 0 and 2 > 0 such that
From (59) follows that
Employing (60) it follows that, for all i ∈ [n], we have
Given that (59) implies 1ȳ1 <ỹ 1 and 2ȳ2 <ỹ 2 , by the inequalities in (53) we have T 1 ( 1ȳ1 , y 2 ) > 1ȳ1 whenever 2ȳ2 ≤ y 2 ≤ỹ 2 , and T 2 (y 1 , 2ȳ2 ) > 2ȳ2 whenever 1ȳ1 ≤ y 1 ≤ỹ 1 . Further application of the inequalities in (53) yields
whenever ( 1ȳ1 , 2ȳ2 ) ≤ (y 1 , y 2 ) ≤ (ỹ 1 ,ỹ 2 ). Then, (58) and (61) show that ( 1ȳ1 , 2ȳ2 ) ≤ T (y 1 , y 2 ) ≤ (ỹ 1 ,ỹ 2 ) whenever ( 1ȳ1 , 2ȳ2 ) ≤ (y 1 , y 2 ) ≤ (ỹ 1 ,ỹ 2 ). By Brouwer's fixed point theorem [44, Theorem 9.3] , there exists at least one fixed point of T (y) in the domain {y = (y 1 , y 2 ) :
Recall that a fixed point of T (y) is equivalent to an equilibrium of (5), hence, by Lemma 6, any fixed point of T (y) must fulfill y 1 + y 2 ≤ 1. In conclusion, system (5) has at least one coexisting equilibrium (ŷ 1 ,ŷ 2 ) 0 in D, such thatŷ 1 +ŷ 2 ≤ 1.
In order to prove Theorem 6, we require the following lemma.
Lemma 8. Consider system (5) under Assumption 1 with m = 2. Suppose that B 1 w and B 2 w are irreducible matrices, that s(B 1 w − D 1 w ) > 0 and s(B 2 w − D 2 w ) > 0, and that 0) and (0,ỹ 2 ) are singlevirus endemic equilibria of (5), thenỹ 1 >ỹ 2 .
Proof: Given thatỹ 1 andỹ 2 are equilibria of (4), and observing that diagonal matrices commute, we obtain
By Lemma 6 we have 0 ỹ 1 1 and 0 ỹ 2 1. Let κ = max i∈[n+1]ỹ 2 i /ỹ 1 i , and thusỹ 2 ≤ κỹ 1 . Since
, we have c 1 = c 2 . Then, by analogous arguments to those in part 2 of the proof of Theorem 3, we know that κ = max i∈[n]ỹ 2 i /ỹ 1 i . Let j be the index in [n] such thatỹ 2 j = κỹ 1 j . Assume, by way of contradiction, that κ > 1, implying y 2 j >ỹ 1 j . Note that since (D 1 w ) −1 B 1 w > (D 2 w ) −1 B 2 w , with both matrices being irreducible and nonnegative, it follows that (D 1 w ) −1 B 1 w y > (D 2 w ) −1 B 2 w y for any y 0. Further, note thatỹ 1 1, ensuring (1 −ỹ 1 j ) > 0. Then, (62) gives us
where (63) follows from (D 1 w ) −1 B 1 w > (D 2 w ) −1 B 2 w , (64) follows fromỹ 2 < κỹ 1 , (65) follows from (62), (66) follows fromỹ 2 j = κỹ 1 j , and (67) follows fromỹ 2 j >ỹ 1 j . Note that (67) is a contradiction, following from our assumption that κ > 1. Therefore, κ ≤ 1, and hence,ỹ 1 ≥ỹ 2 . Assume, by way of contradiction, thatỹ 1 =ỹ 2 =ỹ. Then (I − X(ỹ))(D 1 w ) −1 B 1 w > (I − X(ỹ))(D 2 w ) −1 B 2 w , and it follows from (62) that
It is clear that (68) is a contradiction, following from our assumption thatỹ 1 =ỹ 2 . Therefore we haveỹ 1 >ỹ 2 .
Proof of Theorem 6: Recall that the healthy state is an equilibrium of (5). Since s(B 1 w −D 1 w ) > 0 and s(B 2 w −D 2 w ) > 0, by Proposition 2 there are exactly two single-virus endemic equilibria of system (5), namely (ỹ 1 , 0) and (0,ỹ 2 ), such that 0 ỹ 1 1 and 0 ỹ 2 1. We will now show that with (D 1 w ) −1 B 1 w > (D 2 w ) −1 B 2 w , there are no equilibria in D other than the healthy state, (ỹ 1 , 0) and (0,ỹ 2 ). First, note that by Lemma 6, any additional equilibrium in D must be of the formŷ = (ŷ 1 ,ŷ 2 ), such thatŷ 1 0,ŷ 2 0, andŷ 1 +ŷ 2 1. Assume, by way of contradiction, that such an equilibriumŷ exists. Sinceŷ is an equilibrium of (5) it follows that (I − X(ŷ 1 ) − X(ŷ 2 ))(D 1 w ) −1 B 1 wŷ 1 =ŷ 1 ,
Given thatỹ 1 andỹ 2 are equilibria of (4), and observing that diagonal matrices commute, we obtain
Sinceỹ 1 0, we can define
Thus,ŷ 1 +ŷ 2 ≤ κỹ 1 . Let j ∈ [n + 1] be the index such that (ŷ 1 +ŷ 2 ) j = κỹ 1 j . Assume, by way of contradiction, that κ ≥ 1, implying (ŷ 1 +ŷ 2 ) j ≥ỹ 1 j . Note that since (D 1 w ) −1 B 1 w > (D 2 w ) −1 B 2 w , with both matrices being irreducible and nonnegative, it follows that (D 1 w ) −1 B 1 w y > (D 2 w ) −1 B 2 w y for any y 0. Further, note thatỹ 1 1, ensuring (1 −ỹ 1 j ) > 0. We will now show that κ ≥ 1 leads to contradiction, irrespective of whether j ∈ [n] or j = n + 1. First, consider the case where j ∈ [n]. Then, (69) gives us
where (72) follows from the assumption that (D 1 w ) −1 B 1 w > (D 2 w ) −1 B 2 w , (73) follows fromŷ 1 +ŷ 2 ≤ κỹ 1 , (74) follows from (70), (75) holds since (ŷ 1 +ŷ 2 ) j = κỹ 1 j , and (76) holds due to (ŷ 1 +ŷ 2 ) j ≥ỹ 1 j . Note that (76) is a contradiction, following from the assumption that κ ≥ 1. Now, consider j = n + 1. Then, from (69) we have
where (77) follows from the observation that the matrix X(y) has a 0 in the (n + 1) th position along its diagonal for any y ∈ D k , (78) follows from the assumption that (D 1 w ) −1 B 1 w > (D 2 w ) −1 B 2 w , (79) follows fromŷ 1 +ŷ 2 ≤ κỹ 1 , (80) holds for the same reason as (77), (81) follows from (70), and (82) holds since (ŷ 1 +ŷ 2 ) n+1 = κỹ 1 n+1 . Note that (82) is a contradiction, following from the assumption that κ ≥ 1.
Since the assumption that κ ≥ 1 leads to contradiction for all j ∈ [n + 1], we have κ < 1, implying thatŷ 1 +ŷ 2 ỹ 1 . Now, note that sinceỹ 1 1, andŷ 1 +ŷ 2 1, we know that (I − X(ŷ 1 ) − X(ŷ 2 ))(D 1 w ) −1 B 1 w and (I − X(ỹ 1 ))(D 1 w ) −1 B 1 w are irreducible nonnegative matrices. Then, sinceỹ 1 0, it follows from (70) and item (iii) in Lemma 4 that ρ((I − X(ỹ 1 ))(D 1 w ) −1 B 1 w ) = 1. Likewise,ŷ 1 0, (69) and item (iii) in Lemma 4 give us ρ((I − X(ŷ 1 ) − X(ŷ 2 ))(D 1 w ) −1 B 1 w ) = 1. Following fromŷ 1 +ŷ 2 ỹ 1 , we have
Applying item (iii) in Lemma 4 to (69), (70), and (83) yields
Clearly, (84) is a contradiction following from the assumption thatŷ exists, and hence,ŷ does not exist. Therefore, the only equilibria in D are the healthy state, (ỹ 1 , 0) and (0,ỹ 2 ). It remains to be shown that the healthy state and (0,ỹ 2 ) are unstable, and that (ỹ 1 , 0) is locally exponentially stable. Since, by assumption, s(B 1 w − D 1 w ) > 0 and s(B 2 w − D 2 w ) > 0, instability of the healthy state follows directly. Moreover, since, by assumption, (D 1 w ) −1 B 1 w > (D 2 w ) −1 B 2 w , from Lemma 8, it follows thatỹ 1 >ỹ 2 . Recall that ρ(I − X(ỹ 1 )(D 1 w ) −1 B 1 w ) = 1. Then, sinceỹ 1 >ỹ 2 , item (iv) in Lemma 4 implies that 