# Large population limit for a multilayer _SIR_ model including households and
workplaces

Madeleine Kubasch  Centre de mathématiques appliquées (CMAP), Ecole
Polytechnique, 91128 Palaiseau, France MaIAGE, INRAE, Université Paris-Saclay,
78350 Jouy-en-Josas, France

###### Abstract

We study a multilayer _SIR_ model with two levels of mixing, namely a global
level which is uniformly mixing, and a local level with two layers
distinguishing household and workplace contacts, respectively. We establish
the large population convergence of the corresponding stochastic process. For
this purpose, we use an individual-based model whose state space explicitly
takes into account the duration of infectious periods. This allows to deal
with the natural correlation of the epidemic states of individuals whose
household and workplace share a common infected. In a general setting where a
non-exponential distribution of infectious periods may be considered,
convergence to the unique deterministic solution of a measure-valued equation
is obtained. In the particular case of exponentially distributed infectious
periods, we show that it is possible to further reduce the obtained
deterministic limit, leading to a closed, finite dimensional dynamical system
capturing the epidemic dynamics. This model reduction subsequently is studied
from a numerical point of view. We illustrate that the dynamical system
derived from the large population approximation is a pertinent model reduction
when compared to simulations of the stochastic process or to an alternative
edge-based compartmental model, both in terms of accuracy and computational
cost.

## Introduction

Among the possibilities of non-pharmaceutical epidemic control measures,
social distancing may be used to slow down disease propagation. This has been
emphasized by the recent COVID-19 epidemic, which lead to wide-spread
implementation of non-pharmaceutical interventions. In particular, both
empirical evidence [18] and simulation studies [25] have pointed out that
school closure and teleworking are among the most effective of those measures.
These results are not only limited to COVID-19, as it is equivalently known
that particularly school closure mitigates influenza epidemics [15]. It thus
seems relevant to take into account this kind of control measure in epidemic
models, in order to study their impact on epidemic spread. A difficulty then
arises from the fact that social distancing often disrupts particular
contacts, such as workplace contacts only in the case of teleworking, while
leaving intact other types of contact, for instance within households. Di
Lauro _et al._ [9] have shown through the example of lockdowns that in order
to correctly model control measures impacting preferably specific types of
contact, it is necessary for the model to explicitly distinguish these contact
structures. In addition, even though it is well known that clustered contact
patterns influence epidemic spread [29], precisely understanding the impact on
disease propagation of the way individuals are organized in households and
workplaces is not straightforward [4].

This motivates the study of models which explicitly distinguish different
levels of contact, which as we shall see lead to interesting mathematical
issues due to their multiscale population structure. A minimal model taking
into account households and workplaces has been proposed by Pellis _et al._
[22], referred to as the household-workplace model. This model consists in a
stochastic _SIR_ model with two levels of mixing [7], the _global_ and the
_local_ one, the latter being characterized by the presence of small contact
structures. More precisely, contacts may take place within three different
layers: households, workplaces and the general population. Households and
workplaces correspond to the local level of mixing, while the general
population constitutes the global one. In particular, it is possible to
distinguish different rates at which acquaintances encounter one another
within each of these three layers. This model has already been studied to some
extent, establishing for instance different reproduction numbers, see _e.g._
Ball _et al._ [2] for an extensive overview. In particular, the reproduction
number $R_{I}$ proposed by Pellis _et al._ [22] is interesting, as it is
intimately connected to the proportion of infections taking place within each
layer of the contact network [4]. The $R_{0}$ used throughout this paper will
thus refer to this particular reproduction number, and computations of $R_{0}$
and proportions of infections per layer will make use of the working package
associated to [4].

One drawback of this model is its complexity, both mathematically and
numerically. Indeed, it is not simple to analyse due to correlations arising
as soon as an individual may belong to several small contact structures at
once. Also, simulations require a significant amount of computation time,
especially when considering larger population sizes. As a consequence, it is
of interest to develop reduced models, which may be more prone to theoretical
studies and/or numerical exploration. In a previous work, numerical evidence
indicates that this epidemic model is rather well approached by a
deterministic, uniformly mixing _SIR_ model whose parameters have been
calibrated using the epidemic growth rate of the household-workplace model
[4]. However, even though it captures some key characteristics of the
epidemic, this reduction has the disadvantage of not closely fitting the
epidemic trajectories over time. This paper thus aims at proposing another
reduced model, which more finely captures the epidemic dynamics of the
household-workplace model, and which will arise from the stochastic model as a
rigorously established large population limit.

Indeed, large population approximations of stochastic models have proven
fruitful to achieve such model reductions in many contexts, among which
epidemics on random graphs.

A well-known results of this kind concerns the large graph limit of the _SIR_
model spreading on a configuration graph. The corresponding large population
dynamics have been derived heuristically by Volz [28] and refined by Miller
[20], and convergence of the stochastic model to this limit has since been
established [8, 13]. This model reduction has the particularity of being
especially parsimonious, being a system of very few equations. As we will see,
the heuristic reasoning developed by Volz has given rise to a family of model
reductions in different settings, that are referred to as edge-based
compartmental models (EBCM).

A new difficulty arises as soon as one wants to take into account clustered
networks, containing entirely connected sub-graphs, also called _cliques_.
There exists a variant of the configuration model which allows to take into
account such cliques [21, 19], and Volz _et al._ [29] have proposed an EBCM of
an _SIR_ model spreading on such graphs. Using a different approach, St-Onge
_et al._ [26] propose an approximate master equation for a similar model. In
both cases, the obtained dynamical systems are high-dimensional, contrary to
the case of the classical configuration model. Further, epidemic parameters
may differ between cliques, for instance as (random) functions of clique size,
but these models are not enough to distinguish different levels of contacts.

This has been taken into consideration by models with two levels of mixing.
Originally, the local level would contain only one layer of cliques, often
interpreted as households. For instance, House and Keeling [11] have proposed
a deterministic model where each individual belongs to exactly one household,
and where the global level further is uniformly mixing. Since then, more
general models have been considered, taking into account a general population
represented by a configuration model, using either the EBCM formalism [9], or
the effective degree model [16]. However, in both cases, the local level
contains only one layer, and further within this layer, each individual
belongs to one clique only. Alternatively, Jacobsen _et al._ [12] have
developed an EBCM with an arbitrary number of layers of contacts, but as each
layer corresponds to a dynamic configuration model, cliques are not taken into
account. Finally, Barnard _et al._ [5] have developed an EBCM in the case of a
model with two levels of mixing, where the global level corresponds to a
dynamic configuration model, while at the local level, each individual may
belong to several cliques at once, but without distinguishing different kinds
of cliques.

In this paper, we will establish the large population approximation of the
model with two levels of mixing which explicitly distinguishes households and
workplaces. In order to do so, we will formalize the model in a finite
population as an individual-based stochastic process, and establish that this
process converges when the size of the population grows to infinity. This
allows, on the one hand, to identify our model reduction, and on the other
hand establishes that this model reduction is legitimate since it is
asymptotically exact. Notice here that it is not at all obvious that our model
of interest may be reduced to a closed, finite-dimensional dynamical system,
as each infected individual correlates the epidemic spread in his household
and workplace, being infectious for exactly the same period of time in both
structures. In order to deal with this dependence, the duration of infectious
periods will explicitly be taken into account in the model structure. This
difficulty actually arises as soon as one considers the possibility of an
individual belonging to several cliques at once, whether they are of different
types or not, and we refer to [3] where a similar approach has been developed
for branching approximations. This model formulation amounts to studying a
measure-valued process mixing discrete and continuous components. More
precisely, we establish the convergence of the individual-based process to the
unique solution of an explicit measure-valued equation. This convergence
result is not limited to the case of the Markovian _SIR_ model as the duration
of infectious periods may follow an arbitrary distribution. In the particular
case where this distribution is exponential, it is possible to go one step
further and reduce the epidemic dynamics to a closed, finite dimensional
dynamical system which is similar in spirit to reductions proposed in related
settings [11, 26].

The present paper is structured as follows. Section 1 introduces the
individual-based model, and Section 2.1 subsequently presents the convergence
results in detail. Section 2.2 is devoted to numerical aspects. We first
illustrate that the obtained dynamical system is in good accordance with
stochastic simulations, discuss its implementation and examine its
computational cost in terms of computation time compared to stochastic
simulations. Next, we confront our reduced model to an alternative model
reduction which we obtain using the EBCM formalism. Finally, Section 3
contains the proofs of our results.

Before proceeding, let us introduce some notations that will be used
throughout the paper. For any integers $n\leq m$, we write $\llbracket
n,m\rrbracket=\\{n,\cdots,m\\}$. For a measurable space $(E,\mathcal{E})$, let
$\mathcal{M}_{P}(E)$ be the set of point measures, $\mathcal{M}_{F}(E)$ the
set of finite measures and $\mathcal{M}_{1}(E)$ the set of probability
measures on $E$. We define
$\mathcal{M}_{P,1}(E)=\mathcal{M}_{P}(E)\cap\mathcal{M}_{1}(E)$ and
$\mathcal{M}_{P,F}(E)=\mathcal{M}_{P}(E)\cap\mathcal{M}_{F}(E)$. For a measure
$\mu$ on $E$ and a suitable function $f$ (either non-negative or belonging to
$L^{1}(\mu)$), let $\langle\mu,f\rangle=\int_{E}fd\mu$. Also, for $x\in E$,
$\delta_{x}$ designates the Dirac measure at point $x$. Further, for any
metric space $E$ and any integer $m$, let $\mathcal{C}(E,\mathbb{R}^{m})$ be
the set of continuous functions $f:E\to\mathbb{R}^{m}$. Similarly,
$\mathcal{C}_{b}(E,\mathbb{R}^{m})$ is defined as the subset of bounded
functions $f\in\mathcal{C}(E,\mathbb{R}^{m})$. Finally, the space
$\mathcal{C}^{1}_{b}(E,\mathbb{R}^{m})$ designates the set of bounded
functions $f:E\to\mathbb{R}^{m}$ such that $f$ is differentiable and its
differential is continuous and bounded.

##  1 Presentation of the model

Let us begin by introducing the epidemic model of interest, in two successive
steps. At first, a general model description is yielded, which corresponds to
a more intuitive presentation of the model, before stating the mathematical
model in detail as measure-valued stochastic differential equation.

###  1.1 General presentation of the model

Let us start by describing the population structure of interest. Consider a
population of $K$ individuals. Each individual is part of exactly one
household and one workplace, which are chosen independently from one another,
and independently for each individual.

More precisely, such a population structure can be obtained as described in
[4]. Suppose that households and workplaces are of size at least one and at
most $n_{\max}$. Consider distributions $(\pi^{H}_{j})$ and $(\pi^{W}_{j})$ on
$\llbracket 1,n_{\max}\rrbracket$. Conceptually, these distributions will
correspond to _asymptotic_ household and workplace size distributions, in the
sense that in an infinite population, a proportion $\pi^{H}_{j}$ of households
would be of size $j$, while $\pi^{W}$ would play a similar role for
workplaces. In such an infinite population, the average household and
workplace sizes, respectively $m_{H}$ and $m_{W}$, would be given for
$X\in\\{H,W\\}$ by

$m_{X}=\sum_{j=0}^{n_{\max}}j\pi^{X}_{j}.$

In order to assign each individual of the finite-size population a household,
we proceed as follows. Let $k\in\llbracket 0,K\rrbracket$ be the number of
individuals who are not yet member of a household. While $k>0$, choose a size
$\tilde{n}$ according to $\pi^{H}$, independently from the household sizes
that were chosen during previous steps. The newly uncovered household is then
of size $n=\min(\tilde{n},k)$, and $n$ individuals out of the $k$ remaining
ones are picked uniformly at random to assemble this new household.
Consequently, it remains to update $k$ to $k-n$. The process stops as soon as
$k=0$, as all individuals then belong to a household. Finally, this process is
repeated independently for workplaces, using $\pi^{W}$ instead of $\pi^{H}$.
Throughout the paper, let us consider one particular sequence
$(\mathbf{G}^{K})_{K\geq 1}$ of realisations of this random population
structure. The choice of this sequence being arbitrary among all sequences of
such population structures, this does not lessen the generality of our
results.

It remains to describe the way the disease spreads in the population. The
epidemic model considered here is an extension of the standard _SIR_ model. At
each time, each individual is either _susceptible_ if he has never encountered
the disease and may be contaminated; _infected_ if the individual is currently
infectious, in which case he may transmit the disease to other susceptibles;
or _recovered_ , once the infectious period is over, in which case the
individual has become immune against the disease.

The disease is transmitted among individuals as follows. Any contact between a
susceptible and an infected individual within either level of mixing amounts
in the contamination of the susceptible individual. Such an event is called an
_infectious contact_. Within each household, each workplace and the general
population, uniform mixing is assumed, but the parameterization differs
slightly between the layers. Indeed, for households, we consider a _one-to-
one_ contact rate $\lambda_{H}$, meaning that whenever there are $s$
susceptibles and $i$ infected within a household, the next infectious contact
occurs at rate $\lambda_{H}si$. Similarly, another one-to-one contact rate
$\lambda_{W}$ is associated to workplace contacts. Within the general
population, on the other hand, a _one-to-all_ contact rate $\beta_{G}$ is
considered: when there are $s$ susceptibles and $i$ infected within a
population of size $K$, infectious contacts occur at rate
$\frac{\beta_{G}}{K}si$. Here, contrary to the local level, frequency-
dependence is considered for the global contact rate, as each given encounter
between two individuals chosen at random within the general population is less
likely then within small structures.

Finally, infected individuals remain infectious for a period of time
distributed according to a probability distribution $\nu$ on $\mathbb{R}_{+}$.
Once they recover, they are supposed to be immune against the disease from
there on. In particular, if $\nu$ is an exponential distribution, this
corresponds to the Markovian _SIR_ model.

###  1.2 The epidemic model at the level of households and workplaces

As we aim at investigating the large population limit of this model, we choose
to enrich the population description as to obtain a closed Markov process.
This corresponds to a favourable mathematical setting, as it allows us to use
the associated martingale problem. In order to do so, we will represent the
population in terms of particles which are described by a type. It seems
natural to consider particles which correspond to entire structures, _i.e._
households and workplaces, which are characterized by their size and the
number of susceptible and infected individuals they contain. Indeed, this
point of view has already proven useful for deriving reproduction numbers for
related models [2]. However, this is not enough to obtain a closed system of
Markovian dynamics. The problem is that each infected individual correlates
the spread of the epidemic within his household and his workplace, leading to
an intricate correlation network. In order to circumvent this difficulty,
similarly to [3], we will thus further characterize each structure by the
infectious periods of its infected members.

More precisely, for a population of size $K\geq 1$, let $K_{H}$ be the number
of households and $K_{W}$ the number of workplaces in $\mathbf{G}^{K}$. Label
the $K_{H}$ households in an arbitrary fashion $1,\dots,K_{H}$. Consider the
set

$E=\left\\{(n,s,\tau)\in\llbracket 1,n_{\max}\rrbracket\times\llbracket
0,n_{\max}\rrbracket\times\mathbb{R}^{n_{\max}}:s\leq n;\;\forall
j>n-s,\tau_{j}=0\right\\}.$

Then for $k\in\llbracket 1,K_{H}\rrbracket$, the $k$-th household is
characterized at time $t\geq 0$ by its type

$x^{H}_{k}(t)=(n^{H}_{k},s^{H}_{k}(t),\tau^{H}_{k}(t))\in E.$

The first two components of $x^{H}_{k}$ simply correspond respectively to the
size of the household (which is constant over time), and the number of
susceptible members of the household at time $t$. The third component
$\tau^{H}_{k}$ is a vector containing the _remaining infectious periods_ of
the members of the household and keeps track of infected and removed
individuals. Indeed, at time $t$, there are $n^{H}_{k}-s^{H}_{k}(t)$ infected
or removed individuals within the household. For $j\in\llbracket
1,n^{H}_{k}-s^{H}_{k}(t)\rrbracket$, $\tau^{H}_{k,j}(t)$ corresponds to the
remaining duration of the infectious period at time $t$ of the $j$-th member
of the household who has contracted the disease. If $\tau^{H}_{k,j}(t)>0$, the
individual is still infectious and will remain so for $\tau^{H}_{k,j}(t)$
units of time. Otherwise, if $\tau^{H}_{k,j}(t)\leq 0$, the individual has
recovered, and the recovery has occurred at time $t-|\tau^{H}_{k,j}(t)|$. For
$j>n^{H}_{k}-s^{H}_{k}(t)$, $\tau^{H}_{k,j}(t)$ is set to zero and can be
neglected for interpretation. Notice that, in principle, it would have been
possible to define $\tau^{H}_{k}(t)\in\mathbb{R}^{n^{H}_{k}-s^{H}_{k}(t)}$,
but letting $\tau^{H}_{k}\in\mathbb{R}^{n_{\max}}$ is convenient for
computations.

Similarly, label the $K_{W}$ workplaces in an arbitrary order $1,\dots,K_{W}$.
For $\ell\in\llbracket 1,K_{W}\rrbracket$, the $\ell$-th workplace is
characterized by its type
$x^{W}_{\ell}(t)=(n^{W}_{\ell},s^{W}_{\ell}(t),\tau^{W}_{\ell}(t))$, which is
defined analogously to household types.

By definition, these types evolve over time. On the one hand, for any
$X\in\\{H,W\\}$, for any $k\in\llbracket 1,K_{X}\rrbracket$ and
$j\in\llbracket 1,n_{\max}\rrbracket$, the $j$-th component of $\tau^{X}_{k}$
decreases linearly at unitary rate if it describes the remaining infectious
period of an individual having contracted the disease at some previous time,
and stays constant otherwise:

$\forall j\in\llbracket
1,n_{\max}\rrbracket,\;\frac{d}{dt}\tau^{X}_{k,j}(t)=-\mathbf{1}_{\\{j\leq
n^{X}_{k}-s^{X}_{k}(t)\\}}.$

Let $(e_{j})_{1\leq j\leq n_{\max}}$ denote the canonical basis of
$\mathbb{R}^{n_{\max}}$. Then for any $0\leq t\leq T$, and $x=(n,s,\tau)\in
E$, we may define $\Psi(x,T,t)$ as the type of a structure at time $T$ given
that it was in state $x$ at time $t$, supposing that no infections occurred in
the meantime:

$\Psi(x,T,t)=\left(n,s,\tau-\sum_{j=1}^{n-s}(T-t)e_{j}\right).$

On the other hand, infections within each level of mixing also cause the
modification of the types of the household and the workplace of newly
contaminated individuals. More precisely, consider a contamination occurring
at time $t$. Suppose that the newly infected belongs to the $k$-th household
and $\ell$-th workplace. Let $\sigma$ be the realisation of a random variable
of distribution $\nu$, which is drawn independently for each new infected.
Then $x^{H}_{k}$ and $x^{W}_{\ell}$ jump from $x^{H}_{k}(t-)$ and
$x^{W}_{\ell}(t-)$ to $\mathfrak{j}(x^{H}_{k}(t-),\sigma)$ and
$\mathfrak{j}(x^{W}_{\ell}(t-),\sigma)$ respectively, where for any
$x=(n,s,\tau)\in E$,

$\mathfrak{j}(x,\sigma)=\left(n,s-1,\tau+\sigma e_{n-s+1}\right).$

It remains to describe how one identifies the household $k$ and workplace
$\ell$ the newly infected belongs to. Let $S(t-)$ be the number of
susceptibles in the population previously to the infection event. If it takes
place within the general population, any susceptible individual is chosen with
uniform probability to be contaminated. The newly infected thus belongs to the
$k$-th household with probability $s^{H}_{k}(t-)/S(t-)$, and independently to
the $\ell$-th workplace with probability $s^{W}_{\ell}(t-)/S(t-)$. Similarly,
if the infection occurs within a household, only the workplace of the newly
infected needs to be uncovered, and corresponds to the $\ell$-th workplace
with the same probability as previously. Within-workplace infections are
treated analogously.

We are now ready to properly introduce the stochastic process which represents
the previously described dynamics. We are interested in the stochastic process
$(\zeta^{K}_{t}=(\zeta^{H|K}_{t},\zeta^{W|K}_{t}))_{t\geq 0}$ taking values in
$\mathfrak{M}_{P,1}=\mathcal{M}_{P,1}(E)\times\mathcal{M}_{P,1}(E)$, whose
definition will be detailed hereafter. Here, $\zeta^{H|K}_{t}$ and
$\zeta^{W|K}_{t}$ correspond respectively to the normalized counting measures
associated to the distributions of household and workplace types at time $t$,
_i.e._ for any time $t\geq 0$ and $X\in\\{H,W\\}$,

$\zeta^{X|K}_{t}=\frac{1}{K_{X}}\sum_{k=1}^{K_{X}}\delta_{x^{X}_{k}(t)}.$

Start by noticing that, as both household and workplace sizes are bounded, for
$X\in\\{H,W\\}$, the following inequality holds:

$\frac{K}{n_{\max}}\leq K_{X}\leq K.$

Thus, studying the asymptotic $K\to\infty$ amounts to
$(K_{H},K_{W})\to(\infty,\infty)$. It actually holds that for any Borel set
$A\subset E$, if one were to count at time $t\geq 0$ the number of households
or workplaces of type $x\in A$ in an infinite population, then either there
are none, or there are $O(K)$ of them. The model thus is well suited for
studying the large population limit.

Recall that in the case of infections occurring within the general population,
the rate of infections depends on the number of susceptibles $S(t)$ and
infected $I(t)$ within the whole population at time $t\geq 0$. It is easy to
compute these quantities using $\zeta^{K}$ as follows. For any
$\tau\in\mathbb{R}^{n_{\max}}$, let
$i(\tau)=\sum_{k=1}^{n_{\max}}\mathbf{1}_{\left\\{\tau_{k}>0\right\\}}$. Then
for any $t\geq 0$,

$I^{H}(t)=\frac{1}{K_{H}}\sum_{k=1}^{K_{H}}i(\tau^{H}_{k}(t))$

corresponds to the average number of infected individuals per household at
time $t$. Similarly, one may define $S^{H}(t)$ as the average number of
susceptibles per household at time $t$. Then

$\forall X\in\\{H,W\\},\;S(t)=K_{X}S^{X}(t)\text{ and }I(t)=K_{X}I^{X}(t).$
(1)

Further, let $N^{H}$ be the average household size, which is constant over
time and always equal to $K/K_{H}$. This leads to $I(t)/K=I^{X}(t)/N^{X}$,
which will be of use in computations. Notice that we will need to check that
Equation (1) is well posed, as equalities of the type
$K_{H}S^{H}(t)=K^{W}S^{W}(t)$ technically need to be proven for the stochastic
process formalizing the model.

Finally, let us briefly emphasize that the partition of the population in
households and workplaces is entirely conveyed by
$\zeta^{K}_{0}\in\mathfrak{M}_{P,1}$, as it does not vary over time. In
particular, the proportions of households and workplaces of each size are
supposed to correspond to those observed in $\mathbf{G}^{K}$. Similarly, there
are some natural constraints on $\zeta^{K}_{0}$, as $\zeta^{H|K}$ and
$\zeta^{W|K}$ describe the same population, once dispatched into households,
and once dispatched into workplaces. As these assumptions are intuitive, they
will not all be detailed here, and they will only be emphasized when
necessary. For instance, the total number of members in some epidemic state
(susceptible, infected or recovered) within all households is equal to the
total number of members in this state within all workplaces. Hence
$K_{H}S^{H}(0)=K_{W}S^{W}(0)$ and $K_{H}I^{H}(0)=K_{W}I^{W}(0)$ almost surely.
Further, at time $0$, each infected needs to have the same remaining
infectious period in both his household and workplace. In other words, almost
surely,

$\small\\{\tau^{H}_{k,j}(0):1\leq k\leq K_{H},\,1\leq j\leq
n^{H}_{k}-s^{H}_{k}(0)\\}=\\{\tau^{W}_{\ell,j}(0):1\leq\ell\leq K_{W},\,1\leq
j\leq n^{W}_{\ell}-s^{W}_{\ell}(0)\\}.$ (2)

Forthcoming Lemma 1.2 shows that these conditions are enough to ensure that
the previous characterization of $S(t)$ and $I(t)$ in terms of $S^{X}(t)$ and
$I^{X}(t)$ is legitimate.

Before giving the proper definition of $\zeta^{K}$, let us introduce some
necessary notations. Let $U^{G}_{K}=(\mathbb{R}_{+})^{3}\times\llbracket
1,K_{H}\rrbracket\times\llbracket 1,K_{W}\rrbracket\times\mathbb{R}_{+}$, and
consider the following measure on $U^{G}_{K}$:

$\mu^{G}_{K}(du)=\mu^{G}_{K}(d\bm{\theta},dk,d\ell,d\sigma)=d\bm{\theta}\otimes\mu_{\\#}(dk)\otimes\mu_{\\#}(d\ell)\otimes\nu(d\sigma),$

where $d\bm{\theta}$ and $\mu_{\\#}$ denote the Lebesgue measure on
$\mathbb{R}^{3}$ and the standard counting measure, respectively. For $t\geq
0$ and $u=(\bm{\theta},k,\ell,\sigma)\in U^{G}_{K}$ where
$\bm{\theta}=(\theta_{1},\theta_{2},\theta_{3})$, let us define

$\mathcal{I}_{G}(t,u)=\mathbf{1}_{\left\\{\theta_{1}\leq\frac{\beta_{G}}{K}S(t)I(t),\,\theta_{2}\leq\frac{s^{H}_{\ell}(t)}{S(t)},\,\theta_{3}\leq\frac{s^{W}_{\ell}(t)}{S(t)}\right\\}}.$

The idea is that $\mathcal{I}_{G}$ will yield the correct rate for infection
events in the general population. More precisely, the constraint on
$\theta_{1}$ corresponds to the rate of infectious contacts at that level of
mixing, while the constraints on $\theta_{2}$ and $\theta_{3}$ are related to
the probability that the newly infected belongs to the $k$-th household and
$\ell$-th workplace.

Similarly, let $U_{K}=(\mathbb{R}_{+})^{2}\times\llbracket
1,K_{H}\rrbracket\times\llbracket 1,K_{W}\rrbracket\times\mathbb{R}_{+}$,
endowed with the measure

$\mu_{K}(du)=\mu_{K}(d\bm{\theta},dk,d\ell,d\sigma)=d\bm{\theta}\otimes\mu_{\\#}(dk)\otimes\mu_{\\#}(d\ell)\otimes\nu(d\sigma),$

where with slight abuse of notation, $d\bm{\theta}$ designates the Lebesgue
measure on $\mathbb{R}^{2}$. Then for $t\geq 0$ and
$u=(\bm{\theta},k,\ell,\sigma)\in U_{K}$ where
$\bm{\theta}=(\theta_{1},\theta_{2})$, we further introduce

$\mathcal{I}_{H}(t,u)=\mathbf{1}_{\left\\{\theta_{1}\leq\lambda_{H}s^{H}_{k}(t)i(\tau^{H}_{k}(t)),\,\theta_{2}\leq\frac{s^{W}_{\ell}(t)}{S(t)}\right\\}}.$

This time, the constraints on $\theta_{1}$ and $\theta_{2}$ correspond
respectively to the rate of infection within the $k$-th household and the
probability of the newly infected belonging to the $\ell$-th workplace. Also,
for any $T\geq t\geq 0$ and $u=(\bm{\theta},k,\ell,\sigma)\in U^{G}_{K}\cup
U_{K}$, consider the following quantity, which will allow to keep track of the
change in the household population due to an infection within the $k$-th
household:

$\Delta_{H}(u,T,t)=\delta_{\left(\Psi(\mathfrak{j}(x^{H}_{k}(t-),\sigma),T,t))\right)}-\delta_{\left(\Psi(x^{H}_{k}(t-),T,t))\right)}.$

Finally, define $\mathcal{I}_{W}(t,u)$ and $\Delta_{W}(u,T,t)$ analogously, by
interchanging $H$ and $W$ as well as $k$ and $\ell$:

$\small\mathcal{I}_{W}(t,u)=\mathbf{1}_{\left\\{\theta_{1}\leq\lambda_{W}s^{W}_{\ell}(t)i(\tau^{W}_{\ell}(t)),\,\theta_{2}\leq\frac{s^{H}_{k}(t)}{S(t)}\right\\}}\;\text{and}\;\Delta_{W}(u,T,t)=\delta_{\left(\Psi(\mathfrak{j}(x^{W}_{\ell}(t-),\sigma),T,t))\right)}-\delta_{\left(\Psi(x^{W}_{\ell}(t-),T,t))\right)}.$

In order to simplify notations in the following, for $Y\in\\{H,W\\}$, we let
$U^{Y}_{K}=U_{K}$ and $\mu^{Y}_{K}=\mu_{K}$. We are now ready to state the
main definition of $\zeta^{K}_{t}$.

######  Proposition 1.1.

Define on the same probability space as $\zeta^{K}_{0}$, and independently
from $\zeta^{K}_{0}$, three independent Poisson point measures $Q^{Y|K}$ of
intensity $dt\;\mu^{Y}_{K}(du)$ taking values in $\mathbb{R}_{+}\times
U^{Y}_{K}$, for $Y\in\\{H,W,G\\}$.

Then $\zeta^{K}=(\zeta^{H|K},\zeta^{W|K})$ is defined as the unique strong
solution taking values in $\mathbb{D}(\mathbb{R}_{+},\mathfrak{M}_{P,1})$ of
the following equation. For $X\in\\{H,W\\}$,

$\zeta^{X|K}_{T}=\frac{1}{K_{X}}\left(\sum_{j=1}^{K_{X}}\delta_{\Psi(x^{X}_{j}(0),T,0)}+\sum_{Y\in\\{H,W,G\\}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)\Delta_{X}(u,T,t)Q^{Y|K}(dt,du)\right).$
(3)

The idea behind Equation (3) goes as follows. Let us focus for example on the
distribution $\zeta^{H|K}_{T}$ of household types at time $T$. Each individual
household’s type contributes to the distribution at uniform weight $1/K_{X}$.
If no infection event occurs between times $0$ and $T$, then the state of the
$k$-th household at time $T$ is given by $\Psi(x^{H}_{k}(0),T,0)$. However,
suppose now that before time $T$, at least one initially susceptible member of
the $k$-th household is infected, and let $t$ be the first time at which such
an event occurs. Then $x^{H}_{k}(t)=\mathfrak{j}(x^{H}_{k}(t-),\sigma)$ where
$\sigma$ is distributed according to $\nu$. If no other infections affect this
household up to time $T$, it will be in state
$\Psi(\mathfrak{j}(x^{H}_{k}(t-),\sigma),T,t)$ instead of
$\Psi(x^{H}_{k}(0),T,0)=\Psi(x^{H}_{k}(t-),T,t)$. This reasoning can be
iterated over the whole of $[0,T]$ and is reflected in $\Delta_{H}$. Finally,
the terms $\mathcal{I}_{Y}$ for $Y\in\\{G,H,W\\}$ assure that all infection
events occur at the corresponding rates.

###### Proof.

The proof uses classical arguments, which will only be outlined here. Start by
establishing existence of $\zeta^{K}$. Consider the sequence $(T_{n})_{n\geq
0}$ of successive jump times of $\zeta^{K}$, where we define $T_{0}=0$. Then
using a method similar to rejection sampling (e.g. [27, Section 2.3.1]),
$(T_{n})_{n\geq 0}$ can be obtained as a subsequence of the jump times of a
Poisson process with intensity
$n_{\max}(\lambda_{H}n_{\max}+\lambda_{W}n_{\max}+\beta_{G})K$, whose only
limiting value is $+\infty$. Thus $\lim_{n\to+\infty}T_{n}=+\infty$ almost
surely, ensuring that $\zeta^{K}$ takes values in
$\mathbb{D}(\mathbb{R}_{+},\mathcal{M}_{P,1}(E)\times\mathcal{M}_{P,1}(E))$.

Finally, uniqueness is obtained by an induction argument (e.g. [27,
Proposition 2.2.6]) which proves that for any $n\geq 0$,
$(T_{n},\zeta^{K}_{T_{n}})$ is uniquely determined by
$(\zeta^{K}_{0},Q^{H|K},Q^{W|K},Q^{G|K})$. This obviously is true for $n=0$.
The induction step relies on the observation that $T_{n+1}$ is uniquely
determined by $(T_{n},\zeta^{K}_{T_{n}},Q^{H|K},Q^{W|K},Q^{G|K})$ and
$\zeta^{K}_{T_{n+1}}$ by
$(T_{n+1},T_{n},\zeta^{K}_{T_{n}},Q^{H|K},Q^{W|K},Q^{G|K})$. The induction
hypothesis allows to conclude. ∎

Let us briefly show here that it follows directly from Proposition 1.1 that
Equation (1) is well posed.

######  Lemma 1.2.

Suppose that almost surely, $H_{K}S^{H}(0)=H_{W}S^{W}(0)$ and Equation (2)
holds. Then for any $t\geq 0$, $H_{K}S^{H}(t)=H_{W}S^{W}(t)$ and
$H_{K}I^{H}(t)=H_{W}I^{W}(t)$, almost surely.

###### Proof.

Let $T\geq 0$, and for any $x=(n,s,\tau)\in E$, let $\mathbf{n}(x)=n$,
$\mathbf{s}(x)=s$ and $\mathbf{i}(x)=i(\tau)$. Start by focusing on
$S^{X}(T)=\langle\zeta^{X|K}_{T},\mathbf{s}\rangle$ for $X\in\\{H,W\\}$. It
follows from Equation (3) that

$K_{X}S^{X}(T)=\sum_{j=1}^{K_{X}}\mathbf{s}(\Psi(x^{X}_{j}(0),T,0))+\sum_{Y\in\\{H,W,G\\}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)\langle\Delta_{X}(u,T,t),\mathbf{s}\rangle
Q^{Y|K}(dt,du).$

Notice that on the one hand, for any $x\in E$ and $0\leq t\leq T$,
$\mathbf{s}(\Psi(x,T,t))=\mathbf{s}(x)$. Hence the first term of the right-
hand side equals $K_{X}S^{X}(0)$, and for any $u=(\bm{\theta},k,\ell,\sigma)$,

$\langle\Delta_{H}(u,T,t),\mathbf{s}\rangle=\mathbf{s}(\Psi(\mathfrak{j}(x^{H}_{k}(t-),\sigma),T,t))-\mathbf{s}(\Psi(x^{H}_{k}(t-),T,t))=-1.$

The analogous computation yields
$\langle\Delta_{W}(u,T,t),\mathbf{s}\rangle=-1$. Thus
$K_{H}S^{H}(T)=K_{W}S^{W}(T)$ almost surely.

Let us now turn to $I^{X}(T)=\langle\zeta^{X|K}_{T},\mathbf{i}\rangle$. This
time, for any $u=(\bm{\theta},k,\ell,\sigma)$ and $0\leq t\leq T$, it holds
that
$\langle\Delta_{X}(u,T,t),\mathbf{i}\rangle=\mathbf{1}_{\left\\{\sigma>(T-t)\right\\}}$.
Finally, Equation (2) ensures that

$\sum_{k=1}^{K_{H}}\mathbf{i}(\Psi(x^{H}_{k}(0),T,0))=\sum_{k=1}^{K_{H}}\sum_{j=1}^{n^{H}_{k}-s^{K}_{k}(0)}\mathbf{1}_{\left\\{\tau^{H}_{k,j}(0)>T\right\\}}=\sum_{\ell=1}^{K_{W}}\mathbf{i}(\Psi(x^{W}_{\ell}(0),T,0)).$

The conclusion follows as previously from Equation (3). ∎

##  2 Main results

Let $\mathfrak{M}_{1}=\mathcal{M}_{1}(E)\times\mathcal{M}_{1}(E)$. In this
section, we are going to present our main results on the convergence of
$(\zeta^{K})_{K\geq 1}$ in the Skorokhod space
$\mathbb{D}\left(\mathbb{R}_{+},\mathfrak{M}_{1}\right)$. From now on, for
$X\in\\{H,W\\}$, let $\overline{X}$ be the complementary structure type,
_i.e._ $\overline{X}=W$ if $X=H$ and vice-versa.

###  2.1 Large population approximation of $(\zeta^{K})_{K\geq 1}$

For any $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$, let
$f_{t}(x)=f(t,x)$ for every $(t,x)\in\mathbb{R}_{+}\times E$, and consider the
differential operator $\mathcal{A}$ defined as follows. For any
$x=(n,s,\tau)\in E$,

$\mathcal{A}f_{t}(x)=\partial_{t}f(t,x)-\sum_{k=1}^{n-s}\partial_{\tau_{k}}f(t,x).$

Also, for any $x=(n,s,\tau)\in E$, let $\mathbf{n}(x)=n$, $\mathbf{s}(x)=s$
and $\mathbf{i}(x)=i(\tau)$ be the functions which to a structure in state $x$
associate the corresponding structure size, number of susceptible and number
of infected members, respectively. For instance, for any $X\in\\{H,W\\}$ the
average rage of within-structure infections at time $t$ is given by

$\lambda_{X}\langle\zeta^{X}_{t},\mathbf{s}\mathbf{i}\rangle=\frac{\lambda_{X}}{K_{X}}\sum_{k=1}^{K_{X}}\mathbf{s}(x^{X}_{k}(t))\mathbf{i}(x^{X}_{k}(t)).$

Notice also that as mentioned previously, the average size of structures is
constant over time, hence
$\langle\zeta^{X}_{t},\mathbf{n}\rangle=\langle\zeta^{X}_{0},\mathbf{n}\rangle$
for all $t\geq 0$. We are now ready to state our first result, whose proof is
postponed to Section 3.1.

######  Theorem 2.1.

Suppose that $\left(\zeta^{K}_{0}\right)_{K\geq 1}$ converges in law to
$\eta_{0}\in\mathfrak{M}_{1}$. Then $(\zeta^{K})_{K\geq 1}$ converges in
$\mathbb{D}\left(\mathbb{R}_{+},\mathfrak{M}_{1}\right)$ to
$\eta=(\eta^{H},\eta^{W})$ defined as the unique solution of the following
system of Equations (4). For any $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times
E,\mathbb{R})$, for any $T\geq 0$,

$\displaystyle\langle\eta^{X}_{T},f_{T}\rangle$
$\displaystyle=\langle\eta^{X}_{0},f_{0}\rangle+\int_{0}^{T}\langle\eta^{X}_{t},\mathcal{A}f_{t}\rangle
dt+\lambda_{X}\int_{0}^{T}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt$ (4)
$\displaystyle+\lambda_{\overline{X}}\int_{0}^{T}\frac{\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle}{\langle\eta^{\overline{X}}_{t},\mathbf{s}\rangle}\langle\eta^{X}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt+\beta_{G}\int_{0}^{T}\frac{\langle\eta^{H}_{t},\mathbf{i}\rangle}{\langle\eta^{H}_{0},\mathbf{n}\rangle}\langle\eta^{X}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt,$

where for any $(t,x)\in\mathbb{R}_{+}\times E$,
$f^{\mathcal{I}}_{t}(x)=\langle\nu,f_{t}(\mathfrak{j}(x,\cdot))\rangle$.

From now on, let us assume that $\nu$ is the exponential distribution of
parameter $\gamma$. As we shall see, it then is possible to deduce from
Theorem 2.1 that the proportion of susceptible and infected individuals in the
population converges to the solution of a dynamical system, when the size of
the population grows large.

Let $s(t)$ and $i(t)$ be the proportions of susceptible and infectious
individuals, respectively, in the population at time $t$ according to
distribution $\eta_{t}$. Further, for $N\in\mathbb{N}$, introduce the set

$\mathbb{S}=\left\\{(n-i,i):2\leq n\leq N,0\leq i\leq n-1\right\\}.$

For $(S,I)\in\mathbb{S}$, let $n^{H}_{(S,I)}(t)$ be the proportion of
households containing $S$ susceptible and $I$ infected individuals at time
$t$, according to distribution $\eta^{H}_{t}$. Define $n^{W}_{(S,I)}(t)$
analogously for workplaces. Finally, consider

$\tau_{G}(t)=\beta_{G}i(t),\text{ and
}\tau_{X}(t)=\frac{\lambda_{X}}{m_{X}}\sum_{(S,I)\in\mathbb{S}}SI\;n^{X}_{(S,I)}(t)\text{
for }X\in\\{H,W\\}.$

We assume that at time $0$, a fraction $\varepsilon$ of uniformly chosen
individuals are infected amidst an otherwise susceptible population.
Furthermore, at time $0$, the remaining infectious period of each infected
individual is supposed to be distributed according to $\nu$, independently
from one another. In particular, this setting corresponds to the following
probability distribution
$\eta_{0,\varepsilon}=(\eta^{H}_{0,\varepsilon},\eta^{W}_{0,\varepsilon})\in\mathfrak{M}_{1}$,
where for $X\in\\{H,W\\}$:

$\eta^{X}_{0,\varepsilon}(dn,ds,d\tau)=\pi^{X}_{n}\binom{n}{s}(1-\varepsilon)^{s}\varepsilon^{n-s}\delta_{(n,s)}(dn,ds)\otimes\left(\nu^{\otimes(n-s)}\otimes\delta_{0}^{\otimes\left(n_{\max}-n+s)\right)}(d\tau)\right).$
(5)

It then is possible to describe the epidemic dynamics by a finite, closed set
of ordinary differential equations, as shown in the following result whose
proof is postponed to Section 3.2.

######  Theorem 2.2.

Let $\varepsilon>0$. Suppose that $\nu$ is the exponential distribution of
parameter $\gamma$, and that $\left(\zeta^{K}_{0}\right)_{K\geq 1}$ converges
in law to $\eta_{0,\varepsilon}$. Then the functions
$(s,\,i,\,n^{X}_{(S,I)}:X\in\\{H,W\\},(S,I)\in\mathbb{S})$ are characterized
as being the unique solution of the following dynamical system: for any $t\geq
0$, for any $X\in\\{H,W\\}$, $(S,I)\in\mathbb{S}$,

$\displaystyle\begin{array}[]{l}{\frac{d}{dt}s(t)=-(\tau_{H}(t)+\tau_{W}(t)+\tau_{G}(t)s(t)),}\end{array}$
(6b)
$\displaystyle\begin{array}[]{l}{\frac{d}{dt}i(t)=-\frac{d}{dt}s(t)-\gamma
i(t),}\end{array}$ (6d) $\displaystyle\begin{array}[]{l}{\begin{aligned}
\frac{d}{dt}n^{X}_{(S,I)}(t)&=-\left(\lambda_{X}SI+\tau_{\overline{X}}(t)\frac{S}{s(t)}+\tau_{G}(t)S+\gamma
I\right)n^{X}_{(S,I)}(t)\\\
&\;+\gamma(I+1)n^{X}_{(S,I+1)}(t)\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}}\\\
&\;+\left(\lambda_{X}(S+1)(I-1)+\tau_{\overline{X}}(t)\frac{S+1}{s(t)}+\tau_{G}(t)(S+1)\right)n^{X}_{(S+1,I-1)}(t)\mathbf{1}_{\left\\{I\geq
1\right\\}},\end{aligned}}\end{array}$ (6f)

with initial conditions given by

$s(0)=1-\varepsilon;\;\;i(0)=\varepsilon;\;\;n^{X}_{(S,I)}(0)=\binom{S+I}{I}\pi^{X}_{S+I}(1-\varepsilon)^{S}\varepsilon^{I}.$
(7)

This dynamical system may be understood as follows. Equation (6b) corresponds
to the fact that the proportion of susceptibles decreases whenever a new
infection occurs within the general population, or within a household or
workplace. Similarly, Equation (6d) stems from the fact that newly
contaminated individuals move from the susceptible to the infected state,
which they in turn leave at rate $\gamma$ since the infectious period is
exponentially distributed. It remains to take an interest in Equation (6f).
The first line indicates that a structure of type $(S,I)\in\mathbb{S}$ changes
its composition upon either the infection of one of its susceptible members,
be it an infection within the structure itself or outside of it (within a
structure of the opposite type or in the general population), or upon the
removal of one of its infected members. Simultaneously, a structure of type
$(S,I+1)$ transforms into a structure of type $(S,I)$ whenever one of its
infected members recovers, while a structure of type $(S+1,I-1)$ becomes of
type $(S,I)$ upon infection of a susceptible member. In particular, this
result shows that under the assumptions of Theorem 2.2, in the large
population limit, we may neglect the natural correlation between structures
caused by the fact that infected individuals belong to two structures at once.
This allows to obtain a stronger model reduction than in Theorem 2.1, in the
sense that the model reduces to a finite-dimensional ODE-system instead of a
measure-valued equation.

Before detailing the proofs of Theorems 2.1 and 2.2, let us examine the latter
from a numerical point of view.

###  2.2 Numerical assessment of the limiting dynamical system

The aim of this section is first to portray that the proposed large population
limit, under the form of dynamical system eqs. 6b to 6f, is in good accordance
with the original stochastic model for large population sizes. This
secondarily leads to some practical comments on the implementation of the
dynamical system. Finally, a comparison with another reduced model for
epidemics with two layers of mixing will be established, namely with an Edge-
Based Compartmental Model (EBCM) in the line of work of Volz _et al._ [29].

####  2.2.1 Implementation of the dynamical system and illustration of Theorem
2.2

Figure 1: Household and workplace size distributions $\pi^{H}$ (left) and
$\pi^{W}$ (right) used in simulations.

Let us start by illustrating the result of Theorem 2.2 through numerical
simulations. Using Gillespie’s algorithm, we have performed fifty simulations
of the epidemic within a population of $K=10000$ individuals, where $\pi^{H}$
and $\pi^{W}$ are roughly inspired by the French household and workplace
distributions as observed in 2018 by Insee [4]. These distributions are
represented in Figure 1. Two sets of epidemic parameters have been considered,
leading to either $R_{0}=2.5$ or $R_{0}=1.2$. In both cases, the majority of
contaminations take place at the local level. Indeed, for the first scenario
with $R_{0}=2.5$, 42% and 18% of infections occur within households and
workplaces, respectively. For the second scenario, mean-field infections are
even less common, with 40% of within-household and 40% of within-workplace
infections. Further, the epidemic is started by infecting either $10$ or $100$
individuals chosen uniformly at random at time $0$ ($\varepsilon=0.001$ or
$\varepsilon=0.01$, respectively). For each simulation, we have followed the
evolution of the proportion of susceptible and infected individuals within the
population over time. When $\varepsilon=0.001$, time subsequently is shifted
so that time $0$ corresponds to the first moment when at least 5‰ of the
population are infectious.

The simulation outcomes are presented in Figure 2. For each choice of
parameters, the solutions $s$ and $i$ of dynamical system eqs. 6b to 6f with
initial condition given by (7) are plotted on the same graph. As previously
for $\varepsilon=0.001$, time is subsequently shifted to ensure that
$i(0)=0.005$, in order to synchronize the simulations among themselves and
with their deterministic limit. As expected, one observes good accordance of
the stochastic simulations and the deterministic functions $(s,i)$.

Figure 2: Comparison of the stochastic model with its large population
approximation given by dynamical system eqs. 6b to 6f. Household and workplace
distributions are those of Figure 1. Two sets of epidemic parameters were
considered, namely
$(\beta_{G},\lambda_{H},\lambda_{W},\gamma)=(0.125,1.5,0.00115,0.125)$ and
$(\beta_{G},\lambda_{H},\lambda_{W},\gamma)=(0.03,0.05,0.0015,0.125)$ for the
left and right column respectively ($R_{0}=2.5$ and $R_{0}=1.2$). The initial
conditions are either $\varepsilon=0.001$ in Panels (a) and (b), or
$\varepsilon=0.01$ in (c) and (d). For each of these scenarios, Gillespie’s
algorithm was used to simulate $50$ trajectories of the stochastic model
defined in Proposition 1.1 in a population of $K=10000$ individuals (faint
lines; the same 50 realizations of the population structure were used for each
subplot). For Panels (a) and (b), only trajectories reaching a threshold of
$0.005$ infected were kept, and time was shifted so that time $0$ corresponds
to the moment when this threshold is reached. Finally, the deterministic
solution $(s,i)$ of eqs. 6b to 6f is represented for each scenario (thick
lines). For Panels (a) and (b), the same time shifting procedure as for
simulations is applied.

Before proceeding further, let us briefly emphasize a few aspects of the
implementation of the proposed deterministic model. A potential drawback of
dynamical system eqs. 6b to 6f consists in its large dimension. Indeed, it
holds that $\\#\mathbb{S}=n_{\max}(n_{\max}+1)/2-1$. The number of equations
of dynamical system eqs. 6b to 6f is hence given by $2+2\\#\mathbb{S}$, and
thus of order $O(n_{\max}^{2})$. However, this fast-growing number of
equations actually is manageable, as it is possible to implement the dynamical
system in an automated way, in the sense that each equation does not need to
be written one-by-one by the programmer. We refer to Appendix A.1 for details.

Nevertheless, the large dimension of the dynamical system of interest raises
the question whether it is numerically speaking interesting to actually use it
for numerical explorations when compared to stochastic simulations using
Gillespie’s algorithm. Indeed, if the computation time needed to solve the
dynamical system is not advantageous when compared to stochastic simulations,
one may prefer to use the latter as it contains more information, such as the
fluctuations around the large population limit. In order to address this
question, we have compared the average time needed to either solve once
dynamical system eqs. 6b to 6f, or to simulate one trajectory of the
stochastic model using Gillespie’s algorithm, for different choices of
epidemic parameters. Notice here that in practice, for stochastic simulations,
it is often necessary to compute several individual trajectories in order to
obtain the general behavior of the epidemic. However, as it is possible to
execute these simulations in parallel, comparison to one individual simulation
seemed the most pertinent.

One hundred independent runs of each script were performed and their
computation time assessed. For each run, the runtime of a reference function
was also measured, in order to take into account fluctuations of computation
time which occur among runs of the same script. The procedure and results are
detailed in Appendix A.2. In summary, solving the reduced model is up to one
order of magnitude faster than performing one stochastic simulation for values
of $R_{0}>1$ that are not too close to the critical case $R_{0}=1$. This shows
that the reduced model is pertinent for numerical exploration.

####  2.2.2 Comparison to Edge-Based Compartmental Models

One may notice that the population structure, as described in Section 1.1, can
be regarded as a modification of the well-studied configuration model. Indeed,
our network of household- and workplace-contacts may be seen as a two-layer
graph, where each layer corresponds to a random graph generated as described
in [19, 21], that we shall call clique configuration model (CCM) hereafter.
This random graph model generalizes configuration models to include small,
totally connected sub-graphs referred to as _cliques_. It then is possible to
derive an EBCM for our household-workplace model by reasoning as in [29].
Details are provided in Appendix B.1.

Edge-based compartmental models on CCM variants have been known to be in good
accordance with simulations of the corresponding stochastic epidemic models in
terms of proportions of susceptible and infected individuals within the
population, under the assumption of a very small initial proportion of
infected. In our case, we have confronted the EBCM with dynamical system eqs.
6b to 6f, as well as simulated trajectories of our stochastic model. Results
are reported on Figure B1. The parameters considered are the same as for Panel
(a) of Figure 2, except that several values of the initial proportion of
infected $\varepsilon$ have been explored. As expected, for very small values
of $\varepsilon$, the EBCM and dynamical system eqs. 6b to 6f both yield the
correct epidemic dynamics, whereas for larger values of $\varepsilon$, the
EBCM does not fit the simulated epidemic trajectories. We refer to Appendix
B.1 for details.

Finally, one may notice that the EBCM is of higher dimension than dynamical
system eqs. 6b to 6f. Proceeding like before, we obtain that the number of
equations of the EBCM is of order $O(n_{\max}^{3})$. As expected, this has a
negative impact on computation time, as briefly illustrated in Appendix B.2,
arguing against the applicability of this EBCM for numerical explorations.

To conclude, in the particular case of the household-workplace model studied
in this article, the EBCM seems to be equivalent to the large population
approximation described by dynamical system eqs. 6b to 6f, under the condition
that the initial proportion of infected is very small. However, considering
both the computational cost of its higher dimension and the loss of accuracy
for more general initial conditions of the EBCM, the large population
approximation given by dynamical system eqs. 6b to 6f seems more pertinent in
the case of the epidemic model under consideration.

##  3 Proofs

This section is devoted to establishing Theorems 2.1 and 2.2. As we will see,
the proof of Theorem 2.1 naturally has the intrinsic difficulty of all
convergence results for measure-valued processes, but it allows us to obtain a
deterministic prediction of the dynamics of the structure type distributions
during the course of an epidemic. At this level, the limiting object is rich,
allowing it to convey detailed information on the distribution of remaining
infectious periods within structures. This however comes at the cost of an
infinite-dimensional limiting object, which motivates the interest in trying
to further reduce its dimension by adopting a coarser population description.
In the case where $\nu$ is the exponential distribution, Theorem 2.2 shows
that this actually is possible, the final reduced model taking the form of
dynamical system eqs. 6b to 6f. As mentioned previously, the existence of an
asymptotically exact, closed, finite-dimensional ODE-system capturing the
epidemic dynamics was not obvious from the beginning. Indeed, in order to
obtain this result, we need to eliminate from the state space $E$ the
continuous component $\tau$, which needs to be transformed back into the
simpler distinction between infected and removed states for each member of the
structure. As we shall see in Subsection 3.2, this will demand some effort as
removal events are not explicit in the model formalism adopted in Section 1.2.
This difficulty will become apparent in the proof of forthcoming Proposition
3.13.

###  3.1 Proof of Theorem 2.1

Let us start with the proof of Theorem 2.1. In particular, throughout this
section, $\nu$ is an arbitrary probability distribution on $\mathbb{R}_{+}$.
The proof follows a classical scheme, establishing tightness of
$(\zeta^{K})_{K\geq 1}$, whose limiting values are shown to satisfy Equation
(4). Uniqueness of the solutions of this equation given the initial condition
then ensures the desired convergence result. In particular, the proof follows
[10] and [27].

####  3.1.1 Uniqueness and continuity of the solution of Equation (4).

We are first going to establish a uniqueness result for the solutions of
Equation (4). Notice that we do not need to prove existence of solutions in
this section, as forthcoming Proposition 3.7 constructs such solutions as
limiting values of $(\zeta^{K})_{K\geq 1}.$

Let us start with a technical lemma, whose proof we present for sake of
completeness.

######  Lemma 3.1.

Let $f\in\mathcal{C}_{b}(E,\mathbb{R})$. There exists a sequence
$(f_{k})_{k\geq 1}$ taking values in $\mathcal{C}^{1}_{b}(E,\mathbb{R})$ such
that $f_{k}$ converges simply to $f$ and $\sup_{k\geq 1}\left\lVert
f_{k}\right\rVert_{\infty}\leq\left\lVert f\right\rVert_{\infty}$.

###### Proof.

Consider a mollifier $\psi$, _i.e._
$\psi\in\mathcal{C}^{\infty}(\mathbb{R}^{n_{\max}})$ is compactly supported,
its mass $\int_{\mathbb{R}^{n_{\max}}}\psi(x)dx$ equals $1$, and for $k\geq
1$, the function $\psi_{k}:x\mapsto k^{n_{\max}}\psi(x/k)$ satisfies
$\lim_{k\to\infty}\psi_{k}=\delta_{0}$ in the sense of distributions.

Consider $f\in\mathcal{C}_{b}(E,\mathbb{R})$. Define the sequence
$(f_{k})_{k\geq 1}:E\to\mathbb{R}$ as follows:

$f_{k}:(n,s,\tau)\mapsto f(n,s,\cdot)*\psi_{k}(\tau).$

Then, for any $x\in E$, by definition of $(\psi_{k})_{k\geq 1}$ it holds that
$\lim_{k\to\infty}f_{k}(x)=f(x)$. Further, as $\psi_{k}$ is of integral $1$
for any $k$, it is obvious that for any $k$, $\left\lVert
f_{k}\right\rVert_{\infty}\leq\left\lVert f\right\rVert_{\infty}$. Finally, it
also follows from the usual properties of convolution that for any $k\geq 1$,
$f_{k}$ is smooth with respect to its last variable and all the corresponding
partial derivatives are compactly supported, hence
$f_{k}\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$ for any $k\geq 1$. ∎

We may now turn to the main result of this paragraph.

######  Proposition 3.2.

Let $\eta_{\star}\in\mathfrak{M}_{1}$. Then Equation (4) admits at most one
measure-valued solution $\eta$ which belongs to
$\mathcal{C}(\mathbb{R}_{+},\mathfrak{M}_{1})$, such that
$\eta_{0}=\eta_{\star}$.

From now on, define $n^{H}=\langle\eta^{H}_{0},\mathbf{n}\rangle$,
$s^{H}_{t}=\langle\eta^{H}_{t},\mathbf{s}\rangle$,
$i^{H}_{t}=\langle\eta^{H}_{t},\mathbf{i}\rangle$, and
$s^{W}_{t}=\langle\eta^{W}_{t},\mathbf{s}\rangle$. Let us establish the
proposition.

###### Proof.

First, notice that it follows immediately from Equation (4) that
$\langle\eta^{X}_{T},1\rangle=\langle\eta^{X}_{0},1\rangle$ for
$X\in\\{H,W\\}$, thus $\eta_{0}\in\mathfrak{M}_{1}$ implies that for any
$T\geq 0$, $\eta_{T}\in\mathfrak{M}_{1}$.

Let us show that any solution $\eta$ of Equation (4) belongs to
$\mathcal{C}(\mathbb{R}_{+},\mathfrak{M}_{1})$. In order to do so, it is
enough to show that both
$\eta^{H}\in\mathcal{C}(\mathbb{R}_{+},\mathcal{M}_{1}(E))$ and
$\eta^{W}\in\mathcal{C}(\mathbb{R}_{+},\mathcal{M}_{1}(E))$. We are going to
detail the proof for $\eta^{H}$ only, as $\eta^{W}$ can be handled in the same
way.

Let $T\geq 0$ and $g\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$ such that
$\left\lVert g\right\rVert_{\infty}\leq 1$. Consider the function defined by

$f_{t}(x)=g(\Psi(x,T,t))\;\;\forall(t,x)\in\mathbb{R}\times E,$

and recall from Proposition 1.1 that
$f^{\mathcal{I}}_{t}(x)=\langle\nu,f_{t}(\mathfrak{j}(x,\cdot))\rangle$. Then
by definition, $f_{T}(x)=g(x)$. It follows from the assumption
$g\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$ that
$f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$. The advantage of
this construction is that $t\mapsto\Psi(x,T,t)$ corresponds to a reversal of
time, which cancels out the deterministic dynamics described by the
differential operator $\mathcal{A}$. Indeed, letting $x=(n,s,\tau)$, a brief
computation shows that

$\displaystyle\partial_{t}f_{t}(x)=\sum_{k=1}^{n-s}\delta_{\tau_{k}}g(\Psi(x,T,t))\text{
and }\partial_{\tau_{k}}f_{t}(x)=\delta_{\tau_{k}}g(\Psi(x,T,t)),$

which yields that $\mathcal{A}f_{t}(x)=0$ for all $(t,x)\in\mathbb{R}\times
E$. Using the fact that
$\langle\eta^{H}_{T},g\rangle=\langle\eta^{H}_{T},f_{T}\rangle$, it follows
from Equation (4) that

$\displaystyle\langle\eta^{H}_{T},g\rangle$
$\displaystyle=\langle\eta^{H}_{0},f_{0}\rangle+\lambda_{H}\int_{0}^{T}\langle\eta^{H}_{t},\mathbf{s}\mathbf{i}\left(f^{\mathcal{I}}_{t}-f_{t}\right)\rangle
dt$
$\displaystyle+\lambda_{W}\int_{0}^{T}\frac{1}{s^{W}_{t}}\langle\eta^{W}_{t},\mathbf{s}\mathbf{i}\rangle\langle\eta^{H}_{t},\mathbf{s}\left(f^{\mathcal{I}}_{t}-f_{t}\right)\rangle
dt+\beta_{G}\int_{0}^{T}\frac{i^{H}_{t}}{n^{H}}\langle\eta^{H}_{t},\mathbf{s}\left(f^{\mathcal{I}}_{t}-f_{t}\right)\rangle
dt.$

Recall that $i^{H}_{t}\leq n_{\max}$ and $\frac{i^{H}_{t}}{n^{H}}\leq 1$ since
for any $x\in E$, $\mathbf{i}(x)\leq\mathbf{n}(x)$. We may notice that the
following inequalities hold, as for any $t$, $\left\lVert
f_{t}\right\rVert_{\infty}\leq 1$:

$\begin{aligned}
&\langle\eta^{H}_{t},\mathbf{s}\mathbf{i}\left(f^{\mathcal{I}}_{t}-f_{t}\right)\rangle\leq
2(n_{\max})^{2},\\\
&\frac{1}{s^{W}_{t}}\langle\eta^{W}_{t},\mathbf{s}\mathbf{i}\rangle\leq
n_{\max},\end{aligned}\quad\quad\begin{aligned}
&\langle\eta^{H}_{t},\mathbf{s}\left(f^{\mathcal{I}}_{t}-f_{t}\right)\rangle\leq
2n_{\max},\\\
&\frac{i^{H}_{t}}{n^{H}}\langle\eta^{H}_{t},\mathbf{s}\left(f^{\mathcal{I}}_{t}-f_{t}\right)\rangle\leq
2n_{\max}.\end{aligned}$ (8)

Let
$C=2n_{\max}\left(\lambda_{H}n_{\max}+\lambda_{W}n_{\max}+\beta_{G}\right)$
and let $\epsilon\in\mathbb{R}$. It then follows from inequalities (8) that

$|\langle\eta^{H}_{T}-\eta^{H}_{T+\epsilon},g\rangle|\leq
C|T-(T+\epsilon)|=C|\epsilon|.$ (9)

Consider now $h\in\mathcal{C}_{b}(E,\mathbb{R})$ such that $\left\lVert
h\right\rVert_{\infty}\leq 1$. Lemma 3.1 ensures that there exists a sequence
$(g_{k})_{k\geq 1}$ taking values in $\mathcal{C}^{1}_{b}(E,\mathbb{R})$ which
converges simply to $h$ and such that $\left\lVert
g_{k}\right\rVert_{\infty}\leq 1$. By dominated convergence, this implies that

$|\langle\eta^{H}_{T}-\eta^{H}_{T+\epsilon},h\rangle|\leq
C|T-(T+\epsilon)|=C|\epsilon|.$

As $E$ is a Polish space, it follows from Proposition A.6.1 of [27] that

$||\eta^{H}_{T}-\eta^{H}_{T+\epsilon}||_{TV}=\sup_{h\in\mathcal{C}_{b}(E,\mathbb{R}):\left\lVert
h\right\rVert_{\infty}\leq
1}|\langle\eta^{H}_{T}-\eta^{H}_{T+\epsilon},h\rangle|\leq C|\epsilon|.$ (10)

As this holds for any $\epsilon$ and any $T\geq 0$, the strong continuity of
$\eta^{H}$ is established.

It remains to establish uniqueness of the solution $\eta$ of Equation (4) with
initial condition $\eta_{0}=\eta_{\star}$. We will once more establish
uniqueness component-wise, and focus on $\eta^{H}$ as $\eta^{W}$ is treated in
a similar fashion.

Let $\eta,\overline{\eta}$ be two solutions of Equation (4) with initial
condition $\eta_{\star}$. Let
$\overline{i}^{H}_{t}=\langle\overline{\eta}^{H}_{t},\mathbf{i}\rangle$ and
define $\overline{s}^{W}_{t}$ in analogous manner. As before, let $T\geq 0$.
Consider again $g\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$ such that $\left\lVert
g\right\rVert_{\infty}\leq 1$, and define $f_{t}$ and $f^{\mathcal{I}}_{t}$ as
previously. Then

$\displaystyle|\langle\eta^{H}_{T}-\overline{\eta}^{H}_{T},g\rangle|$
$\displaystyle\leq\lambda_{H}\int_{0}^{T}\left|\langle\eta^{H}_{t}-\overline{\eta}^{H}_{t},\mathbf{s}\mathbf{i}\left(f^{\mathcal{I}}_{t}-f_{t}\right)\rangle\right|dt$
$\displaystyle+\lambda_{W}\int_{0}^{T}\left|\frac{1}{s^{W}_{t}}\langle\eta^{W}_{t},\mathbf{s}\mathbf{i}\rangle\langle\eta^{H}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle-\frac{1}{\overline{s}^{W}_{t}}\langle\overline{\eta}^{W}_{t},\mathbf{s}\mathbf{i}\rangle\langle\overline{\eta}^{H}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle\right|dt$
$\displaystyle+\frac{\beta_{G}}{m_{H}}\int_{0}^{T}\left|i^{H}_{t}\langle\eta^{H}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle-\overline{i}^{H}_{t}\langle\overline{\eta}^{H}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle\right|dt.$

Proceeding similarly as in inequalities (8), we obtain that

$|\langle\eta^{H}_{T}-\overline{\eta}^{H}_{T},g\rangle|\leq\frac{C}{2}\int_{0}^{T}\left|\langle\eta^{H}_{t}-\overline{\eta}^{H}_{t},f^{\mathcal{I}}_{t}-f_{t}\rangle\right|dt\leq\frac{C}{2}\left(\int_{0}^{T}\left|\langle\eta^{H}_{t}-\overline{\eta}^{H}_{t},f^{\mathcal{I}}_{t}\rangle\right|dt+\int_{0}^{T}\left|\langle\eta^{H}_{t}-\overline{\eta}^{H}_{t},f_{t}\rangle\right|dt\right).\\\
$

On the one hand, by definition, $f_{t}\in\mathcal{C}_{b}(E,\mathbb{R})$ and
$\left\lVert f_{t}\right\rVert_{\infty}\leq 1$. Thus, it is obvious that
$\left|\langle\eta^{H}_{t}-\overline{\eta}^{H}_{t},f_{t}\rangle\right|\leq||\eta^{H}_{t}-\overline{\eta}^{H}_{t}||_{TV}$.
On the other hand, it follows from the usual criterion of continuity for
parametric integrals that $f^{\mathcal{I}}_{t}$ is continuous on $E$ and
$\left\lVert f^{\mathcal{I}}_{t}\right\rVert_{\infty}\leq\left\lVert
f_{t}\right\rVert_{\infty}\leq 1$. As a consequence,
$\left|\langle\eta^{H}_{t}-\overline{\eta}^{H}_{t},f^{\mathcal{I}}_{t}\rangle\right|\leq||\eta^{H}_{t}-\overline{\eta}^{H}_{t}||_{TV}$.
Hence

$|\langle\eta^{H}_{T}-\overline{\eta}^{H}_{T},g\rangle|\leq
C\int_{0}^{T}||\eta^{H}_{t}-\overline{\eta}^{H}_{t}||_{TV}dt.$

We then may follow the same steps that allowed to establish Equation (10) from
Equation (9), and obtain that

$||\eta^{H}_{T}-\overline{\eta}^{H}_{T}||_{TV}\leq
C\int_{0}^{T}||\eta^{H}_{t}-\overline{\eta}^{H}_{t}||_{TV}dt.$

Gronwall’s lemma then assures that

$\displaystyle\forall
t\in[0,T],\;\;||\eta^{H}_{t}-\overline{\eta}^{H}_{t}||_{TV}=0.$

Obviously, one obtains the analogous result for $\eta^{W}$ in the same manner.
As $T\geq 0$ is arbitrary, this concludes the proof. ∎

####  3.1.2 Tightness of $(\zeta^{K})_{K\geq 1}$ in
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E),v\otimes
v)\right)$

Let us now turn to the tightness of $(\zeta^{K})_{K\geq 1}$ in
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E),v\otimes
v)\right)$, where $v$ designates the vague topology on $\mathcal{M}_{F}(E)$.
We start by establishing the following preliminary result, whose proof simply
relies on the chain rule.

######  Lemma 3.3.

Let $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$. Then for any
$T\geq t_{0}\geq 0$, for any $x\in E$,

$f\left(T,\Psi(x,T,t_{0})\right)=f\left(t_{0},x\right)+\int_{t_{0}}^{T}\mathcal{A}f\left(t,\Psi(x,t,t_{0})\right)dt.$

###### Proof.

For $x\in E$ and $t_{0}\in\mathbb{R}_{+}$, define
$g_{t_{0},x}:[t_{0},+\infty)\to\mathbb{R},\;\;T\mapsto
f\left(T,\Psi(x,T,t_{0})\right)$. Let us start by noticing that for any
$(t_{0},x)\in\mathbb{R}_{+}\times E$,
$g_{t_{0},x}\in\mathcal{C}^{1}(\mathbb{R}_{+})$. Indeed,
$g_{t_{0},x}=f_{n,s}\circ h_{t_{0},x}$, where

$\displaystyle
f_{n,s}:\mathbb{R}_{+}\times\mathbb{R}^{n_{\max}}\to\mathbb{R},\;\;(u,\mathbf{v})\mapsto
f\left(u,(n,s,\mathbf{v})\right)$ and $\displaystyle
h_{t_{0},x}:\mathbb{R}_{+}\to\mathbb{R}^{1+n_{\max}},\;\;t\mapsto\left(t,\tau-\sum_{k=1}^{n-s}(t-t_{0})e_{k}\right).$

The chain rule and a quick computation of the differentiable of $h_{t_{0},x}$
yields that for every $t\geq t_{0}$,

$g_{t_{0},x}^{\prime}(t)=\partial_{1}f_{n,s}(h_{t_{0},x}(t))-\sum_{k=1}^{n-s}\partial_{k+1}f_{n,s}(h_{t_{0},x}(t)).$

Notice that
$\partial_{1}f_{n,s}(u,\mathbf{v})=\partial_{t}f(u,(n,s,\mathbf{v}))$ and for
$k\geq 2$,
$\partial_{k}f_{n,s}(u,\mathbf{v})=\partial_{\tau_{k-1}}f(u,n,s,\mathbf{v})$.
As a consequence, we have shown that for any $(t_{0},x)\in\mathbb{R}_{+}\times
E$, for every $t\geq t_{0}$, $g_{t_{0},x}$ is differentiable at $t$ and
satisfies

$g_{t_{0},x}^{\prime}(t)=\mathcal{A}f\left(t,\Psi(x,t,t_{0})\right).$

This concludes the proof. ∎

Throughout the section, we use the notation $\mathcal{S}=\\{H,W,G\\}$. Also,
for any $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$, let
$f_{t}(x)=f(t,x)$ for any $(t,x)\in\mathbb{R}_{+}\times E$. Finally, we define
for any continuous bounded function $g:\mathbb{R}_{+}\times E\to\mathbb{R}$,
for any $t\geq 0$ and
$u=(\bm{\theta},k,\ell,\sigma)\in\bigcup_{Y\in\mathcal{S}}U^{Y}_{K}$:

$\displaystyle g^{H}_{t,u}$
$\displaystyle=g(t,\mathfrak{j}(x^{H}_{k}(t),\sigma))-g(t,x^{H}_{k}(t))\quad\text{and}\quad
g^{W}_{t,u}$
$\displaystyle=g(t,\mathfrak{j}(x^{W}_{\ell}(t),\sigma))-g(t,x^{W}_{\ell}(t)).$

######  Proposition 3.4.

Consider $\zeta^{K}$ as introduced in Proposition 1.1. For any
$f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$, $T\geq 0$ and
$X\in\\{H,W\\}$,

$\small\langle\zeta^{X|K}_{T},f_{T}\rangle=\langle\zeta^{X|K}_{0},f_{0}\rangle+\int_{0}^{T}\langle\zeta^{X|K}_{t},\mathcal{A}f_{t}\rangle
dt+\frac{1}{K_{X}}\sum_{Y\in\mathcal{S}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)f^{X}_{t-,u}Q^{Y|K}(dt,du).$

###### Proof.

Let $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$ and
$X\in\\{H,W\\}$. Recall that, by definition, for any bounded function
$g:E\to\mathbb{R}$, for any $T\geq t\geq 0$ and
$u=(\bm{\theta},k,\ell,\sigma)\in\bigcup_{Y\in\mathcal{S}}U^{Y}_{K}$:

$\langle\Delta_{H}(u,T,t),g\rangle=g(\Psi(\mathfrak{j}(x^{H}_{k}(t-),\sigma),T,t))-g(\Psi(x^{H}_{k}(t-),T,t)),$

and $\langle\Delta_{W}(u,T,t),g\rangle$ is defined analogously, by replacing
$H$ by $W$ and $k$ by $\ell$.

From Equation (3), it follows that

$\langle\zeta^{X|K}_{T},f_{T}\rangle=\frac{1}{K_{X}}\sum_{j=1}^{K_{X}}f_{T}\left(\Psi(x^{X}_{j}(0),T,0)\right)+\frac{1}{K_{X}}\sum_{Y\in\mathcal{S}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)\langle\Delta_{X}(u,T,t),f_{T}\rangle
Q^{Y|K}(dt,du).$

Using the result from Lemma 3.3, this becomes:

$\displaystyle\langle\zeta^{X|K}_{T},f_{T}\rangle$
$\displaystyle=\frac{1}{K_{X}}\sum_{j=1}^{K_{X}}\left(f_{0}(x^{X}_{j}(0))+\int_{0}^{T}\mathcal{A}f_{t}(\Psi(x^{X}_{j}(0),t,0))dt\right)$
$\displaystyle+\frac{1}{K_{X}}\sum_{Y\in\mathcal{S}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)\left(\int_{t}^{T}\langle\Delta_{X}(u,z,t),\mathcal{A}f_{z}\rangle
dz\right)Q^{Y|K}(dt,du)$
$\displaystyle+\frac{1}{K_{X}}\sum_{Y\in\mathcal{S}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)f^{X}_{t-,u}Q^{Y|K}(dt,du).$

It follows from the definition of $\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times
E,\mathbb{R})$ that both $f$ and $\mathcal{A}f$ are bounded, hence we may
apply Fubini’s theorem to obtain that

$\displaystyle\langle$
$\displaystyle\zeta^{X|K}_{T},f_{T}\rangle=\frac{1}{K_{X}}\sum_{j=1}^{K_{X}}f_{0}(x^{X}_{j}(0))+\frac{1}{K_{X}}\sum_{Y\in\mathcal{S}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)f^{X}_{t-,u}Q^{Y|K}(dt,du)$
$\displaystyle+\frac{1}{K_{X}}\int_{0}^{T}\left(\sum_{j=1}^{K_{X}}\mathcal{A}f_{z}(\Psi(x^{X}_{j}(0),z,0))+\sum_{Y\in\mathcal{S}}\int_{0}^{z}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)\langle\Delta_{X}(u,z,t),\mathcal{A}f_{z}\rangle
Q^{Y|K}(dt,du)\right)dz.$

The first sum on the right-hand side equals
$\langle\zeta^{X|K}_{0},f_{0}\rangle$. From the second line, one recognizes in
the integrand the definition of
$\langle\zeta^{X|K}_{z},\mathcal{A}f_{z}\rangle$ from Equation (3). This
yields the desired result. ∎

For $Y\in\mathcal{S}$, let

$\widetilde{Q}^{Y|K}(dt,du)=Q^{Y|K}(dt,du)-dt\mu^{Y}_{K}(du)$

be the compensated martingale-measure associated to $Q^{Y|K}$.

It follows that, for $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times
E,\mathbb{R})$ and $X\in\\{H,W\\}$,

$\langle\zeta^{X|K}_{T},f_{T}\rangle=M^{X|K}_{T}(f)+V^{X|K}_{T}(f),$

where we define

$M^{X|K}_{T}(f)=\frac{1}{K_{X}}\sum_{Y\in\mathcal{S}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t-,u)f^{X}_{t-,u}\widetilde{Q}^{Y|K}(dt,du)$

and

$V^{X|K}_{T}(f)=\langle\zeta^{X|K}_{0},f_{0}\rangle+\int_{0}^{T}\langle\zeta^{X|K}_{t},\mathcal{A}f_{t}\rangle
dt+\frac{1}{K_{X}}\sum_{Y\in\mathcal{S}}\int_{0}^{T}\int_{U^{Y}_{K}}\mathcal{I}_{Y}(t,u)f^{X}_{t,u}\mu^{Y}_{K}(du)dt.$

######  Proposition 3.5.

Let $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$ and
$X\in\\{H,W\\}$. Then $(M^{X|K}_{T}(f))_{T\geq 0}$ is a square integrable
martingale. Using the same notations as in Theorem 2.1, its quadratic
variation is given by

$\langle
M^{X|K}(f)\rangle_{T}=\frac{1}{K_{X}}\int_{0}^{T}\langle\zeta^{X|K}_{t},\mathcal{H}^{X}_{t}((f^{2}_{t})^{\mathcal{I}}-2f^{\mathcal{I}}_{t}f_{t}+f^{2}_{t})\rangle
dt,$

where for any $t\geq 0$ and $x\in E$,

$\mathcal{H}^{X}_{t}(x)=\beta_{G}\frac{I^{H}(t)}{N^{H}}\mathbf{s}(x)+\lambda_{X}\mathbf{s}(x)\mathbf{i}(x)+\lambda_{\overline{X}}\frac{\langle\zeta^{\overline{X}|K}_{t},\mathbf{s}\mathbf{i}\rangle}{S^{\overline{X}}_{t}}\mathbf{s}(x).$

###### Proof.

Let $f\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$ and $X\in\\{H,W\\}$. Consider
$M^{X|K}_{T}(f)$, which can be written as

$M^{X|K}_{T}(f)=M^{X|K}_{X,T}(f)+M^{X|K}_{\overline{X},T}(f)+M^{X|K}_{G,T}(f),$

where, for $Y\in\mathcal{S}$ and $T\geq 0$,

$M^{X|K}_{Y,T}(f)=\int_{0}^{T}\int_{U^{Y}_{K}}\frac{1}{K_{X}}\mathcal{I}_{Y}(t-,u)f^{X}_{t,u}\widetilde{Q}^{Y|K}(dt,du).$

Suppose that for any $Y\in\mathcal{S}$,

$\mathbb{E}\left[\int_{0}^{T}\int_{U^{Y}_{K}}\left(\frac{1}{K_{X}}\mathcal{I}_{Y}(t,u)f^{X}_{t,u}\right)^{2}\mu^{Y}_{K}(du)dt\right]<\infty,$

then for all $Y\in\mathcal{S}$, $(M^{X|K}_{Y,T}(f))_{T\geq 0}$ is as square
integrable martingale [17], implying that $(M^{X|K}_{T}(f))_{T\geq 0}$ also is
a square integrable martingale. As $Q^{H|K}$, $Q^{W|K}$ and $Q^{G|K}$ are
independent, it follows that

$\langle M^{X|K}(f)\rangle_{T}=\sum_{Y\in\mathcal{S}}\langle
M^{X|K}_{Y}(f)\rangle_{T}.$

It thus is enough to study $(M^{X|K}_{Y,T}(f))_{T\geq 0}$ for all
$Y\in\mathcal{S}$. In the following, we will detail the necessary computations
in the case $X=H$, the case $X=W$ being similar.

Consider the case $Y=H$. Start by noticing that
$\sum_{\ell=1}^{K_{W}}s^{W}_{\ell}(t)=K_{W}S^{W}(t)$, and that for any
$k\in\llbracket 1,K_{H}\rrbracket$ and $t\in[0,T]$, $s^{H}_{k}(t)$ and
$i^{H}_{k}(t)$ are less then $n_{\max}$. Hence, replacing $S(t)$ by
$K_{W}S^{W}(t)$ in $\mathcal{I}_{H}$,

$\displaystyle\mathbb{E}\left[\langle M^{H|K}_{H}(f)\rangle_{T}\right]$
$\displaystyle=\mathbb{E}\left[\int_{0}^{T}\int_{U^{H}_{K}}\left(\frac{1}{K_{H}}\mathcal{I}_{H}(t,u)f^{H}_{t,u}\right)^{2}\mu^{H}_{K}(du)dt\right]$
$\displaystyle=\mathbb{E}\left[\int_{0}^{T}\frac{1}{K_{H}^{2}}\sum_{k=1}^{K_{H}}\lambda_{H}s^{H}_{k}(t)i(\tau^{H}_{k}(t))\langle\nu,\left(f_{t}(\mathfrak{j}(x^{H}_{k}(t),\cdot))-f_{t}(x^{H}_{k}(t))\right)^{2}\rangle
dt\right]$
$\displaystyle\leq\frac{1}{K_{H}}\lambda_{H}(n_{\max})^{2}4\left\lVert
f\right\rVert_{\infty}^{2}T.$

Since further $K_{H}\geq K/n_{\max}$ and $\left\lVert
f\right\rVert_{\infty}^{2}<\infty$, we obtain that

$\mathbb{E}\left[\langle
M^{H|K}_{H}(f)\rangle_{T}\right]\leq\frac{4}{K}\lambda_{X}(n_{\max})^{3}\left\lVert
f\right\rVert_{\infty}^{2}T<\infty.$ (11)

Thus $(M^{H|K}_{H,t}(f))_{t\geq 0}$ is a square integrable martingale whose
quadratic variations is given by

$\langle
M^{H|K}_{H}(f)\rangle_{T}=\int_{0}^{T}\int_{U^{H}_{K}}\left(\frac{1}{K_{H}}\mathcal{I}_{H}(t,u)f^{H}_{t,u}\right)^{2}\mu^{H}_{K}(du)dt=\frac{\lambda_{H}}{K_{H}}\int_{0}^{T}\langle\zeta^{H|K}_{t},\mathbf{s}\mathbf{i}((f^{2}_{t})^{\mathcal{I}}-2f^{\mathcal{I}}_{t}f_{t}+f^{2}_{t})\rangle
dt,$

using the computations from Equation (3.1.2).

Similarly, $(M^{H|K}_{W,T}(f))_{T\geq 0}$ is a square integrable martingale of
quadratic variation

$\langle
M^{H|K}_{W}(f)\rangle_{T}=\lambda_{W}\int_{0}^{T}\frac{\langle\zeta^{W|K}_{t},\mathbf{s}\mathbf{i}\rangle}{S^{W}(t)}\langle\zeta^{H|K}_{t},\mathbf{s}((f^{2}_{t})^{\mathcal{I}}-2f^{\mathcal{I}}_{t}f_{t}+f^{2}_{t})\rangle
dt.$

Further, for the case $Y=G$, let us use the equalities $S(t)=K_{H}S^{H}(t)$
and $I(t)=K_{H}I^{H}(t)$. As $I(t)/K=I^{H}(t)/N^{H}\leq 1$ almost surely since
$\zeta^{H|K}_{t}\in\mathcal{M}_{1}(E)$, we obtain that

$\mathbb{E}[\left[\langle
M^{H|K}_{G}(f)\rangle_{T}\right]=\mathbb{E}\left[\int_{0}^{T}\int_{U^{G}_{K}}\left(\frac{1}{K_{H}}\mathcal{I}_{G}(t,u)f^{H}_{t,u}\right)^{2}\mu^{G}_{K}(du)dt\right]\leq\frac{4}{K}\beta_{G}(n_{\max})^{2}\left\lVert
f\right\rVert_{\infty}^{2}T.$

As before, $(M^{H|K}_{G,T}(f))_{T\geq 0}$ thus is a square integrable
martingale of quadratic variation given by

$\langle
M^{H|K}_{G}(f)\rangle_{T}=\frac{1}{K_{H}}\beta_{G}\int_{0}^{T}\frac{I^{H}(t)}{N^{H}}\langle\zeta^{H|K}_{t},\mathbf{s}((f^{2}_{t})^{\mathcal{I}}-2f^{\mathcal{I}}_{t}f_{t}+f^{2}_{t})\rangle
dt.$

This yields the desired result for $(M^{H|K}_{T}(f))_{T\geq 0}$, and
proceeding similarly for $(M^{W|K}_{T}(f))_{T\geq 0}$ concludes the proof. ∎

We are now ready to establish the desired result on the tightness of
$(\zeta^{K})_{K\geq 1}$.

######  Proposition 3.6.

Under the assumptions of Theorem 2.1, the sequence $(\zeta^{K})_{K\geq 1}$ is
tight in
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E),v\otimes
v)\right)$.

The proof relies on the fact that in order to establish tightness of
$(\zeta^{K})_{K\geq 1}$, it is enough to show that for any $X\in\\{H,W\\}$,
$\left(\langle\zeta^{X|K},f\rangle\right)_{K\geq 1}$ is tight for a large
enough set of test functions $f$ [24]. This in turn is ensured using the
Aldous [1] and Rebolledo [14] criteria, whose application is straightforward
thanks to the upper bounds established in the previous proof.

###### Proof.

Once more, we will proceed component-wise and show that
$\left(\zeta^{H|K}\right)_{K\geq 1}$ and $\left(\zeta^{W|K}\right)_{K\geq 1}$
are both tight in
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E),v)\right)$.

Let us focus on $\left(\zeta^{H|K}\right)_{K\geq 1}$. According to Theorem 2.1
of [24], it is sufficient to show that for any function $f$ belonging to a
dense subset of

$\mathcal{C}_{0}(E,\mathbb{R})=\left\\{f:E\to\mathbb{R}\text{ continuous
s.t.}\lim_{\left\lVert x\right\rVert_{\infty}\to\infty}|f(x)|=0\right\\},$

the sequence $\left(\langle\zeta^{H|K},f\rangle\right)_{K\geq 1}$ is tight in
$\mathbb{D}(\mathbb{R}_{+},\mathbb{R})$. Notice that by density of
$C^{\infty}_{K}(\mathbb{R}^{n_{\max}})$ in $C_{0}(\mathbb{R}^{n_{\max}})$
endowed with the uniform norm, it follows that
$\mathcal{C}_{0}(E,\mathbb{R})\cap\mathcal{C}^{1}_{b}(E,\mathbb{R})$ is also
dense in $\mathcal{C}_{0}(E,\mathbb{R})$ endowed with the uniform norm. Thus,
let us consider
$f\in\mathcal{C}_{0}(E,\mathbb{R})\cap\mathcal{C}^{1}_{b}(E,\mathbb{R})$.

According to the Aldous [1] and Rebolledo [14] criteria, in order to prove the
tightness of $(\langle\zeta^{H|K},f\rangle)_{K\geq 1}$, it is enough to show
that:

  1. (i)

For any $t$ belonging to a dense subset $\mathcal{T}$ of $\mathbb{R}^{+}$,
both $\left(\langle M^{H|K}(f)\rangle_{t}\right)_{K\geq 0}$ and
$(V^{H|K}_{t}(f))_{K\geq 0}$ are tight in $\mathbb{R}$.

  2. (ii)

For any $T\geq 0$, for any $\epsilon,\alpha>0$, there exist $\delta>0$ and
$K_{0}\in\mathbb{N}$ such that for any two sequences of stopping times
$(S_{K})_{K\geq 1}$ and $(T_{K})_{K\geq 1}$ satisfying $S_{K}\leq T_{K}\leq T$
for all integers $K$,

$\displaystyle\sup_{K\geq K_{0}}\mathbb{P}\left(|\langle
M^{H|K}(f)\rangle_{S_{K}}-\langle
M^{H|K}(f)\rangle_{T_{K}}|\geq\alpha,T_{K}\leq S_{K}+\delta\right)$
$\displaystyle\leq\epsilon$ $\displaystyle\text{and }\sup_{K\geq
K_{0}}\mathbb{P}\left(|V^{H|K}_{S_{K}}(f)-V^{H|K}_{T_{K}}(f)|\geq\alpha,T_{K}\leq
S_{K}+\delta\right)$ $\displaystyle\leq\epsilon.$

Notice that, in order to establish (i), it is enough to show that for any
$t\geq 0$,

$\sup_{K\geq 1}\mathbb{E}\left[|\langle
M^{H|K}(f)\rangle_{t}|\right]<\infty\text{ and }\sup_{K\geq
1}\mathbb{E}\left[\left|V^{H|K}_{t}(f)\right|\right]<\infty.$

Recalling that
$C=2n_{\max}\left(\lambda_{H}n_{\max}+\lambda_{W}n_{\max}+\beta_{G}\right)$,
it follows from Equation (11) that

$\mathbb{E}[|\langle
M^{H|K}(f)\rangle_{t}|]\leq\frac{1}{K}2n_{\max}C\left\lVert
f\right\rVert_{\infty}^{2}t.$

Similar computations yield that

$\mathbb{E}[|V^{H|K}_{t}(f)|]\leq\left\lVert
f\right\rVert_{\infty}+\left\lVert\mathcal{A}f\right\rVert_{\infty}+C\left\lVert
f\right\rVert_{\infty}t.$

As $f\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$, this implies that (i) holds.

It remains to check (ii). Let $\epsilon,\alpha>0$, and consider two sequences
of stopping times $(S_{K})_{K\geq 1}$ and $(T_{K})_{K\geq 1}$ satisfying
$S_{K}\leq T_{K}\leq T$ for all integers $K$. As previously, using Equation
(11), we obtain the following upper bound:

$\displaystyle E\left[|\langle M^{H|K}(f)\rangle_{S_{K}}-\langle
M^{H|K}(f)\rangle_{T_{K}}|\Big{|}T_{K}\leq S_{K}+\delta\right]$
$\displaystyle\leq\mathbb{E}\left[\int_{S_{K}}^{T_{K}}dt\Big{|}T_{K}\leq
S_{K}+\delta\right]\frac{1}{K}2n_{\max}C\left\lVert
f\right\rVert_{\infty}^{2}$
$\displaystyle\leq\frac{\delta}{K}2n_{\max}C\left\lVert
f\right\rVert_{\infty}^{2}.$

Hence, using conditional Markov’s inequality,

$\mathbb{P}\left(|\langle M^{H|K}(f)\rangle_{S_{K}}-\langle
M^{H|K}(f)\rangle_{T_{K}}|\geq\alpha,T_{K}\leq
S_{K}+\delta\right)\leq\frac{\delta}{\alpha K}2n_{\max}C\left\lVert
f\right\rVert_{\infty}^{2}.$ (12)

Proceeding similarly, we also obtain that

$\mathbb{P}\left(|V^{H|K}_{S_{K}}(f)-V^{H|K}_{T_{K}}(f)|\geq\alpha,T_{K}\leq
S_{K}+\delta\right)\leq\frac{\delta}{\alpha}\left(\left\lVert\mathcal{A}f\right\rVert_{\infty}+C\left\lVert
f\right\rVert_{\infty}\right).$ (13)

Equations (12) and (13) imply the existence of $\delta$ and $K_{0}$ such that
(ii) is satisfied. Naturally, $\zeta^{W|K}$ can be handled analogously. This
concludes the proof. ∎

####  3.1.3 Identification of the limiting values of $(\zeta^{K})_{K\geq 1}$

The tightness of $(\zeta^{K})_{K\geq 1}$ in the space
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E),v\otimes
v)\right)$ ensures that from any subsequence of $(\zeta^{K})_{K\geq 1}$, one
may extract a subsubsequence which converges in this space. The limits of
these subsubsequences may be characterized as follows.

######  Proposition 3.7.

Under the assumptions of Theorem 2.1, all limiting values of
$(\zeta^{K})_{K\geq 1}$ in
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E),v\otimes
v)\right)$ are continuous solutions of Equation (4).

###### Proof.

Consider a subsequence $(\zeta^{\varphi(K)})_{K\geq 1}$ of $(\zeta^{K})_{K\geq
1}$ which converges in the Skorokhod space
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E),v\otimes
v)\right)$, and let $\zeta$ be its limit.

Let us start by showing that almost surely,
$\zeta\in\mathcal{C}(\mathbb{R}_{+},\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E))$.
As usually, we will proceed component-wise, and only exhibit the reasoning for
$\zeta^{H}$ as $\zeta^{W}$ can be treated accordingly. Let $T\geq 0$ and let
$f\in\mathcal{C}_{b}(E,\mathbb{R})$ such that $\left\lVert
f\right\rVert_{\infty}\leq 1$. It follows from Proposition 1.1 that

$\left|\langle\zeta^{H|\varphi(K)}_{T},f\rangle-\langle\zeta^{H|\varphi(K)}_{T-},f\rangle\right|\leq\frac{2n_{\max}}{\varphi(K)}\;\;\text{almost
surely}.$

As a consequence, letting $K\to\infty$ and following Proposition A.6.1 of
[27],

$||\zeta^{H}_{T}-\zeta^{H}_{T-}||_{TV}=\sup_{f\in\mathcal{C}_{b}(E,\mathbb{R}):\left\lVert
f\right\rVert_{\infty}\leq
1}\left|\langle\zeta^{H}_{T},f\rangle-\langle\zeta^{H}_{T-},f\rangle\right|=0\;\;\text{almost
surely}.$

It remains to show that $\zeta$ satisfies Equation (4). Let $T\geq 0$ and
$f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$, and consider the
application $\psi^{H}_{T}$ defined by

$\displaystyle\psi^{H}_{T}(\zeta)$
$\displaystyle=\langle\zeta^{H}_{T},f_{T}\rangle-\langle\zeta^{H}_{0},f_{0}\rangle-\int_{0}^{T}\langle\zeta^{H}_{t},\mathcal{A}f_{t}\rangle
dt-\int_{0}^{T}\langle\zeta^{H}_{t},\lambda_{H}\mathbf{s}\mathbf{i}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt$ (14)
$\displaystyle-\int_{0}^{T}\frac{1}{\langle\zeta^{W}_{t},\mathbf{s}\rangle}\langle\zeta^{W}_{t},\lambda_{W}\mathbf{s}\mathbf{i}\rangle\langle\zeta^{H}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt-\int_{0}^{T}\frac{\langle\zeta^{H}_{t},\mathbf{i}\rangle}{\langle\zeta^{H}_{0},\mathbf{n}\rangle}\langle\zeta^{H}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt.$

Start by noticing that
$\psi^{H}_{T}(\zeta^{\varphi(K)})=M^{H|\varphi(K)}_{T}(f)$, as
$K_{H}S^{H}(t)=K_{W}S^{W}(t)$. Using Jensen’s inequality, it follows from
Equation (11) that

$\mathbb{E}[|\psi^{H}_{T}(\zeta^{\varphi(K)})|]^{2}\leq\mathbb{E}[|\psi^{H}_{T}(\zeta^{\varphi(K)})|^{2}]=\mathbb{E}[\langle
M^{H|K}(f)\rangle_{t}]\leq\frac{1}{K}2n_{\max}C\left\lVert
f\right\rVert_{\infty}^{2}T\xrightarrow[K\to\infty]{}0.$

As $\zeta$ is almost surely continuous and
$f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$, $\psi^{H}_{T}$
is continuous at $\zeta$. Hence $(\psi^{H}_{T}(\zeta^{\varphi(K)}))_{K\geq 1}$
converges in law to $\psi^{H}_{T}(\zeta)$. According to Theorem 3.5 of [6], it
then is enough to proof that $(\psi^{H}_{T}(\zeta^{\varphi(K)}))_{K\geq 1}$ is
uniformly integrable to obtain that its expectation converges to the
expectation of $\psi^{H}_{T}(\zeta)$. In our case, uniform integrability is
easily assured as the sequence $(\psi^{H}_{T}(\zeta^{\varphi(K)}))_{K\geq 1}$
is bounded. Indeed, using the fact that for all $T\geq 0$,
$\zeta^{\varphi(K)}_{T}\in\mathfrak{M}_{1}$, we obtain from Equation (14) that

$\left|\psi^{H}_{T}(\zeta^{\varphi(K)})\right|\leq\left((2+C)\left\lVert
f\right\rVert_{\infty}+\left\lVert\mathcal{A}f\right\rVert_{\infty}\right)T.$

We may now conclude that

$\mathbb{E}[|\psi^{H}_{T}(\zeta)|]=\lim_{K\to\infty}\mathbb{E}[|\psi^{H}_{T}(\zeta^{\varphi(K)})|]=0,$

which yields the desired result. ∎

####  3.1.4 Proof of Theorem 2.1

The previous results are sufficient to establish Theorem 2.1. Indeed, it
follows from propositions 3.6 and 3.7 that from every subsequence of
$(\zeta^{K})_{K\geq 1}$, one may extract a subsubsequence converging in
$\mathbb{D}\left(\mathbb{R}_{+},(\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E),v\otimes
v)\right)$ to a solution of Equation (4) which almost surely belongs to
$\mathcal{C}(\mathbb{R}_{+},\mathcal{M}_{F}(E)\times\mathcal{M}_{F}(E))$. As
by assumption, $\zeta^{K}_{0}$ converges in law to
$\eta_{0}\in\mathfrak{M}_{1}$, Proposition 3.2 implies that all of these
subsubsequences converge to the same limit $\eta$, which is the unique
solution of Equation (4) with initial condition $\eta_{0}$. As
$\eta_{0}\in\mathfrak{M}_{1}$, Proposition 3.2 further ensures that
$\eta\in\mathbb{D}(\mathbb{R}_{+},\mathfrak{M}_{1})$. This establishes the
convergence of $(\zeta^{K})_{K\geq 1}$ to $\eta$ in
$\mathbb{D}(\mathbb{R}_{+},\mathfrak{M}_{1})$, as desired.

###  3.2 Proof of Theorem 2.2

This section is devoted to extracting dynamical system eqs. 6b to 6f from the
measure-valued integral equation (4), under the assumption that $\nu$ is the
exponential distribution of parameter $\gamma$. From now on, we suppose that
this assumption is satisfied, and further that $\eta_{0}=\eta_{0,\varepsilon}$
as given by (5).

####  3.2.1 Preliminary study of the dynamical system

Before establishing Theorem 2.2 itself, let us start by showing that the
solution of dynamical system eqs. 6b to 6f endowed with initial condition (7)
is well defined, in the sense that there exists a unique global solution.
Existence actually also follows from the proofs of the forthcoming
subsections, since they construct a solution to the dynamical system. However,
as we will use Cauchy-Lipschitz theorem to establish uniqueness, existence is
obtained at the same time.

For this section, let us rewrite dynamical system eqs. 6b to 6f as follows, in
order to emphasize the associated Cauchy problem. Recall that the dynamical
system is of dimension $d=2+2\\#\mathbb{S}=n_{\max}(n_{\max}+1)$.

Let $y\in\mathcal{C}^{1}(\mathbb{R}_{+},\mathbb{R}^{d})$ and
$f:\mathbb{R}^{d}\to\mathbb{R}^{d}$ be defined such that dynamical system eqs.
6b to 6f amounts to

$y^{\prime}(t)=f(y(t))\;\;\forall t\geq 0.$ (15)

The components of $y$ (and resp. $f$) will be called $s$, $i$ and
$n^{X}_{(S,I)}$ (resp. $f_{s}$, $f_{i}$ and $f_{X;(S,I)}$) for $X\in\\{H,W\\}$
and $(S,I)\in\mathbb{S}$, in order to simplify their identification with the
unknowns of the corresponding dynamical system. More precisely, consider the
applications

$\tau_{X}(y)=-\frac{\lambda_{X}}{m_{X}}\sum_{(S,I)\in\mathbb{S}}SI\;n^{X}_{(S,I)}\text{
for }X\in\\{H,W\\},\text{ and }\tau_{G}(y)=\beta_{G}i.$

Then $f:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is defined as follows, for any
$y=(s,i,n^{X}_{(S,I)}:X\in\\{H,W\\},(S,I)\in\mathbb{S})\in\mathbb{R}^{d}$:

$f_{s}(y)=-(\tau_{H}(y)+\tau_{W}(y)+\tau_{G}(y)s)\;\text{and}\;f_{i}(y)=-f_{s}(y)-\gamma
i,\\\ $

while for all $X\in\\{H,W\\}$ and $(S,I)\in\mathbb{S}$

$\displaystyle f_{X;(S,I)}(y)=$
$\displaystyle-\left[\left(\lambda_{X}I+\frac{\tau_{\overline{X}}(y)}{s}+\tau_{G}(y)\right)S-\gamma
I\right]n^{X}_{(S,I)}+\gamma(I+1)n^{X}_{(S,I+1)}\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}}$
$\displaystyle+\left(\lambda_{X}(I-1)+\frac{\tau_{\overline{X}}(y)}{s}+\tau_{G}(y)\right)(S+1)n^{X}_{(S+1,I-1)}\mathbf{1}_{\left\\{I\geq
1\right\\}}.$

Also, notice that there are some natural constraints that we expect the
solution of dynamical system eqs. 6b to 6f to satisfy. Clearly, $s$, $i$ and
$n^{X}_{(S,I)}$ should belong to $[0,1]$. Also, as the population is
partitioned into susceptible, infected and removed individuals, it follows
that $s+i\leq 1$. Similarly, as all individuals belong to exactly one
household and one workplace, and as $n^{X}_{(S,I)}$ corresponds to the
proportion of structures of type $X$ which contain $S$ susceptible and $I$
infected individuals, we expect that for $X\in\\{H,W\\}$,

$\sum_{(S,I)\in\mathbb{S}}n^{X}_{(S,I)}\leq 1,\text{ and
}\sum_{(S,I)\in\mathbb{S}}Sn^{X}_{(S,I)}\leq m_{X}s.$ (16)

We thus define the following set $V\subset\mathbb{R}^{d}$, which formalizes
these constraints:

$\small V=\left\\{y\in[0,1]^{d}:s+i\leq
1,\sum_{(S,I)\in\mathbb{S}}n^{X}_{(S,I)}\leq 1\text{ and
}m_{X}s-\sum_{(S,I)\in\mathbb{S}}Sn^{X}_{(S,I)}\geq 0\;\;\forall
X\in\\{H,W\\}\right\\}.$

######  Proposition 3.8.

Let $y^{*}\in V$. Then the following assertions hold:

  1. (i)

Suppose that there exists a solution $y$ of the Cauchy problem (15) with
initial condition $y(0)=y^{*}$. Then $y(t)\in V$ for any $t\geq 0$ for which
$y$ is well defined.

  2. (ii)

There exists a unique such solution $y$, which further is global.

  3. (iii)

In particular, for any $\varepsilon>0$, the dynamical system eqs. 6b to 6f
endowed with initial condition (7) admits a unique solution on
$\mathbb{R}_{+}$.

The proof of this proposition is postponed to Appendix C.

####  3.2.2 Absolute continuity

Consider any $n\in\llbracket 1,n_{\max}\rrbracket$ and $s\in\llbracket
0,n\rrbracket$. The aim of this section is to establish that conditionally on
a structure being of size $n$ and containing $s$ susceptibles, the
distribution of the remaining infectious periods of its $n-s$ members who have
contracted the disease at some previous time, is absolutely continuous with
respect to the Lebesgue measure on $\mathbb{R}^{n-s}$. More precisely, we will
see that at any time $t\geq 0$, there exists a function
$\rho^{X,n,s}_{t}:\mathbb{R}^{n-s}\to\mathbb{R}_{+}$ such that for any non-
negative measurable function $f:E\to\mathbb{R}$, for any $X\in\\{H,W\\}$ and
$t\geq 0$,

$\langle\eta^{X}_{t},\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(\bullet)=n\\\
\mathbf{s}(\bullet)=s\end{subarray}\right\\}}f(\bullet)\rangle=\int_{\mathbb{R}^{n-s}}f(\tau)\rho^{X,n,s}(\tau)d\tau.$

In the following, for any $n\geq 1$, let $0_{n}$ be the zero of
$\mathbb{R}^{n}$. Further, for $k\in\llbracket 1,n\rrbracket$ and
$\tau\in\mathbb{R}^{n}$, let $\tau=(\tau_{1},\dots,\tau_{n})$ and
$\tau_{1,k}=(\tau_{1},\dots,\tau_{k})$. Finally, $(e_{k})_{1\leq k\leq
n_{\max}}$ is the canonical basis of $\mathbb{R}^{n_{\max}}$.

######  Lemma 3.9.

For each $X\in\\{H,W\\}$, there exists a family of measurable functions
indexed by $n\in\llbracket 1,n_{\max}\rrbracket$ and $s\in\llbracket
0,n\rrbracket$ such that

$\displaystyle\forall s\in\llbracket
0,n-1\rrbracket,\;\rho^{X,n,s}:\;\;\mathbb{R}_{+}\times\mathbb{R}^{n-s}$
$\displaystyle\to\mathbb{R}_{+}$ $\displaystyle(t,\tau)$
$\displaystyle\mapsto\rho^{X,n,s}_{t}(\tau),$

which verifies that for any $f:E\to\mathbb{R}$ non-negative measurable
function,

$\langle\eta^{X}_{t},f\rangle=\sum_{n=1}^{n_{\max}}\left(f(n,n,0_{n_{\max}})\rho^{X,n,n}_{t}+\sum_{s=0}^{n-1}\int_{\mathbb{R}^{n-s}}f(n,s,\sum_{k=1}^{n-s}\tau_{k}e_{k})\rho^{X,n,s}_{t}(\tau)d\tau\right).$

The proof relies on studying the distribution $\eta^{X}$ conditionally on
structures being of given size and number of susceptibles. Absolute continuity
with regard to the Lebesgue measure is then established by domination with
another, absolutely continuous measure. We refer to Appendix C for details.

Notice that in the case $n=s$, $\rho^{X,n,n}$ depends on $t$ alone, whereas
for $s<n$, $\rho^{X,n,s}$ depends both on $t$ and on
$\tau\in\mathbb{R}^{n-s}$. Throughout the following, with some abuse of
notation, this distinction will mostly be implicit. The adaptation to the case
$s=n$ is generally straightforward. For instance, forthcoming Equation (17)
reduces to

$\frac{d}{dt}\rho^{X,n,n}_{t}=-\Lambda^{X}_{t}n\rho^{X,n,n}_{t},$

where

$\Lambda^{X}_{t}=\lambda_{\overline{X}}\frac{1}{s^{\overline{X}}_{t}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{t}}{n^{H}}.$

When integrating over $\tau$, the case $n=s$ becomes

$\int_{\mathbb{R}^{0}}\rho^{X,n,n}_{t}(\tau)d\tau=\rho^{X,n,n}_{t}.$

Similarly, symbolically, for any $\tau\in\mathbb{R}^{n-s}$, $i(\tau)=0$
whenever $n-s=0$. This makes sense, as in a structure of size $n$ containing
$n$ susceptibles, there are no infected. This leads for instance to

$\displaystyle\int_{\mathbb{R}^{0}}\mathbf{1}_{\left\\{i(\tau)=0\right\\}}\rho^{X,n,n}_{t}(\tau)d\tau=\rho^{X,n,n}_{t}.$

When the adaptation is less clear, the case $n=s$ will be treated separately.

######  Proposition 3.10.

The family $(\rho^{X,n,s})_{1\leq n\leq n_{\max},0\leq s\leq n}$ is a weak
solution to the following system of partial differential equations: for any
$(n,s)\in\llbracket 1,n_{\max}\rrbracket\times\llbracket 0,n\rrbracket$, for
any $\tau\in\mathbb{R}^{n-s}$ and $t\geq 0$,

$\displaystyle\partial_{t}\rho^{X,n,s}_{t}(\tau)-$
$\displaystyle\sum_{k=1}^{n-s}\partial_{\tau_{k}}\rho^{X,n,s}_{t}(\tau)=-\lambda_{X}s(i(\tau)+\Lambda^{X}_{t})\rho^{X,n,s}_{t}(\tau)$
(17) $\displaystyle+\mathbf{1}_{\left\\{s+1\leq
n\right\\}}\lambda_{X}(s+1)i(\tau_{1,n-s-1})\rho^{X,n,s+1}_{t}(\tau_{1,n-s-1})\mathbf{1}_{\left\\{\tau_{n-s}>0\right\\}}\gamma\mathrm{e}^{-\gamma\tau_{n-s}}$
$\displaystyle+\mathbf{1}_{\left\\{s+1\leq
n\right\\}}\Lambda^{X}_{t}(s+1)\rho^{X,n,s+1}_{t}(\tau_{1,n-s-1})\mathbf{1}_{\left\\{\tau_{n-s}>0\right\\}}\gamma\mathrm{e}^{-\gamma\tau_{n-s}},$

with initial condition given by

$\rho^{X,n,n}_{0}=\pi^{X}_{n}\varepsilon^{n}\text{ and }\forall
s<n,\tau\in\mathbb{R}^{n-s},\;\rho^{X,n,s}_{0}(\tau)=\pi^{X}_{n}\binom{n}{s}\varepsilon^{s}(1-\varepsilon)^{n-s}\mathbf{1}_{\left\\{\tau\in(\mathbb{R}_{+}^{*})^{n-s}\right\\}}\mathrm{e}^{-\gamma\sum_{k=1}^{n-s}\tau_{k}}.$

###### Proof.

Let $(n,s)\in\llbracket 1,n_{\max}\rrbracket\times\llbracket 0,n\rrbracket$.
In the case $s=n$, Equation (17) follows directly from Equation (4) applied to
$f(x)=\mathbf{1}_{\left\\{\mathbf{n}(x)=n,\mathbf{s}(x)=s\right\\}}$ for all
$x\in E$ as $\mathcal{A}f=f^{\mathcal{I}}=0$.

Suppose now that $s<n$. Consider $f$ measurable on $E$ such that $f(x)=0$ if
$(\mathbf{n}(x),\mathbf{s}(x))\neq(n,s)$ and $f(x)=f_{n,s}(\tau(x)_{1,n-s})$
otherwise, where $f_{n,s}\in\mathcal{C}^{\infty}(\mathbb{R}^{n-s})$ is
compactly supported. In particular, for any $\sigma>0$ and $x\in E$ such that
$\mathbf{s}(x)\geq 1$, notice that $f(\mathfrak{j}(x,\sigma))$ is zero, unless
$s+1\leq n$ and $(\mathbf{n}(x),\mathbf{s}(x))=(n,s+1)$ in which case it is
equal to $f_{n,s}(\tau(x)_{1,n-s-1}+\sigma e_{n-s})$. Recalling that
$f^{\mathcal{I}}(x)=\langle\nu,f(\mathfrak{j}(x,\cdot)\rangle$, injecting $f$
into Equation (4) and differentiating with regard to $t$ thus yields

$\displaystyle\frac{d}{dt}\int_{\mathbb{R}^{n-s}}f_{n,s}(\tau)\rho^{X,n,s}_{t}(\tau)d\tau=-\sum_{k=1}^{n-s}\int_{\mathbb{R}^{n-s}}\partial_{\tau_{k}}f_{n,s}(\tau)\rho^{X,n,s}_{t}(\tau)d\tau$
$\displaystyle-\int_{\mathbb{R}^{n-s}}(\lambda_{X}si(\tau)+\Lambda^{X}_{t}s)f_{n,s}(\tau)\rho^{X,n,s}_{t}(\tau)d\tau+\mathbf{1}_{\left\\{s=n-1\right\\}}\Lambda^{X}_{t}n\left(\int_{0}^{+\infty}f_{n,n-1}(\sigma
e_{1})\gamma\mathrm{e}^{-\gamma\sigma}d\sigma\right)\rho^{X,n,n}_{t}$
$\displaystyle+\mathbf{1}_{\left\\{s+1<n\right\\}}\int_{\mathbb{R}^{n-s-1}}(\lambda_{X}i(u)+\Lambda^{X}_{t})(s+1)\left(\int_{0}^{+\infty}f_{n,s}(u+\sigma
e_{n-s})\gamma\mathrm{e}^{-\gamma\sigma}d\sigma\right)\rho^{X,n,s+1}_{t}(u)du.$

An integration by parts leads to the following equality, as $f_{n,s}$ is of
compact support which implies that the boundary terms of the integration by
parts vanish:

$-\sum_{k=1}^{n-s}\int_{\mathbb{R}^{n-s}}\partial_{\tau_{k}}f_{n,s}(\tau)\rho^{X,n,s}_{t}(\tau)d\tau=\sum_{k=1}^{n-s}\int_{\mathbb{R}^{n-s}}f_{n,s}(\tau)\partial_{\tau_{k}}\rho^{X,n,s}_{t}(\tau)d\tau.$

Injecting this into the previous computation yields that $\rho^{X,n,s}$ is a
weak solution to Equation (17).

Finally, initial conditions $\rho^{X,n,s}_{0}$ for $n\in\llbracket
1,n_{\max}\rrbracket$ and $s\in\llbracket 0,n\rrbracket$ follow directly from
(5). ∎

####  3.2.3 Proof of Theorem 2.2

Before proceeding, let us introduce the following notation. For
$f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$ and $T\geq t\geq
0$, define

$f_{T,t}:x\in E\mapsto f(T,\Psi(x,T,t))\text{ and }f_{T,t}^{\mathcal{I}}:x\in
E\mapsto\langle\nu,f_{T,t}(\mathfrak{j}(x,\cdot))\rangle.$ (18)

We can now state an intermediate result, which is similar in spirit to
Proposition 3.4.

######  Proposition 3.11.

Let $\eta$ be the unique solution in
$\mathcal{C}(\mathbb{R}_{+},\mathfrak{M}_{1})$ of Equation (4). Then for any
$f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E,\mathbb{R})$, for
$X\in\\{H,W\\}$,

$\displaystyle\langle\eta^{X}_{T},f_{T}\rangle$
$\displaystyle=\langle\eta^{X}_{0},f_{T,0}\rangle+\lambda_{X}\int_{0}^{T}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}(f^{\mathcal{I}}_{T,t}-f_{T,t})\rangle
dt$ (19)
$\displaystyle+\lambda_{\overline{X}}\int_{0}^{T}\frac{1}{s^{\overline{X}}_{t}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle\langle\eta^{X}_{t},\mathbf{s}(f^{\mathcal{I}}_{T,t}-f_{T,t})\rangle
dt+\beta_{G}\int_{0}^{T}\frac{i^{H}_{t}}{n^{H}}\langle\eta^{X}_{t},\mathbf{s}(f^{\mathcal{I}}_{T,t}-f_{T,t})\rangle
dt.$

###### Proof.

The proof follows the exact same lines as the proof of Proposition 3.4,
showing that for any $f\in\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E)$,
Equation (19) leads to Equation (4) using Lemma 3.3. ∎

Throughout this section, let $h$ be the Heaviside step function, _i.e._
$h(z)=\mathbf{1}_{\left\\{z>0\right\\}}$ for any real number $z$. Also,
consider a mollifier $\psi$ on $\mathbb{R}$ in the same sense as in the Proof
of Lemma 3.1. For $\alpha>0$, define $\psi_{\alpha}(z)=\alpha^{-1}\psi(\alpha
z)$ for all $z\in\mathbb{R}$. Let us introduce, for any $\alpha>0$, the
function $h_{\alpha}:\mathbb{R}\to[0,1]$ defined by

$\forall z\in\mathbb{R},\quad h_{\alpha}(z)=h*\psi_{\alpha}(z).$

As $\psi_{\alpha}$ converges to $\delta_{0}$ in the sense of distributions
when $\alpha$ tends to zero, the families of functions
$(h_{\alpha})_{\alpha>0}$ and $(\psi_{\alpha})_{\alpha>0}$ serve as smooth
approximations of $h$ and $\delta_{0}$ respectively. This will be useful in
the proofs of this section.

This allows us to establish the following proposition, which serves as a
starting point of the proof of Theorem 2.2.

######  Proposition 3.12.

Under the assumptions of Theorem 2.2, it holds that for any $X\in\\{H,W\\}$,
$(s^{X}_{t})_{t\geq 0}$ and $(i^{X}_{t})_{t\geq 0}$ satisfy

$\displaystyle\frac{d}{dt}s^{X}_{t}$
$\displaystyle=-n^{X}\left(\frac{\lambda_{X}}{n^{X}}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}\rangle+\frac{\lambda_{\overline{X}}}{n^{\overline{X}}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{t}}{n^{H}}\frac{s^{X}_{t}}{n^{X}}\right),$
(20) $\displaystyle\frac{d}{dt}i^{X}_{t}$
$\displaystyle=-\frac{d}{dt}s^{X}_{t}+\gamma i^{X}_{t}.$

Further,

$s^{X}_{0}=(1-\varepsilon)n^{X}\;\;\text{and}\;\;i^{X}_{0}=\varepsilon n^{X}.$
(21)

###### Proof.

Notice that $\mathbf{s}\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$ is such that
$\mathcal{A}\mathbf{s}=0$ and $\mathbf{s}^{\mathcal{I}}(x)-\mathbf{s}(x)=-1$
for all $x\in E$. It thus follows immediately from Equation (4) that, for any
$X\in\\{H,W\\}$,

$s^{X}_{T}=s^{X}_{0}-\int_{0}^{T}\left(\lambda_{X}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}\rangle+\lambda_{\overline{X}}\frac{s^{X}_{t}}{s^{\overline{X}}_{t}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{t}}{n^{H}}s^{X}_{t}\right)dt.$

Further, since $\eta\in\mathcal{C}(\mathbb{R}_{+},\mathfrak{M}_{1})$, it
follows that for any $t\geq 0$, $(\zeta^{K}_{t})_{K\geq 1}$ converges in law
to $\eta_{t}$. As $\mathbf{s}$ is continuous and bounded on $E$, this implies
that
$\langle\zeta^{X|K}_{t},\mathbf{s}\rangle\to\langle\eta^{X}_{t},\mathbf{s}\rangle$
when $K$ tends to infinity. The analogous result holds for $\mathbf{n}$. Hence
Lemma 1.2 ensures that, for any $t\geq 0$,
$s^{X}_{t}/s^{\overline{X}}_{t}=n^{X}/n^{\overline{X}}$. In other words,

$s^{X}_{T}=s^{X}_{0}-n^{X}\int_{0}^{T}\left(\frac{\lambda_{X}}{n^{X}}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}\rangle+\frac{\lambda_{\overline{X}}}{n^{\overline{X}}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{t}}{n^{H}}\frac{s^{X}_{t}}{n^{X}}\right)dt.$

As $\eta\in\mathcal{C}(\mathbb{R}_{+},\mathfrak{M}_{1})$, and $\mathbf{s}$ and
$\mathbf{i}$ are bounded measurable functions, it follows that the integrand
is continuous with regard to $t$. Thus, the first line of Equation (20) comes
from the fundamental theorem of calculus.

Recall that $\eta_{0}=\eta_{0,\varepsilon}$ as defined in Equation (5). In
particular, we now have $n^{X}=m_{X}$, hence

$s^{X}_{0}=\langle\eta^{X}_{0,\varepsilon},\mathbf{s}\rangle=\sum_{n=1}^{n_{\max}}\pi^{X}_{n}\sum_{s=0}^{n}s\binom{n}{s}(1-\varepsilon)^{s}\varepsilon^{n-s}=(1-\varepsilon)\sum_{n=1}^{n_{\max}}n\pi^{X}_{n}=(1-\varepsilon)n^{X}.$
(22)

This yields the first part of Equation (21).

It remains to take an interest in $i^{X}_{t}$. As $\mathbf{i}$ does not belong
to $\mathcal{C}^{1}_{b}(E,\mathbb{R})$ we cannot proceed in the same way.
Remember that $\mathbf{i}(x)=\sum_{j=1}^{n_{\max}}h(\tau_{j})$ for
$x=(n,s,\tau)\in E$. Let us introduce, for any $\alpha>0$, the function
$\mathbf{i}_{\alpha}:E\to\llbracket 1,n_{\max}\rrbracket$ defined by

$\forall x=(n,s,\tau)\in
E,\quad\mathbf{i}_{\alpha}(x)=\sum_{j=1}^{n_{\max}}h_{\alpha}(\tau_{j}).$

Fix $\alpha>0$. As $\psi_{\alpha}^{\prime}$ is bounded, it follows from the
usual properties of convolution that
$h_{\alpha}\in\mathcal{C}^{1}_{b}(\mathbb{R})$ and hence
$\mathbf{i}_{\alpha}\in\mathcal{C}^{1}_{b}(E,\mathbb{R})$. Thus, we may apply
Proposition 3.11 and obtain that

$\displaystyle\langle\eta^{X}_{T},\mathbf{i}_{\alpha}\rangle$
$\displaystyle=\langle\eta^{X}_{0},(\mathbf{i}_{\alpha})_{T,0}\rangle+\lambda_{X}\int_{0}^{T}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}((\mathbf{i}_{\alpha})^{\mathcal{I}}_{T,t}-(\mathbf{i}_{\alpha})_{T,t})\rangle
dt$
$\displaystyle+\lambda_{\overline{X}}\int_{0}^{T}\frac{1}{s^{\overline{X}}_{t}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle\langle\eta^{X}_{t},\mathbf{s}((\mathbf{i}_{\alpha})^{\mathcal{I}}_{T,t}-(\mathbf{i}_{\alpha})_{T,t})\rangle
dt+\beta_{G}\int_{0}^{T}\frac{i^{H}_{t}}{n^{H}}\langle\eta^{X}_{t},\mathbf{s}((\mathbf{i}_{\alpha})^{\mathcal{I}}_{T,t}-(\mathbf{i}_{\alpha})_{T,t})\rangle
dt.$

Notice that $\left\lVert\mathbf{i}_{\alpha}\right\rVert_{\infty}=n_{\max}$ for
any $\alpha>0$, as well as
$\left\lVert\mathbf{s}\right\rVert_{\infty}=n_{\max}$,
$\left\lVert\mathbf{i}\right\rVert_{\infty}=n_{\max}$,
$\left\lVert\mathbf{s}\mathbf{i}\right\rVert_{\infty}\leq n_{\max}^{2}$ and
$\langle\eta^{\overline{X}},\mathbf{s}\mathbf{i}\rangle/s^{\overline{X}}_{t}\leq
n_{\max}$ for any $t\geq 0$. As further $\eta_{t}\in\mathfrak{M}_{1}$ for any
$t\geq 0$, we may let $\alpha$ go to zero and use dominated convergence to
obtain that

$\displaystyle\langle\eta^{X}_{T},\mathbf{i}\rangle$
$\displaystyle=\langle\eta^{X}_{0},\mathbf{i}_{T,0}\rangle+\lambda_{X}\int_{0}^{T}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}(\mathbf{i}^{\mathcal{I}}_{T,t}-\mathbf{i}_{T,t})\rangle
dt$ (23)
$\displaystyle+\lambda_{\overline{X}}\int_{0}^{T}\frac{1}{s^{\overline{X}}_{t}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle\langle\eta^{X}_{t},\mathbf{s}(\mathbf{i}^{\mathcal{I}}_{T,t}-\mathbf{i}_{T,t})\rangle
dt+\beta_{G}\int_{0}^{T}\frac{i^{H}_{t}}{n^{H}}\langle\eta^{X}_{t},\mathbf{s}(\mathbf{i}^{\mathcal{I}}_{T,t}-\mathbf{i}_{T,t})\rangle
dt.$

Further, for any $\sigma\geq 0$, $T\geq t\geq 0$ and $x\in E$,

$i(\Psi(\mathfrak{j}(x,\sigma),T,t))-i(\Psi(x,T,t))=\mathbf{1}_{\left\\{\sigma>(T-t)\right\\}},$

hence

$\mathbf{i}^{\mathcal{I}}_{T,t}(x)-\mathbf{i}_{T,t}(x)=\nu([T-t,\infty))=\mathrm{e}^{-\gamma(T-t)}.$

Injecting this into Equation (23) and using as before that
$s^{X}_{t}/s^{\overline{X}}_{t}=n^{X}/n^{\overline{X}}$ yields

$i^{X}_{T}=\langle\eta^{X}_{0},\mathbf{i}_{T,0}\rangle+n^{X}\mathrm{e}^{-\gamma
T}\int_{0}^{T}\mathrm{e}^{\gamma
t}\left(\frac{\lambda_{X}}{n^{X}}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}\rangle+\frac{\lambda_{\overline{X}}}{n^{\overline{X}}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{t}}{n^{H}}\frac{s^{X}_{t}}{n^{X}}\right)dt.$
(24)

As $\eta_{0}=\eta_{0,\varepsilon}$, we may compute the first term of the
right-hand side of this equation and obtain that

$\displaystyle\langle\eta^{X}_{0},\mathbf{i}_{T,0}\rangle=\langle\eta^{X}_{0,\varepsilon},\mathbf{i}_{T,0}\rangle$
(25)
$\displaystyle=\sum_{n=1}^{n_{\max}}\pi^{X}_{n}\sum_{s=0}^{n}\binom{n}{s}(1-\varepsilon)^{s}\varepsilon^{n-s}\int_{\mathbb{R}^{n_{\max}}}\sum_{j=1}^{n_{\max}}\mathbf{1}_{\left\\{\tau_{j}-T\mathbf{1}_{\left\\{j\leq
n-s\right\\}}>0\right\\}}\nu^{\otimes(n-s)}\otimes\delta_{0}^{\otimes(n_{\max}-(n-s))}(d\tau_{1},\dots,d\tau_{n_{\max}})$
$\displaystyle=\mathrm{e}^{-\gamma
T}\sum_{n=1}^{n_{\max}}\pi^{X}_{n}\sum_{s=0}^{n}\binom{n}{s}(1-\varepsilon)^{s}\varepsilon^{n-s}(n-s)=\mathrm{e}^{-\gamma
T}\varepsilon n^{X}.$

Using the continuity of the integrand in Equation (24), we may now
differentiate it with regard to $T$:

$\displaystyle\frac{d}{dT}i^{X}_{T}$
$\displaystyle=n^{X}\left(\frac{\lambda_{X}}{n^{X}}\langle\eta^{X}_{T},\mathbf{s}\mathbf{i}\rangle+\frac{\lambda_{\overline{X}}}{n^{\overline{X}}}\langle\eta^{\overline{X}}_{T},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{T}}{n^{H}}\frac{s^{X}_{T}}{n^{X}}\right)$
$\displaystyle-n^{X}\gamma\left(\mathrm{e}^{-\gamma
T}\varepsilon+\mathrm{e}^{-\gamma T}\int_{0}^{T}\mathrm{e}^{\gamma
t}\left(\frac{\lambda_{X}}{n^{X}}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}\rangle+\frac{\lambda_{\overline{X}}}{n^{\overline{X}}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{t}}{n^{H}}\frac{s^{X}_{t}}{n^{X}}\right)dt\right)$
$\displaystyle=-\frac{d}{dT}s^{X}_{T}-\gamma i^{X}_{T}.$

We thus have recovered the second line of Equation (20). Finally, the second
half of Equation (21) is obtained by a computation analogous to (22). This
concludes the proof. ∎

For $(S,I)\in\mathbb{S}$, define the function $f^{(S,I)}:E\to\\{0,1\\}$ by

$f^{(S,I)}(x)=\mathbf{1}_{\left\\{\mathbf{s}(x)=S,\,\mathbf{i}(x)=I\right\\}}.$

For $t\geq 0$, let $n^{X}_{(S,I)}(t)=\langle\eta^{X}_{t},f^{(S,I)}\rangle$,
which defines a continuous function on $\mathbb{R}_{+}$ as $f^{(S,I)}$ is
bounded and measurable and
$\eta\in\mathcal{C}(\mathbb{R}_{+},\mathfrak{M}_{1})$. In words, this
corresponds to the proportion of structures of type $X$ which contain exactly
$S$ susceptible and $I$ infected individuals. Notice that

$\\{x\in E:\mathbf{s}(x)\mathbf{i}(x)>0\\}\subset\\{x\in
E:(\mathbf{s}(x),\mathbf{i}(x))\in\mathbb{S}\\}.$

We may thus rewrite the first line of Equation (20) as follows:

$\small\frac{d}{dt}s^{X}_{t}=-n^{X}\left(\frac{\lambda_{X}}{n^{X}}\sum_{(S,I)\in\mathbb{S}}SIn^{X}_{(S,I)}(t)+\frac{\lambda_{\overline{X}}}{n^{\overline{X}}}\sum_{(S,I)\in\mathbb{S}}SIn^{\overline{X}}_{(S,I)}(t)+\beta_{G}\frac{i^{H}_{t}}{n^{H}}\frac{s^{X}_{t}}{n^{X}}\right).$

Similarly, it holds that

$\Lambda^{X}_{t}=\frac{\lambda_{\overline{X}}}{s^{\overline{X}}_{t}}\langle\eta^{\overline{X}}_{t},\mathbf{s}\mathbf{i}\rangle+\beta_{G}\frac{i^{H}_{t}}{n^{H}}=\frac{\lambda_{\overline{X}}}{s^{\overline{X}}_{t}}\sum_{(S,I)\in\mathbb{S}}SIn^{\overline{X}}_{(S,I)}(t)+\beta_{G}\frac{i^{H}_{t}}{n^{H}},$

which may also be written in terms of the notations of Equation eqs. 6b to 6f:

$\Lambda^{X}_{t}=\left(\frac{s^{\overline{X}}_{t}}{n^{X}}\right)^{-1}\tau_{\overline{X}}(t)+\beta_{G}\frac{i^{H}_{t}}{n^{H}}.$

This motivates a closer study of the functions $n^{X}_{(S,I)}$ for
$(S,I)\in\mathbb{S}$.

######  Proposition 3.13.

Let $X\in\\{H,W\\}$ and $(S,I)\in\mathbb{S}$. Under the assumptions of Theorem
2.2, it holds that

$\displaystyle\frac{d}{dt}n^{X}_{(S,I)}(t)$
$\displaystyle=\gamma\left((I+1)n^{X}_{(S,I+1)}(t)\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}}-In^{X}_{(S,I)}(t)\right)$
(26)
$\displaystyle+\lambda_{X}\left((S+1)(I-1)n^{X}_{(S+1,I-1)}(t)\mathbf{1}_{\left\\{I\geq
1\right\\}}-SIn^{X}_{(S,I)}(t)\right)$
$\displaystyle+\Lambda^{X}_{t}\left((S+1)n^{X}_{(S+1,I-1)}(t)\mathbf{1}_{\left\\{I\geq
1\right\\}}-Sn^{X}_{(S,I)}(t)\right).$

Further

$n^{X}_{(S,I)}(0)=\binom{S+I}{I}\pi^{X}_{S+I}(1-\varepsilon)^{S}\varepsilon^{I}.$
(27)

###### Proof.

Let $X\in\\{H,W\\}$ and $(S,I)\in\mathbb{S}$. It follows from Proposition 3.9
that

$n^{X}_{(S,I)}(t)=\langle\eta^{X}_{t},f^{(S,I)}\rangle=\sum_{n=S+I}^{n_{\max}}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\rho^{(n,S)}_{t}(\tau)d\tau,$

where in the second equality, the sum over $n$ starts at $S+I$ as it is
impossible to observe $S$ susceptible and $I$ infected members in a group of
size $n<S+I$. Using Proposition 3.10, we may deduce that

$\frac{d}{dt}n^{X}_{(S,I)}(t)=\sum_{n=S+I}^{n_{\max}}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\partial_{t}\rho^{(n,S)}_{t}(\tau)d\tau.$

Equation (17) yields that

$\displaystyle\frac{d}{dt}$ $\displaystyle
n^{X}_{(S,I)}(t)=\sum_{n=S+I}^{n_{\max}}\sum_{k=1}^{n-S}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\partial_{\tau_{k}}\rho^{X,n,S}_{t}(\tau)d\tau$
(28)
$\displaystyle-\sum_{n=S+I}^{n_{\max}}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\left(\lambda_{X}SI-\Lambda^{X}_{t}S\right)\rho^{X,n,S}_{t}(\tau)$
$\displaystyle+\sum_{n=S+I}^{n_{\max}}\mathbf{1}_{\left\\{S+1\leq
n\right\\}}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\lambda_{X}(S+1)i(\tau_{1,n-S-1})\rho^{X,n,S+1}_{t}(\tau_{1,n-S-1})\mathbf{1}_{\left\\{\tau_{n-S}>0\right\\}}\gamma\mathrm{e}^{-\gamma\tau_{n-S}}d\tau$
$\displaystyle+\sum_{n=S+I}^{n_{\max}}\mathbf{1}_{\left\\{S+1\leq
n\right\\}}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\Lambda^{X}_{t}(S+1)\rho^{X,n,S+1}_{t}(\tau_{1,n-S-1})\mathbf{1}_{\left\\{\tau_{n-S}>0\right\\}}\gamma\mathrm{e}^{-\gamma\tau_{n-S}}d\tau.$

Let us start by noticing that the second line of Equation (28) may be written
as

$-\left(\lambda_{X}SI-\Lambda^{X}_{t}S\right)\sum_{n=S+I}^{n_{\max}}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\rho^{X,n,S}_{t}(\tau)=-\left(\lambda_{X}SI-\Lambda^{X}_{t}S\right)n^{X}_{(S,I)}(t).$

Next, let us pay attention to the third and fourth line of Equation (28). For
any $\tau\in\mathbb{R}^{n-S}$, it appears that on the event $\tau_{n-S}>0$,
$i(\tau)=I$ if and only if $i(\tau_{1,n-S-1})=I-1$. The third line thus
becomes

$\displaystyle\sum_{n=S+I}^{n_{\max}}$
$\displaystyle\mathbf{1}_{\left\\{S+1\leq
n\right\\}}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau_{1,n-S-1})=I-1\right\\}}\lambda_{X}(S+1)(I-1)\rho^{X,n,S+1}_{t}(\tau_{1,n-S-1})\mathbf{1}_{\left\\{\tau_{n-S}>0\right\\}}\gamma\mathrm{e}^{-\gamma\tau_{n-S}}d\tau$
$\displaystyle=\lambda_{X}(S+1)(I-1)\sum_{n=S+I}^{n_{\max}}\mathbf{1}_{\left\\{S+1\leq
n\right\\}}\int_{\mathbb{R}^{n-S-1}}\mathbf{1}_{\left\\{i(\tau^{\prime})=I-1\right\\}}\rho^{X,n,S+1}_{t}(\tau^{\prime})d\tau^{\prime}$
$\displaystyle=\mathbf{1}_{\left\\{I\geq
1\right\\}}\lambda_{X}(S+1)(I-1)\sum_{n=S+I}^{n_{\max}}\int_{\mathbb{R}^{n-S-1}}\mathbf{1}_{\left\\{i(\tau^{\prime})=I-1\right\\}}\rho^{X,n,S+1}_{t}(\tau^{\prime})d\tau^{\prime}$
$\displaystyle=\mathbf{1}_{\left\\{I\geq
1\right\\}}\lambda_{X}(S+1)(I-1)n^{X}_{(S+1,I-1)}(t).$

In this computation, the indicator function of the event $\\{I\geq 1\\}$
appears as for any $(n,s)$, $\rho^{(n,s)}$ never charges the space
$\\{i(\tau)<0\\}$. As further $S+1\leq n$ whenever $S+I\leq n$ and $I\geq 1$,
the indicator function of the set $\\{S+1\leq n\\}$ can be removed. Similarly,
the fourth line of Equation (28) becomes

$\mathbf{1}_{\left\\{I\geq
1\right\\}}\Lambda^{X}_{t}(S+1)n^{X}_{(S+1,I-1)}(t).$

It remains to focus on the first line of Equation (28). In order to do so, we
will use another expression of $f^{(S,I)}$. For two integers $j\leq n$, let
$\mathbb{B}(n,j)$ be the set of unordered subsets of $j$ elements chosen in
$\llbracket 1,n\rrbracket$. It then holds that for any $(n^{\prime},s,\tau)\in
E$,

$f^{(S,I)}(n^{\prime},s,\tau)=\sum_{n=S+I}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}n^{\prime}=n\\\
s=S\end{subarray}\right\\}}\sum_{\mathbf{v}\in\mathbb{B}(n-S,I)}\prod_{j\in\mathbf{v}}h(\tau_{j})\prod_{j\in\llbracket
1,n-S\rrbracket\setminus\mathbf{v}}(1-h(\tau_{j})),$

where $h$ designates as before the Heaviside step function. The idea behind
this expression is that a structure of type $(n,s,\tau)$ contains $S$
susceptible and $I$ infected members if and only if $s=S$, and further exactly
$I$ out of the $n-S$ first components of $\tau$ (which describe the states of
members of the structure having contracted the disease at some point, _i.e._
infected and removed members) are positive. In other words, there exists at
most one element $\mathbf{v}\in\mathbb{B}(n-S,I)$ for which the term in the
sum is not equal to zero, in which case the set $\mathbf{v}$ corresponds to
the indexes of infectious members, while $\llbracket
1,n-S\rrbracket\setminus\mathbf{v}$ is the set of removed members.

In the following, for any integers $1\leq j\leq n$ and any
$\mathbf{v}\in\mathbb{B}(n,j)$, consider the sets

$A_{\mathbf{v}}=\\{\tau\in\mathbb{R}^{n}:\tau_{k}>0\;\forall
k\in\mathbf{v},\,\tau_{k}\leq 0\;\forall k\notin\mathbf{v}\\}$

and for any $k\in\mathbf{v}$,

$A_{\mathbf{v};k}=\\{(\tau_{1},\dots,\tau_{k-1},\tau_{k+1},\dots,\tau_{n})\in\mathbb{R}^{n-1}:\tau_{j}>0\text{
for }j\in\mathbf{v}\setminus\\{k\\}\text{ and }\tau_{j}\leq 0\text{ for
}j\notin\mathbf{v}\\}.$

With these notations, it holds that

$\displaystyle J_{t}$
$\displaystyle\coloneqq\sum_{n=S+I}^{n_{\max}}\sum_{k=1}^{n-S}\int_{\mathbb{R}^{n-S}}\mathbf{1}_{\left\\{i(\tau)=I\right\\}}\partial_{\tau_{k}}\rho^{X,n,S}_{t}(\tau)d\tau=\sum_{n=S+I}^{n_{\max}}\sum_{k=1}^{n-S}\sum_{\mathbf{v}\in\mathbb{B}(n-S,I)}\int_{A_{\mathbf{v}}}\partial_{\tau_{k}}\rho^{X,n,S}_{t}(\tau)d\tau$
$\displaystyle=\sum_{n=S+I}^{n_{\max}}\sum_{k=1}^{n-S}\sum_{\mathbf{v}\in\mathbb{B}(n-S,I)}(-1)^{\mathbf{1}_{\left\\{k\in\mathbf{v}\right\\}}}\int_{A_{\mathbf{v};k}}\rho^{X,n,S}_{t}(\tau_{1},\dots,\tau_{k-1},0,\tau_{k+1},\dots,\tau_{n})\prod_{\begin{subarray}{c}1\leq
j\leq n\\\ j\neq k\end{subarray}}d\tau_{j},$

where the second line follows from the first by integrating
$\partial_{\tau_{k}}\rho^{X,n,S}_{t}$ over $A_{\mathbf{v}}$ with respect to
$\tau_{k}$. The goal is to show that

$J_{t}=\gamma\left((I+1)n^{X}_{(S,I+1)}(t)\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}}-In_{(S,I)}(t)\right),$
(29)

which will conclude the proof of Equation (26).

In order to do so, let us define the following function for $x\in E$:

$F^{(S,I)}(x)=\sum_{n=S+I}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n\\\
\mathbf{s}(x)=S\end{subarray}\right\\}}\sum_{k=1}^{n-S}\sum_{\mathbf{v}\in\mathbb{B}(n-S,I)}(-1)^{\mathbf{1}_{\left\\{k\in\mathbf{v}\right\\}}}\mathbf{1}_{\left\\{\tau_{k}(x)=0\right\\}}H^{(n,S),\mathbf{v},k}(\tau(x)),$

where for any $\mathbf{v}\subseteq\llbracket 1,n-S\rrbracket$, $k\in\llbracket
1,n-S\rrbracket$ and $\tau\in\mathbb{R}^{n_{\max}}$,

$H^{(n,S),\mathbf{v},k}(\tau)=\prod_{\begin{subarray}{c}j\in\mathbf{v}\\\
j\neq k\end{subarray}}h(\tau_{j})\prod_{\begin{subarray}{c}j\in\llbracket
1,n-S\rrbracket\setminus\mathbf{v}\\\ j\neq k\end{subarray}}(1-h(\tau_{j})).$

It then follows from Proposition 3.9 that
$J_{t}=\langle\eta^{X}_{t},F^{(S,I)}\rangle$.

Furthermore, distinguishing the cases $k\in\mathbf{v}$ and $k\notin\mathbf{v}$
yields that

$F^{(S,I)}=a^{(S,I)}-b^{(S,I)},$

where for any $x\in E$,

$a^{(S,I)}(x)=\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n\\\
\mathbf{s}(x)=S\end{subarray}\right\\}}\sum_{k=1}^{n-S}\mathbf{1}_{\left\\{\tau_{k}(x)=0\right\\}}\sum_{\begin{subarray}{c}\mathbf{v}\in\mathbb{B}(n-S,I)\\\
k\notin\mathbf{v}\end{subarray}}H^{(n,S),\mathbf{v},k}(\tau(x)),$

and

$b^{(S,I)}(x)=\sum_{n=S+I}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n\\\
\mathbf{s}(x)=S\end{subarray}\right\\}}\sum_{k=1}^{n-S}\mathbf{1}_{\left\\{\tau_{k}(x)=0\right\\}}\sum_{\begin{subarray}{c}\mathbf{v}\in\mathbb{B}(n-S,I)\\\
k\in\mathbf{v}\end{subarray}}H^{(n,S),\mathbf{v},k}(\tau(x)).$

The heart of the proof then lies in the following result.

######  Lemma 3.14.

Under the assumptions of Theorem 2.2, for any $X\in\\{H,W\\}$ and
$(S,I)\in\mathbb{S}$, for any $t\geq 0$,

$\langle\eta^{X}_{t},a^{(S,I)}\rangle=\gamma(I+1)n^{X}_{(S,I+1)}(t)\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}},$
(30)

and

$\langle\eta^{X}_{t},b^{(S,I)}\rangle=\gamma In^{X}_{(S,I)}(t).$ (31)

Before establishing Lemma 3.14, let us finish the proof of Proposition 3.13.
It is clear that this lemma implies Equation (29), and thus Equation (26).

In order to obtain the initial condition of Equation (27), notice that by
definition of $\eta_{0,\varepsilon}$,

$\displaystyle\langle\eta^{X}_{0},f^{(S,I)}\rangle$
$\displaystyle=\sum_{n=S+I}^{n_{\max}}\pi^{X}_{n}\binom{n}{S}(1-\varepsilon)^{S}\varepsilon^{n-S}\sum_{\mathbf{v}\in\mathbb{B}(n-S,I)}\int_{\mathbb{R}^{n-S}}\prod_{k\in\mathbf{v}}h(\tau_{k})\prod_{k\in\llbracket
1,n-S\rrbracket\setminus\mathbf{v}}(1-h(\tau_{k}))\nu(d\tau_{1})\dots\nu(d\tau_{n-S})$
$\displaystyle=\sum_{n=S+I}^{n_{\max}}\pi^{X}_{n}\binom{n}{S}(1-\varepsilon)^{S}\varepsilon^{n-S}\sum_{\mathbf{v}\in\mathbb{B}(n-S,I)}\langle\nu,h\rangle^{\\#\mathbf{v}}\langle\nu,1-h\rangle^{(n-S)-\\#\mathbf{v}}.$

Obviously, whenever $(n-S)-\\#\mathbf{v}>0$, the term vanishes as
$\langle\nu,1-h\rangle=0$. Hence only the case $\mathbf{v}=\llbracket
1,n-S\rrbracket$ remains, which in turn corresponds to $n=S+I$. This concludes
the proof. ∎

Let us now turn to Lemma 3.14.

###### Proof of Lemma 3.14.

Throughout the proof, we will continue using notations from the proof of
Proposition 3.13. We will detail the arguments leading to Equation (30), as
Equation (31) is obtained analogously.

Let $(S,I)\in\mathbb{S}$. We may start by noticing that in the case
$S+I=n_{\max}$, it holds that $a^{(S,I)}(x)=0$ for all $x\in E$. Indeed, the
sum over $n$ reduces to the case $n=n_{\max}$, and the set
$\mathbb{B}(n_{\max}-S,I)$ contains exactly one element, namely $\llbracket
1,I\rrbracket$. Hence for any $k\in\llbracket 1,n_{\max}-S\rrbracket$, the set
$\\{\mathbf{v}\in\mathbb{B}(n_{\max}-S,I):k\notin\mathbf{v}\\}$ is empty,
implying that $a^{(S,I)}$ is null on $E$. From now on, we will thus assume
that $S+I<n_{\max}$.

Define the function $g^{(S,I)}$ on $E$ as follows. For any $x\in E$,

$g^{(S,I)}(x)=\sum_{n=S+I}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n\\\
\mathbf{s}(x)=S\end{subarray}\right\\}}\sum_{k=1}^{n-S}h(\tau_{k}(x))\sum_{\begin{subarray}{c}\mathbf{v}\in\mathbb{B}(n-S,I-1)\\\
k\notin\mathbf{v}\end{subarray}}H^{(n,S),\mathbf{v},k}(\tau(x)).$

In particular, it holds by definition of $g^{(S,I)}$ that

$\langle\eta^{X}_{t},g^{(S,I)}\rangle=In^{X}_{(S,I)}(t).$

Since we assume $S+I<n_{\max}$, the state $(S,I+1)$ belongs to $\mathbb{S}$ as
well. Thus, we actually aim to show that for $G^{(S,I+1)}=\gamma
g^{(S,I+1)}-a^{(S,I)}$,

$\forall t\geq 0,\;\langle\eta^{X}_{t},G^{(S,I+1)}\rangle=0.$

In order to achieve this, we will consider the following smooth approximations
of $a^{(S,I)}$ and $g^{(S,I+1)}$ respectively, defined for $\alpha>0$. For any
$x\in E_{X}$,

$\displaystyle a^{(S,I)}_{\alpha}(x)$
$\displaystyle=\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n,\\\
\mathbf{s}(x)=S\end{subarray}\right\\}}\sum_{k=1}^{n-S}\psi_{\alpha}(\tau_{k}(x))H^{n-S,I,k}_{\alpha}(\tau(x)),$
$\displaystyle g^{(S,I+1)}_{\alpha}(x)$
$\displaystyle=\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n,\\\
\mathbf{s}(x)=S\end{subarray}\right\\}}\sum_{k=1}^{n-S}h_{\alpha}(\tau_{k}(x))H^{n-S,I,k}_{\alpha}(\tau(x)),$

where for any integers $J\leq N$ and $k\leq N$, for any
$\tau\in\mathbb{R}^{n_{\max}}$,

$H^{N,J,k}_{\alpha}(\tau)=\sum_{\begin{subarray}{c}\mathbf{v}\in\mathbb{B}(N,J)\\\
k\notin\mathbf{v}\end{subarray}}\prod_{j\in\mathbf{v}}h_{\alpha}(\tau_{j})\prod_{\begin{subarray}{c}j\in\llbracket
1,N\rrbracket\setminus\mathbf{v}\\\ j\neq
k\end{subarray}}(1-h_{\alpha}(\tau_{j})).$

Letting $G^{(S,I+1)}_{\alpha}=\gamma g^{(S,I+1)}_{\alpha}-a^{(S,I)}_{\alpha}$,
dominated convergence ensures that

$\langle\eta^{X}_{t},G^{(S,I+1)}_{\alpha}\rangle\xrightarrow[\alpha\to
0]{}\langle\eta^{X}_{t},G^{(S,I+1)}\rangle.$

Classical properties of convolution imply that both $a^{(S,I)}_{\alpha}$ and
$g^{(S,I+1)}_{\alpha}$ belong to $\mathcal{C}^{1}_{b}(E,\mathbb{R})$, allowing
us to apply Proposition 3.11 to $G^{(S,I+1)}_{\alpha}(x)$. In order to do so,
some preliminary computations are necessary. Indeed, we need to compute
$f_{T,t}$ and $f^{\mathcal{I}}_{T,t}$ as defined in Equation (18) for
$f=a^{(S,I)}_{\alpha}$ and $f=g^{(S,I+1)}_{\alpha}$. Let $T>0$, then for any
$t\in[0,T]$, writing $e_{1,n-S}=(1_{n-S},0_{n_{\max}-(n-S)})$ for the vector
of dimension $n_{\max}$ whose $(n-S)$ first components are set to 1 and the
following all equal 0,

$\left(a^{(S,I)}_{\alpha}\right)_{T,t}(x)=\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n,\\\
\mathbf{s}(x)=S\end{subarray}\right\\}}\sum_{k=1}^{n-S}\psi_{\alpha}(\tau_{k}(x)-(T-t))H^{n-S,I,k}_{\alpha}\left(\tau(x)-(T-t)e_{1,n-S}\right).$

Similarly, for any $\sigma>0$, distinguishing the cases where either $k=n-S$,
or $k\neq n-S$ and $n-S\in\mathbf{v}$, or $k\neq n-S$ and
$n-S\notin\mathbf{v}$ leads to

$\displaystyle\left(a^{(S,I)}_{\alpha}\right)_{T,t}$
$\displaystyle(\mathfrak{j}(x,\sigma))=\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n,\\\
\mathbf{s}(x)-1=S\end{subarray}\right\\}}\Big{[}\psi_{\alpha}(\sigma-(T-t))H^{n-S,I,n-S}_{\alpha}(\tau(x)-(T-t)e_{1,n-S-1})$
$\displaystyle+h_{\alpha}(\sigma-(T-t))\sum_{k=1}^{n-S-1}\psi_{\alpha}(\tau_{k}-(T-t))H^{n-S-1,I-1,k}_{\alpha}\left(\tau(x)-(T-t)e_{1,n-S-1}\right)$
$\displaystyle+(1-h_{\alpha}(\sigma-(T-t)))\sum_{k=1}^{n-S-1}\psi_{\alpha}(\tau_{k}-(T-t))H^{n-S-1,I,k}_{\alpha}\left(\tau(x)-(T-t)e_{1,n-S-1}\right)\Big{]}.$

Define $m^{h}_{T,t}(\alpha)=\langle\nu,h_{\alpha}(\cdot-(T-t))\rangle$ and
$m^{\psi}_{T,t}(\alpha)=\langle\nu,\psi_{\alpha}(\cdot-(T-t))\rangle$. Notice
that the second line is equal to zero if $I=0$, as the set
$\mathbb{B}(n-S-I,-1)$ is empty. Similarly, in the case of the third line,
when $n=S+1+I$ the set $\\{\mathbf{v}\in\mathbb{B}(n-S-1,I)\\}$ reduces to
$\llbracket 1,I\rrbracket$ and hence for any $k\in\llbracket 1,n-S\rrbracket$,
the set $\\{\mathbf{v}\in\mathbb{B}(n-S-1,I):k\neq\mathbf{v}\\}$ is empty. In
particular, the third line vanishes if $S+I+1=n_{\max}$. We may thus recognize
that

$\displaystyle\left(a^{(S,I)}_{\alpha}\right)_{T,t}^{\mathcal{I}}(x)$
$\displaystyle=m^{\psi}_{T,t}(\alpha)\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n,\\\
\mathbf{s}(x)-1=S\end{subarray}\right\\}}H^{n-S,I,n-S}_{\alpha}(\tau(x)-(T-t)e_{1,n-S-1})$
$\displaystyle+m^{h}_{T,t}(\alpha)a^{(S+1,I-1)}_{\alpha;T,t}\mathbf{1}_{\left\\{I\geq
1\right\\}}+(1-m^{h}_{T,t}(\alpha))a^{(S+1,I)}_{\alpha;T,t}\mathbf{1}_{\left\\{S+I+1<n_{\max}\right\\}}.$

Proceeding in the same way also yields

$\displaystyle\left(g^{(S,I+1)}_{\alpha}\right)^{\mathcal{I}}_{T,t}(x)$
$\displaystyle=m^{h}_{T,t}(\alpha)\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n,\\\
\mathbf{s}(x)-1=S\end{subarray}\right\\}}H^{n-S,I,n-S}_{\alpha}(\tau(x)-(T-t)e_{1,n-S-1})$
$\displaystyle+m^{h}_{T,t}(\alpha)g^{(S+1,I)}_{\alpha;T,t}\mathbf{1}_{\left\\{I\geq
1\right\\}}+(1-m^{h}_{T,t}(\alpha))g^{(S+1,I+1)}_{\alpha;T,t}\mathbf{1}_{\left\\{S+I+1<n_{\max}\right\\}}.$

In particular, we thus obtain that

$\displaystyle($ $\displaystyle
G^{(S,I+1)}_{\alpha})^{\mathcal{I}}_{T,t}(x)=m^{h}_{T,t}(\alpha)G^{(S+1,I)}_{\alpha;T,t}\mathbf{1}_{\left\\{I\geq
1\right\\}}+(1-m^{h}_{T,t}(\alpha))G^{(S+1,I+1)}_{\alpha;T,t}\mathbf{1}_{\left\\{S+I+1<n_{\max}\right\\}}$
$\displaystyle+(\gamma
m^{h}_{T,t}(\alpha)-m^{\psi}_{T,t}(\alpha))\sum_{n=S+I+1}^{n_{\max}}\mathbf{1}_{\left\\{\begin{subarray}{c}\mathbf{n}(x)=n,\\\
\mathbf{s}(x)-1=S\end{subarray}\right\\}}H^{n-S,I,n-S}_{\alpha}(\tau(x)-(T-t)e_{1,n-S-1}).$

Notice further that, on the one hand,

$m^{\psi}_{T,t}(\alpha)=\int_{0}^{+\infty}\gamma\mathrm{e}^{-\gamma
x}\psi_{\alpha}(x-(T-t))dx=\int_{0}^{+\infty}\gamma\mathrm{e}^{-\gamma(x+(T-t))}\psi_{\alpha}(x)dx\xrightarrow[\alpha\to
0]{}\gamma\mathrm{e}^{-\gamma(T-t)}.$

On the other hand,

$\displaystyle m^{h}_{T,t}(\alpha)$
$\displaystyle=\int_{0}^{+\infty}\gamma\mathrm{e}^{-\gamma
x}\int_{0}^{+\infty}\psi_{\alpha}(x-(T-t)-z)dzdx=\int_{\mathbb{R}}\int_{y}^{+\infty}\gamma\mathrm{e}^{-\gamma
x}dx\psi_{\alpha}(y-(T-t))dy$
$\displaystyle=\int_{\mathbb{R}}\mathrm{e}^{-\gamma(y+(T-t))}\psi_{\alpha}(y)dy\xrightarrow[\alpha\to
0]{}\mathrm{e}^{-\gamma(T-t)}.$

Thus, for any $T>0$, we may apply Proposition 3.11 to $G^{(S,I+1)}_{\alpha}$,
and let $\alpha$ go to zero in order to obtain by dominated convergence that

$\displaystyle\langle\eta^{X}_{T},G^{(S,I+1)}\rangle=\langle\eta^{X}_{0},G^{(S,I+1)}_{T,0}\rangle$
(32)
$\displaystyle+\lambda_{X}\int_{0}^{T}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}\left(\mathrm{e}^{-\gamma(T-t)}G^{(S+1,I)}_{T,t}\mathbf{1}_{\left\\{I\geq
1\right\\}}+(1-\mathrm{e}^{-\gamma(T-t)})G^{(S+1,I+1)}_{T,t}\mathbf{1}_{\left\\{S+I+1<n_{\max}\right\\}}-G^{(S,I+1)}_{T,t}\right)\rangle
dt$
$\displaystyle+\int_{0}^{T}\Lambda_{t}\langle\eta^{X}_{t},\mathbf{s}\left(\mathrm{e}^{-\gamma(T-t)}G^{(S+1,I)}_{T,t}\mathbf{1}_{\left\\{I\geq
1\right\\}}+(1-\mathrm{e}^{-\gamma(T-t)})G^{(S+1,I+1)}_{T,t}\mathbf{1}_{\left\\{S+I+1<n_{\max}\right\\}}-G^{(S,I+1)}_{T,t}\right)\rangle
dt.$

It is possible to compute $\langle\eta^{X}_{0},G^{(S,I+1)}_{T,0}\rangle$.
Indeed, by definition of $\eta_{0}=\eta_{0,\varepsilon}$,

$\displaystyle\langle\eta^{X}_{0},a^{(S,I)}_{T,0}\rangle=\sum_{n=S+I+1}^{n_{\max}}$
$\displaystyle\pi^{X}_{n}\binom{n}{S}(1-\varepsilon)^{S}\varepsilon^{n-S}$
$\displaystyle\times\Big{[}\sum_{k=1}^{n-S}\langle\nu,\delta_{0}(\cdot-T)\rangle\sum_{\begin{subarray}{c}\mathbf{v}\in\mathbb{B}(n-S,I)\\\
k\notin\mathbf{v}\end{subarray}}\prod_{j\in\mathbf{v}}\langle\nu,h(\cdot-T)\rangle\prod_{\begin{subarray}{c}j\in\llbracket
1,n-S\rrbracket\setminus\mathbf{v}\\\ j\neq
k\end{subarray}}\langle\nu,1-h(\cdot-T)\rangle\Big{]}.$

Of course, the computation of $\langle\eta^{X}_{0},g^{(S,I+1)}_{T,0}\rangle$
is very similar. It suffices to replace the term
$\langle\nu,\delta_{0}(\cdot-T)\rangle$ by $\langle\nu,h(\cdot-T)\rangle$. As
$\langle\nu,\delta_{0}(\cdot-T)\rangle=\gamma\langle\nu,h(\cdot-T)\rangle$, it
follows that

$\langle\eta^{X}_{0},G^{(S,I+1)}_{T,0}\rangle=0.$ (33)

Thus the term $\langle\eta^{X}_{0},G^{(S,I+1)}_{T,0}\rangle$ of Equation (32)
vanishes. Letting
$C=n_{\max}(\lambda_{X}n_{\max}+\lambda_{\overline{X}}n_{\max}+\beta_{G})$, we
obtain that

$\small\lvert\langle\eta^{X}_{T},G^{(S,I+1)}\rangle\rvert\leq
C\int_{0}^{T}\left(\lvert\langle\eta^{X}_{t},G^{(S+1,I)}_{T,t}\rangle\rvert\mathbf{1}_{\left\\{I\geq
1\right\\}}+\lvert\langle\eta^{X}_{t},G^{(S+1,I+1)}_{T,t}\rangle\rvert\mathbf{1}_{\left\\{S+I+1<n_{\max}\right\\}}+\lvert\langle\eta^{X}_{t},G^{(S,I+1)}_{T,t}\rangle\rvert\right)dt.$
(34)

Notice that if $(S,I)\in\mathbb{S}$ is such that $S+I<n_{\max}$, then
$(S^{\prime},I^{\prime})=(S,I+1)$ is such that
$(S^{\prime},I^{\prime})\in\mathbb{S}$, $S^{\prime}+I^{\prime}\geq 3$ and
$I^{\prime}\geq 1$. Similarly, on the event $I\geq 1$,
$(S^{\prime},I^{\prime})=(S+1,I)$ satisfies the same conditions, as well as
$(S^{\prime},I^{\prime})=(S+1,I+1)$ if $S+I+1<n_{\max}$. In order to conclude,
it thus suffices to establish the following lemma.

######  Lemma 3.15.

For any $(S,I)\in\mathbb{S}$ such that $I\geq 1$ and $S+I\geq 3$, for any
$T\geq 0$,

$\forall t\in[0,T],\;\langle\eta^{X}_{t},G^{(S,I)}_{T,t}\rangle=0.$ (35)

The desired result then follows from Equation (34). ∎

###### Proof of Lemma 3.15.

The key argument will be to apply Gronwall’s lemma to
$t\mapsto\langle\eta^{X}_{t},G^{(S,I)}_{T,t}\rangle$. We thus need to show
that for any $(S,I)\in\mathbb{S}$ such that $S+I\geq 3$ and $I\geq 1$, the
function $t\mapsto\langle\eta^{X}_{t},G^{(S,I)}_{T,t}\rangle$ is continuous on
$[0,T]$. In order to do so, notice that we may apply Proposition 3.11 to
$t\mapsto(G^{(S,I)}_{\alpha})_{T,t}$ which belongs to
$\mathcal{C}^{1}_{b}(\mathbb{R}_{+}\times E)$. As for any $z\in\mathbb{R}$ and
$u\in[0,t]$ it holds that $z-(T-t)-(t-u)=z-(T-u)$, it follows that

$\left((G^{(S,I)}_{\alpha})_{T,t}\right)_{t,u}=(G^{(S,I)}_{\alpha})_{T,u}.$

Thus, proceeding like before and letting $\alpha$ go to zero, dominated
convergence allows to conclude that for any $t\in[0,T]$,

$\displaystyle\langle$
$\displaystyle\eta^{X}_{t},G^{(S,I)}_{T,t}\rangle=\langle\eta_{0}^{X},G^{(S,I)}_{T,0}\rangle$
$\displaystyle+\int_{0}^{t}\lambda_{X}\langle\eta^{X}_{u},\mathbf{s}\mathbf{i}\left(\mathrm{e}^{\gamma(T-u)}G^{(S+1,I-1)}_{T,u}\mathbf{1}_{\left\\{I-1\geq
1\right\\}}+(1-\mathrm{e}^{\gamma(T-u)})G^{(S+1,I)}_{T,u}\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}}-G^{(S,I)}_{T,u}\right)\rangle
du$
$\displaystyle+\int_{0}^{t}\Lambda_{u}\langle\eta^{X}_{u},\mathbf{s}\left(\mathrm{e}^{\gamma(T-u)}G^{(S+1,I-1)}_{T,u}\mathbf{1}_{\left\\{I-1\geq
1\right\\}}+(1-\mathrm{e}^{\gamma(T-u)})G^{(S+1,I)}_{T,u}\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}}-G^{(S,I)}_{T,u}\right)\rangle
du.$

The initial condition $t=0$ vanishes according to Equation (33) as it holds
for any $(S,I)\in\mathbb{S}$ satisfying $S+I<n_{\max}$, and indeed