# A pair-based approximation for simplicial contagion

Federico Malizia  Department of Physics and Astronomy, University of Catania,
95125 Catania, Italy  Luca Gallo  Department of Network and Data Science,
Central European University, 1100 Vienna, Austria  Mattia Frasca  Department
of Electrical, Electronics and Computer Science Engineering, University of
Catania, 95125 Catania, Italy  Vito Latora  Department of Physics and
Astronomy, University of Catania, 95125 Catania, Italy  INFN Sezione di
Catania, Via S. Sofia, 64, 95125 Catania, Italy  School of Mathematical
Sciences, Queen Mary University of London, London E1 4NS, UK  Complexity
Science Hub Vienna, A-1080 Vienna, Austria  Giovanni Russo  Department of
Mathematics and Computer Science, University of Catania, 95125 Catania, Italy

###### Abstract

Pairwise interactions alone are often insufficient to characterize contagion
processes, as more complex mechanisms involving groups of three or more
individuals may be at play. Such higher-order interactions can be effectively
modeled using frameworks beyond complex networks, such as simplicial
complexes. The presence of these higher-order interactions has been shown to
play a critical role in shaping the onset and evolution of contagion
processes. However, studying these dynamics can be challenging due to the high
dimensionality of the state space of these structures. To solve this problem,
numerous mean-field models have been developed. Nevertheless, these models
often overlook the correlations between different subsets of nodes, which can
significantly influence the system dynamics. Therefore, more detailed
approximations that account for these correlations are necessary. In this
paper, we present a novel pair-based approximation for studying SIS dynamics
on simplicial complexes. The pair-based approximation takes into consideration
the dynamical correlations that emerge within groups of nodes in a simplicial
complex. Compared to individual-based mean-field approaches, this
approximation yields more accurate predictions of the behavior observed in
stochastic simulations of contagion processes on simplicial complexes.
Specifically, the proposed pair-based approximation provides higher accuracy
in predicting the extent of the region of bistability, the type of the
transition from a disease-free to an endemic state, and the average time
evolution of the fraction of infected individuals. Overall, our findings
highlight the significance of accounting for correlations within groups of
nodes when investigating dynamical processes on simplicial complexes, and
suggest that the pair-based approach can provide valuable insights into the
behavior of such systems.

##  I Introduction

Complex networks have been used to describe a large variety of dynamical
processes involging interacting units, such as epidemic spreading [1, 2],
social contagion [3], random walks [4], and synchronization [5], among many
others. However, networks have inherent limitations as they can only capture
interactions between pairs of units [6, 7]. Consequently, they cannot be
employed to study complex systems characterized by interactions occurring in
groups of three or more units [8]. In such cases, more sophisticated
mathematical structures, such as hypergraphs [9] and simplicial complexes
[10], are required. The inclusion of higher-order interactions has proved that
these interactions give rise to novel collective phenomena in various
dynamical processes [11], including diffusion [12, 13], synchronization [14,
15, 16], and evolutionary games [17].

In particular, the framework of higher-order networks has proven to be
essential in modeling social contagion processes, such as opinion formation,
rumor spreading, or the adoption of novelties, where the exposure to multiple
sources is needed to trigger the transmission [18]. Recent studies of social
contagion in hypergraphs [19] and simplicial complexes [20] have revealed that
higher-order interactions may radically change the characteristics of the
spreading process. In the presence of higher-order interactions, the
transition to an endemic state in SIS models can become discontinuous and a
bistable regime where a disease-free and an endemic state co-exist can appear
[20, 19]. This behavior is closely connected to the microscopic organization
of higher-order interactions as diverse structures may exhibit different types
of transitions [21].

Most models for the analysis of social contagion processes with higher-order
interactions rely on individual-based mean-field approaches, assuming that the
individuals are homogeneously mixed and interact with each other at random.
However, this assumption may lead to an oversimplified approximation, as it
overlooks the dynamic correlations that arise within the higher-order network.
For example, infected individuals are more likely to come into contact with
other infected individuals [2].

In the context of complex networks, to obtain a more precise description of
the spreading process, one can consider more sophisticated approximations. For
instance, pair-based models can be employed [22, 23, 24, 25]. In these models,
the system dynamics is not characterized at the level of nodes, that is, it
does not follow the temporal evolution of the expected number of individuals
in a given state, but it is described at the level of pairs of nodes, namely
it analyzes how the expected number of edges in a given state evolves in time.

However, deriving mean-field models for spreading processes on structures with
higher-order interactions is more challenging, resulting in a smaller number
of available approaches. In this work, we consider the same contagion process
studied in [20], namely an SIS process in the presence of three-body
simplicial interactions, and develop a mean-field model based on a description
at the level of the pairs, which we then analyse in comparison with the
individual-based mean-field model introduced in [20].

In the past years, a few works have investigated the development of mean-field
models for contagion dynamics in the presence of higher-order interactions
[26, 27]. In particular, in [26] a model relying on a discrete-time Markov-
chain approach has been proposed. In this model, the system dynamics is
expressed in terms of joint probabilities of the microscopic states of the
links and the nodes, and thus is described by a set of $N+L$ master equations,
where $N$ is the number of nodes and $L$ the number of links. In our work,
instead, we consider a continuous-time framework yielding a system of
governing equations whose order does not depend on the number of nodes and
links. The model developed in [27] also adopts a continuous-time framework,
but relies on the assumption that adjacent hyperedges share at most one node.
In our work, instead, we tackle the more general case that includes all
possible motifs of three and four nodes that can exist in a simplicial
complex.

The rest of the paper is organized as follows. In Sec. II we provide some
preliminary results, recalling the individual-based and the pair-based
approximation of the SIS model on networks, and then describing the
individual-based approximation for the SIS model on simplicial complexes. In
Sec. III we introduce the pair-based approximation of the simplicial SIS model
and its derivation. In Sec. IV, we compare the predictions of the individual-
based and the pair-based approximation for the simplicial SIS model. Finally,
in Sec. V we summarize and discuss the main results obtained.

##  II Mean-field SIS models

In this section, we first discuss two different mean-field approximations for
the SIS model on networks, namely a contagion process where the disease can be
transmitted exclusively occurs through pairwise interactions. We describe both
the individual-based and pair-based approximations, highlighting the
fundamental mechanisms involved in deriving mean-field models at the
individual- and pair-based levels. Subsequently, we revisit the individual-
based approximation for the simplicial SIS model studied in [20]. This will
lay the groundwork for introducing the pair-based approximation for the
simplicial SIS model in Sec. III.

###  II.1 Individual-based mean-field SIS model

We begin by characterizing the processes governing the transition of an
individual from one state to another in the standard SIS model. In this model,
an individual can be either susceptible (S) or infected/infectious (I), and
can transit from one compartment to another as infection or recovery takes
place. In the former process, a susceptible individual (S) becomes infected
(S$\rightarrow$I) after interacting with an infectious individual (I), which
acts a mediator of the transition. Hence, the infection mechanism is a two-
body nonlinear process. In the second case, an infected individual (I)
recovers after a given period of time, becoming once again susceptible
(I$\rightarrow$S). At variance with the infection mechanism, the recovery is
described by a one-body linear process. Formally, we can express these
transitions in terms of two kinetic equations

$\begin{array}[]{rcl}S+I&\overset{\beta}{\rightarrow}&I+I\\\
I&\overset{\mu}{\rightarrow}&S\end{array},$ (1)

where $\beta$ and $\mu$ represent the transition rates for the infection and
the recovery process, respectively.

Figure 1: Pictorial representation of infection processes among susceptible
(in blue) and infected (in red) individuals in SIS models. (a) Infection of a
node connected to an infected node through a link. This is the only infection
process occurring in the individual-based SIS model. (b) Infection in a pair
of nodes. In a pair-based description of contagion in a SIS model, the disease
can be transmitted to one of the nodes of a pair by the other node of the pair
or by other infected nodes, not belonging to the pair. (c) Simplicial
infection in a triangle by group interaction with two infected nodes. This
infection process characterizes contagion in simplicial SIS models. (d)
Infection in a pair of nodes in the presence of group interactions. The
infection processes depicted in panel (d) arise in the pair-based description
of the simplicial SIS model.

The dynamics of the SIS model on a network can be investigated using different
approaches [2, 1]. A commonly employed strategy is to perform stochastic
simulations of the kinetic equations (1). However, this approach has some
limitations. Firstly, it can be computationally expensive for large networks.
Additionally, it requires a significant amount of data, as both the contact
network and the states of the nodes need to be known empirically. Moreover,
this approach does not facilitate understanding of the relationship between
the model parameters and the emerging collective behavior [25]. For this
reason, a common approach is to make use of a deterministic representation of
the process that focuses on population-level quantities, such as the density
or the expected number of individuals in a given state. This deterministic
approach assumes homogeneous mixing, where the interactions between
individuals are treated as uniform and independent of the network structure,
so that each individual can be considered statistically equivalent to any
other. This hypothesis allows to describe the dynamics of the system in terms
of a so-called mean-field model.

For the SIS model under homogeneous mixing, the dynamics is described by the
following equations

$\begin{array}[]{ccc}\dot{\langle S\rangle}&=&\mu\langle I\rangle-\beta
k\langle SI\rangle\\\ \dot{\langle I\rangle}&=&-\mu\langle I\rangle+\beta
k\langle SI\rangle\end{array},$ (2)

where $\langle S\rangle$ and $\langle I\rangle$ represent the densities of
susceptible and infected individuals, respectively, $\langle SI\rangle$ is the
density of edges in the state $(S,I)$, i.e., a susceptible individual in
contact with an infectious one, while $k$ is the average node degree. The
densities of individuals in state S and I are given by $\langle
S\rangle=\left[S\right]/N$ and $\langle I\rangle=\left[I\right]/N$, where
$\left[S\right]$ and $\left[I\right]$ are the expected number of susceptible
and infected individuals, respectively, and $N=[S]+[I]$ is the total number of
individuals in the population. Consequently, $\langle S\rangle$ and $\langle
I\rangle$ are related through the conservation relation $\langle
S\rangle+\langle I\rangle=1$, so they are not independent. The density
$\langle SI\rangle$ is defined as $\langle SI\rangle=\left[SI\right]/kN$,
where $\left[SI\right]$ represents the expected number of edges in state
$(S,I)$. Note that $kN=2L$, where $L$ is the total number of edges within the
network. The two terms in the equation for $\langle I\rangle$ represent the
decrease in the density of infectious individuals due to recoveries, and the
increase in the same density due to infections, respectively. The first
depends on the transition rate $\mu$ and on the density of infected
individuals, while the second is given by the product of the transition rate
$\beta$ by the average degree $k$ and the edge density $\langle SI\rangle$.
The transition from susceptible to infected through the interaction with an
infected individual is depicted in panel (a) of Fig. 1.

System (2) is exact but not closed, as the equations governing the dynamics of
both $\langle S\rangle$ and $\langle I\rangle$ depend on the quantity $\langle
SI\rangle$. To obtain a closed system of equations, one may apply the law of
mass action, which consists of assuming statistical independence at the level
of individuals. This means to assume that infected individuals are randomly
distributed in the network, so that the probability that a neighbor is
infected is given by $\langle I\rangle$ and does not depend on the state of
the node itself. Under this assumption, by approximating the fraction of edges
in state $(S,I)$ as

$\langle SI\rangle\approx\langle S\rangle\langle I\rangle,$ (3)

and substituting this expression in Eqs. (2), we obtain the following closed
form for the mean-field SIS model:

$\begin{array}[]{ccc}\dot{\langle S\rangle}&=&\mu\langle I\rangle-\beta
k\langle S\rangle\langle I\rangle\\\ \dot{\langle I\rangle}&=&-\mu\langle
I\rangle+\beta k\langle S\rangle\langle I\rangle\\\ \end{array}.$ (4)

###  II.2 Pair-based mean-field SIS model

Assuming that infected individuals are randomly distributed in the network may
lead to an oversimplified approximation, as it does not account for the
dynamic correlations that exist within the contact network (for example,
infected nodes are more likely to come into contact with other infected nodes)
[1]. To provide a more accurate description of the SIS dynamics on a network,
one can consider a pair-based model, which incorporates these dynamic
correlations. In practice, one has to characterize the system dynamics at the
level of pairs of nodes, describing how the expected number of edges in a
given state evolves in time. For the SIS process, the equations describing the
SIS process at the pair-level are:

$\begin{array}[]{ll}\dot{\langle S\rangle}=&\mu\langle I\rangle-\beta k\langle
SI\rangle\\\ \dot{\langle I\rangle}=&-\mu\langle I\rangle+\beta k\langle
SI\rangle\\\ \dot{\langle SS\rangle}=&2\mu\langle SI\rangle-2\beta(k-1)\langle
SSI\rangle\\\ \dot{\langle SI\rangle}=&\mu\langle II\rangle-\mu\langle
SI\rangle+\beta(k-1)\langle SSI\rangle\\\ &-\beta(k-1)\langle
ISI\rangle-\beta\langle SI\rangle\\\ \dot{\langle II\rangle}=&-2\mu\langle
II\rangle+2\beta(k-1)\langle ISI\rangle+2\beta\langle SI\rangle\end{array},$
(5)

where $\langle ISI\rangle$ represents the density of (open and closed)
triplets in state $(I,S,I)$, and is defined as $\langle
ISI\rangle=\left[ISI\right]/(k(k-1)N)$, with $\left[ISI\right]$ being the
expected number of triplets in that state. Note that the term $k(k-1)$
corresponds to the average number of triplets connected to each node. Indeed,
each node $i$ with $k$ neighbors is at the center of $k(k-1)/2$ unordered
triples. $\langle S\rangle$ and $\langle I\rangle$ can be obtained from
$\langle SS\rangle$, $\langle SI\rangle$ and $\langle II\rangle$ through
marginalization, namely $\langle S\rangle=\langle SS\rangle+\langle SI\rangle$
and $\langle I\rangle=\langle SI\rangle+\langle II\rangle$. Additionally,
$\langle SS\rangle$, $\langle SI\rangle$ and $\langle II\rangle$ are not
independent, as they are linked through the conservation relation $\langle
SS\rangle+2\langle SI\rangle+\langle II\rangle=1$. Note that, with respect to
the individual-based approximation, the pair-based SIS model contains a higher
number of terms, reflecting the higher number of possible states and
transitions that a pair has compared to a single node.

To illustrate the terms appearing in Eqs. (5), let us focus on the equation
governing the dynamics of $\langle SI\rangle$. The first two terms are related
to recoveries. Indeed, the first term represents the increase in the density
$\langle SI\rangle$ due to the recovery of one of the nodes of the pairs in
state $(I,I)$. The second term models the decrease in $\langle SI\rangle$ due
to the recovery of the infected nodes of the pairs in state $(S,I)$, i.e., the
transition from $(S,I)$ to $(S,S)$. The remaining three terms encode the
possible transitions of an edge due to infections, as graphically represented
in panel (b) of Fig. 1. In particular, the first encodes the increase in
$\langle SI\rangle$ due to the infection of one of the nodes in the pairs in
$(S,S)$, while the second represents the decrease in $\langle SI\rangle$ due
to the infection of the susceptible nodes in the pairs in state $(S,I)$. Both
transitions are due to the interaction with a third (infectious) node, among
the $(k-1)$ remaining neighbors of the susceptible nodes. Note that the three
nodes involved in the interaction can be arranged in two possible motifs,
namely an open triangle, i.e., a wedge, or a closed one. This point will turn
out to be crucial when considering the closure of the system. The last term
corresponds to the decrease in $\langle SI\rangle$ due to the infection of the
susceptible nodes in the pairs in state $(S,I)$. In this case, however, the
transition is due to the infected nodes of the pair, and not to a third node.

System (5) is exact but not closed. Different approximations of the triplet
densities have been proposed to close the system at the pair level. The most
common and studied one [1] consists of writing the probability that a motif is
in a given state as the product of the probabilities that its edges are in a
certain state, normalized by the probabilities that the nodes in common to the
pairs of edges are in a particular state. To illustrate it, let us first
consider a triple of nodes $i$, $j$, and $z$, in states A, B, and C
respectively, forming a wedge, so that node $j$ is connected to nodes $i$ and
$z$, while $i$ and $z$ are not linked. We can write the density of wedges in
state $(A,B,C)$ as

$\begin{array}[]{l}\langle ABC^{\wedge}\rangle\approx\dfrac{\langle
AB\rangle\langle BC\rangle}{\langle B\rangle},\\\ \end{array}$ (6)

where $\langle AB\rangle=P(AB)$ and $\langle BC\rangle=P(BC)$ are the
probabilities that links $(i,j)$ and $(j,k)$ are in states $(A,B)$ and
$(B,C)$, respectively, while $\langle B\rangle$ is the probability that node
$j$ is in state B. This formula can be obtained as follows.

$\langle ABC^{\wedge}\rangle=P(AB\cap BC)=P(BC)P(AB|BC)$

where $P(AN|BC)$ denotes the conditional probability that $i$ and $j$ are in
states $A$ and $B$, given that $j$ and $z$ are in states $B$ and $C$. Then we
make the approximation

$P(AB|BC)\approx P(AB|B)=P(AB)/P(B)$

Combining these two relations we obtain Eq. (6). When there are few short
cycles, i.e., few closed triangles, in the network, Eq. (6) represents an
accurate approximation of the probability that a triplet of nodes is in state
$(A,B,C)$. However, this closure fails to describe the epidemic dynamics on
networks where the number of closed triplets is not negligible [1, 22], as it
overlooks the correlation between nodes $i$ and $k$. Therefore, in the case of
a triple of nodes $i$, $j$, and $z$ forming a closed triangle, one usually
approximates the density of triangles is in state $(A,B,C)$ using the so-
called Kirkwood superposition [28], namely 111The formula can be understood as
follows: $P\left(AB\cap BC\cap CA\right)=P\left(AB\cap
BC|CA\right)P\left(CA\right)$. Now $P\left(AB\cap
BC|CA\right)=P\left(AB|BC\cap CA\right)\cdot P\left(BC|CA\right)$. We can
write $P\left(AB|BC\cap
CA\right)=\dfrac{P\left(AB\right)}{P\left(A\right)P\left(B\right)}$, and use
the approximation $P\left(AB|BC\cap
CA\right)=\dfrac{P\left(AB\right)}{P\left(A\right)P\left(B\right)}$.
Substituting in the previous expression one has $P(AB\cap
BC|CA)=\dfrac{P\left(AB\right)}{P\left(A\right)P\left(B\right)}\cdot\dfrac{P\left(BC\right)}{P\left(C\right)}$
from which it follows $P\left(AB\cap BC\cap
CA\right)={P(AB)P(BC)P(CA)}/{(P(A)P(B)P(C))}$

$\begin{array}[]{l}\langle ABC^{\Delta}\rangle\approx\dfrac{\langle
AB\rangle\langle BC\rangle\langle AC\rangle}{\langle A\rangle\langle
B\rangle\langle C\rangle}\end{array}.$ (7)

Note that, when we assume the dynamics of nodes $i$ and $z$ to be
uncorrelated, i.e., $\langle AC\rangle=\langle A\rangle\langle C\rangle$, we
find that the approximation for $\langle ABC^{\Delta}\rangle$ recovers the one
for $\langle ABC^{\wedge}\rangle$.

The next step is to calculate the probability that a generic triplet of nodes
is in state $(A,B,C)$, namely $\langle ABC\rangle$, as a function of $\langle
ABC^{\wedge}\rangle$ and $\langle ABC^{\Delta}\rangle$. To this aim, let us
consider the global clustering coefficient $\phi\in[0,1]$, representing the
fraction of closed triangles over all triplets in the structure, and calculate
it in a statistical meaning as

$\begin{array}[]{l}\displaystyle\phi=\dfrac{3\mathbb{E}(\mathcal{T})}{Nk(k-1)},\end{array}$
(8)

where $N$ is the number of nodes, $k$ is the average number of links (i.e.,
1-simplices) per node, and $\mathbb{E}(\mathcal{T})$ represents the expected
number of triangles, which depends on the structure considered. At this point,
we can calculate $\langle ABC\rangle$ as [22, 1]:

$\langle ABC\rangle=(1-\phi)\langle ABC^{\wedge}\rangle+\phi\langle
ABC^{\Delta}\rangle.$ (9)

Note that the densities $\langle ABC^{\wedge}\rangle$ and $\langle
ABC^{\Delta}\rangle$ are defined as $\langle
ABC^{\wedge}\rangle=\left[ABC^{\wedge}\right]/k(k-1)(1-\phi)N$ and $\langle
ABC^{\Delta}\rangle=\left[ABC^{\Delta}\right]/k(k-1)\phi N$, where
$\left[ABC^{\wedge}\right]$ $\left[ABC^{\Delta}\right]$ are the expected
number of wedges and of triangles in state $(A,B,C)$ respectively.

Hereafter, we will consider Eq. (9) to close the compartmental models at the
pair level. In particular, we will approximate the densities $\langle
SSI\rangle$ and $\langle ISI\rangle$ as:

$\begin{array}[]{l}\langle SSI\rangle\approx\left(1-\phi\right)\dfrac{\langle
SS\rangle\langle SI\rangle}{\langle S\rangle}+\phi\dfrac{\langle
SS\rangle\langle SI\rangle^{2}}{\langle S\rangle^{2}\langle I\rangle}\\\
\langle ISI\rangle\approx\left(1-\phi\right)\dfrac{\langle
SI\rangle^{2}}{\langle S\rangle}+\phi\dfrac{\langle SI\rangle^{2}\langle
II\rangle}{\langle S\rangle\langle I\rangle^{2}}\end{array}.$ (10)

Such closure relations are a direct consequence of Eqs. (6)-(9). Eqs. (5)
along with the expressions (10) constitute a closed approximation for the
pair-based mean-field SIS model.

###  II.3 Individual-based mean-field simplicial SIS model

We now discuss the individual-based mean-field simplicial SIS model, presented
in [20]. As done for the SIS model on networks, i.e., in the absence of
higher-order interactions, we first describe the processes ruling the
transition of an individual from one state to another, this time in the
presence of high-order interactions. In this case, in addition to the
infection processes occurring in the SIS model on networks, there is a further
way in which an individual transits from one compartment to another. In fact,
a susceptible individual (S) can become infected (S$\rightarrow$I) through a
three-body interaction, in which the other two individuals are infectious. In
this case, both infected individuals act as mediators of the transition.
Similarly to the network case, we can represent the processes of the
simplicial SIS model in terms of kinetic equations

$\begin{array}[]{rcl}S+I&\overset{\beta}{\rightarrow}&I+I\\\
S+I+I&\overset{\beta_{\Delta}}{\rightarrow}&I+I+I\\\
I&\overset{\mu}{\rightarrow}&S\end{array},$ (11)

where $\beta$ and $\beta_{\Delta}$ are the transition rates for the two-body
and the three-body infection process, respectively, while $\mu$ is the
recovery rate.

We now derive a model that describes the system dynamics in terms of
population-level quantities. Under the homogeneous mixing hypothesis, the
exact equations for the simplicial SIS model are given by

$\begin{array}[]{ll}\dot{\langle S\rangle}=&\mu\langle I\rangle-\beta k\langle
SI\rangle-\beta_{\Delta}k(k-1)\phi\delta\langle ISI^{\Delta}\rangle\\\
\dot{\langle I\rangle}=&-\mu\langle I\rangle+\beta k\langle
SI\rangle+\beta_{\Delta}k(k-1)\phi\delta\langle
ISI^{\Delta}\rangle\end{array},$ (12)

where $\delta\in[0,1]$ is the fraction of triangles that are effectively
$2$-simplices. Compared to Eqs. (2), these equations present an additional
term, i.e., $\beta_{\Delta}k(k-1)\phi\delta\langle ISI^{\Delta}\rangle$, which
takes into account the contribution of the three-body interactions, namely the
infection of a susceptible node due to the simultaneous interaction with two
infected individuals, as shown in panel (c) of Fig. 1.

System (12) may be closed by applying the law of mass action, starting from
the assumption that infected individuals are randomly distributed in the
simplicial complex [20]. Accordingly, the density $\langle
ISI^{\Delta}\rangle$ is approximated as

$\langle ISI^{\Delta}\rangle\approx\langle I\rangle\langle S\rangle\langle
I\rangle,$ (13)

which allows us to write the individual-based mean-field model for the
simplicial SIS model as follows:

$\begin{array}[]{ll}\dot{\langle S\rangle}=&\mu\langle I\rangle-\beta k\langle
S\rangle\langle I\rangle-\beta_{\Delta}k_{\Delta}\langle S\rangle\langle
I\rangle^{2}\\\ \dot{\langle I\rangle}=&-\mu\langle I\rangle+\beta k\langle
S\rangle\langle I\rangle+\beta_{\Delta}k_{\Delta}\langle S\rangle\langle
I\rangle^{2}\end{array},$ (14)

where $k_{\Delta}=k(k-1)\phi\delta$ represents the average number of
$2$-simplices connected to each node.

The system can be analytically investigated in terms of two parameters,
$\lambda=k\beta/\mu$ and $\lambda_{\Delta}=k_{\Delta}\beta_{\Delta}/\mu$,
representing respectively the rescaled infectivity on 1-simplices and on
2-simplices. In particular, this permits to obtain the steady-state solutions
of Eqs. (14) and to study their stability as a function of $\lambda$ and
$\lambda_{\Delta}$. The analysis shows that, when the higher-order
interactions are weak, i.e., when $\lambda_{\Delta}<1$, the system behaves
similarly to the SIS model on networks, described by Eq. (2): for $\lambda<1$
the disease-free equilibrium $\langle I\rangle^{*}=0$ is the only solution,
while for $\lambda>1$ a stable endemic state $\langle I\rangle^{*}\neq 0$
exists. In addition, the phase transition at $\lambda=1$ is continuous.
However, when the higher-order interactions are stronger, i.e.,
$\lambda_{\Delta}>1$, the system shows a different behavior. For $\lambda$
smaller than $\lambda_{c}=2\sqrt{\lambda_{\Delta}}-\lambda_{\Delta}$, only the
disease-free equilibrium exists. For $\lambda_{c}<\lambda<1$, the disease-free
equilibrium and an endemic state coexist in a bistable regime, with the
initial fraction of infected individuals determining whether the system
reaches one equilibrium point or the other: if the initial fraction of
infected individuals is larger than a critical mass, then the system evolves
towards the endemic state; otherwise the disease spreading cannot be
sustained, and the system evolves towards the disease-free equilibrium.
Finally, for $\lambda>1$, the endemic state is the only stable equilibrium.
Differently from the SIS model on networks, when $\lambda_{\Delta}>1$ the
phase transition occurring at $\lambda=\lambda_{c}$ and $\lambda=1$ is
discontinuous.

##  III A pair-based approximation for the simplicial SIS model

In this section, we derive the pair-based mean-field simplicial SIS model. To
this aim, we have to account for all ways a pair of nodes in a simplicial
complex can transit from one state to another. Besides the ones already
discussed in the derivation of the previous models, the transitions
illustrated in panel (d) of Fig. 1 should also be considered. Remarkably, two
of these transitions depend on the state of 4-node motifs, and occur when the
$2$-simplex shares with the edge of interest a single node. Instead, the
remaining one depends of the state of a 3-node motif, and corresponds to the
case where the focal edge belongs to the $2$-simplex. In the first two cases,
an infection can occur when the shared node is susceptible and the other nodes
of the $2$-simplex are infected, while the other node in the pair can be
either susceptible or infected. Hence, both the densities of edges in states
$(S,S)$ and $(S,I)$ are affected by this infection process, with a rate that
is function of the rate $\beta_{\Delta}$ and of the density of 4-node motifs
in the aforementioned states. As shown in Fig. 2, there are three possible
motifs that contribute to this infection process. They differ for the number
of links (zero, one or two) that connect the node external to the $2$-simplex
to the infected nodes of the $2$-simplex.

Figure 2: Graphical representation of the three possible microscopical
configurations of four-node motif states $(I,I,S,S)$ (on top) and $(I,I,S,I)$
(bottom). Square brackets refer to the expected number of the singular
configurations.

Hereby, we denote as $\langle IIS^{\Delta}S\rangle_{0}$, $\langle
IIS^{\Delta}S\rangle_{1}$, and $\langle IIS^{\Delta}S\rangle_{2}$ the
densities of the these motifs in state $(I,I,S,S)$, while we denote as
$\langle IIS^{\Delta}I\rangle_{0}$, $\langle IIS^{\Delta}I\rangle_{1}$, and
$\langle IIS^{\Delta}I\rangle_{2}$ the density of motifs in state $(I,I,S,I)$.
More in general, using the clustering coefficient $\phi$, we can write the
density of 4-node motifs in states $(I,I,S,S)$ and $(I,I,S,I)$ as

$\begin{array}[]{lll}\langle IIS^{\Delta}S\rangle&=&(1-\phi)^{2}\langle
IIS^{\Delta}S\rangle_{0}\\\ &&+2(1-\phi)\phi\langle
IIS^{\Delta}S\rangle_{1}\\\ &&+\phi^{2}\langle
IIS^{\Delta}S\rangle_{2}\\\\[10.0pt] \langle
IIS^{\Delta}I\rangle&=&(1-\phi)^{2}\langle IIS^{\Delta}I\rangle_{0}\\\
&&+(1-\phi)\phi\langle IIS^{\Delta}I\rangle_{1}\\\ &&+\phi^{2}\langle
IIS^{\Delta}I\rangle_{2}\\\ \end{array}$ (15)

In order to quantify the contribution to the model equations of the two
transitions involving the state of 4-node motifs, we also need to calculate
the average number of motifs in which the pair of nodes is not part of the
$2$-simplex. As each node is connected to $k$ links on average, the average
number of motifs composed by a $2$-simplex connected to a link is
$k(k-1)\phi\delta(k-2)$. Indeed, the average number of $2$-simplices connected
to a node is $k(k-1)\phi\delta$, while the fourth node of the motif has to be
chosen among the $(k-2)$ remaining neighbors of the node.

Let us now go back to the infection process that is shown in panel (d) of Fig.
1 and involves the 3-node motifs. In this case, the node pair is part of the
$2$-simplex. Here, the only possible transition for the edge is the one from
state $(S,I)$ to state $(I,I)$. In fact, both the two other nodes of the
$2$-simplex have to be in the infectious state to yield a simplicial
contagion, as expressed by the second kinetic equation (11). Similarly to the
previous case, the overall contribution present in the model will depend on
the infection rate $\beta_{\Delta}$, on the density of $2$-simplices in state
$(I,I,S)$, which we denote as $\langle IIS^{\Delta}\rangle$, and on the
average number of triangles an edge is part of. In particular, the latter is
given by $2(k-1)\phi\delta=2k_{\Delta}/k$, where the factor 2 comes from the
fact that each $2$-simplex has two edges pointing to a node.

Taking into account the previous considerations, we can finally write the
equations governing the simplicial SIS model at the pair-level:

$\begin{array}[]{ll}\dot{\langle S\rangle}=&\mu\langle I\rangle-\beta k\langle
SI\rangle-\beta_{\Delta}k_{\Delta}\langle ISI^{\Delta}\rangle\\\\[5.0pt]
\dot{\langle I\rangle}=&-\mu\langle I\rangle+\beta k\langle
SI\rangle+\beta_{\Delta}k_{\Delta}\langle ISI^{\Delta}\rangle\\\\[5.0pt]
\dot{\langle SS\rangle}=&2\mu\langle SI\rangle-2\beta(k-1)\langle
SSI\rangle\\\\[5.0pt] &-2\beta_{\Delta}\frac{k_{\Delta}}{k}(k-2)\langle
IIS^{\Delta}S\rangle\\\\[5.0pt] \dot{\langle SI\rangle}=&\mu\langle
II\rangle-\mu\langle SI\rangle+\beta(k-1)\langle SSI\rangle\\\\[5.0pt]
&-\beta(k-1)\langle ISI\rangle-\beta\langle SI\rangle\\\\[5.0pt]
&+\beta_{\Delta}\frac{k_{\Delta}}{k}(k-2)\langle
IIS^{\Delta}S\rangle-2\beta_{\Delta}\frac{k_{\Delta}}{k}\langle
ISI^{\Delta}\rangle\\\\[5.0pt]
&-\beta_{\Delta}\frac{k_{\Delta}}{k}(k-2)\langle
IIS^{\Delta}I\rangle\\\\[5.0pt] \dot{\langle II\rangle}=&-2\mu\langle
II\rangle+2\beta(k-1)\langle ISI\rangle+2\beta\langle SI\rangle\\\\[5.0pt]
&+4\beta_{\Delta}\frac{k_{\Delta}}{k}\langle
ISI^{\Delta}\rangle+2\beta_{\Delta}\frac{k_{\Delta}}{k}(k-2)\langle
IIS^{\Delta}I\rangle\end{array}.$ (16)

Here, we remark a property of the model. Analogously to what happens for the
SIS model on networks, where the recovery is a 1-body process, and the
infection a 2-body one, here, in both the individual-based and the pair-based
exact models of simplicial SIS contagion dynamics, we find terms functions of
both $m$-body and $(m+1)$-body quantities. For instance, the dynamics of
$\langle I\rangle$ (1-body variable) depends on $\langle I\rangle$ itself and
on $\langle SI\rangle$ (2-body). Compared the model on networks, however, in
the simplicial SIS model a 3-body infection process is also present.
Consequently, for both the individual-based and the pair-based exact models,
the equations governing the evolution of the $m$-body quantities are
determined by $m$-body, $(m+1)$-body, and $(m+2)$-body quantities. Hence, in
the simplicial SIS model, the dynamics of $\langle I\rangle$ (1-body variable)
depends not only on $\langle I\rangle$ itself and on $\langle SI\rangle$
(2-body), but also on $\langle ISI^{\Delta}\rangle$ (3-body).

Similarly to what we have done for Eqs. (5), we now close Eqs. (16) at the
level of pairs. To this aim, we use the approximations for the triplet
densities $\langle ABC\rangle$ at the pair level given by Eqs. (6), (7), and
(9). For what concerns the densities of the $4$-node motifs, we rely on the
closures introduced in [30] (see Eqs. 15), after re-adapting them to the
simplicial case. Accordingly, we consider

$\begin{array}[]{l}\langle IIS^{\Delta}S\rangle_{0}=\dfrac{\langle
SI\rangle^{2}\langle II\rangle\langle SS\rangle}{\langle S\rangle^{2}\langle
I\rangle^{2}}\\\ \langle IIS^{\Delta}S\rangle_{1}=\dfrac{\langle
SI\rangle^{3}\langle SS\rangle\langle II\rangle}{\langle S\rangle^{3}\langle
I\rangle^{3}}\\\ \langle IIS^{\Delta}S\rangle_{2}=\dfrac{\langle
SS\rangle\langle SI\rangle^{4}\langle II\rangle}{\langle S\rangle^{4}\langle
I\rangle^{4}}\\\ \langle IIS^{\Delta}I\rangle_{0}=\dfrac{\langle
SI\rangle^{3}\langle II\rangle}{\langle S\rangle^{2}\langle I\rangle^{2}}\\\
\langle IIS^{\Delta}I\rangle_{1}=\dfrac{\langle SI\rangle^{3}\langle
II\rangle^{2}}{\langle S\rangle^{2}\langle I\rangle^{4}}\\\ \langle
IIS^{\Delta}I\rangle_{2}=\dfrac{\langle SI\rangle^{3}\langle
II\rangle^{3}}{\langle S\rangle^{2}\langle I\rangle^{6}}\\\ \end{array},$ (17)

where we defined the 4-node motif states densities as a function of the pair
state variables. Summing up, Eqs. (16) along with the expressions (6), (7),
(9), and (17) constitute a closed approximation for the pair-based mean-field
simplicial SIS model.

##  IV Results

In this section we analyse the pair-based approximation of the simplicial SIS
model, and compare it with the corresponding individual-based approximation.
In particular, we will focus on the fraction of infected individuals, denoted
as $\rho(t)=\langle I(t)\rangle$, and on its steady-state value $\rho^{*}$. We
will study $\rho^{*}$ as a function of the model parameters, specifically the
rescaled infection rates $\lambda=k\beta/\mu$ and
$\lambda_{\Delta}=k_{\Delta}\beta_{\Delta}/\mu$, and also compare the
prediction of the models with the results of numerical, stochastic simulations
on a random simplicial complex (RSC) [20]. We will show that pair-based
simplicial SIS model better predicts both the epidemics threshold and the
nature of the phase transitions, i.e., continuous or discontinuous, of the
spreading process.

While for the individual-based model the stationary density $\rho^{*}$ can be
evaluated analytically [20], the pair-based model is not analytically
tractable, and we were not able to find a closed-form expression for the
steady-state equilibrium. For this reason, the stationary density $\rho^{*}$
is obtained by numerically integrating Eqs. (16) for a sufficiently long time
window, such that the system reaches the steady-state.

Since in the individual-based approximation of the simplicial SIS model,
$\rho^{*}$ as a function of $\lambda$ can exhibit a bistable regime [20], to
check whether the same behavior also appears in the pair-based approximation,
for each setting of the parameters we carry out simulations starting from two
different initial conditions of $\rho(0)$, one with a very low fraction of
initially infected individuals, $\rho(0)=10^{-6}$, and another with a higher
fraction, $\rho(0)=0.8$. Once fixed the initial fraction of infected
individuals $\langle I\rangle_{0}$, from which it follows $\langle
S\rangle_{0}=1-\langle I\rangle_{0}$, the initial conditions of the state
variables of the pair-based approximation are set as:

$\begin{array}[]{l}\langle SS\rangle_{0}=\langle S\rangle_{0}^{2}\\\ \langle
SI\rangle_{0}=\langle S\rangle_{0}\langle I\rangle_{0}\\\ \langle
II\rangle_{0}=\langle I\rangle_{0}^{2}\end{array}$ (18)

Finally, to integrate the pair-based approximation, we need to evaluate
$\phi$, namely the fraction of closed triangles over all triplets in the
simplicial complex, which is an input parameter of Eqs. (15). For a RSC with
$N$ nodes, an average number of $1$-simplices equal to $k$, and an average
number of $2$-simplices equal to $k_{\Delta}$, the expected number of closed
triplets is given by:

$\begin{array}[]{l}\displaystyle\mathbb{E}(\mathcal{T})\approx\binom{N}{3}\dfrac{k}{(N-1)}+\binom{N}{3}\dfrac{2k_{\Delta}}{(N-1)(N-2)}.\end{array}$
(19)

After some algebraic manipulations and using the definition in Eq. (8), the
global clustering coefficient $\phi$ for the RSC can be approximated as

$\begin{array}[]{l}\phi\approx\dfrac{k^{2}(N-2)}{(N-1)^{2}(k-1)}+\dfrac{2k_{\Delta}}{k(k-1)}.\end{array}$
(20)

For what concerns the stochastic simulations on the RSC, they are carried out
on a RSC with $N=2000$ nodes, $k\approx 20$ and $k_{\Delta}\approx 6$, setting
$\rho(0)N$ randomly chosen nodes in the infectious state. For each setting of
the model parameters, we perform $M=100$ runs, each one using a different
instance of the RSC model, and we evaluate the stationary state $\rho^{*}$ as
the average over the last $100$ values of $\rho(t)$. We run the simulations
with two different values of $\rho(0)$, namely $\rho(0)=0.001$ and
$\rho(0)=0.8$.

Figure 3: Steady-state behavior, $\rho^{*}$, of the SIS model on a Random
Simplicial Complex (RSC) with $N=2000$, $k=20$ and $k_{\Delta}=6$. (a) Phase
diagrams of both individual (black dashed lines) and pair-based (black solid
lines) models, compared to the stochastic simulations in three different cases
of rescaled infectivies on $2$-simplex values $\lambda_{\Delta}$. The blue
dots represent the values of the stationary densities obtained from $M=100$
iterations of the stochastic simulations with $\rho(0)=0.001$, while the
orange dots represent the stationary densities obtained with $\rho(0)=0.8$.
(b) The average fractions of steady-state infected individuals obtained by
means of stochastic simulations on the RSC in the parameter space ($\lambda$,
$\lambda_{\Delta}$). In the upper panel the solid line represent the epidemic
threshold values predicted by the Pair-based model, while in the bottom the
solid line represent the critical lambda values obtained by the Pair-based
model.

We start illustrating our results by comparing the stationary density
$\rho^{*}$ obtained from the stochastic simulations, the individual-based, and
the pair-based model for different sets of parameters. In particular, Fig.
3(a) shows $\rho^{*}$ as a function of the rescaled infectivity $\lambda$ for
three different values of $\lambda_{\Delta}$. The result of the individual-
based model is represented by a dashed black curve, while the result of the
pair-based model is represented by a continuous black curve. The values
obtained from stochastic simulations on the RSC with $\rho(0)=0.001$ are
depicted as blue dots, and those with an initial density of $\rho(0)=0.8$ are
depicted as yellow dots.

First, we observe that both the individual-based and the pair-based model
correctly predict the existence of a region of values of $\lambda$ for which
the system is bistable, and that the transition from the disease-free
equilibrium, $\rho^{*}=0$, to the endemic equilibrium, $\rho^{*}>0$, is
discontinuous. However, the pair-based approximation generally better predicts
the behavior of the stochastic simulations on the RSC compared to the one
provided by the individual-based approximation. In particular, the individual-
based model tends to overestimate the stationary fraction of infected
individuals, while the pair-based model is in better agreement with the
numerical simulations. Our analysis also shows that the pair-based
approximation better identifies the two epidemic thresholds delimiting the
region of bistability, while the individual-based model overestimates the
width of this region. The individual-based model, in fact, underestimates the
first threshold, which represents the critical point for the simulations
started with a high fraction of infected individuals, i.e., $\rho(0)=0.8$,
and, at the same time, overestimates the second threshold, which represents
the critical point for the simulations started with a low fraction of infected
individuals, i.e. $\rho(0)=0.001$. Indeed, while the individual-based
approximation predicts a transition at $\lambda=1$ for any value of
$\lambda_{\Delta}$, the stochastic simulations display a discontinuity at
values of $\lambda<1$ that depends on $\lambda_{\Delta}$, a behavior that is
well reproduced by the pair-based approximation.

Figure 4: Region of bistability of the SIS model and temporal evolution of the
fraction of infected individuals. (a) Phase diagram, in the parameter space
($\lambda$ , $\lambda_{\Delta}$), highlighting the region of bistability
obtained from the individual-based approximation (in orange) and the one
derived from the pair-based approximation (in blue). (b) Temporal evolution of
the fraction of infected individuals $\rho(t)$ for
$(\lambda,\lambda_{\Delta})=(0.3,3.5)$, which corresponds to the point marked
with a black star in (a), for the individual-based (dashed line) and pair-
based (continuous line) approximations. (b) Temporal evolution of the fraction
of infected individuals $\rho(t)$ for $(\lambda,\lambda_{\Delta})=(0.95,3)$,
which corresponds to the point marked with a black star in (a), for the
individual-based (dashed line) and pair-based (continuous line)
approximations. In both panels (b) and (c), for each mean-field model we have
illustrated the temporal evolution from two distinct initial conditions,
namely $\rho(0)=0.8$ and $\rho(0)=0.001$. The colored curves depict $M=100$
realizations of the spreading process on a RSC, with the colors representing
the two different initial conditions of the simulations, namely $\rho(0)=0.8$
(orange) and $\rho(0)=0.001$ (blue).

We then investigate the behavior of $\rho^{*}$ as a function of both $\lambda$
and $\lambda_{\Delta}$, comparing the results of the stochastic simulations
with the value obtained from the pair-based approximation. Fig. 3(b) shows the
results for $\rho(0)=0.8$ (upper panel) and $\rho(0)=0.001$ (lower panel). In
both cases, the color represents the value of $\rho^{*}$ obtained from the
stochastic simulations on the RSC, while the white curves the critical points
derived from the pair-based approximation. Specifically, in the upper panel,
where the pair-based model is initiated with a high value for the initial
condition $\rho(0)$, the white curve corresponds to the first threshold that
delineates the region of bistability and that we here indicate as
$\lambda_{c}$. Instead, in the lower panel, obtained starting from the smaller
value of the initial condition $\rho(0)$, the white curve corresponds to the
second threshold that we denote as $\lambda^{*}$. As shown in Fig. 3(b), both
$\lambda_{c}$ and $\lambda^{*}$ decrease as a function of $\lambda_{\Delta}$,
meaning that the three-body interactions boost the spreading process by
decreasing the thresholds to observe an endemic behavior. In addition, the
distinct rates of decrease observed in $\lambda_{c}$ and $\lambda^{*}$ as a
function of $\lambda_{\Delta}$ suggest that the width of the bistable region
increases with the intensity of the three-body interactions. Such a dependence
of the thresholds on $\lambda_{\Delta}$ is well reproduced by the pair-based
approximation of the simplicial SIS model.

Next, we compare the region of bistability predicted by the individual-based
and pair-based approximations. Fig. 4(a) shows the model behavior as a
function of $\lambda$ and $\lambda_{\Delta}$, with the orange area indicating
the region of bistability predicted by the individual-based approximation and
the blue area that by pair-based approximation. The darker lines represents
the thresholds associated to the two approximations. Finally, the gray area
represents the region that in both models corresponds to a single equilibrium,
either disease-free $\rho^{*}=0$ or endemic $\rho^{*}>0$. Since the pair-based
approximation well predicts the dependence of the critical points on
$\lambda_{\Delta}$, as shown in Fig. 3(b), then, Fig. 4(a) provides another
indication that the individual-based model overestimates the width of the
region of bistability, by simultaneously underestimating $\lambda_{c}$ and
overestimating $\lambda^{*}$. Overall, these results demonstrate that the
pair-based approximation of the simplicial SIS process provides a better
description of the spreading process on simplicial complexes compared to the
individual-based mean-field approximation.

We conclude this section by illustrating in Fig. 4(b)-(c) the time evolution
of $\rho(t)$ obtained from stochastic simulations on the RSC and from the
individual-based and the pair-based approximations for two settings of
$(\lambda,\lambda_{\Delta})$, namely $(\lambda,\lambda_{\Delta})=(0.3,3.5)$
and $(\lambda,\lambda_{\Delta})=(0.95,3)$. These values correspond to the
points of the phase diagram in Fig. 4(a) marked with a star and a diamond,
respectively. In both cases, the individual-based approximation predicts the
existence of two stable equilibriums, while the pair-based approximation a
single one. Specifically, in the pair-based approximation, the steady-state
equilibrium corresponds to a disease-free state for
$(\lambda,\lambda_{\Delta})=(0.3,3.5)$ and to an endemic state for
$(\lambda,\lambda_{\Delta})=(0.95,3)$. As shown in Fig. 4(b)-(c), the time
evolution of $\rho(t)$ calculated from the stochastic simulations converges to
a single equilibrium point, independently from the system initial conditions,
as correctly predicted by the pair-based approximation. The model is also able
to reproduce the average behavior of the temporal evolution of the epidemic
process.

##  V Conclusions

In this work, we have introduced a pair-based mean-field approximation of the
simplicial SIS model. Firstly, we derived the exact equations for the system,
emphasizing the distinctions between our approach and the individual-based
one. We also introduced suitable closures to obtain a closed form for the
model. Then, we compared the behavior of both approximations (the individual-
based and the pair-based) and contrasted the results with those obtained from
stochastic simulations on random simplicial complexes. Our results demonstrate
that the pair-based approach offers a more accurate approximation of the SIS
contagion process when higher-order interactions are present. This is evident
in several aspects, including the extent of the region exhibiting bistability,
the nature of the transition from a disease-free to an endemic state, and the
average time evolution of the fraction of infected individuals.

The pair-based approximation presented in our work is restricted to cases
where higher-order interactions are represented by simplicial complexes.
However, it can be extended to the more general scenario of hypergraphs by
taking into account that there can be three-body interactions not correlated
with two-body interactions. Another possible direction for future work is to
include both dynamical correlations and degree heterogeneity, so that to
extend the range of the hypernetworks whose behavior can be approximated by
these mean-field models.

## References

  * [1] István Z Kiss, Joel C Miller, Péter L Simon, et al.  Mathematics of epidemics on networks.  Cham: Springer, 598, 2017.
  * [2] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani.  Epidemic processes in complex networks.  Reviews of modern physics, 87(3):925, 2015.
  * [3] Jesus Gómez-Gardenes, L Lotero, SN Taraskin, and FJ Pérez-Reche.  Explosive contagion in networks.  Scientific reports, 6(1):19767, 2016.
  * [4] Jae Dong Noh and Heiko Rieger.  Random walks on complex networks.  Physical review letters, 92(11):118701, 2004.
  * [5] Jesús Gómez-Gardenes and Vito Latora.  Entropy rate of diffusion processes on complex networks.  Physical Review E, 78(6):065102, 2008.
  * [6] Stefano Boccaletti, Vito Latora, Yamir Moreno, Martin Chavez, and D-U Hwang.  Complex networks: Structure and dynamics.  Physics reports, 424(4-5):175–308, 2006.
  * [7] Vito Latora, Vincenzo Nicosia, and Giovanni Russo.  Complex networks: principles, methods and applications.  Cambridge University Press, 2017.
  * [8] Federico Battiston, Giulia Cencetti, Iacopo Iacopini, Vito Latora, Maxime Lucas, Alice Patania, Jean-Gabriel Young, and Giovanni Petri.  Networks beyond pairwise interactions: structure and dynamics.  Physics Reports, 874:1–92, 2020.
  * [9] Claude Berge.  Hypergraphs: combinatorics of finite sets, volume 45.  Elsevier, 1984.
  * [10] Ronald H Atkin.  From cohomology in physics to q-connectivity in social science.  International journal of man-machine studies, 4(2):139–167, 1972\.
  * [11] Timoteo Carletti, Duccio Fanelli, and Sara Nicoletti.  Dynamical systems on hypergraphs.  Journal of Physics: Complexity, 1(3):035006, 2020.
  * [12] Timoteo Carletti, Federico Battiston, Giulia Cencetti, and Duccio Fanelli.  Random walks on hypergraphs.  Physical review E, 101(2):022308, 2020.
  * [13] Michael T Schaub, Austin R Benson, Paul Horn, Gabor Lippner, and Ali Jadbabaie.  Random walks on simplicial complexes and the normalized hodge 1-laplacian.  SIAM Review, 62(2):353–391, 2020.
  * [14] Per Sebastian Skardal and Alex Arenas.  Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes.  Physical review letters, 122(24):248301, 2019.
  * [15] Maxime Lucas, Giulia Cencetti, and Federico Battiston.  Multiorder laplacian for synchronization in higher-order networks.  Physical Review Research, 2(3):033410, 2020.
  * [16] Lucia Valentina Gambuzza, Francesca Di Patti, Luca Gallo, Stefano Lepri, Miguel Romance, Regino Criado, Mattia Frasca, Vito Latora, and Stefano Boccaletti.  Stability of synchronization in simplicial complexes.  Nature communications, 12(1):1255, 2021.
  * [17] Unai Alvarez-Rodriguez, Federico Battiston, Guilherme Ferraz de Arruda, Yamir Moreno, Matjaž Perc, and Vito Latora.  Evolutionary dynamics of higher-order interactions in social networks.  Nature Human Behaviour, 5(5):586–595, 2021.
  * [18] Damon Centola and Michael Macy.  Complex contagions and the weakness of long ties.  American journal of Sociology, 113(3):702–734, 2007.
  * [19] Guilherme Ferraz de Arruda, Giovanni Petri, and Yamir Moreno.  Social contagion models on hypergraphs.  Physical Review Research, 2(2):023032, 2020.
  * [20] Iacopo Iacopini, Giovanni Petri, Alain Barrat, and Vito Latora.  Simplicial models of social contagion.  Nature communications, 10(1):1–9, 2019.
  * [21] Federico Malizia, Santiago Lamata-Otín, Mattia Frasca, Vito Latora, and Jesús Gómez-Gardeñes.  Hyperedge overlap drives explosive collective behaviors in systems with higher-order interactions.  arXiv preprint arXiv:2307.03519, 2023.
  * [22] Matthew J Keeling.  The effects of local spatial structure on epidemiological invasions.  Proceedings of the Royal Society of London. Series B: Biological Sciences, 266(1421):859–867, 1999.
  * [23] Kieran J Sharkey, Istvan Z Kiss, Robert R Wilkinson, and Peter L Simon.  Exact equations for sir epidemics on tree graphs.  Bulletin of mathematical biology, 77:614–645, 2015.
  * [24] Mattia Frasca and Kieran J Sharkey.  Discrete-time moment closure models for epidemic spreading in populations of interacting individuals.  Journal of theoretical biology, 399:13–21, 2016.
  * [25] Federico Malizia, Luca Gallo, Mattia Frasca, Vito Latora, and Giovanni Russo.  Individual-and pair-based models of epidemic spreading: Master equations and analysis of their forecasting capabilities.  Physical Review Research, 4(2):023145, 2022.
  * [26] Joan T Matamalas, Sergio Gómez, and Alex Arenas.  Abrupt phase transition of epidemic spreading in simplicial complexes.  Physical Review Research, 2(1):012049, 2020.
  * [27] Giulio Burgio, Sergio Gómez, and Alex Arenas.  Triadic approximation for contagions on higher-order networks.  arXiv preprint arXiv:2306.11441, 2023.
  * [28] John G Kirkwood.  Statistical mechanics of fluid mixtures.  The Journal of chemical physics, 3(5):300–313, 1935.
  * [29] The formula can be understood as follows: $P\left(AB\cap BC\cap CA\right)=P\left(AB\cap BC|CA\right)P\left(CA\right)$. Now $P\left(AB\cap BC|CA\right)=P\left(AB|BC\cap CA\right)\cdot P\left(BC|CA\right)$. We can write $P\left(AB|BC\cap CA\right)=\genfrac{}{}{}{0}{P\left(AB\right)}{P\left(A\right)P\left(B\right)}$, and use the approximation $P\left(AB|BC\cap CA\right)=\genfrac{}{}{}{0}{P\left(AB\right)}{P\left(A\right)P\left(B\right)}$. Substituting in the previous expression one has $P(AB\cap BC|CA)=\genfrac{}{}{}{0}{P\left(AB\right)}{P\left(A\right)P\left(B\right)}\cdot\genfrac{}{}{}{0}{P\left(BC\right)}{P\left(C\right)}$ from which it follows $P\left(AB\cap BC\cap CA\right)={P(AB)P(BC)P(CA)}/{(P(A)P(B)P(C))}$.
  * [30] Thomas House, Geoffrey Davies, Leon Danon, and Matt J Keeling.  A motif-based approach to network epidemics.  Bulletin of Mathematical Biology, 71(7):1693–1706, 2009.