{'fraction_non_alphanumeric': 0.04974664836904888, 'fraction_numerical': 0.02612688694183469, 'mean_word_length': 4.089724647414372, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 18, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling and Analysis of a Coupled SIS Bi-Virus Model

The paper deals with the setting where two viruses (say virus~1 and virus~2) coexist in a population, and they are not necessarily mutually exclusive, in the sense that infection due to one virus does not preclude the possibility of simultaneous infection due to the other. We develop a coupled bi-virus susceptible-infected-susceptible (SIS) model from a 4n-state Markov chain model, where n is the number of agents (i.e., individuals or subpopulation) in the population. We identify a sufficient condition for both viruses to eventually die out, and a sufficient condition for the existence, uniqueness and asymptotic stability of the endemic equilibrium of each virus. We establish a sufficient condition and multiple necessary conditions for local exponential convergence to the boundary equilibrium (i.e., one virus persists, the other one dies out) of each virus. Under mild assumptions on the healing rate, we show that there cannot exist a coexisting equilibrium where for each node there is a nonzero fraction infected only by virus 1; a nonzero fraction infected only by virus 2; but no fraction that is infected by both viruses 1 and 2. Likewise, assuming that healing rates are strictly positive, a coexisting equilibrium where for each node there is a nonzero fraction infected by both viruses 1 and 2, but no fraction is infected only by virus 1 (resp. virus 2) does not exist. Further, we provide necessary conditions for the existence of certain other kinds of coexisting equilibria. Finally, we illustrate our theoretical findings using an extensive set of in-depth simulations.


Introduction
The phenomenon of spreading processes has been a key facet of human civilization. Several manifestations of this phenomenon are witnessed in the present day too, including the spread of opinions in social networks, diseases in contact networks, of viruses in computer networks, of products in markets, etc. Given the various ramifications of such processes, researchers across diverse disciplines such as physics (Newman et al., 2002), ecology (Munster and Fouchier, 2009), epidemiology (Bailey et al., 1975), computer science (Wang et al., 2003), and economics (Bloom et al., 2018) have devoted significant attention to the same. The focus of this paper is on the spread of viruses in human contact networks. The first model to capture the spread of a virus was proposed by Daniel Bernoulli in the 18 th century to calculate the gain in life expectancy at birth if smallpox were to be eliminated as a cause of death (Bernoulli, 1760). As a discipline in its own right, mathematical epidemiology witnessed enormous growth in the 20 th century, with (Bailey et al., 1975;Hethcote, 2000) being some of the key works.
One of the fundamental research objectives in mathematical epidemiology revolves around comprehending the causes behind a virus becoming persistent (or extinct). Leveraging this comprehension enables the design of mitigation (or eradication) strategies. To this end, the literature abounds with several models. The focus of this paper is on the susceptible-infected-susceptible (SIS) model. In particular, we are interested in networked SIS models. Very briefly, in a networked SIS model, a population of individuals is partitioned into subpopulations, called agents. The interconnection between various agents can be represented by a (possibly) directed graph. Supposing that there is a virus prevalent in the population, it can spread both between individuals in a subpopulation and also between individuals belonging to different subpopulations. If none of the individuals in a subpopulation are infected, then this subpopulation is said to be in the susceptible state; otherwise, it is said to be in the infected state. An otherwise susceptible agent getting infected is conditional on this agent coming into contact with infected neighbors and its own infection rate. An infected agent, depending on its healing rate, recovers from the infection, and becomes susceptible to the virus again. A key feature of the SIS model is that recovering from the virus does not confer immunity. Networked SIS models have been extensively studied in the literature; see, for instance, (Fall et al., 2007;Khanafer et al., 2016;Rami et al., 2014;Peng et al., 2010;Wang et al., 2003;Gracy et al., 2020;Paré et al., 2020).
While all of the aforementioned papers provide insights into how to control an epidemic, they have major limitations. More specifically, they do not account for settings where multiple strains of a virus could be simultaneously active within a population. The dynamics in the multi-virus setting are far richer than those in the singlevirus setting; see, for instance, (Castillo-Chavez et al., 1989, 1999Sahneh and Scoglio, 2014;Liu et al., 2019b;Zhang et al., 2022;Janson et al., 2022;Newman and Ferrario, 2013;Chen et al., 2013;Granell et al., 2013;Ahn et al., 2006;Beutel et al., 2012;Xu et al., 2012;Gao et al., 2016;Zhao et al., 2020). More specifically, suppose that there are two viruses, say virus 1 and virus 2, prevalent, then some of the interaction possibilities are as follows: a) competitive, e.g., leprosy and tuberculosis (Sahneh and Scoglio, 2014;Castillo-Chavez et al., 1989;Ye et al., 2021); b) reinforcing (Newman and Ferrario, 2013;Chen et al., 2013); c) weakening (Granell et al., 2013), and d) asymmetric (Ahn et al., 2006). More specifically, in the competitive regime, an agent can be infected either with virus 1 or with virus 2 or neither. In the reinforcing regime, infection with virus 1 (resp. virus 2) increases the possibility of infection with virus 2 (resp. virus 1); e.g., human immunodeficiency virus (HIV) and syphilis (resp. herpes simplex virus type 2 (HSV-2)). The case of the weakening regime is opposite to that of the reinforcing regime, e.g., epidemic process spreading on a contact network, and an information awareness process diffusing in the network of virtual social contacts between the same individuals. In the asymmetric regime, infection by virus 1 increases the possibility of simultaneous infection by virus 2, but infection by virus 2 decreases the possibility of simultaneous infection by virus 1, e.g., the spread of malicious pathogens and growth of immune cells in living organisms. Due to several interaction possibilities between the viruses, analysing such models becomes mathematically more challenging than the single-virus setting.
Multi-virus models that admit the possibility of an agent being infected with more than one virus at the same time are broadly referred to as coupled bi-virus models. In this paper, our focus is on the development and analysis of a coupled networked SIS bi-virus model. The authors in  proposed such a model, but it accounted only for undirected graphs, and each of the nodes are restricted to have the same healing and infection rates with respect to each of the viruses. A probabilistic coupled bi-virus SIS model was proposed and studied in (Xu et al., 2012), but the authors treated infection and recovery from each virus as probabilities, which implies that those are constrained to stay between zero and one. Different from (Xu et al., 2012), the model we propose in this paper admits nonnegative infection rates and positive healing rates larger than one for each virus. Moreover, we provide a richer analysis of the various equilibria of the co-operative viruses regime. Recently, a coupled bi-virus SIS model that accounts for only a single population node has been proposed in (Gao et al., 2016;Zhao et al., 2020). Overcoming this limitation, the coupled bi-virus model that we propose admits an arbitrary but finite number of nodes. The two viruses may spread through possibly different directed contact graphs. An agent could be either infected by virus 1 or by virus 2 or, at the same time, by both (i.e., viruses 1 and 2) or by neither. We classify the equilibria into the following classes: a) the healthy state (both viruses are eradicated), b) single-virus endemic equilibrium (the endemic equilibrium corresponding to that of virus 1 (resp. virus 2) if only virus 1 (resp. virus 2) were prevalent in the population), c) boundary equilibria (one virus is eradicated, the other one is persistent, fraction of population infected by both is zero), and d) coexisting equilibria (both viruses simultaneously infect possibly same fractions of the population).
Our main contributions are as follows: i) Derivation of the coupled bi-virus model from a 4 nstate Markov chain model; see equations (8)-(10) in Section 2.
ii) We provide a sufficient condition which ensures that, irrespective of the initial state of the network (i.e., healthy or sick), both viruses eventually die out; see Theorem 1.
iii) We provide a sufficient condition for the existence, uniqueness, and asymptotic stability of a single-virus endemic equilibrium; see Theorem 2. iv) We identify a sufficient condition for local exponential stability of the boundary equilibria; see Theorem 3.
v) Assuming that the healing rates are strictly positive, we show that a point in the 3n-dimensional state space where for each node there is a nonzero fraction infected only by virus 1 (resp. virus 2) but no fraction that is infected by both viruses 1 and 2 cannot be an equilibrium point; see Proposition 3. Likewise, under mild assumptions on the healing rates, a point where for each node there is a nonzero fraction infected by both viruses 1 and 2, but no fraction is infected only by virus 1 (resp. virus 2) cannot be an equilibrium; see Proposition 4.
Additionally, we provide a necessary, and sufficient, condition for the healthy state to be the unique equilibrium of the system; see Corollary 1. Assuming both viruses pervade the network, we establish a lower bound on the number of equilibria for the coupled bi-virus system; see Corollary 2. We identify multiple necessary conditions for local exponential convergence to the boundary equilibria; see Proposition 2. We establish a necessary condition for the existence of a coexisting equilibrium wherein the fraction of each node infected by only virus 2 is zero and the rest (i.e., the fraction infected by only virus 1 and the fraction infected by both viruses 1 and 2) are strictly positive; see Proposition 5. Finally, we identify a necessary condition for the existence of an equilibrium where for each node the fraction of subpopulation infected only by virus 1 is strictly positive, the fraction infected only by virus 2 is strictly positive, and the fraction infected by both viruses 1 and 2 is strictly positive; see Proposition 6. Some of the material in this paper was partially presented earlier in (Paré et al., 2018b); the present paper provides a more comprehensive treatment of the work, and considers a more general model. Specifically, the paper provides i) an expansion of the model to the case of virusdependent coupling parameters (i.e., different k for each virus k) that can be greater than 1; ii) complete proofs of all the results; iii) a derivation of the coupled bi-virus model from a 4 n -state Markov chain model; iv) stability results for the boundary equilibria; v) existence results for coexisting equilibria; and vi) additional illustrative simulations in Section 7, none of which were included in (Paré et al., 2018b).
The paper is organized as follows. We conclude this section by introducing all the notations used. The derivation of the coupled bi-virus model from a 4 n -state Markov chain model is provided in Section 2. The problems of interest and standing assumptions are formally stated in Section 3. The main results for the model developed in Section 2 are split across the next three sections: analysis of the DFE is given in Section 4; persistence of a virus in the population is given in Section 5; and analysis of various coexisting equilibria are provided in Section 6. The theoretical findings are illustrated in Section 7. A summary of the results of the paper, and some questions of possible interest to the wider community are given in Section 8.
Notation: For any positive integer n, we use [n] to denote the set {1, 2, . . . , n}. We use 0 and 1 to denote the vectors whose entries all equal 0 and 1, respectively, and I to denote the identity matrix, while the dimensions of the vectors and matrices are to be understood from the context. For any vector x ∈ IR n , we use x to denote its transpose and diag(x) or X to denote the n × n diagonal matrix whose ith diagonal entry equals x i . The notation 1 a=b is used as an indicator function which takes value one if a equals b and zero otherwise. For 1 A=b , where A is a matrix and b is a scalar, the result is a binary matrix of the same dimensions as A with entries 1 aij =b . For any two sets A and B, we use A \ B to denote the set of elements in A but not in B. For any two real vectors . Likewise, for any two real matrices A, B ∈ R n×m , we write For a real square matrix M , we use s(M ) to denote the largest real part among the eigenvalues of M . A real square matrix A is said to be Metzler if all of its off-diagonal entries are nonnegative. A real square matrix A is said to be a Z-matrix if all of its off-diagonal entries are nonpositive. A Z-matrix is an M-matrix if all its eigenvalues have nonnegative real parts. Furthermore, if an M-matrix has an eigenvalue at the origin, then we say that it is singular; if each of its eigenvalues have strictly positive parts, then we say that it is nonsingular.

The Model
Consider two viruses spreading over a network of n agents. Each agent may be infected with either or both viruses at the same time. Specifically, each agent can be infected if one of its neighbors is infected. The neighbor relationships among the n agents are described by an nvertex directed graph. A directed edge from node j to node i means that agent i can be infected by agent j, i.e., agent j is a neighbor of agent i. We use N i to denote the set of neighbors of agent i. The two viruses may spread through different routes in the network. We use N k i to denote the set of neighbors of agent i from which virus k For each virus k ∈ {1, 2}, each agent i has its curing rate δ k i and infection rates β k ji when i ∈ N k j . The former means that if agent i is infected by virus k, it is cured with rate δ k i , and the latter means that if agent i is infected by virus k and its neighbor j is not, agent i can infect agent j at rate β k ji . The two viruses can simultaneously infect the same node, but not independently. Specifically, they are coupled in the following manner. Let i and j be any pair of integers in [n] such that agent j is a neighbor of agent i. If agent j has been infected by only one virus, say virus 1, and agent i is infected by the other virus, virus 2, then, irrespective of whether (or not) it is infected by virus 1, agent i can infect agent j with virus 2 at a rate 1 β 2 ji . See Figure 1 for a depiction of the model. It is worth noting that there is no transition link from state I 1,2 (infected by both viruses) to state S (healthy state), as the probability that the two viruses are cured at the same time is zero. We call k the coupling parameter between the viruses, and assume that each k takes a nonnegative value. If k ∈ (0, 1), that means if a node is infected with virus k, it is less susceptible to the other virus. If k > 1 that means if a node is infected with virus k, it is more likely to become infected with the other virus. If k = 1 for all k, then the two viruses are independent. If k = 0 for all k, then the two viruses are competitive. We will discuss these last two special cases in-depth in Section 4. Let . The spread of the two viruses across the population can be represented by a two-S I 1 Figure 1: The possible states and transitions for node i. State S stands for being healthy but susceptible. States I 1 , I 2 , and I 1,2 stand for being infected by virus 1, virus 2, and both virus 1 and virus 2, respectively. layer graph, where the vertices of the graph correspond to the population nodes. Each layer contains a set of directed edges, E k , specific to virus k; there exists a directed edge from agent j to agent i in E k if, assuming agent j is infected with virus k, it can directly infect agent i. Note that there is a one-to-one correspondence between the no- It is worth emphasizing that the two viruses may spread along different routes, that is, the layers corresponding to B 1 and B 2 are not necessarily the same. We call the layers corresponding to B 1 and B 2 as the spreading graphs of viruses 1 and 2, respectively.
For completeness, we provide here a full description of the 4 n -state Markov model. Each state of the chain, Y k (t), corresponds to a quaternary-valued string s of length n, where s i = 0, s i = 1, s i = 2, or s i = 3 indicate that the ith agent is either susceptible, infected with virus 1, infected with virus 2, or infected with both viruses 1 and 2 (indicated by the number 3), respectively. The generator matrix (Norris, 1998), Q, is defined by for i ∈ [n]. Here virus 1 and virus 2 are propagating over a network whose infection rates are given by β 1 ij and β 2 ij , respectively (nonnegative with β 1 ii = β 2 ii = 0, ∀j), δ 1 i and δ 2 i are the respective healing rates of the ith agent, and, again, s i = 0, s i = 1, s i = 2, or s i = 3 indicate that the ith agent is either susceptible, infected with virus 1, infected with virus 2, or infected with both viruses 1 and 2, respectively. The state vector y(t) is defined as with 4 n k=1 y k (t) = 1. The Markov chain evolves as is the random variable representing whether the ith agent is susceptible or infected with virus 1, or 2, or both. Then where the ith columns of M 1 , M 2 , M 3 indicate the states in the Markov chain where agent i is infected with virus 1, virus 2, and both (all the ternary strings where s i = 1, s i = 2, and s i = 3), respectively, that is, where M ∈ R 4 n ×n has rows of lexicographically-ordered ternary numbers, bit reversed. 1 Therefore, v 1 and v 3 i (t) reflect the summation of all probabilities where s i = 1, s i = 2, and s i = 3. Note that the first state of the chain, which corresponds to s = 0, the healthy state, for δ 1 i , δ 2 i > 0 ∀i, is the absorbing or sink state of the chain. That is, once in the healthy state, the Markov chain will never escape it. Moreover, since the healthy state is the only absorbing state, the system will converge to it with probability one (Norris, 1998).
Remark 1 We consider two special cases of the model. First, let 1 = 2 = 0, and δ 1 i + δ 2 i > 0 for all i ∈ [n]. Then, the system defined by (8)-(10) simplifies tȯ In this case, since the matrix D 1 + D 2 is positive definite, it follows that z(t) converges to 0 exponentially fast, and thus the system will eventually become a competitive bivirus model which has been studied in Santos et al., 2015;Liu et al., 2016Liu et al., , 2019b.
In the second case, we let 1 = 2 = 1. To proceed, we define y 1 i (t) = x 1 i (t) + z i (t) and y 2 i (t) = x 2 i (t) + z i (t) for each i ∈ [n], which represents the total probabilities of agent i being infected by viruses 1 and 2, respectively. From (5)-(7), the dynamics of y 1 i and y 2 i arė which are two independent single SIS dynamics. Therefore, the system defined by (8)-(10) subsumes both the single SIS virus (two single, independent viruses) and the competitive SIS virus models.

Problem Formulation
In this section, we formally state the problems of interest, and the key assumptions needed for ensuring that the model introduced in Section 2 is well defined.

Problem Statements
With respect to the model in (8)-(10), we ask the following seven questions: (i) Identify a sufficient condition under which, irrespective of the initial state, the dynamics converge asymptotically to the healthy state.
(ii) Identify a sufficient condition for virus k, such that for any x k (0) = 0 the dynamics asymptotically converge to the single-virus endemic equilibrium of virus k, for k = 1, 2.
(iii) Identify a sufficient condition for local exponential convergence to the boundary equilibrium of virus k, for k = 1, 2.
(iv) Provide necessary condition(s) for local exponential convergence to the boundary equilibrium of virus k, for k = 1, 2.
(vii) Identify a necessary condition for the existence of the coexisting equilibria (x 1 ,x 2 ,ẑ).

Key Assumptions and Preliminaries
We make the following assumptions on the model to ensure that it is well defined.
Assumption 1 For all i ∈ [n], we have Assumption 2 For all i ∈ [n], we have δ 1 i , δ 2 i ≥ 0. The matrices B 1 and B 2 are nonnegative and irreducible.
Assumption 1 guarantees that the initial infection level with respect to each virus k (k ∈ [2]) in each node i ∈ [n] lies in the set [0, 1], whereas Assumption 2 ensures that the healing and infection rates are nonnegative, and that the spreading graphs for virus 1 and 2 are strongly connected.
Define the set x 1 + x 2 + z ≤ 1} The following lemma establishes that the set D is positively invariant with respect to the system (8)-(10).
Proof: See the Appendix. Since x 1 i , x 2 i , and z i denote the probabilities of sickness of agent i, or fractions of group i, infected by viruses 1, 2, and both 1 and 2 simultaneously, respectively, and 1 − x 1 i − x 2 i − z i denotes the probability of agent i, or fraction of group i that is healthy, it is natural to assume that their initial values are in the interval [0, 1], since otherwise the values will be devoid of any physical meaning for the spread model considered here.

Analysis of the Disease-Free Equilibrium
In this section, we analyze the system defined by (8)-(10). It is easy to see that (0, 0, 0) is an equilibrium of the system defined by (8)-(10). We call it the disease-free equilibrium, or the healthy state. In this section, we focus on identifying conditions under which the healthy state is stable. The following proposition provides a necessary and sufficient condition for local exponential convergence to the healthy state. Proof: From (15), we have Thus, from (Khalil, 2002, Theorem 4.15 and Corollary 4.3), the healthy state is locally exponentially stable if, and only if, s(−D 1 + B 1 ) < 0, s(−D 2 + B 2 ) < 0, and δ 1 i + δ 2 i > 0 for all i ∈ [n]. Note that if s(−D 1 + B 1 ) > 0 or if s(−D 2 + B 2 ) > 0, then s(J(0, 0, 0)) > 0. The claim on instability then follows from (Khalil, 2002, Theorem 4.7).
Note that, on the one hand, the guarantees provided by Proposition 1 are limited in the sense that they concern trajectories that originate in a small neighborhood of the healthy state. On the other hand, no restrictions, besides nonnegativity, are imposed on k , k = 1, 2. Simulations, as we see in Section 7.2, indicate that the region of attraction for the healthy state depends on the choices of k , k = 1, 2. In particular, if the initial state of system (8)-(10) is very close to the healthy state, then, even for a larger value of k , k = 1, 2, the dynamics converge to the healthy state. If the initial state of system (8)-(10) is not too close to the healthy state, then for large values of k , the dynamics do not converge to the healthy state; see Figure 8 in Section 7.2.
The following theorem guarantees global convergence to the healthy state, but with the following caveats: i) the speed of convergence is slower, and ii) more restrictions on k , k = 1, 2, have to be imposed.
Proof: See the Appendix. Theorem 1 answers Question (i) raised in Section 3.1.

Persistence of Viruses
We call an equilibrium (x 1 ,x 2 ,ẑ) an endemic equilibrium if it is not the healthy state, (0, 0, 0). It turns out that if either (or both) of the spectral abscissa conditions in Theorem 1 are violated, then at least one of the viruses pervades the population. We detail the same in the rest of this section.

Existence, Uniqueness and Stability of the Single Virus Endemic Equilibria
We consider the scenario in which either s(B 1 − D 1 ) or s(B 2 − D 2 ) is greater than zero. Without loss of generality, we assume that s(B 1 − D 1 ) > 0, 1 ∈ [0, 1], and s(B 2 − D 2 ) ≤ 0. We denote byx 1 (resp.x 2 ) the singlevirus endemic equilibrium corresponding to virus 1 (resp. virus 2). We have the following result.
Proof: See the Appendix. Theorem 2 establishes the existence, uniqueness, and asymptotic stability of the single-virus endemic equilibrium corresponding to virus 1. An analogous result holds for virus 2; details of which are omitted in the interest of space. Theorem 2 answers Question (ii) raised in Section 3.1.
Combining Theorems 1 and 2, we obtain a necessary and sufficient condition for the healthy state to be the unique equilibrium of the coupled bi-virus system, as stated below.

Both Viruses Pervading the System
Note that Theorem 2 accounts for the case where exactly one of the viruses pervades the system, or, in other words, exactly one of the spectral abscissa conditions in Theorem 1 is violated. However, what happens when both the spectral abscissa conditions in Theorem 1 are violated?

Proof: See the Appendix.
The equilibria of the kind (x 1 , 0, 0) and (0,x 2 , 0) are hereafter referred to as the boundary equilibria. Note that x 1 andx 2 are asymptotically stable in the single virus (i.e., one of the two viruses has died out) systems corresponding to virus 1 and virus 2, respectively; when both viruses pervade the network, the stability of (x 1 , 0, 0) and (0,x 2 , 0), in even the local (let alone global) sense, is not guaranteed. As such, in the rest of this section we will focus on identifying a sufficient condition (resp. some necessary conditions) for local exponential stability of the boundary equilibria.
We need the following assumption, which is slightly stronger than Assumption 2.
Assumption 3 For all i ∈ [n], we have δ 1 i , δ 2 i > 0. The matrices B 1 and B 2 are nonnegative and irreducible.
The following proposition provides necessary conditions for local exponential stability of the equilibrium (x 1 , 0, 0). Proposition 2 Consider system (8)-(10) under Assumptions 1 and 3. Suppose that 1 = 2 = , and that s(−D 1 + B 1 ) > 0. The equilibrium point (x 1 , 0, 0) is locally exponentially stable only if each of the following conditions are satisfied Proof: See the Appendix. Proposition 2 answers Question (iv) raised in Section 3.1.
Remark 2 Note that, in general, there is a gap between the sufficient condition in Theorem 3 and the necessary conditions in Proposition 2. Note, however, that if k = 0 for k = 1, 2, the sufficient condition in Theorem 3 and the necessary conditions in Proposition 2 coincide to yield a necessary and sufficient condition for local exponential convergence to (x 1 , 0, 0). Further, by Assumption 3, δ 1 i > 0 and δ 2 i > 0 for each i ∈ [n]. Hence, if k = 0 for k = 1, 2, then condition ii) in both Theorem 3 and Proposition 2 is always satisfied. As a consequence, condition i) in Theorem 3 becomes a necessary and sufficient condition, which is consistent with (Ye et al., 2022, Theorem 3.10).
Proof: See the Appendix.
Note that by suitably changing the notations of Theorem 3, Proposition 2, and Corollary 3, we can obtain a sufficient condition, necessary conditions, and a condition for instability, respectively, for the boundary equilibrium (0,x 2 , 0).

Analysis of Coexisting Equilibria
Equilibria of the kind (x 1 ,x 2 ,ẑ) with at least any two ofx 1 ,x 2 ,ẑ being non-zero vectors, are referred to as coexisting equilibria. In this section, we show that certain kinds of coexisting equilibria cannot exist, and we focus on identifying necessary conditions for the existence of certain other kinds of coexisting equilibria.

Impossibility of Existence of a Certain Kind of Coexisting Equilibria
It is well known that for the competitive bivirus case (i.e., k = 0 for k = 1, 2), one of the possible equilibria is the so-called coexistence equilibrium, where separate fractions of each node is infected by virus 1 and virus 2 (Liu et al., 2019b;Ye et al., 2022;Janson et al., 2022). More formally, these are equilibria of the kind (x 1 ,x 2 ) with 0 x 1 1, 0 x 2 1 andx 1 +x 2 1. It is natural to wonder whether similar equilibria could exist even when k > 0 for k = 1, 2. It turns out, however, that such coexistence equilibria do not exist for coupled bivirus systems. The following proposition formalizes this. Proposition 3 Consider system (8)-(10) under Assumptions 1 and 2. Suppose that k > 0 for some k ∈ [2]. There does not exist a coexisting equilibrium of the form (x 1 ,x 2 , 0) withx 1 ,x 2 > 0.
Proof: See the Appendix.
In a similar vein to Proposition 3, the following proposition states that the coupled bivirus system cannot have an equilibrium where a fraction of each node is infected by both viruses at the same time, but that no fraction of any node is infected only by one of the viruses. Proposition 4 Consider system (8)-(10) under Assumption 3. There does not exist an equilibrium of the form (0, 0,ẑ), whereẑ > 0.
Proof: See the Appendix.
Taken together, Propositions 3 and 4 answer, in the negative, Question (v) raised in Section 3.1. In so doing, they restrict the set of possible endemic equilibria for the coupled bivirus system.

Necessary Conditions for Existence of Certain Kinds of Coexisting Equilibria
While Subsection 6.1 has dealt with the impossibility of existence of certain kinds of coexisting equilibria, in this subsection we are interested in identifying some necessary conditions for the existence of certain other kinds of coexisting equilibria. We begin by presenting a necessary condition for the existence of a coexisting equilibrium where for each node the fraction infected only by virus 1 is non-zero, by both viruses 1 and 2 is non-zero, but only by virus 2 is zero. Proposition 5 Consider system (8)-(10) under Assumption 3. Suppose that i) ρ((D 1 ) −1 B 1 ) > 1; and ii) 1 = 2 = , with ∈ (0, 1).
Proof: See the Appendix.
Next, we present a necessary condition for the existence of a coexisting equilibrium where each node has i) a fraction that is infected only by virus 1; ii) a fraction that is infected only by virus 2; and iii) a fraction that is infected by both viruses 1 and 2. To this end, we need the following assumption.

Assumption 4
The healing and infection rates are the same for each virus. That is, δ 1 i = δ 2 i for all i ∈ [n], and β 1 ij = β 2 ij for all i = j ∈ [n] and (i, j) ∈ E.

Simulations
In this section, we present a comparison of the 4 n -state Markov model in (1)-(4) to (8)-(10) via simulation, and also provide a set of simulations of various coupled virus models.

Comparison to Full Probabilistic Model
We compare the model in (8)-(10) to the full probabilistic 4 n -state model in (1)-(4) via simulations to illustrate the effectiveness of the approximation. We set 1 = 2 = 3, and use line graphs, star (hub-spoke) graphs, and complete graphs. For examples of each type of graph, see Figure 2. All adjacency matrices for these graphs are symmetric and binary-valued, and both viruses spread over the same graph. In the star graph, the central node is the first agent. Each simulation was run for 10,000 time steps (final time T = 10, 000), with three initial conditions: 1) the first node is infected by virus 1 and the second node is infected by virus 2, 2) the first node is infected by virus 1, the second node is infected by virus 2, and the third node is infected by both virus 1 and virus 2, x 2 (0) = [0 1 0 · · · 0] z(0) = [0 0 1 0 · · · 0] ; and 3) the first node is infected by virus 1, the second node is infected by virus 2, and the remaining nodes are infected by both virus 1 and virus 2, In these tests we explore identical homogeneous viruses, (β, δ) = (β 1 , δ 1 ) = (β 2 , δ 2 ) = (β 3 , δ 3 ). The (β, δ) pairs are [(0.1, 1), (0.215, 1), (0.464, 1), (0.5, 0.5), (1, 0.464), (1, 0.215), (1, 0.1)], and the numbers of agents are n = 4, 6. We limited simulations to these n values since mean field approximations are typically worse for small values of n and there is a computational limitation due to the size of the 4 n -state Markov model.
The results are given in Figures 3-5 in terms of the 2-norm of the difference between the states of (8)-(10) at the final time ([x 1 (T ); x 2 (T ); z(T )]), and the means of the three states in the 4 n -state Markov model at the final time ([v 1 (T ); v 2 (T ); v 3 (T )] as defined by (4)).
The accuracy of the approximation appears to be very similar to the single virus case (Mieghem et al., 2009;Paré et al., 2018a) and to the two-virus case (Liu et al., 2019b). Since the model in (8)-(10) is an upper bounding approximation, the results show that the two models converge to the healthy state for the smaller values of β δ , resulting in small errors between the two models. For many of the larger values of β δ , the model in (8)-(10) again performs quite well since it is at an epidemic state and the 4 n -state Markov model does not appear to reach the healthy state in the finite time considered in the simulations (T = 10, 000). Therefore, for certain values of β δ and certain time scales, the model in (8) appears, in most cases, to be at or close to the healthy state while the model in (8)-(10) is at an epidemic state, resulting in large errors.

Illustrative Examples
Virus 1 is depicted by the color red (r), virus 2 is depicted by the color blue (b), and the state of being infected with both states, z(t), is depicted by the color green (g). For all i ∈ [n], the color at each time t for node i is given by where s i (t) = x 1 i (t) + x 2 i (t) + z i (t). When x 1 i (t) + x 2 i (t) + z i (t) = 0, the color is black, indicating completely healthy, susceptible. These color variations are used to facilitate the depiction of the parallel equilibria (x 1 = αx 2 ), a behavior that is exhibited by the 1 = 2 = 0 case (see (Liu et al., 2016)) and are illustrated by all the nodes converging to the same color. For all i ∈ [n], the diameter of node i is given by with d 0 being the default/smallest diameter and r 0 being the scaling factor depending on the total sickness of node i. Therefore, the color indicates the type of virus(es) each agent has and the diameter indicates the degree of sickness for each agent. The graph structure is as follows: where y i (t) ∈ R 2 is the position of node i, withr = .35.
We consider a network of n = 15 nodes. The binary matrices A 1 and A 2 are populated in correspondence, respectively, to the black and green edges depicted in Figure 6. The initial condition of the network is shown in Figure 6. The elements of diagonal matrices D 1 and D 2 are chosen uniformly at random from [0, 1]. Let b 1 be a vector whose elements are chosen uniformly at random from [0, 1]. Let B 1 = diag(b 1 /6)A 1 . Vector b 2 is chosen analogous to b 1 . Then, let B 2 = diag(b 2 /11)A 2 . Let 1 = 2 = , and fix = 0.5. Given that the matrices involved are of dimension 15 × 15, in the interest of space, we refrain from providing exact values; the exact values used along with the Matlab code used are available here 2 . With these choices of model parameters, it turns out that s(−D 1 + B 1 ) = −0.0037 < 0 and s(−D 2 + B 2 ) = −0.0013 < 0. Consistent with Proposition 1 (and since ∈ (0, 1), also consistent with Theorem 1), both the viruses die out; see Figure 7a. Next, with the same model parameters and initial condition as before, we set = 2. Again, consistent with Proposition 1, the dynamics converge to the healthy state, albeit the rate of convergence is slower than that with = 0.5; see Figure 7b. Once again, with the same model parameters and initial condition as before, we set = 1000. Although s(−D 1 + B 1 ) = −0.0037 < 0 and s(−D 2 + B 2 ) = −0.0013 < 0, as a consequence of the effect of the large value of , the dynamics do not die out; see Figure 8.
For the next simulation, D 1 , D 2 , b 1 , and b 2 are the same as before. Let B 1 = diag(b 1 )A 1 , and let B 2 = diag(b 2 /11)A 2 . Choose = 0.5. With these choices of model parameters, it turns out that s(−D 1 +B 1 ) = 1.4904 > 0 and s(−D 2 + B 2 ) = −0.0013 < 0. Consistent with Theorem 2, virus 2 dies out, virus 1 becomes endemic in the population, and no fraction of any node is infected by both viruses 1 and 2 simultaneously; see Figure 9.
Hence, in line with the result in Theorem 3, the dynamics converge to (x 1 , 0, 0) exponentially fast; see Figure 10.

Conclusion
The paper has tackled the problem of simultaneous infection of an individual (resp. subpopulation) by possibly two viruses. We derived a coupled bi-virus model from a 4 n -state Markov chain model. We identified a condition 2 https://github.com/philpare/coupled bivirus    that leads to the extinction of both viruses; likewise a condition that causes one of the viruses to become endemic in the population. Subsequently, we provided a sufficient condition and two necessary conditions for local exponential convergence of boundary equilibria. With respect to coexistence equilibria, we conclusively ruled out the existence of the following types of coexisting equilibria: i) a point in the state space where for each node there is a non-trivial fraction infected only by virus 1, a non-trivial fraction infected only by virus 2, but no fraction that is infected by both viruses 1 and 2; and ii) a point in the state space where for each node there is a fraction that is infected simultaneously by both viruses 1 and 2, but no fraction is infected only by virus 1 (resp. virus 2). We provided necessary conditions for the existence of certain other kinds of coexisting equilibria. Problems of further interest include a) identifying condition(s) that would lead to global asymptotic stability of the boundary equilibria, and b) stability (local or global) guarantees for various coexisting equilibria. A further extension would be to time-varying graphs.