# Accuracy criterion for mean field approximations of Markov processes on
hypergraphs 111Supported by the ÚNKP-21-1 New National Excellence Program of
the Ministry for Innovation and Technology from the source of the National
Research, Development and Innovation Fund.  222 Partially supported by the ERC
Synergy under Grant No. 810115 - DYNASNET.

Dániel Keliger
Department of Stochastics, Institute of Mathematics,
Budapest University of Technology and Economics
Műegyetem rkp. 3., H-1111 Budapest, Hungary;
Alfréd Rényi Institute of Mathematics, Budapest, Hungary
e-mail<EMAIL_ADDRESS>
Illés Horváth
MTA-BME Information Systems Research Group
e-mail<EMAIL_ADDRESS>

###### Abstract

We provide error bounds for the $N$-intertwined mean-field approximation
(NIMFA) for local density-dependent Markov population processes with a well-
distributed underlying network structure showing NIMFA being accurate when a
typical vertex has many neighbors. The result justifies some of the most
common approximations used in epidemiology, statistical physics and opinion
dynamics literature under certain conditions. We allow interactions between
more than 2 individuals, and an underlying hypergraph structure accordingly.

##  1 Introduction

The analysis of stochastic population processes is an important topic in
several disciplines, such as epidemiology, biology, economics or computer
systems [5, 2, 12, 6, 27]. Such processes consist of a large number of
interacting individuals (agents) that execute random actions based on the
behavior of other individuals.

A widely-used framework is Markov population processes, where each individual
is in a local state from a fixed, finite state space, and can change their
state in a Markovian manner. For such models, the state space increases
exponentially with the population size, making an exact analysis infeasible
even for moderate population sizes, instead raising the question of good
approximations as the next best thing.

The classical result of Kurtz [16, 17] is based on two main assumptions: that
each individual can observe the entire population, and that the Markovian
transition rates of each individual depend on the observation in a density-
dependent manner. The conclusion is that, as the number of individuals
diverges, the evolution of the stochastic system converges to a deterministic
mean-field limit. This limit is straightforward to compute numerically, and
can serve as a good approximation of the stochastic system when the number of
individuals is large. The mean-field limit of Kurtz is referred to as the
_homogeneous mean-field approximation_ in the present paper.

While the density-dependent Markov setting is flexible and covers many
potential applications, the assumption that each individual can observe the
entire population is very restrictive. In many population processes arising
from real-life examples, individuals do not have full information about the
entire population; instead, each individual can observe only a subset of the
population. This information structure can be described by a network topology,
where each individual has interactions only with its neighbors according to
that topology.

The $N$-intertwined mean field approximation (NIMFA) [19] is a quenched mean-
field approximation, where differential equations are considered for each
individual based on their expected evolution. NIMFA is a deterministic process
different from the homogeneous mean-field approximation that incorporates the
network structure naturally, making it a potentially more accurate
approximation. On the flip side, the computational complexity is considerably
increased compared to the homogeneous mean-field approximation; nevertheless,
it remains tractable for population sizes large enough to make it relevant for
practical applications. Unfortunately, unlike for homogeneous systems, the
justification for using NIMFA is poorly understood, mostly relying on
numerical evidence [18, 28] along with a few theoretical results [29, 30, 31,
24].

In the present paper, we focus on a specific class of Markov processes dubbed
_local density-dependent Markov population processes_ , which preserves the
density-dependent assumption of Kurtz, but allows an underlying network
structure that dictates the environments observed by each individual. This
setting covers many of the frequently used stochastic models, such as the SIS
process in epidemiology [7, 13, 14, 3], Glauber dynamics in statistical
physics [10, 22], or the voter model and majority vote in opinion dynamics
[21, 23]. We incorporate interactions between more than 2 vertices into the
model with an underlying hypergraph structure accordingly to reflect on some
recent developments in the theory of higher order interactions.

We provide general error bounds for NIMFA that are strong on well-distributed
networks. Furthermore, under additional homogeneity assumptions, such as
annealed or activity driven networks [11, 26] we show these error bounds to be
small, with the added benefit of further reducing the number of equations to
other well-known approximations, like the _heterogenous mean field
approximation_ [25]. Finally, we elaborate the on the argument given by K.
Devriendt and P. Van Mieghem [9] and show that Szemerédi’s regularity lemma
[32] can be applied to reduce the number of equations (depending on a given
$\varepsilon$ error).

The rest of the paper is structured as follows. Section 2 introduces basic
notation and setup for density-dependent Markov population processes along
with examples of models used in the literature to illustrate these concepts
and their applicability. Section 3 states the main results and also relates
them to the recent work of Sridhar and Kar [30, 31] and Parasnis et al. [24].
Section 4 discusses further reductions of NIMFA to more simple approximations
used throughout the literature. Section 5 contains a summary of this paper
along with the limitations of these results and possible directions for
further research.

Finally, proofs are contained in Section 6.

##  2 Setup

###  2.1 The underlying hypergraph

Let $G$ be a finite hypergraph on $N$ vertices. The vertex set is labeled
$[N]=\\{1,\dots,N\\}$. The hypergraph is not necessarily uniform; edges may
contain up to $M+1$ vertices. The edges are ordered, with the first vertex
being special, and we will usually use the notation $(i,j_{1},\dots,j_{m})$
for an edge where $1\leq m\leq M$ and $i,j_{1},\dots,j_{m}\in[N]$. The idea
behind the distinction of the first vertex in an edge is that
$w^{(m)}_{i,j_{1},\dots,j_{m}}$ will describe the strength of connections
where $j_{1},\dots,j_{m}$ have a joint effect on vertex $i$ (see Figure 1).

Figure 1: Edge (hyperedge) with weight $w^{(m)}_{i,j_{1},\dots,j_{m}}$.

The $M=1$ case corresponds to (directed) graphs.

We allow so-called _secondary loops_ (abbreviated as s. loop), which are
$(i,\underline{j})$ edges with non-distinct vertices among
$j_{1},\dots,j_{m}\in[N]$. Note that traditional loops for the $m=1$ case are
excluded from this definition.

We use the notation $[N]^{m}$ to denote the set of $m$-tuples, and
$\underline{j}$ abbreviates $(j_{1},\dots,j_{m})$.

For unweighted hypergraphs, adjacency indicators
$a^{(m)}_{i,j_{1},\dots,j_{m}}$ (where $1\leq m\leq M$ and
$i,j_{1},\dots,j_{m}\in[N]$)

$\displaystyle a_{i,\underline{j}}^{(m)}=\begin{cases}1\ \text{ if
$i,j_{1},\dots,j_{m}$ are on the same hyperedge}\\\ 0\ \text{else}\end{cases}$

describe the connections between the vertices.

In-degrees for $1\leq m\leq M$ are defined as

$\displaystyle d^{(m)}(i):=$
$\displaystyle\frac{1}{m!}\sum_{\underline{j}\in[N]^{m}}a_{i,\underline{j}}^{(m)},$
(1)

(where $m!$ is included to cancel the re-orderings of $\underline{j}$), and
the average in-degree for each $1\leq m\leq M$ is

$\displaystyle\bar{d}^{(m)}:=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}d^{(m)}(i).$

In the literature, some normalization is usually assumed. In the present
paper, we introduce normalized weights $w^{(m)}_{i,j_{1},\dots,j_{m}}$ and
corresponding normalized in-degrees

$\delta^{(m)}(i):=\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}.$

representing the total weight of $m$-interactions effecting vertex $i\in[N]$.
In the $M=1$ case (classical graphs) we tend to omit the upper index $(m)$ and
write $w_{i,\underline{j}}^{(m)}$ simply as $w_{ij},$ and we also utilize the
matrix notation $W=(w_{ij})_{i,j\in[N]}.$

We have two Conventions for the normalization.

Convention 1:

$\displaystyle
w_{i,\underline{j}}^{(m)}=\frac{a_{i,\underline{j}}^{(m)}}{m!\bar{d}^{(m)}},\qquad\delta^{(m)}(i)=\frac{d^{(m)}(i)}{\bar{d}^{(m)}}.$
(2)

Convention 2:

$\displaystyle
w_{i,\underline{j}}^{(m)}=\frac{a_{i,\underline{j}}^{(m)}}{m!d^{(m)}(i)},\qquad\delta^{(m)}(i)=1.$
(3)

(The same $m!$ from (1) is now included in the conventions.)

For either convention, whenever the denominator would be 0, the numerator will
also be 0, and $w_{i,\underline{j}}^{(m)}$ is simply set to 0 as well.

We set

$\displaystyle w_{\max}=\max_{i,\underline{j},m}w_{i,\underline{j}}^{(m)}.$

We are going to set regularity assumptions for the weights and degrees:

$\displaystyle\delta^{(m)}(i)\leq$ $\displaystyle\,\delta_{\max},$ (4)
$\displaystyle\sum_{\begin{subarray}{c}\underline{j}\in[N]^{m}\\\
\underline{j}\textrm{ s. loop}\end{subarray}}w_{i,\underline{j}}^{(m)}\leq$
$\displaystyle\,R\sqrt{w_{\max}}.$ (5)

For Convention 2, (4) always holds. For Convention 1, we need
$d^{(m)}(i)\leq\delta_{\max}\bar{d}^{(m)}$ (upper regularity of the
hypergraph).

(5) always holds for $M=1$. It also obviously holds if there are no secondary
loops. In other cases, it is an actual restriction on the total weight of
secondary loops.

Symmetry is in general not assumed, that is, the hypergraph may be directed.

For some results concerning classical graphs ($M=1$) with Convention 2, the
extra assumption is needed for _out-_ degrees as well.

$\displaystyle\delta^{\textrm{out}}(j):=\sum_{i\in[N]}w_{ij}\leq\delta_{\textrm{max}}^{\textrm{out}}$
(6)

Assumption (4) and (6) can be understood as a weaker version of double
stochasticity of $W$ assumed in [30, 31].

###  2.2 Local density dependent Markov population process

We define a Markov process on the hypergraph. Each vertex is in a state from a
finite state space $\mathcal{S}$. $\xi_{i,s}(t)$ denotes the indicator that
vertex $i$ is in state $s$ at time $t$; the corresponding vector notation is

$\xi_{i}(t)=\left(\xi_{i,s}(t)\right)_{s\in\mathcal{S}}.$

We also introduce the notation

$\xi^{(m)}_{\underline{i},\underline{s}}(t)=\prod_{k=1}^{m}\xi_{i_{k},s_{k}}(t),$

where $\underline{i}=(i_{1},\dots,i_{m})$ is an edge and
$\underline{s}=(s_{1},\dots,s_{m})$ is a collection of states
($s_{k}\in\mathcal{S},k=1,\dots,m)$.
$\xi^{(m)}_{\underline{i},\underline{s}}(t)$ describes the indicator of
vertices $i_{1},\dots,i_{m}$ being in states $s_{1},\dots,s_{n}$ at time $t$,
respectively.

We define the _$m$ -neighborhood_ of vertex $i$ corresponding to
$\underline{s}=(s_{1},\dots,s_{m})$ as

$\displaystyle\phi^{(m)}_{i,\underline{s}}(t)=\sum_{\underline{j}\in\\{N\\}^{m}}w^{(m)}_{i,\underline{j}}\xi^{(m)}_{\underline{j},\underline{s}}(t).$
(7)

Some explanation is in order. Let $\underline{s}=(s_{1},\dots,s_{m})$ be fixed
for now. According to (7), we consider all edges that include $i$ and $m$
other vertices; for each such edge, we check whether the $m$ other vertices
are exactly according to the configuration of states described by
$\underline{s}$; if yes, their contribution to
$\phi^{(m)}_{i,\underline{s}}(t)$ is $w^{(m)}_{i,\underline{j}}$, otherwise
their contribution is 0.

The $m$-neighborhoods of $i$ consist of $\phi^{(m)}_{i,\underline{s}}(t)$ for
all possible configurations of states $\underline{s}$. The corresponding
vector notation is

$\displaystyle\phi^{(m)}_{i}(t)=\left(\phi^{(m)}_{i,\underline{s}}(t)\right)_{\underline{s}\in\mathcal{S}^{m}},$
(8)

and we may even write

$\displaystyle\phi_{i}(t)=\left(\phi^{(m)}_{i}(t)\right)_{m=1}^{M}$ (9)

for the entire neighborhood of $i$.

In (7), the normalized weights $w^{(m)}_{i,\underline{j}}$ are used; in case
$w^{(m)}_{i,\underline{j}}=0$ for some $\underline{j}$, the corresponding
interaction is simply not present.

Each vertex may transition to another state in continuous time. The transition
rates of a vertex may depend on all of its $m$-neighborhoods for $1\leq m\leq
M$; accordingly, the transition rate from $s^{\prime}$ to $s$ is described by
the function

$q_{ss^{\prime}}:\otimes_{m=1}^{M}\mathbb{R}^{\mathcal{S}^{m}}\to\mathbb{R}$

for each $s^{\prime}\neq s\in\mathcal{S}$.

We assume $q_{ss^{\prime}}$ is locally Lipschitz, and we also require
$q_{ss^{\prime}}(\phi^{(1)},\dots,\phi^{(M)})\geq 0$ for non-negative inputs.

For “diagonal” rates,

$q_{ss}:=-\sum_{s^{\prime}\neq s}q_{s^{\prime}s}$

corresponds to the total outgoing rate from state $s$.

The corresponding transition matrix is
$Q=\left(q_{ss^{\prime}}\right)_{s,s^{\prime}\in\mathcal{S}}.$ We emphasize
that in this convention $q_{ss^{\prime}}$ refers to an $s\leftarrow
s^{\prime}$ transition and not an $s\to s^{\prime}$ one. This ordering allows
us to use column vectors and matrix multiplication from the left.

The dynamics of $\left(\xi_{i}(t)\right)_{i=1}^{N}$ is a continuous-time
Markov chain with state-space $\mathcal{S}^{N}$ where each vertex performs
transitions according to the transition rates $q_{s^{\prime}s}$, independently
from the others. After a transition, vertices update their neighborhood
vectors $\phi_{i}(t)$. We call such dynamics _local-density dependent Markov
processes_.

We define the process $\left(\xi_{i,s}\right)_{i,s}$ formally via Poisson
representation:

$\displaystyle\begin{split}{\xi}_{i,s}(t)=&\xi_{i,s}(0)+\sum_{\begin{subarray}{c}s^{\prime}\in\mathcal{S}\\\
s^{\prime}\neq
s\end{subarray}}\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(t)\right)-\mathcal{N}_{i,s^{\prime}s}\left(\mathcal{H}_{i,s^{\prime}s}(t)\right),\\\
\mathcal{H}_{i,ss^{\prime}}(t)=&\left\\{(\tau,x)\in\mathbb{R}^{2}\left.\right|0\leq\tau\leq
t,\ 0\leq x\leq
q_{ss^{\prime}}\left(\phi_{i}(\tau)\right){\xi}_{i,s^{\prime}}(\tau)\right\\},\end{split}$
(10)

where for each choice of $1\leq i\leq N$ and $s\neq s^{\prime}\in\mathcal{S}$,
$\left(\mathcal{N}_{i,ss^{\prime}}(x,y):x,y\geq 0\right)$ is a 2-dimensional
Poisson-process with density 1, and the processes are independent for
different $(i,s,s^{\prime})$ triples.

(10) is a cumulative formula counting all transitions of the vertex $i$ to and
from state $s$ up to time $t$; $s\leftarrow s^{\prime}$ transitions are
generated using the Poisson points in the 2-dimensional domain
$\mathcal{H}_{i,ss^{\prime}}(t)$ which has area
$\int_{0}^{t}q_{ss^{\prime}}\left(\phi_{i}(\tau)\right){\xi}_{i,s^{\prime}}(\tau)\mathrm{d}\tau$,
ensuring the proper transition rate for $s\leftarrow s^{\prime}$ jumps at time
$\tau$. The second term of the sum corresponds to $s^{\prime}\leftarrow s$
transitions in a similar manner.

###  2.3 N-intertwined mean field approximation

Although the state occupation probabilities of the population process can be
described by the Chapman–Kolmogorov equations, the number of equations is
$\left|\mathcal{S}\right|^{N}$, making it infeasible for numeric or analytic
investigations even for moderate sized populations. To address this issue,
several approximation schemes had been introduced in the literature with
varying complexity.

This chapter discusses the quenched mean field approximation [19], also called
the N-intertwined mean field approximation (NIMFA). NIMFA preserves all
information regarding the graph structure and only neglects dynamical
correlation between vertices. The goal is to derive state occupation
probabilities for each vertex separately, resulting in a total of
$\left|\mathcal{S}\right|N$ equations.

A possible intuition for NIMFA is as follows.

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\left(\xi_{i}(t)\right)=\mathbb{E}\left[Q\left(\phi_{i}(t)\right)\xi_{i}(t)\right]$
(11)

can be derived from (10). To close (11), we apply the approximation
$\phi_{i}(t)\approx\mathbb{E}\left(\phi_{i}(t)\right)$, which is reasonable
when $N$ is large and there is low correlation between vertices:

$\displaystyle\mathbb{E}\left[Q\left(\phi_{i}(t)\right)\xi_{i}(t)\right]\approx\mathbb{E}\left[Q\left(\mathbb{E}\left(\phi_{i}(t)\right)\right)\xi_{i}(t)\right]=Q\left(\mathbb{E}\left(\phi_{i}(t)\right)\right)\mathbb{E}\left(\xi_{i}(t)\right).$

Accordingly, the NIMFA approximation
$z_{i}(t)=(z_{i,s}(t))_{s\in\mathcal{S}},1\leq i\leq N$ is the solution of the
system

$\displaystyle\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}z_{i}(t)=&Q\left(\zeta_{i}(t)\right)z_{i}(t),\\\
\zeta_{i}(t)=&\left(\zeta_{i}^{(m)}(t)\right)_{m=1}^{M},\\\
\zeta_{i}^{(m)}(t)=&\left(\zeta_{i,\underline{s}}^{(m)}(t)\right)_{\underline{s}\in[N]^{m}}=\left(\sum_{\underline{j}\in\mathcal{S}^{m}}w_{i,\underline{j}}^{(m)}z_{\underline{j},\underline{s}}^{(m)}(t)\right)_{\underline{s}\in\mathcal{S}^{m}},\end{split}$
(12)

where $z_{i}(t)$ corresponds to $\xi_{i}(t)$ and $\zeta_{i}(t)$ corresponds to
$\phi_{i}(t)$, and then the approximation used is

$\mathbb{P}\left(\xi_{i,s}(t)=1\right)=\mathbb{E}\left(\xi_{i,s}(t)\right)\approx
z_{i,s}(t).$

The following theorem ensures the existence and uniqueness of the solution of
(12).

######  Theorem 1.

Let $\Delta^{\mathcal{S}}$ denote the set of probability vectors from
$\mathbb{R}^{\mathcal{S}}.$ For any initial condition
$z_{i}(0)\in\Delta^{\mathcal{S}}$ for all $i$ the ODE system (12) has a unique
global solution such that $z_{i}(t)\in\Delta^{S}$ for all $i$ and $t>0$ as
well.

###  2.4 Examples

In this section we give some examples for models covered by the formalism of
Section 2.2.

#### The simplicial SIS model

We will use the simplicial SIS model, also referred to as the contact process
as a running example.

In the $M=1$ case (graphs) the setup is the following: Each vertex can be in
one of two states: susceptible ($S$) and infected ($I$), hence the state space
is $\mathcal{S}=\\{S,I\\}.$ Infected vertices become susceptible at a constant
rate $\gamma\geq 0$ while susceptible vertices receive the illness with rate
proportional to number of its infected neighbhours.

The number of infected neighbhours of vertex $i\in[N]$ at time $t$ equals to

$\sum_{j=1}^{N}a_{ij}\xi_{j,I}(t)$

as $a_{ij}\xi_{j,I}(t)$ the indicator of vertex $j$ is connected to vertex $i$
and that it is infected at time $t$. After normalizing it with $\bar{d}$ or
$d(i)$ depending on our choice of convention 1 or 2 one gets

$\sum_{j=1}^{N}w_{ij}\xi_{j,I}(t)=\phi_{i,I}(t).$

Therefore, the transition rates takes the form $q_{SI}(\phi_{i}(t))=\gamma,\
q_{IS}(\phi_{i}(t))=\beta\phi_{i,I}(t)$ where $\beta\geq 0$ is a suitable
constant factor. In matrix form:

$\displaystyle Q(\phi_{i}(t))=\left[{\begin{array}[]{cc}-\gamma&\gamma\\\
\beta\phi_{i,I}(t)&-\beta\phi_{i,I}(t)\\\ \end{array}}\right]$

For the SIS process NIMFA takes the form:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{i,I}(t)=-\gamma
z_{i,I}(t)+\beta(1-z_{i,I}(t))\sum_{j=1}^{N}w_{ij}z_{j,I}(t).$

Here we used $z_{i,S}(t)=1-z_{i,I}(t)$ which is also the reason why it enough
to write the $I$ components only.

The extension of the SIS model to hypergraphs is called the simplicial SIS
model. The curing rate stays $\gamma$, however the infection dynamics is
modified. A susceptible vertex can be infected via any $(m+1)$-edge if all
other $m$ vertices are infected. The weighted sum of such edges $(m+1)$-edges
is

$\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}\xi_{\underline{j},(I,\dots,I)}^{(m)}(t)=\phi_{i,(I,\dots,I)}^{(m)}(t).$

The infection rates is sum of all the $1\leq m\leq M$ with appropriate
$\beta_{1},\dots,\beta_{M}\geq 0$ factors:

$q_{IS}(\phi_{i}(t))=\sum_{m=1}^{M}\beta_{m}\phi_{i,(I,\dots,I)}^{(m)}(t).$

For the simplicial SIS model NIMFA takes the form

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{i,I}(t)=-\gamma
z_{i,I}(t)+\left(1-z_{i,I}(t)\right)\sum_{m=1}^{M}\beta_{m}\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}z_{\underline{j},(I,\dots,I)}^{(m)}(t),$

#### Glauber dynamics

Glauber dynamics is a stochastic process whose stationary distribution
coincides with the distribution given by a spin system, such as the Ising
model [10].

There are two possible states: $\mathcal{S}=\\{+,-\\}.$ Instead of the
indicators

$\xi_{i,+}(t),\,\xi_{i,-}(t)$

it is customary use the sign variables

$\sigma_{i}(t):=\xi_{i,+}(t)-\xi_{i,-}(t)=2\xi_{i,+}(t)-1.$

In physical systems it is natural to assume $w_{ij}$ is symmetric and
$w_{ii}=0$.

The dynamics is the following:

  * •

At each time step, choose a vertex $i$ uniformly.

  * •

With probability $p_{i}(\sigma)=\frac{e^{\beta S_{i}(\sigma)}}{e^{\beta
S_{i}(\sigma)}+1}$, vertex $i$ switches to state + (else -), where

$S_{i}(\sigma)=\sum_{j=1}^{N}w_{ij}\sigma_{j}.$

Note that $S_{i}(\sigma)$ arises from the reduction of the energy

$\displaystyle H(\sigma):=-\frac{1}{2}\sum_{i<j}w_{ij}\sigma_{i}\sigma_{j}$

when vertex $i$ is flipped from $-$ to $+$. The stationary distribution is
then given by the Gibbs measure

$\displaystyle P(\sigma):=$ $\displaystyle\frac{1}{Z}e^{-\beta H(\sigma)},$
$\displaystyle Z:=$ $\displaystyle\sum_{\sigma}e^{-\beta H(\sigma)}.$

We modify the above dynamics. First, note that, in accordance with (7),

$\displaystyle
S_{i}(\sigma(t))=\sum_{j=1}^{N}w_{ij}\left(\xi_{j,+}(t)-\xi_{j,-}(t)\right)=\phi_{i,+}(t)-\phi_{i,-}(t).$

With a slight abuse of notation, we denote

$\displaystyle
S\left(\phi_{i}(t)\right):=\alpha\phi_{i,+}(t)-\gamma\phi_{i,-}(t),$

allowing the dynamics to have a preferred state.

Furthermore, we turn to the continuous time version instead with transition
rates given by

$\displaystyle q_{+-}(\phi)=$ $\displaystyle e^{\beta S(\phi)},$
$\displaystyle q_{-+}(\phi)=$ $\displaystyle 1.$

Since there are only two states, it is enough to consider the probabilities of
occupying state +. For this, NIMFA gives the following system of ODEs:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{i,+}(t)=$
$\displaystyle(1-z_{i,+}(t))e^{\beta S(\zeta_{i}(t))}-z_{i,+}(t).$ (13)

The equilibrium state is given by the fixed point problem

$\displaystyle z_{i,+}=\frac{e^{\beta S(\zeta_{i})}}{e^{\beta
S(\zeta_{i})}+1}.$ (14)

Assume $\alpha=1,\gamma=-1$ as in the original setting and that the underlying
weighted graph is regular: $\forall i\ \delta(i)=\sum_{j}w_{ij}=1$. Than (14)
reduces to

$\displaystyle\sigma=\tanh\left(\frac{1}{2}\beta\sigma\right),$
$\displaystyle\forall i\ 2z_{i,+}-1=\sigma$

giving back the classical mean field approximation of the Ising model on
lattice. This is not surprising as both NIMFA and the classical mean field
approach is based on the assumption of independence of vertices.

Based on [22], we can generalize the model for hypegraphs via extending
$S(\phi)$ to

$\displaystyle
S(\phi_{i}(t)):=\sum_{m=0}^{M}\alpha_{m}\phi_{i,(+,\dots,+)}^{(m)}(t)-\gamma_{m}\phi_{i,(-,\dots,-)}^{(m)}(t)$

allowing the system to lose even more energy when $3$ or more neighbors have
the same configuration on a hyper-edge.

#### The voter model

The voter model is a conceptually simple stochastic process modeling opinion
dynamics [21]. In the most simple case, there are two possible states:
$\mathcal{S}={0,1}.$

The dynamics can be described the following way: At each time step, we choose
a vertex uniformly. Said vertex chooses an neighbor also uniformly, and copies
its state. Similarly to the Glauber dynamics, we will study the continuous
time version instead.

For vertex $i$, the ratio of neighbors sharing belief $s\in\\{0,1\\}$ is

$\displaystyle\frac{1}{d(i)}\sum_{j=1}^{N}a_{ij}\xi_{j,s}(t)=\phi_{i,s}(t)$

with the choice of Convention 2. Hence, the transition rates take the form

$\displaystyle q_{01}(\phi_{i}(t))=$ $\displaystyle\lambda\phi_{i,0}(t),$
$\displaystyle q_{10}(\phi_{i}(t))=$
$\displaystyle\lambda\phi_{i,1}(t)=\lambda\left(1-\phi_{i,0}(t)\right).$

Using $z_{i,1}(t)=1-z_{i,0}(t)$, NIMFA can be written as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{i,0}(t)=-\lambda(1-\zeta_{i,0}(t))z_{i,0}(t)+\lambda\zeta_{i,0}(t)\left(1-z_{i,0}(t)\right).$

#### A modified majority rule model

Another popular model of opinion dynamics is the majority rule [21]. In this
setting a group of $m+1$ individuals are choosen who update their state
simultaneously to the majority opinion. Ties are usually broke with either a
random choice or setting a preferred opinion to win in this case, say opinion
$1$. For the sake of simplicity, we apply the latter approach.

Due to the continuous time setting we use, we modify the majority rule such
that only one individual updates its opinion during a transition based on the
state of the other vertices (not including its own opinion for the sake of
simplicity).

As it is stated in [21], the hypergraph setting is more suitable for majority
rule. We assume communities have a bounded size $M+1$, while each individual
can be a part of many, possibly overlapping communities.

$a_{i,j_{1},\dots,j_{m}}^{(m)}$ is the indicator of vertices
$i,j_{1},\dots,j_{m}\in[N]$ being in a community. We assume symmetry in the
indices and set $a_{i,j_{1},\dots,j_{m}}^{(m)}=0$ if there are duplicates. We
use a slightly modified version of Convention 1:

$\displaystyle
w_{i,\underline{j}}^{(m)}=\frac{\alpha_{m}a_{i,\underline{j}}}{m!\bar{d}^{(m)}},$

where $\alpha_{m}$ measures how much importance vertices put on communities of
size $m+1$. $w_{\textrm{max}}=\max_{m}\frac{\alpha_{m}}{m!\bar{d}^{(m)}}$ can
be small either due to vertices being part of many communities of size $m+1$
on average or because they put less importance on said communicates.

Introduce the notation $|s|=\sum_{l=1}^{m}s_{l}.$ Vertex $i$ in community
$i,j_{1},\dots,j_{m}$ changes its opinion to the majority of
$j_{1},\dots,j_{m}$ at rate $w_{i,\underline{j}}^{(m)}$. Therefore,

$\displaystyle q_{01}(\phi_{i}(t))=$
$\displaystyle\sum_{m=0}^{M}\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}{\mathds{1}}_{\left\\{0\
\textit{is the majority for $j_{1},\dots,j_{m}$}\right\\}}$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{M}\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}\sum_{|\underline{s}|<\frac{m}{2}}\prod_{l=1}^{m}\xi_{j_{l},s_{l}}(t)=\sum_{m=0}^{M}\sum_{|\underline{s}|<\frac{m}{2}}\phi_{i,\underline{s}}^{(m)}(t),$
$\displaystyle q_{10}(\phi_{i}(t))=$
$\displaystyle\sum_{m=0}^{M}\sum_{|\underline{s}|\geq\frac{m}{2}}\phi_{i,\underline{s}}^{(m)}(t).$

The NIMFA ODEs are

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{i,0}(t)=$
$\displaystyle(1-z_{i,0}(t))\sum_{m=0}^{M}\sum_{|\underline{s}|<\frac{m}{2}}\zeta_{i,\underline{s}}^{(m)}(t)-z_{i,0}(t)\sum_{m=0}^{M}\sum_{|\underline{s}|\geq\frac{m}{2}}\zeta_{i,\underline{s}}^{(m)}(t).$

##  3 Error bounds for NIMFA

In this section we are presenting our main results which bound the error
arising from neglecting the dynamical correlation between vertices.

Recall that (11) was closed by assuming
$\phi_{i}(t)\approx\mathbb{E}\left(\phi_{i}(t)\right)$. We introduce an
auxiliary process where the empirical neighborhood $\phi_{i}(t)$ is replaced
by the approximate $\zeta_{i}(t)$ from (12):

$\displaystyle\begin{split}\hat{\xi}_{i,s}(t)=&\xi_{i,s}(0)+\sum_{\begin{subarray}{c}s^{\prime}\in\mathcal{S}\\\
s^{\prime}\neq
s\end{subarray}}\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{K}_{i,ss^{\prime}}(t)\right)-\mathcal{N}_{i,s^{\prime}s}\left(\mathcal{K}_{i,s^{\prime}s}(t)\right),\\\
\mathcal{K}_{i,ss^{\prime}}(t)=&\left\\{(\tau,x)\in\mathbb{R}^{2}\left.\right|0\leq\tau\leq
t,\ 0\leq x\leq
q_{ss^{\prime}}\left(\zeta_{i}(\tau)\right)\hat{\xi}_{i,s^{\prime}}(\tau)\right\\}.\end{split}$
(15)

The process $\hat{\xi}_{i,s}(t)$ is an indicator process just like
$\xi_{i,s}(t)$, so it takes 0 or 1 values, and
$\sum_{s\in\mathcal{S}}\hat{\xi}_{i,s}(t)=1$ for any $i\in[N]$ and $t\geq 0$.
However, assuming independent initial conditions, $\hat{\xi}_{i}(t)$ remain
independent. Applying total expectation to (15) shows

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\left(\hat{\xi}_{i}(t)\right)=Q\left(\zeta_{i}(t)\right)\mathbb{E}\left(\hat{\xi}_{i}(t)\right),$

which, along with (12), implies that if
$\mathbb{E}\left(\hat{\xi}_{i}(0)\right)=z_{i}(0)$, then
$\hat{\xi}_{i}(t)-z_{i}(t)$ is a martingale and

$\displaystyle\mathbb{E}\left(\hat{\xi}_{i}(t)\right)=z_{i}(t)\quad\forall
t\geq 0.$ (16)

Using the same background Poisson processes $\mathcal{N}_{i,ss^{\prime}}$
provides a coupling between $\xi$ and $\hat{\xi}$ that will be useful later
on.

We aim to give an upper bound for $|\hat{\xi}(t)-\xi(t)|$, as well as for
$|\hat{\xi}(t)-z(t)|$. We start with $|\hat{\xi}(t)-\xi(t)|$ by introducing
the error terms

$\displaystyle D_{i}^{(0)}(t)=$ $\displaystyle\sup_{0\leq\tau\leq
t}\mathbb{E}\left(\sum_{s\in\mathcal{S}}\left|\xi_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right|\right),$
$\displaystyle\tilde{D}_{i}^{(0)}(t)=$
$\displaystyle\mathbb{E}\left(\sup_{0\leq\tau\leq
t}\sum_{s\in\mathcal{S}}\left|\xi_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right|\right).$

Apparently, the only difference between the two is the order in which we take
the supremum in time. $\tilde{D}_{i}^{(0)}(t)$ is more strict as

$\displaystyle D_{i}^{(0)}(t)\leq\tilde{D}_{i}^{(0)}(t).$

Observe that
$\sum_{s\in\mathcal{S}}\left|\xi_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right|$
only has two possible values: $0$ if $\xi_{i}(t)=\hat{\xi}_{i}(t)$, and $2$
otherwise (as there will be two $s\in\mathcal{S}$ indices where
$\xi_{i,s}(t),\hat{\xi}_{i,s}(t)$ differs). This implies

$\displaystyle\sup_{0\leq\tau\leq
t}\mathbb{P}\left(\xi_{i}(\tau)\neq\hat{\xi}_{i}(\tau)\right)=$
$\displaystyle\frac{1}{2}D_{i}^{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{(0)}}}(t),$
$\displaystyle\mathbb{P}\left(\exists\ 0\leq\tau\leq t:\
\xi_{i}(\tau)\neq\hat{\xi}_{i}(\tau)\right)=$
$\displaystyle\frac{1}{2}\bar{D}_{i}^{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{(0)}}}(t)$

We also introduce error terms describing the environments arising from
$\xi_{i}(t)$ and $\hat{\xi}_{i}(t)$:

$\displaystyle D_{i}^{(m)}(t)=$ $\displaystyle\sup_{0\leq\tau\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\phi_{i,\underline{s}}^{(m)}(\tau)-\zeta_{i,\underline{s}}^{(m)}(\tau)\right|\right]\quad(1\leq
m\leq M),$ $\displaystyle\tilde{D}_{i}^{(m)}(t)=$
$\displaystyle\,\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\phi_{i,\underline{s}}^{(m)}(\tau)-\zeta_{i,\underline{s}}^{(m)}(\tau)\right|\right]\quad(1\leq
m\leq M).$

Since the neighborhoods $\phi_{i}(t)$ and $\zeta_{i}(t)$ are constructed from
the indicators $\xi_{i}(t)$ and $\hat{\xi}_{i}(t)$, it is reasonable to expect
$\phi_{i}(t)$ and $\zeta_{i}(t)$ to be close to each other – as long as
$\xi_{i}(t)$ and $\hat{\xi}_{i}(t)$ are also close. To avoid circular
reasoning, we carry on handling these two types of errors together at the same
time. This motivates the introduction of

$\displaystyle D_{\max}^{(m)}(t)=$
$\displaystyle\max_{i\in[N]}D_{i}^{(m)}(t),$ $\displaystyle D_{\max}(t)=$
$\displaystyle\sum_{m=0}^{M}D_{\max}^{(m)}(t),$
$\displaystyle\tilde{D}_{i}(t)=$
$\displaystyle\sum_{m=0}^{M}\tilde{D}_{i}^{(m)}(t).$

The vector notation $\tilde{D}(t)=\left(\bar{D}_{i}(t)\right)_{i\in[N]}$ will
also be utilized.

Now we can go ahead to state the main results of the paper. The idea behind
the staments is when the vertex weights are generally small (the network is
well-distributed) then vertices has low correlation between each other, hence
NIMFA is accurate.

######  Theorem 2.

(Main)

Assume the initial conditions $\xi_{i}(0)$ are independent and (16) is
satisfied. Then for every $t\geq 0$ there is a constant
$C=C\left(t,\delta_{\max},R\right)$ such that

$\displaystyle\max_{i}\sup_{0\leq\tau\leq
t}\mathbb{P}\left(\xi_{i}(\tau)\neq\hat{\xi}_{i}(\tau)\right)\leq\frac{1}{2}D_{\max}(t)\leq$
$\displaystyle C\sqrt{w_{\max}}.$ (17)

Furthermore, if we additionally assume $M=1$ (having $1$-uniform hypergraphs)
then there exist constants
$C_{1}=C_{1}(\delta_{\max}),C_{2}=C_{2}(\delta_{\max})$ such that for all
$t\geq 0$

$\displaystyle\begin{split}\left\|\tilde{D}(t)\right\|\leq&C_{1}\exp\left(C_{2}\left\|W+I\right\|t\right)\|\mu\|,\\\
\mu=&\left(\sqrt{\sum_{j=1}^{N}w_{ij}^{2}}\right)_{i\in[N]},\end{split}$ (18)

where the norm $\|\cdot\|$ is arbitrary, $W=\left(w_{ij}\right)_{i,j=1}^{N}$
and $I$ is the identity matrix.

######  Remark 1.

The reason why we have different results for $M>1$ and $M=1$ is technical in
nature. The main observation is that in the $M=1$ case
$\hat{\xi}_{i,s}(t)-z_{i,s}(t)$ is a martingale making possible to take
$\sup_{0\leq\tau\leq t}$ inside the expectation via Doob’s inequality. It is
no longer the case for $M>1$ where
$\hat{\xi}_{\underline{i},\underline{s}}^{(m)}(t)-z_{\underline{i},\underline{s}}^{(m)}(t)$
is typically not a martingale itself.

(17) is a local result in the sense that it provides a uniform bound, ensuring
that $\hat{\xi}_{i,s}(t)$ and $\xi_{i,s}(t)$ are close for all vertices $i$
simultaneously. For example, in the SIS process it allows us to approximate
infection probabilities for concrete individuals, not just global or
mesoscopic ratios.

(18) will be elaborated on in Theorem 3.

In general, we cannot expect a similar local result for $\hat{\xi}_{i,s}(t)$
and $z_{i,s}(t)$ since $\hat{\xi}_{i,s}(t)$ is an indicator while $z_{i,s}(t)$
is a continuous variable. However, if we average out $\hat{\xi}_{i,s}(t)$ over
a macroscopic set of vertices, a similar result will hold.

In (18) the use of $\ell^{2}$ or $\ell^{\infty}$ is advised. Observe

$\displaystyle\|W\|_{\infty}=\max_{i}\sum_{j}w_{ij}\leq\delta_{\textrm{max}}$
$\displaystyle\|W\|_{2}\leq\sqrt{\|W\|_{1}\|W\|_{\infty}}=\sqrt{\left(\max_{j}\sum_{i}w_{ij}\right)\left(\max_{j}\sum_{i}w_{ij}\right)}\leq\sqrt{\delta_{\textrm{max}}^{\textrm{out}}\delta_{\textrm{max}}},$
(19)

making $\exp\left(C_{2}\left\|W+I\right\|t\right)$ bounded in (18). Note that
(19) is the only step where Assumption (6) regarding
$\delta_{\textrm{max}}^{\textrm{out}}$ is used.

As for $\|\mu\|$:

$\displaystyle\|\mu\|_{\infty}=$ $\displaystyle\max_{1\leq i\leq
N}\sqrt{\sum_{j=1}^{n}w_{ij}^{2}}\leq\max_{1\leq i\leq
N}\sqrt{w_{\textrm{max}}\sum_{j=1}^{n}w_{ij}}\leq\sqrt{w_{\textrm{max}}\delta_{\textrm{max}}},$
$\displaystyle\|\mu\|_{2}=$
$\displaystyle\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}^{2}}.$

Convention 1 works well with the $O\left(\sqrt{w_{\textrm{max}}}\right)$ error
bound as $w_{\textrm{max}}=\frac{1}{\bar{d}}$ holds in that case suggesting
vertices being close to independent when they have a lot of neighbors on
average. Similarly to (17), it also gives a uniform error bound, making it
possible to approximate the probabilities at the individual level. For
Convention 2 on the other hand, $w_{\textrm{max}}=\frac{1}{d_{\min}}$ is
sensitive to even one vertex with a low degree. If we are not attached to
uniform bounds in $i,$ we can provide a more robust on for the error of a
typical vertex, thus, it is possible to describe global or mesoscopic
population statistics.

Let $\iota\sim U\left([N]\right)$ the index of a randomly chosen vertex.

$\displaystyle\mathbb{P}\left(\exists\ \tau\in[0,t]:\
\xi_{\iota}(\tau)\neq\hat{\xi}_{\iota}(\tau)\right)=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\mathbb{P}\left(\exists\ \tau\in[0,t]:\
\xi_{i}(\tau)\neq\hat{\xi}_{i}(\tau)\right)\leq$
$\displaystyle\frac{1}{2N}\sum_{i=1}^{N}\tilde{D}_{i}(t)\leq\sqrt{\frac{1}{4N}\sum_{i=1}^{N}\tilde{D}_{i}^{2}(t)}=O\left(\sqrt{\frac{1}{N}\|\mu\|_{2}^{2}}\right)$

Observe

$\displaystyle\frac{1}{N}\|\mu\|_{2}^{2}=\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}^{2}$
(20)

is the squared and normalized Frobenius norm of the matrix $W.$ We mention
that such bound were used in [31] under more strict assumptions regarding $W$.

Note that for Convention 2

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}^{2}=\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\frac{a_{ij}}{d^{2}(i)}=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{d(i)},$
(21)

meaning the error is small when vertices typically have large degrees.

These observations along with Theorem 2 give the following result:

######  Theorem 3.

For $M=1$ (directed, weighted graphs), there exist constants
$C_{1}=C_{1}(t,\delta_{\textrm{max}})$ and
$C_{2}=C_{2}(t,\delta_{\textrm{max}},\delta_{\textrm{max}}^{\textrm{out}})$
such that

$\displaystyle\max_{i}\mathbb{P}\left(\exists\ \tau\in[0,t]:\
\xi_{i}(\tau)\neq\hat{\xi}_{i}(\tau)\right)\leq$ $\displaystyle
C_{1}\sqrt{w_{\textrm{max}}},$ (22)
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\mathbb{P}\left(\exists\ \tau\in[0,t]:\
\xi_{i}(\tau)\neq\hat{\xi}_{i}(\tau)\right)\leq$ $\displaystyle
C_{2}\sqrt{\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}^{2}}.$ (23)

So far, we have only accounted for the error between $\xi_{i}(t)$ and
$\hat{\xi}_{i}(t)$, however, what we are actually interested in is the
expectation $\mathbb{E}\left(\hat{\xi}_{i}(t)\right)=z_{i}(t)$, the solution
of the ODE system given by NIMFA. Thankfully,
$\left(\hat{\xi}_{i}(t)\right)_{i\in[N]}$ are independent, hence, their
averages must concentrate around the mean:

######  Theorem 4.

Assume (16) holds with independent initial conditions. Then for any $t\geq 0$
and any $1\leq K\leq N$,

$\displaystyle\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\sum_{s\in\mathcal{S}}\left|\frac{1}{K}\sum_{i=1}^{K}\left(\hat{\xi}_{i,s}(\tau)-z_{i,s}(\tau)\right)\right|\right]\leq\frac{2|\mathcal{S}|}{\sqrt{K}}.$
(24)

The most natural application of Theorem 4 is for $K=N$, but it is formulated
in a way so that it can be applied to any convenient subset of vertices (the
fact that the first $K$ vertices are considered has no significance as the
vertices can be reordered arbitrarily).

Together, Theorems 2, 3 and 4 give an error bound for the NIMFA approximation.

######  Theorem 5.

Assume (16) holds with independent initial conditions. Then for any $t\geq 0$,
there exists a constant $C=C(t,\delta_{\max},R)$ such that

$\displaystyle\sup_{0\leq\tau\leq
t}\mathbb{E}\left(\sum_{s\in\mathcal{S}}\left|\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i,s}(\tau)-z_{i,s}(\tau)\right)\right|\right)\leq
C\left(\sqrt{w_{\max}}+\frac{1}{\sqrt{N}}\right).$ (25)

Furthermore, if we additionally assume $M=1$, there exist constants
$C_{1}=C_{1}(t,\delta_{\max}),C_{2}=C_{2}(t,\delta_{\max},\delta_{\textrm{max}}^{\textrm{out}})$
such that

$\displaystyle\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\left(\sum_{s\in\mathcal{S}}\left|\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i,s}(t)-z_{i,s}(t)\right)\right|\right)\right]\leq$
$\displaystyle C_{1}\left(\sqrt{w_{\textrm{max}}}+\frac{1}{\sqrt{N}}\right)$
(26)

and

$\displaystyle\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\left(\sum_{s\in\mathcal{S}}\left|\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i,s}(t)-z_{i,s}(t)\right)\right|\right)\right]\leq$
$\displaystyle
C_{2}\left(\sqrt{\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}w_{ij}^{2}}+\frac{1}{\sqrt{N}}\right)$
(27)

where $\mu$ is the same as for Theorem 2.

Figure 2: The ratio of infected based on the average of $1000$ simulations
(triangles) compared to the estimate of NIMFA (solid line) on an $N=1000$
vertex modified cycle graphs with the closest $10$ (left) and $100$ (right)
neighbors being connected. ($\beta=2,\gamma=1$) As we increase the degrees
NIMFA performs better.

Figure 3: The ratio of infected based on the average of $10$ simulations
(triangles) compared to the estimate of NIMFA (solid line) on an $N=5000$
vertex modified cycle graphs with the closest $10$ (left) and $100$ (right)
neighbors being connected. ($\beta=2,\gamma=1$) As we increase the degrees
NIMFA performs better.

### Related works

In this section we compare our results to the recent independent work of
Sridhar and Kar [30, 31] and Parasnis et al. [24].

In [30] the authors describe how the state densities of certain related
stochastic processes on weighted graphs with doubly symmetric matrix $W$ can
be approximated by a set of $O(N)$ ODEs analogous to NIMFA given that the
normalized Frobenius norm $\frac{1}{N}\sum_{i=1}\sum_{j=1}^{N}w_{ij}^{2}$ is
small and $N$ is large.

Given the conclusions of Theorem 4.2 in [30] and Theorem 5 in the present
paper are very similar in nature, it makes sense to compare the general setup,
the conditions, the conclusions and the technique directly to those in the
present paper.

Setup. Strictly speaking, the stochastic processes discussed in the present
paper and in [30, 31] are different. In our work, time is continuous while
[30] and [31] start from discrete time steps then speed up time. This is a
minor difference though, and with appropriate time scaling, the models in [30,
31] and the present paper define essentially the same object.

Conditions. In the present paper, we require only that the normalized degrees
are bounded. This is more general than the doubly stochastic $W$ assumption of
[30, 31]. Specifically, our result also justifies Example 4.2 in [31].

Via (27), qualitatively the same type of error terms were retained in terms of
the normalized Frobenius norm, but [30, 31] provides an error probability
bound that is exponential in $N$. In the present paper, we do not focus on
this kind of large deviation bound in $N$.

[30, 31] derive bounds for the global average. On the other hand, our results
show more localized, uniform bounds in terms of vertices. This is made
possible by the use of the auxiliary Markov processes $\hat{\xi}_{i}(t)$,
allowing accurate predictions about individual vertices too, not just global
averages.

Our framework also allows higher order interactions, while [30, 31] is
restricted to order 2 interactions (graphs).

In [24] the authors study the SIR process in age-structured populations on
time-varying networks. They show that when $N$ and the rewiring rate is high
the prevalence of the age groups can be described via an ODE system analogous
to the metapopulation NIMFA model (34) in Section 4.2. Note that [24] applies
to cases with fast, but finite rewiring rates as well, while our result only
considers the idealized case of infinite rewiring rates.

##  4 Further reductions to NIMFA

This section relates NIMFA to other approaches from the literature. Although
NIMFA is a major reduction of the exact Kolmogorov-equations, requiring only
$O(N)$ ODEs to be solved, it can be still computationally prohibitive when the
number of vertices is too large. Furthermore, NIMFA requires knowing both the
full network structure and precise initial conditions for all vertices. We
look at further reductions to (12) when additional structure is known for the
network or initial conditions; several of these actually lead to other well-
known models from the literature.

###  4.1 Homogeneous mean field approximation

The homogeneous mean field approximation (HMFA) assumes that the vertices are
_well mixed_ , meaning, every vertex interacts with every other with equal
weights. Formally, this can be this can be described by a complete hypergraph
(with all loops and secondary loops):

$\displaystyle w_{i,\underline{j}}^{(m)}=\frac{1}{N^{m}}.$

This definition may be generalized to include cases when
$w_{i,\underline{j}}^{(m)}=0$ for certain $m$ indices, e.g. $(M+1)$-uniform
hypergraphs. For ease of notation, instead of modifying the definition of
$w_{i,\underline{j}}^{(m)}$, it is also possible to choose the rate functions
$q_{ss^{\prime}}(\phi)$ so that they do not depend on the appropriate
$\phi^{(m)}$ coordinates, making the choice of $w_{i,\underline{j}}^{(m)}$
irrelevant.

It is easy to see that for such networks, $w_{\max}=\frac{1}{N}$ and
$\delta_{\max}=1$. What remains to show is that (5) holds with some bounded
$R$.

$\displaystyle\begin{split}&\sum_{\begin{subarray}{c}\underline{j}\in[N]^{m}\\\
\underline{j}\textit{ is s.\
loop}\end{subarray}}w_{i,\underline{j}}^{(m)}=\frac{1}{N^{m}}\left|\left\\{\left.\underline{j}\in[N]^{m}\right|\underline{j}\textrm{
s.\ loop}\right\\}\right|=\\\
&1-\frac{1}{N^{m}}\left|\left\\{\left.\underline{j}\in[N]^{m}\right|\underline{j}\textrm{
not s. loop}\right\\}\right|=1-\prod_{l=0}^{m-1}\left(1-\frac{l}{N}\right)=\\\
&O\left(\frac{1}{N}\right)\ll\frac{1}{\sqrt{N}}=\sqrt{w_{\max}},\end{split}$
(28)

hence, $R$ can be chosen arbitrarily small for large enough $N$.

Our goal now is to derive a small system of equations for

$\displaystyle u(t):=\frac{1}{N}\sum_{i=1}^{N}z_{i}(t).$

Our strategy is based on the observation that the neighbourhood vectors
$\zeta_{i}(t)$ are the same for all vertices.

$\displaystyle\zeta_{i,\underline{s}}^{(m)}(t)=$
$\displaystyle\frac{1}{N^{m}}\sum_{\underline{j}\in[N]^{m}}\prod_{l=1}^{m}z_{j_{l},s_{l}}(t)=\prod_{l=1}^{m}\left(\frac{1}{N}\sum_{j_{l}=1}^{N}z_{j_{l},s_{l}}(t)\right)=$
$\displaystyle\prod_{l=1}^{m}u_{s_{l}}(t)=:u_{\underline{s}}^{(m)}(t)$

This results in the ODE system:

$\displaystyle\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}u(t)=&Q\left(U(t)\right)u(t),\\\
U(t)=&\left(u^{(m)}(t)\right)_{m=1}^{M},\\\
u^{(m)}(t)=&\left(u_{\underline{s}}^{(m)}(t)\right)_{\underline{s}\in\mathcal{S}^{m}}=\left(\prod_{l=1}^{m}u_{s_{l}}(t)\right)_{\underline{s}\in\mathcal{S}^{m}}.\end{split}$
(29)

For example, the simplicial SIS model (29) takes the form

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}u_{I}(t)=-\gamma
u_{I}(t)+\left(1-u_{I}(t)\right)\sum_{m=1}^{M}\beta_{m}u_{I}^{m}(t).$

which was used in [13].

In this setting, Theorem 5 shows the ratio of vertices in state
$s\in\mathcal{S}$ can be approximated by $u_{s}(t)$ with
$O\left(\frac{1}{\sqrt{N}}\right)$ error. The well known results of Kurtz [16,
17] correspond to the $M=1$ case.

#### Regular hypergraphs

Although (29) is both feasible for analytical and numerical investigations
(due to its finite size) the assumption that the network structure is well-
mixed is quite restrictive. However, as we will see, the well-mixed condition
can be relaxed given uniform initial conditions.

We call a weighted hypergraph _regular_ if

$\displaystyle\forall\ 1\leq i\leq N,\ 1\leq m\leq M\ \ \delta^{(m)}(i)=1.$
(30)

Note that the value $1$ is arbitrary and any other constant value would work
with minor modifications to the rate functions $q_{ss^{\prime}}$.

We note that (30) always holds for Convention 2 hypergraphs. For Convention 1,
it holds when $d^{(m)}(i)=\bar{d}^{(m)}\,\,\forall 1\leq i\leq N,\ 1\leq m\leq
M$ (that is, the hypergraph is regular in the usual sense).

######  Proposition 1.

Assume (30) and

$z_{i}(0)=u(0)\quad\forall\ 1\leq i\leq N$

for some $u(0)\in\Delta^{\mathcal{S}}.$ Then the solution of (12) takes the
form

$\ z_{i}(t)=u(t)\quad\forall\ 1\leq i\leq N$

where $u(t)$ satisfies (29).

We mention that statements similar to Proposition 1 have appeared in the
literature before in certain special cases [15, Proposition 3.18 ]. Combining
Proposition 1 with Theorem 2 ensures the accuracy of the homogeneous mean
field approximation on regular graphs with large degrees and homogeneous
initial conditions disregarding any further network structure.

###### Proof.

(Proposition 1)

Let $u(t)$ be the solution of (29). Set $z_{i}(t)=u(t).$ We have to show that
$z_{i}(t)$ satisfies (12). The initial conditions are satisfied according to
the assumption, and for the derivatives,

$\displaystyle u_{\underline{s}}^{(m)}(t)=$ $\displaystyle
u_{\underline{s}}^{(m)}(t)\delta^{(m)}(i)=u_{\underline{s}}^{(m)}(t)\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}=\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}z_{\underline{j},\underline{s}}^{(m)}(t)=\zeta_{i,\underline{s}}^{(m)}(t),$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{i}(t)=$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}u(t)=Q\left(U(t)\right)u(t)=Q\left(\zeta_{i}(t)\right)z_{i}(t).$

∎

###  4.2 Metapopulation models

As we saw in Section 4.1 , a way to reduce the number of equations is by
grouping vertices together and representing them by a single averaged-out
term. In practice, this approach will only work if the vertices grouped
together are sufficiently homogeneous, which is typically not the case for the
entire population. To mitigate this issue, we may introduce _communities_ ,
inside which we assume homogeneity, then derive the dynamics between
communities. This ”higher resolution” may increase accuracy, at the cost of a
larger ODE system.

In practice, the communities can be chosen by demographic and geographic
criterion such as age and location. Alternatively, it is also possible to
group vertices according to degree, or a third option is the use of community
detection algorithms [1].

We present the general setup for metapopulation models first for graphs in
Section 4.2.1, then for hypergraphs in Section 4.2.2.

For the SIS process on graphs similar results had been derived in [4].

####  4.2.1 Metapopulation models on graphs

First, assume $M=1$. Divide the vertices into a partition $V_{1},\dots,V_{K}$
with size $\left|V_{k}\right|=N_{k}$ such that vertices inside a group are
similar in some sense. The average weight between group $V_{k}$ and $V_{l}$ is

$\displaystyle\tilde{w}_{kl}=\frac{\sum_{i\in V_{k}}\sum_{j\in
V_{l}}w_{ij}}{N_{k}N_{l}}.$ (31)

(In the idealized case of metapopulations, $w_{ij}$ would have the same value
$\tilde{w}_{kl}$ for each $i\in V_{k},j\in V_{l}$ pair.)

Next we derive the dynamics for the averages

$\displaystyle\bar{z}_{k}(t):=\frac{1}{N_{k}}\sum_{i\in V_{k}}z_{j}(t).$ (32)

$\zeta_{i}(t)$ has the same value $\bar{\zeta}_{k}(t)$ for all $i\in V_{k}$:

$\displaystyle\bar{\zeta}_{k}(t)=\zeta_{i}(t)=\sum_{j=1}^{N}w_{ij}z_{j}(t)=\sum_{l=1}^{K}\underbrace{N_{l}\tilde{w}_{kl}}_{\bar{w}_{kl}}\frac{1}{N_{l}}\sum_{j\in
V_{l}}z_{j}(t)=\sum_{l=1}^{K}\bar{w}_{kl}\bar{z}_{l}(t).$ (33)

Therefore, we can derive an ODE system for (32)

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bar{z}_{k}(t)=Q\left(\bar{\zeta}_{k}(t)\right)\bar{z}_{k}(t)$
(34)

which is equivalent to (12) on the graph $\overline{\mathcal{G}}$ with vertex
set $\\{1,\dots,K\\}$ and weights $\left(\bar{w}_{kl}\right)_{k,l=1}^{K}.$

####  4.2.2 Metapopulation models on hypergraphs

For the general metapopulation setting, we assume that for each $m=1,\dots,M$,
the population is partitioned into _local groups_
$V_{1}^{(m)},\dots,V_{K^{(m)}}^{(m)}$. The _type_ of a vertex will be denoted
by $k=\left(k^{(1)},\dots,k^{(M)}\right)$, which means that for each
$m=1,\dots,M$, the given vertex is in the local group $V_{k^{(m)}}^{(m)}$.
Vertices can be partitioned according to their type into
$\prod_{m=1}^{M}K^{(m)}$ _global groups_.

We aim to define a hypergraph on the types, with weights consistent with the
average of weights within each group. That said, with the above setup, this is
easier to do using local groups for each $m=1,\dots,M$.

For a given $m$, $k^{(m)}$ and
$\underline{l}^{(m)}=\left(l_{1}^{(m)},\dots,l_{m}^{(m)}\right)$, the _total
local $m$-weight between $k^{(m)}$ and $\underline{l}^{(m)}$_ is defined as

$\displaystyle W_{k^{(m)},\underline{l}^{(m)}}^{(m)}:=\sum_{i\in
V_{k^{(m)}}^{(m)}}\sum_{j_{1}\in V_{l_{1}^{(m)}}^{(m)}}\dots\sum_{j_{m}\in
V_{l_{m}^{(m)}}^{(m)}}w_{i,\underline{j}}^{(m)}.$ (35)

Then, using the notation

$N_{\underline{l}^{(m)}}:=\prod_{r=1}^{M}N_{l_{r}^{(m)}}^{(m)},$

we define the weight of the edge containing the local groups
$k^{(m)},\underline{l}^{(m)}$ as

$\displaystyle\tilde{w}_{k^{(m)},\underline{l}^{(m)}}^{(m)}:=\frac{W_{k^{(m)},\underline{l}^{(m)}}^{(m)}}{N_{k^{(m)}}N_{\underline{l}^{(m)}}}.$
(36)

Let $k(i)=\left(k^{(1)}(i),\dots,k^{(M)}(i)\right)$ denote the type of $i$.
For easier notation, we will often use $\iota\sim U\left([N]\right)$, which is
a random vertex independent from everything else. Then we define the average
of $z_{i}(t)$ over type $k$ as

$\displaystyle\bar{z}_{k}(t):=\mathbb{E}\left(\left.z_{\iota}(t)\right|k(\iota)=k\right)=\frac{1}{N_{k}}\sum_{i\in
V_{k}}z_{i}(t).$ (37)

In this case as well, $\zeta_{i}(t)$ has the same value for all $i\in V_{k}$;
this common value will be denoted by $\bar{\zeta}_{k}(t).$ Let
$\iota_{1},\dots,\iota_{m}$ denote i.i.d. copies of $\iota.$ Then

$\displaystyle\begin{split}\bar{\zeta}_{k}^{(m)}(t)=&\zeta_{i}^{(m)}(t)=\sum_{j\in[N]^{m}}w_{i,\underline{j}}^{(m)}z_{\underline{j}}^{(m)}(t)=\sum_{\underline{l}^{(m)}}\tilde{w}_{k,\underline{l}^{(m)}}^{(m)}\sum_{j_{1}\in
V_{l_{1}^{(m)}}^{(m)}}\dots\sum_{j_{m}\in
V_{l_{m}^{(m)}}^{(m)}}z_{\underline{j}}^{(m)}(t)\\\
=&\sum_{\underline{l}^{(m)}}\underbrace{N_{\underline{l}^{(m)}}\tilde{w}_{k,\underline{l}^{(m)}}^{(m)}}_{:=\bar{w}_{k^{(m)},\underline{j}^{(m)}}^{(m)}}\mathbb{E}\left(\left.\prod_{r=1}^{m}z_{\iota_{r}}(t)\right|k^{(m)}(\iota_{1})=l_{1}^{(m)},\dots,k^{(m)}(\iota_{m})=l_{m}^{(m)}\right)=\\\
=&\sum_{\underline{l}^{(m)}}\bar{w}_{k^{(m)},\underline{j}^{(m)}}^{(m)}\prod_{r=1}^{m}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{.}}\end{split}$
(38)

This means that the ODE system for (37) is formally the same as (34) (with the
appropriate definition of $\bar{z}_{k}(t)$ and $\bar{\zeta}_{k}(t)$).

Note that $\bar{\zeta}_{k}(t)$ can also be expressed via $\bar{z}_{k}(t)$ as

$\displaystyle\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)=$
$\displaystyle\mathbb{E}\left(\left.\mathbb{E}\left(\left.z_{\iota}(t)\right|k(i)=k\right)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)=$
$\displaystyle\mathbb{E}\left(\left.\bar{z}_{k(\iota)}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right),$

making (34) a closed system.

In the special case when the hypergraph is $(M+1)$-uniform, we can set
$K^{(m)}=1$ for all $m<M$ virtually making the local group $k^{(M)}$ and the
global group $k$ the same (apart from some $1$’s in the first $M-1$
components). In this case, $Q$ only depends on $\bar{\zeta}^{(M)}(t)$ which
can be expressed as

$\displaystyle\bar{\zeta}_{k^{(M)}}^{(M)}=\sum_{\underline{l}^{(m)}}\bar{w}_{k^{(m)},\underline{l}^{(m)}}^{(m)}\prod_{r=1}^{m}\bar{z}_{k^{(M)}(l_{r})}(t).$

###  4.3 Annealed networks

So far, we only focused on the dynamics of the Markov process neglecting the
dynamics of the network itself. When there is a separation of scale between
the speed of the Markov process and the changes to the network itself, two
kinds of idealizations are typically used:

  * •

_quenched networks_ : the speed at which the network changes is much slower
than the Markov process. In this case, the network is assumed constant in
time.

  * •

_annealed networks_ : the speed at which the network changes is much faster
than the Markov process. In this case, we consider the network changes
averaged out for the interactions.

Annealed networks can be modeled by replacing connections
$a_{i,\underline{j}}^{(m)}$ in (2) and (3) with the average $\langle
a\rangle_{i,\underline{j}}^{(m)}$.

In this section, we present a setup for annealed networks generated via the
configuration model [20]. Similar calculations can be made for other models
that include e.g. degree correlation such as equation (93) in [8].

Once again, we start with the graph case.

In the configuration model the degrees $d(1),\dots,d(N)$ are given beforehand,
and vertex $i$ receives $d(i)$ half-edges (_stubs_) initially. Then in each
round, we choose two stubs at random to connect and form an edge, repeating
this procedure until all stubs are paired.

Loops and multiple edges are possible, but their effect will be neglected. The
expected connection between vertices $i$ and $j$ is

$\displaystyle\langle a\rangle_{ij}=\frac{d(i)d(j)}{\bar{d}N}.$

The degree of each vertex $i$ indeed matches the prescribed $d(i)$ as

$\displaystyle\sum_{j=1}^{N}\langle
a\rangle_{ij}=\frac{d(i)}{\bar{d}}\frac{1}{N}\sum_{j=1}^{N}d(j)=d(i).$

$\langle a\rangle_{ij}$ depends only on the degrees of $i$ and $j$, so it can
be interpreted as a metapopulation model where vertices are grouped according
to their degree. (Note that here we also use the index $k=0$ for isolated
vertices if any.) The corresponding weights are

$\displaystyle\tilde{w}_{kl}=\frac{kl}{\bar{d}^{2}N},$

for Convention 1, and

$\displaystyle\tilde{w}_{kl}=\frac{l}{\bar{d}N}.$

for Convention 2.

Let $q_{k}:=\frac{kN_{k}}{\bar{d}N}$ denote the size biased degree
distribution and introduce

$\displaystyle\Theta(t):=\sum_{l=0}^{d_{\max}}q_{l}\bar{z}_{l}(t).$ (39)

Using (33), $\bar{\zeta}_{k}(t)$ can be written as

$\displaystyle\bar{\zeta}_{k}(t)=\frac{k}{\bar{d}}\Theta(t),$

for Convention 1, and

$\displaystyle\bar{\zeta}_{k}(t)=\Theta(t).$

for Convention 2.

For example, the I component of the SIS process assuming Convention 1 is

$\displaystyle\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}\bar{z}_{k,I}(t)&=-\gamma\bar{z}_{k,I}(t)+\frac{\beta}{\bar{\bar{d}}}k\left(1-\bar{z}_{k,I}(t)\right)\Theta_{I}(t){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{,}}\\\
\Theta_{I}(t)&=\sum_{l=0}^{d_{\textrm{max}}}q_{l}\bar{z}_{l,I}(t){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{.}}\end{split}$
(40)

which is the Inhomogeneous Mean Field Approximation (IMFA) studied by Pastor-
Satorras and Vespignani [25].

For Convention 1, to apply the results of the present paper, we need to assume
upper regularity, i.e. $\delta_{\max}=\frac{d_{\max}}{\bar{d}}$ to be bounded.
In many applications, the degree distribution converges to a fixed
distribution, making $\bar{d}$ bounded; in such a setting, we accordingly
require $d_{\max}$ to be bounded as well.

Assuming upper regularity,

$\displaystyle
w_{\max}=\frac{d_{\max}^{2}}{\bar{d}^{2}N}=\frac{1}{N}\delta_{\max}^{2}$

thus Theorem 5 actually provides an $O\left(\frac{1}{\sqrt{N}}\right)$ error
bound.

As for Convention 2, $\delta_{\max}=1$ holds as usual, and

$\displaystyle w_{\max}=\frac{1}{N}\frac{d_{\max}}{\bar{d}}.$

Unfortunately, one can not relax the bound on $d_{\textrm{max}}$ by using (23)
instead of (22) as it requires bounds for the out-degrees:

$\displaystyle\delta^{\textrm{out}}(j)=\sum_{i=1}^{N}w_{ij}=\sum_{k=0}^{d_{\textrm{max}}}N_{k}\frac{d(j)}{\bar{d}N}=\frac{d(j)}{\bar{d}}\leq\frac{d_{\textrm{max}}}{\bar{d}}\leq\delta_{\textrm{max}}^{\textrm{out}}.$

Now we turn to the hypergraph case $M>1$. We generalize the notion of the
configuration model in the following manner: For a fixed $m$, the $m$-degrees
are given as $d^{(m)}(1),\dots,d^{(m)}(N)$ and each vertex receives $m$-stubs
based on their degree. In each round, we choose $m+1$ $m$-stubs at random to
form an $m$-edge, then repeat this procedure until all of the stubs have been
paired. This procedure is performed for each $1\leq m\leq M$ independently.

For distinct $i,j_{1},\dots j_{m}$, the probability of connecting them in a
given round is

$\displaystyle\frac{d^{(m)}(i)\prod_{r=1}^{m}d^{(m)}(j_{r})}{{\bar{d}^{(m)}N\choose
m+1}}\approx\frac{(m+1)!d^{(m)}(i)\prod_{r=1}^{m}d^{(m)}(j_{r})}{\left(\bar{d}^{(m)}N\right)^{m+1}}.$

Since there are $\frac{\bar{d}^{(m)}N}{m+1}$ rounds in total, we set

$\displaystyle\langle
a\rangle_{i,\underline{j}}^{(m)}:=\frac{m!d^{(m)}(i)\prod_{r=1}^{m}d^{(m)}(j_{r})}{\left(\bar{d}^{(m)}N\right)^{m}}.$

For the hypergraph case, we only examine Convention 1, for which

$\displaystyle\tilde{w}_{k^{(m)},\underline{l}^{(m)}}^{(m)}=\frac{k^{(m)}}{\bar{d}^{(m)}}\frac{\prod_{r=1}^{m}l_{r}^{(m)}}{\left(\bar{d}^{(m)}N\right)^{m}}.$

Once again, the resulting hypergraph can be interpreted as a metapopulation
model, where the local groups are given according to the $m$-degrees of the
vertices.

Clearly $\delta^{(m)}(i)=\frac{d^{(m)}(i)}{\bar{d}^{(m)}},$ so we make an
upper regularity assumption in this case as well, from which
$w_{\max}=O\left(\frac{1}{N}\right)$ follows.

For hypergraphs, we also need to check the condition (5).

$\tilde{w}_{k^{(m)},\underline{l}^{(m)}}^{(m)}\leq\frac{\delta_{\max}^{m+1}}{N^{m}},$

so (28) implies

$\displaystyle\sum_{\begin{subarray}{c}\underline{j}\in[N]^{m}\\\
\underline{j}\textit{ is s.\ loop}\end{subarray}}w_{i,\underline{j}}^{(m)}\leq
C\sum_{\begin{subarray}{c}\underline{j}\in[N]^{m}\\\ \underline{j}\textit{ is
s.\
loop}\end{subarray}}\frac{1}{N^{m}}=O\left(\frac{1}{N}\right)\ll\sqrt{w_{\max}},$
(41)

hence arbitrarily small $R$ can be used for large enough $N$.

The next step is to calculate $\bar{\zeta}_{k}(t)$ based on (34). Define

$q_{k^{(m)}}^{(m)}:=\frac{k^{(m)}N_{k^{(m)}}}{\bar{d}^{(m)}N},$

the size-biased degree distribution of the $m$-vertices. Also define

$\displaystyle\Theta^{(m)}(t):=\sum_{l=1}^{d_{\max}^{(m)}}q_{l}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|d^{(m)}(\iota)=l\right),$
(42)

once again using the notation $\iota\sim U\left([N]\right)$.

Using (38),

$\displaystyle\bar{\zeta}_{k}^{(m)}(t)$
$\displaystyle=\sum_{\underline{l}^{(m)}}\bar{w}_{k^{(m)},\underline{j}^{(m)}}^{(m)}\prod_{r=1}^{m}\mathbb{E}\left(\left.z_{\iota}(t)\right|d^{(m)}(\iota)=l_{r}^{(m)}\right)=$
$\displaystyle=\frac{k^{(m)}}{\bar{d}^{(m)}}\sum_{\underline{l}^{(m)}}\prod_{r=1}^{m}q_{l_{r}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|d^{(m)}(\iota)=l_{r}\right)$
$\displaystyle=\frac{k^{(m)}}{\bar{d}^{(m)}}\prod_{r=1}^{m}\sum_{l_{r}=1}^{d_{\textrm{max}}^{(m)}}q_{l_{r}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|d^{(m)}(\iota)=l_{r}\right)=$
$\displaystyle=\frac{k^{(m)}}{\bar{d}^{(m)}}\left(\Theta^{(m)}(t)\right)^{m}.$

Accordingly, e.g. the dynamics for the simplicial SIS model can be written as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bar{z}_{k,I}(t)=-\gamma\bar{z}_{k,I}(t)+(1-\bar{z}_{k,I}(t))\sum_{m=1}^{M}\frac{\beta^{(m)}}{\bar{d}^{(m)}}\left(\Theta_{I}^{(m)}(t)\right)^{m}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{.}}$
(43)

(43) was studied in [14] for the $(M+1)$-uniform case, where
$\mathbb{E}\left(\left.z_{\iota}(t)\right|d^{(M)}(\iota)=l\right)$ simplifies
to $\bar{z}_{k}(t)$ as the global class $k$ and the local class $k^{(M)}$
coincide.

###  4.4 Activity-driven networks

Activity-driven networks were introduced in [26].

Let $a_{1},\dots,a_{K}$ be positive numbers called _activities_ and let $a(i)$
denote the activity of vertex $i$. Instead of a graphs structure, each vertex
chooses a random vertex uniformly with rate $\beta a(i)$ and if they are an SI
pair, the susceptible node becomes infected. Recoveries happen independently
with rate $\gamma.$

The above model corresponds to an SIS process on the weighted graph

$\displaystyle w_{ij}=\frac{a(i)+a(j)}{N}$

since to form the $(i,j)$ pair, either $i$ or $j$ needs to activate, and each
vertex is chosen with probability $\frac{1}{N}.$ The graph is a metapopulation
model, with groups corresponding to the activity values.

We generalize this concept to allow higher order interactions.
$a_{1}^{(m)},\dots,a_{K^{(m)}}^{(m)}$ are the possible $m$-activities and we
assume that vertex $i$ chooses $m$ other vertices at random with rate
$a^{(m)}(i).$ This results in a hypergraph with weights

$\displaystyle
w_{i,\underline{j}}^{(m)}=\frac{1}{N^{m}}\left(a_{i}^{(m)}+\sum_{r=1}^{m}a_{j_{r}}^{(m)}\right).$

Assume the activity rates are bounded from above by some $a_{\max}<\infty.$
Also, introduce

$\bar{a}^{(m)}:=\frac{1}{N}\sum_{i=1}^{N}a^{(m)}(i).$

Then

$\displaystyle\delta^{(m)}(i)=a_{i}^{(m)}+\frac{1}{N^{m}}\sum_{\underline{j}\in[N]^{m}}\sum_{r=1}^{m}a_{j_{r}}^{(m)}=a_{i}^{(m)}+\bar{a}^{(m)}\leq
2a_{\max}$

so (4) is satisfied.

$w_{\max}\asymp\frac{1}{N}$ and (41) is applicable here as well satisfying
(5), hence Theorem 2 applies.

$\bar{\zeta}_{k}(t)$ can also be expressed with the help of (38).

######  Proposition 2.

Let $\iota\sim U([N])$ a random index and $p^{(m)}_{k^{(m)}}$ be the ratio of
vertices in the local group $k^{(m)}.$ Also, define

$\displaystyle\psi^{(m)}(t):=\sum_{l=1}^{K^{(m)}}a_{l}^{(m)}p_{l}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|a^{(m)}(\iota)=l^{(m)}\right).$

Then the neighborhood vectors have the form

$\displaystyle\bar{\zeta}^{(m)}_{k}(t)=\left(a_{k^{m}}^{(m)}\mathbb{E}\left(z_{\iota}(t)\right)+\psi^{(m)}(t)\right)\mathbb{E}^{m-1}\left(z_{\iota}(t)\right).$

The proof of Proposition 2 is given in Section 6.

For activity-driven networks, the simplicial SIS model takes the form

$\displaystyle\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}\bar{z}_{k,I}(t)=&-\gamma\bar{z}_{k,I}(t)+\left(1-\bar{z}_{k,I}(t)\right)\cdot\\\
&\sum_{m=1}^{M}\beta_{m}\mathbb{E}^{m-1}\left(z_{\iota,I}(t)\right)\left(a_{k^{m}}^{(m)}\mathbb{E}\left(z_{\iota,I}(t)\right)+\psi_{I}^{(m)}(t)\right).\end{split}$
(44)

[33] proves that (44) describes the large graph limit correctly when $M=1$.

###  4.5 Dense graphs and Szemerédi’s regularity lemma

We call a hypegraph dense if there is some $0<p_{0}\leq 1$ such that

$\displaystyle\bar{d}^{(m)}\geq p_{0}N^{m}\quad\forall\ 1\leq m\leq M.$ (45)

For Convention 1 graphs,

$\displaystyle\frac{1}{M!\ N}\leq w_{\max}\leq$
$\displaystyle\frac{1}{p_{0}N},$ $\displaystyle\delta_{\max}\leq$
$\displaystyle\frac{1}{p_{0}}$

hold and (41) directly follows, satisfying the conditions for Theorem 2.

We focus on the graph case $M=1$. We assume that the rate functions
$q_{ss^{\prime}}$ are _affine_ , that is, they have the form

$\displaystyle
q_{ss^{\prime}}\left(\phi\right)=q_{ss^{\prime}}^{(0)}+\sum_{r\in\mathcal{S}}q_{ss^{\prime},r}^{(1)}\phi_{r}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{,}}$
(46)

where
$q_{ss^{\prime}}^{(0)},\left(q_{ss^{\prime},r}^{(1)}\right)_{r\in\mathcal{S}}$
are nonnegative constants. Many epidemiological models have this form,
including the SIS process.

As it was pointed out in the preliminary work [9], Szemerédi’s regularity
lemma [32] provides a method to approximate (12) with a finite system up to
arbitrary precision (for large enough $N$).

Roughly speaking, Szemerédi’s regularity lemma states that any large enough
dense graph can be partitioned into finitely many “boxes” (called an
$\varepsilon$-regular partition) which have the same size (except one
remainder box), and besides a few exceptional pairs the edge count between two
boxes behaves as if coming from a randomly mixed graph, with error at most
$\varepsilon$.

We denote an $\varepsilon$-regular partition by $V_{0},V_{1},\dots,V_{K}$,
where $V_{0}$ is the exceptional set.

$\displaystyle e(A,B):=\sum_{i\in A}\sum_{j\in B}a_{ij}$

refers to the number of edges between the vertex sets $A,B$ with the
convention that edges in $A\cap B$ are counted double.

We define the graph $\overline{\mathcal{G}}$ on vertices
$\left(V_{1},\dots,V_{K}\right)$. ($V_{0}$ is neglected.)

The adjacency matrix is replaced by the edge density between $A,B\subseteq[N]$
defined as

$\displaystyle\rho(A,B):=\frac{e\left(A,B\right)}{\left|A\right|\cdot\left|B\right|}$
(47)

It is easy to see that $0\leq\rho\left(A,B\right)\leq 1.$

The adjacency matrix counterpart for $\overline{\mathcal{G}}$ is simply the
edge density between the $V_{1},\dots,V_{K}$ sets. For the average degree we
further define

$\displaystyle p:=$ $\displaystyle\frac{\bar{d}}{N},$ (48)
$\displaystyle\kappa:=$
$\displaystyle\frac{\left|V_{1}\right|}{N}=\dots=\frac{\left|V_{K}\right|}{N}$
(49)

where $p$ is the global edge density of $\mathcal{G}$ and $\kappa$ is the
portion of vertices one box contains. The average degree in
$\overline{\mathcal{G}}$ is $Kp\approx\frac{p}{\kappa},$ motivating the
definition of the weights

$\displaystyle\bar{w}_{kl}:=\frac{\kappa}{p}\rho\left(V_{k},V_{l}\right).$
(50)

The corresponding solution of (12) on the graph $\overline{\mathcal{G}}$ with
weights (50) is denoted by $\left(v_{k}(t)\right)_{k=1}^{K}$ with initial
condition

$\displaystyle v_{k}(0)=\frac{1}{|V_{k}|}\sum_{i\in V_{k}}z_{i}(0).$ (51)

Finally, we define

$\displaystyle\bar{v}(t):=\sum_{k=1}^{K}\frac{\left|V_{k}\right|}{N}v_{k}(t)$
(52)

and the average global density vector

$\displaystyle\bar{z}(t):=\frac{1}{N}\sum_{i=1}^{N}z_{i}(t).$ (53)

######  Theorem 6.

$\forall T>0,\varepsilon>0,p_{0}>0\,\exists K_{\max}\in\mathbb{Z}^{+}$ such
that for any $\mathcal{G}$ simple graph with density parameter $p_{0}$ and
$N\geq K_{\textrm{max}}$, there exists a partition $V_{0},V_{1},\dots,V_{K}$
with $K\leq K_{\max}$ such that

  * •

$\left|V_{1}\right|=\dots=\left|V_{K}\right|{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{,}}$

  * •

$\left|V_{0}\right|\leq\varepsilon N$,

  * •

$\sup_{0\leq t\leq T}\left\|\bar{z}(t)-\bar{v}(t)\right\|_{1}\leq\varepsilon.$

The proof is provided in Section 6.

Szemerédi’s regularity lemma also guarantees that such a partition can be
found in polynomial time [1].

We note that $K_{\max}$ may increase rapidly as $\varepsilon\to 0^{+}$
limiting the applicability of the approach. That said, for networks with extra
community structure, this approach may still be useful.

##  5 Discussion

In this paper we examined the accuracy of the so called N-Intertwined Mean
Field Approximation on hypergraphs. The idea of NIMFA is to assume vertices
are independent from each other, then derive the dynamics of the occupation
probabilities of each vertex. This leaves us with and ODE system of size
$O(N)$ instead of an exponentially increasing system given by the exact
Kolmogorov equations.

Our findings show that when the incoming weights are well distributed – for
example, vertices typically have large degrees – then NIMFA gives an accurate
approximation. Under additional assumptions we showed how the number of ODEs
can be further reduced to give well-known approximation methods from the
literature, such as the heterogenous mean field approximation. Finally, we
showed how Szemerédy’s regularity lemma can be used to reduce the number of
equations to constant order (depending only on the error desired) for large
enough dense graphs.

These results have their limitations. The error bounds work poorly for truly
sparse graphs (with bounded average degrees). Analyzing such systems probably
requires qualitatively different approaches.

The upper regularity condition can be restrictive for certain applications. We
conjecture that the results could be greatly generalized in this direction for
degree distributions with fast decaying tails.

For the reduction for dense graph we applied the strong version of Szemerédy’s
lemma. The weak version of Szemerédy’s lemma, however, has more desirable
algorithmic properties and a smaller bound on the number of ”boxes” one needs
for a given $\varepsilon$. Extending the theorem in this direction might be
beneficial for large, inhomogeneous, dense systems.

Finally, NIMFA has the disadvantage of requiring full knowledge of the network
which is usually not possible in practice. Using metapopulation networks
instead mitigates this problem, and also greatly reduces the number of
equations required. This method, however, relies on the assumption that the
metapopulation dynamics is close enough to the original one. Further research
is needed to understand how well coarse graining performs in terms of
preserving the network dynamics.

##  6 Proofs

###  6.1 General proofs

We state and prove a technical lemma first which will be used throughout other
proofs.

######  Lemma 1.

Let $a_{1},\dots,a_{n}$ and $b_{1},\dots,b_{n}$ two sets of numbers such that
$0\leq\left|a_{i}\right|,\left|b_{i}\right|\leq 1$. Then

$\displaystyle\left|\prod_{i=1}^{n}a_{i}-\prod_{i=1}^{n}b_{i}\right|\leq\sum_{i=1}^{n}\left|a_{i}-b_{i}\right|.$

###### Proof.

(Lemma 1)

The proof is by induction on $n$. The statement is trivial for $n=1$. For
$n>1$,

$\displaystyle\left|\prod_{i=1}^{n}a_{i}-\prod_{i=1}^{n}b_{i}\right|=$
$\displaystyle\left|a_{n}\prod_{i=1}^{n-1}a_{i}-b_{n}\prod_{i=1}^{n-1}b_{i}\right|=$
$\displaystyle\left|\left(a_{n}-b_{n}\right)\prod_{i=1}^{n-1}a_{i}+b_{n}\left(\prod_{i=1}^{n-1}a_{i}-\prod_{i=1}^{n-1}b_{i}\right)\right|\leq$
$\displaystyle\left|a_{n}-b_{n}\right|\prod_{i=1}^{n-1}\left|a_{i}\right|+\left|b_{n}\right|\cdot\left|\prod_{i=1}^{n-1}a_{i}-\prod_{i=1}^{n-1}b_{i}\right|\leq$
$\displaystyle\left|a_{n}-b_{n}\right|+\left|\prod_{i=1}^{n-1}a_{i}-\prod_{i=1}^{n-1}b_{i}\right|\leq\left|a_{n}-b_{n}\right|+\sum_{i=1}^{n-1}\left|a_{i}-b_{i}\right|=$
$\displaystyle\sum_{i=1}^{n}\left|a_{i}-b_{i}\right|.$

∎

Next we show that (12) exhibits a unique global solution.

###### Proof.

(Theorem 1)

The right hand side of (12) is locally Lipschitz, so there is a unique local
solution.

Instead of $q_{ss^{\prime}}$, we use the modified rate functions

$\displaystyle\hat{q}_{ss^{\prime}}(\phi):=$
$\displaystyle\left|q_{ss^{\prime}}(\phi)\right|$ (54)
$\displaystyle\hat{q}_{ss}(\phi)=$ $\displaystyle-\sum_{s^{\prime}\neq
s}\hat{q}_{s^{\prime}s}(\phi)$

which are nonnegative for any input; note that
$\left.\hat{q}_{ss^{\prime}}(\phi)\right|_{\phi\geq
0}=\left.q_{ss^{\prime}}(\phi)\right|_{\phi\geq 0}.$

The modified version of (12) is

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\hat{z}_{i}(t)=\hat{Q}\left(\hat{\zeta}_{i}(t)\right)\hat{z}_{i}(t)$

where
$\hat{Q}(\phi)=\left(\hat{q}_{ss^{\prime}}(\phi)\right)_{s,s^{\prime}\in\mathcal{S}}.$
The local solution uniquely exist in this case as well, and it either extends
to a global solution or blows up at a finite time.

Assume that the local solution blows up at time $t_{0}$. Then
$\hat{\zeta}_{i}(t)$ is well-defined for any $t<t_{0}.$

We construct an auxiliary time-inhomogeneous Markov process on $[0,t_{0})$.
The state space is $\mathcal{S}$ and the transition rates at time $t$ are
given by the matrix $\hat{Q}\left(\hat{\zeta}_{i}(t)\right)$. $p_{s}(t)$
denotes the probability of being in state $s\in\mathcal{S}.$ The Kolmogorov
equations have the form

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(t)=\hat{Q}\left(\hat{\zeta}_{i}(t)\right)p(t).$

Since $\hat{Q}\left(\hat{\zeta}_{i}(t)\right)$ is continuous for $t<t_{0}$,

$\max_{0\leq\tau\leq
t}\left\|\hat{Q}\left(\hat{\zeta}_{i}(\tau)\right)\right\|$

exists and is finite.

Based on Grönwall’s inequality,

$\displaystyle\hat{z}_{i}(t)-p(t)=$
$\displaystyle\hat{z}_{i}(0)-p(0)+\int_{0}^{t}\hat{Q}\left(\hat{\zeta}_{i}(u)\right)\left[\hat{z}_{i}(\tau)-p(\tau)\right]\mathrm{d}\tau,$
$\displaystyle\left\|\hat{z}_{i}(t)-p(t)\right\|=$
$\displaystyle\left\|\hat{z}_{i}(0)-p(0)\right\|+\sup_{0\leq u\leq
t}\left\|\hat{Q}\left(\hat{\zeta}_{i}(\tau)\right)\right\|\int_{0}^{t}\left\|\hat{z}_{i}(\tau)-p(\tau)\right\|\mathrm{d}\tau,$
$\displaystyle\sup_{0\leq\tau\leq
t}\left\|\hat{z}_{i}(\tau)-p(\tau)\right\|\leq$
$\displaystyle\left\|\hat{z}_{i}(0)-p(0)\right\|\exp\left(\sup_{0\leq\tau\leq
t}\left\|\hat{Q}\left(\hat{\zeta}_{i}(\tau)\right)\right\|\cdot t\right).$

Choosing $p(0)=\hat{z}_{i}(0)$ shows that $\hat{z}_{i}(t)=p(t)$ for any $0\leq
t<t_{0}$ as well.

But $p(t)$ is a probability vector, that is,
$\hat{z}_{i}(t)\in\Delta^{\mathcal{S}}$, which contradicts $\hat{z}_{i}(t)$
blowing up as $t\to t_{0}$, so the solution must be global.

Since the solution is on the simplex $\Delta^{S}$, we have
$\hat{q}_{ss^{\prime}}\left(\hat{\zeta}_{i}(t)\right)=q_{ss^{\prime}}\left(\hat{\zeta}_{i}(t)\right)$
(that is, the absolute values in (54) are not necessary). Therefore
$\hat{z}_{i}(t)$ is a solution for the original equation (12) as well. Since
the solution for (12) is unique, $\hat{z}_{i}(t)=z_{i}(t).$ This makes
$z_{i}(t)$ a global solution with values on the simplex $\Delta^{S}.$ ∎

###  6.2 Proof of Theorem 2

The strategy of the proof is to derive an inequality for $D_{\max}(t)$ and
$\tilde{D}_{i}(t)$ such that Grönwall’s inequality could be applied.

In the first step, we are showing an inequality for the error of the
indicators.

######  Lemma 2.

There exists $\tilde{C}_{1}=\tilde{C}_{1}(\delta_{\max})$ such that

$\displaystyle D_{\max}^{(0)}(t)\leq$
$\displaystyle\tilde{C}_{1}\int_{0}^{t}D_{\max}(\tau)\mathrm{d}\tau,$
$\displaystyle\tilde{D}_{i}^{(0)}(t)\leq$
$\displaystyle\tilde{C}_{1}\int_{0}^{t}\tilde{D}_{i}(\tau)\mathrm{d}\tau.$

###### Proof.

(Lemma 2) $\oplus$ denotes symmetric difference.

$\displaystyle\left|\xi_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right|\leq$
$\displaystyle\sum_{\begin{subarray}{c}s^{\prime}\in\mathcal{S}\\\
s^{\prime}\neq
s\end{subarray}}\left|\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(\tau)\right)-\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{K}_{i,ss^{\prime}}(\tau)\right)\right|+\left|\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(\tau)\right)-\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{K}_{i,ss^{\prime}}(\tau)\right)\right|\leq$
$\displaystyle\sum_{\begin{subarray}{c}s^{\prime}\in\mathcal{S}\\\
s^{\prime}\neq
s\end{subarray}}\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(\tau)\oplus\mathcal{K}_{i,ss^{\prime}}(\tau)\right)+\mathcal{N}_{i,s^{\prime}s}\left(\mathcal{H}_{i,s^{\prime}s}(\tau)\oplus\mathcal{K}_{i,s^{\prime}s}(\tau)\right)\leq$
$\displaystyle\sum_{\begin{subarray}{c}s^{\prime}\in\mathcal{S}\\\
s^{\prime}\neq
s\end{subarray}}\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(t)\oplus\mathcal{K}_{i,ss^{\prime}}(t)\right)+\mathcal{N}_{i,s^{\prime}s}\left(\mathcal{H}_{i,s^{\prime}s}(t)\oplus\mathcal{K}_{i,s^{\prime}s}(t)\right)$

In the last step we used the fact that
$\mathcal{H}_{i,ss^{\prime}}(\tau)\oplus\mathcal{K}_{i,ss^{\prime}}(\tau)$ is
an increasing set in $\tau$.

Since the right hand side does not depend on $\tau$, it makes no difference
whether we take $\sup_{0\leq\tau\leq t}$ inside or outside of the expectation.

$\displaystyle D_{i}^{(0)}(t)\leq\tilde{D}_{i}^{(0)}(t)\leq$
$\displaystyle\sum_{s\in\mathcal{S}}\sum_{\begin{subarray}{c}s^{\prime}\in\mathcal{S}\\\
s^{\prime}\neq
s\end{subarray}}\mathbb{E}\left[\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(t)\oplus\mathcal{K}_{i,ss^{\prime}}(t)\right)+\mathcal{N}_{i,s^{\prime}s}\left(\mathcal{H}_{i,s^{\prime}s}(t)\oplus\mathcal{K}_{i,s^{\prime}s}(t)\right)\right]$

The summations with respect to $s$ and $s^{\prime}$ only contribute a constant
factor $\left|\mathcal{S}\right|^{2}$ which will be neglected. Also, the same
bound applies for
$E\left[\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(t)\oplus\mathcal{K}_{i,ss^{\prime}}(t)\right)\right]$
and
$E\left[\mathcal{N}_{i,s^{\prime}s}\left(\mathcal{H}_{i,s^{\prime}s}(t)\oplus\mathcal{K}_{i,s^{\prime}s}(t)\right)\right]$,
so it is enough to keep track of only the first one, with a factor of $2$.

The rate functions are Lipschitz-continuous on a compact domain due to
assumption (4), so they are bounded; their maximum is denoted by $q_{\max}$.

$\displaystyle\mathbb{E}\left[\mathcal{N}_{i,ss^{\prime}}\left(\mathcal{H}_{i,ss^{\prime}}(t)\oplus\mathcal{K}_{i,ss^{\prime}}(t)\right)\right]=$
$\displaystyle\mathbb{E}\left[\int_{0}^{t}\left|q_{ss^{\prime}}\left(\phi_{i}(\tau)\right)\xi_{i,s^{\prime}}(\tau)-q_{ss^{\prime}}\left(\tilde{\phi}_{i}(\tau)\right)\hat{\xi}_{i,s^{\prime}}(\tau)\right|\mathrm{d}\tau\right]\leq$
$\displaystyle\mathbb{E}\left[\int_{0}^{t}q_{\max}\left|\xi_{i,s^{\prime}}(\tau)-\hat{\xi}_{i,s^{\prime}}(\tau)\right|+L_{q}\sum_{m=1}^{M}\sum_{\underline{r}\in\mathcal{S}^{m}}\left|\phi_{i,\underline{r}}^{(m)}(\tau)-\tilde{\phi}_{i,\underline{r}}^{(m)}(\tau)\right|\mathrm{d}\tau\right]\leq$
$\displaystyle\left(q_{\max}+L_{q}\right)\int_{0}^{t}\sum_{m=0}^{M}D_{i}^{(m)}(\tau)\mathrm{d}\tau\leq\left(q_{\max}+L_{q}\right)\int_{0}^{t}\sum_{m=0}^{M}\tilde{D}_{i}^{(m)}(\tau)\mathrm{d}\tau$

Setting
$\tilde{C}_{1}:=2\left(q_{\max}+L_{q}\right)\left|\mathcal{S}\right|^{2}$
yields

$\displaystyle
D_{i}^{(0)}(t)\leq\tilde{D}_{i}^{(0)}(t)\leq\tilde{C}_{1}\int_{0}^{t}\sum_{m=0}^{M}D_{i}^{(m)}(\tau)\mathrm{d}\tau\leq\tilde{C}_{1}\int_{0}^{t}\underbrace{\sum_{m=0}^{M}\tilde{D}_{i}^{(m)}(\tau)}_{=\tilde{D}_{i}(\tau)}\mathrm{d}\tau.$

∎

The second half of the proof of Theorem 2 involves estimating the difference
between the neighbors $\phi_{i}(t)$ and $\zeta_{i}(t)$ via the differences of
the indicators.

$\zeta_{i}(t)$ does not contain the indicators $\hat{\xi}_{i}(t)$ directly,
only their expectation $z_{i}(t)$. To bridge this gap, we introduce
“intermediate neighborhoods”

$\displaystyle\hat{\phi}_{i,\underline{s}}^{(m)}(t)=$
$\displaystyle\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}\hat{\xi}_{\underline{j},\underline{s}}^{(m)}(t).$

Note that under (16) and independent initial conditions,

$\mathbb{E}\left(\hat{\xi}_{\underline{i},\underline{s}}^{(m)}\right)=\mathbb{E}\left(\prod_{l=1}^{m}\hat{\xi}_{i_{l},s_{l}}(t)\right)=\prod_{l=1}^{m}\mathbb{E}\left(\hat{\xi}_{i_{l},s_{l}}(t)\right)=\prod_{l=1}^{m}z_{i_{l},s_{l}}(t)=z_{\underline{i},\underline{s}}^{(m)}$

for non-secondary loop $\underline{i}$ indices. Assumption (5) was made to
ensure secondary loops have low total weight.

$\displaystyle\begin{split}&\left|\mathbb{E}\left(\hat{\phi}_{i,\underline{s}}^{(m)}(t)\right)-\zeta_{i,\underline{s}}^{(m)}(t)\right|=\left|\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}\left[\mathbb{E}\left(\hat{\xi}_{\underline{j},\underline{s}}^{(m)}(t)\right)-z_{\underline{j},\underline{s}}^{(m)}(t)\right]\right|=\\\
&\left|\sum_{\begin{subarray}{c}\underline{j}\in[N]^{m}\\\
\underline{j}\textrm{ s.
loop}\end{subarray}}w_{i,\underline{j}}^{(m)}\left[\mathbb{E}\left(\hat{\xi}_{\underline{j},\underline{s}}^{(m)}(t)\right)-z_{\underline{j},\underline{s}}^{(m)}(t)\right]\right|\leq\sum_{\begin{subarray}{c}\underline{j}\in[N]^{m}\\\
\underline{j}\textrm{ s. loop}\end{subarray}}w_{i,\underline{j}}^{(m)}\leq
R\sqrt{w_{\max}}.\end{split}$ (55)

The next lemma shows that $\hat{\phi}_{i}(t)$ and $\zeta_{i}(t)$ are close.

######  Lemma 3.

Assume (16) holds with independent initial conditions. Then there is a
$\tilde{C}_{2}=\tilde{C}_{2}\left(\delta_{\max},R\right)$ such that for any
$1\leq m\leq M,\ i\in[N]$

$\displaystyle\sup_{0\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\hat{\phi}_{i,\underline{s}}^{(m)}(t)-\zeta_{i,\underline{s}}^{(m)}(t)\right|\right]\leq\tilde{C}_{2}\sqrt{w_{\max}}.$
(56)

If we further assume $M=1$, there exists a $\tilde{C}_{3}$ such that for all
$t\geq 0$,

$\displaystyle\mathbb{E}\left[\sup_{0\leq
t}\sum_{s\in\mathcal{S}}\left|\hat{\phi}_{i,s}(t)-\zeta_{i,s}(t)\right|\right]\leq\tilde{C}_{3}\underbrace{\sqrt{\sum_{j=1}^{n}w_{ij}^{2}}}_{=\mu_{i}}.$
(57)

###### Proof.

(Lemma 3)

We start by applying (55).

$\displaystyle\sup_{0\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\hat{\phi}_{i,\underline{s}}^{(m)}(t)-\zeta_{i,\underline{s}}^{(m)}(t)\right|\right]\leq$
$\displaystyle R\left|\mathcal{S}\right|^{M}\sqrt{w_{\max}}+\sup_{0\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\hat{\phi}_{i,\underline{s}}^{(m)}(t)-\mathbb{E}\left(\hat{\phi}_{i,\underline{s}}^{(m)}(t)\right)\right|\right].$

The first term is of the desired form; we examine the second term.

$\displaystyle\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\hat{\phi}_{i,\underline{s}}^{(m)}(t)-\mathbb{E}\left(\hat{\phi}_{i,\underline{s}}^{(m)}(t)\right)\right|\right]=\sum_{\underline{s}\in\mathcal{S}^{m}}\mathbb{E}\left(\left|\hat{\phi}_{i,\underline{s}}^{(m)}(t)-\mathbb{E}\left(\hat{\phi}_{i,\underline{s}}^{(m)}(t)\right)\right|\right)\leq$
$\displaystyle\sum_{\underline{s}\in\mathcal{S}^{m}}\sqrt{\mathbb{D}^{2}\left(\hat{\phi}_{i,\underline{s}}^{(m)}(t)\right)}=\sum_{\underline{s}\in\mathcal{S}^{m}}\sqrt{\sum_{\underline{j}\in[N]^{m}}\left(w_{i,\underline{j}}^{(m)}\right)^{2}\mathbb{D}^{2}\left(\hat{\xi}_{\underline{j},\underline{s}}^{(m)}(t)\right)}\leq$
$\displaystyle\left|\mathcal{S}\right|^{M}\sqrt{\sum_{\underline{j}\in[N]^{m}}\left(w_{i,\underline{j}}^{(m)}\right)^{2}}\leq\left|\mathcal{S}\right|^{M}\sqrt{\delta_{\max}w_{\max}}.$

The bound is uniform in $t$, so it can be upgraded to $\sup_{0\leq t}$ for
free, and (56) holds with
$\tilde{C}_{2}=\left(R+\sqrt{\delta_{\max}}\right)\left|\mathcal{S}\right|^{M}.$

Next we turn to (57).

$\hat{\xi}_{i,s}(t)-z_{i,s}(t)$ is a martingale, so

$\hat{\phi}_{i,s}(t)-\zeta_{i,s}(t)=\sum_{j=1}^{N}w_{ij}\left[\hat{\xi}_{j,s}(t)-z_{j,s}(t)\right]$

is also a martingale, and Doob’s martingale inequality yields

$\displaystyle\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\sum_{s\in\mathcal{S}}\left|\hat{\phi}_{i,s}(\tau)-\zeta_{i,s}(\tau)\right|\right]\leq\sum_{s\in\mathcal{S}}\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\left|\hat{\phi}_{i,s}(\tau)-\zeta_{i,s}(\tau)\right|\right]\leq$
$\displaystyle\sum_{s\in\mathcal{S}}\sqrt{\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\left|\hat{\phi}_{i,s}(\tau)-\zeta_{i,s}(\tau)\right|^{2}\right]}\leq
2\sum_{s\in\mathcal{S}}\sqrt{\mathbb{E}\left(\left|\hat{\phi}_{i,s}(t)-\zeta_{i,s}(t)\right|^{2}\right)}=$
$\displaystyle
2\sum_{s\in\mathcal{S}}\sqrt{\mathbb{D}^{2}\left(\hat{\phi}_{i,s}(t)\right)}=2\sum_{s\in\mathcal{S}}\sqrt{\sum_{j=1}^{N}w_{ij}^{2}\mathbb{D}^{2}\left(\hat{\xi}_{j,s}(t)\right)}\leq\underbrace{2\left|\mathcal{S}\right|}_{=:\tilde{C}_{3}}\sqrt{\sum_{j=1}^{N}w_{ij}^{2}}.$

∎

Next we show an upper bound for the differences of neighborhood vectors, which
are captured by the values $D^{(m)}_{\max}(t)$.

######  Lemma 4.

Assume (16) and independent initial conditions. Then there exist constants
$\tilde{C}_{4}=\tilde{C}_{5}\left(\delta_{\max}\right)$ such that for any
$t\geq 0$ and $1\leq m\leq M$

$\displaystyle
D^{(m)}_{\max}(t)\leq\tilde{C}_{2}\sqrt{w_{\max}}+\tilde{C}_{4}D_{\max}^{(0)}(t).$

where $\tilde{C}_{2}$ comes from Lemma 3.

If we further assume $M=1$ then

$\displaystyle\tilde{D}^{(1)}(t)\leq\tilde{C}_{3}\mu+W\tilde{D}^{(0)}(t).$

where $\tilde{C}_{3}$ comes from Lemma 3.

###### Proof.

(Lemma 4)

Using Lemma 3, we have

$\displaystyle D_{i}^{(m)}(t)=\sup_{0\leq\tau\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\phi_{i,\underline{s}}^{(m)}(\tau)-\zeta_{i,\underline{s}}^{(m)}(\tau)\right|\right]\leq$
$\displaystyle\tilde{C}_{2}\sqrt{w_{\max}}+\sup_{0\leq\tau\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\phi_{i,\underline{s}}^{(m)}(\tau)-\hat{\phi}_{i,\underline{s}}^{(m)}(\tau)\right|\right]\leq$
$\displaystyle\tilde{C}_{2}\sqrt{w_{\max}}+\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}\left(\sup_{0\leq\tau\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\xi_{\underline{j},\underline{s}}^{(m)}(\tau)-\hat{\xi}_{\underline{j},\underline{s}}^{(m)}(\tau)\right|\right]\right).$

Lemma 1 provides

$\displaystyle\left|\xi_{\underline{j},\underline{s}}^{(m)}(\tau)-\hat{\xi}_{\underline{j},\underline{s}}^{(m)}(\tau)\right|\leq\sum_{l=1}^{m}\left|\xi_{j_{l},s_{l}}(\tau)-\hat{\xi}_{j_{l},s_{l}}(\tau)\right|$
$\displaystyle\sup_{0\leq\tau\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\left|\xi_{\underline{j},\underline{s}}^{(m)}(\tau)-\hat{\xi}_{\underline{j},\underline{s}}^{(m)}(\tau)\right|\right]\leq\sup_{0\leq\tau\leq
t}\mathbb{E}\left[\sum_{\underline{s}\in\mathcal{S}^{m}}\sum_{l=1}^{m}\left|\xi_{j_{l},s_{l}}(\tau)-\hat{\xi}_{j_{l},s_{l}}(\tau)\right|\right]\leq$
$\displaystyle\left|S\right|^{M}\sum_{l=1}^{m}\sup_{0\leq\tau\leq
t}\mathbb{E}\left[\sum_{r\in\mathcal{S}}\left|\xi_{j_{l},r}(\tau)-\hat{\xi}_{j_{l},r}(\tau)\right|\right]\leq\left|S\right|^{M}\sum_{l=1}^{m}D_{j_{l}}^{(0)}(t)\leq
M\left|\mathcal{S}\right|^{M}D_{\max}^{(0)}(t).$

Putting the inequalities together yields

$\displaystyle D_{i}^{(m)}(t)\leq$
$\displaystyle\tilde{C}_{2}\sqrt{w_{\max}}+M\left|\mathcal{S}\right|^{M}D_{\max}^{(0)}(t)\underbrace{\sum_{\underline{j}\in[N]^{m}}w_{i,\underline{j}}^{(m)}}_{=\delta^{(m)}(i)}$
$\displaystyle D_{\max}^{(m)}(t)\leq$
$\displaystyle\tilde{C}_{2}\sqrt{w_{\max}}+\underbrace{M\left|\mathcal{S}\right|^{M}\delta_{\max}}_{=:\tilde{C}_{4}}D_{\max}^{(0)}(t).$

For the second part of Lemma 4, we once again use Lemma 3.

$\displaystyle\tilde{D}_{i}^{(1)}(t)=\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\sum_{s\in\mathcal{S}}\left|\phi_{i,s}(\tau)-\zeta_{i,s}(\tau)\right|\right]\leq$
$\displaystyle\tilde{C}_{3}\mu_{i}+\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\sum_{s\in\mathcal{S}}\left|\phi_{i,s}(\tau)-\hat{\phi}_{i,s}(\tau)\right|\right]\leq$
$\displaystyle\tilde{C}_{3}\mu_{i}+\sum_{j=1}^{N}w_{ij}\left(\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\sum_{s\in\mathcal{S}}\left|\xi_{j,s}(\tau)-\hat{\xi}_{j,s}(\tau)\right|\right]\right)=\tilde{C}_{3}\mu_{i}+\sum_{j=1}^{N}w_{ij}\tilde{D}_{j}^{(0)}(t),$

so

$\displaystyle\tilde{D}^{(1)}(t)\leq\tilde{C}_{3}\mu+W\tilde{D}^{(0)}(t).$

∎

With all the preparations done, we finally turn to proving Theorem 2.

###### Proof.

(Theorem 2)

Using Lemma 2 and 4 and Grönwall’s inequality yields

$\displaystyle
D_{\max}(t)=D_{\max}^{0}(t)+\sum_{m=1}^{M}D_{\max}^{(m)}(t)\leq$
$\displaystyle
M\tilde{C}_{2}\sqrt{w_{\max}}+\left(M\tilde{C}_{4}+1\right)D_{\max}^{(0)}(t)\leq$
$\displaystyle
M\tilde{C}_{2}\sqrt{w_{\max}}+\left(M\tilde{C}_{4}+1\right)\int_{0}^{t}D_{\max}(\tau)\mathrm{d}\tau,$

so

$\displaystyle
D_{\max}(t)\leq\underbrace{M\tilde{C}_{2}e^{\left(M\tilde{C}_{4}+1\right)t}}_{=:C}\sqrt{w_{\max}}.$

Proving the second part is similar.

$\displaystyle\tilde{D}(t)=$
$\displaystyle\tilde{D}^{(0)}(t)+\sum_{m=1}^{M}\tilde{D}^{(m)}(t)\leq\underbrace{M\tilde{C}_{3}}_{=:C_{1}}\mu+M\left(W+I\right)\tilde{D}^{(0)}(t)\leq$
$\displaystyle
C_{1}\mu+\underbrace{\tilde{C_{1}}M}_{=:C_{2}}\int_{0}^{t}\left(W+I\right)\tilde{D}(\tau)\mathrm{d}\tau\Rightarrow$
$\displaystyle\left\|\tilde{D}(t)\right\|\leq$ $\displaystyle
C_{1}\left\|\mu\right\|+C_{2}\left\|W+I\right\|\int_{0}^{t}\left\|\tilde{D}(\tau)\right\|\mathrm{d}\tau,$

so

$\displaystyle\left\|\tilde{D}(t)\right\|\leq$ $\displaystyle
C_{1}e^{C_{2}\left\|W+I\right\|t}\left\|\mu\right\|.$

∎

###  6.3 Proof of Theorems 4 and 5

###### Proof.

(Theorem 4) For a fixed $t$ and $s$, we apply Doob’s inequality for the
martingale $\frac{1}{K}\sum_{i=1}^{K}(\hat{\xi}_{i,s}(t)-z_{i,s}(t))$ and use
independence to get

$\displaystyle\mathbb{E}\left(\sup_{0\leq\tau\leq
t}\left|\frac{1}{K}\sum_{i=1}^{K}\left(z_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right)\right|\right)\leq
2\mathbb{D}\left(\frac{1}{K}\sum_{i=1}^{K}\left(z_{i,s}(t)-\hat{\xi}_{i,s}(t)\right)\right)=$
$\displaystyle
2\left(\frac{1}{K^{2}}\sum_{i=1}^{K}\underbrace{\mathbb{D}^{2}\left(z_{i,s}(t)-\hat{\xi}_{i,s}(t)\right)}_{\leq
1}\right)^{1/2}\leq\frac{2}{\sqrt{K}},$ (58)

and (24) follows by inserting $\sum_{s\in\mathcal{S}}$ on the left hand side
at the cost of an $|\mathcal{S}|$ factor on the right hand side. The bound is
uniform in $t$, so we can upgrade to $\sup_{0\leq t}$. ∎

###### Proof.

(Theorem 5) For (25), we consider $0\leq\tau\leq t$ and use both Theorems 2
and 4:

$\displaystyle\mathbb{E}\left(\sum_{s\in\mathcal{S}}\left|\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i,s}(\tau)-z_{i,s}(\tau)\right)\right|\right)\leq$
$\displaystyle\quad\sum_{s\in\mathcal{S}}\mathbb{E}\left(\left|\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right)\right|+\left|\frac{1}{N}\sum_{i=1}^{N}\left(\hat{\xi}_{i,s}(\tau)-z_{i,s}(\tau)\right)\right|\right)\leq$
$\displaystyle\quad\frac{1}{N}\sum_{i=1}^{N}\underbrace{\sum_{s\in\mathcal{S}}\mathbb{E}\left|\left(\xi_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right)\right|}_{\leq
D_{\max}(t)}+\sum_{s\in\mathcal{S}}\underbrace{\mathbb{E}\left|\frac{1}{N}\sum_{i=1}^{N}\left(\hat{\xi}_{i,s}(\tau)-z_{i,s}(\tau)\right)\right|}_{\leq
2/\sqrt{N}}\leq$ $\displaystyle\quad
D_{\max}(t)+\frac{2|\mathcal{S}|}{\sqrt{N}}\leq
C\left(\sqrt{w_{\max}}+\frac{1}{\sqrt{N}}\right).$

The derivation of (26) is analogous to (25) with the exception of keeping the
$\sup_{0\leq\tau\leq t}$ inside the expectation and using (18) instead of
(17).

For (27), we just note that

$\displaystyle\mathbb{E}\left[\sup_{0\leq\tau\leq
t}\left(\sum_{s\in\mathcal{S}}\left|\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i,s}(\tau)-\hat{\xi}_{i,s}(\tau)\right)\right|\right)\right]$
$\displaystyle\leq\frac{1}{N}\|\tilde{D}(t)\|_{1}$
$\displaystyle\leq\frac{1}{\sqrt{N}}\|\tilde{D}(t)\|_{2}=O\left(\sqrt{\frac{1}{N}\|\mu\|_{2}^{2}}\right),$

and the rest of the argument is essentially identical to the previous one.

∎

###  6.4 Proof of Proposition 2

Let $p_{k^{m}}^{(m)}:=\frac{N_{k^{(m)}}}{N}$ denote the ratio of vertices in
the local group $k^{m}.$

$\displaystyle\bar{\zeta}_{k}^{(m)}(t)$
$\displaystyle=\sum_{\underline{l}^{(m)}}\bar{w}_{k^{(m)},\underline{j}^{(m)}}^{(m)}\prod_{r=1}^{m}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)$
$\displaystyle=\sum_{\underline{l}^{(m)}}\left(\prod_{r=1}^{m}p_{l_{r}}^{(m)}\right)\left(a_{k^{m}}^{(m)}+\sum_{r=1}^{m}a_{l_{r}^{(m)}}^{(m)}\right)\prod_{r=1}^{m}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)$
$\displaystyle=\sum_{\underline{l}^{(m)}}\left(a_{k^{m}}^{(m)}+\sum_{r=1}^{m}a_{l_{r}^{(m)}}^{(m)}\right)\prod_{r=1}^{m}p_{l_{r}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)$
(59)

Observe

$\displaystyle\sum_{l^{m}=1}^{K^{(m)}}p_{l_{r}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|a^{(m)}(\iota)=l_{r}^{(m)}\right)=$
$\displaystyle\mathbb{E}\left(\mathbb{E}\left(\left.z_{\iota}(t)\right|a^{(m)}(\iota)=l_{r}^{(m)}\right)\right)=\mathbb{E}\left(z_{\iota}(t)\right).$

Also introduce

$\displaystyle\psi^{(m)}(t):=\sum_{l=1}^{K^{(m)}}a_{l}^{(m)}p_{l}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|a^{(m)}(\iota)=l^{(m)}\right)$

which is renaissance of an activity biased average.

We expand (59) based on the terms
$a_{k^{m}}^{(m)}+\sum_{r=1}^{m}a_{l_{r}^{(m)}}^{(m)}$. For $a_{k^{m}}^{(m)}$

$\displaystyle
a_{k^{m}}^{(m)}\sum_{\underline{l}^{(m)}}\prod_{r=1}^{m}p_{l_{r}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)=$
$\displaystyle
a_{k^{m}}^{(m)}\left(\sum_{l=1}^{K^{(}m)}p_{l}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l^{(m)}\right)\right)^{m}=$
$\displaystyle a_{k^{m}}^{(m)}\mathbb{E}^{m}\left(z_{\iota}(t)\right).$

For the $a_{l_{r}^{\prime}}^{(m)}$ terms we have

$\displaystyle\sum_{\underline{l}^{(m)}}a_{l_{r^{\prime}}^{(m)}}\prod_{r=1}^{m}p_{l_{r}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)=$
$\displaystyle\underbrace{\sum_{l_{r^{\prime}}=1}^{K^{(m)}}a_{l_{r^{\prime}}^{(m)}}p_{l_{r^{\prime}}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r^{\prime}}^{(m)}\right)}_{\psi^{(m)}(t)}\sum_{\begin{subarray}{c}l_{r}^{(m)}=1\\\
r\neq r^{\prime}\end{subarray}}^{K^{(m)}}\prod_{\begin{subarray}{c}r=1\\\
r\neq
r^{\prime}\end{subarray}}^{m}p_{l_{r}}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l_{r}^{(m)}\right)=$
$\displaystyle\psi^{(m)}(t)\left(\sum_{l=1}^{K^{(}m)}p_{l}^{(m)}\mathbb{E}\left(\left.z_{\iota}(t)\right|k^{(m)}(\iota)=l^{(m)}\right)\right)^{m-1}=\psi^{(m)}(t)\mathbb{E}^{m-1}\left(z_{\iota}(t)\right).$

Therefore, (59) reduces to

$\displaystyle\bar{\zeta}^{(m)}_{k}(t)$
$\displaystyle=a_{k^{m}}^{(m)}\mathbb{E}^{m}\left(z_{\iota}(t)\right)+\psi^{(m)}(t)\mathbb{E}^{m-1}\left(z_{\iota}(t)\right)$
$\displaystyle=\left(a_{k^{m}}^{(m)}\mathbb{E}\left(z_{\iota}(t)\right)+\psi^{(m)}(t)\right)\mathbb{E}^{m-1}\left(z_{\iota}(t)\right).$

###  6.5 Proof of Theorem 6

Recall (47). We call the sets $X,Y\subset[N]$ $\varepsilon$-regular if for all
$A\subseteq X,\ B\subseteq Y$ such that
$\left|A\right|>\varepsilon\left|X\right|,\
\left|B\right|>\varepsilon\left|Y\right|$ one has

$\displaystyle\left|\rho\left(A,B\right)-\rho\left(X,Y\right)\right|<\varepsilon.$

We use Szemerédi’s regularity lemma.

######  Lemma.

(Szemerédi’s regularity lemma)

For every $\varepsilon>0,\ K_{\min}\in\mathbb{Z}^{+}$ there is a $K_{\max}$
such that if $N\geq K_{\max}$ there is a partition $V_{0},V_{1},\dots,V_{K}$
such that

$\displaystyle\left|V_{0}\right|<\varepsilon N,$
$\displaystyle\left|V_{1}\right|=\dots=\left|V_{K}\right|,$ $\displaystyle
K_{\min}\leq K\leq K_{\max}$

and there are at most $\varepsilon{K\choose 2}$ pairs of
$\left(V_{k},V_{l}\right),\ 1\leq k<l\leq K$ such that they are not
$\varepsilon$-regular.

Fix a $\varepsilon^{\prime}>0$ and a $K_{\min}$ such that

$\displaystyle K_{\min}>\frac{1}{\varepsilon^{\prime}}.$

This choice ensures that there are enough boxes such that most of the vertices
are between boxes and not within them. This is a fairly common approach in the
context of Szemerédi’s regularity lemma [32].

Using Szemerédi’s regularity lemma for $\varepsilon^{\prime}$, we obtain a
partition denoted by $V_{0},V_{1},\dots,V_{K}.$

For $p$ and $\kappa$, as defined in (48) and (49), we have the following
inequalities:

$\displaystyle
p=\frac{\bar{d}}{N}\geq\frac{p_{0}(N-1)}{N}\geq\frac{p_{0}}{2}>0$
$\displaystyle
1=\sum_{k=0}^{K}\frac{\left|V_{k}\right|}{N}\geq\sum_{k=1}^{K}\frac{\left|V_{k}\right|}{N}=K\kappa\Longrightarrow\kappa\leq\frac{1}{K}\leq\frac{1}{K_{\min}}<\varepsilon^{\prime}$

where we used $N\geq 2$.

Introduce the notations

$\displaystyle\bar{z}_{k}(t):=$
$\displaystyle\frac{1}{\left|V_{k}\right|}\sum_{i\in V_{k}}z_{i}(t),$
$\displaystyle\psi(t):=$
$\displaystyle\sum_{k=1}^{K}\frac{\left|V_{k}\right|}{N}\left\|\bar{z}_{k}(t)-v_{k}(t)\right\|_{1}=\kappa\sum_{k=1}^{K}\left\|\bar{z}_{k}(t)-v_{k}(t)\right\|_{1}.$

If $V_{0}=\emptyset,$ we use the convention $z_{0}(t)\equiv 0.$

From (53) and (52), we have

$\displaystyle\bar{z}(t)=\frac{1}{N}\sum_{i=1}^{N}z_{i}(t)=\sum_{k=0}^{K}\frac{\left|V_{k}\right|}{N}\frac{1}{\left|V_{k}\right|}\sum_{i\in
V_{k}}z_{i}(t)=\sum_{k=0}^{K}\frac{\left|V_{k}\right|}{N}\bar{z}_{k}(t)$
$\displaystyle\left\|\bar{z}(t)-\bar{v}(t)\right\|_{1}=\left\|\frac{\left|V_{0}\right|}{N}\bar{z}_{0}(t)+\sum_{k=1}^{K}\frac{\left|V_{k}\right|}{N}\left[\bar{z}_{k}(t)-v_{k}(t)\right]\right\|_{1}\leq$
$\displaystyle\frac{\left|V_{0}\right|}{N}\left\|\bar{z}_{0}(t)\right\|_{1}+\sum_{k=1}^{K}\frac{\left|V_{k}\right|}{N}\left\|\bar{z}_{k}(t)-v_{k}(t)\right\|_{1}\leq\varepsilon^{\prime}+\psi(t)$

where in the last step we used $\left|V_{0}\right|<\varepsilon^{\prime}N$ and

$\displaystyle\left\|\bar{z}_{0}(t)\right\|_{1}\leq\frac{1}{\left|V_{0}\right|}\sum_{i\in
V_{0}}\underbrace{\left\|z_{i}(t)\right\|_{1}}_{=1}=1.$

Going forward, it is enough to examine $\psi(t).$

Next we calculate the derivative of $\bar{z}_{k}(t).$ As $M=1,$ (29) takes the
form

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{i,s}(t)=$
$\displaystyle\sum_{s^{\prime}\in\mathcal{S}}q_{ss^{\prime}}\left(\zeta_{i}(t)\right)z_{i,s^{\prime}}(t)=$
$\displaystyle\sum_{s^{\prime}\in\mathcal{S}}q_{ss^{\prime}}^{(0)}z_{i,s^{\prime}}(t)+\sum_{s^{\prime}\in\mathcal{S}}\sum_{r\in\mathcal{S}}q_{ss^{\prime},r}^{(1)}\zeta_{i,r}(t)z_{i,s^{\prime}}(t)=$
$\displaystyle\sum_{s^{\prime}\in\mathcal{S}}q_{ss^{\prime}}^{(0)}z_{i,s^{\prime}}(t)+\sum_{s^{\prime}\in\mathcal{S}}\sum_{r\in\mathcal{S}}q_{ss^{\prime},r}^{(1)}\left[\sum_{j=1}^{N}\underbrace{\frac{a_{ij}}{\bar{d}}}_{w_{ij}^{(m)}}z_{i,s^{\prime}}(t)z_{j,r}(t)\right]$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bar{z}_{k,s}(t)=$
$\displaystyle\sum_{s^{\prime}\in\mathcal{S}}q_{ss^{\prime}}^{(0)}\bar{z}_{k,s^{\prime}}(t)+\sum_{s^{\prime}\in\mathcal{S}}\sum_{r\in\mathcal{S}}q_{ss^{\prime},r}^{(1)}\left[\frac{1}{\left|V_{k}\right|}\sum_{i\in
V_{k}}\sum_{j=1}^{N}\frac{a_{ij}}{\bar{d}}z_{i,s^{\prime}}(t)z_{j,r}(t)\right]$

Similarly,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}v_{k,s}(t)=$
$\displaystyle\underbrace{\sum_{s^{\prime}\in\mathcal{S}}q_{ss^{\prime}}^{(0)}v_{k,s^{\prime}}(t)+\sum_{s^{\prime}\in\mathcal{S}}\sum_{r\in\mathcal{S}}q_{ss^{\prime},r}^{(1)}\sum_{l=1}^{K}{\bar{w}_{kl}}v_{k,s^{\prime}}(t)v_{l,r}(t)}_{=:f_{k,s}\left(V(t)\right)}$

where $V(t):=\left(v_{k,s}(t)\right)_{k\in[K],\ s\in\mathcal{S}}$ and
$\overline{Z}(t)=\left(\bar{z}_{k,s}(t)\right)_{k\in[K],\ s\in\mathcal{S}}$
analogously.

Next we show a Lipschitz-type inequality for
$f_{k}=\left(f_{k,s}\right)_{s\in\mathcal{S}}.$

$\displaystyle\left|\bar{z}_{k,s^{\prime}}(t)\sum_{l=1}^{K}\bar{w}_{kl}\bar{z}_{l,r}(t)-v_{k,s^{\prime}}(t)\sum_{l=1}^{K}\bar{w}_{kl}v_{l,r}(t)\right|\leq$
$\displaystyle\left|\bar{z}_{k,s^{\prime}}(t)-v_{k,s^{\prime}}(t)\right|\underbrace{\sum_{l=1}^{K}\bar{w}_{kl}\bar{z}_{l,r}(t)}_{\leq\sum_{k=1}^{K}\bar{w}_{kl}\leq\frac{2}{p_{0}K}K=\frac{2}{p_{0}}}+\underbrace{v_{k,s^{\prime}}(t)}_{\leq
1}\sum_{l=1}^{K}\underbrace{\bar{w}_{kl}}_{\leq\frac{2\kappa}{p_{0}}}\left|\bar{z}_{l,r}(t)-v_{l,r}(t)\right|\leq$
$\displaystyle\frac{2}{p_{0}}\left(\left|\bar{z}_{k,s^{\prime}}(t)-v_{k,s^{\prime}}(t)\right|+\kappa\sum_{l=1}^{K}\left|\bar{z}_{l,r}(t)-v_{l,r}(t)\right|\right),$

so

$\displaystyle\left|f_{k,s}\left(\bar{Z}(t)\right)-f_{k,s}\left(V(t)\right)\right|\leq
q_{\max}\sum_{s^{\prime}\in\mathcal{S}}\left|\bar{z}_{k,s^{\prime}}(t)-v_{k,s^{\prime}}(t)\right|+$
$\displaystyle\frac{2q_{\max}}{p_{0}}\sum_{s^{\prime}\in\mathcal{S}}\sum_{r\in\mathcal{S}}\left(\left|\bar{z}_{k,s^{\prime}}(t)-v_{k,s^{\prime}}(t)\right|+\kappa\sum_{l=1}^{K}\left|\bar{z}_{l,r}(t)-v_{l,r}(t)\right|\right)=$
$\displaystyle
q_{\max}\left(1+\frac{2\left|\mathcal{S}\right|}{p_{0}}\right)\left\|\bar{z}_{k}(t)-v_{k}(t)\right\|_{1}+\frac{2q_{\max}\left|\mathcal{S}\right|}{p_{0}}\psi(t).$

Summation for $s\in\mathcal{S}$ results only in an extra $\mathcal{S}$ factor,
so there exists a constant $L_{f}$ such that

$\displaystyle\left\|f_{k}\left(\overline{Z}(t)\right)-f_{k}\left(V(t)\right)\right\|_{1}\leq
L_{f}\left(\left\|\bar{z}_{k}(t)-v_{k}(t)\right\|_{1}+\psi(t)\right).$ (60)

Next we look to replace the right hand side of
$\frac{\mathrm{d}}{\mathrm{d}t}\bar{z}_{k,s}(t)$ with
$f_{k,s}\left(\overline{Z}(t)\right).$ The corresponding error term is

$\displaystyle
g_{k,s}(t):=\sum_{s^{\prime}\in\mathcal{S}}\sum_{r\in\mathcal{S}}q_{ss^{\prime},r}^{(1)}\left[\frac{1}{\left|V_{k}\right|}\sum_{i\in
V_{k}}\sum_{j=1}^{N}\frac{a_{ij}}{\bar{d}}z_{i,s^{\prime}}(t)z_{j,r}(t)-\sum_{l=1}^{K}{\bar{w}_{kl}}\bar{z}_{k,s^{\prime}}(t)\bar{z}_{l,r}(t)\right],$
(61)

and from
$\frac{\mathrm{d}}{\mathrm{d}t}\bar{z}_{k}(t)=g_{k}(t)+f_{k}\left(\overline{Z}(t)\right)$,
we have

$\displaystyle\bar{z}_{k}(t)=$
$\displaystyle\bar{z}_{k}(0)+\int_{0}^{t}g_{k}(\tau)\mathrm{d}\tau+\int_{0}^{t}f_{k}\left(\overline{Z}(\tau)\right)\mathrm{d}\tau.$

Using $\bar{z}_{k}(0)=v_{k}(0)$, $\psi(t)$ can be bounded from above by

$\displaystyle\psi(t)=$
$\displaystyle\kappa\sum_{k=1}^{K}\left\|\bar{z}_{k}(t)-v_{k}(t)\right\|_{1}\leq$
$\displaystyle t\cdot\sup_{0\leq\tau\leq
t}\kappa\sum_{k=1}^{K}\left\|g_{k}(\tau)\right\|_{1}+\int_{0}^{t}\kappa\sum_{k=1}^{K}\left\|f_{k}\left(\overline{Z}(\tau)\right)-f_{k}\left(V(\tau)\right)\right\|_{1}\mathrm{d}\tau\leq$
$\displaystyle t\cdot\sup_{0\leq\tau\leq
t}\kappa\sum_{k=1}^{K}\left\|g_{k}(\tau)\right\|_{1}+L_{f}\int_{0}^{t}\kappa\sum_{k=1}^{K}\left(\left\|\bar{z_{k}}(\tau)-v_{k}(\tau)\right\|_{1}+\psi(\tau)\right)\mathrm{d}\tau\leq$
$\displaystyle t\cdot\sup_{0\leq\tau\leq
t}\kappa\sum_{k=1}^{K}\left\|g_{k}(\tau)\right\|_{1}+2L_{f}\int_{0}^{t}\psi(\tau)\mathrm{d}\tau,$

so from Grönwall’s inequality,

$\displaystyle\sup_{0\leq t\leq T}\psi(t)\leq$
$\displaystyle\left(T\cdot\sup_{0\leq t\leq
T}\kappa\sum_{k=1}^{K}\left\|g_{k}(t)\right\|_{1}\right)e^{2L_{f}T}.$

Therefore it is enough to show that $\sup_{0\leq t\leq
T}\kappa\sum_{k=1}^{K}\left\|g_{k}(t)\right\|_{1}=O\left(\varepsilon^{\prime}\right)$,
and with an appropriate choice of $\varepsilon=C\varepsilon^{\prime}$ we can
conclude

$\displaystyle\sup_{0\leq t\leq
T}\left\|\bar{z}(t)-\bar{v}(t)\right\|_{1}\leq\varepsilon.$

$\displaystyle\kappa\sum_{l=1}^{K}\left\|g_{k}(t)\right\|_{1}=$
$\displaystyle\kappa\sum_{s\in\mathcal{S}}\sum_{k=1}^{K}\left|\sum_{s^{\prime},r\in\mathcal{S}}q_{ss^{\prime},r}^{(1)}\left[\frac{1}{\left|V_{k}\right|}\sum_{i\in
V_{k}}\sum_{j=1}^{N}\frac{a_{ij}}{\bar{d}}z_{i,s^{\prime}}(t)z_{j,r}(t)-\sum_{l=1}^{K}{\bar{w}_{kl}}\bar{z}_{k,s^{\prime}}(t)\bar{z}_{l,r}(t)\right]\right|\leq$
$\displaystyle\kappa
q_{\max}\sum_{s,s^{\prime},r\in\mathcal{S}}\sum_{k=1}^{K}\sum_{l=0}^{K}\left|\frac{1}{|V_{k}|}\sum_{i\in
V_{k}}\sum_{j\in
V_{l}}\frac{a_{ij}}{\bar{d}}z_{i,s^{\prime}}(t)z_{j,r}(t)-\bar{w}_{kl}\bar{z}_{k,s^{\prime}}(t)\bar{z}_{l,r}(t)\right|$

$\sum_{s,s^{\prime},r\in\mathcal{S}}(\dots)$ only contributes a factor of
$\left|\mathcal{S}\right|^{3}$ which we can include in the constant factor
along with $q_{\max}.$ The remaining terms are

$\displaystyle\kappa\sum_{k=1}^{K}\sum_{l=0}^{K}\left|\frac{1}{|V_{k}|}\sum_{i\in
V_{k}}\sum_{j\in
V_{l}}\frac{a_{ij}}{\bar{d}}z_{i,s^{\prime}}(t)z_{j,r}(t)-\bar{w}_{kl}\bar{z}_{k,s^{\prime}}(t)\bar{z}_{l,r}(t)\right|.$
(62)

In the next step we shall get rid of the diagonal $(k,l)$ terms and also the
terms with $l=0$. We have

$\displaystyle\frac{1}{|V_{k}|}\sum_{i\in V_{k}}\sum_{j\in
V_{l}}\frac{a_{ij}}{\bar{d}}z_{i,s^{\prime}}(t)z_{j,r}(t)\leq\frac{1}{\left|V_{k}\right|\bar{d}}\sum_{i\in
V_{k}}\sum_{j\in
V_{k}}1=\frac{\left|V_{l}\right|}{\bar{d}}\leq\frac{\varepsilon^{\prime}}{p}\leq\frac{2\varepsilon^{\prime}}{p_{0}},$
$\displaystyle\bar{w}_{kl}\bar{z}_{k,s^{\prime}}(t)\bar{z}_{l,r}(t)\leq\frac{\kappa}{p}\leq\frac{2\varepsilon^{\prime}}{p_{0}},$

so each term in the sum of (62) is $O\left(\varepsilon^{\prime}\right).$ There
are $O(K)$ pairs which are either diagonal or $l=0$, so their overall
contribution to the sum is $O\left(\kappa
K\varepsilon^{\prime}\right)=O\left(\varepsilon^{\prime}\right),$ hence we can
neglect them and what we are left with is

$\displaystyle\kappa\
\sum_{(k,l)\in\mathcal{I}}\left|\frac{1}{|V_{k}|}\sum_{i\in V_{k}}\sum_{j\in
V_{l}}\frac{a_{ij}}{\bar{d}}z_{i,s^{\prime}}(t)z_{j,r}(t)-\bar{w}_{kl}\bar{z}_{k,s^{\prime}}(t)\bar{z}_{l,r}(t)\right|.$
(63)

where $\mathcal{I}=\\{(k,l)|k,l\in[K],k\neq l\\}$.

In order to have an upper bound for (63) we want to use the properties of the
$\varepsilon^{\prime}$-regular partition. However, Szemerédi’s regularity
lemma uses subsets of $[N]$, or in other words, $0-1$ valued indicators of
vertices compared to $z_{i,s}(t)$ which may take any value from $[0,1].$

To account for this problem, we introduce $N$ independent homogeneous Markov
processes taking values from $\mathcal{S}$. Each process makes Markov
transitions according to the transition rate matrix
$Q\left(\zeta_{i}(t)\right)$ and its initial distribution is given by
$\left(z_{i,s}(0)\right)_{s\in\mathcal{S}}.$ Let $\eta_{i,s}(t)$ be an
indicator of the $i$’th such process is at state $s$ at time $t$. We also
apply the notations

$\displaystyle\eta_{i}(t)=$
$\displaystyle\left(\eta_{i,s}(t)\right)_{s\in\mathcal{S}},$
$\displaystyle\bar{\eta}_{k}(t):=$
$\displaystyle\frac{1}{\left|V_{k}\right|}\sum_{i\in V_{k}}\eta_{i}(t).$

It is easy to see that $\mathbb{E}\left(\eta_{i}(t)\right)=z_{i}(t).$ Also,
since $i\in V_{k}$ and $j\in V_{l}$, $i$ and $j$ are different for $k\neq l$,
hence the corresponding processes are independent, so

$\displaystyle z_{i,s^{\prime}}(t)z_{j,r}(t)=$
$\displaystyle\mathbb{E}\left(\eta_{i,s^{\prime}}(t)\eta_{j,k}(t)\right),$
$\displaystyle\bar{z}_{k,s^{\prime}}(t)\bar{z}_{l,r}(t)=$
$\displaystyle\mathbb{E}\left(\bar{\eta}_{k,s^{\prime}}(t)\bar{\eta}_{l,r}(t)\right).$

Therefore, (63) can be bounded from above by

$\displaystyle\mathbb{E}\left[\kappa\
\sum_{(k,l)\in\mathcal{I}}\left|\frac{1}{|V_{k}|}\sum_{i\in V_{k}}\sum_{j\in
V_{l}}\frac{a_{ij}}{\bar{d}}\eta_{i,s^{\prime}}(t)\eta_{j,r}(t)-\bar{w}_{kl}\bar{\eta}_{k,s^{\prime}}(t)\bar{\eta}_{l,r}(t)\right|\right].$
(64)

The upper bound we aim to obtain does not depend on the artificial randomness
just introduced, hence the expectation is ignored.

We make some algebraic manipulation to end up with edge densities needed for
Szemerédi’s regularity lemma. We use the notation

$\displaystyle V_{k,s}(t):=\left\\{\left.i\in
V_{k}\right|\eta_{i,s}(t)=1\right\\}.$

Then

$\displaystyle\frac{1}{|V_{k}|}\sum_{i\in V_{k}}\sum_{j\in
V_{l}}\frac{a_{ij}}{\bar{d}}\eta_{i,s^{\prime}}(t)\eta_{j,r}(t)=\frac{1}{\left|V_{k}\right|\bar{d}}e\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)=$
$\displaystyle\frac{\left|V_{l}\right|}{\bar{d}}\rho\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)\frac{\left|V_{k,s^{\prime}}(t)\right|}{\left|V_{k}\right|}\frac{\left|V_{l,r}(t)\right|}{\left|V_{l}\right|}=\frac{\kappa}{p}\rho\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)\bar{\eta}_{k,s^{\prime}}(t)\bar{\eta}_{l,k}(t).$

By recalling (50), the inside of (64) can be rewritten as

$\displaystyle\frac{\kappa^{2}}{p}\sum_{(k,l)\in\mathcal{I}}\left|\rho\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)-\rho\left(V_{k},V_{l}\right)\right|\bar{\eta}_{k,s^{\prime}}(t)\bar{\eta}_{l,r}(t).$
(65)

Note that the summands of (65) are $O(1)$.

Using Szemerédi’s lemma to (65) is relatively straightforward from now on. We
still have to deal with non-$\varepsilon^{\prime}$-regular $k,l$ pairs, and
pairs where either
$\left|V_{k,s^{\prime}}(t)\right|\leq\varepsilon^{\prime}\left|V_{k}\right|$
or $\left|V_{l,r}(t)\right|\leq\varepsilon^{\prime}\left|V_{l}\right|$. The
former set of pairs are denoted by $\mathcal{I}_{1}$ and the latter by
$\mathcal{I}_{2}$, and
$\mathcal{I}_{3}:=\mathcal{I}\setminus\left(\mathcal{I}_{1}\cup\mathcal{I}_{2}\right)$
denotes the non-problematic pairs.

Then from $\left|\mathcal{I}_{1}\right|\leq\varepsilon^{\prime}{K\choose
2}\leq\varepsilon^{\prime}K^{2}$ we have

$\displaystyle\frac{\kappa^{2}}{p}\sum_{(k,l)\in\mathcal{I}_{1}}\left|\rho\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)-\rho\left(V_{k},V_{l}\right)\right|\bar{\eta}_{k,s^{\prime}}(t)\bar{\eta}_{l,r}(t)=O\left(\varepsilon^{\prime}\kappa^{2}K^{2}\right)=O\left(\varepsilon^{\prime}\right).$

$(k,l)\in\mathcal{I}_{2}$ is equivalent with
$\bar{\eta}_{k,s^{\prime}}(t)\leq\varepsilon^{\prime}$ or
$\bar{\eta}_{l,k}(t)\leq\varepsilon^{\prime}$, yielding

$\displaystyle\frac{\kappa^{2}}{p}\sum_{(k,l)\in\mathcal{I}_{2}}\left|\rho\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)-\rho\left(V_{k},V_{l}\right)\right|\bar{\eta}_{k,s^{\prime}}(t)\bar{\eta}_{l,r}(t)\leq$
$\displaystyle\frac{\varepsilon^{\prime}\kappa^{2}}{p}\sum_{(k,l)\in\mathcal{I}_{2}}1=O\left(\varepsilon^{\prime}\kappa^{2}K^{2}\right)=O\left(\varepsilon^{\prime}\right).$

Finally, $(k,l)\in\mathcal{I}_{3}$ gives

$\displaystyle\left|\rho\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)-\rho\left(V_{k},V_{l}\right)\right|<\varepsilon^{\prime}\Rightarrow$
$\displaystyle\frac{\kappa^{2}}{p}\sum_{(k,l)\in\mathcal{I}_{3}}\left|\rho\left(V_{k,s^{\prime}}(t),V_{l,r}(t)\right)-\rho\left(V_{k},V_{l}\right)\right|\bar{\eta}_{k,s^{\prime}}(t)\bar{\eta}_{l,r}(t)\leq$
$\displaystyle\frac{\varepsilon^{\prime}\kappa^{2}}{p}\sum_{(k,l)\in\mathcal{I}_{2}}1=O\left(\varepsilon^{\prime}\kappa^{2}K^{2}\right)=O\left(\varepsilon^{\prime}\right).$

This ensures that $\sup_{0\leq t\leq
T}\kappa\sum_{k=1}^{K}\left\|g_{k}(t)\right\|_{1}=O\left(\varepsilon^{\prime}\right)$
indeed holds, concluding the proof of Theorem 6.

∎

## References

  * [1] N. Alon, R. A. Duke, H. Lefmann, V. Rödl, and R. Yuster.  The algorithmic aspects of the regularity lemma.  Proceedings, 33rd Annual Symposium on Foundations of Computer Science, pages 473–481, 1992.
  * [2] R. Bakhshi, L. Cloth, W. Fokkink, and B. R. Haverkort.  Mean-field framework for performance evaluation of push-pull gossip protocols.  Performance Evaluation, 68(2):157–179, Feb. 2011.
  * [3] Á. Bodó, G. Y. Katona, and P. L. Simon.  SIS epidemic propagation on hypergraphs.  Bulletin of Mathematical Biology, 78:713–735, 2016.
  * [4] S. Bonaccorsi, S. Ottaviano, D. Mugnolo, and F. De Pellegrini.  Epidemic outbreaks in networks with equitable or almost-equitable partitions.  SIAM Journal on Applied Mathematics, 75:2421–2443, 11 2015.
  * [5] D. Bruneo, M. Scarpa, A. Bobbio, D. Cerotti, and M. Gribaudo.  Markovian agent modeling swarm intelligence algorithms in wireless sensor networks.  Performance Evaluation, 2011.
  * [6] G. Caravagna.  Formal Modeling and Simulation of Biological Systems with Delays.  Ph. D. thesis, Università di Pisa, 2011.
  * [7] P. Cisneros-Velarde and F. Bullo.  Multi-group SIS epidemics with simplicial and higher-order interactions, 2020.  https://arxiv.org/abs/2005.11404.
  * [8] G. F. de Arruda, F. A. Rodrigues, and Y. Moreno.  Fundamentals of spreading processes in single and multilayer complex networks.  Physics Reports, 756:1–59, 2018.
  * [9] K. Devriendt and P. Van Mieghem.  Unified mean-field framework for susceptible-infected-susceptible epidemics on networks, based on graph partitioning and the isoperimetric inequality.  Physical Review E, 96:052314, Nov 2017.
  * [10] R. J. Glauber.  Time‐dependent statistics of the Ising model.  Journal of Mathematical Physics, 4(2):294–307, 1963.
  * [11] B. Guerra and J. Gómez-Gardeñes.  Annealed and mean-field formulations of disease dynamics on static and adaptive networks.  Physical Review E, Statistical, nonlinear, and soft matter physics, 82:035101, 09 2010.
  * [12] R. A. Hayden, I. Horváth, and M. Telek.  Mean field for performance models with generally-distributed timed transitions.  In G. Norman and W. Sanders, editors, Quantitative Evaluation of Systems, volume 8657 of Lecture Notes in Computer Science, pages 90–105. Springer International Publishing, 2014.
  * [13] I. Iacopini, G. Petri, A. Barrat, and V. Latora.  Simplicial models of social contagion.  Nature Communications, 10:2485, 06 2019.
  * [14] B. Jhun, M. Jo, and B. Kahng.  Simplicial SIS model in scale-free uniform hypergraph.  arXiv, 2019.  https://arxiv.org/abs/1910.00375.
  * [15] I. Kiss, J. Miller, and P. Simon.  Mathematics of Epidemics on Networks, volume 46.  Springer, 01 2017.
  * [16] T. Kurtz.  Solutions of ordinary differential equations as limits of pure jump Markov processes.  Journal of Applied Probability, 7:49–58, 04 1970.
  * [17] T. G. Kurtz.  Strong approximation theorems for density dependent Markov chains.  Stochastic Processes and their Applications, 6(3):223 – 240, 1978\.
  * [18] C. Li, R. van de Bovenkamp, and P. Van Mieghem.  Susceptible-infected-susceptible model: A comparison of $n$-intertwined and heterogeneous mean-field approximations.  Physical Review E, 86:026116, Aug 2012.
  * [19] P. Mieghem.  The N-intertwined SIS epidemic network model.  Computing, 93:147–169, 12 2011.
  * [20] M. Molloy and B. Reed.  A critical point for random graphs with a given degree sequence.  Random Structures & Algorithms, 6(2‐3):161–180, 1995.
  * [21] J. Noonan and R. Lambiotte.  Dynamics of majority rule on hypergraphs.  Physical Review E, 104:024316, Aug 2021.
  * [22] M. Ostilli and F. Mukhamedov.  Continuous- and discrete-time Glauber dynamics. First- and second-order phase transitions in mean-field Potts models.  EPL (Europhysics Letters), 101(6):60008, 2013.
  * [23] N. Papanikolaou, G. Vaccario, E. Hormann, R. Lambiotte, and F. Schweitzer.  Consensus from group interactions: An adaptive voter model on hypergraphs.  Physical Review E, 105:054307, May 2022.
  * [24] R. Parasnis, R. Kato, A. Sakhale, M. Franceschetti, and B. Touri.  Usefulness of the age-structured SIR dynamics in modelling COVID-19, 2022.  https://arxiv.org/abs/2203.05111.
  * [25] R. Pastor-Satorras and A. Vespignani.  Epidemic spreading in scale-free networks.  Physical Review Letters, 86:3200–3203, Apr 2001.
  * [26] N. Perra, B. Gonçalves, R. Pastor-Satorras, and A. Vespignani.  Activity driven modeling of time varying networks.  Scientific reports, 2, 03 2012.
  * [27] R. Schlicht and G. Winkler.  A delay stochastic process with applications in molecular biology.  Journal of Mathematical Biology, 57(5):613–48, Nov. 2008.
  * [28] D. H. Silva, S. C. Ferreira, W. Cota, R. Pastor-Satorras, and C. Castellano.  Spectral properties and the accuracy of mean-field approaches for epidemics on correlated power-law networks.  Physical Review Research, 1:033024, Oct 2019.
  * [29] P. L. Simon and I. Z. Kiss.  On bounding exact models of epidemic spread on networks, 2017.  https://arxiv.org/abs/1704.01726.
  * [30] A. Sridhar and S. Kar.  Mean-field approximation for stochastic population processes in networks under imperfect information, 01 2021.  https://arxiv.org/abs/2101.09644.
  * [31] A. Sridhar and S. Kar.  On the accuracy of deterministic models for viral spread on networks.  ArXiv, 2021.  https://arxiv.org/abs/2104.04913.
  * [32] E. Szemeredi.  Regular partitions of graphs.  Proceedings of the C.N.R.S. International Colloquium, 260:11, 04 1975.
  * [33] L. Zino, A. Rizzo, and M. Porfiri.  An analytical framework for the study of epidemic models on activity driven networks.  Journal of Complex Networks, 5(6):924–952, 11 2017.