key: cord-0175252-a1xedm8e
authors: Ye, Mengbin; Anderson, Brian D.O.; Liu, Ji
title: Convergence and Equilibria Analysis of a Networked Bivirus Epidemic Model
date: 2021-11-15
journal: nan
DOI: nan
sha: 6ef0d3ec6033fd4db91e8dc0626db3e571629b89
doc_id: 175252
cord_uid: a1xedm8e

This paper studies a networked bivirus model, in which two competing viruses spread across a network of interconnected populations; each node represents a population with a large number of individuals. The viruses may spread through possibly different network structures, and an individual cannot be simultaneously infected with both viruses. Focusing on convergence and equilibria analysis, a number of new results are provided. First, we show that for networks with generic system parameters, there exist a finite number of equilibria. Exploiting monotone systems theory, we further prove that for bivirus networks with generic system parameters, then convergence to an equilibrium occurs for all initial conditions, except possibly for a set of measure zero. Given the network structure of one virus, a method is presented to construct an infinite family of network structures for the other virus that results in an infinite number of equilibria in which both viruses coexist. Necessary and sufficient conditions are derived for the local stability/instability of boundary equilibria, in which one virus is present and the other is extinct. A sufficient condition for a boundary equilibrium to be almost globally stable is presented. Then, we show how to use monotone systems theory to generate conclusions on the ordering of stable and unstable equilibria, and in some instances identify the number of equilibria via rapid simulation testing. Last, we provide an analytical method for computing equilibria in networks with only two nodes, and show that it is possible for a bivirus network to have an unstable coexistence equilibrium and two locally stable boundary equilibria.

1. Introduction. Mathematical models of epidemic spreading processes have been of interest to the broad scientific community for decades [23] , and have recently come into the limelight as a result of the ongoing COVID-19 pandemic [35] . In the context of infectious diseases, such models are of interest to predict the dynamics of the disease and the course of an outbreak. One may seek to draw conclusions on whether the disease will eventually disappear or become endemic, examine the impact of key characteristics such as infection and recovery rates in shaping the epidemic, and design control strategies to reduce or stop the spread [21, 22, 23, 18, 39] .

Among the different modeling frameworks, compartmental models have become especially popular; different compartments indicate particular health states and each individual in a large population may move between compartments due to the infectious nature of the disease. A classical compartmental model is the Susceptible-Infected-Susceptible (SIS) model, in which each individual is either healthy and susceptible (S) to the disease, or infected (I) but can recover with no immunity [23] . If immunity is permanently acquired after recovery, then the Susceptible-Infected-Removed (SIR) framework is used [23] . For a particular compartmental framework, both stochastic and deterministic models exist, and are often related by mean-field approximation [21] .

Although not as realistic in capturing the probabilistic nature of infectious disease transmissions, deterministic models remain popular as they offer a balance in terms of analytical tractability, simulation implementation and modeling accuracy.

More recently, attention has grown on models that study the spread of two or more diseases/viruses within the same population (the literature most commonly refers to these as bi-or multi-virus models). Depending on the application of interest, the viruses may be assumed to compete against, reinforce, or weaken one another [26, 20, 6] . Competitive virus models assume that if an individual is infected by one virus, then they cannot be infected by any other virus, and are especially popular. Different variants of SIS-or SIR-like models also have broad applications beyond infectious diseases, including competing ideas, decisions, and internet memes [26, 36, 13] . This paper focuses on the "bivirus model" in the SIS framework, which considers two competing viruses, called virus 1 and virus 2, spreading across a network of interconnected populations on possibly different network structures [17, 24, 26, 37, 28] . Such a scenario may reflect two strains of a disease spreading in a population, such as gonorrhea and a strain of partially drug-resistant gonorrhea [2, 3] . To discuss existing works, we first introduce the reproduction number of virus 1 and virus 2 as R 1 and R 2 , respectively. In general, the infection and recovery rate parameters are assumed to be heterogeneous among the nodes of the network, and then R i for i ∈ {1, 2} is a nontrivial but computable function of the recovery and infection rates. For homogeneous parameters, R i simplifies to the ratio of the infection and recovery rate for virus i, multiplied by the spectral radius of the adjacency matrix of the virus i spreading network. In the case that only virus i, for i ∈ {1, 2}, is present in the population, R i ≤ 1 or R i > 1 determines if virus i asymptotically disappears from or persists in the population, respectively. Much of the analysis in existing works focuses on the assumption that R i ≤ 1 for at least one of i = 1 or i = 2; in this case, either one or both viruses will be eliminated asymptotically [17, 37, 28] . The works of [26, 28, 29] consider R i > 1 for both i = 1, 2. Assuming homogeneous infection and recovery rates among individuals for each virus, but allowing generic network structures, [26] deals only with local stability of equilibria. In contrast, [28, 29] assume general infection rates but homogeneous recovery rates, and focus on global stability of specific equilibria of interest. This paper considers the deterministic bivirus model, with the two viruses spreading on possibly different network structures across the same set of nodes (populations). We allow heterogeneous infection and recovery parameters, which yields new dynamical phenomena, such as the possibility of an unstable coexistence equilibrium, where both viruses are present. We establish several novel results on the system dynamics and associated equilibria.

1. We formally prove that a bivirus system with generic values of infection and recovery parameters has a finite number of equilibria. Prior to our work, the finiteness of equilibria has only been proved for special parameter values [17, 12] , or assumed to be true with no rigorous arguments [26] . We also prove that one can conduct equilibria analysis by assuming unity recovery rates without loss of generality, thus reducing the parameter space dimension. The work of [29, 28] made a similar assumption but without justification. 2. We use a coordinate transformation to demonstrate that the bivirus system has monotone trajectories, simplifying the proof in [28] . Going beyond [28] , we explicitly connect the bivirus dynamics to the monotone systems literature [8, 31, 32] . We significantly extend the existing results by providing a general conclusion on the limiting dynamical behaviour of the bivirus net-worked system: a generic bivirus system converges to an equilibrium among the finite number of equilibria for all initial conditions except possibly a set of measure zero. Thus, no chaos is possible, and limit cycles (if they exist) are nonattractive. Our result differs from [28, 29, 17] , which establish sufficient conditions for a particular equilibrium to be (almost) globally attractive. By further exploiting the literature, we establish several conclusions on the ordering of stable and unstable equilibria. 3. We identify a class of system parameters that yields a connected set containing an infinite number of equilibria, comprising an interval of a straight line, in which both viruses coexist, covering a much broader set of parameters than [17] . We term such a connected set a "line of equilibria" for brevity, and show that this line of equilibria is locally exponentially attractive. Then, simple necessary and sufficient conditions are given for checking whether boundary equilibria, in which one virus is present and the other is absent, are locally stable or unstable. A sufficient condition for one of the boundary equilibria to be almost globally attractive is provided. 4. For networks with only two nodes, an analytic method to compute coexistence equilibria. We then provide numerical examples showcasing the different stability and instability configurations for the boundary and coexistence equilibria. We report a system which has an unstable coexistence equilibrium, and two locally stable boundary equilibria with nontrivial regions of attraction, highlighting the nontrivial equilibria and convergence properties of networked bivirus systems. The paper is organized as follows. Section 2 details the mathematical preliminaries and presents the bivirus model. The main results are detailed in section 3, and conclusions are presented in section 4.

In this section, we introduce notation, and specify the nature of the system under consideration. The set-up is drawn from that in [17] , and we also review the key conclusions of [17, 26, 28, 29] . This section also includes a review of some key results in monotone systems theory [32, 31] , a tool we shall use in this paper.

2.1. Notation and related. Given a natural number n, define the index set [n] = {1, 2, . . . , n}. For real vectors x, y ∈ R n , with entries x i , y i , i ∈ [n], we write x ≥ y iff x i ≥ y i ∀i, we write x > y iff x i ≥ y i ∀i and x = y, and we write x y iff x i > y i ∀i. For matrices A, B of the same dimension, the notation A ≥ B, A > B, A B mean the same thing as the corresponding inequalities relating vec(A), vec(B); we say that A is nonnegative if A ≥ 0 n×n . The n-vectors of all 1's and 0's are denoted by 1 n , 0 n respectively, and the n-dimensional identity matrix is denoted I n . The sets R n ≥0 and R n >0 denote {x ∈ R n : x ≥ 0 n } and {x ∈ R n : x 0 n } respectively. The set {x ∈ R n , 0 n ≤ x ≤ 1 n } will be denoted by Ξ n .

Suppose A is a square matrix. Then ρ(A) and s(A) will denote the spectral radius of A and the greatest real part of any eigenvalue of A, respectively. It is Hurwitz if s(A) < 0. We say that A is reducible iff there is a permutation matrix P such that P AP is block upper triangular; otherwise A is said to be irreducible. When a nonnegative A is irreducible, and Ax = y for x > 0 n , y > 0 n , y cannot have a zero entry in every position where x has a zero entry; equivalently, A is reducible if there exist x > 0 n , y > 0 n with the same set of zero entries. A square matrix A is a Metzler matrix if all off-diagonal entries are nonnegative. For an irreducible Metzler A, and by an extension of the Perron-Frobenius theorem [9] , s(A) is a simple eigenvalue and the only eigenvalue with this real part, and the corresponding eigenvector can be taken to have all positive entries, while the eigenvectors corresponding to other eigenvalues do not have this property. A square matrix A is an M -matrix if −A is Metzler and all eigenvalues of A have positive real parts except for any at the origin; we say A is a singular or nonsingular M -matrix if it has an eigenvalue at the origin or if its eigenvalues have strictly positive real parts, respectively [9] . Further key properties, detailed in [ 

. . , n} as the vertex set, E ⊆ V × V the edge set and A a nonnegative n × n square matrix. Further, a ij > 0 if and only if (j, i) ∈ E, which connotes the existence of a directed edge from node j to node i. Strong connectivity of G is the property that there exists a path connecting any two nodes, and is equivalent to A being irreducible [1] .

To begin, we review the classical Susceptible-Infected-Susceptible (SIS) network model [14] . We consider n > 1 populations of individuals, with each population being of large and constant size. Each individual has two possible health states, being healthy but susceptible (S), or infected (I) with the virus. Individuals who recover are assumed to do so with no immunity, being immediately susceptible again to infection. Defining the fraction x i ∈ [0, 1] of population i ∈ [n] who are infected, the SIS model posits that

where δ i > 0 is the recovery rate of population i and β ij ≥ 0 is the transmission rate of the virus from infected individuals in population j to susceptible individuals in population i. With x = [x 1 , x 2 , . . . , x n ] , X = diag(x), B = (β ij ), and D = diag(δ 1 , δ 2 , . . . , δ n ), the system equation becomes

The forward invariance of the set Ξ n for (2.2) is well known [14] , and this guarantees that x i , for all i ∈ [n], retains its physical meaning. The graph G = (V, E, B) captures the network structure over which the virus spreads, and the standard assumption that B is irreducible is equivalent to strong connectivity of G; the virus can reach any population i from any other population j through a path on G.

Given the irreducibility of B, a complete convergence characteristic can be determined by the reproduction number R ρ(D −1 B). If R ≤ 1, which is equivalent to s(−D + B) ≤ 0, the only equilibrium of (2.2) is the healthy equilibrium 0 n and it is globally asymptotically stable for x(0) ∈ Ξ n (exponentially stable if s(−D + B) < 0). If R > 1, i.e. s(−D + B) > 0, there is in addition to the healthy equilibrium (which is unstable) a unique nonzero/endemic 1 equilibriumx 0 n which is globally attractive for x(0) ∈ Ξ n \ 0 n [14, 18] . The equation (2.2) is the mean-field approximation of a stochastic discrete-time system with 2n states. Detailed discussion of the relationship between (2.2) and the stochastic model can be found in [21, 19, 33] , and we do not explore this aspect further.

This SIS framework can be extended to study the scenario in which there are circulating two viruses, termed virus 1 and virus 2, as in [17, 26, 28] . They are competing, in that an individual infected with one virus has immunity to an infection from the other virus, but like the single virus framework, an infected individual can recover with no immunity and immediately become susceptible again to infection from either virus. Similar to (2.2), the bivirus dynamics presented below is a mean-field approximation of a 3 n -state Markov chain model. Due to space limitations, we do not discuss the details of this approximation. The reader is referred to [17, Section II] and also the comprehensive Chapter 5 in the PhD thesis [27] for two different treatments.

Neighbor relationships for each virus are modeled by a directed graph between the populations (nodes), corresponding to the graph vertices, with an edge from node j to node i denoting the direction in which virus transmission can occur; the nonnegative infection rates β 1 ij , β 2 ij for i, j ∈ [n] capture the transmission rates for virus 1 and virus 2, respectively. Each group i ∈ [n] also has associated with it positive δ 1 i , δ 2 i corresponding to the rate of recovery from virus 1 and virus 2, respectively. With virus 1 and virus 2, we distinguish the vectors of fractions of infected individuals as x 1 (t) ∈ R n and x 2 (t) ∈ R n , respectively. The corresponding system equations becomeẋ

In order for x 1 , x 2 to have physical meaning, observe that in addition to the requirement that 0 n ≤ x i (t) ≤ 1 n for i = 1, 2, we require that x 1 (t) + x 2 (t) ≤ 1 n , the quantity on the left corresponding to the vector of fractions of individuals infected by either virus. Following a similar proof to the single virus case, an invariant set for the bivirus model can be identified which satisfies such requirements.

Lemma 2.1 ([17, Lemma 8] ). With the above notation, suppose that the initial conditions for (2.3) satisfy 0 n ≤ x i (0) ≤ 1 n for i = 1, 2, and x 1 (0) + x 2 (0) ≤ 1 n . Then for all t > 0, there holds 0 n ≤ x i (t) ≤ 1 n for i = 1, 2 and x 1 (t) + x 2 (t) ≤ 1 n .

In order to discuss existing results in the literature such as [17, 26, 37, 28] , and gain an appreciation of some open questions that remain unanswered, we introduce a key assumption to hold in the rest of the paper that parallels one typically imposed when analyzing the single virus case 2 .

Standing Assumption. The matrices B 1 , B 2 are nonnegative and irreducible. The matrices D 1 , D 2 are positive diagonal.

Positive diagonal D 1 , D 2 implies that for any population (node) i ∈ V, the recovery rate against each of the two viruses is strictly positive, but may differ between the viruses and between populations. By associating B i with the graph G i = (V, E i , B i ), for i = 1, 2, one can view G 1 and G 2 as the graphs capturing the network over which virus 1 and virus 2 spread, respectively. The assumption on irreducibility is equivalent to both G 1 and G 2 being strongly connected, but we allow the edges and associated infection rates to differ between the two graphs.

We first explain some reasoning behind several existing results of [17, 28, 37] , as opposed to reproducing the proofs. Consider action of a single virus, say virus 1, and the effect on x 1 in two cases: when virus 2 is not present, and when virus 2 is present. First, let us define the reproduction numbers of virus 1 and virus 2 as R 1 ρ((D 1 ) −1 B 1 ) and R 2 ρ((D 2 ) −1 B 2 ), respectively. Examination of (2.1) and (2.3a) shows that the presence of x 2 ≥ 0 n serves to slow down the increase in x 1 due to infection, because there necessarily holds (I n − X 1 (t) − X 2 (t))B 1 x 1 (t) ≤ (I n − X 1 (t))B 1 x 1 (t). This means that if in the absence of virus 2, the healthy state x 1 = 0 n is an attractive equilibrium state for virus 1, i.e. any nonzero x 1 (0) yields x 1 (t) → 0 n , then x 1 = 0 n must a fortiori remain an attractive equilibrium in the presence of virus 2. This implies then that if R 1 ≤ 1, any equilibrium (x 1 ,x 2 ) of the bivirus system (2.3) will necessarily be of the form (0 n ,x 2 ), and consequently in this casex 2 is necessarily an equilibrium of the single virus system applicable to virus 2 alone due to (2.3b). And then, there will further holdx 2 = 0 n if and only if R 2 ≤ 1.

(Of course, one can interchange the roles of the two viruses). These observations effectively cover [17, Theorems 1 to 3], [37] , and part of [28] . The remaining theorems in [17] focus on equilibria analysis of systems (2.3) with nongeneric parameter matrices D i , B i , i = 1, 2. As identified in [17] and as we will further explore in the sequel, stability and equilibria of nongeneric parameters are not always indicative of what happens for generic parameters.

The astute reader will recognize that the observations above also clarify that the condition R 1 > 1 and R 2 > 1 remains one of great interest, since this is not covered by a collapse to two single virus problems. Note that when both these conditions hold, there will still necessarily be three equilibria (0 n , 0 n ), (x 1 , 0 n ) and (0 n ,x 2 ); we term the first the healthy equilibrium and the latter two boundary equilibria of the bivirus system. We remark thatx 1 0 n andx 2 0 n are separately the unique endemic equilibria for each of the two single virus systems. This is evident since for x 2 (0) = 0 n , the underlying equations (2.3) imply x 2 (t) = 0 n for all t, and the evolution of x 1 (t) from some x 1 (0) is the same as what would occur with the corresponding single virus system in (2.1). This also implies that (x 1 , 0 n ) is the only equilibrium where virus 1 is present and virus 2 is not present. The same holds correspondingly with virus 1 and 2 interchanged.

The works [26, 28, 29] consider R 1 > 1 and R 2 > 1, and allow G 1 and G 2 to differ. In [26] , infection and recovery parameters are assumed homogeneous among the nodes for any one virus, but can be different between the two viruses. The existence and local stability of equilibria is studied using linearization and perturbation methods, coupled with extensive numerical simulations [26] . In [28] and [29] , differing sufficient conditions on infection parameters but with homogeneous recovery parameters are presented which guarantee the boundary equilibrium (0 n ,x 2 ) is attractive for all initial conditions satisfying x 2 (0) > 0 n (and thus a similar condition exists for (0 n ,x 2 ) to be attractive). However, a number of important open questions remain for generic bivirus systems. For instance, i) are there a finite number of equilibria? and ii) can a dichotomy be established for the typical limiting behaviour? In this paper, we answer these questions and more, extensively and rigorously expanding on existing results.

As background for the results in later sections, we outline the concept of monotone systems [31, 32] and recall some key results. Indeed, in the sequel, we use a standard method to demonstrate that the bivirus system is a monotone system, offering an alternative and simpler proof than the approach in [28] . The discussion in this subsection however is in terms of a general system

We assume henceforth that conditions on F exist which guarantee existence and uniqueness of solutions for all time, and indeed the Jacobian dF x exists at every point on allowed trajectories. Let φ t (x 0 ) denote the solution x(t) of (2.4) at time t when x(0) = x 0 ∈ U . There is special interest in the behavior of the difference between trajectory pairs when confined to particular orthants of R n , and we begin with the central definition.

Let m = (m 1 , m 2 , . . . , m n ), with m i ∈ {0, 1} for all i ∈ [n], be a prescribed sequence, and associate with it an orthant

We say x ≤ Km y and x Km y if y − x ∈ K m and y − x ∈ int(K m ) respectively. Here, int(·) denotes the interior of a set.

There is a straightforward condition for the monotone property, in terms of dF x , see [ 

Given this matrix characterization, it makes sense to define an irreducible monotone system as one for which P m dF x P m is irreducible for all x ∈ U .

One of the most important results concerning irreducible monotone systems is that only certain forms of limiting behavior are permitted, provided that there are known to be only a finite number of equilibria, see [ The above result establishes that the monotone system cannot exhibit chaos, and any limit cycles in M must be nonattractive. In other words, convergence to an equilibrium occurs for all initial conditions in M except possibly a set Y of measure zero. Note that if x i is unstable, then int(B(x i )) = ∅. This implies that convergence will occur to a stable equilibrium, except possibly for a set Z of initial conditions of measure zero which would yield convergence to an unstable equilibrium if one exists (a saddle point or a source). Clearly Z ⊂ Y, if in fact the two subsets are not empty.

3. Main results. In this section, we will present a series of new findings on the bivirus system. In order to keep the focus on the new results and to aid exposition and discussion, all proofs are presented in the Appendix. Moreover, and unless explicitly stated otherwise, the Standing Assumption is assumed to hold throughout.

In much of this section we will appeal to the equilibrium equations associated with (2.3). With an overbar denoting an equilibrium, these equations are

Since the outcomes when x i (0) = 0 n for some i are fully understood from the single virus dynamics, we focus on initial conditions in the set

Thus, when we refer to a globally attractive equilibrium, it is with respect to initial conditions in ∆, and an almost globally attractive equilibrium excludes initial conditions in a set of measure zero that is a subset of ∆.

3.1. General result on equilibria and convergence. In this subsection, we establish general properties of the equilibria of (2.3), and provide a general convergence result. First, (3.1) can be used to establish the following.

We term any equilibrium (x 1 ,x 2 ) withx 1 0 n andx 2 0 n a coexistence equilibrium. This lemma, whose proof is given in Appendix A.1, restricts substantially the equilibria which can lie on the boundary of the set

There are corresponding restrictions on the trajectories, as expressed in the following, which strengthens Lemma 2.1, and is essentially obtained as a result of taking into account the Standing Assumption (see Appendix A.2 for the proof). 3) has initial conditions satisfying (x 1 (0), x 2 (0)) ∈ ∆. Then for all finite 3 t > 0, there holds 0 n x i (t) 1 n for i = 1, 2 and x 1 (t) + x 2 (t) 1 n .

It is well known that ifx is an equilibrium of (2.4), then the (local) stability of x can often be determined through examination of the eigenvalues of the Jacobian dF x (x). If the linear systemż = dF x (x)z is (exponentially) asymptotically stable or is unstable, the same is true of (2.4), at least locally aroundx [30, Section 5.8] .

In the case of the bivirus system, the Jacobian provides information not only concerning equilibria. Denoting by J the Jacobian of the bivirus system, it is straightforward to verify that, withB i (

The Jacobian at a general point on a trajectory is the same as that in (3.2), save thatx i ,X i are replaced by x i , X i . Let P diag(I n , −I n ). Then observe (now with

It is immediately clear that this matrix is a Metzler matrix. However, more is true, as the following result illustrates (the proof is given in Appendix A.3):

Then the matrix P J(x 1 (t), x 2 (t))P evaluated at an arbitrary point along the trajectories of (2.3) with t < ∞, and the matrix P J(x 1 ,x 2 )P for an equilibrium satisfyingx 1 0 n ,x 2 0 n , are irreducible.

Thus the bivirus system given the Standing Assumption is an irreducible monotone system of type K m where m = (0 n , 1 n ). While monotonicity of the bivirus trajectories with homogeneous recovery rates was established in [28, Theorem 18] , the proof via Lemma 3.3 is significantly more direct. We also explicitly make a connection to the monotone systems literature 4 and competitive species models [11, 15, 3] . Indeed, the bivirus system is a cooperative system in the terminology of [31, 8] (not to be confused with competitive species models). A key contribution of this paper is to show how to leverage the monotone systems literature, including the seminal contributions by Morris M. Hirsh in the 1980s, to establish a number of new results for networked bivirus systems. Our results also extend the work of [3] , which examined a three-node bivirus network with a specific tree structure.

To begin, a general convergence result can be established using Lemma 2.3 if the bivirus system can be shown to have a finite number of equilibria. As we will shortly prove, given generic parameter matrices D i , B i , i = 1, 2, the bivirus system in (2.3) has a finite number of equilibria. We first define what we mean by "'free parameters" and "generic". Off-diagonal entries of D i are not free parameters, being always zero. If D i is constrained to be the identity matrix (as occurs later), the diagonal entries are also not free parameters. All entries of B i are free parameters. The free parameters of D i and B i take values from the nonnegative real interval, with those of D i required to be strictly positive while those of B i are only assumed to be nonnegative. At the very least, the existence of such a set would need to be demonstrated to conclude that generic values exist. Often, such a set can be characterized. For the bivirus system, the particular algebraic set is identified as part of the proof of the following result, and is presented in Appendix B.2.

Theorem 3.6. For generic parameter matrices D i , B i , i = 1, 2, the bivirus equation set (2.3) has a finite number of equilibria. If D i = I n , i = 1, 2, then for generic parameter matrices B i , i = 1, 2 the same conclusion holds. Moreover, for all initial conditions (x 1 (0), x 2 (0)) ∈ ∆, except possibly for a set of measure zero, the system (2.3) will converge to an equilibrium. If the system does not converge to an equilibrium, then it is on a nonattractive limit cycle.

This thoroughly answers the two questions posed at the end of subsection 2.2.

Concerning the first, and as detailed at the end of subsection 2.2, when R i ≤ 1 for some i = 1, 2, there are at most 3 equilibria; these 3 equilibria remain present when R i > 1 for i = 1, 2. Sahneh and Scoglio [26] showed then that for some parameter choices, there exist coexistence equilibria, in which both virus are present. However, [26] never proves (or even discusses) if the coexistence equilibrium is unique or indeed whether there are a finite number of them. Simulations suggest that one typically expects a finite number of equilibria, but a special scenario that yields a connected set containing an infinite number of equilibria was identified in [17] .

To keep the exposition clear, the precise algebraic set mentioned above is given in Appendix B.2, and we remark that for any n nodes, the algebraic set is always easily identified. The important point to notice is that the existence of such a set with measure zero (understood in the usual sense for real numbers) establishes the notion of "generic" parameter matrices D i and B i .

Concerning the second, previous works such as [17, 37, 28, 29, 3] establish sufficient conditions on the D i , B i parameter matrices to ensure a specific equilibrium among the healthy and two boundary equilibria is globally attractive for all initial conditions in ∆. In contrast, Theorem 3.6 establishes that the dynamical system itself is "almost globally stable", in the sense that the typical outcome is convergence to an equilibrium, not necessarily the healthy and two boundary ones as illustrated in subsection 3.4. This provides a general conclusion on the limiting behaviour of the dynamical system, as opposed to a general condition for a particular equilibrium to be globally attractive.

Our result does not specify which equilibrium the system converges to and this is deliberate. Identifying conditions on D i , B i that yield precise conclusions on the number of equilibria and associated regions of attraction is often difficult [12, 28, 29] . Indeed, there are instances of the bivirus dynamics with no global attractor equilibrium [3] . In the sequel, we provide an explicit example of such a network, whereby the two boundary equilibria each have a region of attraction with nonzero measure and the boundary between the two regions of attraction forms the stable manifold of an unstable coexistence equilibrium.

The remainder of this paper seeks to explore the properties of the equilibria, and we establish a number of different conclusions i) based on conditions for the D i and B i parameter matrices, and ii) exploiting the monotonicity of the bivirus trajectories.

3.2. Properties of the Equilibria. Because the Jacobian under transformation in (3.3) is a Metzler matrix, we can identify different collections of differential equations that give rise to the same equilibria that have the same local stability properties. We start with the following result. This proposition, with proof given in Appendix A.4, provides an important conclusion: for many theoretical investigations, there will be no loss of generality in assuming D 1 = D 2 = I n , which can simplify computations and reduce parametrization of the system to just the two matrices B 1 , B 2 defining the infection rates between nodes. To this end, note that the existence of (D i ) −1 ensures that (D i ) −1 B i is irreducible if and only if B i is irreducible. So the Standing Assumption 'flows through' if the simplification is undertaken. We further remark that the lemma makes no assertion that global stability (should it be present) in one system implies the same property for the other. Nonetheless, this property is indeed true, as we demonstrate with a simple extension of the lemma in the sequel.

Remark 3.8. This result allows us to examine scenarios in which the two viruses have different time-scales, by associating S andŜ in Lemma 3.7 with the quadruples B 1 , D 1 , B 2 , D 2 andB 1 = (D 1 ) −1 B 1 ,D 1 = I n ,B 2 = ( D 2 ) −1 B 2 ,D 2 = I n , respectively, for a positive that captures the relative rate of evolution of the two viruses. The location and the stability of equilibria remain invariant as varies in magnitude. While the invariance of the equilibria locations may seem intuitive, it is certainly not a trivial conclusion that the stability property of each equilibrium is preserved when the relative time-scale of a coupled nonlinear system is changed [30] . On the other hand, varying can certainly change the region of attraction of equilibria.

3.2.1. Nongeneric networks with an infinite number of equilibria. We will now show that with a generic B 1 , there are an infinite number of B 2 , satisfying a specific functional form, which yield a connected set containing an infinite number of equilibria, and this set comprises an interval of a straight line. We term this "a line of equilibria" for convenience.

For convenience, but with no loss of generality as demonstrated by Lemma 3.7, take D 1 = D 2 = I n and B 1 an arbitrary irreducible matrix. We require that s(−I n + B 1 ) > 0, thereby ruling out the possibility of the healthy equilibrium (x 1 ,x 2 ) = (0 n , 0 n ) being attracting for the bivirus system for any B 2 when (x 1 (0), x 2 (0)) ∈ ∆. Let 1 n z 0 n with Z = diag(z), satisfy

That is, z is the unique endemic equilibrium of the single virus system (2.2) with D = I n . Now let C be any other nonnegative irreducible matrix for which z is also an eigenvector with unity eigenvalue. Obviously it is straightforward to find such a matrix, and there is an infinite number. The Perron-Frobenius Theorem [10] implies ρ(C) = 1 and z is the only positive eigenvector of C (up to a scaling). Define Proposition 3.9. With D 1 = D 2 = I n , with B 1 an arbitrary nonnegative irreducible matrix and with z and B 2 defined using (3.4) and (3.5), a set of equilibrium points of the bivirus equilibrium equations (3.1) is given by (αz, (1 − α)z), with α ∈ [0, 1]. This equilibrium set is locally exponentially attractive, i.e., for all allowed initial conditions of the bivirus system sufficiently close to the set, the trajectory will approach the set exponentially fast.

The proof is given in Appendix A.5.

In [17] , a line of equilibria was shown to exist if D 1 = kD 2 , B 1 = kB 2 for some positive scalar k (and at one point even more specialized conditions). Here, we have relaxed this condition by identifying that a line of equilibria can exist under a broader type of nongenericity; the result in [17] is covered by Proposition 3.9 after applying Lemma 3.7. Whether it is possible to obtain an infinite number of equilibria not lying on a straight line is at the moment an open question.

It is often appealing to consider special examples of systems where the choice of parameters makes the derivation of analytic results possible. There is a risk however that nongenericity of the examples leads to nongeneric conclusions. The existence of nongeneric values of matrices defining a bivirus problem which simplify computations but give rise to an infinite number of equilibrium points is an illustration of the possibility.

In the situation described above, the connected set of equilibria forms an interval of a straight line, and each end of the interval corresponds to each of the single virus equilibria associated with the two viruses being active one at a time. In effect, a bifurcation produces this result. Suppose that B 2 above were replaced by µ(I n − Z) −1 C where µ is a scalar positive parameter. If µ < 1, one can verify that (x 1 ,x 2 ) = (z, 0 n ) is a locally exponentially stable equilibrium, while when µ > 1, one can check that it is a saddle point equilibrium, with divergent trajectories. When µ = 1, the Jacobian acquires a zero eigenvalue and the bifurcation occurs. A similar sort of analysis applies if B 2 is as above, but B 1 , after determination of z, is then replaced by µB 1 and one considers then the equilibrium at (x 1 ,x 2 ) = (0 n , z) with µ varying from less than to greater than one. This aspect will now be further explored.

We now explore the stability properties of the boundary equilibria (x 1 , 0 n ) and (0 n ,x 2 ) wherex 1 0 n andx 2 0 n are separately the unique endemic equilibria for each of the two single virus systems. Whilex 1 andx 2 are almost globally stable in the single virus systems, the stability of (x 1 , 0 n ) and (0 n ,x 2 ), local or global, in the bivirus case are not guaranteed. The above result, with proof appearing in Appendix A.6, provides necessary and sufficient conditions for local stability of the boundary equilibria, which can be determined by examining two separate single virus systems, with iterative algorithms available for computingx 1 andx 2 , e.g. [19, Theorem 5] and [18, Theorem 4.3] . Some insightful sufficient conditions can be obtained as follows (see Appendix A.7).

Corollary 3.11. With notation as above, the following statements hold:

1. If B 2 > B 1 , (x 1 , 0 n ) is unstable, and (0 n ,x 2 ) is locally stable, and there is no coexistence equilibrium (x 1 ,x 2 ) ∈ ∆. 2. If b 2 min i n j=1 β 2 ij >b 1 max k n j=1 β 1 kj , then (x 1 , 0 n ) is unstable, and (0 n ,x 2 ) is locally stable, and there is no coexistence equilibrium (x 1 ,x 2 ) ∈ ∆. 3. Ifx 2 >x 1 then (x 1 , 0 n ) is unstable, and (0 n ,x 2 ) is locally stable.

Item 1 and Item 2 of the corollary give different sufficient conditions on (x 1 , 0 n ) and (0 n ,x 2 ) to be unstable and locally stable, respectively. Since there are no coexistence equilibria (x 1 ,x 2 ), the only other equilibrium is the healthy equilibrium. Later, we present Corollary 3.16: when there is no coexistence equilibria, one of the two boundary equilibria is in fact stable for all initial conditions in ∆ (and in this instance, it would be (0 n ,x 2 )). Notice that Item 1 conditions are entry-wise on the B 1 and B 2 , while Item 2 concerns row sums. Neither subsumes the other, and it is possible to find examples of B i that satisfy one condition but not the other. Item 2 first appeared in [29, Corollary 4] (which refines [28, Theorem 20] ), and here we give an alternative and short proof. The sufficient condition in Item 3 of the corollary is implied by the condition in Item 1, but the converse is not necessarily true. That is,

x 1 , as detailed in [17, Proposition 4] . All three configurations of boundary equilibria stability properties can occur: i) both boundary equilibria are locally exponentially stable, ii) both boundary equilibria are unstable, and iii) one boundary equilibria is locally stable and the other unstable. In subsection 3.4 below, we give examples of each configuration for an n = 2 network.

Knowing the stability configuration of the two boundary equilibria, one can exploit the properties of monotone systems to obtain simple counting results that lower bound the number of interior equilibria, see [31, Theorem 2.8] . When there are no interior equilibria, a global stability result for one of the boundary equilibria follows. We explore this in further detail in the next subsection.

In this subsection, we argue that for a given system, all trajectories with x i (0) 0 n , i = 1, 2 are bounded by two special trajectories, starting from special corners of the allowed set of initial conditions. Other trajectory bounding results have appeared in [28, 29] , but by using two special trajectories, we are able to go further by demonstrating that from this result flows several important conclusions on equilibria properties. For instance, the results provide a simulation tool involving construction of two trajectories which on occasions can be expected to exclude the possibility of any other stable equilibria than those associated with the two trajectories, and on other occasions, to flag the presence of an unstable equilibrium. We draw on the conclusion, recorded after Lemma 3.3, that the bivirus system is type K m monotone with m = (0 n , 1 n ).

Theorem 3.12. Consider the equation set (2.3), and in particular consider the trajectories x A (t), x B (t) defined for arbitrarily small but positive η by the initial conditions x 1

. Suppose x C (t) is a trajectory beginning at any initial condition satisfying x 1

, and x 1 C (0) + x 2 C (0) < 1 n . Then the following inequalities hold for all t > 0

Appendix A.8 provides a very short proof based on known properties from the monotone systems literature, while [28] obtains a similar result using computations tailored to the bivirus dynamics.

The theorem also allows us to directly extend Lemma 3.7 to deal with bivirus systems with a globally stable equilibrium, as stated in the following result (with proof in Appendix A.9).

Lemma 3.13. Assume the same hypotheses as Lemma 3.7, with the two systems S andŜ. If one system has an equilibrium which is globally stable for all initial conditions in ∆, then the global stability property holds for the second system.

We further expand the theorem via a second corollary dealing with the limit points of the trajectories x A (t), x B (t), x C (t). Before stating it, we make an important observation, based on Lemma 2.3 and the discussions below it. Specifically, we know that if convergence to a stable equilibrium point does not occur for one of the initial conditions used in Theorem 3.12, then convergence will occur to a stable equilibrium for a perturbation of that initial condition within a ball of arbitrarily small radius. We abbreviate this idea below with the words 'perturbed if necessary'.

Corollary 3.14. With the same hypothesis as Theorem 3.12, with values for the matrices D i , B i assuring the number of equilibria is finite, and with initial conditions Fig. 1 : Illustration of Corollary 3.14. The triangles indicate initial conditions, with the system trajectory shown converging to the circles, which indicate the equilibria. The dotted rectangle shows the intersection of the hyperrectangle W with the plane defined by the ith coordinates of x 1 and x 2 .

perturbed if necessary, assume that the trajectories

with m = (0 n , 1 n ) as above, denote the closed hyperrectangle 5 whose edges are axis-parallel and which has a principal diagonal joining the pointsx A andx B . Then 1. If x C (0) lies outside W, the associated trajectory converges to a limitx C within W, obeying the constraints implied by Theorem 3.12. 2. If x C (0) lies within W, the associated trajectory either converges to a limit

x C again obeying the constraints implied by Theorem 3.12 or lies on a nonattractive limit cycle, which requires that W has dimension at least 2.

The above corollary is illustrated in Figure 1 , and proved in Appendix A.10. The corollary effectively offers a simulation-type test for establishing whether there is a single equilibrium for whichx i 0 n , i = 1, 2: one simply computes two trajectories, x A (t) and x B (t), from initial conditions as close as possible to opposite corners of the region of interest, and checks that they approach a common limit. What happens if the two trajectories do not approach a common limit? The answer, again appealing to monotone systems theory, is that an unstable equilibrium must lie in W. The following is an amalgam of Theorem 2.8 and Proposition 2.9 in [31] , and thus no proof is presented. x D x B , then there are two unstable trajectories emanating from the equilibriumx C that are tangent to the eigenvector associated 5 In some cases, one or more coordinates ofx A andx B may assume the same value, (though if

is possible by Lemma 3.1), and W is then degenerate.

with the unstable eigenvalue of the associated Jacobian matrix J(x C ), satisfying the following properties. The first trajectory satisfiesẋ 1 0 n and −ẋ 2 0 n along its entirety and its limiting point isx B . The second trajectory satisfies −ẋ 1 0 n anḋ x 2 0 n along its entirety and its limiting point isx A .

Finally, if the particular system has no coexistence equilibria, then, as recorded in the next corollary the properties established above imply that one of the boundary equilibria is globally stable, being attractive for all initial conditions in ∆. Sufficient conditions on D i , B i for there to be no coexistence equilibria include Corollary 3.11 and [12, Theorem 6] . The proof is in Appendix A.11. Corollary 3.16. For the bivirus system in (2.3) with generic parameter matrices, suppose there are no coexistence equilibria. Then precisely one of the boundary equilibria is an attractive equilibrium, and the set ∆ is in its region of attraction.

3.4. Equilibria for a low order system. We first explain how, for generic values of system parameters and with n = 2, one can compute using nothing more than a single quadratic equation solution followed by linear equation solutions, any equilibrium for whichx 1 0 n ,x 2 0 n . Then we illustrate a series of possible outcomes which provides a comprehensive account of the possible limiting behaviours for the bivirus system.

The solution procedure is straightforward. Without loss of generality, assume D 1 = D 2 = I 2 . Setting α = (x 1 2 /x 1 1 ) and γ = (x 2 2 /x 2 1 ), the equilibrium equations in (3.1) yields after some manipulation

One can eliminate γ and obtain a quadratic equation for α. For each solution of the quadratic equation, it is then possible using only linear equations to obtain values for thex i j . Note that these have no constraints on the signs of their entries. The fact that a quadratic equation underpins the algortihm means that there can be at most two interior equilibria. Exploiting properties of monotone systems, an argument centred around [31, Proposition 2.9] will lead to the conclusion that if there are two interior equilibria, one must be stable and the other unstable. In the examples below, we report either no interior equilibria, or the interior equilibrium is unique.

With D 1 = D 2 = I 2 , we fix B 1 = 1.6 1 1 1.6

and vary the B 2 matrix according to Table 1 . The equilibria (x 1 , x 2 ), with

, for each case as computed using the above solution procedure are reported, excluding the healthy equilibrium (0 2 , 0 2 ).

Case 1: A line of coexistence equilibria (which is locally exponentially attractive by Proposition 3.9) joins the boundary equilibrium (0.615 · 1 2 , 0 2 ) to the other boundary equilibrium (0 2 , 0.615 · 1 2 ). There are no other coexistence equilibria other than those in the line, and extensive simulations suggest the line of equilibria is in fact globally attractive for ∆. results on regions of attraction for stable equilibria [4] , it follows that the regions of attraction for (x 1 , 0 2 ) and (0 2 ,x 2 ) encompass all of ∆ except for a set of measure zero, being the attractive manifold of (x 1 ,x 2 ) on which there are no two distinct points y, z with y < Km z. In fact, this manifold is part of the boundary of the regions of attraction for (x 1 , 0 2 ) and (0 2 ,x 2 ). Remark 3.17. Case 2 has a particularly notable outcome, because in the context of a "survival-of-the-fittest" battle between the two viruses, different initial conditions in ∆ can result in either virus surviving. A similar outcome was presented in [3] , but for a three-node network with specific tree structure. Our recent work has proposed a method to systematically construct bivirus networks, with an arbitrary number of nodes and arbitrary topology, where either virus can win the battle [38] . Such outcomes are distinct from existing results such as [16, 26, 28, 29] , which identify parameter regimes for a specific virus to win the battle independent of initial conditions.

In this paper, we have analyzed a deterministic networked bivirus model. We have proved the bivirus system with generic parameters has a finite number of equilibria. Using monotone systems theory, and assuming generic parameters, we established convergence to an equilibrium for (almost) all initial conditions. The properties of equilibria have been further explored, using a mixture of monotone systems theory, matrix theory, and algebraic analysis. Future work will focus on the implications of stable/unstable boundary equilibria on the number and stability of interior equilibria, and control strategies that leverage one virus to eliminate the other, perhaps by optimal design of parameters to ensure a specific boundary equilibrium is globally stable.

Appendix A. Proofs of main results.

A.1. Proof of Lemma 3.1. As we know, there are equilibria (0 n , 0 n ), (x 1 , 0 n ) and (0 n ,x 2 ) for vectorsx 1 ,x 2 which are separately equilibria of single virus systems. Becausex 1 0 n is the unique endemic equilibrium associated with the dynamics (2.2) if only virus 1 is considered, it is clear that there cannot exist an equilibrium (x 1 , 0 n ) with 0 n ≤x 1 =x 1 . By the same reasoning, there cannot exist an equilibrium (0 n ,x 2 ) with 0 n ≤x 2 =x 2 .

If and only if R i > 1, or s(−D i + B i ) > 0 for i = 1, 2, it becomes possible for further equilibria to exist withx i > 0 n . To complete the first part of the proof, we show that there can be no equilibrium (x 1 ,x 2 ) in whichx i = 0 n , butx i 0 n fails. Observe first that there exists no j for which 1−x 1 j −x 2 j = 0. For if there were such a j, the j-th row of the two equilibrium equations (3.1) would imply x 1 j = x 2 j = 0, a contradiction. Hencex 1 +x 2 1 n and I n −X 1 −X 2 is a nonsingular positive diagonal matrix, so that (I n −X 1 −X 2 )B i and indeed (D i ) −1 (I n −X 1 −X 2 )B i for i = 1, 2 are then both irreducible and nonnegative. This means that for any y > 0 n , there is at least one positive entry of the vector (D i ) −1 (I n −X 1 −X 2 )B i y in those rows where y has a zero entry (see subsection 2.1). Since howeverx i = (D i ) −1 (I n −X 1 −X 2 )B ixi , identifying y withx i would yield a contradiction, unlessx i 0 n . Lastly, observe that

and becausex 1 0 n and B 1x1 0 n it follows thatX 2 =X 2 orx 2 =x 2 as required.

It is straightforward to derive thaṫ

Consider the second equation. If z i (t) = 0 for some i, 1 − x 1 i (t) − x 2 i (t) = 0, so that one at least of x 1 i (t) and x 2 i (t) is nonzero, i.e.ż i is positive. In fact, there exists , such that if z i (t) ∈ [0, ], thenż i (t) > 0; this implies that z i (t) can never be zero for t > 0, i.e. z(t) 0 n ∀t > 0. Obviously this also establishes that x i (t) 1 n for t > 0. Next, suppose that at some time T , k entries of x 1 (t) in positions, i 1 , i 2 , . . . i k say, are zero, with k < n. Then since B 1 is irreducible, B 1 x 1 must have at least one entry in one of the positions i 1 , i 2 , . . . i k which is nonzero, as detailed in subsection 2.1. Suppose it is position i 1 . Then from the above equations,ẋ 1 i1 (t) > 0. It follows that for t > T with t − T sufficiently small, there are fewer than k entries of x i which are zero. An extension of this argument in fact shows that for all t > 0, all entries of x i must be nonzero, i.e. x i (t) 0 n for t > 0.

A.3. Proof of Lemma 3.3. Under the hypothesis of the lemma, Lemma 3.2 guarantees that I n − X 1 − X 2 is nonsingular, and so (I n − X 1 − X 2 )B i is irreducible. Moreover,B i , i = 1, 2 are nonsingular positive diagonal matrices since x i 0 n for i = 1, 2. Evidently, it is enough then to show that the nonnegative matrix

is irreducible, which is equivalent to the graph G = (V, E, H) being strongly connected. The graph G is comprised of two strongly connected subgraphs, call them G 1 and G 2 , associated with irreducible matrices B 1 and B 2 , respectively. The matrixB 1 captures edges from nodes in G 2 to nodes in G 1 , whileB 2 captures edges from nodes in G 1 to nodes in G 2 . SinceB i are positive diagonal for i = 1, 2, there exists a path from any node k in G 1 to any node j in G 2 , and vice versa. It follows that G is strongly connected, and hence H is irreducible. The proof for the equilibrium value P J(x 1 ,x 2 )P is effectively identical.

A.4. Proof of Lemma 3.7. When the two equations of (3.1) corresponding to S are multiplied on the left by (D 1 ) −1 and (D 2 ) −1 , respectively, equilibrium equations forŜ result (and the converse is trivial). The proves that the equilibrium sets are identical. The equivalence of equilibrium stability properties requires slightly more work. Begin with system S. BecauseJ S =: P J S (x 1 ,x 2 )P is Metzler, local exponential stability of the equilibrium (x 1 ,x 2 ) implies that −J S is an M -matrix. As the Mmatrix property is preserved under multiplication by a positive diagonal matrix, it follows on multiplyingJ S by diag (D 1 ) −1 , (D 2 ) −1 that the result is another Mmatrix. But this multiplication yieldsJŜ = −P JŜ P . Thus the stability property for the equilibrium is the same. An identical argument works if −J S is a singular M -matrix and by negation if it is not an M -matrix. . We now consider the Jacobian at such an equilibrium. LetB i α denote diag(B ixi α ). The Jacobian at any equilibrium on this interval is

Now the equilibrium equations give immediately diag(B 1x1 α ) = α(I n − Z) −1 Z and diag(B 2x2 α ) = (1 − α)(I n − Z) −1 Z. Then, J is seen to be similar tō

using the transformation matrix P = diag(I n , −I n ). The matrixJ is an irreducible Metzler matrix. Further, [z , z ] 0 n is a nullvector ofJ for all α. It follows that all other eigenvectors of the matrix have negative real parts. It also follows that the bivirus equations defining the line of equilibria define a one-dimensional center manifold along which the Jacobian is singular. By standard center manifold theory, see e.g. [30, Section 7.6] , the eigenvalue properties of J then imply the center manifold is exponentially attractive, i.e. for initial conditions sufficiently close to the manifold, convergence occurs to some point on the manifold exponentially fast.

A.6. Proof of Theorem 3.10. We establish the result for the boundary equilibrium (x 1 , 0 n ), with the proof for (0 n ,x 2 ) being identical after adjustment of certain indices. First, we give some relevant results concerning the single virus model (2.2). Let R > 1, and thus the unique endemic equilibriumx of (2.2) satisfies

The positive vectorx 0 n can be seen as the nullvector of the matrix P = −D + (I n −X)B, which is an irreducible Metzler matrix. From the discussions in section 2, it follows that s(P ) = 0, which in turn implies that −P is a singular irreducible Mmatrix. It is known that an irreducible singular M -matrix plus a nonnegative diagonal matrix with at least one positive diagonal element yields an irreducible nonsingular M -matrix [25, Theorem 4.31] . WithB = diag(Bx), it follows that −P +B = D − (I n −X)B +B is a nonsingular M -matrix, and thus s(−P +B) < 0.

The Jacobian J(x 1 , 0 n ) is upper block triangular, as seen from (3.2). From the arguments immediately above, one can conclude that the upper diagonal block matrix −I n + (I n −X 1 )B 1 −B 1 is the negative of a nonsingular M -matrix, and is therefore Hurwitz. The lower diagonal block matrix I n + (I n −X 1 )B 2 is an irreducible Metzler matrix, and note thatX 1 is uniquely determined by B 1 and is therefore independent of B 2 . As outlined in section 2, one has that s − I n + (I n − Z 1 )B 2 < 0 ⇔ ρ (I n − Z 1 )B 2 < 1. The condition for instability of (x 1 , 0 n ) can be similarly proved. This completes the proof.

A.7. Proof of Corollary 3.11. Item 1: As detailed below (A.4), one has that s(−I n +(I n −X 1 )B 1 ) = 0, which in turn implies that ρ((I n −X 1 )B 1 ) = 1, according to section 2. Since B 2 > B 1 , it follows that there exists at least one entry of (I n −X 1 )B 2 is strictly greater than the corresponding entry of (I n −X 1 )B 1 . Then, [34, Theorem 2.7] establishes that ρ((I n −X 1 )B 2 ) > ρ((I n −X 1 )B 1 ) = 1, which in conjunction with Theorem 3.10 delivers the claim on instability. The argument for local stability of (0 n ,x 2 ) is similar, starting with the observation s(−I n + (I n −X 2 )B 2 ) = 0.

To show there is no equilibrium of the form (x 1 ,x 2 ) withx 1 0 n ,x 2 0 n , let us assume to the contrary that such a (x 1 ,x 2 ) exists. Then, (3.1a) indicates thatx 1 is a positive eigenvector associated with the simple eigenvalue 0 of the irreducible Metzler matrix −I + (I −X 1 −X 2 )B 1 . Let y be the associated positive left eigenvector, normalized to satisfy y x 1 = 1. Set C = (I −X 1 −X 2 )(B 2 − B 1 ). Observe that C > 0 n×n . Further, (3.1b) yields [−I+(I−X 1 −X 2 )B 1 ]x 2 +Cx 2 = 0 n . Premultiplying the left by y gives a contradiction, since C > 0 n×n , y 0 n ,x 2 0 n . Item 2: We first show there are no coexistence equilibria, and then show the instability of (x 1 , 0 n ). Suppose, to obtain a contradiction, that there is an interior equilibrium (x 1 ,x 2 ). Let J, K be indices such thatx 1 i ≤x 1 ary equilibrium. Suppose that there exists an interior point x(0) = (x 1 (0), x 2 (0)) for which the associated trajectory does not converge to the stable boundary equilibrium (perhaps on a nonattractive limit cycle). Let N be a suitably small neighbourhood of x(0). Then, there exist x A (0) ∈ N and x B (0) ∈ N with x B (0) x(0) x A (0) such that the trajectories starting from x A (0) and x B (0) converge to the stable boundary equilibrium, and Corollary 3.14 implies that the trajectory beginning at x(0) also converges to the common limit, delivering the contradiction.

Appendix B. Applying algebraic geometry to the proof of Theorem 3.6. The main focus is to prove that for generic D i , B i , the bivirus system in (2.3) has a finite number of equilibria. To do this, we shall first argue that there can only be an infinite number of equilibria for values of the free parameters in D i , B i lying in an algebraic set, i.e. a set defined by setting a multivariate polynomial in the parameters to zero. The existence of this set will be demonstrated by algebraic geometry. If the multivariate polynomial is the zero polynomial, the algebraic set comprises the whole space. If it is not zero, then the algebraic set has measure zero and for almost all parameter values, there will be a finite number of equilibria. Then we prove that for any given n, there exists a specific choice of D i and B i for which the bivirus system has a finite number of equilibria, implying that the algebraic set cannot be the whole set and so necessarily has measure zero.

B.1. Background on algebraic geometry. We first provide relevant background on algebraic geometry, the main reference being [5, Chapter 3] .

With x denoting an n-vector of unknowns, a monomial is a term of the form x α1 1 x α2 2 . . . x αn n , where the α i are nonnegative integers. A single multivariate polynomial equation is obtained by taking a linear combination of such monomials and setting it to zero. Algebraic geometry enables examination of the solvability question for n such equations in n unknowns. They can be written in the form P (x, β) = 0 where β is the set of free coefficients of the different monomials. (In any one equation, some monomials may be absent, and others again may have a fixed coefficient, such as 1; the coefficients of the remaining monomials make up the entries of β). We focus on the case of real polynomial equations, so that β can be regarded as a real vector, the first block of entries corresponding to monomials in the first equation, the second block to monomials in the second equation, and so on.

In order to deal with a phenomenon (so-called 'solutions at infinity') not encountered in working with a purely linear equation set, a technical device is introduced.

A single multivariate polyonomial equation is termed homogeneous when the sum of the powers in each monomial in the equation is the same for all terms. If one has a nonhomogeneous equation, it can be made homogeneous by introducing a further variable, x 0 , with a monomial x α1 1 x α2 2 . . . x αn n being replaced by x α0 0 x α1 1 x α2 2 . . . x αn n with α 0 chosen to ensure all monomials in the single polynominal equation have the same total degree n i=0 α i . Consider in fact n such homogeneous polynomial equations in the n + 1 unknowns x 0 , x 1 , . . . x n . Note that while separately homogeneous, the individual equation degrees are not necessarily the same. The equations may well have been obtained by making homogeneous the inhomogeneous set P (x, β) = 0. We shall write the homogeneous set asP (x 0 , x, β) = 0. This is related to P byP (1, x, β) = P (x, β). The key result from algebraic geometry is the following, [5, pg. 86, Theorem 2.3].

Theorem B.1. Consider a set of n homogeneous polynomial equations, denoted P (x 0 , x, β) = 0 with free parameters in the vector β and with unknowns x 0 , x 1 , . . . , x n .

Then there exists an expression R(β) (termed the resultant) that is polynomial in the entries of β such that if R(β) = 0 for particular values of the β entries, there are a finite number of nonzero solutions (possibly complex) toP (x 0 , x, β) = 0, disregarding scaling. If R(β) = 0, there are either no solutions or an infinity of solutions (disregarding scaling).

Evidently, if (x 0 ,x 1 , . . . ,x n ) is a solution to a homogeneous set of equations P (x 0 , x, β) = 0, then (λx 0 , λx 1 , . . . λx n ) for any nonzero real or complex λ is also a solution. Ifx 0 = 0, the choice λ = (x 0 ) −1 recovers a solution to the inhomogeneous set P (x, β) = 0. However, ifx 0 = 0, a solution to P (x, β) = 0 cannot be recovered in this way; such solutions are however sometimes termed 'solutions at infinity' of the set P (x, β) = 0, [5, pg. 115 ].

Corollary B.2. Assume the same hypotheses as Theorem B.1. If for a particular choice of values for the entries of β, call itβ, the resultant takes a nonzero value, then for almost all choices of β, the resultant will be nonzero.

The proof is immediate: it is obvious that if a polynomial in a single variable is nonzero for some value of that variable, it is nonzero for almost all values, i.e. nonzero everywhere except on a set of measure zero (which in the case of a polynomial in a single variable is a finite set). The first property is clearly also true for any (nonzero) multivariate polynomial such as the resultant.

Thus if the resultant takes a nonzero value for a particular choice of valuesβ for the free parameters of the equation setP (x 0 , x, β), then for almost all values assumed by the entries of β, the resultant polynomial R(β) will evaluate as a nonzero number. Further, if it takes a nonzero value for a particular choice of values for the free parameters in the equation setP (x 0 , x, β) = 0 and there is an associated solution with x 0 = 0, i.e. there is a solution of P (x, β) = 0, the same will hold true for almost all values of the free parameters.

B.2. Proof of Theorem 3.6. As noted just prior to Definition 3.4, evidently the key part of the theorem we have to prove is that the number of equilibria is finite.

The equilibrium equations for the bivirus system, given in (3.1), are a set of polynomial equations in thex j i . Suppose that we take D 1 = D 2 = I n , and B 1 , B 2 as positive diagonal matrices with no two corresponding entries equal. The equation set then becomes decoupled, and the i-th entry of each of the vector equations in (3.1) is

ii ]x 2 i = 0 The solvability of these equations yielding a finite number of solution pairs (x 1 i ,x 2 i ) is easily checked. The solutions in fact are (0, 0), (0, 1 − (1/β 2 ii )), (1 − (1/β 1 ii ), 0). Of course, the other entries of thex i can be treated in the same way due to the decoupling. Since there is a particular choice for the entries of B 1 , B 2 and D 1 , D 2 (the latter particular choice being the identity matrix) yielding a finite number of solutions, the algebraic geometry arguments presented above show that for almost all choices of the free parameters in B i , and the diagonal entries of D i , a finite number of solutions to the equilibrium equations exist.

Due to Lemma 3.2, it follows that for any x(0) ∈ ∆, we have x(τ ) ∈∆ for some finite time τ ≥ 0, with∆ an open strict subset of ∆ that is positively invariant. With finiteness of the equilibria assured, we apply Lemma 2.3 to the bivirus system (2.3), relating M to∆. Note, we can relate M to ∆ {x 1 , x 2 |0 n ≤ x i ≤ 1 n for i = 1, 2, and x 1 + x 2 ≤ 1 n }.

Nonnegative Matrices in the Mathematical Sciences, Computer Science and Applied Mathematics

Competitive exclusion in gonorrhea models and other sexually transmitted diseases

Competitive exclusion and coexistence of multiple strains in an SIS STD model

Stability Regions of Nonlinear Dynamical Systems: Theory, Estimation, and Applications

Using algebraic geometry

Dynamical Interplay between Awareness and Epidemic Spreading in Multiplex Networks

Systems of differential equations which are competitive or cooperative: I. limit sets

Systems of differential equations that are competitive or cooperative ii: Convergence almost everywhere

Topics in matrix analysis

Matrix Analysis

Competitive Exclusion and Coexistence for Competitive Systems on Ordered Banach Spaces

Networked multivirus spread with a shared resource: Analysis and mitigation strategies

Competing epidemics on complex networks

A Deterministic Model for Gonorrhea in a Nonhomogeneous Population

A remark on the global dynamics of competitive systems on ordered banach spaces

On the Analysis of a Continuous-Time Bi-Virus Model

Analysis and Control of a Continuous-Time Bi-Virus Model

On the dynamics of deterministic epidemic propagation over networks

Virus spread in networks

Interacting Epidemics and Coinfection on Contact Networks

Analysis and Control of Epidemics: A Survey of Spreading Processes on Complex Networks

Epidemic processes over time-varying networks

Epidemic processes in complex networks

Winner takes all: Competing viruses or ideas on fair-play networks

Cooperative Control of Dynamical Systems: Applications to Autonomous Vehicles

Competitive epidemic spreading over arbitrary multilayer networks

Bi-virus Epidemics over Large-scale Networks: Emergent Dynamics and Qualitative Analysis

Bi-virus sis epidemics over networks: Qualitative analysis

Sufficient Condition for Survival of the Fittest in a Bi-virus Epidemics

Nonlinear systems: analysis, stability, and control

Systems of Ordinary Differential Equations Which Generate an Order Preserving Flow. A Survey of Results

Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems

Accuracy criterion for the mean-field approximation in susceptible-infected-susceptible epidemics on networks

Matrix Iterative Analysis

Modelling COVID-19

Competing memes propagation on networks: A network science perspective

A Bi-Virus Competing Spreading Model with Generic Infection Rates

Competitive epidemic networks with multiple survival-of-the-fittest outcomes

Applications of the Poincaré-Hopf Theorem: Epidemic Models and Lotka-Volterra Systems

Acknowledgments. We would like to acknowledge the anonymous referees and their comments and suggestions for improving the manuscript.

For any irreducible nonnegative matrices A 1 , A 2 the condition A 1 > A 2 implies ρ(A 1 ) > ρ(A 2 ) and ρ(A 1 ) ≥ min i n j=1 a ij . This means that ρ[(I −X 1 )B 2 ] ≥ ρ[(1 − x 1 J )B 2 ] ≥ (1−x 1 J )b 2 > 1 and the instability claim is complete. The proof that (0 n ,x 2 ) is locally stable, by establishing ρ[(I −X 2 )B 1 ] < 1, is the same, mutatis mutandis. Item 3: Sincex 2 >x 1 , it immediately follows that (I n −X 1 )B 2 > (I n −X 2 )B 2 , which according to [34, Theorem 2.7] implies that ρ((I n −X 1 )B 2 ) > ρ((I n −X 2 )B 2 ) = 1, thus establishing the instability of (x 1 , 0 n ) by the result of Theorem 3.10. Similarly, we have that (I n −X 1 )B 1 > (I n −X 2 )B 1 , and it follows from [34, Theorem 2.7] that ρ((I n −X 2 )B 1 ) < ρ((I n −X 1 )B 1 ) = 1, and hence (0 n ,x 2 ) is locally stable.A.8. Proof of Theorem 3.12. It is evident that with m = (0 n , 1 n ) the initial conditions specified in the theorem statement obey 0)) holds for all t. These inequalities yield the first two lines of (3.6). The inequality of the last line was established in Lemma 3.2.A.9. Proof of Lemma 3.13. Without loss of generality, suppose that S has a globally stable equilibrium atx. To obtain a contradiction, suppose there is an initial state x(0) ∈ ∆ forŜ for which the trajectory x(t) does not converge tox. Let N ⊂ ∆ be a suitably small neighborhood of x(0). Then there exist x A (0) ∈ N and x B (0) ∈ N with x B (0)x(0) x A (0) such that, by Theorem 3.6, the trajectories x A (t) and x B (t) starting from x A (0) and x B (0) converge to a stable equilibrium ofŜ. But by Lemma 3.7, the systemŜ cannot have a stable equilibrium other thanx, since the system S has no stable equilibrium point other thanx. In other words, convergence of trajectories in the systemŜ must occur tox. Now use the key result of Theorem 3.12 onŜ. The inequality for the initial conditions implies x B (t)x(t) x A (t) ∀t and taking the limit as t → ∞ yields lim t→∞ x(t) =x. Note we can always select x A (0) and x B (0) to satisfy Theorem 3.12; due to Lemma 3.2, any x(0) ∈ ∆ will satisfy 0 n x i (t) 1 n , i = 1, 2, after a finite time t.A.10. Proof of Corollary 3.14. The underlying system is an irreducible monotone system and by assumption, almost all trajectories converge to a limit point (which may be locally stable, or a saddle point); the remaining trajectories are nonattractive limit cycles. The inequalities established in Theorem 3.12 imply that the limiting trajectory resulting from x C (0) necessarily lies within the hyperrectangle W, and this can only be a nonattractive limit cycle if x C (0) itself lies in W and W is neither a single point nor a one-dimensional interval. These observations imply one of the three conclusions in the corollary statement must hold.A.11. Proof of Corollary 3.16. The argument relies on the fact that all trajectories of the bivirus system approach a locally stable equilibrium, except possibly from a set of initial conditions of measure zero, as per Theorem 3.6 and Lemma 2.3. If there were two unstable boundary equilibria, a contradiction is immediate because almost all trajectories approach a stable equilibrium; boundary equilibria are excluded through their instability and a coexistence equilibrium is not present by hypothesis. If there were two stable boundary equilibria, Corollary 3.14 yields existence of an interior unstable equilibrium, another contradiction. Hence there is precisely one stable bound-