<figure><img src="image_1.png" /></figure>

HAL
open science

Networked Competitive Bivirus SIS spread with Higher
Order Interactions

Sebin Gracy, Brian D.O. Anderson, Mengbin Ye, César A. Uribe

► To cite this version:

Sebin Gracy, Brian D.O. Anderson, Mengbin Ye, César A. Uribe. Networked Competitive Bivirus SIS spread with Higher Order Interactions. 2025. hal-04669193v3

HAL Id: hal-04669193

https://hal.science/hal-04669193v3

Preprint submitted on 12 Sep 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
NETWORKED COMPETITIVE BIVIRUS SIS SPREAD WITH
HIGHER ORDER INTERACTIONS*

SEBIN GRACY†, BRIAN D.O. ANDERSON ‡, MENGBIN YE§, AND CÉSAR A URIBE¶

**Abstract.** The paper studies the simultaneous spread of two competing viruses over a network of population nodes by also accounting for the possibility of higher-order interactions (HOI). To this end, we consider a continuous-time time-invariant competitive bivirus networked susceptible-infected-susceptible (SIS) HOI system. First, we show that, under the assumption that the hypergraph associated with the system is strongly connected, the model is strongly monotone. Subsequently, using the Parametric Transversality Theorem of differential topology, we show that, for generic parameter choices, the system admits only a finite number of equilibria, and that the Jacobian, evaluated at any equilibrium, is a nonsingular matrix. The aforementioned two findings together guarantee that the typical behavior of the model is convergence to some stable equilibrium point. The equilibria of this system are i) the disease-free equilibrium (DFE), ii) single-virus endemic equilibria, and iii) coexistence equilibria (where both viruses are present in separate fractions of a population node). We identify a parameter regime that admits the possibility of three equilibria (namely, the DFE, and two single-virus endemic equilibria) being simultaneously stable. We then provide sufficient conditions for the existence of a coexistence equilibrium, both for the same parameter regime as mentioned above, and for a different one. Thereafter, we identify a necessary condition for the existence of a coexistence equilibrium. Finally, we illustrate our results using several numerical examples.

**Key words.** Epidemic processes, competing viruses, finiteness of equilibria, coexistence equilibrium.

**MSC codes.** 34D05, 37C75, 92D30

**1. Introduction.** Epidemics have massively influenced the trajectories of civilizations and societies. Disease outbreaks may accelerate certain changes in society, such as the large-scale adoption of digital technologies in daily lives, as witnessed during the COVID-19 pandemic [25]; permanent shutting down of many small businesses [26]; exacerbation of social unrest, and the possible weakening and eventual collapse of existing governments [35]. Given the consequential effects of disease outbreaks, it has increasingly drawn the attention of several scientific communities besides the medical community, such as physics [42], mathematics [6], computer science [36], and economics [5]. The basic quest behind all these efforts is to understand when a disease spreads in a community and leverage this understanding for the effective design of disease mitigation/eradication strategies.

Fundamental to such a pursuit has been the development of compartmental models; some of the most popular models include susceptible-infected-recovered (SIR) [27];

<footnote>*Submitted to the editors on August 1, 2024. A preliminary version of this paper has been accepted for publication in the Proceedings of the 2024 American Control Conference, [14].</footnote>

<footnote>**Funding:** Part of this material is based upon work supported by the National Science Foundation under Grants #2211815 and No. #2213568. The work of MY is funded by the Western Australian Government through the Premier's Science Fellowship Program, and by the Australian Government through the Office of National Intelligence (NI240100203).</footnote>

<footnote>†Department of Electrical Engineering and Computer Science, South Dakota School of Mines and Technology, Rapid City, SD, USA. (sebin.gracy@sdsmt.edu)</footnote>

<footnote>‡School of Engineering, Australian National University, Canberra, Australia. (brian.anderson@anu.edu.au)</footnote>

<footnote>§Centre for Optimisation and Decision Science, Curtin University, Australia. (mengbin.ye@curtin.edu.au).</footnote>

<footnote>¶Department of Electrical and Computer Engineering, Rice University, Houston, TX, USA. (cauribe@rice.edu)</footnote>
susceptible-exposed-infected-recovered (SEIR) [31]; susceptible-infected (SI) [23]; and susceptible-infected-susceptible (SIS) [29, 42]. The reasons why mathematical models have gained traction across multiple scientific communities are threefold: First, they provide strict guarantees regarding (some of) the various outcomes associated with disease outbreaks (in the model framework, of course). Second, the conditions involved in providing said guarantees can be checked. Lastly, because of the low cost of simulations, they can generate insightful observations relevant to a real rather than model epidemic. The present paper deals with a networked SIS model, with each node in the network being a proxy for a large population, and the interconnection among the nodes indicating the pathways for the virus to spread. The interest in networked SIS models is not new; these have been studied since the classical paper [29].

Classic epidemic models pertain to the scenario where there is a *single* virus circulating in a metapopulation (i.e., a network of population nodes). It is, however, not unusual for settings where multiple viruses are *simultaneously* circulating in a metapopulation. In such settings, the viruses to be competitive [32] or cooperative [15]. In the competitive case, infection with one virus rules out the possibility of *simultaneous* infection with another virus. Examples of such competitive viruses include simultaneous circulation of a) multiple strains of influenza A virus [34]; b) gonorrhea and a strain of partially drug-resistant gonorrhea [7]; and c) influenza A virus and human respiratory syncytial virus [8]. Other than disease outbreaks, the notion of competitive viruses also finds resonance in, among other scenarios, product adoption in a marketplace [2], and the spread of opinions in social networks [33].

The dynamics of competitive networked multi-virus SIS models are more involved than their single-virus counterparts; see, for instance, [32, 45, 13]. Specifically, the complexity stems from the fact that competitive networked multi-virus SIS models, in sharp contrast to single-virus SIS models, exhibit phenomena such as coexistence (i.e., multiple viruses infecting separate fractions of a population node) and competitive exclusion (i.e., one (or more) virus(es) pushing out the remaining virus(es)). Motivated by these challenges, these models have received a lot of attention lately; see [44, 16, 1]. Nonetheless, classical multi-virus SIS models suffer from the following limitation: they do not account for the possibility that an individual might *simultaneously* interact with more than one other individual, thus possibly giving rise to behavior that would not be observed in classical multi-virus SIS models. Consequently, to model such higher-order interactions (HOI), as has been correctly pointed out in [4], there is a need to use hypergraphs, <fnref n="1" />, i.e., graphs where an edge can connect more than two nodes. The present paper focuses on the spread of two competing SIS viruses over a (possibly) directed hypergraph.

Indeed, the topic being studied in this paper is not without precedent. Specifically, an SIS model on a hypergraph, albeit with certain restrictions imposed on the network structure, has been studied in [12]; this model relies on the approach proposed in [22], which is a scalar model. The findings in [22] have been extended to the network case (i.e., vector case) in [9]; for global convergence to the DFE, a sufficient condition, simpler than the one in [9, Theorem 5.1, statement ii]), has been provided in [11, Corollary 3]. Note that the papers [22, 12, 9, 11] do not account for the possibility of simultaneous circulation of multiple viruses. This limitation has been partly overcome with a bivirus competitive networked SIS model that admits the possibility of HOI in [10]. Observe that [10, 9] account for the case when HOIs involve simultaneous interactions with at most two other individuals (so groups of three individuals), since,

<footnote n="1">¹Simplicial networks (see [18]) have also been used for studying HOI, see [4].</footnote>
as noted in [9, Section VI], the analysis for the case where there are more than three
individuals in a group is analogous. Therefore, the present paper will also obey the
same restriction as in the setups of [9, 10].

**Paper Contributions.** The present paper analyzes a competitive networked bivirus SIS HOI model. Our analysis significantly extends the existing literature on such models, by establishing general properties of its equilibria and convergence, as well as specific conditions for various equilibria to exist (or not). Our main contribu- tions are as follows:

i) We show that, under the assumption that the aforementioned hypergraph is strongly connected, the model is strongly monotone; see Theorem 3.1.

ii) We show, by means of the Parametric Transversality Theorem from differ-
ential topology, that for generic parameter choices, the system has a finite
number of equilibria, and that the associated vector field zero is nondegener-
ate; see Theorem 3.2. Leveraging Theorems 3.1 and 3.2, it turns out that, for
generic parameter choices, the system converges to some stable equilibrium
point; see Theorem 3.3.

iii) We identify a parameter regime that admits three equilibria (a single-virus endemic equilibrium corresponding to each of the two viruses, and the DFE) that are simultaneously locally exponentially stable; see Proposition 4.1.

iv) A parameter regime that admits at least one coexistence equilibrium, when the single-virus endemic equilibria corresponding to each of the two viruses are unstable, has been identified in [10, Theorem 5.4]. We provide a novel proof for a special case of [10, Theorem 5.4]; see proof of Proposition 5.1. For the same parameter regime, we show that there exists an unstable coexistence equilibrium even when the single-virus endemic equilibria corresponding to each of the two viruses are stable; see Theorem 5.2. Further, we show that for the same parameter regime as covered by Proposition 4.1, there exists an unstable coexistence equilibrium; see Proposition 5.3.

v) We identify a necessary condition for the existence of a coexistence equilib-
rium; see Proposition 5.4.

As a byproduct, we also make the following auxiliary contributions. For the
directed hypergraph associated with our system, we provide two algebraic character-
izations for strong connectivity; see Proposition 2.2 and Corollary 2.4. We identify
a sufficient condition for two single-virus endemic equilibrium (corresponding to each
of the two viruses) to be simultaneously unstable; see Proposition 4.3. For the special
case where there is only one virus circulating in the population, we secure a condition
for local asymptotic stability; see Theorem 4.4.

Some of the results in this paper have appeared in the Proceedings of the 2024 American Control Conference [14]. Note that Proposition 2.2, Corollary 2.4, Corollary 2.5, Lemma 2.9, Theorem 3.1 (and hence Lemma 2.11), Theorem 3.2 (and hence Lemma 2.6), Theorems 3.3 and 4.4, and Propositions 4.1, 5.3 and 5.4 were not included in [14].

**Structure of the Paper.** The paper is organized as follows. We gather all the notations used in the sequel, and recall certain graph-theoretic notions of pertinence in the rest of this section. We introduce the competitive networked bivirus SIS model with HOI, and formally state the problems of interest in Section 2. Our main results are split across Section 3, 4 and 5. Simulations to illustrate our theoretical findings are provided in Section 6. Finally, we summarize our paper, and list a few problems of interest to the wider community in Section 7.
**1.1. Notations.** We denote the set of real numbers by $\mathbb{R}$, and the set of non-negative real numbers by $\mathbb{R}_+$. For any positive integer $n$, we use $[n]$ to denote the set $\{1, 2, ..., n\}$. The $i^{th}$ entry of a vector $x$ is denoted by $x_i$. The element in the $i^{th}$ row and $j^{th}$ column of a matrix $M$ is denoted by $M_{ij}$. We use **0** and **1** to denote the vectors whose entries all equal 0 and 1, respectively, and use $I$ to denote the identity matrix, while the sizes of the vectors and matrices are to be understood from the context. For a vector $x$ we denote the diagonal square matrix with $x$ along the diagonal by $\text{diag}(x)$. For any two real vectors $a, b \in \mathbb{R}^n$ we write $a \ge b$ if $a_i \ge b_i$ for all $i \in [n]$, $a > b$ if $a_i > b_i$ for all $i \in [n]$. Likewise, for any two real matrices $A, B \in \mathbb{R}^{n \times m}$, we write $A \ge B$ if $A_{ij} \ge B_{ij}$ for all $i \in [n], j \in [m]$, and $A > B$ if $A \ge B$ and $A \neq B$. For a square matrix $M$, we use $\sigma(M)$ to denote the spectrum of $M$, $\rho(M)$ to denote the spectral radius of $M$, and $s(M)$ to denote the largest real part among the eigenvalues of $M$, i.e., $s(M) = \max\{\text{Re}(\lambda) : \lambda \in \sigma(M)\}$. We denote a subset by $P \subseteq Q$, a proper subset by $P \subset Q$, and a set difference by $P \setminus Q$.

A real square matrix $A$ is called Metzler if all its off-diagonal entries are nonnegative. If $A(=[a_{ij}]_{n \times n})$ is a nonnegative matrix, then $\rho(A)$ decreases monotonically with a decrease in $a_{ij}$ for any $i, j \in [n]$. The matrix $A$ is reducible if, and only if, there is a permutation matrix $P$ such that $P^\top AP$ is block upper triangular; otherwise, $A$ is said to be irreducible. If a nonnegative $A$ is irreducible, and $Ax = y$ for $x > 0$, then $y > 0$, and $y$ cannot have a zero in every position where $x$ has a zero.

**1.2. Graph theoretical notions.** The notion of hypergraphs will play a crucial role in our analysis. Here, we recall certain concepts related to hypergraphs that will be needed in the sequel.

Note that the hypergraphs of interest in this paper have two important characteristics: they are directed, and the *head vertex set* of any hyperedge is always a single vertex; the latter is a consequence of the restriction imposed in this paper that an individual in node $i$ can simultaneously interact with individuals belonging to at most two different nodes, say node $j$ and node $\ell$. For such a hypergraph $\mathcal{H} = (\mathcal{V}, \mathcal{E})$, where $\mathcal{V}$ is the vertex set, and $\mathcal{E}$ is the set of hyperedges, each hyperedge $e \in \mathcal{E}$ is defined by an ordered pair $e = (S, v)$ where $S$ is a set of *tail vertices* and $v \in \mathcal{V}$ is a *head vertex*. In fact, each hyperedge consists of either a simple edge joining two vertices or a pair of edges leaving from two possibly different vertices but terminating in a single vertex. For the special case where all tail vertex sets have one element, the notions of directed hypergraph and directed graph coincide.

Following the nomenclature in [3], we say that a *hyperpath* from a set of vertices $V_1$ to a single vertex $t$ is an alternating sequence of vertices and hyperedges $\mathcal{P} = \{S, e_1, v_2, e_2, v_3, \dots, e_k, t\}$ such that

i) For each $j = 1, 2, \dots, \ell$, there holds $e_j = (S_j, v_{j+1})$ for some $S_j$ satisfying $S_j \subseteq S \cup \{v_2, v_3, \dots, v_{j-1}\}$, and

ii) For each $j = 1, 2, \dots, \ell + 1$, there holds $v_j \in \cup_{j=2}^\ell S_j$

Thus, each hyperedge used in the hyperpath must be a hyperedge of $\mathcal{H}$, all the tail vertices of a given hyperedge must have been reached previously in the hyperpath, and the head vertex is in a tail vertex of another hyperedge. A directed hypergraph is strongly connected when, between any two vertices, there is a hyperpath.

It is clear that this is an equivalence relation for the vertices, and paralleling the situation applicable for a directed graph, the vertices of a hypergraph can be partitioned into subsets, with each subset constituting vertices that are maximally strongly connected. Further, from such a partition, a condensation graph can be constructed, with
each vertex corresponding to one maximal strongly connected set of vertices of the original hypergraph, and with an edge directed from one vertex of the condensation graph, call it *a*, to another, call it *b*, when the corresponding vertex sets $S_a$ and $S_b$ of the original hypergraph have the property that there is a hyperedge from $S_a$ to any vertex in $S_b$. The condensation graph is an acyclic directed graph.

As such, and assuming the original hypergraph is not strongly connected, the condensation graph has more than one node, and the condensation graph has at least one leaf node. Consequently, in the hypergraph, there is one maximally strongly connected set of vertices from which no other vertex, outside that set, is reachable. If the vertex set is denoted by $S_0$, there is no vertex $v$ outside $S_0$ such that for $H_0$ a subset of $S_0$, ($H_0, v$) is an edge of the hypergraph.

## 2. Problem Formulation.

**2.1. Model.** Consider a network of $n$ nodes. A node represents a well-mixed<fnref n="2" /> population of individuals; the size of the population is fixed. We suppose that there are two viruses, say virus 1 and virus 2, spreading over such a network. Throughout this paper, we will assume that the two aforementioned viruses are competing with each other, which implies that no individual can be simultaneously infected with virus 1 and with virus 2. Through pairwise or higher-order interactions (HOI) (which we describe in more detail below), a healthy individual in node $i$ gets infected with virus 1 (resp. virus 2) due to contact with either other individuals who are infected with virus 1 (resp. virus 2). When, at a given time instant, a single interaction is involved (i.e., between two individuals in node $i$ or between an individual in node $i$ and an individual in node $j$), we say that the infection is due to *pairwise interactions*. An individual in node $i$ might also be infected with virus 1 (resp. virus 2) due to *simultaneous interactions* with infected individuals in nodes $j$ and $\ell$, where, with $j$ being possibly equal to $\ell$ either a) $j = i$, and/or $\ell = i$, or b) $j, \ell$ are neighbors of $i$. Such interactions are referred to as *higher-order interactions* (HOI).

We assume that the pairwise infection (resp. HOI) rate with respect to virus $k$ is the same for all nodes, i.e., $\beta_1^k$ (resp. $\beta_2^k$) for all $i \in [n]$ and $k \in [2]$. An individual infected with virus $k$, depending on its healing rate $\delta_i^k$, recovers from said infection and becomes susceptible to again being infected by virus 1 or by virus 2. For a given node, all individuals in it have the same healing rate with respect to virus $k$; individuals in different nodes possibly have different healing rates. Node $i$ is said to be healthy if all individuals in node $i$ are healthy; otherwise, we say it is infected. Note that, within the same node, it is possible for there to exist simultaneously i) a fraction that is infected with virus 1, and ii) a different fraction that is infected with virus 2.

The spread of the two competing viruses, with the possibility of HOI factored in, can be modeled using a 2-layer hypergraph, say $\mathcal{G}$. The vertices of $\mathcal{G}$ represent the population nodes. The $k^{th}$ layer denotes the contact hypergraph for the spread of virus $k$, where $k = 1, 2$. More specifically, for each layer of $\mathcal{G}$, there exists a hyperedge corresponding to i) each pairwise interaction between individuals in nodes $j$ and $i$, and ii) each HOI jointly between individuals in node $j$ and individuals in node $\ell$, and individuals in node $i$. Let $E^k$ denote the hyperedge set corresponding to the $k^{th}$ layer of $\mathcal{G}$. Let $A^k$ (where $a_{ij}^k \ge 0$) denote the matrix, the entries of which capture the strength of pairwise interactions with respect to virus $k$. For each node $i$, let $B_i^k = [b_{ij}^k]$ denote the matrix whose entries capture the strength of HOI with respect

<footnote n="2">²Well-mixed means that the probability of any two individuals in a node interacting with each other is the same.</footnote>
to virus $k$. The rows and columns of $B_i^k$ are indexed by $j$ and $\ell$, which are the tail vertices of the hyperedge linking those vertices to vertex $i$. Clearly, the matrix must be symmetric. It is, however, easier when $j \neq \ell$, to take each of $(B_i^k)_{j\ell}$ and $(B_i^k)_{\ell j}$ to be one-half the total strength of the interaction attributable to the hyperedge. To summarize, for the $k^{th}$ layer of the hypergraph $G$, there exists a directed hyperedge from node $j$ to node $i$ or from nodes $j, \ell$ to node $i$ if, and only if, $a_{ij}^k > 0$ or $b_{ij\ell}^k > 0$.

Note that $\beta_i^k$ can be considered the base infection rate for virus $k$ for pairwise ($i=1$) and HOI ($i=2$) interactions, while $a_{ij}^k$ and $b_{ij\ell}^k$ are parameters that scale this base rate. We choose this notation for convenience and to be consistent with past literature concerning HOI models ([9, 10, 11]), but point out that our formulation allows complete heterogeneity in infection rates for each node, in line with classical SIS bivirus literature [44].

We use $x_i^k(t)$ to denote the fraction of individuals infected with virus $k$ in agent $i$ at time instant $t$. The evolution of this fraction can, therefore, be represented by the following scalar differential equation [10, Section 5]. Specifically, for $i = 1, 2, \dots, n$, we have

$$ (2.1) \qquad \begin{aligned} \dot{x}_i^1 &= -\delta_i^1 x_i^1 + \beta_1^1(1-x_i^1-x_i^2) \sum_{j=1}^n a_{ij}^1 x_i^1 + \beta_2^1(1-x_i^1-x_i^2) \sum_{j,\ell=1}^n b_{ij\ell}^1 x_j^1 x_\ell^1 \\ \dot{x}_i^2 &= -\delta_i^2 x_i^2 + \beta_1^2(1-x_i^1-x_i^2) \sum_{j=1}^n a_{ij}^2 x_i^2 + \beta_2^2(1-x_i^1-x_i^2) \sum_{j,\ell=1}^n b_{ij\ell}^2 x_j^2 x_\ell^2 \end{aligned} $$

Define $D^1 = \text{diag}(\delta_i^1)$, where $i = 1, 2, \dots, n$, and define $D^2$ analogously. Define $X^1 = \text{diag}(x_i^1)$, where $i = 1, 2, \dots, n$, and define $X^2$ analogously. The matrices $A^k$ and $B_i^k$ for each $i \in [n]$ and $k \in [2]$ are as defined previously. Let $x^k = [x_1^k \quad x_2^k \quad \dots \quad x_n^k]^\top$ for $k = 1, 2$.

Therefore, in vector form, equation (2.1) can be written as:

$$ (2.2) \qquad \begin{aligned} \dot{x}^1 &= -D^1 x^1 + \beta_1^1 (I - X^1 - X^2) A^1 x^1 + \\ & \qquad \beta_2^1 (I - X^1 - X^2) ((x^1)^\top B_1^1 x^1, (x^1)^\top B_2^1 x^1, \dots, (x^1)^\top B_n^1 x^1)^\top \\ \dot{x}^2 &= -D^2 x^2 + \beta_1^2 (I - X^1 - X^2) A^2 x^2 + \\ & \qquad \beta_2^2 (I - X^1 - X^2) ((x^2)^\top B_1^2 x^2, (x^2)^\top B_2^2 x^2, \dots, (x^2)^\top B_n^2 x^2)^\top \end{aligned} $$

Note that $(x^\top B_1 x, x^\top B_2 x, \dots, x^\top B_n x)^\top = \begin{bmatrix} x^\top B_1 x \\ x^\top B_2 x \\ \vdots \\ x^\top B_n x \end{bmatrix}$. We will need the Jacobian of the right side of (2.2), and we note for this purpose, using the symmetric nature of $B_i$, that

$$ \frac{\partial}{\partial x_j} (x^\top B_i x) = 2x^\top B_i e_j $$

and so

$$ (2.3) \qquad \nabla \begin{bmatrix} x^\top B_1 x \\ x^\top B_2 x \\ \vdots \\ x^\top B_n x \end{bmatrix} = 2 \begin{bmatrix} x^\top B_1 \\ x^\top B_2 \\ \vdots \\ x^\top B_n \end{bmatrix}. $$

We note that system (2.2) is a special case of [10, system 5.5] in the following sense: System (2.2) only accounts for the case where virus 1 (resp. virus 2) spread only due to contact with the infected individuals. In contrast, the model in [10] (see [10, system 5.5]) allows for the possibility of the viruses to spread through additional mediums such as a water distribution network, a public transit network, etc.
Remark 2.1. Note that setting $\beta_2^k = 0$ for $k = 1, 2$ results in system (2.2) coinciding with the classic networked bivirus SIS model, see [44] and references therein. Setting, for all $t \in \mathbb{R}_{\ge 0}$, $x^1(t) = \mathbf{0}$ (resp. $x^2(t) = \mathbf{0}$) results in system (2.2) coinciding with the model used for studying the spread of a single virus over hypergraphs in [9].

The Jacobian of system (2.2) evaluated at an arbitrary point, $(x^1, x^2)$, in the state space is as given in (2.4).

$$ (2.4) \qquad J(x^1, x^2) = \begin{bmatrix} J_{11} & J_{12} \\ J_{21} & J_{22} \end{bmatrix}, $$

where

$$ (2.5) \qquad J_{11} = -D^1 + \beta_1^1(I - X^1 - X^2)A^1 - \operatorname{diag}(\beta_1^1 A^1 x^1) + \beta_2^1(I - X^1 - X^2)O_1(x^1) - \beta_2^1 \operatorname{diag}((x^1)^\top B_i^1 x^1)_{i=1,2,...,n} $$

$$ (2.6) \qquad J_{12} = -\operatorname{diag}(\beta_1^1 A^1 x^1) - \beta_2^1 \operatorname{diag}((x^1)^\top B_i^1 x^1)_{i=1,2,...,n} $$

$$ (2.7) \qquad J_{21} = -\operatorname{diag}(\beta_1^2 A^2 x^2) - \beta_2^2 \operatorname{diag}((x^2)^\top B_i^2 x^2)_{i=1,2,...,n} $$

$$ J_{22} = -D^2 + \beta_1^2(I - X^1 - X^2)A^2 - \operatorname{diag}(\beta_1^2 A^2 x^2) + \beta_2^2 \operatorname{diag}((x^2)^\top B_i^2 x^2)_{i=1,2,...,n} $$

$$ (2.8) \qquad \beta_2^2(I - X^1 - X^2)O_2(x^2) - \beta_2^2 \operatorname{diag}((x^2)^\top B_i^2 x^2)_{i=1,2,...,n} $$

The terms $O_1(x^1)$ and $O_2(x^2)$ are as given in (2.9), being obtained from (2.3):

$$ (2.9) \qquad O_1(x^1) = 2 \begin{bmatrix} (x^1)^\top B_1^1 \\ (x^1)^\top B_2^1 \\ \vdots \\ (x^1)^\top B_n^1 \end{bmatrix} \qquad O_2(x^2) = 2 \begin{bmatrix} (x^2)^\top B_1^2 \\ (x^2)^\top B_2^2 \\ \vdots \\ (x^2)^\top B_n^2 \end{bmatrix} $$

We need the following assumption so as to ensure that the model is well-defined.

ASSUMPTION 1. The matrix $D^k$, for $k = 1, 2$, is a positive diagonal matrix. The matrix $A^k$, for $k = 1, 2$, and $B_i^k$ for $i \in [n]$ and $k \in [2]$ is nonnegative.

We define the set $\mathcal{D}$ as follows:

$$ (2.10) \qquad \mathcal{D} := \{(x^1, x^2) \mid x^k \ge \mathbf{0}, k = 1, 2, \sum_{k=1}^2 x^k \le \mathbf{1}\}. $$

It is known that the origin is always an equilibrium for system (2.2), and, under Assumption 1, the set $\mathcal{D}$ is positively invariant; see [10, Lemma 5.1]. The fact that $\mathcal{D}$ is positively invariant guarantees that the state values $x_i^k, k \in [2]$, $i \in [n]$, stay in the $[0, 1]$ interval for all time instants. Given that the states represent fractions of a population node that is infected, values outside the $[0, 1]$ interval do not correspond to physical reality.

The model in System (2.2) has three kinds of equilibria, viz. healthy state or disease-free equilibrium (DFE), $(\mathbf{0}, \mathbf{0})$; single-virus endemic equilibria corresponding to virus $k$ of the form $(\bar{x}^k, \mathbf{0})$, where $\mathbf{0} \ll \bar{x}^k \ll \mathbf{1}$ for $k = 1, 2$; and coexisting equilibria, $(\bar{x}^1, \bar{x}^2)$, where, as we will show in Lemma 2.6, $\mathbf{0} \ll \bar{x}^1, \bar{x}^2 \ll \mathbf{1}$, and $\bar{x}^1 + \bar{x}^2 \ll \mathbf{1}$. We will use the terms “single-virus endemic equilibria” and “boundary equilibria” interchangeably. Note that it is not known if the single-virus endemic equilibria corresponding to virus $k$ are necessarily unique; the papers [9, 10] identify conditions for the existence of a single-virus endemic equilibria corresponding to virus $k$, but remain silent on the question of uniqueness.
**2.2. Strong connectivity for hypergraphs.** Recall that a directed hypergraph is strongly connected when, between any two vertices, there is a hyperpath. In this subsection, we will provide an algebraic characterization of the notion of strong connectedness in hypergraphs, and state and prove several properties of the matrices associated with a hypergraph that would be useful in proving some of the main results of this paper. For ease of exposition, we recall the single-virus networked SIS HOI model, which, as noted previously, is a special case of the model in (2.2).

$$ (2.11) \qquad \begin{aligned} \dot{x}(t) = & -Dx(t) + \beta_1(I - X)Ax(t) + \\ & \beta_2(I - X)(x(t)^\top B_1x(t), x(t)^\top B_2x(t), \dots, x(t)^\top B_nx(t)) \end{aligned} $$

Let $G$ denote the hypergraph associated with system (2.11), and note that the hypergraph $G$ is the same as the 2-layer hypergraph of system (2.2) if the latter consists of a single layer.

With the graph $G$ defined, the following proposition gives an algebraic characterization of the notion of strong connectedness in hypergraphs.

**PROPOSITION 2.2.** *Consider system (2.11). The following two conditions are equivalent:*

i) The hypergraph $G$ is not strongly connected.

ii) There exists a strict subset $S_0$ of the vertex set $V$, assumed to have $|S_0| = n_0 < n$, such that if the vertices are ordered so that all vertices in $S_0$ appear before the remaining $n_1 = n - n_0$ vertices, then the matrices $A$ and $B_i$ have the following structural properties:

$$ (2.12) \qquad A = \begin{bmatrix} A_{00} & A_{01} \\ 0 & A_{11} \end{bmatrix} $$
and for all $i \notin S_0$, i.e., for all $i \in \{n_0 + 1, \dots, n\}$,

$$ (2.13) \qquad B_i = \begin{bmatrix} 0 & (B_i)_{01} \\ (B_i)_{10} & (B_i)_{11} \end{bmatrix} $$

*Proof.* Suppose that the hypergraph $G$ is not strongly connected. Construct $S_0$ as described in subsection 1.2. If the $A_{10}$ (i.e., 2-1) block of $A$ had a nonzero element, this would indicate the existence of a simple edge from one vertex in $S_0$ to a vertex in $V \setminus S_0$, which would contradict the defining property of $S_0$ that there is no join of a single vertex in $S_0$ to a vertex outside $S_0$. Similarly, if the $(B_i)_{00}$ or 1-1 block of $B_i$ had a single nonzero entry for some $i \notin S_0$, this would indicate the existence of a vertex pair of $S_0$ with a hyperedge to a vertex $i$ outside $S_0$. Hence the zero entries of $A$ and $B_i$ are a necessary consequence of the lack of strong connectedness of the hypergraph. The argument is very straightforward to reverse, to complete the demonstration of the equivalence of the two conditions. □

**Proposition 2.2** is a generalization to hypergraphs of the well-known connection between the notion of irreducibility of the adjacency matrix and strong connectivity of a directed graph. There is an additional property of nonnegative irreducible matrices that is useful to generalize to the hypergraph context. Recall first that a nonnegative matrix $A$ is irreducible if and only if for any nonnegative and nonzero vector $x$, the vector $Ax$ has at least one nonzero entry in the same position as where $x$ has a zero entry (see Section 1.1). The following is a generalization of this idea.

**LEMMA 2.3.** *Consider system (2.11). Let $x$ be an arbitrary nonnegative but non-zero $n$-vector and suppose that $j_1, j_2, \dots, j_{n_1}$ denote the indices of the zero entries of $x$. Then under the condition that the entries $(Ax)_{jk}$ of $Ax$ and the scalars $x^\top B_{jk} x$ for*
all $k = 1, 2, \dots, n_1$ are all zero, the hypergraph is not strongly connected. Conversely, if the hypergraph is not strongly connected, there exists such a nonnegative $x$ such that the nominated quantities are zero.

*Proof.* Suppose there exists a nonnegative but nonzero $x$ for which the condition of the Lemma statement is fulfilled. Reorder the entries of $x$ and correspondingly the rows and columns of $A$ and the $B_i$, including the indexing of this matrix set, so that the nonzero entries of $x$ occur before all the zero entries. If the last $n_1$ entries of $Ax$ are all zero, it follows that $A$ has the form of (2.12). Next, notice that $x^\top B_i x$ is nonzero precisely if the subblock matrix of $B_i$ formed from its first $n_0$ rows and columns is nonzero. So if $x^\top B_j x$ is zero for all $j = n_0+1, n_0+2, \dots, n$, equation (2.13) holds for all such $B_j$, i.e., the lack of strong connectivity is established. Conversely, suppose there is a lack of strong connectivity. With the structure defined in Proposition 2.2, any $x$ with the last $n_1$ entries zero will satisfy the conditions of the lemma. $\square$

The following corollary is almost immediate by negation of Lemma 2.3.

**COROLLARY 2.4.** *Consider system (2.11). The hypergraph $G$ is strongly connected if, and only if, for any nonnegative, nonzero $x$ for which entries $j_1, \dots, j_{n_1}$ constitute the set of zero entries, either the corresponding entries of $Ax$ must contain one that is nonzero and/or one of the $x^\top B_{jk} x$ must be nonzero*

Next, we turn our attention to showing that the irreducibility of a matrix including $A$ and the $B_i$ (defined as in system (2.11)) is a consequence of the strong connectivity of the hypergraph. Indeed, we have the following corollary:

**COROLLARY 2.5.** *Consider system (2.11) and suppose that the hypergraph $G$ is strongly connected. Let each $w_i$, for $i = 1, 2, \dots, n$, be a strictly positive vector. Then*

for any positive $\beta_1, \beta_2$, the matrix $\beta_1 A + \beta_2 \begin{bmatrix} w_1^\top B_1 \\ w_2^\top B_2 \\ \vdots \\ w_n^\top B_n \end{bmatrix}$ is irreducible

*Proof.* We shall prove the negation of the claim, by showing that if irreducibility fails, then the hypergraph is not strongly connected. Accordingly, suppose there exists a nonzero $x > 0$ such that if the zero entries of $x$ are in positions $j_1, j_2, \dots, j_{n_1}$, then the corresponding entries $(Ax)_{jk}$ of $Ax$ are zero as are all $w_{jk}^\top B_{jk} x$. Without loss of generality, reorder the vertex numbering so that vertices $j_1, j_2, \dots, j_{n_1}$ become vertices $n_0+1, n_0+2, \dots, n$. This implies that the matrix $A$ satisfies (2.12), while, for $i = n_0+1, n_0+2, \dots, n$, the matrix $B_i$ satisfies (2.13). Therefore, by Proposition 2.2, it follows that the hypergraph is not strongly connected. $\square$

Notice that if either or both of the two matrices appearing in the statement of Corollary 2.5 are multiplied on the left by a diagonal matrix with positive diagonal entries, the conclusion of Corollary 2.5 remains true.

In light of the above, we adopt the following assumption in the present paper.

**ASSUMPTION 2.** *The hypergraph $G^k$ for each $k \in [2]$ is strongly connected.*

**2.3. Problem Statements.** With respect to system (2.2), we aim to answer the following questions in this paper:

i) What is the typical behavior that the trajectories exhibit as time goes to infinity? (For example, can chaos or limit cycles occur, or is there always convergence to an equilibrium?)
ii) Can a parameter regime be identified such that multiple (three, to be precise) equilibria are simultaneously (locally) stable?

iii) Supposing that $x^2(t) = 0$ for all $t$ (resp. $x^1(t) = 0$ for all $t$), provide a sufficient condition for local exponential convergence to a single-virus endemic equilibrium of virus 1 (resp. virus 2).

iv) Supposing that the reproduction number of each virus is greater than one, identify sufficient conditions for the existence of a coexistence equilibrium. Furthermore, establish stability properties of such an equilibrium based on knowledge of the stability properties of the boundary equilibria.

v) Identify a sufficient condition which guarantees that no coexistence equilibrium can exist.

**2.4. Preliminary Lemmas.** In this subsection, we secure certain preliminary results on the nature of equilibria of system (2.1), and recall certain known results on irreducible matrices. These results will play a crucial role in establishing the main findings of this paper.

The following lemma restricts the kinds of equilibria that system (2.2) may possess. The proof is partly inspired from [45, Lemma 3.1].

**LEMMA 2.6.** Consider system (2.2) under Assumptions 1 and 2. If $\bar{x} = (\bar{x}^1, \bar{x}^2) \in D$ is an equilibrium of (2.2), then, for each $k \in [2]$, either $\bar{x}^k = 0$, or $0 \ll \bar{x}^k \ll 1$. Moreover, $\sum_{k=1}^2 \bar{x}^k \ll 1$.

*Proof.* It is clear that $(0,0)$ is an equilibrium of (2.2). Therefore, in the rest of the proof, we will show that any non-zero equilibrium $\bar{x} = (\bar{x}^1, \bar{x}^2)$ of (2.2) must satisfy, for each $k \in [2]$, $0 \ll \bar{x}^k \ll 1$ and $\sum_{k=1}^2 \bar{x}^k \ll 1$. We start off by showing that $\bar{x}^1 + \bar{x}^2 \ll 1$. For any $i \in [n]$, observe that the following is satisfied:

$$ (2.14) \qquad \begin{aligned} \dot{\bar{x}}_i^1 + \dot{\bar{x}}_i^2 = & - \delta_i^1 \bar{x}_i^1 - \delta_i^2 \bar{x}_i^2 + \beta_1^1 (1 - \bar{x}_i^1 - \bar{x}_i^2) \sum_{j=1}^n a_{ij} \bar{x}_j^1 \\ & + \beta_2^1 (1 - \bar{x}_i^1 - \bar{x}_i^2) \sum_{j,l=1}^n b_{ijl} \bar{x}_j^1 \bar{x}_l^1 \\ & + \beta_1^2 (1 - \bar{x}_i^1 - \bar{x}_i^2) \sum_{j=1}^n a_{ij}^2 \bar{x}_j^2 \\ & + \beta_2^2 (1 - \bar{x}_i^1 - \bar{x}_i^2) \sum_{j,l=1}^n b_{ijl}^2 \bar{x}_j^2 \bar{x}_l^2 \end{aligned} $$

Suppose that, for some $i \in [n]$, $\bar{x}_i^1 + \bar{x}_i^2 = 1$. Therefore, since, by Assumption 1, $\delta_i^k > 0$ for $k = 1, 2$, and since $\bar{x}_i^k \in D$, from (2.14), it is clear that $\dot{\bar{x}}_i^1 + \dot{\bar{x}}_i^2 < 0$. However, since by assumption $\bar{x} = (\bar{x}^1, \bar{x}^2)$ is an equilibrium, it must be that $\dot{\bar{x}}_i^1 + \dot{\bar{x}}_i^2 = 0$, which is a contradiction. Therefore, for all $i \in [n]$, $\bar{x}_i^1 + \bar{x}_i^2 < 1$, which implies that $\sum_{k=1}^2 x^k \ll 1$; thus guaranteeing that $\bar{x}^k \ll 1$ for $k = 1, 2$.

We are left to show that $\bar{x}^k \gg 0$ for $k = 1, 2$. To this end, suppose that $\bar{x}^1 > 0$ is an equilibrium point for which there exists at least one (but possibly more) $i \in [n]$ such that $\bar{x}_i^1 = 0$. Note that the equilibrium version of the first line of equation (2.2) yields the following:

$$ (2.15) \qquad \begin{aligned} \dot{\bar{x}}^1 = & - D^1 \bar{x}^1 + \beta_1^1 (I - \bar{X}^1 - \bar{X}^2) A^1 \bar{x}^1 + \\ & \beta_2^1 (I - \bar{X}^1 - \bar{X}^2) ((\bar{x}^1)^\top B_1^1 \bar{x}^1, (\bar{x}^1)^\top B_2^1 \bar{x}^1, \dots, (\bar{x}^1)^\top B_n^1 \bar{x}^1)^\top \end{aligned} $$

By noting that $\bar{x}^1$ is an equilibrium point, and by a suitable rearrangement of terms, we obtain:

$$ (2.16) \qquad \bar{x}^1 = S \bar{x}^1, $$
where

$$ (2.17) \quad S = (D^1)^{-1} \beta_1^1 (I - \bar{X}^1 - \bar{X}^2) A^1 + (D^1)^{-1} \beta_2^1 (I - \bar{X}^1 - \bar{X}^2) ((\bar{x}^1)^\top B_1^1, \dots, (\bar{x}^1)^\top B_n^1)^\top. $$

By Assumption 2, we know that the hypergraph is strongly connected. Therefore, since, by assumption, $\bar{x}^1 > 0$, from Corollary 2.4, it follows that there exists $i \in [n]$ such that if $\bar{x}_i^1 = 0$, then either a) $(A^1 \bar{x}_i^1)_i > 0$ and/or b) $(\bar{x}^1)^\top B_i^1 \bar{x}^1 > 0$. Therefore, since $S$ is a nonnegative matrix (which follows from Assumption 1 and the fact that $I - \bar{X}^1 - \bar{X}^2$ is a positive diagonal matrix), there exists $i \in [n]$ such that $\bar{x}_i^1 = 0$ yet $(S \bar{x}^1)_i > 0$. Note that this contradicts (2.16). Hence, if $\bar{x}^1 > 0$ is an equilibrium point, then it must be that $\bar{x}^1 \gg 0$. By an analogous argument, it can be shown that if $\bar{x}^2 > 0$ is an equilibrium point, then it must be that $\bar{x}^2 \gg 0$, thus completing the proof. □

The following lemmas will be needed for establishing (in)stability of different kinds of equilibria of system (2.2).

**LEMMA 2.7.** [43, Lemma 2.3] Suppose that $M$ is an irreducible Metzler matrix. Then $r = s(M)$ is a simple eigenvalue of $M$, with an eigenvector $\zeta \gg 0$.

**LEMMA 2.8.** [43, Theorem 2.7] Suppose that $N$ is an irreducible nonnegative matrix. Then,

(i) $r = \rho(N)$ is a simple eigenvalue of $N$.

(ii) There is an eigenvector $\zeta \gg 0$ corresponding to the eigenvalue $r$.

(iii) $x > 0$ is an eigenvector only if $Nx = rx$ and $x \gg 0$.

(iv) If $A$ is a nonnegative matrix such that $A < N$, then $\rho(A) < \rho(N)$.

**LEMMA 2.9.** Suppose that $\Lambda$ is a negative diagonal matrix and $N$ is a nonnegative matrix. Let $M$ be the Metzler matrix such that $M = \Lambda + N$. Then, $s(M) < 0$ if and only if $\rho(-\Lambda^{-1}N) < 1$, and $s(M) > 0$ if and only if, $\rho(-\Lambda^{-1}N) > 1$.

*Proof.* Let $Q$ be a matrix of all 1's, and set $N_\epsilon = N + \epsilon Q$. Suppose that $\rho(-\Lambda^{-1}N) < 1$. Then for small enough $\epsilon > 0$, $\rho(-\Lambda^{-1}N_\epsilon) < 1$, due to continuity of the spectral radius with respect to the entries of the matrix. Since $N_\epsilon$ is irreducible, then [32, Proposition 1] implies that $s(M_\epsilon) < 0$, where $M_\epsilon = \Lambda + N_\epsilon$. The spectral abscissa is a monotone function of $\epsilon$ and, since $M \le M_\epsilon$, it must be that $s(M) \le s(M_\epsilon)$ and hence $s(M) < 0$. For the other direction, suppose that $s(M) < 0$. Since $M < M_\epsilon$, we have that $s(M) \le s(M_\epsilon)$. By choosing $\epsilon$ small enough, we obtain $s(M_\epsilon) \le 0$. We separately consider the case where $s(M_\epsilon) < 0$ and the case where $s(M_\epsilon) = 0$. Suppose that $s(M_\epsilon) < 0$, which, given that $M_\epsilon$ is irreducible Metzler, from [32, Proposition 1] implies that $\rho(-\Lambda^{-1}N_\epsilon) < 1$. This, from Lemma 2.8, further implies that $\rho(-\Lambda^{-1}N) < 1$. Suppose that $s(M_\epsilon) = 0$. Therefore, from [32, Proposition 1], it must be that $\rho(-\Lambda^{-1}N_\epsilon) = 1$, which, from Lemma 2.8, further implies that $\rho(-\Lambda^{-1}N) < 1$.

Suppose that $\rho(-\Lambda^{-1}N) > 1$. By noting that $-\Lambda^{-1}N_\epsilon > -\Lambda^{-1}N$, and by applying Lemma 2.8 statement iv), it is clear that $\rho(-\Lambda^{-1}N_\epsilon) > 1$. Since $-\Lambda^{-1}N_\epsilon$ is nonnegative irreducible, from [32, Proposition 1] it must be that $s(M_\epsilon) > 0$. Then for small enough $\epsilon$, $s(M) > 0$. By reversing the arguments, it can be shown that $s(M) > 0$ implies $\rho(-\Lambda^{-1}N) > 1$. □

Lemma 2.9 generalizes [32, Proposition 1], in that it does not insist that the matrix $N$ be irreducible.

**LEMMA 2.10.** [37, Proposition 2] Let $A \in \mathbb{R}^{n \times n}$ be Metzler. Then, $A$ is Hurwitz if, and only if, there exists an $x \in \mathbb{R}^n$ such that $x \gg 0$ and $Ax \ll 0$.
In classical SIS bivirus systems, there are certain limitations on the trajectories;
those starting on the boundary of $\mathcal{D}$ are always inward pointing, see [45, Lemma 3.2].
Hence, it is natural to ask if there are similar restrictions on the trajectories of sys-
tem (2.2). We will address this question in the rest of this subsection. Moreover, we
will be doing so *without* insisting on the matrices $A^k$, for $k = 1, 2$, being irreducible
- this is a sharp contrast from the approach adopted in [45, Lemma 3.2]. We do,
however, invoke an assumption that the hypergraphs $\mathcal{G}^k$ are strongly connected.

LEMMA 2.11. Consider system (2.2) under Assumptions 1 and 2. Suppose that the matrices $B_i^k$ for $i \in [n]$ and $k \in [2]$ are nonnegative. Suppose that the initial conditions satisfy i) $x^k(0) > 0$ for $k \in [2]$, and ii) $(x^1(0), x^2(0)) \in \mathcal{D}$. Then for all finite $t > 0$, $0 \ll x^k(t) \ll 1$ for $k \in [2]$, and $x^1(t) + x^2(t) \ll 1$.

*Proof.* Define $z := \mathbf{1} - x^1 - x^2$ (and $Z = \text{diag}(I - X^1 - X^2)$), and $\hat{B}^i := \text{diag}(\beta_1^i A^i x^i)$. Therefore, (2.2) can be rewritten as:

$$
\begin{align}
\dot{x}^i(t) &= -D^i x^i(t) + \hat{B}^i z(t) + \beta_2^i Z((x^i)^\top B_1^i x^i, (x^i)^\top B_2^i x^i, \dots, (x^i)^\top B_n^i x^i)^\top, \quad i=1,2 \nonumber \\
\dot{z}(t) &= D^1 x^1(t) + D^2 x^2(t) - [\hat{B}^1 + \hat{B}^2]z(t) \nonumber \\
&\quad - \beta_2^1 Z((x^1)^\top B_1^1 x^1, (x^1)^\top B_2^1 x^1, \dots, (x^1)^\top B_n^1 x^1)^\top \nonumber \\
&\quad - \beta_2^2 Z((x^2)^\top B_1^2 x^2, (x^2)^\top B_2^2 x^2, \dots, (x^2)^\top B_n^2 x^2)^\top \tag{2.18}
\end{align}
$$

Suppose that for some $\tau \in \mathbb{R}_{\ge 0}$, and for some $i \in [n]$, $z_i(\tau) = 0$. Since $z_i(\tau) = 1-x_i^1(\tau)-x_i^2(\tau)$, it is immediate that either $x_i^1(\tau) \ne 0$ and/or $x_i^2(\tau) \ne 0$. Therefore, a) since $D^1$ and $D^2$ are positive diagonal matrices, and b) since, by assumption, $z_i(\tau) = 0$, from the second line in equation (2.18), we have that $\dot{z}_i(\tau) > 0$. This implies that, for $t > \tau$ and $t - \tau$ sufficiently small, $z_i(t) > 0$. Note that in the above calculation, the choice of node $i$ is arbitrary; hence, the argument above applies for all $i \in [n]$. Hence, we have that $z_i(t) > 0$ for all $t$ and $i \in [n]$. As a result, we obtain $z(t) \gg 0$, which implies that $x^1(t) + x^2(t) \ll 1$. Hence, we have $x^k(t) \ll 1$ for $k \in [2]$.

We now turn to establishing that $0 \ll x^k(t)$ for all finite $t > 0$. Suppose that for some $t \in \mathbb{R}_{\ge 0}$, $\ell$ entries in the vector $x^1(t)$ equal zero. We label these entries as $i_1, i_2, \dots, i_\ell$. By Assumption 2, we know that the hypergraph is strongly connected. Therefore from Corollary 2.4, since $x^1(t) \ge 0$ for all $t \in \mathbb{R}_{\ge 0}$ (see [10, Lemma 5.1]), it follows that there exists $i \in [n]$ such that $[x^1(t)]_{i_1} = 0$ yet either $[A^1 x^1(t)]_{i_1} > 0$ and/or $(x^1)^\top B_i x > 0$. We assume that this happens for the $i_1^{th}$ entry in $x^1$, i.e., $[x^1(t)]_{i_1} = 0$ yet $[A^1 x^1(t)]_{i_1} > 0$ and/or $((x^1(t))^\top B_i^1 x^1(t) > 0$. The evolution of the infection level for virus 1 in node $i_1$ is as follows:

$$
\begin{equation}
\begin{aligned}
\dot{x}_{i_1}^1(t) = & [-D^1 x^1(t)]_{i_1} + [\beta_1^1 (I - X^1 - X^2) A^1 x^1(t)]_{i_1} \\
& + \beta_2^1 (1 - x^1(t)_{i_1} - x^2(t)_{i_1}) x^1(t)^\top B_{i_1}^1 x^1(t) \\
& > 0,
\end{aligned}
\tag{2.19}
\end{equation}
$$

where the inequality in (2.19) follows by noting that $x_{i_1}(t) = 0$ by assumption, and either $[\beta_1^1 A^1 x^1(t)]_{i_1} > 0$ and/or $x^1(t)^\top B_{i_1}^1 x^1(t) > 0$ as discussed above. This means that there must exist some time instant $t'$, with $t' - t$ not too large, such that $x^1(t')$ has fewer than $\ell$ zero entries. Repeating the argument for all the other zero entries in the vector $x^1(t')$ (indeed, this can be done since the choice of node $i_1$ was arbitrary), we have that $x^1(t) \gg 0$. Analogously, we can prove that $x^2(t) \gg 0$, thus concluding the proof. $\square$
**3. Monotone dynamical systems and competitive bivirus networked SIS models with HOI.** In this section, we establish that the bivirus SIS HOI system is strongly monotone, and provide a general conclusion on the limiting dynamical behavior of the same.

Monotone dynamical systems (MDS) are a class of non-linear systems that are highly appealing in the context of mathematical epidemiology; for an overview of MDS, see [40, 20]. This is mainly because of the fact that MDS, assuming generic parameter values,<fnref n="3" /> have a finite number of equilibria, and converge to a (stable) equilibrium point for all initial conditions except those lying on a measure-zero set. The definition of MDS is as follows. Suppose that $(x_A^1(0), x_A^2(0))$ and $(x_B^1(0), x_B^2(0))$ are two initial conditions in int($\mathcal{D}$) satisfying i) $x_A^1(0) > x_B^1(0)$ and ii) $x_A^2(0) < x_B^2(0)$. Then, if for all $t \in \mathbb{R}_{\ge 0}$ (resp. $t \in \mathbb{R}_{>0}$), i) $x_A^1(t) \ge x_B^1(t)$ (resp. $x_A^1(t) \gg x_B^1(t)$) and ii) $x_A^2(t) \le x_B^2(t)$ (resp. $x_A^2(t) \ll x_B^2(t)$), then system (2.2) is monotone (resp. strongly monotone) [20]. In simpler terms, MDS are endowed with a trajectory ordering property.

**3.1. The bivirus HOI system is strongly monotone.** The following is the first main result of this paper.

**THEOREM 3.1.** Consider system (2.2) under Assumptions 1 and 2. Then, system (2.2) is strongly monotone.

*Proof.* First, recall that a system is monotone (resp. strongly monotone) if the associated Jacobian (or a suitably transformed version of it) is Metzler (resp. irreducible Metzler); see [19, p. 424]. Consider the permutation matrix $P := \text{diag}(I_n, -I_n)$, and observe that $PJ(x^1, x^2)P = [\frac{J_{11}}{-J_{21}} - \frac{J_{12}}{J_{22}}]$. By Assumption 1 and Lemma 2.11, matrices $-J_{12}$ and $-J_{21}$, and the off-diagonal terms in $J_{11}$ and $J_{22}$ are nonnegative; thus, implying that $PJ(x^1, x^2)P$ is Metzler. We are left to show that $PJ(x^1, x^2)P$ is irreducible in the interior of the region of interest.

As a connectivity argument reveals, it is intuitively reasonable, and was argued in [45], that the required irreducibility property will hold provided that $J_{11}$ and $J_{22}$ are separately irreducible, and $J_{12}$ and $J_{21}$ are both nonzero. Note that the Jacobian is as given in (2.4). By disregarding certain diagonal summands in the expression for $J_{11}$, it is evident that the irreducibility of $J_{11}$ is equivalent to the irreducibility of $\hat{J}_{11}$, with

$$ (3.1) \qquad \hat{J}_{11} = \beta_1^1(I - X^1 - X^2)A^1 + 2\beta_2^1(I - X^1 - X^2) \begin{bmatrix} (x^1)^\top B_1^1 \\ (x^1)^\top B_2^1 \\ \vdots \\ (x^1)^\top B_n^1 \end{bmatrix}. \quad \square $$

In the interior of $\mathcal{D}$, the matrix $I - X^1 - X^2$ is nonsingular with positive diagonal elements. Hence the irreducibility follows using Corollary 2.5, The irreducibility of $J_{22}$ follows in the same way.

To see that $J_{12}$ (and by the same argument $J_{21}$) is nonzero, suppose that to the contrary, $J_{12} = 0$. This would imply that for some $x^1 \gg 0$, there holds $A^1 x^1 = 0$ and $(x^1)^\top B_i^1 x^1 = 0$ for all $i$. Together with Lemma 2.3, this implies that all entries of $A^1$ and all entries of all $B_i^1$ are zero, which contradicts Assumption 2.

It is known that under Assumptions 1 and the assumption that the pairwise interaction matrices $A^k$ are irreducible (see Assumption 3 in Section 4 for a formal statement), system (2.2) is strongly monotone; see [10, Theorem 5.5]. Indeed, if there are no

<footnote n="3">³The term “generic” means the following: for all but a set of parameter values that has measure zero. An algebraic or semi-algebraic set defines this set of exceptional values.</footnote>
HOI (i.e., $\beta_2^k = 0$ for $k = 1, 2$), the assumption that $A^k$ is irreducible plays a central role in establishing that classical bi-virus systems are strongly monotone; see [45, Theorem 3.6]. Theorem 3.1 shows that for bivirus systems with HOI the irreducibility of $A^k$ is not required for establishing that the system is strongly monotone.

**3.2. Finiteness of the number of equilibria for generic parameter choices.** In this section, we establish that, for generic parameter choices, system (2.2) has a finite number of equilibria, and furthermore that the associated Jacobian is nonsingular. Note that generic finiteness of the number of equilibria has already been established, by means of an algebraic geometry approach, in [10, Theorem 5.5], but [10, Theorem 5.5] does not establish that, for generic parameter choices, the associated Jacobian is nonsingular. The gap is filled here by use of a different and efficient approach for proving finiteness of equilibria; the tool is the Parametric Transversality Theorem from differential topology, see [30, see p. 145] and [17, see p.68]. We use arguments very much like those in [1]. Essentially because the healthy equilibrium and the single-virus boundary equilibria can be conveniently studied using single-virus techniques, it is easily established that there are no continua of equilibria confined to any boundary, i.e., any continuum of equilibria necessarily includes a continuum of coexistence equilibria. We focus, therefore, on showing that for generic parameter values, such equilibria cannot exist. The main result is as follows:

**THEOREM 3.2.** *Consider system (2.2) under Assumptions 1 and 2. With any fixed matrices $A^k$ and nonnegative $B_i^k$, and the exclusion of a set of values for the entries of $D^1$, $D^2$ of measure zero, the number of coexistence equilibrium points is finite, and the associated vector field zero is nondegenerate, i.e., the associated Jacobian is nonsingular. Similarly, with any fixed $D^1$, $D^2$, and the exclusion of a set of values for the entries of $A^k$, $B^k$ of measure zero, the same properties of equilibrium points hold.*

*Proof:* See Appendix.

**3.3. Generic convergence.** Theorems 3.1 and 3.2 together guarantee that all the assumptions of Hirsch's generic convergence theorem [19] are fulfilled. This allows us to draw conclusions regarding the limiting behavior of system (2.2); we have the following result.

**THEOREM 3.3.** *Consider system (2.2) under Assumptions 1 and 2. For all initial conditions $(x^1(0), x^2(0)) \in D$ except possibly for a set of measure zero, the system (2.2) will converge to an equilibrium. If the system does not converge to an equilibrium, then it is on a nonattractive limit cycle.*

In words, Theorem 3.3 establishes that the typical behavior of system (2.2) is convergence to some equilibrium; this could be healthy, or (one of the possibly many) single-virus boundary equilibria, or a coexistence equilibrium. It further says that limit cycles, if any, are nonattractive. Note that no further complicated behavior is allowed; in particular, chaos can be ruled out, see [41]. Thus, Theorem 3.3 answers question i) raised in Section 2.3. Theorem 3.3 (resp. Theorem 3.2) is more general than [14, Theorem 2] (resp. [14, Theorem 1]), since [14, Theorem 2] (resp. [14, Theorem 1]) require the pairwise interaction matrices $A^k$ to be irreducible.

**4. Existence and local stability of boundary equilibria.** In this section, we identify a parameter regime that permits three equilibria of system (2.2) to be simultaneously locally exponentially stable. Thereafter, for a different parameter regime, we identify a condition for the existence and instability of a boundary equilibrium. Finally, for the case when there is only one virus circulating in the network, we identify
a condition for the existence and local exponential stability of an endemic equilibrium. As we will see in Remark 4.2, the results in this section offer significantly more insight into the spread of viruses over hypergraphs than has hitherto been available.

**PROPOSITION 4.1.** Consider system (2.2) under Assumption 1, and $B_i^k \ge 0$ for all $i \in [n]$ and $k \in [2]$. Define, for $k = 1, 2$, $\mathbf{1}_{B^k} \in \{0, 1\}^n$ by $(\mathbf{1}_{B^k})_i = 1$ if $B_i^k \neq \mathbf{0}$; otherwise $(\mathbf{1}_{B^k})_i = 0$. Suppose that the following conditions are fulfilled for $k = 1, 2$:

a) $\rho(\beta_1^k(D^k)^{-1}A^k) < 1$, and

b) $\min_{i \in [n] \text{ s.t. } B_i^k \neq \mathbf{0}} \left( \frac{\beta_i^k}{\delta_i^k} (A^k \mathbf{1}_{B^k})_i + \frac{\beta_i^k}{2\delta_i^k} \mathbf{1}_{B^k}^\top B_i^k \mathbf{1}_{B^k} \right) > 2$.

Then, the following statements are true:

i) The DFE is locally exponentially stable.

ii) Under Assumption 2, there exist equilibria $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$ where for $k = 1, 2$, there holds $\bar{x}^k \gg \mathbf{0}$, with also $\bar{x}_i^k > \frac{1}{2}$ for $k = 1, 2$, for any $i$ such that $B_i^k \neq \mathbf{0}$.

iii) Under Assumption 2, any such equilibrium point $(\bar{x}^1 \mathbf{0})$ is locally exponentially stable.

iv) Under Assumption 2, any such equilibrium point $(\mathbf{0}, \bar{x}^2)$ is locally exponentially stable.

*Proof. Proof of statement i)*: Note that the Jacobian evaluated at the DFE is as follows:

$$ (4.1) \qquad J(\mathbf{0}, \mathbf{0}) = \begin{bmatrix} -D^1 + \beta_1^1 A^1 & \mathbf{0} \\ \mathbf{0} & -D^2 + \beta_1^2 A^2 \end{bmatrix}. $$

By assumption, $\beta_1^k \rho((D^k)^{-1}A^k) < 1$, for $k = 1, 2$. Therefore, from Lemma 2.9, it must be that $s(-D^k + \beta_1^k A^k) < 0$ for $k = 1, 2$, which, since $J(\mathbf{0}, \mathbf{0})$ is a block diagonal matrix, and since the matrices $-D^1 + \beta_1^1 A^1$ and $-D^2 + \beta_1^2 A^2$ are the only blocks along the main diagonal, implies that $s(J(\mathbf{0}, \mathbf{0})) < 0$. Local exponential stability of the DFE, then, follows from [28, Theorem 4.15 and Corollary 4.3].

*Proof of statement ii)*: Observe that the point $(\bar{x}^1, \mathbf{0})$, with $\bar{x}^1 \gg \mathbf{0}$, is an equilibrium of (2.2) (in fact, a single-virus endemic equilibrium of virus 1) if and only if $\bar{x}^1$ is an endemic equilibrium of (2.11). The same is true when we consider $(\mathbf{0}, \bar{x}^2)$. Thus, we focus on (2.11), and drop the superscript $k$ from $\bar{x}$, $A$, $B_i$, $D$, and $\mathbf{1}_B$ for clarity. Our approach is the same as that in the proof of [9, Theorem 5.1, statement iv)]; the crucial difference is that we establish the claim using Assumption 2, whereas [9, Theorem 5.1, statement iv)] requires a more restrictive assumption, namely that the matrix $A$ be irreducible. First, we introduce the following functions: a) $h_+(z) = \frac{z}{1+z}$ for any $z \in \mathbb{R}_{\ge 0}$; and b) $H_+(y) = (h_+(y_1), h_+(y_2), \dots, h_+(y_n))^\top$ for $y \ge \mathbf{0}$. Define the map $H(x) := H_+(\bar{A}x + (x^\top \bar{B}_1x, x^\top \bar{B}_2x, \dots, x^\top \bar{B}_nx))$, where $\bar{A} := \beta_1 D^{-1}A$, and, for each $i \in [n]$, $\bar{B}_i := \frac{\beta_2}{\delta_i} B_i$. Clearly, a fixed point of map $H(\cdot)$ is an equilibrium of (2.11). Consider the vector $\mathbf{1}_B^\top$, as given in the statement of the proposition but particularized for the single-virus case. Define $Y := \{y \in [0, 1]^n | \frac{1}{2}\mathbf{1}_B \le y \le \mathbf{1}\}$. It turns out that $H(y) > \frac{1}{2}\mathbf{1}_B$. To see this, consider the following: For any $y \in Y$,

$$ (4.2) \qquad \begin{aligned} H(y) &= H_+(\bar{A}y + (y^\top \bar{B}_1y, y^\top \bar{B}_2y, \dots, y^\top \bar{B}_ny)^\top) \\ &\ge H_+(\bar{A}\frac{1}{2}\mathbf{1}_B + \frac{1}{4}((\mathbf{1}_B)^\top \bar{B}_1, \dots, (\mathbf{1}_B)^\top \bar{B}_n)^\top \mathbf{1}_B), \end{aligned} $$

where inequality (4.2) comes from the fact that a) $y \ge \frac{1}{2}\mathbf{1}_B$, and b) the function $h_+(\cdot)$
is monotonic in its argument. Observe that

$$
\begin{equation} \tag{4.3}
\begin{aligned}
& \frac{1}{2}(\bar{A}\mathbf{1}_B)_i + \frac{1}{4}(\mathbf{1}_B)^\top \bar{B}_i \mathbf{1}_B = \frac{1}{2}\frac{\beta_1}{\delta_i} \sum_{j=1}^n a_{ij}(\mathbf{1}_B)_j + \frac{1}{4}\frac{\beta_2}{\delta_i}(\mathbf{1}_B)^\top B_i \mathbf{1}_B \\
& \qquad = \frac{1}{4}\alpha_i,
\end{aligned}
\end{equation}
$$

where $\alpha_i := 2 \frac{\beta_1}{\delta_i} \sum_{j=1}^n a_{ij}(\mathbf{1}_B)_j + \frac{\beta_2}{\delta_i} (\mathbf{1}_B)^\top B_i \mathbf{1}_B$. Notice that (4.3) holds for every $i \in [n]$, and that the left hand side of (4.3) is bounded from below by zero. By statement b) of the proposition hypothesis, it is clear that, for any $i$ such that $B_i \neq 0$, $\alpha_i > 4$. This implies that $h_+(\frac{1}{4}\alpha_i) > \frac{1}{2}$. Consequently, we have that $H(y) \ge H_+(\frac{1}{4} \text{diag}(\alpha_i)_{i=1,2,\dots,n} \mathbf{1}_B) > \frac{1}{2} \mathbf{1}_B$.

Thanks to Assumption 1, we can be assured that $(\bar{A}y)_i + y^\top \bar{B}_i y \ge 0$, which implies that $h_+((\bar{A}y)_i + y^\top \bar{B}_i y) \le 1$. This, in turn, implies that $H(y) \le \mathbf{1}$. This, coupled with the fact that $H(y) > \frac{1}{2} \mathbf{1}_B$ as established above, means that $H: Y \to Y$. Therefore, using Brouwer's fixed-point theorem [38, Theorem 4.5], we conclude that there exists a fixed point $\bar{y} \in Y$ such that $H(\bar{y}) = \bar{y}$.

Next, we show that $\bar{y} \gg 0$. To this end, note that any fixed point of $H$, say $\bar{y}$, must satisfy the following

$$
(4.4) \qquad \frac{(\bar{A}\bar{y})_i + \bar{y}^\top \bar{B}_i \bar{y}}{1 + (\bar{A}\bar{y})_i + \bar{y}^\top \bar{B}_i \bar{y}} = \bar{y}_i, \quad i = 1, 2, \dots, n
$$

Suppose, to obtain a contradiction, that $\bar{y}_{j_1}, \bar{y}_{j_2}, \dots, \bar{y}_{j_k}$ are the zero entries of $\bar{y}$.
Then the left side of (4.4) is zero for $i = j_1, j_2, \dots, j_k$ and there must then hold

$$
\bar{y}^{\top} \bar{B}_{j_1} \bar{y} = \bar{y}^{\top} \bar{B}_{j_2} \bar{y} = \dots = \bar{y}^{\top} \bar{B}_{j_k} \bar{y} = 0
$$

and also

$$
(\bar{A}\bar{y})_{j_1} = (\bar{A}\bar{y})_{j_2} = \dots = (\bar{A}\bar{y})_{j_k} = 0
$$

By Corollary 2.4, it then follows that the hypergraph is not strongly connected, which contradicts the proposition hypothesis (namely, that Assumption 2 holds). This implies that $\bar{y} \gg 0$.

We will now show that $\bar{y}_i > \frac{1}{2}$ for $i$ such that $B_i \neq 0$. To this end, recall that by statement b) of the proposition hypothesis, we have the following:

$$
(4.5) \qquad \min_{i \in [n] \text{ s.t. } B_i \neq 0} \left( \frac{\beta_1}{\delta_i} (A \mathbf{1}_B)_i + \frac{\beta_2}{2\delta_i} \mathbf{1}_B^\top B_i \mathbf{1}_B \right) > 2
$$

Now the mapping $y \to H(y)$ maps any $y \in Y$ with $\frac{1}{2}\mathbf{1}_B \le y \le \mathbf{1}_n$ into $Y$, and the
equilibrium point $\bar{y}$ is contained in $Y$. Invoking the lower limit, we observe that for
such a point,

$$
(4.6) \qquad (\bar{A}\bar{y})_i + \bar{y}^\top \bar{B}_i \bar{y} \ge (\bar{A}\frac{1}{2}\mathbf{1}_B)_i + \frac{1}{4}\mathbf{1}_B^\top \bar{B}_i \mathbf{1}_B
$$

But in view of (4.5), the inequality in (4.6) means that $(\bar{A}\bar{y})_i + \bar{y}^\top \bar{B}_i \bar{y} > 1$; this, in combination with (4.4), ensures that $\bar{y}_i > \frac{1}{2}$.

*Proof of statement iii)*: Our approach is similar to that in the proof of [9, Theorem 5.1, statement v)]. Consider the equilibrium point $(\bar{x}^1, \mathbf{0})$, and observe that the Jacobian
evaluated at this equilibrium is as follows:

$$ (4.7) \qquad J(\bar{x}^1, \mathbf{0}) = \begin{bmatrix} \bar{J}_{11} & \bar{J}_{12} \\ \mathbf{0} & \bar{J}_{22} \end{bmatrix}, $$

where

$$ \begin{align*} \bar{J}_{11} &= -D^1 + \beta_1^1(I - \bar{X}^1)A^1 - \operatorname{diag}(\beta_1^1 A^1 \bar{x}^1) + \\ & \beta_2^1(I - \bar{X}^1)O_1(\bar{x}^1) - \beta_2^1 \operatorname{diag}((\bar{x}^1)^\top B_i^1 \bar{x}^1)_{i=1,2,\dots,n} \\ \bar{J}_{12} &= -\operatorname{diag}(\beta_1^1 A^1 \bar{x}^1) - \beta_2^1 \operatorname{diag}((\bar{x}^1)^\top B_i^1 \bar{x}^1)_{i=1,\dots,n} \\ \bar{J}_{22} &= -D^2 + \beta_1^2(I - \bar{X}^1)A^2. \end{align*} $$

The term $O_1(\bar{x}^1)$ is as defined in (2.9). We will establish the exponential stability of the 11 and 22 blocks (i.e., $\bar{J}_{11}$ and $\bar{J}_{22}$) separately. Observe that

$$ \bar{J}_{11} = -D^1 + \beta_1^1(I - \bar{X}^1)A^1 - \operatorname{diag}(\beta_1^1 A^1 \bar{x}^1) + \beta_2^1(I - \bar{X}^1)O_1(\bar{x}^1) - \beta_2^1 \begin{bmatrix} (\bar{x}^1)^\top B_1^1 \bar{x}^1 & \cdots & (\bar{x}^1)^\top B_n^1 \bar{x}^1 \end{bmatrix}. $$

Define the summands

$$ \begin{align*} Q_1 &:= -D^1 + \beta_1^1(I - \bar{X}^1)A^1 + \beta_2^1(I - \bar{X}^1) \begin{bmatrix} (\bar{x}^1)^\top B_1^1 \\ \vdots \\ (\bar{x}^1)^\top B_n^1 \end{bmatrix}, \quad \text{and} \\ Q_2 &:= \beta_2^1(I - \bar{X}^1) \begin{bmatrix} (\bar{x}^1)^\top (B_1^1)^\top \\ \vdots \\ (\bar{x}^1)^\top (B_n^1)^\top \end{bmatrix} - \operatorname{diag}(\beta_1^1 A^1 \bar{x}^1) - \beta_2^1 \begin{bmatrix} (\bar{x}^1)^\top B_1^1 \bar{x}^1 \\ \vdots \\ (\bar{x}^1)^\top B_n^1 \bar{x}^1 \end{bmatrix}. \end{align*} $$

It is immediate that $\bar{J}_{11} = Q_1 + Q_2$, which implies that $\bar{J}_{11} \bar{x}^1 = Q_1 \bar{x}^1 + Q_2 \bar{x}^1$. Since $\bar{x}^1$ is a single-virus endemic equilibrium corresponding to virus 1, by taking recourse to the equilibrium version of the first line of equation (2.2), it is clear that $Q_1 \bar{x}^1 = \mathbf{0}$. Hence, $\bar{J}_{11} \bar{x}^1 = Q_2 \bar{x}^1$.

Note that

$$ (4.8) \qquad Q_2 \bar{x}^1 = \beta_2^1 (I - \bar{X}^1) \begin{bmatrix} (\bar{x}^1)^\top (B_1^1)^\top \bar{x}^1 \\ \vdots \\ (\bar{x}^1)^\top (B_n^1)^\top \bar{x}^1 \end{bmatrix} - \operatorname{diag}(\beta_1^1 A^1 \bar{x}^1) \bar{x}^1 - \beta_2^1 \begin{bmatrix} (\bar{x}^1)^\top B_1^1 \bar{x}^1 \\ \vdots \\ (\bar{x}^1)^\top B_n^1 \bar{x}^1 \end{bmatrix} \bar{x}^1. $$

Denote by $(Q_2 \bar{x}^1)_i$ the $i^{th}$ entry of the vector $Q_2 \bar{x}^1$. Therefore, in view of (4.8), we have the following:

$$ (4.9) \qquad (Q_2 \bar{x}^1)_i = -\beta_1^1 \left( \sum_{j=1}^n a_{ij}^1 \bar{x}_j^1 \right) \bar{x}_i^1 + \beta_2^1 (1 - 2\bar{x}_i^1) ((\bar{x}^1)^\top B_i^1 \bar{x}^1) $$

We consider its sign, under two circumstances.

**Case 1:** Suppose that $B_i^1 = \mathbf{0}$. Note that by Assumption 2 the hypergraph $\mathcal{G}$ is strongly connected, which, from Corollary 2.5 coupled with Lemma 2.11, implies that

the matrix $\beta_1^1 A^1 + \begin{bmatrix} (x^1)^\top B_1^1 \\ (x^1)^\top B_2^1 \\ \vdots \\ (x^1)^\top B_n^1 \end{bmatrix}$ is irreducible. Consequently, the hypothesis that $B_i^1 = \mathbf{0}$

implies that the $i^{th}$ row of matrix $A^1$ has a positive sum. This, since $\bar{x}_i^1 > 0$, implies
that $\beta_1^1(\sum_{j=1}^n a_{ij}\bar{x}_j^1)\bar{x}_i^1 > 0$. Thus, $(Q_2\bar{x}^1)_i < 0$.

**Case 2:** Suppose secondly that $B_i^1 \neq 0$. Therefore, from statement ii), it follows that $\bar{x}_i^1 > \frac{1}{2}$. Consequently, $1 - 2\bar{x}_i^1 < 0$; thus implying that $(Q_2\bar{x}^1)_i < 0$.

Note that, for the analysis above, the choice of index $i$ was arbitrary; therefore, again, we have $(Q_2\bar{x}^1)_i < 0$ for all $i \in [n]$. Hence, since $\bar{J}_{11}\bar{x}^1 = Q_2\bar{x}^1$, it follows that $(\bar{J}_{11}\bar{x}^1)_i < 0$ for all $i \in [n]$. Note that Assumptions 1 and 2, together with Corollary 2.5 and Lemma 2.11 guarantee that the matrix $\bar{J}_{11}$ is an irreducible Metzler matrix; the reasons are as detailed in the proof of Theorem 3.1. Therefore, from Lemma 2.10, it must be that the matrix $\bar{J}_{11}$ is Hurwitz.

Turning to the 22 block of the Jacobian (i.e., matrix $\bar{J}_{22}$), consider the matrices $\beta_1^2(D^2)^{-1}A^2$ and $\beta_1^2(D^2)^{-1}(I - \bar{X}^1)A^2$. From Assumption 1, it is clear that $\beta_1^2(D^2)^{-1}A^2$ is a nonnegative matrix. Since, from statement ii), $\bar{x}^1$ satisfies $0 \ll \bar{x}^1 \ll 1$, it is also clear that $\beta_1^2(D^2)^{-1}(I - \bar{X}^1)A^2$ is a nonnegative matrix. Furthermore, we immediately obtain $\beta_1^2(D^2)^{-1}(I - \bar{X}^1)A^2 < \beta_1^2(D^2)^{-1}A^2$. Therefore, since the spectral radius of a nonnegative matrix decreases monotonically with a decrease in any entry of said matrix (see Section 1.1), it follows that $\rho(\beta_1^2(D^2)^{-1}(I - \bar{X}^1)A^2) \le \rho(\beta_1^2(D^2)^{-1}A^2)$. By assumption, $\rho(\beta_1^2(D^2)^{-1}A^2) < 1$, which implies that $\rho(\beta_1^2(D^2)^{-1}(I - \bar{X}^1)A^2) < 1$, and consequently, from Lemma 2.9, we have that $s(-D^2 + (I - \bar{X}^1)A^2) < 0$. Therefore, since $J(\bar{x}^1, 0)$ is block upper triangular, and since we have already established that $\bar{J}_{11}$ is Hurwitz, it follows that $s(J(\bar{x}^1, 0)) < 0$. Local exponential stability of $(\bar{x}^1, 0)$, then, follows from [28, Theorem 4.15 and Corollary 4.3].

*Proof of statement iv):* The proof is same as that of statement iii), up to a suitable adjustment of notation. $\square$

Proposition 4.1 answers question ii) raised in Section 2.3. We have the following remark.

**Remark 4.2.** Proposition 4.1 highlights an interesting phenomenon that one observes in bivirus spread over hypergraph but not in classical competitive bivirus settings: identification of a parameter regime that permits three equilibria, namely the DFE, a single-virus endemic equilibrium corresponding to virus 1, and a single-virus endemic equilibrium corresponding to virus 2, to be simultaneously locally exponentially stable. In the classical bivirus system, it is impossible for the DFE and an endemic equilibrium to be simultaneously stable; see [44, Section E].

Observe that Proposition 4.1 guarantees existence of boundary equilibria for the case when $\rho(\beta_1^k(D^k)^{-1}A^k) < 1$. It is quite natural to ask if one is assured of existence even if the spectral radius inequality was flipped. In order to answer this, we require a stronger assumption, which is as follows:

**ASSUMPTION 3.** The matrix $A^k$, for $k = 1, 2$, is irreducible.

Notice that, due to Proposition 2.2, it is clear Assumption 3 implies Assumption 2; the converse, however, is not necessarily true. With Assumption 3 in place, we have the following proposition.

**PROPOSITION 4.3.** Consider system (2.2) under Assumptions 1 and 3. Suppose that, for all $k \in [2]$, $\rho(\beta_k^k(D^k)^{-1}A^k) > 1$. Then system (2.2) has at least three equilibria, namely the DFE, which is unstable; a single-virus endemic equilibrium corresponding to virus 1 $(\bar{x}^1, 0)$; and a single-virus endemic equilibrium corresponding to virus 2 $(0, \bar{x}^2)$. Furthermore, if $s(-D^i + \beta_1^i(I - \bar{X}^k)A^i) > 0$ for each $i, k \in [2]$ such
that $i \neq k$, then the equilibrium points $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$ are unstable.

*Proof.* Observe that the DFE is always an equilibrium of system (2.2). Suppose that for some $k \in [2]$, $\rho(\beta_1^k (D^k)^{-1} A^k) > 1$, then from (4.1) it is clear that $s(J(\mathbf{0}, \mathbf{0})) > 0$. This, in view of [28, Theorem 4.7, statement ii)], guarantees that the DFE is unstable. Further, from [9, Theorem 5.1, statement vii)] we know that there exists an endemic equilibrium, $\bar{x}^k$, where $\mathbf{0} \ll \bar{x}^k \ll \mathbf{1}$. Since, by assumption, $\rho(\beta_1^k (D^k)^{-1} A^k) > 1$ for all $k \in [2]$, it is also immediate that there exist equilibria, $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$, where $\mathbf{0} \ll \bar{x}^k \ll \mathbf{1}$, for $k=1,2$.
Observe that the Jacobian evaluated at the point $(\bar{x}^1, \mathbf{0})$ is as given in (4.7). Further, note that $J(\bar{x}^1, \mathbf{0})$ is block upper triangular. By assumption, $s(-D^1+\beta_1^1(I-\bar{X}^2)A^1) > 0$, which implies that $s(J(\bar{x}^1, \mathbf{0})) > 0$. Consequently, instability of $(\bar{x}^1, \mathbf{0})$ follows from [28, Theorem 4.7, statement ii)]. The instability of $(\mathbf{0}, \bar{x}^2)$ can be shown analogously, thus completing the proof. $\square$

Note that Proposition 4.3 does not say anything about the situation when one of the spectral radii is greater than one, and the other is less than one. We discuss this here. Suppose that $s(-D^2 + \beta_1^2(I - \bar{X}^1)A^2) > 0$ and $s(-D^1 + \beta_1^1(I - \bar{X}^2)A^1) < 0$. Then, since by assumption $s(-D^2 + \beta_1^2(I - \bar{X}^1)A^2) > 0$, from Proposition 4.3 it is clear that $(\bar{x}^1, \mathbf{0})$ is unstable. By assumption, $s(-D^1 + \beta_1^1(I - \bar{X}^2)A^1) < 0$. However, whether or not $s(-D^1 + \beta_1^1(I - \bar{X}^2)A^1) < 0$ guarantees (local exponential) stability of $(\mathbf{0}, \bar{x}^2)$ is an open question. So no conclusions can be drawn regarding the stability of $(\mathbf{0}, \bar{x}^2)$.

While Proposition 4.3 assures existence of single-virus endemic equilibria given satisfaction of a spectral radius condition, without further information it remains silent on two aspects of said equilibria. First, it is not clear if there is a unique single-virus endemic equilibrium for a given virus $k$ (for the classic bivirus network SIS model, it is known that the single-virus endemic equilibrium is unique). Second, the stability of said equilibria, even in a local sense, is not known. For the single-virus case, a sufficient condition for global stability of the endemic equilibrium has been identified in [9, Theorem 5.1, statement viii)]. Regardless, the following caveat with respect to said result should be noted: [9, Theorem 5.1, statement viii] is reliant on that fact that for a given $\beta_1$, there exists a small enough $\beta_2$ such that for the resulting single-virus system the endemic equilibrium is globally exponentially stable. This is clearly *not* the case if one were to seek a sufficient condition where both the parameters $\beta_1$ and $\beta_2$ are arbitrary. It turns out, however, that one can identify a sufficient condition for local exponential attractivity of the endemic equilibrium in the single-virus case, even when the parameters $\beta_1$ and $\beta_2$ are arbitrary. We have the following theorem.

**THEOREM 4.4.** *Consider system (2.11) under Assumptions 1 and 3. Suppose that the parameters associated with system (2.11) are generic. If $\beta_1\rho(D^{-1}A) > 1$, then there exists an equilibrium point $\bar{x}$ such that $\bar{x} \gg \mathbf{0}$. Furthermore, the point $\bar{x}$ is asymptotically stable.*

*Proof.* First, observe that [10, Theorem 4.7] implies that system (2.11) is an irreducible monotone system, and that, for generic parameter choices, for almost all initial conditions, the dynamics converge to a stable equilibrium point. Note that for system (2.11) the origin is an equilibrium, and there are no other equilibria on the boundary of the set $[0, 1]^n$; see [9, Lemma 5.1, statement iii)]. By assumption, $\beta_1\rho(D^{-1}A) > 1$, which, from Lemma 2.9 implies that $s(-D+\beta_1A) > 0$. Consequently, from the analysis in the proof of Proposition 4.1, it is clear that the origin is unstable. Since $\beta_1\rho(D^{-1}A) > 1$, it follows from [9, Theorem 5.1, statement vii)] that there exists
an equilibrium $\bar{x}$ such that $\bar{x} \gg 0$. Thus, all the hypotheses of [40, Theorem 2.8] are satisfied, and consequently we conclude that $s(J(\bar{x})) \le 0$.

**Theorem 3.2** implies that, for generic parameter matrices, the Jacobian of system (2.11), at any equilibrium point, is nonsingular. This means that $J(\bar{x})$ cannot have an eigenvalue with zero real part and zero imaginary part. It can be easily verified that $J(\bar{x})$ is an irreducible Metzler matrix. Therefore, from Lemma 2.7, it follows that the eigenvalue with the largest real part does not have any imaginary parts. Consequently, the possibility that, for generic parameter choices, $s(J(\bar{x})) = 0$ can be ruled out. It follows that [40, Theorem 2.8, statement iv)], coupled with the fact that the set $D$ is positively invariant (see [10, Lemma 5.1]), guarantees that the equilibrium point $\bar{x}$ attracts at least all initial conditions $x$ such that $0 \le x \ll \bar{x}$. Therefore, it is clear that $\bar{x}$ is locally asymptotically stable. $\square$

**5. (Non)existence of Coexistence equilibria.** This section deals with the establishment of sufficient conditions for the existence (resp. nonexistence) of coexistence equilibria. To this end, we consider two parameter regimes, viz. for $k = 1, 2$, i) $s(-D^k + \beta_i^k A^k) > 0$, and ii) $s(-D^k + \beta_i^k A^k) < 0$, and show that for each of these, there exists at least one coexistence equilibrium.

**PROPOSITION 5.1.** *Consider system (2.2) under Assumptions 1 and 3. Let $(\bar{x}^1, 0)$ and $(0, \bar{x}^2)$ denote a single-virus endemic equilibrium corresponding to virus 1 and virus 2, respectively. Suppose that the following conditions are satisfied:*

i) $s(-D^1 + \beta_1^1 A^1) > 0$;

ii) $s(-D^2 + \beta_1^2 A^2) > 0$;

iii) $s(-D^1 + \beta_1^1 (I - \bar{X}^2) A^1) > 0$; and

iv) $s(-D^2 + \beta_1^2 (I - \bar{X}^1) A^2) > 0$.

Then, there exists at least one equilibrium of the form $(\hat{x}^1, \hat{x}^2)$ such that $0 \ll \hat{x}^1, \hat{x}^2 \ll 1$ and $\hat{x}^1 + \hat{x}^2 \ll 1$.

Prior to proving the claim in Proposition 5.1, we need the following background material. In line with the terminology of [21], given an equilibrium point of system (2.2), it is classified as saturated or unsaturated. An equilibrium is saturated (resp. strictly saturated) if the diagonal block of the Jacobian corresponding to the zero entries of said equilibrium has no eigenvalues with positive real part (resp. has every eigenvalue to be strictly less than zero) and unsaturated otherwise [21]. A boundary equilibrium of (2.2) is saturated if and only if said boundary equilibrium is locally exponentially stable; this follows immediately by noting the structure of the Jacobian matrix, evaluated at a boundary equilibrium, see (4.7). The definition also implies that every fixed point in the interior of $D$, irrespective of its stability properties, is saturated [21]; therefore, from Lemma 2.6, we have that every coexistence equilibrium of system (2.2) is saturated.

**Proof.** Assumptions i) and ii) of Proposition 5.1 guarantee existence of boundary equilibria, $(\bar{x}^1, 0)$ and $(0, \bar{x}^2)$; see Proposition 4.3. Observe that [10, Lemma 5.1] guarantees that, for each $k \in [2]$, $x^k(0) \ge 0$ implies that $x^k(t) \ge 0$ for all $t \in \mathbb{R}_{>0}$, and that the set $D$ (which is compact) is forward invariant. Therefore, from [21, Theorem 2], it follows that system (2.2) has at least one saturated fixed point. There are two cases to consider.

**Case 1:** Suppose that the aforementioned saturated fixed point is in the interior of $D$. Note that any fixed point in the interior of $D$ is of the form $(\hat{x}^1, \hat{x}^2)$, where $0 \ll (\hat{x}^1, \hat{x}^2) \ll 1$, thus implying that $(\hat{x}^1, \hat{x}^2)$ is a coexistence equilibrium. From Lemma 2.6, it must necessarily satisfy $0 \ll (\hat{x}^1, \hat{x}^2) \ll 1$, and $\hat{x}^1 + \hat{x}^2 \ll 1$.
Case 2: Suppose, but we will demonstrate a contradiction, that there are no fixed points in the interior of $\mathcal{D}$. This implies that there must be a saturated fixed point on the boundary of $\mathcal{D}$ [21]. Therefore, at least one of the single-virus boundary equilibria is saturated.

However, from Proposition 4.3, it is clear that assumptions iii) and iv) guarantee that the boundary equilibria are unstable; thus implying that they are unsaturated, and the contradiction is obtained. □

Proposition 5.1 is implied by [10, Theorem 5.4], which, assuming $\beta_2^k = 0$ for $k = 1, 2$, is the same as [24, Theorem 5], [13, Theorem 5.4]. The proof of [10, Theorem 5.4] is lengthy, since it primarily relies on fixed point mapping, Perron-Frobenius theory, etc. Our proof is significantly shorter, since we exploit the advantages that the notion of MDS endows our system with.

Note further that Proposition 5.1 makes no comment on the number of coexistence equilibria. This gap may be partially addressed as follows. Thanks to Theorem 3.2, one could perhaps leverage [21, Theorem 2] so as to obtain a lower bound on the number of coexistence equilibria; this is very much the situation for classic bivirus networked SIS models; see [1, Corollary 3.10, statement 1]. Thereafter, since system (2.2) is monotone, the properties of MDS can be exploited to possibly conclude that there must exist a locally exponentially stable coexistence equilibrium. A detailed treatment of this issue is beyond the scope of the present paper.

For the same parameter regime as in Proposition 5.1, we now present a different condition that ensures the existence of a coexistence equilibrium point $(\hat{x}^1, \hat{x}^2)$, but, differently from the situation in Proposition 5.1, we assume that both boundary equilibria are stable.

**THEOREM 5.2.** Consider system (2.2) under Assumptions 1 and 3. Suppose that the parameters associated with system (2.2) are generic. Let $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$ denote a single-virus endemic equilibrium corresponding to virus 1 and virus 2, respectively. Suppose that the following conditions are satisfied:

i) $s(-D^1 + \beta_1^k A^1) > 0$;

ii) $s(-D^2 + \beta_1^k A^2) > 0$;

Suppose that both $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$ are locally exponentially stable. Then there exists at least one equilibrium of the form $(\hat{x}^1, \hat{x}^2)$ such that $\mathbf{0} \ll \hat{x}^1, \hat{x}^2 \ll \mathbf{1}$ and $\hat{x}^1 + \hat{x}^2 \ll \mathbf{1}$, such that $(\hat{x}^1, \hat{x}^2)$ is unstable.

*Proof.* By assumption, $s(-D^k + \beta_1^k A^k) > 0$ for $k = 1, 2$. Therefore, from Proposition 4.3, it follows that there exists a single-virus endemic equilibrium corresponding to virus 1, $\bar{x}^1 \gg \mathbf{0}$, and a single-virus endemic equilibrium corresponding to virus 2, $\bar{x}^2 \gg \mathbf{0}$. By assumption, both $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$ are locally exponentially stable. The condition $s(-D^1 + \beta_1^k A^1) > 0$ implies that the origin is unstable; this can be observed from the proof of statement i) in Proposition 4.1. We are left to show that the stable manifold of the origin does not lie in the interior of $\mathcal{D}$. We will do so by relying on the proof technique of [1, Lemma 3.8]. It suffices to show that for the (linear) system

$$ (5.1) \qquad \begin{bmatrix} \dot{x}^1 \\ \dot{x}^2 \end{bmatrix} = \begin{bmatrix} -D^1 + \beta_1^k A^1 & \mathbf{0} \\ \mathbf{0} & -D^2 + \beta_1^k A^2 \end{bmatrix} \begin{bmatrix} x^1 \\ x^2 \end{bmatrix} $$

no trajectory starting in the interior of $\mathcal{D}$ converges to the origin. First, consider $x^1(t)$. Let $w^\top$ be the positive left eigenvector associated with $s(-D^1 + \beta_1^k A^1)$ such that all its entries sum to one. Define $z := w^\top x^1$, and observe that $\dot{z} = w^\top \dot{x}^1$, which, from (5.1), further implies that $\dot{z} = s(-D^1 + \beta_1^k A^1)z$. Since, by assumption, $s(-D^1 + \beta_1^k A^1) > 0$, it is clear that the projection onto $w$ (which is a positive vector)
of the points of (5.1) in the interior of $\mathcal{D}$ is away from $x^1 = \mathbf{0}$. An analogous argument can be made for $x^2(t)$, since, by assumption, $s(-D^2+\beta_1^2 A^2) > 0$. Therefore, the stable manifold of the origin does not lie in the interior of $\mathcal{D}$. Consequently, since we know that system (2.2) is monotone and the monotone condition $\bar{x}^1 \gg \mathbf{0}$, $\bar{x}^2 \gg \mathbf{0}$ relates the two exponentially stable equilibrium points, from [40, Proposition 2.9] it follows that there exists an equilibrium point of the form $(\hat{x}^1, \hat{x}^2)$ such that $\mathbf{0} \ll \hat{x}^1, \hat{x}^2 \ll \mathbf{1}$ and $\hat{x}^1 + \hat{x}^2 \ll \mathbf{1}$. Furthermore, the point $(\hat{x}^1, \hat{x}^2)$ satisfies $s(J(\hat{x}^1, \hat{x}^2)) \ge 0$. From Theorem 3.2, we know that the Jacobian, at any equilibrium point, is nonsingular. This means that $J(\hat{x}^1, \hat{x}^2)$ cannot have an eigenvalue with zero real part and zero imaginary part. Observe that $J(\hat{x}^1, \hat{x}^2)$ is an irreducible Metzler matrix. Therefore, from Lemma 2.7, it follows that the eigenvalue with the largest real part does not have any imaginary parts. Consequently, $s(J(\hat{x}^1, \hat{x}^2)) > 0$, which implies that the equilibrium point $(\hat{x}^1, \hat{x}^2)$ is unstable. $\square$

**Proposition 5.1** and **Theorem 5.2** partially answer question iv) raised in Section 2.3.

We now turn our attention to a different parameter regime, namely $s(-D^k + \beta_1^k A^k) < 0$ for $k \in [2]$, and provide a sufficient condition which guarantees the existence of an unstable coexistence equilibrium. We have the following result.

**PROPOSITION 5.3.** Consider system (2.2) with generic parameter matrices, under Assumptions 1 and 2. Suppose further that for each $k \in [2]$ there exists some $i \in [n]$ such that $B_i^k \neq \mathbf{0}$. Define, for $k = 1, 2$, $\mathbf{1}_{B^k} \in \{0, 1\}^n$ by $(\mathbf{1}_{B^k})_i = 1$ if $B_i^k \neq \mathbf{0}$; otherwise $(\mathbf{1}_{B^k})_i = 0$. Suppose that the following conditions are fulfilled for $k = 1, 2$:

a) $\rho(\beta_1^k (D^k)^{-1} A^k) < 1$, and

b) $\min_{i \in [n] \text{ s.t. } B_i^k \neq \mathbf{0}} \left( \frac{\beta_i^k}{\delta_i^k} (A^k \mathbf{1}_{B^k})_i + \frac{\beta_i^k}{2\delta_i^k} \mathbf{1}_{B^k}^\top B_i \mathbf{1}_{B^k} \right) > 2$.

Then there exists at least one equilibrium of the form $(\hat{x}^1, \hat{x}^2)$ such that $\mathbf{0} \ll \hat{x}^1, \hat{x}^2 \ll \mathbf{1}$ and $\hat{x}^1 + \hat{x}^2 \ll \mathbf{1}$ that is unstable.

*Proof.* Suppose that the conditions in Proposition 5.3 are fulfilled. Then, it follows that there exists boundary equilibria $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$, and that both are locally exponentially stable; see statements iii) and iv) in Proposition 4.1. Therefore, since we know that system (2.2) is monotone, from [40, Proposition 2.9] it follows that there exists (at least) one equilibrium point of the form $(\hat{x}^1, \hat{x}^2)$ such that $\mathbf{0} \ll \hat{x}^1, \hat{x}^2 \ll \mathbf{1}$ and $\hat{x}^1 + \hat{x}^2 \ll \mathbf{1}$. The rest of the proof (i.e., showing that $(\hat{x}^1, \hat{x}^2)$ is unstable) is similar to the approach in the proof of Theorem 5.2. $\square$

**Proposition 5.3** assures us of the existence of at least one coexistence equilibrium. Recall that system (2.2) is monotone. Therefore, since Theorem 3.2 says that, for generic parameter matrices, for each of the equilibrium points the associated Jacobian is nonsingular, the conditions in Proposition 5.3 can be shown to guarantee the existence of an odd number of coexistence equilibria, each of which must be unstable. The proof follows from a Brouwer degree argument; see proof of [39, Proposition 3.7]. In fact, for classical bivirus SIS systems, for the same stability configuration as in Theorem 5.2 and Proposition 5.3, a lower bound on the number of coexistence equilibria has been obtained in [1, Corollary 3.10, statement 2].

We now present a necessary condition for the existence of a coexistence equilibrium that involves a comparison (in the sense of matrix inequalities) between the spread parameters of the two viruses. We define for $k = 1, 2$ the following matrices: $R^k := [B_1^k \mathbf{1} \ B_2^k \mathbf{1} \ ... \ B_n^k \mathbf{1}]^\top$ We have the following proposition.

**Proposition 5.4.** Consider system (2.2) under Assumptions 1 and 2. Suppose
that, for all $k \in [2]$, $\rho(\beta_1^k(D^k)^{-1}A^k) > 1$. Suppose that $D^k = I$ for $k = 1,2$. If $\beta_1^1A^1 + \beta_2^1R^1 < \beta_1^2A^2$, then there exists no coexistence equilibrium.

*Proof.* Suppose that $\beta_1^1A^1 + \beta_2^1R^1 < \beta_1^2A^2$. Assume by way of contradiction that there exists a coexisting equilibrium, $(\hat{x}^1, \hat{x}^2)$ such that $0 \ll (\hat{x}^1, \hat{x}^2) \ll 1$, and $\hat{x}^1 + \hat{x}^2 \ll 1$. Hence, since by assumption $D^k = I$ for $k = 1,2$, by considering the equilibrium version of equation (2.2), we have the following:

$$ (5.2) \qquad -\hat{x}^1 + \beta_1^1(I - \hat{X}^1 - \hat{X}^2)A^1\hat{x}^1 + \beta_2^1(I - \hat{X}^1 - \hat{X}^2) \begin{bmatrix} (\hat{x}^1)^\top B_1^1 \hat{x}^1 \\ (\hat{x}^1)^\top B_2^1 \hat{x}^1 \\ \vdots \\ (\hat{x}^1)^\top B_n^1 \hat{x}^1 \end{bmatrix} = 0 $$

$$ (5.3) \qquad -\hat{x}^2 + \beta_1^2(I - \hat{X}^1 - \hat{X}^2)A^2\hat{x}^2 + \beta_2^2(I - \hat{X}^1 - \hat{X}^2) \begin{bmatrix} (\hat{x}^2)^\top B_1^2 \hat{x}^2 \\ (\hat{x}^2)^\top B_2^2 \hat{x}^2 \\ \vdots \\ (\hat{x}^2)^\top B_n^2 \hat{x}^2 \end{bmatrix} = 0 $$

From (5.2) and (5.3), we have the following:

$$ (5.4) \qquad (\beta_1^1(I - \hat{X}^1 - \hat{X}^2)A^1 + \beta_2^1(I - \hat{X}^1 - \hat{X}^2)) \begin{bmatrix} (\hat{x}^1)^\top B_1^1 \\ (\hat{x}^1)^\top B_2^1 \\ \vdots \\ (\hat{x}^1)^\top B_n^1 \end{bmatrix} \hat{x}^1 = \hat{x}^1 $$

$$ (5.5) \qquad (\beta_1^2(I - \hat{X}^1 - \hat{X}^2)A^2 + \beta_2^2(I - \hat{X}^1 - \hat{X}^2)) \begin{bmatrix} (\hat{x}^2)^\top B_1^2 \\ (\hat{x}^2)^\top B_2^2 \\ \vdots \\ (\hat{x}^2)^\top B_n^2 \end{bmatrix} \hat{x}^2 = \hat{x}^2 $$

By Assumption 1, we know that the matrices $A^1$ and $A^2$ are nonnegative and $\beta_1^1, \beta_2^1, \beta_1^2$ and $\beta_2^2$ are all strictly positive, and $(I - \hat{X}^1 - \hat{X}^2)$ is a positive diagonal matrix (because by Lemma 2.6, $\hat{x}^1 + \hat{x}^2 \ll 1$), whereas, because of Assumption 2, from Corollary 2.5

we have that the matrix $\beta_1^k(I - \hat{X}^1 - \hat{X}^2)A^k + \beta_2^k(I - \hat{X}^1 - \hat{X}^2) \begin{bmatrix} (\hat{x}^k)^\top B_1^k \\ (\hat{x}^k)^\top B_2^k \\ \vdots \\ (\hat{x}^k)^\top B_n^k \end{bmatrix}$ is irreducible.

Since $\hat{x}^1, \hat{x}^2 \gg 0$, from Lemma 2.8, it follows that

$$ (5.6) \qquad \rho(\beta_1^1(I - \hat{X}^1 - \hat{X}^2)A^1 + \beta_2^1(I - \hat{X}^1 - \hat{X}^2)) \begin{bmatrix} (\hat{x}^1)^\top B_1^1 \\ (\hat{x}^1)^\top B_2^1 \\ \vdots \\ (\hat{x}^1)^\top B_n^1 \end{bmatrix} = 1 $$

$$ (5.7) \qquad \rho(\beta_1^2(I - \hat{X}^1 - \hat{X}^2)A^2 + \beta_2^2(I - \hat{X}^1 - \hat{X}^2)) \begin{bmatrix} (\hat{x}^2)^\top B_1^2 \\ (\hat{x}^2)^\top B_2^2 \\ \vdots \\ (\hat{x}^2)^\top B_n^2 \end{bmatrix} = 1 $$

By assumption, $\beta_1^1A^1 + \beta_2^1R^1 < \beta_1^2A^2$. Since $\hat{x}^1 \ll 1$, it must be that $\beta_1^1A^1 +

$$ \beta_2^1 \begin{bmatrix} (\hat{x}^1)^\top B_1^1 \\ (\hat{x}^1)^\top B_2^1 \\ \vdots \\ (\hat{x}^1)^\top B_n^1 \end{bmatrix} < \beta_1^1A^1 + \beta_2^1R^1 < \beta_1^2A^2. \text{ Since } (I - \hat{X}^1 - \hat{X}^2) \text{ is positive diagonal, we} $$

obtain the following:

$$ (5.8) \qquad \beta_1^1(I - \hat{X}^1 - \hat{X}^2)A^1 + \beta_2^1(I - \hat{X}^1 - \hat{X}^2) \begin{bmatrix} (\hat{x}^1)^\top B_1^1 \\ (\hat{x}^1)^\top B_2^1 \\ \vdots \\ (\hat{x}^1)^\top B_n^1 \end{bmatrix} < \beta_1^2(I - \hat{X}^1 - \hat{X}^2)A^2 $$
Applying the condition in statement iv) of Lemma 2.8 to (5.8) yields:

$$
(5.9) \quad \rho(\beta_1^1(I - \hat{X}^1 - \hat{X}^2)A^1 + \beta_2^1(I - \hat{X}^1 - \hat{X}^2)) \begin{bmatrix} (\hat{x}^1)^\top B_1^1 \\ (\hat{x}^1)^\top B_2^1 \\ \vdots \\ (\hat{x}^1)^\top B_n^1 \end{bmatrix} < \rho(\beta_1^2(I - \hat{X}^1 - \hat{X}^2)A^2)
$$

From (5.6) and (5.9), $\rho(\beta_1^2(I - \hat{X}^1 - \hat{X}^2)A^2) > 1$. Therefore, since $\beta_2^2 > 0$, $(I - \hat{X}^1 - \hat{X}^2)$

is positive diagonal and the matrix
$$
\begin{bmatrix}
(\hat{x}^2)^\top B_1^2 \\
(\hat{x}^2)^\top B_2^2 \\
\vdots \\
(\hat{x}^2)^\top B_n^2
\end{bmatrix}
$$
is nonnegative, it is clear that

$$
(5.10) \quad \rho(\beta_1^2(I - \hat{X}^1 - \hat{X}^2)A^2 + \beta_2^2(I - \hat{X}^1 - \hat{X}^2)) \begin{bmatrix} (\hat{x}^2)^\top B_1^2 \\ (\hat{x}^2)^\top B_2^2 \\ \vdots \\ (\hat{x}^2)^\top B_n^2 \end{bmatrix} > 1,
$$

which contradicts (5.7). Therefore, there does not exist a coexistence equilibrium ($\hat{x}^1, \hat{x}^2$), where $0 \ll (\hat{x}^1, \hat{x}^2) \ll 1$ and $\hat{x}^1 + \hat{x}^2 \ll 1$; thus completing the proof. □

Proposition 5.4 answers question v) raised in Section 2.3. It can also be shown that the condition in Proposition 5.4 implies that a) the DFE is unstable, and b) the boundary equilibrium, ($\bar{x}^1, \mathbf{0}$) is unstable.

**6. Numerical Examples.** In this section, we present a series of simulations of the bivirus system with HOIs to highlight interesting phenomena that can emerge when HOIs are incorporated. The network has $n = 4$ nodes, and we set $D^1 = D^2 = I$. The pairwise interactions are captured by the following matrices

$$
(6.1) \qquad A^1 = \begin{bmatrix} 2 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}, \quad A^2 = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 \\ 1 & 1 & 2 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}.
$$

The HOI are captured by the following nonzero $B_i^k$ matrices (all other $B_i^k = 0$):

$$
B_1^1 = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{bmatrix}, B_4^1 = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}, B_1^2 = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}, B_4^2 = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}.
$$

First, note that $A^k$ is not irreducible for $k = 1, 2$, and so the graph associated with just the pairwise interactions is not strongly connected (Assumption 3 is not satisfied). Nonetheless, the hypergraph associated with $A^k$ and $B_i^k$, $i = 1, 2, \dots, n$ is strongly connected for $k = 1, 2$ (Assumption 2 is satisfied).

In our simulations, we randomly sample $x_i^k(0)$ from a uniform distribution $(0, 1)$, and then normalize the vectors $x^1(0)$ and $x^2(0)$ to ensure that $(x^1(0), x^2(0)) \in \text{int}(\mathcal{D})$. The $\beta_i^k$ are varied to yield different stability properties for the system in (2.2).

**Example 1:** We set $\beta_1^1 = \beta_2^1 = 0.2$ and $\beta_1^2 = \beta_2^2 = 5$. This ensures the inequalities of both conditions for Proposition 4.1 are satisfied. For initial conditions close to the DFE (Fig. 1a), the trajectories converge to the locally exponentially stable DFE, $(x^1 = \mathbf{0}, x^2 = \mathbf{0})$. When the initial conditions are further in the interior of $\mathcal{D}$, and depending on the particular initial condition, we observe convergence to one of two boundary equilibrium (both of which are locally exponentially stable), $(\bar{x}^1, \mathbf{0})$ or $(\mathbf{0}, \bar{x}^2)$ for some positive $\bar{x}^1 > 0.5 \times \mathbf{1}$ and $\bar{x}^2 > 0.5 \times \mathbf{1}$, see Figs. 1b and 1c, respectively. It is important to note that, without HOIs, it is impossible for a bivirus system to have the DFE, $(\bar{x}^1, \mathbf{0})$, and $(\mathbf{0}, \bar{x}^2)$ all locally exponentially stable [44, Section E].
**Example 2:** We set $\beta_1^1 = \beta_1^2 = 2$, and $\beta_2^1 = \beta_2^2 = 0.5$. The resulting bivirus system satisfies the first two inequality conditions in Proposition 5.1. The system also has two boundary equilibria $(\bar{x}^1, \mathbf{0})$ and $(\mathbf{0}, \bar{x}^2)$, which satisfy the last two inequality conditions in Proposition 5.1. Proposition 4.3 indicates that the two boundary equilibria are unstable, and this can indeed be verified. As noted below Proposition 4.3, there is no theoretical proof that $(\bar{x}^1, \mathbf{0})$ (resp. $(\mathbf{0}, \bar{x}^2)$) is the unique boundary equilibrium for virus 1 (resp. virus 2) but simulations suggest the uniqueness property holds true. Figure 1d indicates that the system converges to a coexistence equilibrium, $(\hat{x}^1, \hat{x}^2)$ for some $\hat{x}^1 \gg \mathbf{0}$ and $\hat{x}^2 \gg \mathbf{0}$; the presence of such a coexistence equilibrium is predicted by Proposition 5.1 (though the Proposition is silent on the stability properties of any coexistence equilibrium). Additional simulations suggest that this coexistence equilibrium is the unique attractor of all initial conditions in the interior of $\mathcal{D}$.

Next, we increase the HOI interaction strength for virus 2, by changing $\beta_2^2 = 0.7$. Fig. 1e and 1f show convergence to different equilibria for different initial conditions (i.e., there are multiple locally exponentially stable equilibria). Importantly, of the two stable equilibria observed, one is a boundary equilibria $(\mathbf{0}, \bar{x}^2)$, where virus 2 is endemic and virus 1 extinct, while the other is a coexistence equilibrium. Note that the inequality $s(-D^2 + \beta_1^2(I - \bar{X}^1)A^2) > 0$ of Proposition 4.3 is violated for the boundary equilibrium $(\mathbf{0}, \bar{x}^2)$, but as noted below the proposition, this is not enough to theoretically establish the local exponential stability of the boundary equilibrium $(\mathbf{0}, \bar{x}^2)$; this contrasts with the bivirus system without HOIs [45].

**Example 3:** We set $\beta_1^1 = 2.5$, $\beta_1^2 = 2$, $\beta_2^1 = 2$ and $\beta_2^2 = 2.5$. The inequalities of Theorem 5.2 are satisfied. Figs. 2a and 2b reveal that there are two locally exponentially stable boundary equilibria $(\bar{x}^1, \mathbf{0})$ or $(\mathbf{0}, \bar{x}^2)$, with convergence to either depending on the initial conditions. However, the DFE is unstable, and no trajectories in $\mathcal{D}$ converge to the DFE.

We conclude with several remarks about the simulations. First, the majority of the propositions and theorems in Sections 4 and 5 are stated under Assumption 3, and yet in the simulations, we find that the same conclusions regarding stability and existence of various equilibria continue to hold under the more relaxed Assumption 2. This suggests that an important future direction is to formally prove the results in Sections 4 and 5 will hold under Assumption 2. Second, we have a number of example systems that demonstrate multi-stability, viz., having multiple attractive (locally exponentially stable) equilibria with different regions of attraction. In the classical bivirus systems *without* HOIs, multi-stable systems are more difficult to come by. They are often constructed by special design of the two-layer network to create strong asymmetry in the spread of virus 1 and virus 2 [1, 46]. Here, all three examples appear to suggest that the presence of HOIs helps to facilitate a multi-stability property. This suggests a second line of future work, viz., further analysis of how HOIs can lead to novel dynamical phenomena in SIS models.

**7. Conclusion.** In this paper, we analyzed a competitive networked bivirus SIS HOI model. We obtained general results on equilibria and convergence under milder connectivity conditions compared to the literature. We then derived new conditions for the stability and existence (or otherwise) of various equilibria. Several lines of research are worthy of further attention, as we have noted throughout the paper. We mention a few additional ones as well: First, the paper has not commented on conditions for global stability of single-virus endemic equilibria. Hence, one line of investigation may center around that. Second, provision of conditions for uniqueness and local stability of the coexistence equilibria discussed in Section 5 remains an open
<figure><img src="image_2.png" /><figcaption>FIG. 1. Trajectories of the system (2.2), for different simulation parameters and initial conditions. 1a, 1b and 1c correspond to different initial conditions for Example 1. 1d, 1e, and 1f correspond to Example 2, for different HOI interaction strengths for virus 2, and different initial conditions.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Trajectories of the system (2.2), for different simulation parameters and initial conditions. 2a and 2b correspond to different initial conditions for Example 3.</figcaption></figure>
problem. Finally, devising control strategies for mitigating the spread of competing viruses with(out) accounting for HOI is another promising research direction.

REFERENCES

[1] B. D. ANDERSON AND M. YE, *Equilibria Analysis of a networked bivirus epidemic model using Poincaré-Hopf and Manifold Theory*, SIAM Journal on Applied Dynamical Systems, 22 (2023), pp. 2856–2889.

[2] S. ARAL AND D. WALKER, *Creating social contagion through viral product design: A randomized trial of peer influence in networks*, Management science, 57 (2011), pp. 1623–1639.

[3] G. AUSIELLO AND L. LAURA, *Directed hypergraphs: Introduction and fundamental algorithms—a survey*, Theoretical Computer Science, 658 (2017), pp. 293–306.

[4] C. BICK, E. GROSS, H. A. HARRINGTON, AND M. T. SCHAUB, *What are higher-order networks?*, SIAM Review, 65 (2023), pp. 686–731.

[5] D. E. BLOOM, D. CADARETTE, AND J. SEVILLA, *Epidemics and economics*, Finance & Development, 55 (2018).

[6] C. CASTILLO-CHAVEZ, H. W. HETHCOTE, V. ANDREASEN, S. A. LEVIN, AND W. M. LIU, *Epidemiological models with age structure, proportionate mixing, and cross-immunity*, Journal of Mathematical Biology, 27 (1989), pp. 233–258.

[7] C. CASTILLO-CHAVEZ, W. HUANG, AND J. LI, *Competitive exclusion in gonorrhea models and other sexually transmitted diseases*, SIAM Journal on Applied Mathematics, 56 (1996), pp. 494–508.

[8] K. F. CHAN, L. A. CAROLAN, D. KORENKOV, J. DRUCE, J. MCCAW, P. C. READING, I. G. BARR, AND K. L. LAURIE, *Investigating viral interference between influenza A virus and human respiratory syncytial virus in a ferret model of infection*, Journal of Infectious Diseases, 218 (2018), pp. 406–417.

[9] P. CISNEROS-VELARDE AND F. BULLO, *Multigroup SIS epidemics with simplicial and higher order interactions*, IEEE Transactions on Control of Network Systems, 9 (2021), pp. 695–705.

[10] S. CUI, F. LIU, H. JARDÓN-KOJAKHMETOV, AND M. CAO, *General SIS diffusion process with indirect spreading pathways on a hypergraph*, arXiv preprint arXiv:2306.00619, (2023).

[11] S. CUI, G. ZHANG, H. JARDÓN-KOJAKHMETOV, AND M. CAO, *On metzler positive systems on hypergraphs*, arXiv preprint arXiv:2401.03652, (2024).

[12] G. F. DE ARRUDA, G. PETRI, AND Y. MORENO, *Social contagion models on hypergraphs*, Physical Review Research, 2 (2020), p. 023032.

[13] V. DOSHI, S. MALLICK, ET AL., *Convergence of Bi-Virus Epidemic Models With Non-Linear Rates on Networks – A Monotone Dynamical Systems Approach*, IEEE/ACM Transactions on Networking, (2022).

[14] S. GRACY, B. ANDERSON, M. YE, AND C. A. URIBE, *Competitive networked bivirus SIS spread over hypergraphs*, in 2024 American Control Conference (ACC), IEEE, 2024, pp. 1–1, https://arxiv.org/pdf/2309.14230.pdf. Note: In Press.

[15] S. GRACY, P. E. PARÉ, J. LIU, H. SANDBERG, C. L. BECK, K. H. JOHANSSON, AND T. BAŞAR, *Modeling and analysis of a coupled SIS bi-virus model*, arXiv:2207.11414, (2023).

[16] S. GRACY, M. YE, B. D. ANDERSON, AND C. A. URIBE, *Towards understanding the endemic behavior of a competitive tri-virus SIS networked model*, SIAM Journal on Applied Dynamical Systems, 23 (2024), pp. 1372–1410.

[17] V. GUILLEMIN AND A. POLLACK, *Differential Topology*, vol. 370, American Mathematical Soc., 2010.

[18] A. HATCHER, *Algebraic Topology*, 2005.

[19] M. W. HIRSCH, *Systems of differential equations that are competitive or cooperative ii: Convergence almost everywhere*, SIAM Journal on Mathematical Analysis, 16 (1985), pp. 423–439.

[20] M. W. HIRSCH, *Stability and convergence in strongly monotone dynamical systems*, Journal for die reine und angewandte Mathmatik, (1988).

[21] J. HOFBAUER, *An index theorem for dissipative semiflows*, The Rocky Mountain Journal of Mathematics, 20 (1990), pp. 1017–1031.

[22] I. IACOPINI, G. PETRI, A. BARRAT, AND V. LATORA, *Simplicial models of social contagion*, Nature Communications, 10 (2019), p. 2485.

[23] J. A. JACQUEZ AND C. P. SIMON, *The stochastic SI model with recruitment and deaths I. Comparison with the closed SIS model*, Mathematical Biosciences, 117 (1993), pp. 77–125.

[24] A. JANSON, S. GRACY, P. E. PARÉ, H. SANDBERG, AND K. H. JOHANSSON, *Networked multi-virus spread with a shared resource: Analysis and mitigation strategies*, arXiv:2011.07569,
(2020).

[25] F. JAUMOTTE, M. OIKONOMOU, C. PIZZINELLI, AND M. TAVARES, *How pandemic accelerated digital transformation in advanced economies*, International Monetary Fund, 21 (2023).

[26] S. KALOGIANNIDIS, *Covid impact on small business*, International Journal of Social Science and Economics Invention, 6 (2020), pp. 387–391.

[27] W. O. KERMACK AND A. G. MCKENDRICK, *Contributions to the mathematical theory of epidemics. ii.—the problem of endemicity*, Proceedings of the Royal Society of London. Series A, containing papers of a mathematical and physical character, 138 (1932), pp. 55–83.

[28] H. KHALIL, *Nonlinear Systems*, Prentice Hall, 2002.

[29] A. LAJMANOVICH AND J. A. YORKE, *A deterministic model for gonorrhea in a nonhomogeneous population*, Mathematical Biosciences, 28 (1976), pp. 221–236.

[30] J. M. LEE, *Introduction to Smooth Manifolds*, Springer, 2013.

[31] M. Y. LI AND J. S. MULDOWNEY, *Global stability for the SEIR model in epidemiology*, Mathematical biosciences, 125 (1995), pp. 155–164.

[32] J. LIU, P. E. PARÉ, A. NEDIĆ, C. Y. TANG, C. L. BECK, AND T. BAŞAR, *Analysis and control of a continuous-time bi-virus model*, IEEE Transactions on Automatic Control, 64 (2019), pp. 4891–4906.

[33] S. A. MYERS AND J. LESKOVEC, *Clash of the contagions: Cooperation and competition in information diffusion*, in Proceedings of the 12th International Conference on Data Mining, IEEE, 2012, pp. 539–548.

[34] M. NUNO, Z. FENG, M. MARTCHEVA, AND C. CASTILLO-CHAVEZ, *Dynamics of two-strain influenza with isolation and partial cross-immunity*, SIAM Journal on Applied Mathematics, 65 (2005), pp. 964–982.

[35] J. PIRET AND G. BOIVIN, *Pandemics throughout history*, Frontiers in Microbiology, 11 (2021), p. 631736.

[36] B. A. PRAKASH, A. BEUTEL, R. ROSENFELD, AND C. FALOUTSOS, *Winner takes all: competing viruses or ideas on fair-play networks*, in Proceedings of the 21st International Conference on World Wide Web, 2012, pp. 1037–1046.

[37] A. RANTZER, *Distributed control of positive systems*, in Proceedings of the 50th IEEE Conference on Decision and Control, 2011, pp. 6608–6611.

[38] J. H. SHAPIRO AND J. H. SHAPIRO, *Nash equilibrium: Fixed points in game theory*, A Fixed-Point Farrago, (2016), pp. 51–64.

[39] H. L. SMITH, *Competing subcommunities of mutualists and a generalized kamke theorem*, SIAM Journal on Applied Mathematics, 46 (1986), pp. 856–874.

[40] H. L. SMITH, *Systems of ordinary differential equations which generate an order preserving flow. a survey of results*, SIAM Review, 30 (1988), pp. 87–113.

[41] E. D. SONTAG, *Monotone and near-monotone biochemical networks*, Systems and Synthetic Biology, 1 (2007), pp. 59–87.

[42] P. VAN MIEGHEM, J. OMIC, AND R. KOOIJ, *Virus spread in networks*, IEEE/ACM Transactions On Networking, 17 (2008), pp. 1–14.

[43] R. VARGA, *Matrix Iterative Analysis*, Springer Series in Computational Mathematics, Springer Berlin Heidelberg, 1999.

[44] M. YE AND B. D. ANDERSON, *Competitive epidemic spreading over networks*, IEEE Control Systems Letters, 7 (2022), pp. 545–552.

[45] M. YE, B. D. ANDERSON, AND J. LIU, *Convergence and equilibria analysis of a networked bivirus epidemic model*, SIAM Journal on Control and Optimization, 60 (2022), pp. S323–S346.

[46] M. YE, B. D. O. ANDERSON, A. JANSON, S. GRACY, AND K. H. JOHANSSON, *Competitive epidemic networks with multiple survival-of-the-fittest outcomes*, Systems & Control Letters, (2023), arXiv:2111.06538. Note: Under Review.

**Appendix: Proof of Theorem 3.2.** We treat first the case where the $A^k$ and $B_i^k$ are fixed. Consider (2.2) in the shorthand form $\dot{x} = f(x)$, and to emphasise the possibility of variations in the entries of $D^k$, show this explicitly by writing the equations as $\dot{x} = f_\delta(x, \delta)$, where $\delta$ is a vector containing all the diagonal entries of the $D^k$. It is easily established that

$$ (7.1) \qquad \frac{\partial f_{\delta}(x^1, x^2, D^1, D^2)}{\partial \delta} = \frac{\partial f_{\delta}(x^1, x^2, D^1, D^2)}{\partial (D^1, D^2)} = \begin{bmatrix} X^1 & 0 \\ 0 & X^2 \end{bmatrix} $$
At a coexistence equilibrium, from Lemma 2.6, it is clear that this matrix has full row rank. Consequently, the matrix $\frac{\partial f_{\delta}(x, \delta)}{\partial(x, \delta)}$ (which is obtained simply by adding further columns) also has full row rank at a coexistence equilibrium.

Similarly, in preparation for considering the second claim of the theorem, we can consider $f$ as depending explicitly on the possibly variable $A^k$ and $B^k$. It is helpful to introduce a rewriting of the key system equation (2.2). Define for $k = 1, 2$ the matrices

$$ (7.2) \qquad C^k(x^k) = \beta_1^k A^k + \beta_2^k \begin{bmatrix} (x^k)^\top B_1^k \\ (x^k)^\top B_2^k \\ \vdots \\ (x^k)^\top B_n^k \end{bmatrix} $$

This means that the system equations can be written as

$$ (7.3) \qquad \begin{aligned} \dot{x}^1 &= -D^1 x^1 + (I - X^1 - X^2) C^1(x^1) x^1 \\ \dot{x}^2 &= -D^2 x^2 + (I - X^1 - X^1) C^2(x^2) x^2 \end{aligned} $$

There are a finite number of patterns of zero entries for the $A^k$ and $B^k$ for which the hypergraph is strongly connected. When the hypergraph is strongly connected, by Corollary 2.5 for any positive $x^k$ the matrices $C^k$ are irreducible. For each matrix separately, arbitrarily choose any one such pattern, and let $\gamma^k$ denote the vector of nonzero entries of the $C^k$; also set $\gamma = [(\gamma^1)^\top (\gamma^2)^\top]^\top$. Denote the associated $f$ as $f_\gamma(x, \gamma)$. Then, one can show that

$$ (7.4) \qquad \frac{\partial f_\gamma(x, \gamma)}{\partial \gamma} = \begin{bmatrix} I - X^1 - X^2 & 0 \\ 0 & I - X^1 - X^2 \end{bmatrix} \begin{bmatrix} \frac{\partial(C^1(x^1)x^1)}{\partial \gamma^1} & 0 \\ 0 & \frac{\partial(C^2(x^2)x^2)}{\partial \gamma^2} \end{bmatrix} $$

A detailed calculation as set out in [1] shows that the second matrix in the product (on the right-hand side of (7.4)) has full row rank. (The irreducibility of each $C^k(x^k)$ for any positive $x^k$ ensures in particular that every row is nonzero, and this is the key property which guarantees the conclusion.) We also know from Lemma 2.6 that no equilibrium of (2.2) can cause any diagonal entry of $I - X^1 - X^2$ to be zero, which means that the first matrix of the product is nonsingular. Hence, the matrix on the left has full row rank at any equilibrium. The matrix $\frac{\partial f_\gamma(x, \gamma)}{\partial(x, \gamma)}$, obtained by the addition of further columns to $\frac{\partial f_\gamma(x, \gamma)}{\partial \gamma}$, then has full row rank.

We now complete the proof for the first part of the theorem. (Proof of the second part is virtually identical). Without loss of generality, suppose that for some $\bar{\delta}$, all diagonal entries $\delta_i$ of the $D^k$ satisfy $0 < \bar{\delta} < \delta_i < \bar{\delta}^{-1}$. Denote the set of such $\delta$ as $\mathcal{D}$. Call the set of $x$ in the interior of the region of interest $\mathcal{X}$. Then $\mathcal{X} \times \mathcal{D}$ is a manifold, being a product of open sets, and $f_\delta(x, \delta)$ is a mapping of the manifold $\mathcal{X} \times \mathcal{D}$. Let the image be denoted by $\mathcal{Y}$, and let $\mathcal{Z}$ be the submanifold of $\mathcal{Y}$ consisting of the single point $\mathbf{0}_{2n}$. (The set $f_\delta^{-1}(\mathcal{Z})$ is the set of zeros of $f$.) Since the Jacobian of $f_\delta$ with respect to $(x, \delta)$ has full row rank on this set, this means that the map $f_\delta$ is transversal to $\mathcal{Z}$. By the Parametric Transversality Theorem, this means that for all choices of $\delta = \delta^*$ excluding a set of measure zero, the Jacobian $\frac{\partial f_{\delta^*}(x, \delta^*)}{\partial x}$ will be nonsingular at a preimage of $\mathcal{Z}$, i.e at a zero. This is equivalent to saying that the zeros are nondegenerate. The boundedness of the region $\mathcal{X}$ then guarantees that there can only be a finite number of such zeros.