hypothesisHypothesis claimClaim General SIS on HypergraphsS. Cui, F. Liu, H.
Jardón-Kojakhmetov, and M. Cao ex_supplement

# General SIS diffusion process with indirect spreading pathways on a
hypergraph

Shaoxuan Cui Bernoulli Institute, Faculty of Science and Engineering,
University of Groningen, The Netherlands<EMAIL_ADDRESS>Fangzhou Liu
Research Institute of Intelligent Control and Systems, School of Astronautics,
Harbin Institute of Technology, China<EMAIL_ADDRESS>Hildeberto
Jardón-Kojakhmetov Bernoulli Institute, Faculty of Science and Engineering,
University of Groningen, The Netherlands<EMAIL_ADDRESS>Ming
Cao ENTEG, Faculty of Science and Engineering, University of Groningen, The
Netherlands<EMAIL_ADDRESS>
###### Abstract

Conventional graphs captures pairwise interactions; by contrast, higher-order
networks (hypergraphs, simplicial complexes) are introduced to the
interactions involving more parties, which has been rapidly applied to
characterize a growing number of complex real systems. It becomes apparent
that higher-order interactions give rise to new challenges for the rigorous
analysis of dynamics on a higher-order network. In this paper, we study a
series of Susceptible-infected-susceptible-type (SIS-type) diffusion processes
with both indirect and direct pathways on a directed hypergraph. We start with
a concrete setup where we choose a polynomial interaction function to describe
how several agents influence one another over a hyperedge. Then, under this
interaction function, we further extend the system and propose a bi-virus
competing model on a directed hypergraph by coupling two single-virus models
together. Finally, the most general model in this paper considers a general
interaction function for both the single-virus and bi-virus cases. For the
single-virus model, we provide a comprehensive characterization of the healthy
state and endemic equilibrium. For the bi-virus setting, we further give the
analysis of the existence and stability of the healthy state, dominant endemic
equilibria, and coexisting equilibria. All theoretical results are supported
additionally by some numerical examples.

###### keywords:

Susceptible-infected-susceptible (SIS) model, Vector transmission, Networked
systems, Hypergraphs, Dynamical systems

05C65, 34D05, 34C12, 37N25, 92D30

##  1 Introduction

Mathematical modeling of epidemic spreading processes is crucial in
epidemiology and has a long history of scientific study. Especially during the
recent Covid-19 period, epidemic models have demonstrated even greater
importance [59, 54]. The precise and timely prediction [1], parameter
identification [51, 63] and design of control strategies [24, 4, 25, 65] using
epidemic models have provided us with strong guidance to fight against the
pandemic. Besides the scenario of the spreading of a virus, epidemic models
also show great potential to describe very general diffusion processes,
including information diffusion in social networks [27, 67, 29], computer
virus dissemination [34, 21, 62] and others [17, 48].

Among various epidemic models, compartmental models are especially popular
[38, 68, 40]. In such a framework, each compartment indicates a state, and
transmissions between compartments show changes between the states of
infections. For example, usually, $S$ denotes the susceptible, $I$ is the
infected, and $R$ is the recovered states. By considering various combinations
of state transitions, epidemic models, such as SIS, SIR, and SIRS, have been
proposed. In the early stage of the development of compartmental epidemic
models, a relatively simple modeling framework, called _scalar models_ , was
proposed and studied. In this framework, a single population is studied and is
assumed to be well-mixed, i.e., everyone within the population has the same
chance to get infected. This simple model is detailed in, e.g., the survey
[15].

It is evident that such a simple framework with its assumption on homogeneity
provides only a relatively rough approximation of the real complex diffusion
process and fails to capture the heterogeneity in the real world (for example,
age structures, spatial diversity, and social behaviors). In order to overcome
this drawback, the population can be divided into several sub-populations, and
a networked structure (hence referred to as _network models_) is further
introduced to describe the interactions among the sub-populations. For
instance, the contributions to the networked SIS models include those in
continuous-time [58, 57, 23] and others in discrete-time [42, 46, 30]. By
further considering that two or more viruses spread over the network, the bi-
or multi-virus networked SIS models have also been considered, for example in
[64, 32, 14, 43, 44]. All these above-mentioned models assume that the
interaction between sub-populations is linear and thus additive. Such an
assumption may be violated under some situations since the infection process
sometimes may have the saturation effect and thus slow down [33, 60, 66, 49].
With the introduction of non-linear interaction functions, the networked bi-
virus SIS model [61] is able to capture the saturation effect of the infection
process.

_Networked epidemic models_ still use the conventional graphs, one edge of
which just connects one node to another. This leads to the restriction that
dynamical systems on conventional graphs only take pairwise interactions into
account. In contrast, it can be frequently observed that people not only
interact with each other one one one but also as a group. Recently, it has
been shown that higher-order networks, including _hypergraphs_ and _simplicial
complexes_ , provide a powerful modeling tool to characterize effects beyond
pairwise interactions [2]. In this higher-order network structure, one
hyperedge can contain more than two nodes and it is thus suitable to model
group-wise multi-body interactions. Along this line, an SIS model on a
simplicial complex is proposed and further studied via a mean-field
approximation in [16]. In each simplex, the healthy agent can only get
infected via the hyperedge when all other agents in the hyperedge are
infected. Followed by this modeling framework and by further modifications, a
similar SIS model on a hypergraph is proposed and studied in [7]. By
calculating the probabilities of infection and using epidemic link equations
[36], another type of SIS models on a simplicial complex is studied in [37].
In [9], an SIS model on a hypergraph is proposed, and most analytical results
are achieved via adopting a special network structure (hyperstar or
homogeneous hypergraph). Furthermore, the model they studied is a threshold
model and one agent may get infected through the hyperedge when more than a
certain amount of the agents in the hyperedge are infected. By considering
that one may react differently in different surroundings, an appropriate model
is further proposed and studied in [53]. In [28], a bi-virus SIS model on a
simplicial complex has been studied. However, under the mean-field
approximation and homogeneous assumption, the model in [28] is approximated by
a planar system, and the coexistence results are established purely via
simulations, which restricts the model from wider applications.

Apart from discrepancies in network topology, all the aforementioned models
only describe direct human-to-human spreading processes. However, in a real
epidemic process, especially those concerning waterborne or insect-borne
viruses, the _pathogen_ also spreads through different kinds of media.
Analogously, for information diffusion on social networks, a post on Twitter
or a video on Tiktok both facilitates information spreading. For this reason,
it is also of great significance to capture the indirect spreading process
through different types of media. This is usually realized by coupling the
aforementioned epidemic models with extra pathogen dynamics on resources (or
media). With this idea, the SIWS111Here, the $W$ stands for “water” but it is
only indicative of the medium, and is by no means restricted to it.222In this
paper, SIWS means susceptible-infected-water-susceptible, and the use of $W$
is different from waning immunity, which has also been studied in [6, 26, 20].
model on a graph, which is an SIS model coupled with pathogen dynamics, is
proposed and studied by, for example, [13, 19, 18, 41, 31, 8]. However, the
analysis of SIWS models on higher-order networks is, although extremely
relevant, still missing.

The main purpose of this paper is to propose and study a general framework for
spreading processes on hypergraphs. It is then very natural to think about an
SIWS model on a higher-order network potentially with non-linear and non-
polynomial interaction functions.

Briefly, we list the contributions and organization of our paper as follows.
Firstly, in section 3, making use of a specific choice of interaction
function, we propose a single-virus SIWS model on a directed hypergraph, which
incorporate both a direct human-to-human pathway and an indirect human-
resource-human pathway. For the modeling, we introduce interaction functions,
the spreading mechanism of different orders, and higher-order network
structures (hypergraphs, simplicial complex).

In section 4, we give a full analysis of the model proposed in section 3. We
give the results regarding the healthy state and the endemic equilibrium.
Especially, we find out that bistability may occur due to the introduction of
higher-order interactions related to the hypergraph structure. We also confirm
that the proposed model is an irreducible monotone system and presents no
chaos.

In section 5, we further consider a bi-virus competing scenario and provide
abundant analytical results. Particularly, we give the results regarding the
healthy state, dominant endemic equilibrium, and coexisting equilibrium. The
bi-virus system is also an irreducible monotone system and presents no chaos.

In section 6, to further generalize our framework by considering an abstract
interaction function and prove similar results, for both the single-virus and
the bi-virus cases. This general system remains irreducible monotone.

Finally, in section 7, some numerical examples are given to highlight the
theoretical contributions.

##  2 Preliminaries

A tensor $T\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{k}}$ is a
multidimensional array, where the order is the number of its dimension $k$ and
each dimension $n_{i}$, $i=1,\cdots,k$ is a mode of the tensor. A tensor is
cubical if every mode has the same size, that is $T\in\mathbb{R}^{n\times
n\times\cdots\times n}$. We further write $k$-th order cubical tensor as
$T\in\mathbb{R}^{n\times n\times\cdots\times n}=\mathbb{R}^{[n,k]}$. A cubical
tensor $T$ is called supersymmetric if $T_{j_{1}j_{2}\ldots j_{k}}$ is
invariant under any permutation of the indices.

The definition of a hypergraph that we consider in this paper follows from
[11]. Furthermore, we introduce a set of several tensors to collect the
information of the weights. A weighted and directed hypergraph is a triplet
$\mathbf{H}=(\mathcal{V},\mathcal{E},A)$. The set $\mathcal{V}$ denotes a set
of vertices and $\mathcal{E}=\\{E_{1},E_{2},\cdots\\}$ is the set of
hyperedges. A hyperedge is an ordered pair $E=(\mathcal{X},\mathcal{Y})$ of
disjoint subsets of vertices; $\mathcal{X}$ is the tail of $E$ while
$\mathcal{Y}$ is the head. According to the physical meaning of the scenario
of a diffusion process, we assume that each hyperedge has only one tail but
one or multiple ($\geq 1$) heads. The set of tensors
$A=\\{A_{2},A_{3},\cdots\\}$ collect the weights of all hyperedges, where
$A_{2}=[A_{ij}]$ collects the weights of all second-order hyperedges,
$A_{3}=[A_{ijk}]$ collects the weights of all third-order hyperedges, and so
on. For instance, $A_{ijkl}$ denotes the weight of the hyperedge where $i$ is
the tail and $j,k,l$ are the heads. For simplicity, in this paper, we also use
the weight (for example, $A_{\bullet}$) to denote the corresponding hyperedge.
If all hyperedges only have one tail and one head, then the network is a
standard directed and weighted graph. As a special case, a weighted and
undirected hypergraph is a triplet $\mathbf{H}=(\mathcal{V},\mathcal{E},A)$,
where $\mathcal{E}$ is now a finite collection of non-empty subsets of
$\mathcal{V}$ [2]. An undirected hypergraph can be represented as a
supersymmetric tensor. For simplicity, we can use the notation $\pi(\bullet)$
, where $\pi(\bullet)$ denotes any permutation of $\bullet$, to indicate the
undirected hyperedge. For example, an undirected hyperedge can be written as
$A_{\pi(ijk)}$ since $A_{ijk},A_{ikj},A_{kij},A_{kji},A_{jik},A_{jki}$ take
the same value in an undirected graph. An undirected simplicial complex is a
special type of undirected hypergraph that contains all nonempty subsets of
hyperedges as hyperedges. If $A_{\pi(ijk)}$ is a hyperedge of the simplicial
complex then $A_{\pi(ij)},A_{\pi(ik)},A_{\pi(jk)}$ are also hyperedges of the
simplicial complex. In order to capture the heterogeneity of the spreading
process, we extend the definition of the simplicial complex into a directed
simplicial complex. That is, a simplicial complex is a special type of
directed hypergraph such that for all
$E=(\mathcal{X},\mathcal{Y})\in\mathcal{E}$ and
$\tilde{E}=(\mathcal{\tilde{X}},\mathcal{\tilde{Y}}),\mathcal{\tilde{X}}\subset\mathcal{{X}},\mathcal{\tilde{Y}}\subset\mathcal{{Y}}$,
then $\tilde{E}\in\mathcal{E}$; For example, if $A_{ijk}$ is a hyperedge of
the directed simplicial complex then $A_{ij},A_{ik}$ are also hyperedges of
the simplicial complex.

_Notation:_ $\mathbb{R}$ and $\mathbb{N}$ denote the set of real numbers and
nonnegative integers, respectively. Given a square matrix
$M\in\mathbb{R}^{n\times n}$, $\rho(M)$ is the spectral radius of $M$, which
is the largest absolute value of the eigenvalues of $M$. By $s(M)$, we refer
to the the largest real part among the eigenvalues of $M$. For a matrix
$M\in\mathbb{R}^{n\times r}$ and a vector $a\in\mathbb{R}^{n}$, $M_{ij}$ and
$a_{i}$ denote the element in the $i$th row and $j$th column and the $i$th
entry, respectively. For any two vectors $a,b\in\mathbb{R}^{n}$, $a\gg(\ll)b$
represents that $a_{i}>(<)b_{i}$, for all $i=1,\ldots,n$; $a>(<)b$ means that
$a_{i}\geq(\leq)b_{i}$, for all $i=1,\ldots,n$ and $a\neq b$; and
$a\geq(\leq)b$ means that $a_{i}\geq(\leq)b_{i}$, for all $i=1,\ldots,n$ or
$a=b$. These component-wise comparisons are also applicable to matrices with
the same dimension. The vector $\mathbf{1}$ ($\mathbf{0}$) represents the
column vector or matrix of all ones (zeros) with appropriate dimensions. The
matrix $I$ represents the identity matrix with the appropriate dimension.

##  3 Model Set-up

In this section, we introduce our model set-up and some important concepts
(spreading mechanisms, interaction functions, indirect spreading and so on).

###  3.1 General spreading mechanism and SIS model set-up

We first consider the directed simplicial complex structure and explain the
physical meaning of this network structure. The directed hypergraph structure
allows us to study more general dynamics compared with the undirected one. We
adopt the following setting based on the physical interpretation of a
spreading process. Firstly, we recall that there is only one tail agent and
the other agents are all head agents in each hyperedge. The hyperedge
describes all head agents’ joint influence on the tail agent as all agents
gather in the group. For example, a weight of a hyperedge of a $2$-simplex is
denoted as $A_{ijk}$ and is interpreted as the joint influence of $j$ and $k$
on $i$. Here, we always set the first index in the subscript as the tail agent
and all other are heads. We also assume that if the agents appear in a
hyperedge, then the hyperedge with other combinations of these agents must
also exist. For example, if $A_{ijk}$ exists, then $A_{jik}$, $A_{kij}$,
$A_{ikj}$, $A_{jki}$, $A_{kji}$ must also exist. Since usually $A_{ijk}$ and
$A_{ikj}$ represent the same influence, we set them equal. Note also that the
physical meaning of $A_{jik}$ and $A_{ikj}$ are usually different, they will
take usually different values. Figure 1 shows examples of fourth-order
hyperedges.

$i$$j$$k$$l$(a)$i$$j$$k$$l$(b) Figure 1: Graphical representations of two
fourth-order hyperedges. (a) A directed hyperedge $A_{kijl}$ in a simplicial
complex. If $A_{kijl}$ is in the simplicial complex, then the (directed)
$3$-hyperedges $A_{kjl}$, $A_{klj}$, $A_{kil}$, $A_{kli}$, $A_{kij}$ and
$A_{kji}$ are also in the simplicial complex. Moreover, $A_{kil}$ further
indicates that the directed $2$-hyperedges $A_{ki}$ and $A_{kl}$ are in the
simplicial complex (which are, of course, conventional directed edges). The
directions of $A_{jl}$, $A_{li}$, and $A_{ij}$ are not relevant for this
example. (b) A (general) directed hyperedge $A_{lijk}$. We remind the reader
that we only consider hyperedges with one tail, hence in this paper the first
index always denotes the tail.

Next, in order to specify our model, we need to know about which mechanism the
spreading process obeys. For example, a pairwise spreading mechanism is clear
from e.g. [8, 42, 58]. An agent can only get infected only if its neighbor in
the graph is infected. As for higher-order interactions, in [16], it is the
rule that if all the other (head) agents in the simplex are infected, then the
(tail) agent also has a chance to get infected. However, inspired by some
threshold models similar to [9], we consider a more general rule of the
spreading mechanism. That is, for an $n$-simplex (a hyperedge of $n$ agents),
the tail agent may get infected if and only if at least $m$ of the $n-1$ head
agents are infected, where $m\leq n-1$. We say the rule described above is an
$m$-th-order rule of an $n$-body interaction. If $m=n-1$, we say that the
interaction is a full-order $n$-body interaction.

Now, after knowing the interaction rule, we can specify the interaction
function of the spreading process. Firstly, we consider the simple case of
pairwise interactions. If there is a directed edge with the weight $A_{ij}$,
then the influence of $j$ on $i$ is described as an interaction function
$\beta_{i}A_{ij}X_{j}(t)$, where $\beta_{i}$ is the intrinsic parameter of
agent $i$ to describe how likely $i$ is infected when the spreading rule is
satisfied; the weight $A_{ij}$ describes how close the contact between $i$ and
$j$ is; $X_{j}(t)$ is a binary variable denoting whether $j$ is infected
($X_{j}=1$) or not ($X_{j}=0$). One can observe that this interaction is well-
defined because it is non-negative and upper bounded by $1$, which means that
the contribution of $j$ to the possibility of $i$ getting infected is in the
interval $[0,1]$. Furthermore, as long as $j$ is healthy, $j$ will never
infect $i$.

When it comes to an $m$-th order $n$-body interaction, we define the
interaction function as

(1) $\beta_{in}A_{ijk\dots}\frac{\Sigma^{m}_{n}}{C_{n-1}^{m}},$

where the subscript $n$ in $\beta_{in}$ denotes that the parameter is w.r.t
$n$-body interaction, $\Sigma^{m}_{n}$ is the sum of all the polynomial
combinations (ignoring the order) of $m$ binary variables of the head agents
and $C_{n-1}^{m}$ is the combination $m$ out of $n-1$. For example, if we
consider a second-order $4$-body interaction, we get the interaction function
of

(2) $\beta_{i4}A_{ijkl}\frac{X_{j}X_{k}+X_{j}X_{l}+X_{k}X_{l}}{C_{3}^{2}}.$

In the case of the full-order $(m+1)$-body interaction, the interaction
function is $\beta_{i(m+1)}A_{ijk\dots}X_{j}X_{k}\cdots$. Thus, one can use
multiple lower-order hyperedges where the interaction is full-order to replace
the higher-order hyperedge where the interaction is not full-order.

######  Lemma 3.1.

Any not-full-order ($l$-order, $l<m$) ($m+1$)-body interaction (1) can be
replaced by a composition of multiple full-order $(l+1)$-body interactions,
where $l<m$.

######  Proof 1.

Firstly, consider the equation (1). If the interaction is not in full order,
$\Sigma^{l}_{n}$ is a composition of several products among $l$ variables.
Each of these products corresponds to a full-order $(l+1)$-body interaction.
By rearranging the weights, any not-full-order ($m+1$)-body interactions (1)
can be replaced by a composition of multiple full-order $(l+1)$-body
interactions, where $l<m$.

For instance, $A_{ijkl}$ in the interaction function (2) can be replaced by a
combination of $A_{ijk}$, $A_{ikj}$, $A_{ikl}$, $A_{ilk}$, $A_{ijl}$,
$A_{ilj}$. Due to the simplicial complex structure, if $A_{ijkl}$ exists, all
the links $A_{ijk}$, $A_{ikj}$, $A_{ikl}$, $A_{ilk}$, $A_{ijl}$, $A_{ilj}$
exist too. By choosing again a proper weight, one only needs to consider full-
order interactions under this choice of the interaction function. Note that,
by doing so, the model originally based on a simplicial complex is now
constructed on a hypergraph instead, because some lower-order interaction may
be absent while the corresponding higher-order interaction may still exist
after the replacement.

The choice of interaction function is not unique. In order to make sure it is
well-defined, one only needs to ensure that

  * •

it is zero when the corresponding rule is not satisfied;

  * •

otherwise, it is positive but upper bounded by $1$.

In this paper, we firstly focus on a concrete system where we use the
polynomial form $\beta_{in}A_{ijk\dots}\frac{\Sigma}{C_{n}^{m}}$ as
interaction functions. Afterwards, we further study a more general case when
the interaction function takes an abstract form. If other types of
interactions are considered, the equivalence between one higher-order
interaction and multiple lower-order interactions may not be valid anymore,
depending on the function.

By adopting the interaction function (1) and inspired by [16], we model an
epidemic process as a continuous-time Markov chain. For simplicity, we only
consider the full-order interaction by adopting the interaction function
introduced above; however, bear in mind that it indeed contains all other
lower-order interactions and is just replaced by the full-order interaction on
lower-order hyperedges. Although our model is firstly developed on a
simplicial complex structure, after the replacement, the model is built on a
hypergraph instead. If every hyperedge originally incorporates only the full-
order interaction, then the simplicial complex structure remains unchanged.

Our model is derived by considering all these interactions (1) with different
orders. The state transitions in the SIS model can be characterized by a
Markov chain with $2^{n}$ states, where $n$ is the number of nodes in the
simplicial complex. Specifically, we denote the random variable $X_{i}(t)$ as
the state of node $i$ at time $t$. Moreover, $X_{i}=0$ if node $i$ is in state
$S$ (suceptible) and $X_{i}=1$ if node $i$ is in state $I$ (infected). The
infection process is described as the interaction above; whereas, the curing
process is assumed to be passive with the rate $\delta_{i}$. The local nodal
dynamics can be modeled as a two-state Markov chain in continuous time. For a
sufficiently small time $\Delta t$, it holds for the infection process that

(3) $\begin{split}&\textbf{Prob}[X_{i}(t+\Delta
t)=1|X_{i}(t)=0]=\sum_{j\in\mathbf{N}_{2}^{\text{i}}}\beta_{i}A_{ij}X_{j}(t)\Delta
t\\\
&+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{i3}A_{ijk}X_{j}(t)X_{k}(t)\Delta
t+\cdots+\sum_{j,k,l,\cdots\in\mathbf{N}_{n}^{\text{i}}}\beta_{in}A_{ijkl\cdots}X_{j}(t)X_{k}(t)X_{l}(t)\cdots\Delta
t\end{split}$

where the notation $\mathbf{N}_{2}^{\text{i}}$ denotes the set including all
the conventional edges where the node $i$ is the tail, the notation
$\mathbf{N}_{n}^{\text{i}}$ describes the set including all the $2$-simplexes
where the node $i$ is the tail, and the $n$-body interaction within the
$n$-simplex is of the full order.

For the curing process, it holds that

(4) $\textbf{Prob}[X_{i}(t+\Delta t)=0|X_{i}(t)=1]=\delta_{i}X_{i}(t)\Delta
t.$

We now apply the mean-field approximation, leading to

$\sum_{j\in\mathbf{N}_{2}^{\text{i}}}\beta_{i}A_{ij}X_{j}(t)=\sum_{j\in\mathbf{N}_{2}^{\text{i}}}\beta_{i}A_{ij}\mathbf{E}(X_{j}(t)),$
$\sum_{j,k\in\mathbf{N}_{3}^{\text{i2}}}\beta_{i3}A_{ijk}X_{j}(t)X_{k}(t)=\sum_{j,k\in\mathbf{N}_{3}^{\text{i2}}}\beta_{i3}A_{ijk}\mathbf{E}(X_{j}(t)X_{k}(t)).$

This approximation is valid and accurate if the states of the neighbors are
sufficiently independent and the number of the in-neighbors of node $i$ is
large.

Thus, the central limit theorem [35] applies and it follows that when $\Delta
t\rightarrow 0$, by adopting (3), one gets

(5)
$\begin{split}&\frac{d\mathbf{E}(X_{i}(t))}{dt}=\mathbf{E}\Big{(}\Big{[}\sum_{j\in\mathbf{N}_{2}^{\text{i}}}\beta_{i}A_{ij}X_{j}(t)+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{i3}A_{ijk}X_{j}(t)X_{k}(t)\\\
&+\sum_{j,k,l,\cdots\in\mathbf{N}_{n}^{\text{i}}}\beta_{in}A_{ijkl\cdots}X_{j}(t)X_{k}(t)X_{l}(t)\Big{]}(1-X_{i}(t))-\delta_{i}X_{i}(t)\Big{)}.\end{split}$

Notice that
$\mathbf{E}(X_{i}(t)X_{j}(t))=\mathbf{E}(X_{i}(t))\mathbf{E}(X_{j}(t))-\mathbf{cov}(X_{i}(t),X_{j}(t))$.
For the large-scale network, we adopt the independence assumption [29, 39] so
that
$\mathbf{cov}(X_{i}(t),X_{j}(t)))=0$. The article [van2015accuracy] indicates
that the accuracy of the approximation increases with the network scale by
using an SIS model on a conventional graph. It shows that this approximation
is generally sufficiently accurate on large-scale networks and the error is
usually tolerable on a network with several hundred nodes. In [16], they show
that a mean-field model captures quantitative features of the origin
stochastic process using a simplicial complex setting. By ignoring all
covariance and letting $x_{i}(t)=\textbf{Prob}[X_{i}(t)=1]$, we get

(6)
$\begin{split}\dot{x}_{i}=&(1-x_{i})\left(\sum_{j\in\mathbf{N}_{2}^{\text{i}}}\beta_{i}A_{ij}x_{j}\right)+(1-x_{i})\left(\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{i3}A_{ijk}x_{j}x_{k}\right)+\cdots\\\
&+(1-x_{i})\left(\sum_{j,k,l,\cdots\in\mathbf{N}_{n}^{\text{i}}}\beta_{in}A_{ijkl\cdots}x_{j}x_{k}x_{l}\cdots\right)-\delta_{i}x_{i},\end{split}$

which is similar to the general higher-order SIS model in [7] but now the
parameter of $\beta_{i},\beta_{i3},\cdots,\beta_{in}$ can be heterogeneous,
and our model is actually based on a directed hypergraph. In contrast to [7],
we now provide an explanation, at the microscopic level, of the physical
meaning of our model, the mechanism behind the modeling framework, and how it
relates to the stochastic epidemic model. Furthermore, we explore the
connection between the epidemic model on a general hypergraph and on a
simplicial complex. The spreading process on a simplicial complex may become a
process on a hypergraph when one assumes that the different spreading
mechanisms occur on different hyperedges.

In [16], the authors use a simplicial structure but with undirected and
equally-weighted hyperedges to construct the stochastic process. They only
consider a full-order interaction. If one assumes that the weights of the
hypergraph are the same in a group of the same agents and chooses the same
interaction function above (1), then our model is equivalent to the SCM
(Social Contagion Model) model proposed in [16]. It is worthwhile mentioning
that our model does consider any-order interaction, but because of the
equivalence result we establish, they are implicit in the lower-order
hyperedge. Because diverse kinds of interactions are considered, our model is
eventually equivalent to a model on a general hypergraph according to the
Lemma 3.1. In [16], it shows that the mean-field model is indeed
quantitatively comparable to the original stochastic process using simplicial
complexes with several hundred nodes.

Since achieving analytical results for the arbitrary higher-order case is
challenging and tedious, we first consider the hypergraph structure up to
hyperedges of degree 2 (three nodes involved in a hyperedge). In such a case,
the dynamics of (6) read as

(7)
$\begin{split}\dot{x}_{i}=&(1-x_{i})\left(\sum_{j\in\mathbf{N}_{2}^{\text{i}}}\beta_{i}A_{ij}x_{j}\right)+(1-x_{i})\left(\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{i3}A_{ijk}x_{j}x_{k}\right)-\delta_{i}x_{i},\end{split}$

which is similar to the simplicial SIS model in [7] but here the parameters of
$\beta_{i}$ and $\beta_{i3}$ can be heterogeneous. Just to recall, we now only
consider the hyperedge up to 3 bodies. However, the results based on a
hypergraph with even higher-order hyperedges will be still provided later when
we consider a general case in section 6.

Finally, although the framework so far describes a direct spreading process
reasonably well, it cannot characterize indirect spreading. For example, one
can get infected not only because one gets in contact with infected neighbors
but also because one uses a resource that carries the pathogen of the disease.
A common way to overcome this drawback is to introduce coupled pathogen
dynamics with a resource. In the next subsection, we will introduce how to
couple the system that we have proposed so far with a pathogen evolution.

###  3.2 Higher-order SIWS Model

Consider an SIS-type dynamics over a human contact hypergraph consisting of
$n$ groups of individuals (agents) and the pathogen diffuses over $m$
different types of resources that are utilized by some agents in the human
contact network. An individual can be infected by the disease either through
its infected neighbor in the human contact network or the polluted resource
the individual utilizes. In turn, the resource can be contaminated by the
infected agents using the resource.

By adding a $W$ (Water) compartment, we propose a higher-order SIWS model. We
emphasize that the compartment $W$ can refer to any other media.

In this setting, the epidemic process and pathogen dynamics in the resource
evolve as

(8)
$\displaystyle\dot{x}_{i}=-\delta_{i}x_{i}+(1-x_{i})\left(\sum_{j=1}^{n}\beta_{ij}x_{j}+\sum_{j=1}^{m}\beta_{ij}^{w}w_{j}\right)$
$\displaystyle\qquad+(1-x_{i})\left(\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{ijk}x_{i}x_{j}+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{ijk}^{wx}w_{j}x_{k}+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{ijk}^{xw}x_{j}w_{k}\right),$
(9)
$\displaystyle\dot{w}_{j}=-\delta_{j}^{w}w_{j}+\sum_{k=1}^{n}c_{jk}^{w}x_{k}+\sum_{i,k\in\mathbf{N}_{3}^{\text{i}}}c_{jik}^{xx}x_{i}x_{k},$

where $i=1,\ldots,n$, $j=1,\ldots,m$, and all notations and parameters are
defined in table 1. We remind the reader that for now we consider interactions
of up to 3 bodies.

Table 1: Notations and parameters in equations (8)-(9) $x_{i}=x_{i}(t)$ |  Infection probability of agent $i$
---|---
|  or infection proportion of group $i$ at time $t$
$w_{j}=w_{j}(t)$ |  Pathogen concentration in the $j$-th resource at time $t$
$\delta_{i}$ |  Healing rate of node $i$
$\delta_{j}^{w}$ |  Decaying rate of the pathogen in resource $j$
$A_{ij}$ | Weight of the conventional edges (first-order hyperedges)
|  in the corresponding $n$-node simplicial $\mathcal{G}(A)$
$A_{ijk}$ | Weight of the second-order hyperedges
|  in the corresponding $n$-node simplicial $\mathcal{G}(A)$
$\beta_{i}$ |  First-order infection rate of node $i$
$\beta^{xx}_{i}$ |  Second-order infection rate of node $i$
$\beta_{ij}:=\beta_{i}A_{ij}$ |  First-order effective transmission rate from node $j$ to node $i$
$\beta_{ijk}:=\beta_{i}^{xx}A_{ijk}$ | Second-order effective transmission rate
|  from node $j$ and $k$ to node $i$
$c_{jk}^{w}$ |  First-order person-resource contact rate of agent $j$
|  with resource $k$
$c_{ijk}^{xx}$ | Second-order person-resource contact rate
|  of agent $i$ and $k$ with resource $j$
$\beta_{ij}^{w}$ |  First-order resource-person transmission rate of agent $i$
|  with resource $j$
$\beta_{ijk}^{wx},\beta_{ijk}^{xw}$ | Second-order human-infrastructure co-infection rate
|  of node $j$ and $k$ to node $i$, the superscript $xw$ ($wx$) indicates $j$
($k$)
|  is a human index and $k$ ($j$) is a resource index

For (8)-(9), we make the assumption that no direct interaction exists between
resources because we mainly focus on the pathogen (information) diffusing over
different types of media. The inner interaction inside the same type of medium
can be ignored when we model the same type of medium as a whole. Since we also
assume that the different types of media are sufficiently independent, we do
not need to consider the higher-order interaction purely induced by resources.
In the case of the higher-order interaction where both resource and agent are
present, we assume that the higher-order interaction takes a similar form as
when only agents are present, i.e., there are no $xw$-cross-terms, and thus
the resource remains independent of each other.

Due to the fact that the indirect spreading in the resource is a macroscopic
process and the concentration cannot be simply modeled as a binary variable,
it is not suitable to model the concentration in the Markov chain. However,
the idea of coupling the epidemic dynamics in the human contact network with
the pathogen dynamics is widely adopted and supported in e.g. [41, 19, 45,
55].

It is worthwhile mentioning that the SIWS model usually describes a process of
water-borne disease or any other virus that can spread over a medium other
than the human-to-human pathway. So in the above, we also explain the
Simplicial SIWS model in this context. However, notice that a real virus
spreading usually obeys the rule that one may get infected if and only if one
of the agents in the (hyper)-edge is infected, and thus all the interactions
are of the first order (pairwise). That is to say, the conventional graph
structure is sufficient to study such a case if we assume the interaction
takes a polynomial form. Interested readers may refer to the literature about
SIWS models on the conventional graphs in discrete time [8] or in continuous
time [41, 19, 45].

In the context of _information diffusion_ , things are different. Firstly, the
higher-order interaction really counts. A real spreading process is a hybrid
process of interactions with all orders. One may obey different rules under
different situations. For instance, one will tend to believe the information
if half of the neighbors in the hyperedge believe the information, since all
the neighbours are his close friends. However, if he is less familiar with the
neighbors, he may tend to believe the information only unless all of the
neighbors in the hyperedge believe the information. Secondly, the resource
compartment also plays a role in information diffusion. While the direct
information spreading usually corresponds to private chat (pairwise) or group
chat (higher-order) via Whatsapp or some similar apps, the indirect spreading
are usually the posts on Twitter or videos on Tiktok. The posts or videos also
influence the agents to an extent depending on how popular the posts or videos
are. Thus, such posts or videos are resources, like water, and their
popularity and influence behave like the pathogen concentration correlated to
the infection level. Once the infection level is high, the post is getting
popular and in turn promote the spreading of the information. The higher-order
interaction is also necessary when the resource gets involved here. For
example, when one sees an interesting video on Tiktok, he or she may directly
forward it to a friend and have a discussion. Both the content of the video
and his or her conversation with a friend will jointly influence how strongly
one believes in the information.

To summarize, in this section, we have proposed a general model describing a
spreading process on a hypergraph up to 3-body interactions. Our model
considers any higher-order interaction and a general rule for spreading. The
model we have proposed can be used, for example, in the context of the
physical virus spreading and information diffusion. Furthermore, we take into
account both, the direct and indirect spreading processes. The rest of the
paper aims at achieving analytical results for the higher-order SIWS model and
its potential extension to deepen our understanding of a general spreading
process.

##  4 Analytical Results

Let us now consider a higher-order SIWS model on a hypergraph of $n$
populations and $m$ resources.

Component-wise, the system we consider is given by (8) and (9). We can further
derive it in a matrix form as

(10) $\dot{z}=-D_{f}z+(I-Z)(B_{f}z+H(z));$

where

(11) $\displaystyle D_{f}:=\left[\begin{matrix}D&\mathbf{0}\\\
\mathbf{0}&D^{w}\\\
\end{matrix}\right],\;Z=\left[\begin{matrix}\operatorname*{\mathsf{diag}}(x)&0\\\
0&0\\\ \end{matrix}\right],\;B_{f}=\left[\begin{matrix}B&B_{w}\\\ C_{w}&0\\\
\end{matrix}\right],$ $\displaystyle
H(z)=[z^{\top}B_{f1}z,\cdots,z^{\top}B_{fn+m}z]^{\top};$

and

(12)
$\begin{split}D&=\operatorname*{\mathsf{diag}}((\delta_{1},\cdots,\delta_{n})^{\top})\\\
D_{w}&=\operatorname*{\mathsf{diag}}((\delta^{w}_{1},\cdots,\delta^{w}_{m})^{\top})\\\
B_{w}&=[\beta^{w}_{ij}]_{n\times m}\\\ B&=[\beta_{ij}]_{n\times n}\\\
C_{w}&=[c^{w}_{jk}]_{m\times n}\end{split}$

Furthermore, if $i\leq n$, the matrix $B_{fi}$ is defined by

(13) $\begin{split}B_{fi}&=\left[\begin{matrix}[\beta_{ijk}]_{n\times
n}&[\beta^{xw}_{ijk}]_{n\times m}\\\ [\beta^{wx}_{ijk}]_{m\times
n}&\mathbf{0}\\\ \end{matrix}\right],\qquad\text{ if }i\leq n\\\
B_{fi}&=\left[\begin{matrix}[c^{xx}_{ijk}]_{n\times n}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\\\ \end{matrix}\right],\qquad\text{ if }i\geq
n+1.\end{split}$

Now, we present the analytical results of the general SIWS on a hypergraph
(10). Firstly, we need the following assumptions to ensure that our model is
well-defined.

######  Assumption 1.

The initial condition of the infection level satisfies $x_{i}(0)\in[0,1]$ for
all $i=1,\cdots,n$.

######  Remark 4.1.

Since $x_{i}(0)\in[0,1]$, according to (9),

(14) $w_{j}\leq
w_{j\max}\coloneqq\frac{\sum_{k=1}^{n}c_{jk}^{w}+\sum_{i,k\in\mathbf{N}_{3}^{\text{i}}}c_{jik}^{xx}}{\delta_{j}^{w}}.$

######  Assumption 2.

The matrices $B_{f}$, $B_{fi}$ are non-negative, $i=1,\cdots,n+m$. The
diagonal entries of $D_{f}$ are positive.

######  Remark 4.2.

We define later the reproduction number of the higher-order spreading as
$\rho(D_{f}^{-1}B_{f})$. To ensure the reproduction number is well-defined,
$D_{f}$ must be positive.

In the following Lemma, we describe the general properties of (10).

######  Remark 4.3.

For simplicity, we define the notation $F(z)=B_{f}z$,
$f_{i}=F(z)_{i}=(B_{f}z)_{i}$, $h_{i}=H(z)_{i}=z^{\top}B_{fi}z$. Then, the
concrete system (10) is a special case of the general system (37), which is
discussed later. Moreover, the concrete interaction function defined by
$F(z)=B_{f}z,\,H(z)_{i}=z^{\top}B_{fi}z$ satisfies Assumptions 8-11 for the
general system (37).

######  Lemma 4.4 (General properties).

Consider the system (10). Let $w_{j\max}$ be as in (14),
$\mathbf{D}=\left\\{z=(x,w)^{\top}=(x_{1},\ldots,x_{n},w_{1},\ldots,w_{n})^{\top}\,|\,x_{i}\in(0,1),\,w_{j}\in(0,w_{j\max})\right\\}$
and $\mathbf{\bar{D}}$ be its closure. If Assumptions 1-2 hold, then

  * i)

the set $\mathbf{\bar{D}}$ is positively invariant;

  * ii)

the set $\mathbf{D}$ is also positively invariant;

  * iii)

the origin is always an equilibrium of the model and there is no other
equilibrium on the boundary of $\mathbf{D}$.

######  Proof 4.5.

This lemma is a special case of Lemma 6.1. Thus, the proof is omitted here and
we provide the proof of a more general case in Lemma 6.1.

######  Remark 4.6.

The general properties of Lemma 4.4 show that the system is well-defined. The
set $\mathbf{\bar{D}}$ is the system’s domain for (10).

Now, we need one more assumption to achieve analytical results.

######  Assumption 3.

The matrix $B_{f}$ is irreducible.

######  Remark 4.7.

Assumption 3 implies that the graph of pairwise interactions is
strongly connected. Our model is originally based on a simplicial complex
structure. If any two nodes in the simplicial complex can be connected with a
sequence of edges (or hyperedges), then the graph of pairwise interaction is
strongly connected. Since we have performed some replacements introduced in
Section 3.1, the structure of the simplicial complex may break. However, since
the graph of pairwise interactions is at the lowest level, the property of its
strong connectivity does not break after the replacements. Hence, the
assumption is rather natural.

It is easy to check that the system (10) has an all-zero equilibrium, which is
the so-called _healthy state_.

######  Theorem 4.8 (Healthy-state behaviour).

Consider the system (10). If Assumptions 1-3 hold, then the healthy state
(zero equilibrium) always exists and is:

  1. 1.

locally stable if $\rho(D_{f}^{-1}B_{f})<1$,

  2. 2.

globally exponentially stable if
$\rho(D^{-1}_{f}B_{f}+D^{-1}_{f}[u^{\top}B_{f1},\cdots,u^{\top}B_{fn+m}]^{\top})<1$,
where $u=[1,\cdots,1,w_{1\max},\cdots,w_{m\max}]^{\top}$

  3. 3.

unstable if $\rho(D^{-1}_{f}B_{f})>1$.

######  Proof 4.9.

The local stability and instability of the healthy state is a special case of
Theorem 6.4. The proof is omitted here and we provide a proof for a more
general case in Theorem 6.4.

Next, we prove global stability. It is straightforward to check that

(15)
$D^{-1}_{f}B_{f}+D^{-1}_{f}[u^{\top}B_{f1},\cdots,u^{\top}B_{fn+m}]^{\top}$

is an irreducible non-negative matrix. Let $v$ be the strictly positive left
Perron-Frobenius eigenvector of the matrix (15) and $\lambda$ be the
corresponding Perron-Frobenius eigenvalue. Define the Lyuapunov function
$V=v^{\top}D_{f}^{-1}z$. We observe that $V\geq 0$ and $V=0$ if and only if
$z=\mathbf{0}$. Let $\delta=(\delta_{1},\cdots,\delta_{n})^{\top}$ and
$\delta^{w}=(\delta^{w}_{1},\cdots,\delta^{w}_{m})^{\top}$. Then, we get

(16) $\begin{split}\dot{V}&=-v^{\top}D_{f}^{-1}(D_{f}z+(I-Z)B_{f}z\\\
&+(I-Z)[z^{\top}B_{f1}z,\cdots,z^{\top}B_{fn+m}z]^{\top})\\\
&\leq-v^{\top}z+v^{\top}(D_{f}^{-1}B_{f}z\\\
&+D^{-1}_{f}[u^{\top}B_{f1}z,\cdots,u^{\top}B_{fn+m}z)^{\top}]\\\
&\leq(-1+\lambda)v^{\top}D_{f}D^{-1}_{f}z\leq(-1+\lambda)\min(\delta,\delta^{w})V\end{split}$

In view of the comparison lemma [22], we have $0\leq V(t)\leq
V(0)e^{(-1+\lambda)\min(\delta,\delta^{w})t}$, which further implies that
$x_{i}\leq\frac{\max(\delta^{-1},(\delta^{w})^{-1})v^{\top}z(0)}{\min(\delta^{-1},(\delta^{w})^{-1})v_{i}}e^{(-1+\lambda)\min(\delta,\delta^{w})t}$
and
$w_{j}\leq\frac{\max(\delta^{-1},(\delta^{w})^{-1})v^{\top}z(0)}{\min(\delta^{-1},(\delta^{w})^{-1})v_{j+n}}e^{(-1+\lambda)\min(\delta,\delta^{w})t}.$
Thus, we complete the proof.

######  Remark 4.10.

In epidemiology, the reproduction number (usually denoted as $R_{0}$) is
defined to be the expected number of cases directly generated by one infected
case in a population. Generally, if $R_{0}<1$, then the virus gradually dies
out, while if $R_{0}>1$, the pandemic breaks out. As for our model (10), we
suggest directly using $R_{0}=\rho(D_{f}^{-1}B_{f})$ as the reproduction
number, since it indicates whether the healthy state is locally stable or not.
Naturally, the behavior of the higher-order system (10) is different from the
SIWS model on a conventional graph studied by [41, 19, 45]. In the following,
we show that for the system (10), even if $R_{0}<1$, the pandemic may still
break out.

Other than the zero equilibrium, the system also processes a strictly positive
equilibrium, which is called _endemic equilibrium_.

######  Theorem 4.11.

(Existence of an endemic equilibrium with high infection levels) Consider the
system (10). If Assumptions 1-3 hold, then there exists an endemic equilibrium
$z^{*}\gg\mathbf{0}$ such that $x_{i}^{*}\geq\frac{2}{3}$ and
$w_{i}^{*}\geq\frac{2w_{i\max}}{3}$ for any $i$ such that
$B_{fi}\neq\mathbf{0}$ whenever
$\theta=\min_{i}\\{(D^{-1}_{f}B_{f})_{i}+[D^{-1}_{f}(\tilde{u}^{\top}B_{f1},\cdots,\tilde{u}^{\top}B_{fn+m})^{\top}]_{i}\\}\geq
4.5$, where
$\tilde{u}=[\tilde{u}_{1},\tilde{u}_{2},\cdots,\tilde{u}_{n+m}]^{\top}$ is
given by

$\tilde{u}_{i}=\begin{cases}1,&B_{fi}\neq\mathbf{0},\qquad i=1,\ldots,n\\\
w_{i\max},&B_{fi}\neq\mathbf{0},\qquad i=n+1,\ldots,n+m\\\
0,&B_{fi}=\mathbf{0}\end{cases}$

######  Proof 4.12.

Define the map

(17) $\begin{split}\hat{T}(z)&=T(D^{-1}_{f}F(z)+D^{-1}_{f}H(z));\\\
\end{split}$

where $T_{i}=\frac{p_{i}}{1+p_{i}}$for $i=1,\ldots,n$ and $T_{j}=p_{i}$ for
$i=n+1,\ldots,n+m$ with $p_{i}=(D^{-1}_{f}F(z)+D^{-1}_{f}H(z))_{i}$.

Next, we calculate the endemic equilibrium
$z^{*}=(x^{*},w^{*})^{\top}=(x_{1}^{*},\ldots,x_{n}^{*},w_{1}^{*}\ldots,w_{m}^{*})^{\top}$
by setting $\dot{x}_{i}$ and $\dot{w}_{j}$ to zero. This leads to
$x_{i}^{*}=\frac{p_{i}(z^{*})}{1+p_{i}(z^{*})}$ and
$w^{*}_{j}=p_{j+n}(z^{*})$, which implies that the endemic equilibrium $z^{*}$
is a fixed point of $\hat{T}$, i.e., $\hat{T}(z^{*})=z^{*}$. We know that both
$\frac{p}{1+p}$ and $p$ are monotonically increasing in $z\geq 0$. Thus, the
maps $\hat{T}(p(z))$ and $T(p(z))$ are monotonically increasing in $z\geq 0$.
Moreover, if $\alpha>1$, $p\leq 1-\frac{1}{\alpha}$ with $p\geq 0$ if and only
if $\alpha\geq\frac{1}{1-p}$. If $\alpha>1$ and $p>0$, $\alpha p>p$. Thus, we
have $T(\alpha p(z))\geq p(z)$ with $p(z)\geq\mathbf{0}$ if and only if
$p_{i}\leq 1-\frac{1}{\alpha}$ for $i\leq n$.

Then, we have that

$\begin{split}\hat{T}(z)&=T(D^{-1}_{f}B_{f}z+D^{-1}_{f}[z^{\top}B_{f1}z,z^{\top}B_{f2}z,\cdots]^{\top})\\\
&\geq
T\left(\frac{2}{3}D^{-1}_{f}B_{f}\tilde{u}+\frac{4}{9}D_{f}^{-1}[\tilde{u}^{\top}B_{f1}\tilde{u},\tilde{u}^{\top}B_{f2}\tilde{u},\cdots]^{\top}\right)\\\
&\geq
T\left(\frac{4}{9}\theta\tilde{u}\right)\geq\frac{2}{3}\tilde{u}.\end{split}$

By using the monotonicity, we derive that $u$ must be the upper bound. In
other words, $\hat{T}$ maps the set

(18)
$\begin{split}\Big{\\{}z&=(x,w)^{\top}=(x_{1},\cdots,x_{n},w_{1},\cdots,w_{m})^{\top}\,|\\\
&\frac{2}{3}\tilde{u}_{i}\leq x_{i}\leq
1,\frac{2w_{j\max}}{3}\tilde{u}_{j+n}\leq w_{j}\leq w_{j\max},\\\
&i=1,\cdots,n,\,j=1,\cdots,m\Big{\\}}\end{split}$

into itself. By using Brouwer’s Fixed-Point Theorem [52], we complete the
proof.

Furthermore, we can show that the endemic equilibrium in Theorem 4.11 is
locally stable.

######  Theorem 4.13.

(Bi-stability may occur) Consider the model (10). If Assumptions 1-3 hold and
all $B_{fi}\neq\mathbf{0}$, any endemic equilibriunm in the form of
$z^{*}\gg\mathbf{0}$ such that $\frac{2}{3}\leq x_{i}^{*}\leq 1$ and
$\frac{2w_{\max}}{3}\leq w_{j}^{*}\leq w_{\max}$ is locally stable.

######  Proof 4.14.

Firstly, we calculate the Jacobian for the whole system:

$\begin{split}J&=-D_{f}+(I-{Z})\Big{(}B_{f}-\hat{I}\operatorname*{\mathsf{diag}}([z^{\top}B_{f1}z,\cdots,z^{\top}B_{fn+m}z]^{\top}\\\
&+B_{f}z+[z^{\top}(B_{f1}+B_{f1}^{\top}),\cdots,z^{\top}(B_{fn+m}+B_{fn+m}^{\top})]^{\top})\Big{)}.\end{split}$

It is easy to check that $J$ is an irreducible Metzler matrix.

Next, we define the vector $\tilde{z}^{*}=(\frac{x^{*}}{2},w^{*})^{\top}$.
Then,

$\begin{split}J\tilde{z}^{*}&=-D_{f}\tilde{z}^{*}+(I-{Z})(B_{f}\tilde{z}^{*}\\\
&+[z^{\top}(B_{f1}+B_{f1}^{\top})\tilde{z}^{*},\cdots,z^{\top}(B_{fn+m}+B_{fn+m}^{\top})\tilde{z}^{*}]^{\top})\\\
&-\hat{I}(\operatorname*{\mathsf{diag}}([z^{\top}B_{f1}z,\cdots,z^{\top}B_{fn+m}z]^{\top}+B_{f}z)\tilde{z}^{*})\end{split}$

with $\hat{I}=\left[\begin{matrix}I&\mathbf{0}\\\ \mathbf{0}&\mathbf{0}\\\
\end{matrix}\right]$.

Element-wise, we get the following. For $n+1\leq i\leq n+m$,

$\begin{split}(J(z^{*})\tilde{z}^{*})_{i}&=-\delta^{w}_{i-n}w^{*}_{i-n}+\sum_{j}c^{w}_{(i-n)j}\frac{x^{*}_{j}}{2}\\\
&+2\sum_{i,j}c^{xx}_{(i-n)jk}x^{*}_{j}\frac{x^{*}_{k}}{2}\\\
&=-\sum_{j}c^{w}_{(i-n)j}\frac{x^{*}_{j}}{2}<-d_{i}w^{*}_{i-n}=-d_{i}z_{i}^{*},\end{split}$

with some $d_{i}>0$.

On the other hand, for $i\leq n$,

$\begin{split}(J(z^{*})\tilde{z}^{*})_{i}&=-\frac{x_{i}^{*}}{2}(\sum_{j}\beta_{ij}x^{*}_{j})+(\frac{1}{2}-x^{*}_{i})\sum_{j}\beta_{ij}^{w}\frac{w^{*}_{j}}{2}\\\
&+\frac{1}{2}(2-3x^{*}_{i})(\sum_{j,k}\beta^{wx}_{ijk}x^{*}_{j}w^{*}_{k}+\sum_{j,k}\beta^{wx}_{ijk}x^{*}_{k}w^{*}_{j})\\\
&+(\frac{1}{2}-x^{*}_{i})\sum_{j,k}\beta^{xx}_{ijk}x^{*}_{j}x^{*}_{k}.\end{split}$

Since $2-3x_{i}^{*}<0$ and $\frac{1}{2}-x^{*}_{i}<0$, we observe that
$(J\tilde{z}^{*})_{i}<-d_{i}\frac{z_{i}^{*}}{2}$ for $i\leq n$ with some
$d_{i}>0$. By combining the results from the case $n+1\leq i\leq n+m$, we have
$J\tilde{z}^{*}<-d\tilde{z}^{*}$ with some $d>0$. According to Theorem 10.14
in [3], we confirm that $J$ is Hurwitz and the endemic equilibrium we consider
here is locally stable.

Notice that the condition guaranteeing the existence of an endemic equilibrium
in the form of Theorem 4.11 and its local stability can hold together with the
condition that the reproduction number is less than $1$. This is an
interesting phenomenon of the higher-order SIWS model on hypergraphs because
we see that the system may have a multi-stability domain when the reproduction
number is smaller than $1$. In such a domain, the healthy state and the
endemic equilibrium in the form of Theorem 4.11 are locally stable. Different
from the higher order system, the SIWS on a conventional graph has one unique
globally stable equilibrium, healthy state, when the reproduction number is
less than $1$. We will give some numerical examples of this case later in the
section on numerical examples (Figure 3).

######  Theorem 4.15.

(Existence of an endemic equilibrium when $R_{0}>1$) Consider the model (10),
if Assumptions 1-3 hold, there exists an endemic equilibrium
$z^{*}\gg\mathbf{0}$ if $R_{0}=\rho=\rho(D^{-1}_{f}B_{f})>1$.

######  Proof 4.16.

Define the set $\tilde{Y}=\\{z=(x,w)^{\top}|c_{i}\leq x_{i}\leq 1,c_{j+n}\leq
w_{j}\leq w_{j\max}\\}$, where $c=(c_{1},\cdots,c_{n+m})^{\top}=\alpha v$, $v$
is the positive Perron-eigenvector associated with
$\rho=\rho(D^{-1}_{f}B_{f})$ and $\alpha$ is small enough such that $c=\alpha
v\leq 1-\frac{1}{\rho}$. We use again the map (17). For any $z\in\tilde{Y}$,
we have
$\hat{T}(z)=T(D^{-1}_{f}B_{f}z+D^{-1}_{f}[z^{\top}B_{f1}z,z^{*\top}B_{f2}z,\cdots]^{\top})\geq
T(D^{-1}_{f}B_{f}z)\geq T(\alpha\rho v)\geq T(\rho c)\geq c$. Notice that the
map is upper bounded by $u$ when $z$ is in the system domain. Using again the
Brouwer’s fixed-Point Theorem [52], one confirms that there exists an endemic
equilibrium $z^{*}\gg\mathbf{0}$ if $\rho(D^{-1}_{f}B_{f})>1$.

Now, we show that the endemic equilibrium is unique and globally stable if
$\rho(D^{-1}_{f}B_{f})>1$ and the $\beta^{w}_{ij}$ terms in $B_{f}$ and all
the higher order parameters in all $B_{fi}$ are sufficiently small.

######  Theorem 4.17.

(Global stability of the endemic equilibrium) Consider the system (10), if
Assumptions 1-3 hold and $R_{0}=\rho(D^{-1}_{f}B_{f})>1$, then the endemic
equilibrium is unique and globally exponentially stable if the
$\beta^{w}_{ij}$ terms in $B_{f}$ and all the higher order parameters in all
$B_{fi}$ are sufficiently small.

######  Proof 4.18.

We recall that $\tilde{Y}=\\{z=(x,w)^{\top}|c_{i}\leq x_{i}\leq 1,c_{j+n}\leq
w_{j}\leq w_{j\max}\\}$, where $c=(c_{1},\cdots,c_{n+m})^{\top}=\alpha v$, $v$
is the positive Perron-eigenvector associated with
$\rho=\rho(D^{-1}_{f}B_{f})$ and $\alpha$ is small enough such that $c=\alpha
v\leq 1-\frac{1}{\rho}$. It is easy to check that $\tilde{Y}$ is a forward
invariant set. By using the condition $c_{i}=\alpha v_{i}$, it yields that
$\dot{z}_{i}>0$ when $z_{i}=c_{i}$.

Then, the idea is to rewrite (10) in the form of $\dot{z}=D(z,z^{*})(z-z^{*})$
and apply Theorem 3.1 (Exponential convergence from Coppel’s inequality) in
[7]. More precisely, we can rewrite (10) as $\dot{z}=D(z,z^{*})(z-z^{*}),$
where

$\begin{split}&D(z,z^{*}):=-D_{f}+(I-{Z}^{*})B_{f}-\operatorname*{\mathsf{diag}}(\left[\begin{matrix}B&B_{w}\\\
\mathbf{0}&\mathbf{0}\\\ \end{matrix}\right]z)\\\
&+(I-{Z}^{*})\left[\begin{matrix}\frac{1}{2}(z+z^{*})^{\top}B_{f1}+\frac{1}{2}(z+z^{*})^{\top}B_{f1}^{\top}\\\
\vdots\\\
\frac{1}{2}(z+z^{*})^{\top}B_{fn+m}+\frac{1}{2}(z+z^{*})^{\top}B_{fn+m}^{\top}\end{matrix}\right]\\\
&-\operatorname*{\mathsf{diag}}[(z^{\top}B_{f1}z,\cdots,z^{\top}B_{fn}z,0,\cdots,0)^{\top}]\end{split}$

is an irreducible Metzler matrix.

Again, we define $\tilde{z}^{*}=(\frac{x^{*}}{2},w^{*})^{\top}$. For $i\geq
n+1$, we have

$\begin{split}(D\tilde{z}^{*})_{i}=-\frac{1}{2}\delta^{w}_{i-n}w_{i-n}^{*}+\frac{1}{4}\sum_{j,k}c^{xx}_{(i-n)jk}(x_{i}^{*}x_{k}+x_{i}x_{k}^{*}).\end{split}$

On the other hand, for $i\leq n$ we have

$\begin{split}(D\tilde{z}^{*})_{i}&=(1-x_{i}^{*})\Big{(}\sum_{j}\frac{1}{2}\beta^{w}_{ij}w_{j}^{*}\Big{)}\\\
&-\frac{1}{2}x^{*}_{i}\Big{(}\sum_{j}\beta_{ij}x_{j}+\sum_{j}\beta^{w}_{ij}w_{j}+\sum_{j,k}\beta_{ijk}x_{i}x_{j}\\\
&+\sum_{j,k}\beta_{ijk}^{wx}w_{j}x_{k}+\sum_{j,k}\beta_{ijk}^{xw}x_{j}w_{k}\Big{)}\\\
&+(1-x_{i}^{*})\Big{(}\frac{1}{4}\sum_{j,k}\beta^{xw}_{ijk}(x^{*}_{j}w_{k}+x_{j}w^{*}_{k}+x^{*}_{j}w^{*}_{k})\\\
&+\sum_{j,k}\frac{1}{4}\beta^{wx}_{ijk}(w^{*}_{j}x^{*}_{k}+w^{*}_{j}x_{k}+w_{j}x^{*}_{k})\\\
&+\sum_{j,k}\frac{1}{4}\beta^{xx}_{ijk}(x^{*}_{j}x_{k}+x_{j}x^{*}_{k})\Big{)}.\end{split}$

One can confirm that if the $\beta^{w}_{ij}$ terms in $B_{f}$ and all the
higher order parameters in all $B_{fi}$
($\beta^{xw}_{ijk},\beta^{xx}_{ijk},\beta^{wx}_{ijk},c^{xx}_{ijk}$) are
sufficiently small, then $D\tilde{z}^{*}$ is upper bounded by some negative
constant. Thus, according to Theorem 3.1 (Exponential convergence from
Coppel’s inequality) in [7], we complete the proof.

Next, we show that the system (10) is actually an irreducible monotone system.

######  Theorem 4.19 (Irreducible monotone system).

If Assumption 1-3 hold, the system (10) is an irreducible monotone system in
$\mathbf{{D}}$ (defined in Lemma 4.4). Furthermore, the model (10) has a
finite number of equilibria in $\mathbf{\bar{D}}$, which is the closure of
$\mathbf{{D}}$ (defined in Lemma 4.4), for a generic choice of parameters.
Moreover, the solution of (10) converges to an equilibrium for almost all
initial conditions. That is to say, the set of initial conditions, such that
the model does not converge to an equilibrium, is a set of Lebesgue measure
zero.

######  Proof 4.20.

Firstly, we calculate the Jacobian of the system (10):

$\displaystyle
J=-D_{f}+(I-Z)(B_{f}+J_{H}(z))-Z\operatorname*{\mathsf{diag}}(B_{f}z+H(z));$

where $J_{H}(z)$ is the Jacobian of $H(z)$. Note that
$B_{f}+J_{H}(z)>\mathbf{0}$. Furthermore, $D_{f}$ is diagonal. Thus, the
Jacobian is an irreducible Metzler matrix at an arbitrary point in
$\mathbf{{D}}$ along the trajectories of (10). This confirms that the system
(10) is an irreducible monotone system by Lemma 2.2 of [64].

We just need to further confirm that the system (10) has a finite number of
equilibria in the closure of $\mathbf{\bar{D}}$ for a generic choice of
parameter. We know that the healthy state is globally exponentially stable if
$\rho(D^{-1}_{f}B_{f}+D^{-1}_{f}[u^{\top}B_{f1},\cdots,u^{\top}B_{fn+m}]^{\top})<1$,
which indicates that the healthy state is a unique equilibrium. Another
example could be to choose the infection rates (some are zeros, the rest take
positive values) and recovery rates (set to be 1) such that it yields a
partially-decoupled equation of
$\\{-1+(1-x_{i})[\beta_{ii}+\beta_{iii}x_{i}]\\}x_{i}=0$ and
$-w_{i}+[c_{ij}x_{j}+c_{ijj}(x_{j})^{2}]=0$. It is straightforward to see that
such an equation system has a finite number of solutions. We indeed see some
particular examples. Thus, we found out that the model (10) has a finite
number of equilibria in the closure of $\mathbf{\bar{D}}$ with a particular
choice of parameter. We can also see that the equation $\dot{z}=\mathbf{0}$ is
actually a polynomial equation. Thus, the model has a finite number of
equilibria in the closure of $\mathbf{\bar{D}}$ if the parameters do not lie
on a certain algebraic set of measure zero by Theorem B.1 and Corollary B.2 in
[64].

Moreover, by Lemma 2.3 of [64], the solutions of (10) converge to an
equilibrium with some domain of attraction; by the same lemma, the set of
initial conditions, such that the model does not converge to an equilibrium,
is a set of measure zero.

######  Remark 4.21.

As a consequence of (10) being an irreducible monotone system and according to
the Theorem 4.19, the convergence to equilibrium occurs for almost all initial
conditions. Once the system processes a unique locally stable equilibrium, it
is then globally stable. If the healthy state is unstable, then the
convergence to an endemic equilibrium happens for almost all initial
conditions. However, there may be multiple endemic equilibria. In the case of
multiple locally stable equilibria, then multi-stability is induced. However,
we have not determined the domain of attraction of each locally stable
equilibrium yet. This remains future work.

##  5 Extension to a bi-virus competing SIWS model with higher-order
interactions

In most scenarios, not just one virus (or information) spreads over the
network. Usually, several different viruses spread simultaneously and compete
with each other. Concrete examples of this competition in real life include
competing opinions on a social network [50] and competing rumors spreading in
politics [56]. Here, we study the simplest case of a bi-virus competition as a
typical example. Since we concentrate on a case of extreme competition, we
assume both viruses are not compatible with each other and one can only get
infected with one virus at a time. We want to emphasize that this competitive
assumption and coupling method are widely adopted among related works in
epidemic models on a conventional graph [32, 8]. Then, a bi-virus model can be
obtained in at least two ways. One can directly couple two separate single-
virus models (10) into one system. Alternatively, one can first model a Markov
chain for the bi-virus case and repeat the procedure in the section 3.1 and
3.2.

The model (10) can be straightforwardly extended into a bi-virus competing
SIWS model, based on the assumption of no simultaneous co-infection, which
leads to

$\displaystyle\dot{x}^{1}_{i}$
$\displaystyle=-\delta_{i}^{1}x_{i}^{1}+\left[(1-x_{i}^{1}-x_{i}^{2})(\sum_{j=1}^{n}\beta^{1}_{ij}x_{j}^{1}+\sum_{j=1}^{m}\beta_{ij}^{w1}w_{j}^{1})\right]$
(19)
$\displaystyle+\left[(1-x_{i}^{1}-x_{i}^{2})\left(\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta^{xx1}_{ijk}x_{i}^{1}x_{j}^{1}+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{ijk}^{wx1}w_{j}^{1}x_{k}^{1}+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{ijk}^{xw1}x_{j}^{1}w_{k}^{1}\right)\right],$
(20) $\displaystyle\dot{w}_{j}^{1}$
$\displaystyle=-\delta_{j1}^{w}w_{j}^{1}+\sum_{k=1}^{n}c_{jk}^{w1}x_{k}^{1}+\sum_{i,k\in\mathbf{N}_{3}^{\text{i}}}c_{ijk}^{xx1}x^{1}_{i}x^{1}_{k},$
$\displaystyle\dot{x}^{2}_{i}$
$\displaystyle=-\delta_{i}^{2}x_{i}^{2}+\left[(1-x_{i}^{1}-x_{i}^{2})(\sum_{j=1}^{n}\beta_{ij}^{2}x_{j}^{2}+\sum_{j=1}^{m}\beta_{ij}^{w2}w_{j}^{2})\right]$
(21)
$\displaystyle+\left[(1-x_{i}^{1}-x_{i}^{2})\left(\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta^{xx2}_{ijk}x_{i}^{2}x_{j}^{2}+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{ijk}^{wx2}w_{j}^{2}x_{k}^{2}+\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta_{ijk}^{xw2}x_{j}^{2}w_{k}^{2}\right)\right],$
(22) $\displaystyle\dot{w}_{j}^{2}$
$\displaystyle=-\delta_{j}^{w2}w_{j}^{2}+\sum_{k=1}^{n}c_{jk}^{w2}x_{k}^{2}+\sum_{i,k\in\mathbf{N}_{3}^{\text{i}}}c_{ijk}^{xx2}x^{2}_{i}x^{2}_{k},$

where all notations remain the same as Table 1 and the superscript $1,2$,
denotes the first or the second virus. Without loss of generality and for
simplicity, we refer to $\nu=1$ as the first virus and to $\nu=2$ as the
second. We recall that, for now, we deal with the interaction of up to 3
bodies.

For $\nu=1,2$, we can rewrite our system (19)-(22) as

(23)
$\begin{split}\dot{z}^{\nu}&=-D_{f}^{\nu}z^{\nu}+(I-Z^{1}-Z^{2})(B_{f}^{\nu}z^{\nu}+H^{\nu}(z^{\nu})).\end{split}$

In the following, we present the bi-virus version of assumptions 1, 2, and 3.

######  Assumption 4.

The initial condition of infection level satisfies $x^{\nu}_{i}(0)\in[0,1]$
for all for all $\nu=1,2$ and $i=1,\ldots,n$. Furthermore,
$x^{1}_{i}(0)+x^{2}_{i}(0)\leq 1$ holds for all $i=1,\ldots,n$.

######  Remark 5.1.

We can now obtain the upper bound of each pathogen dynamic as

(24)
$w^{\nu}_{j\max}=\frac{\sum_{k=1}^{n}c_{jk}^{w\nu}+\sum_{i,k\in\mathbf{N}_{3}^{\text{i}}}c_{jik}^{xx\nu}}{\delta_{j}^{w}}.$

######  Assumption 5.

The matrices $B_{f}^{\nu}$, $B_{fj}^{\nu}$ are non-negative for $\nu=1,2$. The
diagonal entries of $D_{f}^{\nu}$ are positive for $\nu=1,2$.

######  Assumption 6.

The matrix $B^{\nu}_{f}$ is irreducible for $\nu=1,2$.

For simplicity, we define the notation $F(z)^{\nu}=B_{f}^{\nu}z^{\nu}$,
$f_{i}^{\nu}=(B_{f}^{\nu}z^{\nu})_{i}$, $h_{i}^{\nu}=(H^{\nu})_{i}$. Then, the
concrete system (23) is a special case of the general system (46), which is
discussed later.

######  Lemma 5.2 (General properties).

Consider the bi-virus model (23). Let $w^{k}_{j\max}$ be as defined in (24),

(25)
$\begin{split}\mathbf{D}&=\left\\{z|z=(x^{1},w^{1},x^{2},w^{2})^{\top}=(x^{1}_{1},\cdots,x^{1}_{n},w^{1}_{1},\cdots,w^{1}_{m},x^{2}_{1},\cdots,x^{2}_{n},w^{2}_{1},\cdots,w^{2}_{m})^{\top},\right.\\\
&\qquad\left.x^{\nu}_{i}\in(0,1),\,w^{\nu}_{j}\in(0,w^{\nu}_{j\max}),\,\nu=1,2\right\\},\end{split}$

and $\mathbf{\bar{D}}$ be its closure. If Assumptions 4-5 hold, then

  * i)

the set $\mathbf{\bar{D}}$ is positively invariant;

  * ii)

the set $\mathbf{D}$ is also positively invariant;

  * iii)

the origin is always an equilibrium of the model;

  * iv)

the dynamics of each virus in (23) are upper bounded by its single virus
counterpart in (10) with the same system parameters.

######  Proof 5.3.

This lemma is a special case of Lemma 6.17. Thus, the proof is omitted here
and we provide the proof of a more general case in Lemma 6.17.

######  Remark 5.4.

Notice that now we study a bi-virus system (23). Hence, the system domain
$\mathbf{D}$ is different from the system domain of a single-virus system in
the former section (10). We keep the same notation because we are referring to
the equivalent physical domain.

######  Theorem 5.5 (Healthy-state behaviour).

Consider the system (23). If Assumptions 4-6 hold, then the healthy state
(zero equilibrium) always exists and is

  1. 1.

locally stable if $\rho((D^{\nu}_{f})^{-1}B_{f}^{\nu})<1$ for $\nu=1,2$,

  2. 2.

globally exponentially stable if

$\rho((D^{\nu}_{f})^{-1}B^{\nu}_{f}+(D^{\nu}_{f})^{-1}[u^{\nu\top}B^{\nu}_{f1},\cdots,u^{\nu\top}B^{\nu}_{fn+m}]^{\top})<1,$

with
${u}^{\nu}=[\mathbf{1}^{n},w^{\nu}_{1\max},\cdots,w^{\nu}_{m\max}]^{\top}$,
and for $\nu=1,2$,

  3. 3.

unstable if $\rho((D^{\nu}_{f})^{-1}B_{f}^{\nu})>1$ for either $\nu=1$ or
$\nu=2$.

######  Proof 5.6.

It is clear that the healthy state (zero equilibrium) always exists. The
Jacobian matrix of the model at the origin reads as

$J(0,0)=\left[\begin{matrix}-D^{1}_{f}+B_{f}^{1}&\mathbf{0}\\\
\mathbf{0}&-D^{2}_{f}+B_{f}^{2}\\\ \end{matrix}\right].$

Thus, the healthy state is locally stable if
$\rho((D^{\nu}_{f})^{-1}B_{f}^{\nu})<1$ for $\nu=1,2$ and unstable if
$\rho((D^{\nu}_{f})^{-1}B_{f}^{\nu})>1$ for either $\nu=1$ or $\nu=2$.

The global exponential stability when

$\rho((D^{\nu}_{f})^{-1}B^{\nu}_{f}+(D^{\nu}_{f})^{-1}[u^{\nu\top}B^{\nu}_{f1},\cdots,u^{\nu\top}B^{\nu}_{fn+m}]^{\top})<1$

for $\nu=1,2$, is a direct consequence of statement iv) in Lemma 5.2, Theorem
4.8 and the comparison principle.

######  Theorem 5.7 (Dominant endemic equilibrium).

Consider the bi-virus model (23), if Assumptions 4-6 hold, we have the
following statements. If $z^{*}=(x^{*},w^{*})^{\top}$ is an endemic
equilibrium of the single-virus counterpart (10) of the first virus, then
$(z^{*},0,0)^{\top}$ is the dominant endemic equilibrium of (23). The dominant
endemic equilibrium is locally stable when $z^{*}$ is locally stable and
$\rho(-D^{2}_{f}+(I-Z^{*})B_{f}^{2})<0$ with respect to the single-virus
system and is unstable when $z^{*}$ is unstable or
$\rho(-D^{2}_{f}+(I-Z^{*})B_{f}^{2})>0$, where
$Z^{*}=\left[\begin{matrix}x^{*}&0\\\ 0&0\\\ \end{matrix}\right]$.
Furthermore, if $\rho((D^{1}_{f})^{-1}B^{1}_{f})>1$,
$\rho((D^{2}_{f})^{-1}B^{2}_{f}+(D^{2}_{f})^{-1}[u^{2\top}B_{f1},\cdots,u^{2\top}B_{fn+m}]^{\top})<1$,
the $\beta^{w1}_{ij}$ terms in $B_{f}$ and all the higher order parameters in
all $B^{1}_{fi}$ are sufficiently small, then the dominant endemic equilibrium
is globally stable with domain of attraction
$\bar{\mathbf{D}}\setminus\mathcal{H}$, where
$\mathcal{H}=\\{[z^{1},z^{2}]^{\top}\,|\,z^{2}\in\bar{\mathbf{D}},z^{1}=\mathbf{0}\\}$.

######  Proof 5.8.

If we set $z^{2}=\mathbf{0}$, the system (23) is reduced into the single virus
counterpart of the first virus. Thus, if $z^{*}$ is an endemic equilibrium of
the single-virus counterpart (10) of the first virus, then
$(z^{*},0,0)^{\top}$ is the dominant endemic equilibrium of the system (23).

The Jacobian matrix of the model at the dominant endemic equilibrium is

$J(z^{*},\mathbf{0})=\left[\begin{matrix}\hat{J}(z^{*})&\mathbf{Res}\\\
\mathbf{0}&-D^{2}_{f}+(I-Z^{*})B_{f}^{2}\\\ \end{matrix}\right],$

where $\hat{J}(z^{*})$ is the Jacobian of the single-virus counterpart at the
endemic equilibrium $z^{*}$ and $\mathbf{Res}$ is the matrix of some residual
which we do not need to deal with.

By observing the Jacobian $J(z^{*},\mathbf{0})$ and notice that
$\hat{J}(z^{*})$ and $-D^{2}_{f}+(I-Z^{*})B_{f}^{2}$ is irreducible
nonnegative, we can directly obtain the results of local stability. That is,
the dominant endemic equilibrium is locally stable when $z^{*}$ is locally
stable and $\rho(-D^{2}_{f}+(I-Z^{*})B_{f}^{2})<0$ with respect to the single-
virus system and is unstable when $z^{*}$ is unstable or
$\rho(-D^{2}_{f}+(I-Z^{*})B_{f}^{2})>0$.

For the case of global stability, since

$\rho((D^{2}_{f})^{-1}B^{2}_{f}+(D^{2}_{f})^{-1}[u^{2\top}B_{f1},\cdots,u^{2\top}B_{fn+m}]^{\top})<1$

holds, the single-virus counterpart of the second virus converges to the
healthy state. By the comparison principle, the dynamics of the second virus
converge to zero. The dynamic of the first virus is thus its single virus
counterpart plus a vanishing perturbation. Therefore, the whole system
converges to the dominant endemic equilibrium.

######  Remark 5.9.

The analytical results of the Theorem 5.7 are still true if we interchange the
superscript $1$ (first virus) and $2$ (second virus). Bistability may also
occur in the bi-virus system because the condition for the local stability of
the healthy state and a dominant endemic equilibrium may be satisfied at the
same time.

######  Theorem 5.10 (Coexisting equilibrium).

Consider the model (23), if Assumptions 4-6 hold, we have the following
statements. If $\rho((D^{\nu}_{f})^{-1}B^{\nu}_{f})>1$ for all $\nu=1,2$ and
$\rho(-D^{1}_{f}+(I-Z^{2*})B^{1}_{f})>0$ as well as
$\rho(-D^{2}_{f}+(I-Z^{1*})B^{2}_{f})>0$, where $z^{1*}$ and $z^{2*}$ are the
endemic equilibrium of the single-virus counterpart of the first and second
virus respectively and the notation $Z^{1*}$ and $Z^{2*}$ sticks to the form
of (11), then there exists at least one coexisting equilibrium
$(\hat{z}^{1},\hat{z}^{2})\gg\mathbf{0}$ such that
$\hat{x}^{1}+\hat{x}^{2}\leq\mathbf{1}$, where
$\hat{z}^{\nu}=(\hat{x}^{\nu},\hat{w}^{\nu})^{\top}$.

######  Proof 5.11.

Let us define the map $M:\mathbb{R}^{2n+2m}\rightarrow\mathbb{R}^{2n+2m}$ as
follows: component-wise, $M$ is given by

(26) $\displaystyle M^{1}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=\frac{\left(1-z_{i}^{2}\right)\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}z^{1}+\left(D_{f}^{1}\right)^{-1}H^{1}(z^{1})\right)_{i}}{1+\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}z^{1}+\left(D_{f}^{1}\right)^{-1}H^{1}(z^{1})\right)_{i}}$
$\displaystyle M^{2}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=\frac{\left(1-z_{i}^{1}\right)\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}z^{2}+\left(D_{f}^{2}\right)^{-1}H^{2}(z^{2})\right)_{i}}{1+\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}z^{2}+\left(D_{f}^{2}\right)^{-1}H^{2}(z^{2})\right)_{i}},$

for $i\leq n$. Furthermore, for $i\geq n+1$, we let

(27) $\displaystyle M^{1}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=(D_{f}^{1})^{-1}(B_{f}^{1}z^{1}+H^{1}(z^{1}))_{i},$
$\displaystyle M^{2}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=(D_{f}^{2})^{-1}(B_{f}^{2}z^{2}+H^{2}(z^{2}))_{i}.$

It is straightforward to check that the fixed point of the map
$M(z)=M(z^{1},z^{2})=(M^{1}(z^{1},z^{2}),M^{2}(z^{1},z^{2}))^{\top}$ yields
the equilibrium point of the model (23). Furthermore, the map $M^{1}$ is a
non-increasing function of $z^{2}$ and a non-decreasing function of $z^{1}$.
Similarly, the map $M^{2}$ is a non-decreasing function of $z^{2}$ and a non-
increasing function of $z^{1}$.

Suppose $z^{1*}$ and $z^{2*}$ are the endemic equilibrium of the single-virus
counterpart (10) of the first and second virus respectively. It follows that
$M_{1}(z^{1*},\mathbf{0})=z^{1*}$ and $M_{2}(\mathbf{0},z^{2*})=z^{2*}$. Then,
we look at the map $M_{1}$. If $z^{2}\geq\mathbf{0}$, then it holds
$M_{1}(z^{1*},z^{2})\leq z^{1*}$. Moreover, if $z^{1}\leq z^{1*}$, then
$M_{1}(z^{1},z^{2})\leq z^{1*}$. Similarly, we have $M_{2}(z^{1},z^{2})\leq
z^{2*}$ if $z^{1}\geq\mathbf{0}$ and $z^{2}\leq z^{2*}$. Combining these two
results, we get $M(z^{1},z^{2})\leq M(z^{1*},z^{2*})$ if
$(z^{1},z^{2})\leq(z^{1*},z^{2*})$.

Now, we construct another map $\hat{M}$. Component-wise, the expression for
each row is given by,

(28) $\displaystyle\hat{M}^{1}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=\frac{\left(1-z_{i}^{2}\right)\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}z^{1}\right)_{i}}{1+\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}z^{1}\right)_{i}}\leq
M^{1}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle\hat{M}^{2}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=\frac{\left(1-z_{i}^{1}\right)\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}z^{2}\right)_{i}}{1+\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}z^{2}\right)_{i}}\leq
M^{2}\left(z^{1},z^{2}\right)_{i};$

for $i\leq n$. These inequalities are true because the function
$\frac{s}{1+s}$ is increasing as $s$ increases. Then, similar proof techniques
from Theorem 5 of [19] can be further applied.

Then, for $i\geq n+1$,

(29) $\displaystyle\hat{M}^{1}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=(D_{f}^{1})^{-1}(B_{f}^{1}z^{1})\leq
M^{1}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle\hat{M}^{2}\left(z^{1},z^{2}\right)_{i}$
$\displaystyle=(D_{f}^{2})^{-1}(B_{f}^{2}z^{2})\leq
M^{2}\left(z^{1},z^{2}\right)_{i}.$

Hence, it follows that $\hat{M}(z)\leq M(z)$. Due to the condition
$\rho(-D^{1}_{f}+(I-Z^{2*})B^{1}_{f})>0$ as well as
$\rho(-D^{2}_{f}+(I-Z^{1*})B^{2}_{f})>0$ and also considering that
$-D^{1}_{f}+(I-Z^{2*})B^{1}_{f}$ and $-D^{2}_{f}+(I-Z^{1*})B^{2}_{f}$ are both
irreducible Metzler matrices, we can obtain
$\rho((I-Z_{3-\nu}^{*})(D^{\nu}_{f})^{-1}B^{\nu}_{f})>1$, for $\nu=1,2$. Note
that $(I-Z_{3-\nu}^{*})(D^{\nu}_{f})^{-1}B^{\nu}_{f}$ is irreducible
nonnegative. Let $\lambda^{\nu}>1$ and $y^{\nu}$ be its corresponding simple
eigenvalue and positive eigenvector. Since $(D^{\nu}_{f})^{-1}B^{\nu}_{f}$ is
also irreducible nonnegative, we have
$((D^{\nu}_{f})^{-1}B^{\nu}_{f}y^{\nu})_{k}>0$. Hence, there must exist
$\epsilon^{1}>0$ and $\epsilon^{2}>0$ such that

(30)
$\begin{array}[]{r}\epsilon^{1}<\min\left\\{\frac{\lambda^{1}-1}{\max_{i\in[n+1]}\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}{{y}}^{1}\right)_{i}},\min_{i\in[n+1]}\frac{({z}_{1}^{*})_{i}}{{{y}}_{i}^{1}}\right\\}\\\
\epsilon^{2}<\min\left\\{\frac{\lambda^{2}-1}{\max_{i\in[n+1]}\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}{{y}}^{2}\right)_{i}},\min_{i\in[n+1]}\frac{({z}_{2}^{*})_{i}}{{{y}}_{i}^{2}}\right\\}\end{array}$

It follows that

(31) $\displaystyle
1+\max_{i\in[n]}\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}\epsilon^{1}{{y}}^{1}\right)_{i}<\lambda^{1}$
$\displaystyle
1+\max_{i\in[n]}\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}\epsilon^{2}{{y}}^{2}\right)_{i}<\lambda^{2}$

Thus, for $i\leq n$,

$\displaystyle\hat{M}_{i}^{1}\left(\epsilon^{1}{{y}}^{1},z^{2*}\right)$
$\displaystyle=\frac{\left(\left(I-Z^{2*}\right)\left(D_{f}^{1}\right)^{-1}B_{f}^{1}\epsilon^{1}y^{1}\right)_{i}}{1+\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}\epsilon^{1}{y}^{1}\right)_{i}}$
$\displaystyle=\frac{\lambda^{1}\epsilon^{1}{y}_{i}^{1}}{1+\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}\epsilon^{1}{{y}}^{1}\right)_{i}}>\epsilon^{1}{y}_{i}^{1}$

Similarly, we have
$\hat{M}_{i}^{2}\left({z}^{*}_{1},\epsilon^{2}{{y}}^{2}\right)>\epsilon^{2}{y}_{i}^{2}$
for $i\leq n$. On the other hand, for $i\geq n+1$,

$\displaystyle\hat{M}_{i}^{1}\left(\epsilon^{1}{{y}}^{1},z^{2*}\right)$
$\displaystyle=\left((D_{f}^{1})^{-1}(B_{f}^{1}\epsilon^{1}{{y}}^{1})\right)_{i}$
$\displaystyle=\epsilon^{1}\lambda^{1}{y}_{i}^{1}>\epsilon^{1}{y}_{i}^{1},$

and analogously, we have
$\hat{M}_{i}^{2}\left(z^{1*},\epsilon^{2}{{y}}^{2}\right)>\epsilon^{2}{y}_{i}^{2}$
for $i\geq n+1$.

By using the monotonicity of the map ($M,\hat{M}$) and noticing that
$M\geq\hat{M}$, it follows that

${M}^{1}\left(z^{1},z^{2}\right)\geq\hat{M}^{1}\left(z^{1},z^{2}\right)>\epsilon^{1}{{y}}^{1},\quad
M^{2}\left(z^{1},z^{2}\right)\geq\hat{M}^{2}\left(z^{1},z^{2}\right)>\epsilon^{2}{{y}}^{2},$

if $\epsilon^{1}{{y}}^{1}\leq z^{1}\leq z^{1*}$ and $\epsilon^{2}{{y}}^{1}\leq
z^{2}\leq z^{2*}$. Thus, the map $M$ maps
$\\{(z^{1},z^{2})|\epsilon^{1}{{y}}^{1}\leq z^{1}\leq
z^{1*},\quad\epsilon^{2}{{y}}^{1}\leq z^{2}\leq z^{2*}\\}$ into itself.
According to the Brouwer fixed point theorem [52], we finish the proof.

######  Remark 5.12.

In [28], the authors study a very similar mean-field bi-virus planar SIS model
in a simplicial complex. Such a model can be seen as a simplified version of
our bi-virus model (23), as it does not account for the indirect spreading via
resources. In [28], it is shown that the mean-field bi-virus model also
possesses some coexisting equilibria. Further simulation results suggest that
the coexistence equilibrium is unstable. In contrast, we have provided
rigorous proof of the existence of a coexisting equilibrium in an even more
general network structure and modeling setting. Its stability remains a future
research question. However, from the numerical study performed by us, such a
coexisting equilibrium seems to always be unstable.

######  Theorem 5.13 (Irreducible monotone system).

If Assumptions 4-6 hold, the bi-virus model (23) is an
$(\mathbf{0}_{n+m},\mathbf{1}_{n+m})-$type irreducible monotone system in
$\mathbf{{D}}$. The model has a finite number of equilibria in
$\mathbf{{\bar{D}}}$ for a generic choice of parameters. The model converges
to a equilibrium for almost all initial conditions. That is to say, the set of
the initial conditions, such that the model does not converge to an
equilibrium, is a set of Lebesgue measure zero.

######  Proof 5.14.

We recall that the system (23) can be formulated as

$\begin{split}\dot{z}^{\nu}&=-D_{f}^{\nu}z^{\nu}+(I-Z^{1}-Z^{2})(B_{f}^{\nu}z^{\nu}+H^{\nu}(z^{\nu})),\nu=1,2.\end{split}$

Thus, the Jacobian is

$J=\left[\begin{matrix}J^{1}&A^{1}\\\ A^{2}&J^{2}\\\ \end{matrix}\right];$

where

$\begin{split}J^{\nu}&=-D^{\nu}_{f}+(I-{Z}^{1}-Z^{2})(B^{\nu}_{f}+[z^{\nu\top}(B^{\nu}_{f1}+B_{f1}^{\nu\top}),\cdots,z^{\nu\top}(B^{\nu}_{fn+m}+B_{fn+m}^{\nu\top})]^{\top})\\\
&\qquad-\hat{I}([z^{\nu\top}B^{\nu}_{f1}z,\cdots,z^{\nu\top}B^{\nu}_{fn+m}z^{\nu}]^{\top}+B_{f}^{\nu}z^{\nu})\end{split}$

and $A^{\nu}=-(B_{f}^{\nu}z^{\nu}+H(z^{\nu}))$. We notice that $J^{\nu}$ is an
irreducible Metzler matrix and $A^{\nu}$ is a non-positive matrix as long as
the system variables are in $\mathbf{\bar{D}}$. Let
$P=\operatorname*{\mathsf{diag}}((1,\cdots,1,-1,\cdots,-1)^{\top})$, then
$PJP$ is an irreducible Metzler matrix. Thus, the bi-virus model (23) is an
$(\mathbf{0}_{n+m},\mathbf{1}_{n+m})-$type irreducible monotone system in
$\mathbf{{D}}$.

Furthermore, by Theorem B.1 and Corollary B.2 in [64], the system (23) has a
finite number of equilibria in the closure of $\mathbf{{D}}$ if the parameters
do not lie on a certain algebraic set of measure zero because we can observe
that there are some cases where the healthy state is the unique equilibrium
and is globally stable and isolated. Alternatively, we can set parameters so
that the equation system is partially decoupled, similar to the proof of
Theorem 4.19. We further see it yields a finite number of soultions.

By Lemma 2.3 of [64], the solutions of (23) converge to an equilibrium with
some domain of attraction; By the same Lemma, the set of the initial
conditions, such that the model does not converge to an equilibrium, is a set
of measure zero.

##  6 A general SIWS process on a hypergraph

Recall that the single-virus system (10) and the bi-virus system (23) studied
so far are based on a polynomial interaction function. In reality, the
interaction functions can be any general smooth function. Thus, in this
section, we propose a generalized SIWS model, where interaction functions
remain abstract. Note also that we only deal with up to 3-body interactions
with polynomial functions in sections 3-5. The results in this section will be
also valid for a SIWS model beyond 3-body interactions with either polynomial
or non-polynomial but smooth functions.

Let us now consider the general SIWS model on a hypergraph of $n$ populations
and $m$ resources. Component-wise, we propose a system given by

(32) $\displaystyle\dot{z}_{i}$
$\displaystyle=\dot{x}_{i}=-\delta_{i}x_{i}+(1-x_{i})\left(f_{i}(z)+h_{i}(z)\right),\qquad
1\leq i\leq n;$ (33) $\displaystyle\dot{z}_{i}$
$\displaystyle=\dot{w}_{j}=-\delta^{w}_{j}w_{j}+f_{i}(x)+h_{i}(x),\qquad\qquad
n+1\leq i=j+n\leq n+m,$

where the system variable is $z=(x,w)^{\top}\in\mathbb{R}^{n+m}$,
$f_{i}(z),i=1,\cdots,n$ denotes the pairwise interaction
$f_{i}(x,w)=\sum_{j=1}^{n+m}A_{ij}f_{ij}(z_{j})$ on a conventional graph,
$h_{i}(z)$ denotes the higher-order interaction

$h_{i}(z)=\sum_{j,k=1}^{n+m}A_{ijk}f_{ijk}(z_{j},z_{k})+\sum_{j,k,l=1}^{n+m}A_{ijkl}f_{ijkl}(z_{j},z_{k},z_{l})+\cdots$

among those higher-order hyperedges on a hypergrpah. For the network
structure, we assume that the resources are independent of each other and each
resource is connected with at least one population node. For this reason, the
interaction function in $w_{j}$ is reduced to $f_{i}(x)+h_{i}(x)$, where
$f_{i},h_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$, instead of
$f_{i}(z)+h_{i}(z)$, where
$f_{i},h_{i}:\mathbb{R}^{n+m}\rightarrow\mathbb{R}$. The interpretation of
this model is that we take further heterogeneity into account and consider
different kinds of interactions (not necessarily to be a polynomial).

We can further compactify (32)-(33) as

(34) $\displaystyle\dot{x}$
$\displaystyle=-Dx+(I-\operatorname*{\mathsf{diag}}(x))\left(F_{x}(z)+H_{x}(z)\right),$
(35) $\displaystyle\dot{w}$ $\displaystyle=-D_{w}w+F_{w}(z)+H_{w}(x),$

where

(36)
$\begin{split}D&=\operatorname*{\mathsf{diag}}((\delta_{1},\cdots,\delta_{n})^{\top})\\\
D_{w}&=\operatorname*{\mathsf{diag}}((\delta^{w}_{1},\cdots,\delta^{w}_{m})^{\top})\\\
F_{x}&=(f_{1}(z),\cdots,f_{n}(z))^{\top}\\\
F_{w}&=(f_{n+1}(z),\cdots,f_{n+m}(z))^{\top}\\\
H_{x}&=(h_{1}(z),\cdots,h_{n}(z))^{\top}\\\
F_{w}&=(h_{n+1}(z),\cdots,h_{n+m}(z))^{\top}\end{split}$

We can rewrite our system (34)-(35) as

(37) $\dot{z}=-D_{f}z+(I-Z)F(z)+(I-Z)H(z);$

where

(38) $\displaystyle D_{f}:=\left[\begin{matrix}D&\mathbf{0}\\\
\mathbf{0}&D^{w}\\\
\end{matrix}\right],Z=\left[\begin{matrix}\operatorname*{\mathsf{diag}}(x)&0\\\
0&0\\\ \end{matrix}\right],F(z)=(F_{x}(z),F_{w}(z))^{\top},$ $\displaystyle
H(z)=(H_{x}(z),H_{w}(x))^{\top}.$

Now, we present the analytical results of the general SIWS on a hypergraph
(37). Firstly, we need the following assumptions to make sure our model is
well-defined.

######  Assumption 7.

[Natural bound of probability] The initial condition of infection level
satisfies $x_{i}(0)\in[0,1]$ for all $i=1,\cdots,n$.

######  Assumption 8.

[Healthy state as a fixed point] The interaction process must satisfy $F(0)=0$
and $H(0)=0$.

######  Assumption 9.

[Infected neighbours increases the chance of infection] The interaction
function satisfies $\frac{\partial f_{i}(z)}{\partial z_{j}}>0$ for all $i$
and $j$ such that there is a pairwise connection $A_{ij}>0$ (conventional
edge) and $\frac{\partial h_{i}(z)}{\partial z_{j}}\geq 0$ for all
$i=1,\cdots,n+m$ and $j=1,\cdots,n$.

######  Assumption 10.

[Higher order interaction is in higher order] The higher-oder interaction
function satisfies $\frac{\partial h_{i}(\mathbf{0})}{\partial z_{j}}=0$ for
all $i$ and $j$.

######  Assumption 11.

[Strong connectivity] The matrix $[A_{ij}]^{(n+m)\times(n+m)}$ is
irreducible.

In the following Lemma, we describe the general properties of the model (37).

######  Lemma 6.1 (General properties).

Consider the model (37). Let
$w_{j\max}=\frac{f_{j}(1)+h_{j}(1)}{\delta_{j}^{w}}$,

$\mathbf{D}=\left\\{z=(x,w)^{\top}=(x_{1},\ldots,x_{n},w_{1},\ldots,w_{m})^{\top}\,|\,x_{i}\in(0,1),\,w_{j}\in(0,w_{j\max})\right\\},$

and $\mathbf{\bar{D}}$ be its closure. If Assumptions 8-10 hold, then

  * i)

the set $\mathbf{\bar{D}}$ is positively invariant;

  * ii)

the set $\mathbf{D}$ is also positively invariant;

  * iii)

the origin is always an equilibrium of the model and there is no other
equilibrium on the boundary of $\mathbf{D}$.

######  Proof 6.2.

According to Assumptions 8-10, we have the following arguments. If $x_{i}=0$,
we have $\dot{x}_{i}=f_{i}(z)+h_{i}(z)\geq 0$. If $w_{j}=0$,
$\dot{w}_{j}=f_{j+n}(x)+h_{j+n}(x)\geq 0$. On the boundary (i.e., at least one
$z_{i}=0$), we have that $x=\mathbf{0}$ and $w=\mathbf{0}$ if and only if
$\dot{x}=\mathbf{0}$ and $\dot{w}=\mathbf{0}$.

If $x_{i}=1$, we get $\dot{x}_{i}=-\delta_{i}<0$. If $w_{j}=w_{j\max}$, we
have
$\dot{w}_{j}=-f_{n+j}(\mathbf{1})-h_{n+j}(\mathbf{1})+f_{n+j}(x)+h_{n+j}(x)\leq
0$, since the pairwise interaction is increasing in $x$ and the higher-order
interaction function is non-decreasing in $x$. If $w_{j}=w_{j\max}$ and
$x\neq\mathbf{1}$, $\dot{w}_{j}<0$.

All these ensure that statements i) ii) and iii) must be true.

######  Remark 6.3.

These general properties show that the system is well-defined. The set
$\mathbf{\bar{D}}$ is the system domain of the model (37).

It is obvious that the system has an all-zero equilibrium, which is the so-
called healthy state. Let us define the notation $J_{F}$, which denotes the
Jacobian of the corresponding function $F$.

######  Theorem 6.4 (Healthy-state behaviour).

Consider the model (37), if Assumptions 7-11 hold, we have the following
statements. The healthy state (zero equilibrium)

  * 1)

is locally stable if $\rho(D^{-1}_{f}J_{F}(0))<1$;

  * 2)

is globally exponentially stable if $D^{-1}_{f}F(z)+D^{-1}_{f}H(z)\leq dz$ for
some $0<d<1$;

  * 3)

is unstable if $\rho(D^{-1}_{f}J_{F}(0))>1$.

######  Proof 6.5.

The Jacobian matrix of the model at the origin reads

$J(0,0)=-D_{f}+J_{F}(0).$

Since the matrix $J_{F}(0)$ is irreducible nonnegative from Assumption 11 and
$-D_{f}$ is a negative diagonal matrix, $\rho(D^{-1}_{f}J_{F}(0))>1$ is
equivalent to $s(-D_{f}+B_{f})>0$ and $\rho(D^{-1}_{f}B_{f})<1$ is equivalent
to $s(-D_{f}+B_{f})<0$ by Proposition 1 of [32]. This ensures that the origin
is locally stable if $\rho(D^{-1}_{f}J_{F}(0))<1$ and it is unstable if
$\rho(D^{-1}_{f}J_{F}(0))>1$.

Next, we show global stability. Define the Lyapunov function $V=D_{f}^{-1}z$.
We observe that $V\geq\mathbf{0}$ and $V=\mathbf{0}$ if and only if
$z=\mathbf{0}$.

Then, for some $0<d<1$, we get

(39)
$\begin{split}\dot{V}&=D_{f}^{-1}\left(-D_{f}z+(I-Z)F(z)+(I-Z)H(z)\right)\\\
&\leq-z+D_{f}^{-1}(F(z)+H(z))\\\ &\leq-z+D_{f}^{-1}dz=(-1+d)z\\\
&=(-1+d)D_{f}D_{f}^{-1}z=(-1+d)\min(\delta,\delta^{w})V,\end{split}$

The rest of the proof is similar to the proof of Theorem 4.8.

The condition $D^{-1}_{f}F(z)+D^{-1}_{f}H(z)\leq dz$ generally requires that
the infection interaction process doesn’t increase sharply and is upper
bounded by a line with the slope smaller than $1$. Next, we deal with the
existence of an endemic equilibrium.

######  Theorem 6.6 (Existence of an endemic equilibrium).

Consider the system (37), if Assumptions 7-11 hold, there exists an endemic
equilibrium $z^{*}\gg\mathbf{0}$ such that $x_{i}^{*}\geq\frac{k-1}{k}$ and
$w^{*}_{i}\geq\frac{(k-1)w_{i\max}}{k}$ for any $i$ such that
$B_{fi}\neq\mathbf{0}$ if
$\theta=\min_{i}(D^{-1}_{f}F(\frac{k-1}{k}\tilde{u})+D^{-1}_{f}H(\frac{k-1}{k}\tilde{u}))_{i}\geq\frac{(k-1)^{2}}{k}\tilde{u}$,
where

$\tilde{u}_{i}=\begin{cases}1,&H(z)_{i}\not\equiv 0,\qquad i=1,\ldots,n\\\
w_{i\max},&H(z)_{i}\not\equiv 0,\qquad i=n+1,\ldots,n+m\\\ 0,&H(z)_{i}\equiv
0\end{cases}$

######  Proof 6.7.

We use the same map (17) and it follows the same monotone property. It is easy
to check that the fixed point of the map (17) is also the equilibrium of (37).
Now, we define the set
$Y=\\{z=(x,w)^{\top}\in\mathbb{R}^{n+m}|\frac{1}{2}\tilde{u}\leq z\leq u\\}$.

We have that

$\begin{split}\hat{T}(z)&=T(D^{-1}_{f}F(z)+D^{-1}_{f}H(z))\\\ &\geq
T\left(D^{-1}_{f}F(\frac{1}{2}\tilde{u})+D^{-1}_{f}H(\frac{1}{2}\tilde{u})\right)\\\
&\geq T\left(\theta\right)\geq\frac{k-1}{k}\tilde{u}\end{split}$

Notice that $\frac{p_{i}}{1+p_{i}}\leq 1$ for all $i\leq n$ and $p_{j+n}\leq
w_{j\max}$ for all $j\leq m$. Since $\hat{T}(z)$ is a continuous map and it
maps the set $Y$ to itself. The Brouwer Fixed-Point Theorem [52] ensures that
there exists at least one $z^{*}\gg\mathbf{0}$ of the given form.

######  Theorem 6.8 (Further results on the existence of an endemic
equilibrium).

Consider the model (37), if Assumptions 7-11 hold and there exists a $c\neq 0$
such that $c\ll u$ and $D^{-1}_{f}F(c)\geq c$, then there exists an endemic
equilibrium $z^{*}\gg\mathbf{0}$.

######  Proof 6.9.

Define the set $\tilde{Y}=\\{z=(x,w)^{\top}|c_{i}\leq x_{i}\leq 1,c_{j+n}\leq
w_{j}\leq w_{j\max}\\}$, where $c=(c_{1},\cdots,c_{m+n})^{\top}$.

We use again the map (17). For any $z\in\tilde{Y}$, we have
$\hat{T}(z)=T(D^{-1}_{f}F(z)+D^{-1}_{f}H(z))\geq T(D^{-1}_{f}F(z))\geq
T(D^{-1}_{f}F(c))\geq c$. We use again the Brouwer Fixed-Point Theorem [52]
and confirm that there exists an endemic equilibrium $z^{*}\gg\mathbf{0}$ if
$D^{-1}_{f}F(c)\geq c$.

######  Remark 6.10.

The condition $D^{-1}_{f}F(c)\geq c$ is not strict. We consider $c$ now very
close to the all-zero vector. We expand the function $F$ around $0$ and get
$D_{f}^{-1}(F(0)+J_{F}(0)c)\geq c$. Since $F(0)=0$, we finally get
$J_{F}(0)c\geq D_{f}c$, which implies the influence of infection rates around
zero is stronger than the influence of healing rates.

Next, we can further show that the general system (37) is also an irreducible
monotone system.

######  Theorem 6.11 (Irreducible monotone system).

If Assumption 7-11 hold,
the model (37) is an irreducible monotone system in $\mathbf{{D}}$.
Furthermore, if the model (37) has a finite number of equilibria in the
closure of $\mathbf{\bar{D}}$ for a generic choice of parameters, then the
model converges to an equilibrium for almost all initial conditions. That is,
the set of initial conditions, such that the model does not converge to an
equilibrium, is a set of Lebesgue measure zero.

######  Proof 6.12.

Firstly, we calculate the Jacobian

$\displaystyle
J=-D_{f}+(I-Z)(J_{F}(z)+J_{H}(z))-Z\operatorname*{\mathsf{diag}}(F(z)+H(z));$

Note that $J_{F}(z)+J_{H}(z)>\mathbf{0}$. Furthermore, $D_{f}$ is diagonal.
Thus, the Jacobian is an irreducible Metzler matrix at an arbitrary point in
$\mathbf{{D}}$ along the trajectories of (37). This confirms that the system
(37) is an irreducible monotone system by Lemma 2.2 of [64].

By assumption, the equilibrium set is finite. Thus, by Lemma 2.3 of [64], the
model converges to an equilibrium with some domain of attraction. The set of
the initial conditions, such that the model does not converge to an
equilibrium, is a set of measure zero.

######  Remark 6.13.

If Assumption 7-10 hold, and furthermore the pairwise interaction and the
higher order interaction is a polynomial function, the model (37) under this
setting has a finite number of equilibria in the closure of $\mathbf{\bar{D}}$
for a generic choice of parameter, then the model converges to an equilibrium
for almost all initial conditions. That is, the set of the initial conditions,
such that the model does not converge to an equilibrium, is a set of Lebesgue
measure zero.

Furthermore, if $\rho(D^{-1}_{f}B_{f})>1$, the system with pairwise
interaction is similar to the conventional system (10) and without the higher-
order interaction, the system is reduced to the conventional continuous-time
SIWS model on a graph [41]. Now, the system on a graph reads as

(40) $\dot{z}=-D_{f}z+(I-{Z})B_{f}z.$

######  Proposition 1.

Consider the system (40) with $D_{f}$ positive diagonal and $B_{f}$
irreducible non-negative. If $\rho(D^{-1}_{f}B_{f})>1$, then the system has a
unique endemic equilibrium (strictly positive equilibrium) and it is
hyperbolic and globally stable in $\bar{\mathbf{D}}\setminus\\{\mathbf{0}\\}$.

######  Proof 6.14.

The existence, uniqueness, and global stability of the endemic equilibrium are
reported in theorem 2 in [41]. So, we only show that the endemic equilibrium
is hyperbolic, which means that the Jacobian at the endemic equilibrium has no
eigenvalues with a zero real part.

So we can see that the Jacobian is $J=-D_{f}+B_{f}-{Z^{*}}B_{f}-B_{f}{Z^{*}}$.
We notice that
$Jz^{*}=-D_{f}z^{*}+B_{f}z^{*}-{Z^{*}}B_{f}z^{*}-B_{f}{Z^{*}}z^{*}=-B_{f}{Z^{*}}z^{*}<-dz^{*}$
with some positive $d$. Thus, $J$ is Hurwitz and contains only eigenvalues
with a negative real part. Thus, the endemic equilibrium is hyperbolic.

Let us now consider the system (40) with $D_{f}$ positive diagonal and $B_{f}$
irreducible non-negative with some small higher-order term as $\epsilon H(z)$
satisfying Assumption 8-10. That is,

(41) $\dot{z}=-D_{f}z+(I-{Z})(B_{f}z+\epsilon H(z)).$

######  Theorem 6.15.

Consider the perturbed system (41). If Assumptions 7-11 hold, and
$\rho(D^{-1}_{f}B_{f})>1$, then there is at least one endemic equilibrium,
which is hyperbolic and locally stable, for $0<\epsilon\ll 1$ sufficiently
small.

######  Proof 6.16.

We know from [41] that if $\rho(D^{-1}_{f}B_{f})>1$, the unperturbed system
has a unique endemic equilibrium and it is hyperbolic. We can rewrite the
unperturbed system as $g(z)$ and the perturbed system as $G(z,\epsilon)$. Let
$z^{*}$ be the unique endemic equilibrium of the unperturbed system. By
definition, $G(z^{*},0)=0$ and $\frac{\partial G}{\partial
x}(z^{*},0)=\frac{\partial f}{\partial x}(z^{*})$. Since $z^{*}$ is
hyperbolic, $\frac{\partial G}{\partial x}(z^{*},0)=\frac{\partial f}{\partial
x}(z^{*})$ has a nonvanishing determinant. By the implicit function theorem,
there is a unique equilibrium in the neighborhood of $z^{*}$ for sufficiently
small $\epsilon$. This equilibrium is also hyperbolic by continuous dependence
of the eigenvalues of $\frac{\partial G}{\partial x}$ on $\epsilon$. Thus, the
local stability of the equilibrium also persists.

The model (37) can be easily extended into a competing bi-virus SIWS model,
which reads as

(42) $\displaystyle\dot{z}^{1}_{i}$
$\displaystyle=\dot{x}^{1}_{i}=-\delta^{1}_{i}x^{1}_{i}+(1-x^{1}_{i}-x^{2}_{i})\left(f^{1}_{i}(z^{1})+h^{1}_{i}(z^{1})\right),$
$\displaystyle\quad 1\leq i\leq n,$ (43) $\displaystyle\dot{z}^{1}_{i}$
$\displaystyle=\dot{w}^{1}_{j}=-\delta^{w1}_{j}w^{1}_{j}+f^{1}_{i}(x^{1})+h^{1}_{i}(x^{1}),$
$\displaystyle\quad n+1\leq i=j+n\leq n+m,$ (44)
$\displaystyle\dot{z}^{2}_{i}$
$\displaystyle=\dot{x}^{2}_{i}=-\delta^{2}_{i}x^{2}_{i}+(1-x^{1}_{i}-x^{2}_{i})\left(f^{2}_{i}(z^{2})+h^{2}_{i}(z^{2})\right),$
$\displaystyle\quad 1\leq i\leq n$ (45) $\displaystyle\dot{z}^{2}_{i}$
$\displaystyle=\dot{w}^{1}_{j}=-\delta^{w2}_{j}w^{2}_{j}+f^{2}_{i}(x^{2})+h^{2}_{i}(x^{2}),$
$\displaystyle\quad n+1\leq i=j+n\leq n+m,$

where all notations remain the same as in (36), (38) and the superscript $1,2$
denotes the first or the second virus.

For $\nu=1,2$, we can rewrite (42)-(45) as

(46)
$\begin{split}\dot{z}^{\nu}&=-D_{f}^{\nu}z^{\nu}+(I-Z^{1}-Z^{2})(F^{\nu}(z^{\nu})+H^{\nu}(z^{\nu})).\end{split}$

Throughout the rest of this paper, we assume that if Assumption 8-10 hold for
the system (46), then Assumptions 8-10 hold for each
$F^{\nu}(z^{\nu}),H^{i}(z^{\nu}),A^{\nu},\nu=1,2$.

For the bi-virus system, we have a new assumption on the initial condition.

######  Assumption 12.

The initial condition of infection level satisfies $x^{\nu}_{i}(0)\in[0,1]$
for all $i$ and $\nu$. Furthermore, $x^{1}_{i}(0)+x^{2}_{i}(0)\leq 1$ holds
for all $i=1,\cdots,n$.

######  Lemma 6.17 (General properties).

Consider the bi-virus model (46) and let
$w^{\nu}_{j\max}=\frac{f^{\nu}_{j}(1)+h^{\nu}_{j}(1)}{\delta_{j}^{w\nu}}$,

(47)
$\begin{split}\mathbf{D}&=\left\\{z=(x^{1},w^{1},x^{2},w^{2})^{\top}=(x^{1}_{1},\cdots,x^{1}_{n},w^{1}_{1},\cdots,w^{1}_{m},x^{2}_{1},\cdots,x^{2}_{n},w^{2}_{1},\cdots,w^{2}_{m})^{\top}|\right.\\\
&\qquad\left.x^{\nu}_{i}\in(0,1),w^{\nu}_{j}\in(0,w^{\nu}_{j\max}),\nu=1,2\right\\},\end{split}$

and $\mathbf{\bar{D}}$ be its closure. If Assumptions 8-12 hold, then

  * i)

the set $\mathbf{\bar{D}}$ is positively invariant;

  * ii)

the set $\mathbf{D}$ is also positively invariant;

  * iii)

the origin is always an equilibrium of the model;

  * iv)

the dynamics of each virus (46) are upper bounded by its single virus
counterpart (37) with the same system parameters.

######  Proof 6.18.

If $x^{\nu}_{i}=0$, we have
$\dot{x}^{\nu}_{i}=(1-x^{3-\nu}_{i})(f^{\nu}_{i}(z)+h^{\nu}_{i}(z))\geq 0$. If
$w^{\nu}_{j}=0$, $\dot{w}^{\nu}_{j}=f^{\nu}_{j+n}(x)+h^{\nu}_{j+n}(x)\geq 0$.
If $x^{\nu}=\mathbf{0}$ and $w^{\nu}=\mathbf{0}$, we have $\dot{x}=\mathbf{0}$
and $\dot{w}=\mathbf{0}$, which shows that the origin is always an
equilibrium.

If $x^{\nu}_{i}=1$, we get $\dot{x}^{\nu}_{i}=-\delta_{i}<0$. If
$w_{j}^{\nu}=w^{\nu}_{j\max}$, we have
$\dot{w}^{\nu}_{j}=-f^{\nu}_{n+j}(\mathbf{1})-h^{\nu}_{n+j}(\mathbf{1})+f^{\nu}_{n+j}(x^{\nu})+h^{\nu}_{n+j}(x^{k})\leq
0$, since the pairwise interaction is increasing in $x$ and the higher-order
interaction function is non-decreasing in $x$. If
$w^{\nu}_{j}=w^{\nu}_{j\max}$ and $x^{\nu}\neq\mathbf{1}$, then
$\dot{w}^{\nu}_{j}<0$.

If $x^{1}_{i}+x^{2}_{i}=1$, then $x^{\nu}_{i}=-\delta_{i}x^{\nu}_{i}<0$. All
these arguments show that statements i) and ii) are true.

The statement iv) is straightforward. It suffices to compare the dynamics of
each virus (46) with its single virus counterpart with the same system
parameters (37).

######  Theorem 6.19 (Healthy-state behaviour).

Consider the model (46), if Assumptions 8-12 hold, we have the following
statements. The healthy state (zero equilibrium)

  * 1)

is locally stable if $\rho((D^{\nu}_{f})^{-1}J_{F^{\nu}}(0))<1$ for $\nu=1,2$;

  * 2)

is globally exponentially stable if
$(D^{\nu}_{f})^{-1}F^{\nu}(z^{\nu})+(D^{\nu}_{f})^{-1}H^{\nu}(z^{\nu})\leq
dz^{\nu}$, $\nu=1,2$ for some $0<d<1$;

  * 3)

is unstable if $\rho((D^{\nu}_{f})^{-1}J_{F^{\nu}}(0))>1$ for either $\nu=1$
or $\nu=2$.

######  Proof 6.20.

It is already clear that the healthy state (zero equilibrium) always exists.
The Jacobian matrix of the model at the origin reads as

$J(0,0)=\left[\begin{matrix}-D^{1}_{f}+J_{F^{1}}(0)&\mathbf{0}\\\
\mathbf{0}&-D^{2}_{f}+J_{F^{2}}(0)\\\ \end{matrix}\right].$

Thus, the healthy state is locally stable if
$\rho((D^{\nu}_{f})^{-1}J_{F^{\nu}}(0))<1$ for $\nu=1,2$ and unstable if
$\rho((D^{\nu}_{f})^{-1}J_{F^{\nu}}(0))>1$ for either $\nu=1$ or $\nu=2$.

The global exponential stability when
$(D^{\nu}_{f})^{-1}F^{\nu}(z^{\nu})+(D^{\nu}_{f})^{-1}H^{\nu}(z^{\nu})\leq
dz^{\nu}$ for some $0<d<1$ and $\nu=1,2$ is the direct consequence of
statement iv) in Lemma 6.17 and Theorem 6.4.

######  Theorem 6.21 (Dominant endemic equilibrium).

Consider the bi-virus model
(46). If Assumptions 8-12 hold, we have the following statements. If
$z^{*}=(x^{*},w^{*})^{\top}$ is the endemic equilibrium of the single-virus
counterpart (10) of the first virus, then $(z^{*},0,0)^{\top}$ is the dominant
endemic equilibrium of the system (46). The dominant endemic equilibrium is
locally stable when $z^{*}$ is locally stable and
$\rho(-D^{2}_{f}+(I-Z^{*})J_{F^{2}}(0))<0$ with respect to the single-virus
system and is unstable when $z^{*}$ is unstable or
$\rho(-D^{2}_{f}+(I-Z^{*})J_{F^{2}}(0))>0$, where
$Z^{*}=\left[\begin{matrix}x^{*}&0\\\ 0&0\\\ \end{matrix}\right]$.

######  Proof 6.22.

If we set $z^{2}=\mathbf{0}$, the system is reduced into the single virus
counterpart of the first virus. Thus, if $z^{*}$ is the endemic equilibrium of
the single-virus counterpart (10) of the first virus, then
$(z^{*},0,0)^{\top}$ is the dominant endemic equilibrium of the system (46).

The Jacobian matrix of the model at the dominant endemic equilibrium is

$J(z^{*},\mathbf{0})=\left[\begin{matrix}\hat{J}(z^{*})&\mathbf{Res}\\\
\mathbf{0}&-D^{2}_{f}+(I-Z^{*})J_{F^{2}}(0)\\\ \end{matrix}\right],$

where $\hat{J}(z^{*})$ is the Jacobian of the single-virus counterpart at the
endemic equilibrium $z^{*}$ and $\mathbf{Res}$ is the matrix of some residual
which we do not need to deal with.

By observing the Jacobian and noticing that $\hat{J}(z^{*})$ and
$-D^{2}_{f}+(I-Z^{*})J_{F^{2}}(0)$ is irreducible nonnegative, we can directly
obtain the results.

######  Theorem 6.23 (Coexisting equilibrium).

Consider the model (46). If Assumptions 8-12 hold, we have the following
statements. If there exist $(D^{\nu}_{f})^{-1}F^{\nu}(c)\geq c$ and $c\neq 0$,
$c\ll u^{i}$, $F^{\nu}(z^{\nu})\geq B^{\nu}_{f}z^{\nu}$ for both viruses
$\nu=1,2$, and $\rho(-D^{1}_{f}+(I-Z^{2*})B^{1}_{f})>0$ as well as
$\rho(-D^{2}_{f}+(I-Z^{1*})B^{2}_{f})>0$, where $z^{1*}$ and $z^{2*}$ are the
endemic equilibrium of the single-virus counterpart of the first and second
virus respectively and always exist, then there exists at least one coexisting
equilibrium $(\hat{z}^{1},\hat{z}^{2})\gg\mathbf{0}$ such that
$\hat{x}^{1}+\hat{x}^{2}\leq\mathbf{1}$.

######  Proof 6.24.

First of all, we need to confirm that $z^{1*}$ and $z^{2*}$, the endemic
equilibrium of the single-virus counterpart of the first and second virus,
always exist under the conditions. Since there must exist
$(D^{\nu}_{f})^{-1}F^{\nu}(c)\geq c$ for both viruses $\nu=1,2$, according to
the Theorem 6.8, $z^{1*}$ and $z^{2*}$ must always exist.

Now, define the map $M$ as follows. Component-wise, it is given by

(48) $\displaystyle M_{i}^{1}\left(z^{1},z^{2}\right)$
$\displaystyle=\frac{\left(1-z_{i}^{2}\right)\left(\left(D_{f}^{1}\right)^{-1}F^{1}(z^{1})+\left(D_{f}^{1}\right)^{-1}H^{1}(z^{1})\right)_{i}}{1+\left(\left(D_{f}^{1}\right)^{-1}F^{1}(z^{1})+\left(D_{f}^{1}\right)^{-1}H^{1}(z^{1})\right)_{i}}$
$\displaystyle M_{i}^{2}\left(z^{1},z^{2}\right)$
$\displaystyle=\frac{\left(1-z_{i}^{1}\right)\left(\left(D_{f}^{2}\right)^{-1}F^{2}(z^{2})+\left(D_{f}^{2}\right)^{-1}H^{2}(z^{2})\right)_{i}}{1+\left(\left(D_{f}^{2}\right)^{-1}F^{2}(z^{2})+\left(D_{f}^{2}\right)^{-1}H^{2}(z^{2})\right)_{i}}$

for $i\leq n$. Furthermore, for $i\geq n+1$,

(49) $\displaystyle M_{i}^{1}\left(z^{1},z^{2}\right)$
$\displaystyle=(D_{f}^{1})^{-1}(F^{1}(z^{1})+H^{1}(z^{1}))$ $\displaystyle
M_{i}^{2}\left(z^{1},z^{2}\right)$
$\displaystyle=(D_{f}^{2})^{-1}(F^{2}(z^{2})+H^{2}(z^{2})).$

Similarly, by using the property of the monotonicity of $M$, we get
$M(z^{1},z^{2})\leq M(z^{1*},z^{2*})$ if $(z^{1},z^{2})\leq(z^{1*},z^{2*})$.

Now, we construct another map $\hat{M}$. Component-wise, it is given by, for
$i\leq n$,

(50) $\displaystyle\hat{M}_{i}^{1}\left(z^{1},z^{2}\right)$
$\displaystyle=\frac{\left(1-z_{i}^{2}\right)\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}z^{1}\right)_{i}}{1+\left(\left(D_{f}^{1}\right)^{-1}B_{f}^{1}z^{1}\right)_{i}}\leq
M_{i}^{1}\left(z^{1},z^{2}\right)$
$\displaystyle\hat{M}_{i}^{2}\left(z^{1},z^{2}\right)$
$\displaystyle=\frac{\left(1-z_{i}^{1}\right)\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}z^{2}\right)_{i}}{1+\left(\left(D_{f}^{2}\right)^{-1}B_{f}^{2}z^{2}\right)_{i}}\leq
M_{i}^{2}\left(z^{1},z^{2}\right);$

these inequalities are true because the function $\frac{s}{1+s}$ is increasing
as $s$ increases, and furthermore $F^{\nu}(z^{\nu})+H^{\nu}(z^{\nu})\geq
F^{\nu}(z^{\nu})\geq B_{f}^{\nu}z^{\nu}$.

Furthermore, for $i\geq n+1$,

(51) $\displaystyle\hat{M}_{i}^{1}\left(z^{1},z^{2}\right)$
$\displaystyle=(D_{f}^{1})^{-1}(B_{f}^{1}z^{1})\leq
M_{i}^{1}\left(z^{1},z^{2}\right)$
$\displaystyle\hat{M}_{i}^{2}\left(z^{1},z^{2}\right)$
$\displaystyle=(D_{f}^{2})^{-1}(B_{f}^{2}z^{2})\leq
M_{i}^{2}\left(z^{1},z^{2}\right).$

The rest proof remains the same as for Theorem 5.10.

######  Remark 6.25.

The condition $F^{\nu}(z^{\nu})\geq B^{\nu}_{f}z^{\nu}$ shows that the
pairwise interaction is lower bounded by a linear function. Thus, the rate of
the increase of the pairwise interaction must exceed the slope of the linear
function. Alternatively, the general higher-order model can be seen as lower-
bounded by the conventional SIWS model on a graph, since the linear function
denotes the infection process of the conventional SIWS model on a graph.

######  Theorem 6.26 (Irreducible monotone system).

If Assumptions 8-12 hold, the model (46) is an
$(\mathbf{0}_{n+m},\mathbf{1}_{n+m})-$type irreducible monotone system in
$\mathbf{{D}}$. Furthermore, if the model (46) has a finite number of
equilibria in the closure of $\mathbf{\bar{D}}$ for a generic choice of
parameter, then the model converges to an equilibrium for almost all initial
conditions. That is to say, the set of initial conditions, such that the model
does not converge to an equilibrium, is a set of Lebesgue measure zero.

######  Proof 6.27.

The proof is similar to the Theorem 5.13. Notice that the Jacobian can be
permutated into an irreducible Metzler matrix via the permutation matrix
$P=\operatorname*{\mathsf{diag}}((1,\cdots,1,-1,\cdots,-1)^{\top})$. The rest
of the proof remains the same.

######  Remark 6.28.

If the pairwise and the higher-order interactions are both polynomial
functions, then the model (46) has a finite number of equilibria in the
closure of $\mathbf{\bar{D}}$ for a generic choice of parameter. As a
consequence, the model converges to an equilibrium for almost all initial
conditions.

##  7 Numerical examples

In this section, we carry out numerical simulations to illustrate our
analytical results.

As for the simulation setup, we consider a hypergraph with 5 population nodes
and 2 resource nodes. The recovery rate and the decay rate are randomly picked
up from $[0,1]$. Furthermore, the infection rate and the contamination rate
are randomly picked up from $[0,0.2]$. By this setting, the Assumptions 1-3
hold already. We observe that the parameters of this setting usually
reasonably provide a reproduction number. Otherwise, for example, if we allow
the maximum infection rate and the maximum contamination rate to exceed $0.2$,
then it is very likely to observe a super big reproduction number. For
simplicity, each figure of the simulation represents the mean value of the
infection (or contamination) level of all population (or resource) nodes. In
each figure, the short name _IC_ denotes initial condition and _V_ means
virus.

###  7.1 Conventional higher-order single- and bi-virus system

Firstly, we make a simulation on the single-virus system (10) and the bi-virus
system (23). The simulation results of figure 2-8 are in line with the
analytical results of this paper. Especially, we see that bi-stability occurs
for both single-virus and bi-virus systems. Although we confirm the existence
of coexisting equilibrium under some conditions, the simulation results
suggest that it is unstable.

Figure 2: The simulation of the single-virus case when the reproduction number
equals $0.0644$. From 3 different random initial conditions, the model always
converges to zero. Figure 3: The simulation of the single-virus case when the
reproduction number equals $0.4578$. From 3 different random initial
conditions, the model either converges to the origin or to an endemic
equilibrium. Figure 4: The simulation of the single-virus case when the
reproduction number equals $1.7052$. From 3 different random initial
conditions, the model always converges to the same endemic equilibrium. Figure
5: The simulation of the bi-virus case when the reproduction numbers equal
$0.3632$ and $0.4563$ respectively. From 3 different random initial
conditions, the model always converges to the healthy state. Figure 6: The
simulation of the bi-virus case when the reproduction numbers equal $2.3867$
and $0.5647$ respectively. From 3 different random initial conditions, the
model always converges to the same dominant endemic equilibrium. Figure 7: The
simulation of the bi-virus case when the reproduction numbers equal $0.6640$
and $0.4685$ respectively. From 3 different random initial conditions, the
model either converges to a dominant endemic equilibrium or to the healthy
state. Figure 8: The simulation of the bi-virus case when the reproduction
numbers are the same and equal $1.4575$. From 3 different random initial
conditions, the model always converges to a dominant endemic equilibrium.

###  7.2 Bi-virus system with HOIs up to 4-body interactions

Here, we simulate the case when we consider the HOIs up to 4-body
interactions. That is to consider bi-virus system (46) with pairwise
interactions $F^{l}(z)=B^{l}_{f}z$ and HOIs
$h^{l}_{i}(z)=\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta^{l}_{i3}A^{l}_{ijk}z^{l}_{j}z^{l}_{k}+\sum_{j,k,l\in\mathbf{N}_{4}^{\text{i}}}\beta_{i4}A_{ijkl}z_{j}z_{k}z_{l},l=1,2$.
From figure 9 and 10, we see the system still processes bistability. These
numerical results support our analytical results of a generalized SIWS model.
Moreover, the numerical examples further suggest that any coexisting
equilibrium is unstable since we don’t observe any solution converging to any
coexisting equilibrium by calibrating system parameters and initial
conditions.

Figure 9: The simulation of the single-virus case. From 3 different random
initial conditions, the model either converges to an endemic equilibrium or to
the healthy state. Figure 10: The simulation of the bi-virus case. From 3
different random initial conditions, the model converges either to a dominant
endemic equilibrium or to the healthy state.

###  7.3 Bi-virus system with non-linear pairwise interactions and polynomial
HOIs

In this subsection, we simulate the case when we consider pairwise
interactions are non-linear and the HOIs are polynomial functions. That is to
consider bi-virus system (46) with pairwise interactions
$f^{l}_{i}(z)=\sum_{j,k\in\mathbf{N}_{2}^{\text{i}}}\beta^{l}_{i2}A^{l}_{ij}\log(1+0.2z^{l}_{j})$
and HOIs
$h^{l}_{i}(z)=\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta^{l}_{i3}A_{ijk}z^{l}_{j}z^{l}_{k},l=1,2$.
These pairwise interactions are the same as the form suggested by [10] and the
HOIs are the same form as the conventional single-virus system (10) and the
bi-virus system (23). From figure 11 and 12, we see the system still processes
bistability. These numerical results again support our analytical results of a
generalized SIWS model. Furthermore, by calibrating system parameters and
initial conditions, we don’t observe that any solution converges to any
coexisting equilibrium. Coexisting equilibrium seems always unstable.

Figure 11: The simulation of the single-virus case. From 3 different random
initial conditions, the model either converges to an endemic equilibrium or to
the healthy state. Figure 12: The simulation of the bi-virus case. From 3
different random initial conditions, the model converges either to a dominant
endemic equilibrium or to the healthy state.

###  7.4 Bi-virus system with non-linear pairwise interactions and non-
polynomial HOIs

In this subsection, we simulate the case when we consider pairwise
interactions are non-linear and the HOIs are in a similar form. That is to
consider bi-virus system (46) with pairwise interactions
$f^{l}_{i}(z)=\sum_{j,k\in\mathbf{N}_{2}^{\text{i}}}\beta^{l}_{i2}A^{l}_{ij}\log(1+0.2z^{l}_{j})$
and HOIs
$h^{l}_{i}(z)=\sum_{j,k\in\mathbf{N}_{3}^{\text{i}}}\beta^{l}_{i3}A^{l}_{ijk}z^{l}_{j}\log(1+0.2z^{l}_{k}),l=1,2$.
From figure 13 and 14, we see the system processes bistability. These
numerical results are still in line with our analytical results of a
generalized SIWS model. Furthermore, by calibrating system parameters and
initial conditions, we don’t observe that any solution converges to any
coexisting equilibrium. Coexisting equilibrium seems always unstable.

Figure 13: The simulation of the single-virus case. From 3 different random
initial conditions, the model either converges to an endemic equilibrium or to
the healthy state. Figure 14: The simulation of the bi-virus case. From 3
different random initial conditions, the model converges either to a dominant
endemic equilibrium or to the healthy state.

##  8 Conclusion and future works

In this paper, we study a general SIS-type diffusion process with both
indirect and direct pathways on a directed hypergraph. Based on a polynomial
interaction function, although the interaction among agents is firstly based
on a simplicial complex structure, we show that such a model is eventually
equivalent to a model on a general hypergraph if we adopt a spreading
mechanism of different orders. From a formal level, this has allowed us to
study a spreading process on a rather general hypernetwork setting. We further
give a full analysis of our proposed model. The system analysis mainly focuses
on the healthy-state and the endemic behaviors. Importantly, induced by the
higher order interaction, we find out that bistability may even occur when the
reproduction number is smaller than $1$. When the reproduction number is
greater than $1$, if the higher-order terms are sufficiently small, the unique
endemic equilibrium is globally stable. Then, we continue to propose a
competing bi-virus model on a directed hypergraph and the system analysis is
analogously carried out. The bistability may also occur in the bi-virus system
which is a direct consequence of the bistability of a single-virus system.
Particularly, we give a condition for the system having a coexisting
equilibrium. However, further numerical studies show that coexisting
equilibrium seems always unstable. To further generalize the proposed
framework, we consider an abstract interaction function, which in real-life
processes may account for interactions that do not follow a polynomial rule.
The corresponding single-virus and bi-virus models are proposed and studied.
The global system behavior is generally similar to the counterpart with a
polynomial interaction function. In addition, we show that all proposed
systems in this paper belong to irreducible monotone systems. Under the
further assumption of a finite number of equilibria, solutions of these
systems always converge to an equilibrium. The theoretical results are finally
supported by some numerical examples.

Several future research directions stem from the work presented here; these
include, but are not limited to, proving the instability of a coexisting
equilibrium, providing further results on the endemic behavior (stability) of
the general system, and designing control strategies to regulate the spreading
process (e.g. adding or eliminating a node or a hyperedge from a hypergraph).
For the last item, we may need to investigate tensor spectral theory [5, 47],
since a hypergraph can be abstracted as tensors, and manipulation of such a
hypergraph may change the eigenvalues of the tensor, which may further
influence the system’s behavior.

## References

  * [1] M. A. Achterberg, B. Prasse, L. Ma, S. Trajanovski, M. Kitsak, and P. Van Mieghem, Comparing the accuracy of several network-based covid-19 prediction algorithms, International journal of forecasting, 38 (2022), pp. 489–504.
  * [2] C. Bick, E. Gross, H. Harrington, and M. Schaub, What are higher order networks?, SIAM Review, (2022).
  * [3] F. Bullo, Lectures on Network Systems, Kindle Direct Publishing, 1.6 ed., 2022, http://motion.me.ucsb.edu/book-lns.
  * [4] R. Carli, G. Cavone, N. Epicoco, P. Scarabaggio, and M. Dotoli, Model predictive control to mitigate the covid-19 outbreak in a multi-region scenario, Annual Reviews in Control, 50 (2020), pp. 373–393.
  * [5] K. Chang, L. Qi, and T. Zhang, A survey on the spectral theory of nonnegative tensors, Numerical Linear Algebra with Applications, 20 (2013), pp. 891–912.
  * [6] S. S. Chaves, P. Gargiullo, J. X. Zhang, R. Civen, D. Guris, L. Mascola, and J. F. Seward, Loss of vaccine-induced immunity to varicella over time, New England Journal of Medicine, 356 (2007), pp. 1121–1129.
  * [7] P. Cisneros-Velarde and F. Bullo, Multi-group sis epidemics with simplicial and higher-order interactions, IEEE Transactions on Control of Network Systems, (2021).
  * [8] S. Cui, F. Liu, H. Jardón-Kojakhmetov, and M. Cao, Discrete-time layered-network epidemics model with time-varying transition rates and multiple resources, arXiv preprint arXiv:2206.07425, (2022).
  * [9] G. F. de Arruda, G. Petri, and Y. Moreno, Social contagion models on hypergraphs, Physical Review Research, 2 (2020), p. 023032.
  * [10] 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).
  * [11] G. Gallo, G. Longo, S. Pallottino, and S. Nguyen, Directed hypergraphs and applications, Discrete applied mathematics, 42 (1993), pp. 177–201.
  * [12] L. Gallo, R. Muolo, L. V. Gambuzza, V. Latora, M. Frasca, and T. Carletti, Synchronization induced by directed higher-order interactions, Communications Physics, 5 (2022), p. 263.
  * [13] S. Gracy, I. C. Morarescu, V. S. Varma, and P. E. Pare, Analysis and on/off lockdown control for time-varying sis epidemics with a shared resource, in 2022 European Control Conference (ECC), IEEE, 2022, pp. 1660–1665.
  * [14] S. Gracy, M. Ye, B. Anderson, and C. A. Uribe, On the endemic behavior of a competitive tri-virus sis networked model, arXiv preprint arXiv:2209.11826, (2022).
  * [15] H. W. Hethcote, The mathematics of infectious diseases, SIAM review, 42 (2000), pp. 599–653.
  * [16] I. Iacopini, G. Petri, A. Barrat, and V. Latora, Simplicial models of social contagion, Nature communications, 10 (2019), pp. 1–9.
  * [17] R. Illner and J. Ma, An sis-type marketing model on random networks, Communications in mathematical sciences, 14 (2016), pp. 1723–1740.
  * [18] A. Janson, S. Gracy, P. E. Paré, H. Sandberg, and K. H. Johansson, Analysis of a networked sis multi-virus model with a shared resource, IFAC-PapersOnLine, 53 (2020), pp. 797–802.
  * [19] 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 preprint arXiv:2011.07569, (2020).
  * [20] H. Jardón-Kojakhmetov, C. Kuehn, A. Pugliese, and M. Sensi, A geometric analysis of the SIR, SIRS and SIRWS epidemiological models, Nonlinear Analysis: Real World Applications, 58 (2021), p. 103220.
  * [21] J. O. Kephart and S. R. White, Directed-graph epidemiological models of computer viruses, in Computation: the micro and the macro view, World Scientific, 1992, pp. 71–102.
  * [22] H. K. Khalil, Nonlinear systems third edition, Patience Hall, 115 (2002).
  * [23] A. Khanafer, T. Başar, and B. Gharesifard, Stability properties of infected networks with low curing rates, in 2014 American Control Conference, IEEE, 2014, pp. 3579–3584.
  * [24] J. Köhler, L. Schwenkel, A. Koch, J. Berberich, P. Pauli, and F. Allgöwer, Robust and optimal predictive control of the covid-19 outbreak, Annual Reviews in Control, 51 (2021), pp. 525–539.
  * [25] R. M. Kovacevic, N. I. Stilianakis, and V. M. Veliov, A distributed optimal control model applied to covid-19 pandemic, SIAM Journal on Control and Optimization, 60 (2022), pp. S221–S245.
  * [26] J. S. Lavine, A. A. King, and O. N. Bjørnstad, Natural immune boosting in pertussis dynamics and the potential for long-term vaccine failure, Proceedings of the National Academy of Sciences, 108 (2011), pp. 7259–7264.
  * [27] M. Li, X. Wang, K. Gao, and S. Zhang, A survey on information diffusion in online social networks: Models and methods, Information, 8 (2017), p. 118.
  * [28] W. Li, X. Xue, L. Pan, T. Lin, and W. Wang, Competing spreading dynamics in simplicial complex, Applied Mathematics and Computation, 412 (2022), p. 126595.
  * [29] F. Liu and M. Buss, Optimal control for heterogeneous node-based information epidemics over social networks, IEEE Transactions on Control of Network Systems, 7 (2020), pp. 1115–1126, https://doi.org/10.1109/TCNS.2019.2963488.
  * [30] F. Liu, C. Shaoxuan, X. Li, and M. Buss, On the stability of the endemic equilibrium of a discrete-time networked epidemic model, IFAC-PapersOnLine, 53 (2020), pp. 2576–2581.
  * [31] J. Liu, P. E. Paré, E. Du, and Z. Sun, A networked sis disease dynamics model with a waterborne pathogen, in 2019 American Control Conference (ACC), IEEE, 2019, pp. 2735–2740.
  * [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] W.-m. Liu, H. W. Hethcote, and S. A. Levin, Dynamical behavior of epidemiological models with nonlinear incidence rates, Journal of mathematical biology, 25 (1987), pp. 359–380.
  * [34] A. L. Lloyd and R. M. May, How viruses spread among computers and people, Science, 292 (2001), pp. 1316–1317.
  * [35] E. Lukacs, Probability and mathematical statistics: an introduction, Academic Press, 2014.
  * [36] J. T. Matamalas, A. Arenas, and S. Gómez, Effective approach to epidemic containment using link equations in complex networks, Science advances, 4 (2018), p. eaau4212.
  * [37] J. T. Matamalas, S. Gómez, and A. Arenas, Abrupt phase transition of epidemic spreading in simplicial complexes, Physical Review Research, 2 (2020), p. 012049.
  * [38] W. Mei, S. Mohagheghi, S. Zampieri, and F. Bullo, On the dynamics of deterministic epidemic propagation over networks, Annual Reviews in Control, 44 (2017), pp. 116–128.
  * [39] P. E. Pare, Virus spread over networks: Modeling, analysis, and control, PhD thesis, University of Illinois at Urbana-Champaign, 2018.
  * [40] P. E. Paré, C. L. Beck, and T. Başar, Modeling, estimation, and analysis of epidemics over networks: An overview, Annual Reviews in Control, 50 (2020), pp. 345–360.
  * [41] P. E. Paré, A. Janson, S. Gracy, J. Liu, H. Sandberg, and K. H. Johansson, Multi-layer sis model with an infrastructure network, IEEE Transactions on Control of Network Systems, (2022).
  * [42] P. E. Paré, J. Liu, C. L. Beck, B. E. Kirwan, and T. Başar, Analysis, estimation, and validation of discrete-time epidemic processes, IEEE Transactions on Control Systems Technology, 28 (2018), pp. 79–93.
  * [43] P. E. Paré, J. Liu, C. L. Beck, A. Nedić, and T. Başar, Multi-competitive viruses over time-varying networks with mutations and human awareness, Automatica, 123 (2021), p. 109330.
  * [44] P. E. Paré, D. Vrabac, H. Sandberg, and K. H. Johansson, Analysis, online estimation, and validation of a competing virus model, in 2020 American Control Conference (ACC), IEEE, 2020, pp. 2556–2561.
  * [45] P. E. Paré, J. Liu, H. Sandberg, and K. H. Johansson, Multi-layer disease spread model with a water distribution network, in 2019 IEEE 58th Conference on Decision and Control (CDC), 2019, pp. 8335–8340, https://doi.org/10.1109/CDC40024.2019.9029305.
  * [46] B. Prasse and P. Van Mieghem, The viral state dynamics of the discrete-time nimfa epidemic model, IEEE Transactions on Network Science and Engineering, 7 (2019), pp. 1667–1674.
  * [47] L. Qi and Z. Luo, Tensor analysis: spectral theory and special tensors, SIAM, 2017.
  * [48] H. S. Rodrigues and M. J. Fonseca, Can information be spread as a virus? viral marketing as epidemiological model, Mathematical methods in the applied sciences, 39 (2016), pp. 4780–4786.
  * [49] S. Ruan and W. Wang, Dynamical behavior of an epidemic model with a nonlinear incidence rate, Journal of differential equations, 188 (2003), pp. 135–163.
  * [50] F. D. Sahneh and C. Scoglio, Competitive epidemic spreading over arbitrary multilayer networks, Physical Review E, 89 (2014), p. 062817.
  * [51] T. Sauer, T. Berry, D. Ebeigbe, M. M. Norton, A. J. Whalen, and S. J. Schiff, Identifiability of infection model parameters early in an epidemic, SIAM Journal on Control and Optimization, 60 (2021), pp. S27–S48.
  * [52] J. H. Shapiro, A fixed-point farrago, Springer, 2016.
  * [53] G. St-Onge, H. Sun, A. Allard, L. Hébert-Dufresne, and G. Bianconi, Universal nonlinear infection kernel from heterogeneous exposure on higher-order networks, Physical review letters, 127 (2021), p. 158301.
  * [54] L. Stella, A. P. Martínez, D. Bauso, and P. Colaneri, The role of asymptomatic infections in the covid-19 epidemic via complex networks and stability analysis, SIAM Journal on Control and Optimization, 60 (2022), pp. S119–S144.
  * [55] J. H. Tien and D. J. Earn, Multiple transmission pathways and disease dynamics in a waterborne pathogen model, Bulletin of mathematical biology, 72 (2010), pp. 1506–1533.
  * [56] D. Trpevski, W. K. Tang, and L. Kocarev, Model for rumor spreading over networks, Physical Review E, 81 (2010), p. 056102.
  * [57] P. Van Mieghem and J. Omic, In-homogeneous virus spread in networks, arXiv preprint arXiv:1306.2588, (2013).
  * [58] P. Van Mieghem, J. Omic, and R. Kooij, Virus spread in networks, IEEE/ACM Transactions On Networking, 17 (2008), pp. 1–14.
  * [59] A. Vespignani, H. Tian, C. Dye, J. O. Lloyd-Smith, R. M. Eggo, M. Shrestha, S. V. Scarpino, B. Gutierrez, M. U. Kraemer, J. Wu, et al., Modelling covid-19, Nature Reviews Physics, 2 (2020), pp. 279–281.
  * [60] L.-X. Yang and X. Yang, The impact of nonlinear infection rate on the spread of computer virus, Nonlinear Dynamics, 82 (2015), pp. 85–95.
  * [61] L.-X. Yang, X. Yang, and Y. Y. Tang, A bi-virus competing spreading model with generic infection rates, IEEE Transactions on Network Science and Engineering, 5 (2017), pp. 2–13.
  * [62] X. Yang and L.-X. Yang, Towards the epidemiological modeling of computer viruses, Discrete Dynamics in Nature and Society, 2012 (2012).
  * [63] L. Ye, P. E. Paré, and S. Sundaram, Parameter estimation in epidemic spread networks using limited measurements, SIAM Journal on Control and Optimization, 60 (2021), pp. S49–S74.
  * [64] 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.
  * [65] Y. Yi, L. Shan, P. E. Paré, and K. H. Johansson, Edge deletion algorithms for minimizing spread in sir epidemic models, SIAM Journal on Control and Optimization, 60 (2022), pp. S246–S273.
  * [66] H. Yuan, G. Liu, and G. Chen, On modeling the crowding and psychological effects in network-virus prevalence with nonlinear epidemic model, Applied Mathematics and Computation, 219 (2012), pp. 2387–2397.
  * [67] X. Zhu, Y. Nie, and A. Li, Demographic prediction of online social network based on epidemic model, in Web Technologies and Applications: APWeb 2014 Workshops, SNA, NIS, and IoTS, Changsha, China, September 5, 2014. Proceedings 16, Springer, 2014, pp. 93–103.
  * [68] L. Zino and M. Cao, Analysis, prediction, and control of epidemics: A survey from scalar to dynamic network models, IEEE Circuits and Systems Magazine, 21 (2021), pp. 4–23.