# Theoretical and numerical considerations of the assumptions behind triple
closures in epidemic models on networks

Nicos Georgiou  Nicos Georgiou
University of Sussex
Department of Mathematics
Falmer Campus
Brighton BN1 9QH
UK<EMAIL_ADDRESS>http://www.sussex.ac.uk/profiles/329373 ,  István
Z. Kiss  Istvan Kiss
University of Sussex
Department of Mathematics
Falmer Campus
Brighton BN1 9QH
UK<EMAIL_ADDRESS>http://www.sussex.ac.uk/profiles/211073 and  Péter
Simon  Péter Simon
Institute of Mathematics, Eötvös Loránd University Budapest, and Numerical
Analysis and Large Networks Research Group, Hungarian Academy of Sciences
Hungary<EMAIL_ADDRESS>http://simonp.web.elte.hu/index_en.html

###### Abstract.

Networks are widely used to model the contact structure within a population
and in the resulting models of disease spread. While networks provide a high
degree of realism, the analysis of the exact model is out of reach and even
numerical methods fail for modest network size. Hence, mean-field models (e.g.
pairwise) focusing on describing the evolution of some summary statistics from
the exact model gained a lot of traction over the last few decades. In this
paper we revisit the problem of deriving triple closures for pairwise models
and we investigate in detail the assumptions behind some of the well-known
closures as well as their validity. Using a top-down approach we start at the
level of the entire graph and work down to the level of triples and combine
this with information around nodes and pairs. We use our approach to derive
many of the existing closures and propose new ones and theoretically connect
the two well-studied models of multinomial link and Poisson link selection.
The theoretical work is backed up by numerical examples to highlight where the
commonly used assumptions may fail and provide some recommendations for how to
choose the most appropriate closure when using graphs with no or modest degree
heterogeneity.

###### Key words and phrases:

closures, multinomial link model, epidemic models, SIS

###### 2000 Mathematics Subject Classification:

92D30, 90B10, 90B15

NG was partially supported by the EPSRC First Grant EP/P021409/1. IZK was
partially supported by the Leverhulme Trust, Research Project Grant
RPG-2017-370. PS was partially supported by the Hungarian Scientific Research
Fund, OTKA, Grant no. 115926.

##  1\. Introduction

Many complex systems are forbiddingly high dimensional and one efficient way
to deal with such a challenge is to focus on some summary statistics or
moments of the full-system. However, more often than not, the evolution of a
moment (usually a coarse grained view or quantity that is computable from the
full model) depends on higher-order moments. In [5], the author summarises the
four main steps of a typical moment closure based approach. These are: (a)
select the moments and their hierarchy, (b) write down evolution equations for
the moments, (c) derive and apply the moment-closure and, finally, (d) justify
and validate the choice of moments and closures and perform further numerical
tests.

Unfortunately, there is no single or precise way in which moment hierarchies
and moment closures are derived. Often these rely on empirical or numerical
observations, may only work for particular problems or in specific contexts
and rigorous mathematical proofs for closures are rare. As a result in this
paper, we revisit some of the existing closures in the context of epidemic
models and show that these can be derived by a hybrid method combining a top-
down approach (counting over the entire graph/network) and a bottom-up
approach where assumptions about the states of the neighbours of nodes are
made.

This is a well-know and well-studied area with some initial important results
by [6, 3]. In particular they showed that closures for some epidemic models
are possible by starting at the level of nodes and considering their degrees
as a random variable with a given mean plus some random error with zero mean.
Coupling this with knowledge at the graph-level about the counts of nodes,
links and triples in various states, they managed to derive two distinct
closures which we will investigate in detail. One of the main assumption made
in their approach is about the distribution of the states of nodes around a
given node (i.e. Multinomial or Poisson).

In what follows we start from the entire graph and aim to re-derive some of
the existing closures and propose two new ones, but with the main contribution
being that we provide a rather general approach and we illustrate it by using
a hypothetical (i.e. uniform) form of the distribution of the sates of the
nodes around a give node.

The paper is structured as follows. In Section 2 we give the main ingredients
of the model including the network and epidemic dynamics on it together with
the unclosed pairwise model for SIS dynamics. In Section 3 we present the top-
down derivation of an exact expression for the expected value of $k$-tuples of
states in the epidemics via probabilistic considerations. In Section 4 we re-
derive some of the classical closures used in the literature and present two
new ones, that performs well in numerical test. We also highlight when the
Multinomial and Poisson assumptions work and when they are hard to
distinguish. The final part of the section explains how to rigorously derive
the Poisson model from the Multinomial link model, both widely used in the
literature. In Section 5 we provide a large number of numerical tests of the
closures and show that the pairwise model with the new closure leads to an
epidemic threshold which is well-known in the literature. Further results and
applications of our method are presented in Section 6. Finally, Section 7 is
dedicated to a discussion and summary of our findings.

##  2\. Model: network, epidemic dynamics and mean-field models

The starting point is to model the contact structure of a population of $N$
individuals as an undirected and unweighted network
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, so that $|\mathcal{V}|=N$. Such a
graph can also be represented by adjacency matrix
$G=(g_{ij})_{i,j=1,2,\dots,N}$, where $g_{ij}=1$ if nodes $i$ and $j$ are
connected and zero otherwise. Self-loops are excluded, so $g_{ii}=0$ and
$g_{ij}=g_{ji}$ for all $i,j=1,2,\dots N$. While we will consider a general
epidemic model where nodes can be in an arbitrary number of discrete states
(i.e. $m$ different states $A_{1},A_{2},\ldots A_{m}$), our examples will
focus on the standard susceptible-infected-susceptible (SIS) epidemic dynamics
on a network. The SIS dynamics is driven by two processes: (a) infection and
(b) recovery from infection. Infection spreads from an infected node (I) to
any of its susceptible neighbours (S) and this is modelled as a Poisson point
process with per-link infection rate $\tau$. Infectious nodes recover from
infection at constant rate $\gamma$, independently of the network, and become
susceptible again. The resulting model is a continuous-time Markov Chain over
a state space with $2^{N}$ elements. This consists of all arrangements of
length $N$ with each entry being either S or I independently. While this is
easy to generalise and write down theoretically, the numerical integration of
the system becomes intractable even for modest values of $N$ [4].

One way to deal with such a high-dimensional model is to derive mean-field
approximations for some of the summary statistics of the exact process.
Probably of most interest is the expected number of infected nodes over time.
There are many different approaches that can help achieve this [4], but almost
all rely on starting at ‘node’ or ‘node and its neighbourhood-level’ and
proceed by writing down differential equations for their evolution. This
immediately leads to a dependency on higher-order moments since whatever
quantity we focus on, its evolution will depend on the type and states of the
neighbouring nodes, e.g. singles depend on pairs and pairs depend on triple.
One straightforward mean-field model, also extensively used in this paper, is
the pairwise model which is given below,

(2.1) $\displaystyle\dot{[S]}$ $\displaystyle=-\tau[SI]+\gamma[I],$ (2.2)
$\displaystyle\dot{[I]}$ $\displaystyle=\tau[SI]-\gamma[I],$ (2.3)
$\displaystyle\dot{[SS]}$ $\displaystyle=-2\tau[SSI]+2\gamma[SI],$ (2.4)
$\displaystyle\dot{[SI]}$
$\displaystyle=\tau([SSI]-[ISI]-[SI])-\gamma([SI]-[II]),$ (2.5)
$\displaystyle\dot{[II]}$ $\displaystyle=2\tau([ISI]+[SI])-2\gamma[II].$

Here, $[\cdot]$ stands for the expected number of the respective quantities.
This unclosed model was derived from the full system of master equations and
it was proved to be exact [8].

The equations above are straightforward to interpret [4]. The evolution
equation $[I]$, see (2.2), has two terms: (a) a positive term which is
proportional to the expected number of S-I links ($[SI]$) and represents the
incoming flux of new infections, and (b) a negative term which is proportional
to $[I]$ and stands for the recovery of infected nodes. More importantly, we
notice that singles depend on pairs and this hierarchy of dependency continues
whereby $[SI]$ links are created at rate $\tau[SSI]$, meaning that one of the
S nodes in an S-S link can be infected by an external node, e.g. S-S-I.
Equally, the $[SI]$ links are depleted due to events within pair, i.e. I can
infect S or I can recover, or events from outside the pair, i.e. an external
node infecting S, that is I-S-I. It is now clear that there is hierarchy of
dependency on ever higher moments.

Keeping the hierarchy would lead to a large number of more and more
complicated equations. But what if higher order moments, say triples, can be
approximated by lower-order ones, such as singles and pairs? This process is
referred to as ‘closure’ and it is fundamental to derive mean-field models. In
what follows we start from the entire graph and aim to re-derive some of the
existing closures and propose new ones, but with special focus on staring at
the graph-level and combining this with information at the local or node-
level.

##  3\. Triple counts from neighbourhood distribution

Deriving evolution equation for some of the summary statistics from the exact
model will require to find reasonable approximations to expressions of the
form

$[A_{i_{1}}A_{i_{2}}\ldots A_{i_{k}}],\text{ where }k\in\mathbb{N}\text{ and
}i_{j}\in\\{1,\ldots,m\\}\text{ for all }j.$

All quantities depend on time, so we assume it fixed and do not encumber
notation with an extra $t$ index. Also keep $k$ fixed for the moment. Define
the set of paths of size $k$ as

$\Pi_{k}=\\{(v_{1},v_{2},\ldots,v_{k}):(v_{1},v_{2},\ldots,v_{k})\in\mathcal{V}^{k},v_{i}\neq
v_{j},(v_{i},v_{i+1})\in\mathcal{E}\text{ for }1\leq i\leq k-1\\}.$

Admissible and non-admissible paths are shown in Figure 1. Irrespective of the
Markov chain which happens on the network $\mathcal{G}$, the process will
induce a probability mass function on the $k$-tuples of $\Pi_{k}$. At a given
time every element of the state space has a well-defined probability mass
which is given by the Markov chain itself, and obviously depends also on the
initial state of the system. Then the probabilities of the states determine
the expected value of the elements of $\Pi_{k}$.

$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$ Figure
1. Admissible (left) and non-admissible (right) ordered network paths going
through 5 nodes. The left path visits each of its 5 nodes exactly once. The
right one is inadmissible in our calculations, as the second and fourth node
on it coincide.

An expression of the expected value of the $k-$tuple can be computed using the
law of total probability, by conditioning on the type of site $v_{j}$:

$\displaystyle\mathbb{E}(A_{i_{1}}$ $\displaystyle A_{i_{2}}\ldots
A_{i_{k}})=\mathbb{E}\Big{(}\sum_{(v_{1},\ldots,v_{k})\in\Pi_{k}}\mathbbm{1}\\{v_{1}=A_{i_{1}},v_{2}=A_{i_{2}},\ldots,v_{k}=A_{i_{k}}\\}\Big{)}$
(3.1)
$\displaystyle=\sum_{(v_{1},\ldots,v_{k})\in\Pi_{k}}\mathbb{P}\\{v_{1}=A_{i_{1}},v_{2}=A_{i_{2}},\ldots,v_{k}=A_{i_{k}}\\}$
(3.2)
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A_{i_{j}}\\}\\!\\!\\!\\!\sum_{(v_{1},\ldots,v_{j}=v,\ldots,v_{k})\in\Pi_{k}}\\!\\!\\!\\!P\\{v_{1}=A_{i_{1}},\ldots,v_{k}=A_{i_{k}}|v_{j}=v=A_{i_{j}}\\}.$

This expression is based on the fact that we are computing the expectation of
a specific arrangement of states on a non-self-intersecting paths. Different
expressions involving conditional quantities for expected values of more
complicated sets of $k$-tuples will hold in the case of non-admissible paths;
the applications below are for $k=3$ involving open triples and therefore we
do not need to consider more complicated paths.

For the theoretical part of this article, we consider a multi-type epidemic on
$\mathcal{G}$ and assume a node can be in any of $m$ different states, say
$A_{1},\ldots,A_{m}$. We will first compute theoretical expressions for
$[A_{i}AA_{j}]$ using (3.2).

For a given node $v$ denote by

(3.3) $N_{i}^{(v)}=\text{card}\\{\text{neighbours of $v$ that take the value
$A_{i}$}\\},\quad A_{i}\in\\{A_{1},\ldots,A_{m}\\}.$

Keep in mind that $\displaystyle\sum_{i=1}^{m}N_{i}^{(v)}=\deg v$.

######  Proposition 3.1.

Consider a multi-type epidemic on $\mathcal{G}$ and assume a node can be on
any of $m$ different states $A_{1},\ldots,A_{m}$. Then, the following two
formulas hold:

For $A_{i}\neq A_{j}$

(3.4)
$[A_{i}AA_{j}]=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}N_{i}^{(v)}N_{j}^{(v)}\Big{|}v=A\Big{]},$

for $A_{i}=A_{j}$

(3.5)
$[A_{i}AA_{i}]=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}N_{i}^{(v)}(N_{i}^{(v)}-1)\Big{|}v=A\Big{]}.$

###### Proof of proposition 3.1.

We first assume that the types $A_{i}\neq A_{j}$, and we start computing from
equation (3.2):

$\displaystyle[A_{i}AA_{j}]=\mathbb{E}(A_{i}AA_{j})$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}\sum_{(v_{1},v,v_{3})\in\Pi_{3}}P\\{v_{1}=A_{i},v_{3}=A_{j}|v=A\\}\text{
by \eqref{eq:condmass}}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}\sum_{v_{1}:v_{1}\sim
v}\sum_{v_{3}:v_{3}\sim
v}\mathbbm{1}\\{v_{1}=A_{i},v_{3}=A_{j}\\}\Big{|}v=A\Big{]}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}N_{i}^{(v)}N_{j}^{(v)}\Big{|}v=A\Big{]}.$

Similarly, we repeat the computation for when $A_{i}=A_{j}$. The change is
during the second equality above, when we compute the double sum.

$\displaystyle[A_{i}AA_{i}]=\mathbb{E}(A_{i}AA_{i})$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}\sum_{(v_{1},v,v_{3})\in\Pi_{3}}P\\{v_{1}=A_{i},v_{3}=A_{i}|v=A\\}\text{
by \eqref{eq:condmass}}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}\sum_{v_{1}:v_{1}\sim
v}\sum_{v_{3}:v_{3}\sim v,v_{3}\neq
v_{1}}\mathbbm{1}\\{v_{1}=A_{i},v_{3}=A_{i}\\}\Big{|}v=A\Big{]}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}\sum_{v_{1}:v_{1}\sim
v}\mathbbm{1}\\{v_{1}=A_{i}\\}\sum_{v_{3}:v_{3}\sim v,v_{3}\neq
v_{1}}\mathbbm{1}\\{v_{3}=A_{i}\\}\Big{|}v=A\Big{]}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}N_{i}^{(v)}(N_{i}^{(v)}-1)\Big{|}v=A\Big{]}.$

∎

In fact fixing an arbitrary network and disease parameters allows us to
extract the expected values of triples that drive the epidemic (e.g. for an
SIS epidemic these are $[SIS](t)$ and $[ISI](t)$) directly from sampling
multiple realisations of the stochastic epidemic. Sampling across a discrete
time set can then give the empirical conditional distribution of

$(N_{1}^{(v)},N_{2}^{(v)},\ldots,N_{m}^{(v)})|v=A$

and also a way to estimate $\mathbb{P}\\{v=A\\}$ on these times, which can
then be used in Proposition 3.1. Networks with enough symmetries (i.e. if
nodes are exchangeable) and a large number of nodes can give a good
approximation of the empirical distribution with just one simulation.

These expected triple counts in time can then be fed into the pairwise model
(2.1)-(2.5) and hence its numerical integration is possible without a closure.
In fact, this leads to a very accurate mean-field model as shown in Figure 2
for three different network models. Of course, the applicability of this
method is limited as it relies on output from the simulation. Nevertheless, it
shows that being able to accurately approximate triples is key to derive
accurate mean-field models.

Figure 2. Individual realisations of the stochastic SIS model (multicoloured
noisy lines) on regular (left), Erdős-Rényi (middle) and bi-modal (right)
networks. The average of all realisations is also shown (thick continuous
line) along with the solution of the pairwise model (2.1)-(2.5) (thin dashed
line) with the triples $[SSI]$ and $[ISI]$ taken as the expected value
computed across all 200 individual epidemic simulations. In all figures
$N=200$, $\tau=0.5$ and $\gamma=1$. Regular and Erdős-Rényi networks have
$\langle k\rangle=4$ while the bi-modal network has half of the nodes with
degree 4 and the the other half with degree 8. Each epidemic starts with 10
infected nodes chosen uniformly at random and only epidemics that reach 20
infected nodes count towards the average with time reset to $t=0$ when the
state with twenty infected nodes is first reached.

######  Remark 3.2.

Note that the formulas in Proposition 3.1 have the sum over all vertices of
$\mathcal{G}$. This is just the convenient way to write the formula for the
sequel; the terms that need to be estimated are then the $\mathbb{P}\\{v=A\\}$
and $E\Big{[}N_{i}^{(v)}N_{j}^{(v)}\Big{|}v=A\Big{]}$. We explain in the next
sections how these quantities can be approximated.

The conditional expectation indicates that the important quantities clearly
depend on the state of $v$ being $A$ and, per realisation, it is those sites
that guide the value of the triple. Here is a different way that highlights
this, assuming $A_{i}\neq A_{j}$, we get:

$\displaystyle\\#\\{A_{i}AA_{j}\\}$
$\displaystyle=\sum_{v\in\mathcal{G}}\sum_{w:w\sim v}\sum_{z:z\sim
v}\mathbbm{1}\\{v=A\\}\mathbbm{1}\\{w=A_{i}\\}\mathbbm{1}\\{z=A_{j}\\}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbbm{1}\\{v=A\\}N_{i}^{(v)}N_{j}^{(v)}=\sum_{v:v=A}N_{i}^{(v)}N_{j}^{(v)}.$

This is an exact count, and the sum goes through the (random) set of nodes in
state A. What we need for the mean-field model in this article is the expected
value. Taking expectations cannot simplify the right-hand side any further as
the sum is over a random index and as such it does not commute with the
expectation operation. This is the reason why we opt to use the formulas in
Proposition 3.1 rather than this more intuitive one.

##  4\. Closures based on the distribution of states of the neighbours

Proposition 3.1 can be useful in certain situations, as can be seen in the
following cases. In all that follow, a theoretical ansatz is made on the
distribution of $(N_{1}^{(v)},N_{2}^{(v)},\ldots,N_{m}^{(v)})|v=A$ that can be
used to further develop the formulas in Proposition 3.1. As a consequence, we
can recover directly several closure formulas that are used in the literature,
and showcase some further examples.

###  4.1. Closures for the multinomial distribution of states of the
neighbours

Consider a a multi-type epidemic on a network, with possible states
$A_{1},\ldots,A_{m}$. For the Multinomial link model we assume that for any
given vertex $v$, the distribution of the $A_{i}$’ s of its neighbouring
vertices depends on the state of $v$, (say $v=A$) and conditional on the state
of $v$, each neighbouring vertex takes a value $A_{i}$ from
$\\{A_{1},\ldots,A_{m}\\}$ with probability $p^{(v,A)}_{i}$, independently of
everything else. The only condition required at this point is that

(4.1) $\sum_{i=1}^{m}p^{(v,A)}_{i}=1.$

This description is equivalent to the following

######  Assumption 4.1.

Conditional on $v=A$, the conditional distribution

(4.2) $(N_{1}^{(v)},N_{2}^{(v)},\ldots,N_{m}^{(v)})|v=A\sim{\rm Mult}(\deg
v,m;p^{(v,A)}_{1},\ldots,p^{(v,A)}_{m}),$

i.e. the distribution of types around a node of type $A$ is multinomial with
$m$ possible outcomes of $\text{deg}\,v$ independent experiments.

######  Remark 4.2.

When $\deg v$ is constant (and not random), the marginal distributions of
$N_{i}^{(v)}$ are binomial and they were considered before. They were used to
derive a triple closure for simple epidemics (see [6], chapter 4, p.104). In
that notation [6], this distribution is denoted by $Q_{v}(i|j)$. Correlations
were introduced to take into account the fact that the degree creates
constraints, turning it into the multinomial model. The authors also discuss
Poisson closures as viable approximations and we also do that at the end of
the section, re-deriving their formulas.

As we mentioned, Assumption 4.1 implies that the (conditional) marginal
distributions of the coordinates are binomial and we can immediately obtain
that

$E(N_{i}^{(v)}|v=A)=p_{i}^{(v,A)}\deg v,\text{ for all }i\in\\{1,\ldots,m\\}.$

Moreover, a direct computation using (4.2) gives

$E\Big{[}N_{i}^{(v)}N_{j}^{(v)}\Big{|}v=A\Big{]}=p_{i}^{(v,A)}p_{j}^{(v,A)}\deg
v(\deg v-1).$

The conditional distribution is multinomial, therefore the conditional
expectation is the expectation of products of binomially distributed
coordinates of the multinomial distribution (4.2). A different way to write
$E[N_{i}^{(v)}N_{j}^{(v)}|v=A]=\text{Cov}_{v=A}(N_{i}^{(v)},N_{j}^{(v)})+E[N_{i}^{(v)}|v=A]\,E[N_{j}^{(v)}|v=A]$,
which corresponds to the approach in [6]. Therefore when we substitute in
(3.4) for the case where $A_{i}\neq A_{j}$ we get that

(4.3) $\displaystyle[A_{i}AA_{j}]$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}p_{i}^{(v,A)}p_{j}^{(v,A)}\deg
v(\deg v-1).$

The final step is to derive as accurate as possible approximations for the
probabilities/quantities in (4.3).

######  Remark 4.3.

Consider SIS epidemics on a regular network of degree $n$. There are only two
types, so the distribution of types around (say) S nodes is completely
characterised by the number (and therefore the distribution) of type I around
S nodes. If the multinomial model assumption is accurate, then that
distribution of I’s should be Binomial with $n$ trials. Indeed this can be
numerically verified. Figure 3 shows histograms of the marginal distributions
of the number of I nodes around susceptible nodes for an SIS epidemic
simulated on a regular network. The probability of the binomial distribution
is estimated using formula (4.4), and it fits well to the simulated data. This
particular fit suggests that Assumption 4.1 is not unrealistic in many
situations.

Figure 3. The distribution of infected nodes around susceptible nodes at
different points in the epidemic on a Regular network with $N=200$ nodes with
$\langle k\rangle=6$, $\gamma=1$ and $\tau=0.75$. This is base on one single
realisation of the epidemic on one single realisation of the network. The
continuous and dashed lines stand for the Binomial (with $\langle k\rangle=6$
trials and probability of success $[SI]/(\langle k\rangle[S])$ ) and Poisson
distributions respectively (with rate $[SI]/[S]$), where the pair and single
are taken from simulation. Notice how the binomial mass function is a better
fit for this data than the Poisson one.

At this point we carefully explain the remaining assumptions that lead to the
closures we are concerned with.

######  Assumption 4.4.

$\mathbb{P}\\{v=A\\}\sim\frac{[A]}{N},\quad{\rm
for\,\,\,all\quad}v\in\mathcal{G}.$

We emphasise that this is not an assumption that applies to all networks, as
it suggests that at any given time $t$, that any node has the same probability
of being (say) infected; intuitively this fails when in networks with high
degree heterogeneity, as well-connected or high-degree nodes have higher
probability of being in a certain state. This is a reasonable assumption on
networks with low degree heterogeneity.

The second assumption is about the values of $p_{i}^{(v,A)}$ in the
conditional multinomial distribution.

######  Assumption 4.5.

We assume that

$p_{i}^{(v,A)}\equiv p_{i}^{(A)}.$

In other words, the parameters of the multinomial distribution do not depend
on $v$ in any way, but only on the state $v$ is in.

Under Assumptions 4.1, 4.4 and 4.5 we can first derive approximations for
$p_{i}^{(A)}$:

$\displaystyle[A_{i}A]$
$\displaystyle=\mathbb{E}(A_{i}A)=\mathbb{E}\Big{(}\sum_{v\in\mathcal{G}}\sum_{w:w\sim
v}\mathbbm{1}\\{w=A_{i},v=A\\}\Big{)}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E(N_{i}^{(v)}|v=A)\approx\frac{[A]}{N}\sum_{v\in\mathcal{G}}\text{deg}\,v\,p_{i}^{(A)}$
$\displaystyle=[A]\,p_{i}^{(A)}\langle\text{deg}\,v\rangle.$

Above we defined $\langle\text{deg}\,v\rangle$ as the average degree of the
network. The above computation implies

(4.4)
$p_{i}^{(A)}\approx\frac{[A_{i}A]}{\langle\text{deg}\,v\rangle[A]},\text{ for
all }i\in\\{1,2,\ldots,m\\}.$

We can now conclude computation (4.3) under these approximations:

$\displaystyle[A_{i}AA_{j}]$
$\displaystyle\approx\sum_{v\in\mathcal{G}}\frac{(\text{deg}v-1)\text{deg}v}{\langle\text{deg}v\rangle^{2}}\frac{1}{N}\frac{[A_{i}A][AA_{j}]}{[A]}$
$\displaystyle=\frac{[A_{i}A][AA_{j}]}{[A]\langle\text{deg}\,v\rangle^{2}}\Big{(}\frac{1}{N}\sum_{v\in\mathcal{G}}(\text{deg}v-1)\text{deg}v\Big{)}$
$\displaystyle=\frac{[A_{i}A][AA_{j}]}{[A]\langle\text{deg}\,v\rangle^{2}}\langle\deg
v(\deg v-1)\rangle.$

Now in the case where $A_{i}=A_{j}$ we use (3.5) and have

$\displaystyle[A_{i}AA_{i}]$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}(N_{i}^{(v)})^{2}-N_{i}^{(v)}\Big{|}v=A\Big{]}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}\Big{(}\text{Var}(N_{i}^{(v)}|v=A)+(\mathbb{E}N_{i}^{(v)}|v=A)^{2}-\mathbb{E}(N_{i}^{(v)}|v=A)\Big{)}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}(p_{i}^{(A}(1-p_{i}^{(A)})\deg
v+(p_{i}^{(A)}\deg v)^{2}-p_{i}^{(A)}\deg v)$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}(p_{i}^{(A)})^{2}\deg
v(\deg v-1).$

Substitute in the approximations for $\mathbb{P}\\{v=A\\}$ and $p_{i}^{(A)}$
to obtain

$[A_{i}AA_{i}]\approx\frac{[A]}{[N]}\frac{[A_{i}A]^{2}}{\langle\deg
v\rangle^{2}[A]^{2}}\sum_{v\in\mathcal{G}}\deg v(\deg
v-1)=\frac{[A_{i}A]^{2}}{[A]\langle\deg v\rangle^{2}}\langle\deg v(\deg
v-1)\rangle.$

Thus starting from Proposition 3.1 and using the above assumptions we have
proved the following.

######  Theorem 4.6.

Under Assumptions 4.1, 4.4 and 4.5 the expected value of triples can be given
as

(4.5)
$[A_{i}AA_{j}]=\frac{[A_{i}A][AA_{j}]}{[A]\langle\text{deg}\,v\rangle^{2}}\langle\deg
v(\deg v-1)\rangle$

and

(4.6) $[A_{i}AA_{i}]=\frac{[A_{i}A]^{2}}{[A]\langle\deg
v\rangle^{2}}\langle\deg v(\deg v-1)\rangle.$

Applying this theorem to SIS dynamics we obtain

######  Corollary 4.7.

For SIS epidemics, closures (4.5), (4.6) become

(4.7) $[SSI]\approx\frac{[SS][SI]}{[S]\langle\deg\,v\rangle^{2}}\langle\deg
v(\deg v-1)\rangle,$

and

(4.8) $[ISI]\approx\frac{[SI]^{2}}{[S]\langle\deg\,v\rangle^{2}}\langle\deg
v(\deg v-1)\rangle.$

######  Remark 4.8.

The formulas in Corollary 4.7 are strikingly similar to the ones in article
[7], which was obtained by decomposing triples based on the degree of the
susceptible nodes. However, there is an a-priori significant difference, as
the degree moments here are over all nodes in the network, not just those in
state S at time $t$. If the nodes in state S at any given time in the network
is a statistically significant sample of nodes in the network and the nodes
are exchangeable (e.g. low-degree heterogeneity), the moments here and the
moments used in [7] are very close to each other and the two closures will
perform in a similar way.

######  Example 1.

($n$-regular networks $\mathcal{G}$) If $\mathcal{G}$ is $n$-regular, i.e. all
nodes have degree $n$, then $\deg v=n=\langle\deg v\rangle$. Assume for the
moment that Assumption 4.5 does not necessarily hold. We may still proceed
with the computation of $[A_{i}A]$, but now we would have

$[A_{i}A]\approx\frac{[A]}{N}\sum_{v\in\mathcal{G}}p_{i}^{(v,A)}\deg v,$

and it can be pulled out of the sum, and after an elementary algebraic
manipulation we have that

$\frac{1}{N}\sum_{v\in\mathcal{G}}p_{i}^{(v,A)}\approx\frac{[A_{i}A]}{n[A]}.$

Therefore the statistical average of $p_{i}^{(v,A)}$ (on the left hand side)
is only a function of the state (right-hand side) which coincides with the
expression in (4.4).

For $\mathcal{G}$ $n$-regular, equation (4.5) becomes the classical binomial
closure

(4.9) $[A_{i}AA_{j}]\approx\frac{n-1}{n}\frac{[A_{i}A][AA_{j}]}{[A]},$

which can be found in [6, 3].

Equation (4.4) is only an approximation. As it is stated, it violates the
consistency condition (4.1). That can be easily accounted for by defining a
normalised version. There are several ways to do this, each one leading to a
different closure. We conclude this section with one such closure for SIS
epidemics. Since there are only two possible states one way to normalise
probabilities is to assume that

(4.10) $p^{(S)}_{I}=\frac{[SI]}{[S]\langle\deg v\rangle}\text{ and
}p^{(S)}_{S}=1-p^{(S)}_{I}.$

With this, we have the following.

######  Corollary 4.9.

Under approximation (4.10)

(4.11) $[SSI]\approx[S]p^{(S)}_{I}(1-p^{(S)}_{I})\langle\deg v(\deg
v-1)\rangle,\quad[ISI]\approx[S](p^{(S)}_{I})^{2}\langle\deg v(\deg
v-1)\rangle.$

Note that while we may also assume

(4.12) $p^{(S)}_{S}=\frac{[SS]}{[S]\langle\deg v\rangle}\text{ and
}p^{(S)}_{I}=1-p^{(S)}_{S},$

this and the resulting closure are not driven by the important quantity in the
epidemics, namely the number of $[SI]$ links since the quantity does not show
up in the closure. As such, it is expected that a normalised closure using
(4.12) would perform badly. Indeed this is supported by numerical evidence.
Furthermore, we also tried to normalise by simply dividing each ‘probability’
with the sum of the two but this led to an ill-behaved closed pairwise model
and thus we omitted it here.

###  4.2. Closures for the Poisson distribution of states of the neighbours

In the regime where $n$ is large and the epidemic level is not high, we assume
that conditional on $v=A$ the collection $\\{N_{i}^{(v)}|v=A\\}_{1\leq i\leq
m-1}$ is a collection of independent Poisson random variables; one type is the
prevalent one which we omitted from the collection - for example in the
initial stages of SIR epidemic, we can use this to model $N_{R}$ and $N_{I}$
since $S$ is prevalent. We assume that the parameters of the Poisson
distribution do not depend on specific $v$ and we formalise this as follows.

######  Assumption 4.10.

(4.13) $N_{i}^{(v)}|v=A\sim{\rm Poisson}(\lambda_{i}^{(A)}).$

For example if $\lambda_{i}^{(v,A)}<n^{\alpha}<<n$ ($\alpha<1$), the
independence assumption is not restrictive, since Poisson random variables are
sharply concentrated around their mean, and the error of the approximation
below can be controlled for types $1\leq i\neq j\leq m-1$. In Figure 4 we show
an example where the Poisson distribution does fit well the marginal of the
number of infected nodes around susceptible ones. The independence assumption
is harder to verify numerically, but we present at the end of this section how
the Poisson model naturally arises in Erdős-Rényi graphs with low connection
probability and prove the independence assumption.

Figure 4. The distribution of infected nodes around susceptible nodes at
different points in the epidemic on a Erdős-Rényi network with $N=200$ nodes
with $\langle k\rangle=6$, $\gamma=1$ and $\tau=0.25$. Maximum degree in the
network is twelve. This is based on one single realisation of the epidemic on
one single realisation of the network. The dashed line line stands for the
Poisson distribution with mean given by $[SI]/[S]$ taken from simulation.

Assumption 4.10 implies that

$E(N_{i}^{(v)}|v=A)=\lambda_{i}^{(A)},\text{ for all }i\in\\{1,\ldots,m\\}.$

Moreover, the independence assumption on the coordinates implies the equality
below,

$E\Big{[}N_{i}^{(v)}N_{j}^{(v)}\Big{|}v=A\Big{]}=E\Big{[}N_{i}^{(v)}\Big{|}v=A\Big{]}E\Big{[}N_{j}^{(v)}\Big{|}v=A\Big{]}=\lambda_{i}^{(A)}\lambda_{j}^{(A)}.$

Under Assumption 4.10 we can first derive an approximation for
$\lambda_{i}^{(A)}$:

$\displaystyle[A_{i}A]$
$\displaystyle=\mathbb{E}(A_{i}A)=\mathbb{E}\Big{(}\sum_{v\in\mathcal{G}}\sum_{w:w\sim
v}\mathbbm{1}\\{w=A_{i},v=A\\}\Big{)}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E(N_{i}^{(v)}|v=A)=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}\lambda_{i}^{(A)}=[A]\,\lambda_{i}^{(A)}.$

Thus we have

(4.14) $\lambda_{i}^{(A)}\approx\frac{[A_{i}A]}{[A]},\text{ for all
}i\in\\{1,2,\ldots,m\\}.$

Starting again from (3.4) we obtain,

$[A_{i}AA_{j}]=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}N_{i}^{(v)}N_{j}^{(v)}\Big{|}v=A\Big{]}=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}\lambda_{i}^{(A)}\lambda_{j}^{(A)}\approx[A]\lambda_{i}^{(A)}\lambda_{j}^{(A)}.$

We note that this relation holds also for $i=j$. Combining this equation with
(4.14) leads to the following theorem.

######  Theorem 4.11.

Under Assumption 4.10 the expected value of triples can be given as

(4.15) $[A_{i}AA_{j}]=\frac{[A_{i}A][A_{j}A]}{[A]},$

which can also be found in [6].

Focus for the moment on Poisson closures in SIS epidemics.

######  Corollary 4.12.

For SIS epidemics, closure (4.15) yields

(4.16)
$[SSI]\approx\frac{[SS][SI]}{[S]},\qquad[ISI]\approx\frac{[SI]^{2}}{[S]}.$

The Poisson model is not as ad-hoc as it may initially seem, but it actually
arises from the multinomial model in the case of Erdős-Reńyi graphs with low
average degree. For the remainder of this section we restrict the discussion
to Erdős-Rényi graphs of $N+1$ nodes, but with low average degree $\langle
k\rangle=D$.

In this case, we assume that the probability of a link being present in the
network is $\frac{D}{N}$. The degree of each node $v$ is a random Binomial
number $\deg v\sim\text{Bin}(N,\frac{D}{N})$. The law of rare events suggests
that for $N$ large, the degree distribution is approximately Poisson$(D)$, so
for this particular discussion we assume the network has marginal degree
distribution

$\deg v\sim\text{Poisson}(D).$

For clarity of the exposition, focus only on $SIS$ epidemics on the network.
We still operate under the multinomial link distribution (Assumption 4.1),
that is,

$(N_{S}^{(v)},N_{I}^{(v)})|v=S\sim\text{Mult}(\deg
v,2;p_{S}^{(S)},p_{I}^{(S)}).$

However, here the number of trials of the multinomial is a random Poisson
number $\deg v$. Given that the degree is $x\in\mathbb{N}$, then the
distribution of types is multinomial. We first use the law of total
probability to find the marginal distribution of $N_{I}^{(v)}$:

$\displaystyle\mathbb{P}\\{N_{I}^{(v)}=k\\}$
$\displaystyle=\sum_{x=0}^{\infty}P\\{N_{I}^{(v)}=k|\deg
v=x\\}\mathbb{P}\\{\deg v=x\\}$ $\displaystyle=\sum_{x=k}^{\infty}{x\choose
k}(p_{I}^{(S)})^{k}(1-p_{I}^{(S)})^{x-k}e^{-D}\frac{D^{x}}{x!}$
$\displaystyle=e^{-D}\frac{D^{k}(p_{I}^{(S)})^{k}}{k!}\sum_{x=k}^{\infty}\frac{(D(1-p_{I}^{(S)}))^{x-k}}{(x-k)!}=e^{-p_{I}^{(S)}D}\frac{D^{k}(p_{I}^{(S)})^{k}}{k!},$

and therefore the conditional marginal
$[N_{I}^{(v)}|v=S]\sim\text{Poisson}(Dp_{I}^{(S)})$. Similarly, the marginal
distribution of $N_{S}^{(v)}|v=S$ is that of a Poisson$(Dp_{S}^{(S)}))$ and
conveniently, we may write

$(N_{S}^{(v)}|v=S)=(\deg v-N_{I}^{(v)}|v=S)=\deg v-\text{Bin}(\deg
v,p_{I}^{(S)})\sim\text{Poisson$(Dp_{S}^{(S)}))$ }.$

Moreover, these two marginals are in fact independent. Here is the short
calculation. Since $\deg v$ can take values in $\mathbb{N}$, as is Poisson
distributed, the same holds for both $N_{S}^{(v)}$ and $N_{I}^{(v)}$. Then

$\displaystyle\mathbb{P}\\{N_{S}^{(v)}=\ell,N_{I}^{(v)}=m\\}$
$\displaystyle=\mathbb{P}\\{\deg
v-N_{I}^{(v)}=\ell,N_{I}^{(v)}=m\\}=\mathbb{P}\\{\deg
v=m+\ell,N_{I}^{(v)}=m\\}$ $\displaystyle=P\\{N_{I}^{(v)}=m|\deg
v=m+\ell\\}\mathbb{P}\\{\deg v=m+\ell\\}$ $\displaystyle={{\ell+m}\choose
m}(p_{I}^{(S)})^{m}(1-p_{I}^{(S)})^{\ell}e^{-D}\frac{D^{\ell+m}}{(\ell+m)!}$
$\displaystyle=e^{-p_{I}^{(S)}D}\frac{(p_{I}^{(S)}D)^{m}}{m!}\cdot
e^{-(1-p_{I}^{(S)})D}\frac{(D(1-p_{I}^{(S)}))^{\ell}}{\ell!}$
$\displaystyle=\mathbb{P}\\{N_{I}^{(v)}=m\\}\mathbb{P}\\{N_{S}^{(v)}=\ell\\}.$

To summarise, on Erdős-Rényi graphs with low average degree but many nodes,
the fact that the marginal degree is Poisson distributed reduces the
multinomial link model to the Poisson model. For the Erdős-Rényi case, the
assumption that the marginal distributions of types are Poisson, with rates
given by $\lambda_{S}^{(S)}=\langle k\rangle p_{S}^{(S)}$ and
$\lambda_{I}^{(S)}=\langle k\rangle p_{I}^{(S)}$ respectively, and the fact
that they are independent, is proven.

Because of this, we expect that when the Poisson closure is used for the mean-
field epidemic model on Erdős-Rényi graphs, the approximation performs at
least as well as the one using the moments closure given by Corollary 4.7.
This is also suggested by the simulations in the middle row of Figure 5, where
the two curves lie on top of each other, and both approximate the semi-stable
regime better than the binomial closure, when the infection rate is not too
small.

##  5\. Comparison of the closed pairwise model to stochastic simulations and
further results

We performed a fair number of numerical test where we compared the accuracy of
four closures (i.e. binomial, Poisson, moments and $p_{I}$ given by equations
(4.9), (4.16), (4.7)-(4.8) and (4.11), respectively) by comparing the expected
proportion of infected nodes in time from the closed system with that
resulting from explicit stochastic simulation of the SIS dynamics on regular,
Erdős-Rényi and bi-modal networks. It is well known that agreement between
closed pairwise systems and simulation tend to improve for large epidemics. In
this case this is equivalent to either having dense networks or high rates of
transmission, when assuming that the recovery rate is fixed. However, from the
viewpoint of deriving closures, it is well-known that the most interesting
behaviour may be for networks with low average degree. Due to this our tests
are conducted for networks with an average degree equal to six. The recovery
rate is also fixed and the value of the transmission rate $\tau$ is varied to
go from small to large epidemics. The results are shown in Figure 5.

Several observations can be made. As expected, for large epidemics the various
closures lead to very similar output and the agreement with simulation is
good. However, the Poisson closure leads to an over estimation in all cases,
while the binomial closure does better especially for regular and bi-modal
networks. For regular networks, the moments and the binomial closure coincide.
In general, the binomial closure does better on regular networks. The new
moments closure in many cases leads to results which are close to those based
on the Poisson closure. One such case is when the graph is Erdős-Rényi with
low average degree where the multinomial link model becomes the Poisson model,
as we discussed in the previous section.

In arbitrary networks of low average degree the independence assumption for
the Poisson marginals is unrealistic. This is because the rates add up to the
average degree but the probability of the sum of two independent Poisson to be
equal to the degree of a node, for all nodes, with small fluctuations is
small. In other words, it is expected that the binomial or moments
approximations would perform better as the total degree of a node is built
into the modelling. This can be seen from the closure’s performance in Figure
5 on the regular and bi-modal degree networks when the average degree is low
(first and last row of Figure 5).

When the infection rate $\tau$ of the SIS epidemics is low, the probability
$p_{I}^{(S)}$ is low as well. Then the mass function of the binomial
distribution with probability of success $p_{I}^{(S)}$ and $\langle k\rangle$
trials and the mass function of the Poisson distribution with rate $\langle
k\rangle p_{I}^{(S)}$ are very close. This implies that the information needed
for judging which closures are suitable for the epidemic model is not
contained only in the marginal distributions of types around nodes.

In Figure 5 the first column shows the estimated mean field epidemic model
using all closures for low $\tau$, and the moment closure (cyan line) does
estimate better than, or as good as, the Poisson one, while the empirical
distribution of infected nodes is well-approximated by both Binomial and
Poisson. However, when the value of $\tau$ becomes high, the marginal
distributions do contain more information. For example in Figure 3, we see
that the binomial mass function (black line) approximates the empirical
distribution better than the Poisson (dashed magenta line) and this is
reflected in the better approximation that can be seen in the first row of
Figure 5.

Figure 5. The average of 200 individual realisations (thick grey line)
together with output from four mean-field models (binomial (black), Poisson
(dashed magenta), moments (cyan) and $p_{I}$(yellow) given by equations (4.9),
(4.16), (4.7)-(4.8) and (4.11), respectively). The regular, Erdős-Rényi and
bi-modal networks are shown on the top, middle and bottom row, respectively.
All networks have an average degree of $\langle k\rangle=6$ with the bi-modal
having half of the nodes with degree four and the the other half with degree
eight. In all figures $N=200$ and $\gamma=1$. The rate of infection in the
left, middle and right columns is $\tau=0.25$, $\tau=0.5$ and $\tau=0.75$,
respectively. Each epidemic starts with five infected nodes chosen uniformly
at random and only epidemics that reach ten infected nodes count towards the
average with time reset to $t=0$ when the state with ten infected nodes is
first reached.

Applying the slightly more general closures (4.7)-(4.8) in the unclosed
pairwise system (2.1)-(2.5) leads to following equation for the expected
number of $[SI]$ pairs,

(5.1)
$\dot{[SI]}=\tau\mathcal{D}\frac{[SS][SI]}{[S]}-\tau\mathcal{D}\frac{[SI][SI]}{[S]}-\tau[SI]-\gamma[SI]+\gamma[II],$

where

(5.2) $\mathcal{D}=\frac{\langle\deg v(\deg v-1)\rangle}{\langle\deg
v\rangle^{2}}.$

Since the epidemic is driven by the number and growth rate of $[SI]$, we can
analyse the equation above at the disease free steady state, i.e.
$([S],[I],[SS],[SI],[II])=(N,0,\langle\deg v\rangle N,0,0)$. This allows us to
work out the rate of growth by looking at all the positive terms in (5.1) when
evaluated at time $t=0$. This leads to

(5.3) $\dot{[SI]}\sim\tau\mathcal{D}\frac{\langle\deg v\rangle
N}{N}[SI]=\frac{\tau\langle\deg v(\deg v-1)\rangle}{\langle\deg
v\rangle}[SI].$

However, the lifetime of an S-I link is simply $1/(\tau+\gamma)$ and hence,
the average number of S-I links produced by a typical S-I link during its
lifetime is

(5.4) $\mathcal{R}=\frac{\tau}{\tau+\gamma}\frac{\langle\deg v(\deg
v-1)\rangle}{\langle\deg v\rangle},$

which is a well know quantity in the epidemics on networks mathematical theory
[2]. This really reassuring as it suggests that our newly derived closure is
based on sound assumptions.

##  6\. Further extensions of the method

###  6.1. Conditioning on a link

There are many ways to count triples in a network. In the multinomial link
example we counted triples $[ABC]$ by focusing on the state of the middle
node. A different way would be to count them by looking at the neighbours of
node $v_{2}$, given that the link $(v_{2},v_{3})=(BC)$. In order to count this
way we will use equation (3.1) by conditioning on the values of $v_{2}$ and
$v_{3}$.

Let $\mathcal{G}$ be an $n$-regular graph. Assume we want to compute closures
of the form $[A_{i}AB]$, so we are given that the state of $v_{2}=A$ and state
of $v_{3}=B$, where $A,B\in\\{A_{1},\ldots,A_{m}\\}$. Homogeneity of the graph
implies that for any triple $(v_{1},v_{2},v_{3})\in\Pi_{3}$ the conditional
probabilities $P\\{v_{1}=A_{i}|(v_{2},v_{3})=(A,B)\\}$ do not depend of the
choice of particular triple $(v_{1},v_{2},v_{3})$ and they are only a function
of $A_{i},A$ and $B$. Denote their common value

$P\\{v_{1}=A_{i}|(v_{2},v_{3})=(A,B)\\}=p_{A_{i}|A-B}.$

Now, we compute from equation (3.1) to obtain

$\displaystyle[A_{i}AB]$
$\displaystyle=\sum_{(v_{1},v_{2},v_{3})\in\Pi_{3}}P\\{v_{1}=A_{i}|(v_{2},v_{3})=(A,B)\\}\mathbb{P}\\{(v_{2},v_{3})=(A,B)\\}$
$\displaystyle=\sum_{(v_{1},v_{2},v_{3})\in\Pi_{3}}p_{A_{i}|A-B}\mathbb{P}\\{(v_{2},v_{3})=(A,B)\\}$
$\displaystyle=\sum_{(v_{2},v_{3})=(A,B)}\mathbb{P}\\{(v_{2},v_{3})=(A,B)\\}\sum_{v_{1}:v_{1}\sim
v_{2},v_{1}\neq v_{3}}p_{A_{i}|A-B}$
$\displaystyle=\sum_{(v_{2},v_{3})=(A,B)}\mathbb{E}(\mathbbm{1}\\{(v_{2},v_{3})=(A,B)\\})(n-1)p_{A_{i}|A-B}$
$\displaystyle=[AB](n-1)p_{A_{i}|A-B}.$

Hence, in order to compute the closure $[A_{i}AB]$, one needs a good
approximation for the probabilities $p_{A_{i}|A-B}$.

Note that, again, under Assumption 4.1 of the multinomial link model of
Section 4.1, one can immediately see that

$p_{A_{i}|A-B}=p_{i}^{(v_{2},A)}=p_{i}^{(A)}=\frac{[A_{i}A]}{n[A]},$

which leads to the same result as given in Theorem 4.6 and in [8, 1]. The
multinomial link model does not take into account the extra information given
by $B$, so we expect that more refined closures can be discovered by taking
this information into account.

###  6.2. Closures for the Uniform distribution of states of the neighbours

We reiterate the arguments of Section 4.1 in an example where the conditional
distribution of types in neighbouring nodes is uniform. For simplicity we are
assuming that the network is $n$-regular.

######  Example 2.

(Uniform link distribution around a node.) Consider $m$ possible states for a
node on a network $\mathcal{G}$ that is $n$-regular. Furthermore, we know that
given the state $A$ of a node $v$ there exists an integer number
$k_{i,j}^{(A)}$, $0\leq k_{i,j}^{(A)}\leq n$ the conditional distribution

(6.1)
$(N_{i}^{(v)},N_{j}^{(v)})|v=A\sim\text{Uniform}[\Delta_{k_{i,j}^{(A)}}],$

where
$\Delta_{\ell}=\\{(x_{1},x_{2}):x_{i}\in\mathbb{Z}_{+},x_{1}+x_{2}\leq\ell\\}$.
Using this, and starting from equation (3.4), we calculate

$\displaystyle[A_{i}AA_{j}]$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E\Big{[}N_{i}^{(v)}N_{j}^{(v)}\Big{|}v=A\Big{]}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}\sum_{x=0}^{k_{i,j}^{(A)}}\sum_{y=0}^{k_{i,j}^{(A)}-x}xyP\\{N_{i}^{(v)}=x,N_{j}^{(v)}=y\\}$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}\sum_{x=0}^{k_{i,j}^{(A)}}\sum_{y=0}^{k_{i,j}^{(A)}-x}xy\frac{2}{(k_{i,j}^{(A)}+1)(k_{i,j}^{(A)}+2)}.$

The last line comes from the fact the conditional joint distribution is
uniform on the simplex $\Delta_{k_{i,j}^{(A)}}$. Now, the two inner sums can
be directly computed to be

$\displaystyle\sum_{x=0}^{k_{i,j}^{(A)}}\sum_{y=0}^{k_{i,j}^{(A)}-x}xy$
$\displaystyle=\frac{1}{2}\sum_{x=1}^{k_{i,j}^{(A)}}\big{\\{}x((k_{i,j}^{(A)})^{2}+k_{i,j}^{(A)})-x^{2}(2k_{i,j}^{(A)}+1)+x^{3}\big{\\}}$
$\displaystyle=\frac{1}{24}k_{i,j}^{(A)}(k_{i,j}^{(A)}+1)(k_{i,j}^{(A)}-1)(k_{i,j}^{(A)}+2).$

Substitute in the calculation for $[A_{i}AA_{j}]$ to obtain

(6.2) $\displaystyle[A_{i}AA_{j}]$
$\displaystyle=\frac{1}{12}\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}k_{i,j}^{(A)}(k_{i,j}^{(A)}-1)\approx\frac{[A]}{12}k_{i,j}^{(A)}(k_{i,j}^{(A)}-1).$

The last approximation holds when $\mathbb{P}\\{v=A\\}\approx[A]N^{-1}$ and
when the number $k_{i,j}^{(A)}$ does not depend on $v$. It may be necessary to
approximate $k_{i,j}^{(A)}$ using the network; one such approximation can be
obtained by taking expected values

$\displaystyle[AA_{i}]$
$\displaystyle=\sum_{v\in\mathcal{G}}\mathbb{P}\\{v=A\\}E[N_{i}^{(v)}|v=A]\approx\sum_{v\in\mathcal{G}}\frac{[A]}{N}\frac{k_{i,j}^{(A)}}{2}=[A]\frac{k_{i,j}^{(A)}}{2}.$

This gives that $k_{i,j}^{(A)}\approx 2\frac{[AA_{j}]}{[A]}.$ A similar
calculation would give that $k_{i,j}^{(A)}=2\frac{[AA_{i}]}{[A]}$. We can
combine the two to obtain two different estimates for $k_{i,j}^{(A)}$, namely

$\hat{k}_{i,j}^{(A)}=\frac{2}{[A]}\sqrt{[AA_{i}][AA_{j}]}\quad\text{ or
}\quad\tilde{k}_{i,j}^{(A)}=\frac{1}{[A]}([AA_{i}]+[AA_{j}]).$

A final substitution of these estimates in (6.2) yields two different closures
for $[A_{i}AA_{j}]$. In fact, there are infinitely many closures implied; if
one introduces a parameter $\alpha\in[0,1]$ then any convex combination
$\alpha\hat{k}_{i,j}^{(A)}+(1-\alpha)\tilde{k}_{i,j}^{(A)}$ is a different
closure.

##  7\. Discussion

A mean-field approximation to Markovian epidemics is widely used in various
scientific disciplines. Its strength relies on a drastic reduction of the
number of equations, which are a priori exponentially (in the number of nodes)
many, but in the mean-field model are polynomially many. Usually, mean-field
models are made possible by using some kind of closure which are
approximations of higher-order moments by lower-order ones, e.g. approximating
the expected number of triples by the expected number of singles and pairs.
The earlier such approximations are performed the more information is lost.
For example, closing the pairs is simpler than approximating triples but the
accuracy of approximation will be worse. For the standard SIS and SIR
epidemics one typically closes (or approximates) triples [6, 3, 4].

In this article we presented a top-down probabilistic approach to obtain a
rigorous formula for the expected number of $k$-tuples (and in particular
triples) in multi-type epidemics (Prop. 3.1). No assumptions were necessary
for the calculations, so the proposition works in all types of networks. We
then proceeded by approximating the expected value of triples in various ways.
Each method of approximating gave rise to a closure - either in a theoretical
or a numerical/statistical way.

Our main contributions outside of the robust theoretical framework mentioned
above, are summarised below:

  1. (1)

Derivation of a “moments” closure that performs well in SIS epidemics on
regular (or networks with low degree heterogeneity), Erdős-Rényi and bi-modal
graphs.

  2. (2)

Derivation of closures already used in the literature. We provide a careful
list of background assumptions that are necessary for each approximation to
work, and which are sufficient to derive these classical closures.

  3. (3)

Derivation of a few new closures which either bare similarities with existing
ones or allowed us to obtain epidemic threshold results that are well known in
the literature.

  4. (4)

A derivation of the Poisson link model from the multinomial link model, that
naturally arises for Erdős-Rényi graphs of low connection probability.

  5. (5)

Numerical verification of theoretical assumptions, particularly the marginal
distribution of types of nodes around susceptible ones.

The new closures show some promise but need some more testing and better
understanding when these agree with known ones, on which network do they work
best and in what way are they different from the existing ones. The new
closure based on assuming a uniform distribution is a good example of how our
method can be extended and applied beyond epidemic models. Finally, we hope
that these new insight may stimulate more research and may lead to some more
rigorous results for closures and closed systems in whatever context or
application.

## References

  * [1] R.C. Barnard, L. Berthouze, P.L. Simon and I.Z. Kiss, Journal of mathematical biology, online, 1–38, 2019.
  * [2] O. Diekmann and J.A.P. Heesterbeek, Mathematical epidemiology of infectious diseases: model building, analysis and interpretation, John Wiley & Sons, 2000.
  * [3] M.J. Keeling, Proc. R. Soc. Lond. Ser. B: Biol. Sci. 266, 859–867, 1999.
  * [4] I.Z. Kiss, J.C. Miller and P.L. Simon, Mathematics of epidemics on networks: from exact to approximate models, Springer, 2017.
  * [5] K. Christian, Moment closure a brief review, in Control of self-organising nonlinear systems, 253–271, Springer, 2016.
  * [6] D.A. Rand, Correlation Equations and Pair Approximations for Spatial Ecologies, in Advanced ecological theory: principles and applications. 100– 142\. Blackwell Science, Oxford, 1999.
  * [7] P.L. Simon and I.Z. Kiss, J. of Compl. Networks 4, 187–200, 2015.
  * [8] M. Taylor, P.L. Simon, P. L., D.M. Green, T. House and I.Z. Kiss, J. Math. Biol. 64, 1021–1042 (2012).