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

Madeleine Kubasch

<figure><img src="image_1.png" /></figure>

## To cite this version:

Madeleine Kubasch. Large population limit for a multilayer SIR model including households and workplaces. 2023. hal-04106800v3

HAL Id: hal-04106800

https://hal.science/hal-04106800v3

Preprint submitted on 13 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
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

October 13, 2023

**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 specifies the remaining infectious period length for each infected. 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.

**Keywords.** Measure-valued process, SDE with jumps, large population limit, model re-
duction, epidemic process, household-workplace models, two layers of mixing.

Code availability. https://github.com/m-kubasch/household-workplace-model

# 1 Introduction

Epidemic spread depends by essence on the way individuals interact with one another. As a consequence, models have been developed which take into account main features of real-life contacts, for instance through contact networks [23]. In particular, some attention has been drawn to studying clustered networks, as clustering has a strong impact on epidemic spread [15, 42] which is intimately related to the way clustering is achieved within the network [17]. A particular form of clustering consists in the presence of entirely connected small social structures, such as households or workplaces, which exist in addition to random contacts in the general population. This kind of population structure is captured by models with two levels of mixing [8, and references therein], which are related to efficient control measures. Indeed, both COVID-19 and influenza epidemics illustrate the pertinence of teleworking and school closures [29, 37, 26], and models explicitly distinguishing different contact types are well suited to simulate the impact of these measures [11]. However, precisely understanding the impact on disease propagation of the way individuals are organized in households and workplaces is not straightforward [5].
This motivates the study of models with several levels of contact, which as we shall see
lead to interesting mathematical issues due to their multiscale population structure. Pellis et
al. [33] have proposed a model with two levels of mixing structured in three layers of contacts:
households, workplaces and the general population. This household-workplace model has already
been studied to some extent, establishing for instance the epidemic growth rate [34, 5] and several
reproduction numbers [3]. In particular, the R₀ used throughout this paper was introduced in
[33], and computations of R₀ and proportions of infections per contact layer will make use of
the working package associated to [5].

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 particular, large population approximations of stochastic models have
proven fruitful to achieve such model reductions in many contexts, among which epidemics on
random graphs.

Historically, the standard SIR model developed by Kermack and McKendrick itself corre-
sponds to the large population limit of the uniformly mixing stochastic SIR model. In the
Markovian setting, the convergence of the stochastic model to its deterministic limit can be
established using classical results on the convergence of finite type density-dependent Markov
jump processes, e.g. [2]. When infectious periods are not restricted to being exponentially dis-
tributed, the large population convergence of the stochastic model to the unique deterministic
solution of a system of integral equations can also be obtained [24, 13]. For more complex con-
tact networks however, it often is challenging to propose closed systems of equations correctly
describing the epidemic dynamics.

A well-understood case is the *SIR* model on the configuration graph, for which a reduced
model referred to as edge-based compartmental model (EBCM) [41, 31] has been proven to be
the large population limit of the underlying stochastic model [9, 19]. Since then, the equivalence
with other reduced models has been established under appropriate assumptions [17, 43, 18, 22],
and the EBCM formalism has been extended to related models [36, 18]. The configuration graph
is a favourable setting for this analysis thanks to the absence of clustering in the large graph
limit, which however also constitutes a major limitation, prohibiting for instance the existence
of household-like structures.

Some attention has thus been drawn to models where each individual belongs to a random
number of fully connected subgraphs (*cliques*) of the same type, hence being closely related
to the household-workplace model. Several reduced models have been proposed, including an
EBCM [42, 38, 15]. To our knowledge, the convergence of the underlying stochastic model to
the proposed reduced model has not been established in any of these settings. Notice that these
models share a major common point, which will also hold true in our setting: they focus on the
epidemic at the level of structures, as they keep track of the proportions of cliques containing a
certain number of susceptibles and infected, leading to high-dimensional dynamical systems for
larger clique sizes.

When considering two levels of mixing, the first model for which a large population limit
has been determined is the household model, which assumes a uniformly mixing general popu-
lation and that each individual belongs to exactly one household [16], thus being a special case
of the household-workplace model. Here, the stochastic model can again be formalized as a
finite type density-dependent Markov jump process, ensuring the large population convergence
to the deterministic model. If either of these two assumptions is relaxed, e.g. considering a
configuration graph at the global level [11, 27] or individuals belonging to several households [6],
reduced models have been proposed, but without rigorous derivation from stochastic models.
In particular, the case of the household-workplace model is not covered, and the only reduced models proposed so far approach the epidemic dynamics using well calibrated uniformly mixing models [5, 10]. While these are capable of capturing some key characteristics of the epidemic, such as the epidemic peak size and final size, they do not allow for an accurate prediction of the epidemic dynamic over time.

As a consequence, in this paper, we will study the large population limit of the multilayer SIR model with 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 sequence of processes converges in law when the size of the population grows to infinity. This allows to identify a new model reduction, and establishes that it is asymptotically exact. Besides, it paves the way for more quantitative estimates on this approximation.

Notice here that 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 mathematical representation of the model. This difficulty actually arises as soon as one considers the probability of an individual belonging to several cliques at once, whether they are of different types or not, and we refer to [4] where a similar approach has been developed for branching approximations. Let us emphasize that this model formulation allows to immediately consider a wide range of infectious period length distributions instead of being restricted to the Markovian case, which is a pertinent generalisation for many epidemic models [36, 13, 12, 25].

The model will be represented by 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. 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 [16, 38].

The present paper is structured as follows. Section 2 introduces the individual-based model, and Section 3.1 subsequently presents the convergence results in detail. Section 3.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 4 contains the proofs of our results.

Before proceeding, let us introduce some notations that will be used throughout the paper. For any integers $n \le m$, we write $\llbracket n, m \rrbracket = \{n, \dots, m\}$. For a measurable space $(E, \mathcal{E})$, let $M_P(E)$ be the set of point measures, $M_F(E)$ the set of finite measures and $M_1(E)$ the set of probability measures on $E$. We define $M_{P,1}(E) = M_P(E) \cap M_1(E)$ the set of punctual probability measures on $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 f d\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 $C(E, \mathbb{R}^m)$ be the set of continuous functions $f : E \to \mathbb{R}^m$. Similarly, $C_b(E, \mathbb{R}^m)$ is defined as the subset of bounded functions $f \in C(E, \mathbb{R}^m)$. Finally, the space $C_b^1(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.

## 2 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 using a measure-valued stochastic differential equation.

## 2.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 [5]. Suppose that households and workplaces are of size at least one and at most $n_{max}$. Consider distributions $(\pi_j^H)$ and $(\pi_j^W)$ on $[1, n_{max}]$. These distributions correspond to the large population limit of household and workplace size distributions, in the sense that in an infinite population, a proportion $\pi_j^H$ 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_j^X.$$

On a probability space $(\Omega_G, \mathbb{P}_G, \mathcal{F}_G)$, we construct a sequence $(\mathbf{G}^K)_{K \ge 1}$ of this random population structure as follows. For $K \ge 1$, let $k \in [0, K]$ 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$.

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. 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, 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$. Indeed, each given encounter within the general population becomes less likely when the population size $K$ increases.

Finally, infected individuals remain infectious for a period of time which is independent from $(\mathbf{G}_K)_{K \ge 1}$ and distributed according to a probability distribution $\nu$ on $\mathbb{R}_+$, which we assume to be absolutely continuous with regard to the Lebesgue measure. 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.
## 2.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 [3], as well as the epidemic growth rate of the household-workplace model [34]. 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 [4], we will thus further characterize each structure by the infectious periods of its infected members. Adopting this point of view is key, as it allows to handle both the progressive discovery of the graph as the epidemic spreads, and the correlations arising from infected individuals, without explicitly keeping in memory the discovered graph.

Let $\omega \in \Omega_G$. For a population of size $K \ge 1$, let $K_H$ be the number of households and $K_W$ the number of workplaces in $G^K(\omega)$. While $K_H$ and $K_W$ depend on $\omega$, this dependency is not specified explicitly for readability. This will also apply to the forthcoming notations. Label the $K_H$ households in an arbitrary fashion $1, \dots, K_H$. Consider the set

$$E = \{(n, s, \tau) \in [1, n_{\max}] \times [0, n_{\max}] \times \mathbb{R}^{n_{\max}} : s \le n; \forall j > n - s, \tau_j = 0\}.$$

Then for $k \in [1, K_H]$, the k-th household is characterized at time $t \ge 0$ by its type

$$x_k^H(t) = (n_k^H, s_k^H(t), \tau_k^H(t)) \in E.$$

The first two components of $x_k^H$ 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_k^H$ is a vector containing the remaining infectious periods of the members of the household. Indeed, at time $t$, there are $n_k^H - s_k^H(t)$ infected or removed individuals within the household. For $j \in [1, n_k^H - s_k^H(t)]$, if $\tau_{k,j}^H(t) > 0$, the individual is still infectious and will remain so for $\tau_{k,j}^H(t)$ units of time. Otherwise, if $\tau_{k,j}^H(t) \le 0$, the individual has recovered, and the recovery has occurred at time $t - |\tau_{k,j}^H(t)|$. For $j > n_k^H - s_k^H(t)$, $\tau_{k,j}^H(t)$ has no interpretation, and is set to zero for convenience in computations. In other words, the infectiousness of a previously contaminated individual with remaining infectious period $\tau$ is given by $\mathbf{1}_{\{\tau>0\}}$.

Similarly, label the $K_W$ workplaces in an arbitrary order $1, \dots, K_W$. For $\ell \in [1, K_W]$, the $\ell$-th workplace is characterized by its type $x_\ell^W(t) = (n_\ell^W, s_\ell^W(t), \tau_\ell^W(t))$, which is defined analogously to household types.

Notice that all of these quantities depend on the population size $K$, but this dependency is omitted to simplify notations.

By definition, these types evolve over time. On the one hand, for any $X \in \{H, W\}$, for any $k \in [1, K_X]$ and $j \in [1, n_{\max}]$, the j-th component of $\tau_k^X$ 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 [1, n_{\max}], \quad \frac{d}{dt} \tau_{k,j}^X(t) = -\mathbf{1}_{\{j \le n_k^X - s_k^X(t)\}}.$$
Let $(e_j)_{1 \le j \le n_{\max}}$ denote the canonical basis of $\mathbb{R}^{n_{\max}}$. Then for any $0 \le t \le 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_k^H$ and $x_\ell^W$ jump from $x_k^H(t-)$ and $x_\ell^W(t-)$ to $j(x_k^H(t-), \sigma)$ and $j(x_\ell^W(t-), \sigma)$ respectively, where for any $x = (n, s, \tau) \in E$,

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

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_k^H(t-)/S(t-)$, and independently to the $\ell$-th workplace with probability $s_\ell^W(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 introduce the stochastic process $(\zeta_t^K = (\zeta_t^{H|K}, \zeta_t^{W|K}))_{t \ge 0}$ taking values in $\mathfrak{M}_{P,1} = \mathcal{M}_{P,1}(E) \times \mathcal{M}_{P,1}(E)$. $\zeta_t^{H|K}$ and $\zeta_t^{W|K}$ 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 \ge 0$ and $X \in \{H, W\}$,

$$ \zeta_t^{X|K} = \frac{1}{K_X} \sum_{k=1}^{K_X} \delta_{x_k^X(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}} \le K_X \le K. $$

Thus, studying the asymptotic $K \to \infty$ amounts to $(K_H, K_W) \to (\infty, \infty)$.

Observe that the number of infected individuals in a household of type $(n, s, \tau)$ is given by $i(\tau) = \sum_{k=1}^{n_{\max}} \mathbf{1}_{\{\tau_k > 0\}}$. Then for any $t \ge 0$,

$$ I_H(t) = \frac{1}{K_H} \sum_{k=1}^{K_H} i(\tau_k^H(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$, as well as the workplace-related quantities $I_W(t)$ and $S_W(t)$. Then

$$ \forall X \in \{H, W\}, S(t) = K_X S_X(t) \text{ and } I(t) = K_X I_X(t). \quad (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. Notice that for $X \in \{H, W\}$, $S_X$ and $I_X$ actually depend on the population size $K$, which is omitted in notations for readability.

Finally, let us briefly emphasize that the partition of the population in households and workplaces is entirely conveyed by $\zeta_0^K \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 $G^K(\omega)$. Similarly, there are some natural constraints on $\zeta_0^K$, as $\zeta^{|K|}$ and $\zeta^{|W|}$ describe the same population, once dispatched into households, and once dispatched into workplaces. 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 or recovered needs to have the same remaining infectious period in both his household and workplace. In other words, almost surely,

$$ \{\tau_{k,j}^H(0): 1 \le k \le K_H, 1 \le j \le n_k^H - s_k^H(0)\} = \{\tau_{\ell,j}^W(0): 1 \le \ell \le K_W, 1 \le j \le n_\ell^W - s_\ell^W(0)\}. \quad (2) $$

Forthcoming Lemma 2.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. Again, $K$-dependency is not specified explicitly in order to simplify notations. Let

$$ U_G = (\mathbb{R}_+)^3 \times [1, K_H] \times [1, K_W] \times \mathbb{R}_+, \quad (3) $$

and consider the following measure on $U_G$:

$$ \mu_G(du) = \mu_G(d\theta, dk, d\ell, d\sigma) = d\theta \otimes \mu_\#(dk) \otimes \mu_\#(d\ell) \otimes \nu(d\sigma), $$

where $d\theta$ and $\mu_\#$ denote the Lebesgue measure on $\mathbb{R}^3$ and the standard counting measure, respectively. For $t \ge 0$ and $u = (\theta, k, \ell, \sigma) \in U_G$ where $\theta = (\theta_1, \theta_2, \theta_3)$, let us define

$$ \mathcal{I}_G(t, u) = \mathbf{1}\left\{\theta_1 \le \frac{\beta_G}{K} S(t) I(t), \theta_2 \le \frac{s_k^H(t)}{S(t)}, \theta_3 \le \frac{s_\ell^W(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 = (\mathbb{R}_+)^2 \times [1, K_H] \times [1, K_W] \times \mathbb{R}_+, \quad (4) $$

endowed with the measure

$$ \mu(du) = \mu(d\theta, dk, d\ell, d\sigma) = d\theta \otimes \mu_\#(dk) \otimes \mu_\#(d\ell) \otimes \nu(d\sigma), $$

where with slight abuse of notation, $d\theta$ designates the Lebesgue measure on $\mathbb{R}^2$. Then for $t \ge 0$ and $u = (\theta, k, \ell, \sigma) \in U$ where $\theta = (\theta_1, \theta_2)$, we further introduce

$$ \mathcal{I}_H(t, u) = \mathbf{1}\left\{\theta_1 \le \lambda_H s_k^H(t) i(\tau_k^H(t)), \theta_2 \le \frac{s_\ell^W(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 \ge t \ge 0$ and $u = (\theta, k, \ell, \sigma) \in U_G \cup U$, 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(\Psi(j(x_k^H(t-), \sigma), T, t)) - \delta(\Psi(x_k^H(t-), T, t)). $$

Finally, define $\mathcal{I}_W(t, u)$ and $\Delta_W(u, T, t)$ analogously:

$$ \mathcal{I}_W(t, u) = \mathbf{1}\left\{\theta_1 \le \lambda_W s_{\ell}^W(t) i(\tau_{\ell}^W(t)), \theta_2 \le \frac{s_{\ell}^H(t)}{S(t)}\right\}, \\ \text{and} \\ \Delta_W(u, T, t) = \delta(\Psi(j(x_{\ell}^W(t-), \sigma), T, t)) - \delta(\Psi(x_{\ell}^W(t-), T, t)). $$

We are now ready to yield the main characterization of $\zeta_t^K$, as inspired by [39].

**Proposition 2.1.** Define on the same probability space as $\zeta_0^K$, and independently from $\zeta_0^K$, three independent Poisson point measures $Q_Y$ on $\mathbb{R}_+ \times U_Y$ with intensity $d\mu_Y(du)$, 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\}$ and $T \ge 0$,

$$ \zeta_T^{X|K} = \frac{1}{K_X} \left( \sum_{j=1}^{K_X} \delta_{\Psi(x_j^X(0), T, 0)} + \sum_{Y \in \{H, W, G\}} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) \Delta_X(u, T, t) Q_Y(dt, du) \right), \quad (5) $$

where $U_G$ and $U_H = U_W = U$ are defined by Equations (3) and (4), respectively.

The idea behind Equation (5) goes as follows. Let us focus for example on the distribution $\zeta_T^{H|K}$ of household types at time $T$. Each household's type contributes to the distribution at uniform weight $1/K_H$. 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_k^H(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_k^H(t) = j(x_k^H(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(j(x_k^H(t-), \sigma), T, t)$ instead of $\Psi(x_k^H(0), T, 0) = \Psi(x_k^H(t-), T, t)$. This reasoning is reflected in $\Delta_H$, and can be iterated over the whole of $[0, T]$. 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 \ge 0}$ of successive jump times of $\zeta^K$, where we define $T_0 = 0$. Then using a method similar to rejection sampling, $(T_n)_{n \ge 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}_+, \mathfrak{M}_{P,1}(E))$.

Finally, uniqueness is obtained by an induction argument which proves that for any $n \ge 0$, $(T_n, \zeta_{T_n}^K)$ is uniquely determined by $(\zeta_0^K, (Q_Y)_{Y \in S})$ where $S = \{G, H, W\}$. 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_{T_n}^K, (Q_Y)_{Y \in S})$ and $\zeta_{T_{n+1}}^K$ by $(T_{n+1}, T_n, \zeta_{T_n}^K, (Q_Y)_{Y \in S})$. The induction hypothesis allows to conclude. $\square$

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

**Lemma 2.2.** Suppose that almost surely, $K_H S_H(0) = K_W S_W(0)$ and Equation (2) holds. Then for any $t \ge 0$, $K_H S_H(t) = K_W S_W(t)$ and $K_H I_H(t) = K_W I_W(t)$, almost surely.
Proof. Let $T \ge 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_T^{X|K}, \mathbf{s} \rangle$ for $X \in \{H, W\}$. It follows from Equation (5) that

$$K_X S_X(T) = \sum_{j=1}^{K_X} \mathbf{s}(\Psi(x_j^X(0), T, 0)) + \sum_{Y \in \{H, W, G\}} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) \langle \Delta_X(u, T, t), \mathbf{s} \rangle Q_Y(dt, du).$$

Notice that on the one hand, for any $x \in E$ and $0 \le t \le 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 = (\theta, k, \ell, \sigma)$,

$$\langle \Delta_H(u, T, t), \mathbf{s} \rangle = \mathbf{s}(\Psi(j(x_k^H(t-), \sigma), T, t)) - \mathbf{s}(\Psi(x_k^H(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_T^{X|K}, \mathbf{i} \rangle$. This time, for any $u = (\theta, k, \ell, \sigma)$ and $0 \le t \le T$, it holds that $\langle \Delta_X(u, T, t), \mathbf{i} \rangle = \mathbf{1}_{\{\sigma>(T-t)\}}$. Finally, Equation (2) ensures that

$$\sum_{k=1}^{K_H} \mathbf{i}(\Psi(x_k^H(0), T, 0)) = \sum_{k=1}^{K_H} \sum_{j=1}^{n_k^H - s_k^K(0)} \mathbf{1}_{\{\tau_{k,j}^H(0)>T\}} = \sum_{\ell=1}^{K_W} \mathbf{i}(\Psi(x_\ell^W(0), T, 0)).$$

The conclusion follows as previously from Equation (5). □

### 3 Main results

In this section, we are going to present our main results on the convergence of $(\zeta^K)_{K \ge 1}$ in the Skorokhod space $D(\mathbb{R}_+, M_1(E))^2$. For $X \in \{H, W\}$, let $\bar{X}$ be the complementary structure type, i.e. $\bar{X} = W$ if $X = H$ and vice-versa.

As we are interested in studying the limit $K \to \infty$, an important ingredient will be the asymptotic behavior of the sequence of random graphs on which the epidemic spreads. Let $\pi^{H|K}$ and $\pi^{W|K}$ be the household and workplace size distributions observed in $\mathbf{G}^K$. The law of large numbers ensures that $(\pi^{H|K}, \pi^{W|K})_{K \ge 1}$ converges $\mathbb{P}_{\mathbf{G}}$-almost everywhere to $(\pi^H, \pi^W)$. We hence define

$$\Omega_G^* = \{ \omega \in \Omega_G : (\pi^{H|K}(\omega), \pi^{W|K}(\omega)) \xrightarrow[K \to \infty]{} (\pi^H, \pi^W) \},$$

and our main results will hold for $\omega \in \Omega_G^*$.

Further, the following assumption on the sequence of initial conditions $(\zeta_0^K)_{K \ge 1}$ will be required from now on.

**Assumption 3.1.** For any $X \in \{H, W\}$ and $T \ge 0$, suppose that:

$$(i) \quad \lim_{N \to \infty} \sup_{K \ge 1} \mathbb{E} \left[ \sup_{0 \le t \le T} \frac{1}{K_X} \sum_{k=1}^{K_X} \sum_{i=1}^{n_{\max}} \mathbf{1}_{\{n_k^X - s_k^X(0) \ge i, |\tau_{k,i}^X(0) - t| \ge N\}} \right] = 0.$$

$$(ii) \quad \text{For any } c \in \mathbb{R}, \text{ for any } i \in [1, n_{\max}],$$

$$\lim_{\epsilon \to 0} \sup_{K \ge 1} \mathbb{E} \left[ \frac{1}{K_X} \sum_{k=1}^{K_X} \mathbf{1}_{\{n_k^X - s_k^X(0) \ge i, |(\tau_{k,i}^X(0) - T) - c| \le \epsilon\}} \right] = 0.$$

Briefly, the first assumption allows to control the impact of the initial condition on the queues of the distribution of remaining infectious periods at each time, while the second condition is
related to aspects of absolute continuity. These conditions are for instance satisfied if for any $K \ge 1$, at time 0, the remaining infectious periods of infected individuals are i.i.d. of law $\nu$, while those of recovered individuals are set to be equal to zero. Notice that this choice for recovered individuals does not represent a loss of generality, as it does not affect the epidemic spread and initially recovered individuals will remain recognizable at any time $T$ as the only ones whose remaining infectious periods equal $-T$.

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

For any $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$, let $f_t(x) = f(t, x)$ and $f_t^I(x) = \langle \nu, f_t(j(x, \cdot)) \rangle$ for every $(t, x) \in \mathbb{R}_+ \times E$. 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 rate of within-structure infections at time $t$ is given by

$$ \lambda_X \langle \zeta_t^{X|K}, \mathbf{s}\mathbf{i} \rangle = \frac{\lambda_X}{K_X} \sum_{k=1}^{K_X} \mathbf{s}(x_k^X(t)) \mathbf{i}(x_k^X(t)). $$

Notice also that as mentioned previously, the average size of structures is constant over time, hence $\langle \zeta_t^{X|K}, \mathbf{n} \rangle = \langle \zeta_0^{X|K}, \mathbf{n} \rangle$ for all $t \ge 0$. Finally, let $\mathfrak{M}_1 = \mathcal{M}_1(E)^2$. We are now ready to state our first result, whose proof is postponed to Section 4.1.

**Theorem 3.2.** Let $\omega \in \Omega_G^*$. Suppose that $(\zeta_0^K)_{K \ge 1}$ satisfies Assumption 3.1 and converges in law to $\eta_0 \in \mathfrak{M}_1$. Then $(\zeta^K)_{K \ge 1}$ converges in $\mathbb{D}(\mathbb{R}_+, \mathcal{M}_1(E))^2$ to $\eta = (\eta^H, \eta^W)$ defined as the unique solution of the following system of Equations (6). For any $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$, for any $T \ge 0$,

$$ \begin{aligned} \langle \eta_T^X, f_T \rangle = \langle \eta_0^X, f_0 \rangle &+ \int_0^T \langle \eta_t^X, \mathcal{A}f_t \rangle dt + \lambda_X \int_0^T \langle \eta_t^X, \mathbf{s}\mathbf{i}(f_t^I - f_t) \rangle dt \\ &+ \lambda_{\bar{X}} \int_0^T \frac{\langle \bar{\eta}_t^X, \mathbf{s}\mathbf{i} \rangle}{\langle \bar{\eta}_t^X, \mathbf{s} \rangle} \langle \eta_t^X, \mathbf{s}(f_t^I - f_t) \rangle dt + \beta_G \int_0^T \frac{\langle \eta_0^H, \mathbf{i} \rangle}{\langle \eta_0^H, \mathbf{n} \rangle} \langle \eta_t^X, \mathbf{s}(f_t^I - f_t) \rangle dt. \end{aligned} \quad (6) $$

This measure-valued equation can further be related to a system of PDEs. Indeed, it is possible to establish an absolute continuity result for the marginals of $\eta^X$ conditioned on the structure's size and number of susceptible members. The associated densities can be shown to satisfy, in the sense of distributions, a system of differential equations related to non-linear nonlocal transport equations.

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 3.2 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 introduce the set

$$ \mathcal{S} = \{(n - i, i) : 2 \le n \le n_{\max}, 0 \le i \le n - 1\}. $$
For $(S, I) \in \mathbb{S}$, let $n_{S,I}^H(t)$ be the proportion of households containing $S$ susceptible and $I$ infected individuals at time $t$, according to distribution $\eta_t^H$. Define $n_{S,I}^W(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_{S,I}^X(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. Let us emphasize here that actually, only this second assumption is crucial for the results to hold, while the original distribution of infected individuals does not need to be uniform (in which case forthcoming Equations (7) and (9) need to be adapted). We have chosen this particular initial condition as it has been previously considered in the literature, and refer to the Discussion for further comments.

In practice, this setting corresponds to the following probability distribution $\eta_{0,\varepsilon} = (\eta_{0,\varepsilon}^H, \eta_{0,\varepsilon}^W) \in \mathfrak{M}_1$ characterized for $X \in \{H, W\}$ as follows. For any $n \in [1, n_{\max}]$ and $s \in [0, n]$:

$$ \eta_{0,\varepsilon}^X(n, s, d\tau) = \pi_n^X \binom{n}{s} (1-\varepsilon)^s \varepsilon^{n-s} \left( \nu^{\otimes(n-s)} \otimes \delta_0^{\otimes(n_{\max}-n+s)} \right) (d\tau). \quad (7) $$

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 4.2.

**Theorem 3.3.** Let $\varepsilon > 0$ and suppose that $\eta$ satisfies Equation (6) with $\eta_0 = \eta_{0,\varepsilon}$. Then the functions $(s, i, n_{S,I}^X : X \in \{H, W\}, (S, I) \in \mathbb{S})$ are characterized as being the unique solution of the following dynamical system: for any $t \ge 0$, $X \in \{H, W\}$ and $(S, I) \in \mathbb{S}$,

$$ \frac{d}{dt}s(t) = -(\tau_H(t) + \tau_W(t) + \tau_G(t)s(t)), \quad (8a) $$

$$ \frac{d}{dt}i(t) = -\frac{d}{dt}s(t) - \gamma i(t), \quad (8b) $$

$$ \begin{aligned} \frac{d}{dt} n_{S,I}^X(t) = & - \left( \lambda_X SI + \tau_{\bar{X}}(t) \frac{S}{s(t)} + \tau_G(t)S + \gamma I \right) n_{S,I}^X(t) \\ & + \gamma(I+1)n_{S,I+1}^X(t)\mathbf{1}_{\{S+I<n_{\max}\}} \\ & + \left( \lambda_X(S+1)(I-1) + \tau_{\bar{X}}(t) \frac{S+1}{s(t)} + \tau_G(t)(S+1) \right) n_{S+1,I-1}^X(t)\mathbf{1}_{\{I \ge 1\}}, \end{aligned} \quad (8c) $$

with initial conditions given by

$$ s(0) = 1 - \varepsilon; \quad i(0) = \varepsilon; \quad n_{S,I}^X(0) = \binom{S+I}{I} \pi_{S+I}^X (1-\varepsilon)^S \varepsilon^I. \quad (9) $$

This dynamical system may be understood as follows. Equation (8a) 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 (8b) is due to newly contaminated individuals moving from the susceptible to the infected state, which they in turn leave at rate $\gamma$. It remains to take an interest in Equation (8c). 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 which may occur in any layer of the graph, or upon the removal of one of its infected members. Simultaneously, a structure of type $(S, I+1)$ transforms into a structure of
<figure><img src="image_3.png" /><figcaption>Figure 1: Household and workplace size distributions $\pi^H$ (left) and $\pi^W$ (right) used in simulations.</figcaption></figure>

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 3.3, 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 3.2, 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 3.2 and 3.3, let us examine the latter from a numerical point of view.

## 3.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 (8a-c), 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. [42].

### 3.2.1 Implementation of the dynamical system and illustration of Theorem 3.3

Let us start by illustrating the result of Theorem 3.3 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 [5]. 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, and these proportions both equal 40% in the second scenario. Further, the epidemic is started by infecting either 10 or 100 individuals chosen uniformly at random at time 0. For each simulation, we have followed the evolution of the proportion of susceptible and infected individuals within the population over time.

The simulation outcomes are presented in Figure 2. For each choice of parameters, the solutions $s$ and $i$ of dynamical system (8a-c) with initial condition given by (9) are plotted on
the same graph. As expected, one observes good accordance of the stochastic simulations and the deterministic functions $(s, i)$.

Before proceeding further, let us briefly emphasize a few aspects of the implementation of the proposed deterministic model. A potential drawback of dynamical system (8a-c) consists in its large dimension. Indeed, it holds that $\#\mathcal{S} = n_{\max}(n_{\max} + 1)/2 - 1$. The number of equations of dynamical system (8a-c) is hence 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. We have compared the average time needed to either solve once dynamical system (8a-c), or to simulate one trajectory of the stochastic model using Gillespie's algorithm, for different choices of epidemic parameters. 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. 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.

### 3.2.2 Comparison to edge-based compartmental models

One may notice that the population structure, as described in Section 2.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 [30, 32], 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 [42]. Details are provided in Appendix B.

Edge-based compartmental models on CCM variants have been known to be in good accordance with simulations of the corresponding stochastic epidemic models, under the assumption of a very small initial proportion of infected. In our case, we have confronted the EBCM with dynamical system (8a-c), as well as simulated trajectories of our stochastic model. As expected, for very small values of $\varepsilon$, the EBCM and dynamical system (8a-c) 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 for details.

Finally, proceeding like before, we obtain that the number of equations of the EBCM is of order $O(n_{\max}^3)$. This has a strong negative impact on computation time, as briefly illustrated in Appendix B, 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 (8a-c), 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 (8a-c) seems more pertinent in the case of the epidemic model under consideration.
<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison of the stochastic model with its large population approximation given by dynamical system (8a-c). Household and workplace distributions are those of Figure 1. Two sets of epidemic parameters are 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 is used to simulate 50 trajectories of the stochastic model defined in Proposition 2.1 in a population of $K = 10000$ individuals (faint lines). For Panels (a) and (b), only trajectories reaching a threshold proportion of 0.005 infected are kept, and time is shifted so that time 0 corresponds to the moment when this threshold is reached. Finally, the deterministic solution $(s, i)$ of (8a-c) is represented for each scenario (thick lines). For Panels (a) and (b), the same time shifting procedure as for simulations is applied.</figcaption></figure>
# 4 Proofs

This section is devoted to establishing Theorems 3.2 and 3.3. As we will see, the proof of Theorem 3.2 has the intrinsic difficulty of all convergence results for measure-valued processes, with some technical difficulty arising from the infectiousness being a discontinuous function of the remaining infectious period. This will become apparent in the proof of forthcoming Proposition 4.8. Nevertheless, 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 3.3 shows that this actually is possible, the final reduced model taking the form of dynamical system (8a-c). 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 show that structuring the population of previously contaminated individuals by remaining infectious period is not necessary to handle the correlation of epidemic states of structures sharing a common infected. While this result fundamentally relies on the memory-less property of the exponential distribution, it will demand some effort, as illustrated in forthcoming Propositions 4.12 and 4.14.

## 4.1 Proof of Theorem 3.2

Let us start with the proof of Theorem 3.2. It follows a classical scheme, establishing tightness of $(\zeta^K)_{K \ge 1}$, whose limiting values are shown to satisfy Equation (6). Uniqueness of the solutions of this equation given the initial condition then ensures the desired convergence result. In particular, the proof is inspired by [14] and [39].

### 4.1.1 Uniqueness and continuity of the solution of Equation (6).

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

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

**Lemma 4.1.** Let $f \in C_b(E, \mathbb{R})$. There exists a sequence $(f_k)_{k \ge 1}$ taking values in $C_b^1(E, \mathbb{R})$ such that $f_k$ converges simply to $f$ and $\sup_{k \ge 1} \|f_k\|_\infty \le \|f\|_\infty$.

*Proof.* Consider a mollifier $\psi$, i.e. $\psi \in C^\infty(\mathbb{R}^{n_{\max}})$ is compactly supported, its mass $\int_{\mathbb{R}^{n_{\max}}} \psi(x)dx$ equals 1, and for $k \ge 1$, the function $\psi_k: x \mapsto k^{n_{\max}}\psi(kx)$ satisfies $\lim_{k \to \infty} \psi_k = \delta_0$ in the sense of distributions. Let $f \in C_b(E, \mathbb{R})$. Define the sequence $(f_k)_{k \ge 1}$ as follows:

$$f_k: (n, s, \tau) \in E \mapsto f(n, s, \cdot) * \psi_k(\tau).$$

Then, for any $x \in E$, by definition of $(\psi_k)_{k \ge 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$, $\|f_k\|_\infty \le \|f\|_\infty$. Finally, it also follows from the usual properties of convolution that for any $k \ge 1$, $f_k$ is smooth with respect to its last variable and the corresponding partial derivatives are bounded, hence $f_k \in C_b^1(E, \mathbb{R})$. $\square$

We may now turn to the main result of this paragraph. With slight abuse of notation, for an element $\eta = (\eta_1, \eta_2) \in \mathfrak{M}_1$, we define its total variation norm by $\|\eta\|_{TV} = \|\eta_1\|_{TV} \vee \|\eta_2\|_{TV}$.

**Proposition 4.2.** Let $\eta_* \in \mathfrak{M}_1$. Then Equation (6) admits at most one measure-valued solution $\eta$ which belongs to $C(\mathbb{R}_+, (\mathfrak{M}_1, \|\cdot\|_{TV}))$, such that $\eta_0 = \eta_*$. □
From now on, for $X \in \{H, W\}$, define $n_X = \langle \eta_0^X, \mathbf{n} \rangle$, $s_X(t) = \langle \eta_t^X, \mathbf{s} \rangle$ and $i_X(t) = \langle \eta_t^H, \mathbf{i} \rangle$. Let us establish the proposition.

*Proof.* First, notice that it follows immediately from Equation (6) that $\langle \eta_T^X, 1 \rangle = \langle \eta_0^X, 1 \rangle$ for $X \in \{H, W\}$, thus $\eta_0 \in \mathfrak{M}_1$ implies that for any $T \ge 0$, $\eta_T \in \mathfrak{M}_1$.

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

Let $T \ge 0$ and $g \in C_b^1(E, \mathbb{R})$ such that $\|g\|_\infty \le 1$. Consider the function defined by

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

and recall that $f_t^I(x) = \langle \nu, f_t(j(x, \cdot)) \rangle$. Then by definition, $f_T(x) = g(x)$. It follows from the assumption $g \in C_b^1(E, \mathbb{R})$ that $f \in C_b^1(\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

$$ \partial_t f_t(x) = \sum_{k=1}^{n-s} \partial_{\tau_k} g(\Psi(x, T, t)) \quad \text{and} \quad \partial_{\tau_k} f_t(x) = \partial_{\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_T^H, g \rangle = \langle \eta_T^H, f_T \rangle$, it follows from Equation (6) that

$$ \langle \eta_T^H, g \rangle = \langle \eta_0^H, f_0 \rangle + \lambda_H \int_0^T \langle \eta_t^H, \mathbf{si} (f_t^I - f_t) \rangle dt \\ + \lambda_W \int_0^T \frac{1}{s_W(t)} \langle \eta_t^W, \mathbf{si} \rangle \langle \eta_t^H, \mathbf{s} (f_t^I - f_t) \rangle dt + \beta_G \int_0^T \frac{i_H(t)}{n_H} \langle \eta_t^H, \mathbf{s} (f_t^I - f_t) \rangle dt. $$

Recall that $i_H(t) \le n_{\max}$ and $\frac{i_H(t)}{n_H} \le 1$ since for any $x \in E$, $\mathbf{i}(x) \le \mathbf{n}(x)$. We may notice that the following inequalities hold, as for any $t$, $\|f_t\|_\infty \le 1$:

$$ \begin{aligned} \langle \eta_t^H, \mathbf{si} (f_t^I - f_t) \rangle &\le 2(n_{\max})^2, & \langle \eta_t^H, \mathbf{s} (f_t^I - f_t) \rangle &\le 2n_{\max}, \\ \frac{1}{s_W(t)} \langle \eta_t^W, \mathbf{si} \rangle &\le n_{\max}, & \frac{i_H(t)}{n_H} \langle \eta_t^H, \mathbf{s} (f_t^I - f_t) \rangle &\le 2n_{\max}. \end{aligned} \quad (10) $$

Let $C = 2n_{\max}(\lambda_H n_{\max} + \lambda_W n_{\max} + \beta_G)$ and let $\epsilon \in \mathbb{R}$. It then follows from Inequalities (10) that

$$ |\langle \eta_T^H - \eta_{T+\epsilon}^H, g \rangle| \le C|T - (T + \epsilon)| = C|\epsilon|. \quad (11) $$

Consider now $h \in C_b(E, \mathbb{R})$ such that $\|h\|_\infty \le 1$. Lemma 4.1 ensures that there exists a sequence $(g_k)_{k \ge 1}$ taking values in $C_b^1(E, \mathbb{R})$ which converges simply to $h$ and such that $\|g_k\|_\infty \le 1$. By dominated convergence, this implies that

$$ |\langle \eta_T^H - \eta_{T+\epsilon}^H, h \rangle| \le C|T - (T + \epsilon)| = C|\epsilon|. $$

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

$$ \|\eta_T^H - \eta_{T+\epsilon}^H\|_{TV} = \sup_{h \in C_b(E, \mathbb{R}): \|h\|_\infty \le 1} |\langle \eta_T^H - \eta_{T+\epsilon}^H, h \rangle| \le C|\epsilon|. \quad (12) $$

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

It remains to establish uniqueness of the solution $\eta$ of Equation (6) 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, \bar{\eta}$ be two solutions of Equation (6) with initial condition $\eta_*$. Let $\bar{i}_H(t) = \langle \bar{\eta}_t^H, \mathbf{i} \rangle$ and define $\bar{s}_W(t)$ in analogous manner. As before, let $T \ge 0$. Consider again $g \in C_b^1(E, \mathbb{R})$ such that $\|g\|_\infty \le 1$, and define $f_t$ and $f_t^I$ as previously. Then

$$ 
\begin{aligned}
|\langle \eta_T^H - \bar{\eta}_T^H, g \rangle| &\le \lambda_H \int_0^T |\langle \eta_t^H - \bar{\eta}_t^H, \mathbf{s}\mathbf{i} (f_t^I - f_t) \rangle| dt \\
&\quad + \lambda_W \int_0^T \left| \frac{1}{s_W(t)} \langle \eta_t^W, \mathbf{s}\mathbf{i} \rangle \langle \eta_t^H, \mathbf{s}(f_t^I - f_t) \rangle - \frac{1}{\bar{s}_W(t)} \langle \bar{\eta}_t^W, \mathbf{s}\mathbf{i} \rangle \langle \bar{\eta}_t^H, \mathbf{s}(f_t^I - f_t) \rangle \right| dt \\
&\quad + \frac{\beta_G}{m_H} \int_0^T |i_H(t) \langle \eta_t^H, \mathbf{s}(f_t^I - f_t) \rangle - \bar{i}_H(t) \langle \bar{\eta}_t^H, \mathbf{s}(f_t^I - f_t) \rangle| dt.
\end{aligned}
$$

Proceeding similarly as in Inequalities (10), we obtain that

$$ 
\begin{aligned}
|\langle \eta_T^H - \bar{\eta}_T^H, g \rangle| &\le \frac{C}{2} \int_0^T |\langle \eta_t^H - \bar{\eta}_t^H, f_t^I - f_t \rangle| dt \\
&\le \frac{C}{2} \left( \int_0^T |\langle \eta_t^H - \bar{\eta}_t^H, f_t^I \rangle| dt + \int_0^T |\langle \eta_t^H - \bar{\eta}_t^H, f_t \rangle| dt \right).
\end{aligned}
$$

On the one hand, by definition, $f_t \in C_b(E, \mathbb{R})$ and $\|f_t\|_\infty \le 1$. Thus $|\langle \eta_t^H - \bar{\eta}_t^H, f_t \rangle| \le \|\eta_t^H - \bar{\eta}_t^H\|_{TV}$. On the other hand, it follows from the usual criterion of continuity for parametric integrals that $f_t^I$ is continuous on $E$ and $\|f_t^I\|_\infty \le \|f_t\|_\infty \le 1$. As a consequence, $|\langle \eta_t^H - \bar{\eta}_t^H, f_t^I \rangle| \le \|\eta_t^H - \bar{\eta}_t^H\|_{TV}$. Hence

$$ |\langle \eta_T^H - \bar{\eta}_T^H, g \rangle| \le C \int_0^T \|\eta_t^H - \bar{\eta}_t^H\|_{TV} dt. $$

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

$$ \|\eta_T^H - \bar{\eta}_T^H\|_{TV} \le C \int_0^T \|\eta_t^H - \bar{\eta}_t^H\|_{TV} dt. $$

Gronwall's lemma then assures that

$$ \forall t \in [0, T], \quad \|\eta_t^H - \bar{\eta}_t^H\|_{TV} = 0. $$

One obtains the analogous result for $\eta^W$ in the same manner. As $T \ge 0$ is arbitrary, this concludes the proof. $\square$

### 4.1.2 Tightness of $(\zeta^K)_{K \ge 1}$ in $D(\mathbb{R}_+, (M_F(E), w))^2$

Let us now turn to the tightness of $(\zeta^K)_{K \ge 1}$ in $D(\mathbb{R}_+, (M_F(E), w))^2$, where $w$ designates the weak topology on $M_F(E)$. We start by establishing the following preliminary result, whose proof relies on the chain rule.

**Lemma 4.3.** Let $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$. Then for any $T \ge t_0 \ge 0$, for any $x \in E$,

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

*Proof.* For $x = (n, s, \tau) \in E$ and $t_0 \in \mathbb{R}_+$, define

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

$$f_{n,s} : \mathbb{R}_+ \times \mathbb{R}^{n_{\max}} \to \mathbb{R}, \quad (u, \mathbf{v}) \mapsto f(u, (n, s, \mathbf{v}))$$
and $h_{t_0,x} : \mathbb{R}_+ \to \mathbb{R}^{1+n_{\max}}, \quad 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 \ge t_0$,

$$\frac{d}{dt} g_{t_0,x}(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 on the one hand, $\partial_1 f_{n,s}(u, \mathbf{v}) = \partial_t f(u, (n, s, \mathbf{v}))$ and on the other, $k \ge 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 \ge t_0$, $g_{t_0,x}$ is differentiable at $t$ and satisfies

$$\frac{d}{dt} g_{t_0,x}(t) = \mathcal{A} f(t, \Psi(x, t, t_0)).$$

This concludes the proof.

Throughout the section, we use the notation $\mathcal{S} = \{H, W, G\}$. Also, for any $f \in C_b^1(\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 \ge 0$ and $u = (\theta, k, \ell, \sigma) \in \bigcup_{Y \in \mathcal{S}} U_Y$:

$$g_{t,u}^H = g(t, j(x_k^H(t), \sigma)) - g(t, x_k^H(t)) \quad \text{and} \quad g_{t,u}^W = g(t, j(x_\ell^W(t), \sigma)) - g(t, x_\ell^W(t)).$$

**Proposition 4.4.** Consider $\zeta^K$ as introduced in Proposition 2.1. For any $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$, $T \ge 0$ and $X \in \{H, W\}$,

$$\langle \zeta_T^{X|K}, f_T \rangle = \langle \zeta_0^{X|K}, f_0 \rangle + \int_0^T \langle \zeta_t^{X|K}, \mathcal{A} f_t \rangle dt + \frac{1}{K_X} \sum_{Y \in \mathcal{S}} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) f_{t-,u}^X Q_Y(dt, du).$$

*Proof.* Let $f \in C_b^1(\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 \ge t \ge 0$ and $u = (\theta, k, \ell, \sigma) \in \bigcup_{Y \in \mathcal{S}} U_Y$:

$$\langle \Delta_H(u, T, t), g \rangle = g(\Psi(j(x_k^H(t-), \sigma), T, t)) - g(\Psi(x_k^H(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 (5), it follows that

$$\langle \zeta_T^{X|K}, f_T \rangle = \frac{1}{K_X} \left( \sum_{j=1}^{K_X} f_T \left( \Psi(x_j^X(0), T, 0) \right) + \sum_{Y \in \mathcal{S}} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) \langle \Delta_X(u, T, t), f_T \rangle Q_Y(dt, du) \right).$$

Using the result from Lemma 4.3, this becomes:

$$\langle \zeta_T^{X|K}, f_T \rangle = \frac{1}{K_X} \sum_{j=1}^{K_X} \left( f_0(x_j^X(0)) + \int_0^T \mathcal{A} f_t(\Psi(x_j^X(0), t, 0)) dt \right) \\ + \frac{1}{K_X} \sum_{Y \in \mathcal{S}} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) \left( \int_t^T \langle \Delta_X(u, z, t), \mathcal{A} f_z \rangle dz \right) Q_Y(dt, du) \\ + \frac{1}{K_X} \sum_{Y \in \mathcal{S}} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) f_{t-,u}^X Q_Y(dt, du).$$
It follows from the definition of $C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$ that both $f$ and $Af$ are bounded, hence we may apply Fubini's theorem to obtain that

$$ \begin{aligned} \langle \zeta_T^{X|K}, f_T \rangle = & \frac{1}{K_X} \sum_{j=1}^{K_X} f_0(x_j^X(0)) + \frac{1}{K_X} \sum_{Y \in S} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) f_{t-, u}^X Q_Y(dt, du) \\ & + \frac{1}{K_X} \int_0^T \left( \sum_{j=1}^{K_X} \mathcal{A}f_z(\Psi(x_j^X(0), z, 0)) + \sum_{Y \in S} \int_0^z \int_{U_Y} \mathcal{I}_Y(t-, u) \langle \Delta_X(u, z, t), \mathcal{A}f_z \rangle Q_Y(dt, du) \right) dz. \end{aligned} $$

The first sum on the right-hand side equals $\langle \zeta_0^{X|K}, f_0 \rangle$. From the second line, one recognizes in the integrand the definition of $\langle \zeta_z^{X|K}, \mathcal{A}f_z \rangle$ from Equation (5). This yields the desired result. $\square$

For $Y \in S$, let

$$ \tilde{Q}_Y(dt, du) = Q_Y(dt, du) - dt\mu_Y(du) $$

be the compensated martingale-measure associated to $Q_Y$.

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

$$ \langle \zeta_T^{X|K}, f_T \rangle = M_T^{X|K}(f) + V_T^{X|K}(f), $$

where we define

$$ M_T^{X|K}(f) = \frac{1}{K_X} \sum_{Y \in S} \int_0^T \int_{U_Y} \mathcal{I}_Y(t-, u) f_{t-, u}^X \tilde{Q}_Y(dt, du) $$

and

$$ V_T^{X|K}(f) = \langle \zeta_0^{X|K}, f_0 \rangle + \int_0^T \langle \zeta_t^{X|K}, \mathcal{A}f_t \rangle dt + \frac{1}{K_X} \sum_{Y \in S} \int_0^T \int_{U_Y} \mathcal{I}_Y(t, u) f_{t, u}^X \mu_Y(du) dt. $$

**Proposition 4.5.** Let $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$ and $X \in \{H, W\}$. Then $(M_T^{X|K}(f))_{T \ge 0}$ is a square integrable martingale. Using the same notations as in Theorem 3.2, its quadratic variation is given by

$$ \langle M^{X|K}(f) \rangle_T = \frac{1}{K_X} \int_0^T \langle \zeta_t^{X|K}, \mathcal{H}_t^X((f_t^2)^\mathcal{I} - 2f_t^\mathcal{I}f_t + f_t^2) \rangle dt, $$

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

$$ \mathcal{H}_t^X(x) = \beta_G \frac{I_H(t)}{N_H} \mathbf{s}(x) + \lambda_X \mathbf{s}(x) \mathbf{i}(x) + \lambda_X \frac{\langle \zeta_t^{\overline{X}|K}, \mathbf{s}\mathbf{i} \rangle}{S_{\overline{X}}(t)} \mathbf{s}(x). $$

*Proof.* Let $f \in C_b^1(E, \mathbb{R})$ and $X \in \{H, W\}$. Consider $M_T^{X|K}(f)$, which can be written as

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

where, for $Y \in S$ and $T \ge 0$,

$$ M_{Y,T}^{X|K}(f) = \int_0^T \int_{U_Y} \frac{1}{K_X} \mathcal{I}_Y(t-, u) f_{t, u}^X \tilde{Q}_Y(dt, du). $$

Suppose that for any $Y \in S$,

$$ \mathbb{E} \left[ \int_0^T \int_{U_Y} \left( \frac{1}{K_X} \mathcal{I}_Y(t, u) f_{t, u}^X \right)^2 \mu_Y(du) dt \right] < \infty, $$
then for all $Y \in S$, $(M_{Y,T}^{X|K}(f))_{T \ge 0}$ is as square integrable martingale [28], implying that $(M_T^{X|K}(f))_{T \ge 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 S} \langle M_Y^{X|K}(f) \rangle_T. $$

It thus is enough to study $(M_{Y,T}^{X|K}(f))_{T \ge 0}$ for all $Y \in 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_{\ell}^W(t) = K_W S_W(t)$, and that for any $k \in [1, K_H]$ and $t \in [0, T]$, $s_k^H(t)$ and $i_k^H(t)$ are less then $n_{\max}$, almost surely. Hence, replacing $S(t)$ by $K_W S_W(t)$ in $\mathcal{I}_H$,

$$ \begin{aligned} \mathbb{E} \left[ \langle M_H^{H|K}(f) \rangle_T \right] &= \mathbb{E} \left[ \int_0^T \int_{U_H} \left( \frac{1}{K_H} \mathcal{I}_H(t, u) f_{t,u}^H \right)^2 \mu_H(du) dt \right] \\ &= \mathbb{E} \left[ \int_0^T \frac{1}{K_H^2} \sum_{k=1}^{K_H} \lambda_H s_k^H(t) i(\tau_k^H(t)) \langle \nu, (f_t(j(x_k^H(t), \cdot)) - f_t(x_k^H(t)))^2 \rangle dt \right] \\ &\le \frac{1}{K_H} \lambda_H (n_{\max})^{24} \|f\|_{\infty}^2 T. \end{aligned} $$

Since further $K_H \ge K/n_{\max}$ and $\|f\|_{\infty}^2 < \infty$, we obtain that

$$ \mathbb{E} \left[ \langle M_H^{H|K}(f) \rangle_T \right] \le \frac{4}{K} \lambda_X (n_{\max})^3 \|f\|_{\infty}^2 T < \infty. \quad (13) $$

Thus $(M_{H,t}^{H|K}(f))_{T \ge 0}$ is a square integrable martingale whose quadratic variation is given by

$$ \begin{aligned} \langle M_H^{H|K}(f) \rangle_T &= \int_0^T \int_{U_H} \left( \frac{1}{K_H} \mathcal{I}_H(t, u) f_{t,u}^H \right)^2 \mu_H(du) dt \\ &= \frac{\lambda_H}{K_H} \int_0^T \langle \zeta_t^{H|K}, \mathbf{s}\dot{i}((f_t^2)^{\mathcal{I}} - 2f_t^{\mathcal{I}} f_t + f_t^2) \rangle dt, \end{aligned} $$

using the computations from Equation (4.1.2).

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

$$ \langle M_W^{H|K}(f) \rangle_T = \lambda_W \int_0^T \frac{\langle \zeta_t^{W|K}, \mathbf{s}\dot{i} \rangle}{S_W(t)} \langle \zeta_t^{H|K}, \mathbf{s}((f_t^2)^{\mathcal{I}} - 2f_t^{\mathcal{I}} f_t + f_t^2) \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 \le 1$ almost surely, we obtain that

$$ \mathbb{E} \left[ \langle M_G^{H|K}(f) \rangle_T \right] = \mathbb{E} \left[ \int_0^T \int_{U_G} \left( \frac{1}{K_H} \mathcal{I}_G(t, u) f_{t,u}^H \right)^2 \mu_G(du) dt \right] \le \frac{4}{K} \beta_G (n_{\max})^2 \|f\|_{\infty}^2 T. $$

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

$$ \langle M_G^{H|K}(f) \rangle_T = \frac{1}{K_H} \beta_G \int_0^T \frac{I_H(t)}{N_H} \langle \zeta_t^{H|K}, \mathbf{s}((f_t^2)^{\mathcal{I}} - 2f_t^{\mathcal{I}} f_t + f_t^2) \rangle dt. $$

This yields the desired result for $(M_T^{H|K}(f))_{T \ge 0}$. Proceeding in the same way for $(M_T^{W|K}(f))_{T \ge 0}$ concludes the proof. $\square$
We are now ready to focus on the tightness of $(\zeta^K)_{K \ge 1}$, endowing $\mathcal{M}_F(E)$ with the vague topology $v$ as a first step.

**Proposition 4.6.** Under the assumptions of Theorem 3.2, the sequence $(\zeta^K)_{K \ge 1}$ is tight in $D(\mathbb{R}_+, (\mathcal{M}_F(E), v))^2$.

The proof relies on the fact that in order to establish tightness of $(\zeta^K)_{K \ge 1}$, it is enough to show that for any $X \in \{H, W\}$, $(\langle \zeta^{X|K}, f \rangle)_{K \ge 1}$ is tight for a large enough set of test functions $f$ [35]. This in turn is ensured using the Aldous [1] and Rebolledo [20] 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 $(\zeta^{H|K})_{K \ge 1}$ and $(\zeta^{W|K})_{K \ge 1}$ are both tight in $D(\mathbb{R}_+, (\mathcal{M}_F(E), v))$.

Let us focus on $(\zeta^{H|K})_{K \ge 1}$. According to Theorem 2.1 of [35], 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_{\|x\|_{\infty} \to \infty} |f(x)| = 0 \right\}, $$

the sequence $(\langle \zeta^{H|K}, f \rangle)_{K \ge 1}$ is tight in $D(\mathbb{R}_+, \mathbb{R})$. Notice that by density of the set of smooth compactly supported functions in $C_0(\mathbb{R}^{n_{\max}})$ endowed with the uniform norm, it follows that $C_0(E, \mathbb{R}) \cap C_b^1(E, \mathbb{R})$ is also dense in $C_0(E, \mathbb{R})$ endowed with the uniform norm. Thus, let us consider $f \in C_0(E, \mathbb{R}) \cap C_b^1(E, \mathbb{R})$.

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

(i) For any $t$ belonging to a dense subset $\mathcal{T}$ of $\mathbb{R}^+$, $(\langle M^{H|K}(f) \rangle_t)_{K \ge 0}$ and $(V_t^{H|K}(f))_{K \ge 0}$ are tight in $\mathbb{R}$.

(ii) For any $T \ge 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 \ge 1}$ and $(T_K)_{K \ge 1}$ satisfying $S_K \le T_K \le T$ for all integers $K$,

$$ \sup_{K \ge K_0} \mathbb{P} \left( |\langle M^{H|K}(f) \rangle_{S_K} - \langle M^{H|K}(f) \rangle_{T_K}| \ge \alpha, T_K \le S_K + \delta \right) \le \epsilon $$

$$ \text{and } \sup_{K \ge K_0} \mathbb{P} \left( |V_{S_K}^{H|K}(f) - V_{T_K}^{H|K}(f)| \ge \alpha, T_K \le S_K + \delta \right) \le \epsilon. $$

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

$$ \sup_{K \ge 1} \mathbb{E} \left[ |\langle M^{H|K}(f) \rangle_t| \right] < \infty \text{ and } \sup_{K \ge 1} \mathbb{E} \left[ |V_t^{H|K}(f)| \right] < \infty. $$

Recalling that $C = 2n_{\max}(\lambda_H n_{\max} + \lambda_W n_{\max} + \beta_G)$, it follows from Equation (13) that

$$ \mathbb{E}[|\langle M^{H|K}(f) \rangle_t|] \le \frac{1}{K} 2n_{\max}C \|f\|_{\infty}^2 t. $$

Similar computations yield that

$$ \mathbb{E}[|V_t^{H|K}(f)|] \le \|f\|_{\infty} + \|Af\|_{\infty} + C \|f\|_{\infty} t. $$

As $f \in C_b^1(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 \ge 1}$ and $(T_K)_{K \ge 1}$ satisfying $S_K \le T_K \le T$ for all integers $K$. As previously, using Equation (13), we obtain the following upper bound:

$$ \mathbb{E} \left[ |\langle M^{H|K}(f) \rangle_{S_K} - \langle M^{H|K}(f) \rangle_{T_K}| \mid T_K \le S_K + \delta \right] \le \mathbb{E} \left[ \int_{S_K}^{T_K} dt \mid T_K \le S_K + \delta \right] \frac{2}{K} n_{\max} C \|f\|_{\infty}^2 \\ \le \frac{\delta}{K} 2 n_{\max} C \|f\|_{\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}| \ge \alpha, T_K \le S_K + \delta \right) \le \frac{\delta}{\alpha K} 2 n_{\max} C \|f\|_{\infty}^2. \quad (14) $$

Proceeding similarly, we also obtain that

$$ \mathbb{P} \left( |V_{S_K}^{H|K}(f) - V_{T_K}^{H|K}(f)| \ge \alpha, T_K \le S_K + \delta \right) \le \frac{\delta}{\alpha} (\|Af\|_{\infty} + C \|f\|_{\infty}). \quad (15) $$

Equations (14) and (15) 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. $\square$

Finally, this result on the tightness of $(\zeta^K)_{K \ge 1}$ in $\mathbb{D}(\mathbb{R}_+, (\mathcal{M}_F(E), v))^2$ lets us establish the main result of this subsection.

**Proposition 4.7.** Under the assumptions of Theorem 3.2, the sequence $(\zeta^K)_{K \ge 1}$ is tight in $\mathbb{D}(\mathbb{R}_+, (\mathcal{M}_F(E), w))^2$.

*Proof.* Let $X \in \{H, W\}$. Tightness in $\mathbb{D}(\mathbb{R}_+, (\mathcal{M}_F(E), w))$ of $(\zeta^{X|K})_{K \ge 1}$ will be shown using Theorem 1.1.8 from [40], which we state in our setting for the sake of completeness. Let $\Phi: z \in \mathbb{R} \mapsto 6z^2 - 15z^4 + 10z^3$ and for $N \ge 1$, define smooth approximations of $x \in E \mapsto \mathbf{1}_{\{\|\tau(x)\|_{\infty} \ge N\}}$ by:

$$ \forall x \in E, \forall N \ge 1, \phi_N(x) = \Phi(0 \lor (\|\tau(x)\|_{\infty} - (N-1)) \land 1). $$

Then in order to ensure the tightness of $(\zeta^{X|K})_{K \ge 1}$ in $\mathbb{D}(\mathbb{R}_+, (\mathcal{M}_F(E), w))$, it is sufficient to show that for any $T \ge 0$, the following conditions hold:

(i) There exists a family of functions $F$ which is dense in $C_0(E, \mathbb{R})$ and stable under addition, such that for any $f \in F \cup \{x \in E \mapsto 1\}$, the sequence $(\langle\zeta^{X|K}, f\rangle)_{K \ge 1}$ is tight in $\mathbb{D}(\mathbb{R}_+, \mathbb{R})$.

$$ \lim_{N \to \infty} \limsup_{K \to \infty} \mathbb{E}[\sup_{t \le T} \langle \zeta_t^{X|K}, \phi_N \rangle] = 0. $$

(iii) Any limiting value of $(\zeta^{X|K})_{K \ge 1}$, if it exists, belongs to $C([0, T], (\mathcal{M}_F(E), w))$.

The proof hence consists in checking those assumptions. Let $T \ge 0$, and consider $X = H$, as the case $X = W$ can be treated similarly. We may see that (i) is satisfied, as we have shown in the proof of Proposition 4.6 that for any $f \in C_0(E, \mathbb{R}) \cap C_b^1(E, \mathbb{R})$, $(\langle\zeta^{H|K}, f\rangle)_{K \ge 1}$ is tight, and further for any $K \ge 1$, for any $T \ge 0$, $\langle\zeta_T^{H|K}, 1\rangle = 1$ almost surely.

Let us now turn our attention to (ii). Start by noticing that for any $N \ge 1$ and $x \in E$,

$$ \phi_N(x) \le \mathbf{1}_{\{\|\tau(x)\|_{\infty} \ge N-1\}} \le f_{N-1}(x) := \sum_{i=1}^{n_{\max}} f_{N-1,i}(x) $$

where $f_{N-1,i}(x) = \mathbf{1}_{\{n(x)-s(x) \ge i, |\tau_i(x)| \ge N-1\}}$.
Let $t \in [0, T]$. For any $N-1 \ge t$, $x \in E$, $z \in [0, t]$ and $\sigma \ge 0$, it holds by definition that

$$f_{N-1,i}(\Psi(j(x, \sigma), t, z)) - f_{N-1,i}(\Psi(x, t, z)) = \mathbf{1}_{\{n(x)-s(x)=i-1, |\sigma-(t-z)|>N-1\}} \le \mathbf{1}_{\{\sigma>N-1\}}.$$

Hence, using Proposition 2.1 and the above upper bounds, it follows that almost surely,

$$\langle \zeta_t^{H|K}, \phi_N \rangle \le \frac{1}{K_H} \sum_{k=1}^{K_H} f_{N-1}(\Psi(x_k^H(0), t, 0)) + \frac{n_{\max}}{K_H} \sum_{Y \in S} \int_0^t \int_{U_Y} \mathcal{I}_Y(z-, u) \mathbf{1}_{\{\sigma>N\}} Q_Y^K(dz, du).$$

Defining as previously $C = 2n_{\max}(\lambda_H n_{\max} + \lambda_W n_{\max} + \beta_G)$, this leads to the following upper bound:

$$\mathbb{E}[\sup_{t \le T} \langle \zeta_t^{H|K}, \phi_N \rangle] \le \mathbb{E}[\sup_{t \le T} \frac{1}{K_H} \sum_{k=1}^{K_H} f_{N-1}(\Psi(x_k^H(0), t, 0))] + \frac{CT}{2} n_{\max} \nu((N-1, +\infty)).$$

As a consequence, Assumption 3.1 (i) ensures that (ii) is satisfied. Notice that this assumption could actually be a little bit relaxed here, as it would be enough if the supremum over $K$ were replaced by the limit superior over $K \to \infty$.

In order to check that condition (iii) holds, we will follow the arguments presented in [21]. Suppose that $\eta^H$ is a limiting value of $(\zeta^{H|K})_{K \ge 1}$. By definition,

$$\sup_{t \in [0, T]} \sup_{f \in L^\infty, \|f\|_\infty \le 1} |\langle \zeta_t^{H|K}, f \rangle - \langle \zeta_{t-}^{H|K}, f \rangle| \le \frac{1}{K_H}.$$

As the application $\mu \mapsto \sup_{t \in [0, T]} |\langle \mu_t, f \rangle - \langle \mu_{t-}, f \rangle|$ is continuous on $\mathbb{D}([0, T], (\mathcal{M}_F(E), v))$ for any $f$ in a measure-determining countable set, it follows that $\eta^H \in \mathcal{C}([0, T], (\mathcal{M}_F(E), v))$.

Let us now introduce $\phi_{N,M} = \phi_N(1 - \phi_M)$, which serves as a smooth and compactly supported approximation of $x \in E \mapsto \mathbf{1}_{\{N \le \|\tau(x)\|_\infty \le M\}}$, for $N \le M$. As, on the one hand, $\mu \mapsto \sup_{t \in [0, T]} \langle \mu_t, \phi_{N,M} \rangle$ is continuous on $\mathbb{D}([0, T], (\mathcal{M}_F(E), v))$, and on the other hand, for any $K \ge 1$, $\sup_{t \in [0, T]} \langle \zeta_t^{X|K}, \phi_{N,M} \rangle \le 1$, it follows that:

$$\mathbb{E}[\sup_{t \in [0, T]} \langle \eta_t^H, \phi_{N,M} \rangle] = \lim_{K \to \infty} \mathbb{E}[\sup_{t \in [0, T]} \langle \zeta_t^{X|K}, \phi_{N,M} \rangle] \le \limsup_{K \to \infty} \mathbb{E}[\sup_{t \in [0, T]} \langle \zeta_t^{X|K}, \phi_N \rangle].$$

Letting $M$ go to infinity in the left hand side, dominated convergence ensures that

$$\mathbb{E}[\sup_{t \in [0, T]} \langle \eta_t^H, \phi_N \rangle] \le \limsup_{K \to \infty} \mathbb{E}[\sup_{t \in [0, T]} \langle \zeta_t^{X|K}, \phi_N \rangle] \xrightarrow{N \to \infty} 0,$$

where the convergence of the right hand side is achieved as in the proof of (ii). In particular, it thus is possible to extract a subsequence from $(\sup_{t \in [0, T]} \langle \eta_t^H, \phi_N \rangle)_N$ which converges almost surely to zero. This implies that for any $\epsilon$, there exists $N$ such that almost surely,

$$1 - \sup_{t \in [0, T]} \langle \eta_t^H, \mathbf{1}_{\{\|\tau(\cdot)\|_\infty \le N\}} \rangle \le \sup_{t \in [0, T]} \langle \eta_t^H, \phi_N \rangle < \epsilon$$

Thus $(\eta_t^H)_{t \in [0, T]}$ is almost surely tight.

Let $g \in \mathcal{C}_b(E)$, and let $g_N = g(1 - \phi_N)$. It then holds that for $h$ small so that $t + h \in [0, T]$,

$$|\langle \eta_{t+h}^H, g \rangle - \langle \eta_t^H, g \rangle| \le |\langle \eta_{t+h}^H, g - g_N \rangle| + |\langle \eta_{t+h}^H, g_N \rangle - \langle \eta_t^H, g_N \rangle| + |\langle \eta_t^H, g - g^N \rangle|.$$

Let $\epsilon > 0$. As $|g - g_N| \le \|g\|_\infty \phi_N$, there exists $N_0$ such that $\sup_{t \in [0, T]} \langle \eta_t^H, g - g_{N_0} \rangle < \epsilon/3$. Further, as $\eta^H \in \mathcal{C}([0, T], (\mathcal{M}_F(E), v))$, for $h$ small enough, $|\langle \eta_{t+h}^H, g_{N_0} \rangle - \langle \eta_t^H, g_{N_0} \rangle| < \epsilon/3$. This allows to conclude that $\eta^H \in \mathcal{C}([0, T], (\mathcal{M}_F(E), w))$, establishing (iii) and finally tightness of $\eta^H$ in $\mathbb{D}(\mathbb{R}_+, (\mathcal{M}_F(E), w))$. $\square$
### 4.1.3 Identification of the limiting values of $(\zeta^K)_{K \ge 1}$

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

**Proposition 4.8.** Under the assumptions of Theorem 3.2, all limiting values of $(\zeta^K)_{K \ge 1}$ in $D(\mathbb{R}_+, (\mathcal{M}_F(E), w))^2$ are continuous with regard to the total variation norm, and solutions of Equation (6).

Before proceeding to the proof of this proposition, let us emphasize that there is some technical difficulty due to infectiousness being a discontinuous function of an individual's remaining infectious period. Indeed, tightness of $(\zeta^K)_{K \ge 1}$ in $D(\mathbb{R}_+, (\mathcal{M}_F(E), w))^2$ allows us to extract a subsequence $(\zeta^{\varphi^{(K)}})_{K \ge 1}$ which converges in law in this space to some limiting value $\eta$, and our aim is to show that $\eta$ satisfies Equation (6). However, convergence in law in $D(\mathbb{R}_+, (\mathcal{M}_F(E), w))$ is not enough to ensure that $\langle \eta_t^X, \mathbf{i}f \rangle$ converges in law to $\langle \eta_t^X, \mathbf{i}f \rangle$ for $f \in C_b^1(E)$, as $\mathbf{i}$ is discontinuous on $E$. This leads to forthcoming Proposition 4.9.

*Proof.* Consider a subsequence $(\zeta^{\varphi^{(K)}})_{K \ge 1}$ of $(\zeta^K)_{K \ge 1}$ which converges in law in the space $D(\mathbb{R}_+, (\mathcal{M}_F(E), w))^2$, and let $\eta$ be its limit.

Notice that it follows from the Proof of Proposition 4.7 that $\eta \in C(\mathbb{R}_+, (\mathcal{M}_F(E), w))^2$ almost surely. Hence, following Proposition A.6.1 of [39], for any $X \in \{H, W\}$,

$$ \|\eta_T^X - \eta_{T-}^X\|_{TV} = \sup_{f \in C_b(E, \mathbb{R}): \|f\|_\infty \le 1} |\langle \eta_T^X, f \rangle - \langle \eta_{T-}^X, f \rangle| = 0 \quad \text{almost surely.} $$

It remains to show that $\eta$ satisfies Equation (6). Let $T \ge 0$ and $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$, and consider the application $\psi_T^H$ defined by

$$ \begin{aligned} \psi_T^H(\eta) = & \langle \eta_T^H, f_T \rangle - \langle \eta_0^H, f_0 \rangle - \int_0^T \langle \eta_t^H, A f_t \rangle dt - \int_0^T \langle \eta_t^H, \lambda_H \mathbf{s} \langle \mathbf{i}(f_t^I - f_t) \rangle dt \\ & - \int_0^T \frac{1}{\langle \eta_t^W, \mathbf{s} \rangle} \langle \eta_t^W, \lambda_W \mathbf{s} \rangle \langle \eta_t^H, \mathbf{s}(f_t^I - f_t) \rangle dt - \beta_G \int_0^T \frac{\langle \eta_t^H, \mathbf{i} \rangle}{\langle \eta_0^H, \mathbf{n} \rangle} \langle \eta_t^H, \mathbf{s}(f_t^I - f_t) \rangle dt. \end{aligned} \quad (16) $$

Start by noticing that $\psi_T^H(\zeta^{\varphi^{(K)}}) = M_T^{H|\varphi^{(K)}}(f)$, as $K_H S_H(t) = K_W S_W(t)$. Using Jensen's inequality, it follows from Equation (13) that

$$ \mathbb{E}[|\psi_T^H(\zeta^{\varphi^{(K)}})|]^2 \le \mathbb{E}[|\psi_T^H(\zeta^{\varphi^{(K)}})|^2] = \mathbb{E}[\langle M^{H|K}(f) \rangle_t] \le \frac{1}{K} 2n_{\max} C \|f\|_{\infty}^2 T \underset{K \to \infty}{\longrightarrow} 0. $$

Suppose that $(\psi_T^H(\zeta^{\varphi^{(K)}}))_{K \ge 1}$ converges in law to $\psi_T^H(\zeta)$. According to Theorem 3.5 of [7], it then is enough to proof that $(\psi_T^H(\zeta^{\varphi^{(K)}}))_{K \ge 1}$ is uniformly integrable to obtain that its expectation converges to the expectation of $\psi_T^H(\zeta)$. In our case, uniform integrability is easily assured as the sequence $(\psi_T^H(\zeta^{\varphi^{(K)}}))_{K \ge 1}$ is bounded. Indeed, using the fact that for all $T \ge 0$, $\zeta_T^{\varphi^{(K)}} \in \mathfrak{M}_1$, we obtain from Equation (16) that

$$ |\psi_T^H(\zeta^{\varphi^{(K)}})| \le ((2 + C) \|f\|_{\infty} + \|Af\|_{\infty}) T. $$

We may now conclude that

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

which yields the desired result.
It thus suffices to show that $(\psi_T^H(\zeta^{\varphi^{(K)}}))_{K \ge 1}$ converges in law to $\psi_T^H(\zeta)$. According to Skorokhod's representation theorem, there exists a probability space $\Omega$ on which one may define $(\tilde{\zeta}^K)_{K \ge 1}$ and $\tilde{\eta}$ equal in law to $(\zeta^{\varphi^{(K)}})_{K \ge 1}$ and $\eta$, respectively, such that $(\tilde{\zeta}^K)_{K \ge 1}$ converges almost surely in $D(\mathbb{R}_+, (\mathcal{M}_F(E), w))^2$ to $\tilde{\eta}$ on $\Omega$. In particular, it holds that

$$ \forall T \ge 0, \forall X \in \{H, W\}, \forall g \in \mathcal{C}_b(E), \quad \langle \tilde{\zeta}_T^{X|K}, g \rangle \xrightarrow{K \to \infty} \langle \tilde{\eta}_T^X, g \rangle \text{ almost surely.} $$

It follows immediately that for any $t \in [0, T]$ and $X \in \{H, W\}$, almost surely,

$$ (\langle \tilde{\zeta}_t^{X|K}, f_t \rangle, \langle \tilde{\zeta}_t^{X|K}, A f_t \rangle) \xrightarrow{K \to \infty} (\langle \tilde{\eta}_t^X, f_t \rangle, \langle \tilde{\eta}_t^X, A f_t \rangle). \quad (17) $$

Since $|\langle \tilde{\zeta}_t^{X|K}, A f_t \rangle| \le 2 \|Df\|_{\infty}$, where $Df$ designates the differential of $f$, dominated convergence ensures that

$$ \int_0^T \langle \tilde{\zeta}_t^{X|K}, A f_t \rangle dt \xrightarrow{K \to \infty} \int_0^T \langle \tilde{\eta}_t^X, A f_t \rangle dt \text{ almost surely.} \quad (18) $$

In order to establish the desired convergence of the last three terms of $\psi_T^H(\zeta)$, we will make use of the following proposition, whose proof is postponed. In this context, a $d$-dimensional rectangle is a set defined as the product of $d$ intervals of $\mathbb{R} \cup \{-\infty, +\infty\}$.

**Proposition 4.9.** For any $n \in [1, n_{\max}]$, for any $s \in [0, n]$, consider $m = m(n, s) < \infty$, a set $(A_k^{n,s})_{k \le m}$ of pairwise disjoint $(n-s)$-dimensional rectangles and a set $(\varphi_k^{n,s})_{k \le m}$ of functions belonging to $C_b^1(\mathbb{R}^{n-s})$. For any $\tau \in \mathbb{R}^{n_{\max}}$, let $\tau_{1,n-s} = (\tau_1, \dots, \tau_{n-s})$. Define the function $\phi: E \to \mathbb{R}$ by

$$ \forall x = (n, s, \tau) \in E, \quad \phi(x) = \sum_{k=1}^{m(n,s)} \mathbf{1}_{A_k^{n,s}}(\tau_{1,n-s}) \varphi_k^{n-s}(\tau_{1,n-s}). $$

Then for any $X \in \{H, W\}$ and $T \ge 0$, it holds that

$$ \langle \tilde{\zeta}_T^{X|K}, \phi \rangle \xrightarrow{K \to \infty} \langle \tilde{\eta}_T^X, \phi \rangle \text{ in } L^1. $$

Let us focus on the second-to-last term, representing infection events occurring within work-places, as the other two can be treated similarly.

The application $\phi(x) = \lambda_W \mathbf{s}(x)\mathbf{i}(x)$ is of the form described in Proposition 4.9, hence for any $t \in [0, T]$, $\langle \tilde{\zeta}_t^{W|K}, \lambda_W \mathbf{s}\mathbf{i} \rangle$ converges in $L^1$ to $\langle \tilde{\eta}_t^W, \lambda_W \mathbf{s}\mathbf{i} \rangle$ as $K$ tends to infinity. Also, notice that as $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$, it follows that for any $t \in [0, T]$, $f_t^I \in C_b^1(E)$. Thus Proposition 4.9 ensures that for any $t \in [0, T]$, $\langle \tilde{\zeta}_t^{H|K}, \mathbf{s}(f_t^I - f_t) \rangle$ converges in $L^1$ to $\langle \tilde{\eta}_t^H, \mathbf{s}(f_t^I - f_t) \rangle$ as $K$ tends to infinity.

In particular, the following convergence holds in probability:

$$ X_t^K := \langle \tilde{\zeta}_t^{W|K}, \lambda_W \mathbf{s}\mathbf{i} \rangle \langle \tilde{\zeta}_t^{H|K}, \mathbf{s}(f_t^I - f_t) \rangle \xrightarrow{K \to \infty} X_t := \langle \tilde{\eta}_t^W, \lambda_W \mathbf{s}\mathbf{i} \rangle \langle \tilde{\eta}_t^H, \mathbf{s}(f_t^I - f_t) \rangle. $$

Letting $c = 2\lambda_W n_{\max}^2 \|f\|_{\infty}$ and $D = \{(x, y) : |x| \le cy^2\}$, then $(X_t^K, \langle \tilde{\zeta}_t^{W|K}, \mathbf{s} \rangle)$ and $(X_t, \langle \tilde{\eta}_t^W, \mathbf{s} \rangle)$ belong almost surely to $D$, for any $K \ge 1$. As $\langle \tilde{\zeta}_t^{W|K}, \mathbf{s} \rangle$ converges almost surely to $\langle \tilde{\eta}_t^W, \mathbf{s} \rangle$, and the application $(x, y) \mapsto (x/y)\mathbf{1}_{\{y \ne 0\}}$ is continuous on $D$, we deduce the following convergence in probability:

$$ Y_t^K := \frac{\langle \tilde{\zeta}_t^{W|K}, \lambda_W \mathbf{s}\mathbf{i} \rangle}{\langle \tilde{\zeta}_t^{W|K}, \mathbf{s} \rangle} \langle \tilde{\zeta}_t^{H|K}, \mathbf{s}(f_t^I - f_t) \rangle \xrightarrow{K \to \infty} Y_t := \frac{\langle \tilde{\eta}_t^W, \lambda_W \mathbf{s}\mathbf{i} \rangle}{\langle \tilde{\eta}_t^W, \mathbf{s} \rangle} \langle \tilde{\eta}_t^H, \mathbf{s}(f_t^I - f_t) \rangle. $$
In addition, for any $K \ge 1$ and any $t \in [0, T]$, $|Y_t^K| \le 2\lambda_W n_{\max} \|f\|_{\infty}$. Thus using twice dominated convergence, we first obtain that the above convergence of $(Y_t^K)_{K \ge 1}$ to $Y_t$ also holds in $L^1$, and subsequently the following convergence holds in $L^1$:

$$ \int_0^T \frac{\langle \tilde{\zeta}_t^{W|K}, \lambda_W \mathbf{s} \rangle}{\langle \tilde{\zeta}_t^{W|K}, \mathbf{s} \rangle} \langle \tilde{\zeta}_t^H|K, \mathbf{s}(f_t^I - f_t) \rangle dt \xrightarrow{K \to \infty} \int_0^T \frac{\langle \tilde{\eta}_t^W, \lambda_W \mathbf{s} \rangle}{\langle \tilde{\eta}_t^W, \mathbf{s} \rangle} \langle \tilde{\eta}_t^H, \mathbf{s}(f_t^I - f_t) \rangle dt. \quad (19) $$

Reasoning in a similar manner, one also obtains:

$$ \int_0^T \frac{\langle \tilde{\zeta}_t^H, \lambda_H \mathbf{s} \rangle (f_t^I - f_t) dt}{\langle \tilde{\zeta}_t^H, \mathbf{i} \rangle \langle \tilde{\zeta}_t^H, \mathbf{n} \rangle} \xrightarrow{K \to \infty} \int_0^T \frac{\langle \tilde{\eta}_t^H, \lambda_H \mathbf{s} \rangle (f_t^I - f_t) dt}{\langle \tilde{\eta}_t^H, \mathbf{n} \rangle} \langle \tilde{\eta}_t^H, \mathbf{s}(f_t^I - f_t) dt \quad \text{in } L^1. \quad (20) $$

Thus, Equations (17-20) imply that all the terms on the right hand side of the definition of $\psi_T^H(\tilde{\zeta}^K)$ as stated in Equation (16) converge in probability, and thus their linear combination converges in probability to the linear combination of their limits. In other words, $\psi_T^H(\tilde{\zeta}^K)$ converges in probability to $\psi_T^H(\tilde{\eta})$, which ensures as desired that $\psi_T^H(\zeta^K)$ converges in law to $\psi_T^H(\eta)$. This concludes the proof. $\square$

In order to conclude, we only need to show that Proposition 4.9 holds.

**Proof of Proposition 4.9. Step 1.** Recall that a $d$-dimensional rectangle is a set $A$ defined as the product of $d$ intervals of $\mathbb{R} \cup \{-\infty, +\infty\}$. If all $d$ intervals are included in $\mathbb{R}$, the rectangle will further be said finite.

Let $X \in \{H, W\}, n \in [1, n_{\max}]$ and $s \in [0, n]$. We start by showing that for any $T \ge 0$ and any finite $(n-s)$-dimensional rectangle $B$,

$$ \langle \tilde{\zeta}_T^{X|K}, \mathbf{1}_{\{(n,s)\}\times B} \rangle \xrightarrow{K \to \infty} \langle \tilde{\eta}_T^X, \mathbf{1}_{\{(n,s)\}\times B} \rangle \quad \text{almost surely.} \quad (21) $$

As $\tilde{\eta}^X \in \mathcal{C}(\mathbb{R}_+, (M_F(E), w))$, it follows that for any $T$, $\tilde{\zeta}_T^{X|K}$ converges almost surely to $\tilde{\eta}_T^X$ in $(M_F(E), w)$. Thus, in order to establish the desired result, it is sufficient to show that $\mathbf{B} = \{(n, s)\} \times B$ is a $\tilde{\eta}_T^X$-continuity set, in which case the Portmanteau theorem allows to conclude.

For any set $A$, let $\partial A$ be the boundary of $A$. Then $\partial \mathbf{B} = \{(n, s)\} \times \partial B$. As $B$ is a $(n-s)$-dimensional rectangle, there exist $a_i < b_i \in \mathbb{R}$ for $i \in [1, n-s]$ such that $B$ can be written as the product of intervals (potentially open, closed or half-open) delimited by $a_i < b_i$, for $i \in [1, n-s]$. Thus

$$ \partial B = \bigcup_{i=1}^{n-s} \bigcup_{c \in \{a_i, b_i\}} \left( \prod_{j=1}^{i-1} [a_j, b_j] \times \{c\} \times \prod_{k=i+1}^{n-s} [a_k, b_k] \right). $$

Consider any $i \in [1, n-s]$ and $c \in \mathbb{R}$. We are going to prove that

$$ \langle \tilde{\eta}_T^X, \mathbf{1}_{\{(n,s)\}} \times (\prod_{j=1}^{i-1} [a_j, b_j] \times \{c\} \times \prod_{k=i+1}^{n-s} [a_k, b_k]) \rangle = 0, \quad (22) $$

which will be enough to conclude. In order to achieve this, let us introduce a mollifier $\psi \in \mathcal{C}^{\infty}(\mathbb{R})$ in the same sense as in the proof of Lemma 4.1, with compact support in $[-1, 1]$. For $\varepsilon > 0$, define the function $\varphi_{\varepsilon} : x \mapsto \varepsilon^{-1}\psi(x/\varepsilon)$, whose support lies in $[-\varepsilon, \varepsilon]$ and which converges to $\delta_0$ in the sense of distributions, when $\varepsilon$ goes to zero. For any $x = (n, s, \tau) \in E$, let

$$ \phi_{\varepsilon}(x) = \mathbf{1}_{\{\mathbf{n}(x)=n, \mathbf{s}(x)=s\}} \left( \prod_{\substack{j=1 \\ j \ne i}}^{n-s} \mathbf{1}_{[a_j, b_j]} * \varphi_{\varepsilon}(\tau_j) \right) \mathbf{1}_c * \varphi_{\varepsilon}(\tau_i). $$
As $\phi_{\varepsilon} \in C_b^1(E)$ and $\tilde{\zeta}_T^{X|K}$ converges almost surely to $\tilde{\eta}_T^X$ in $(M_F(E), w)$, dominated convergence implies that

$$ \mathbb{E}[\langle \tilde{\eta}_T^X, \phi_\varepsilon \rangle] = \lim_{K \to \infty} \mathbb{E}[\langle \tilde{\zeta}_T^{X|K}, \phi_\varepsilon \rangle]. $$

Notice that

$$ \phi_{\varepsilon}(x) \le \mathbf{1}_{\{\mathbf{n}(x)-\mathbf{s}(x)>i\}} \mathbf{1}_c * \varphi_{\varepsilon}(\tau_i(x)). $$

Hence proceeding as in the proof of Proposition 4.7, it follows that

$$ \mathbb{E}[\langle \tilde{\zeta}_T^{X|K}, \phi_\varepsilon \rangle] \le \mathbb{E} \left[ \frac{1}{K_X} \sum_{k=1}^{K_X} \mathbf{1}_{\{n_k^X - s_k^X(0) \ge i, |(\tau_{k,i}^X(0)-T)-c| \le \epsilon\}} \right] \\ + \frac{C}{2} \int_0^T \nu([c+(T-t)-\varepsilon, c+(T-t)+\varepsilon])dt. $$

Absolute continuity of $\nu$ with regard to the Lebesgue measure and Assumption 3.1 ensure that the right hand side is dominated by a function $c(\varepsilon)$ which does not depend on $K$, and which goes to zero with $\varepsilon$. Thus $\mathbb{E}[\langle \tilde{\eta}_T^X, \phi_\varepsilon \rangle] \le c(\varepsilon)$. In particular, one may construct a sequence $(\varepsilon_n)_{n \ge 1}$ which converges to 0 and satisfies $\sum_{n \ge 0} c(\varepsilon_n) < \infty$. Then on the one hand, the Borel-Cantelli lemma ensures that $\langle \tilde{\eta}_T^X, \phi_\varepsilon \rangle$ converges almost surely to 0 as $n$ tends to infinity. On the other hand, by dominated convergence, $\langle \tilde{\eta}_T^X, \phi_\varepsilon \rangle$ converges almost surely to the left-hand side of Equation (22) as $\varepsilon$ tends to zero, hence Equation (22) is proven to be true.

As a consequence, we conclude that $\langle \tilde{\eta}_T^X, \partial \mathbf{B} \rangle = 0$, and thus Equation (21) holds.

**Step 2.** Consider now a function $\phi$ as described in the proposition. For any $N$, let us introduce a partition of $\mathbb{R}^{n-s}$ whose elements consist in (partially open) hypercubes of side length $2^{-N}$. For every $k \le m(n, s)$, a partition $(B_{k,j}^N)_{j \ge 1}$ of $A_k^{n,s}$ is obtained by taking the intersection of $A_k^{n,s}$ with those hypercubes. As $A_k^{n,s}$ is a rectangle itself, the family $(B_{k,j}^N)_{j \ge 1}$ consists of rectangles of side length at most $2^{-N}$. For every $j$, consider a point $z_{k,j}^N$ belonging to $B_{k,j}^N$. Finally, define the set $J_N(k) = \{j \ge 0 : \sup\{\|x\|_{\infty} : x \in B_{k,j}^N\} \le N\}$, which contains only a finite number of elements. Then we can define the following approximation of $\phi$:

$$ \forall x \in E, \quad \phi_N(x) = \sum_{\substack{1 \le n \le n_{\max} \\ 0 \le s \le n}} \mathbf{1}_{\{\mathbf{n}(x)=n, \mathbf{s}(x)=s\}} \sum_{k=1}^{m(n,s)} \sum_{j \in J_N(k)} \varphi_k^{n-s}(z_{k,j}^N) \mathbf{1}_{B_{k,j}^N}(\tau_{1,n-s}(x)). $$

Using our result from the first step, for every $(n, s)$ such that $1 \le n \le n_{\max}$ and $0 \le s \le n$, for every $k \le m(n, s)$ and $j \le J_N(k)$, we obtain that

$$ \lim_{K \to \infty} \langle \tilde{\zeta}_T^{X|K}, \phi_N \rangle = \langle \tilde{\eta}_T^X, \phi_N \rangle \quad \text{almost surely.} \quad (23) $$

Notice that for any $x \in E$ such that $\|\tau(x)\|_{\infty} > N$, $\phi_N(x) = 0$. Hence for any $x = (n, s, \tau)$, one obtains the following inequality:

$$ |\phi_N(x) - \phi(x)| \le |\phi(x)| \mathbf{1}_{\{\|\tau\|_{\infty} > N\}} + \sum_{k=1}^{m(n,s)} \sum_{j \in J_N(k)} |\varphi_k^{n-s}(z_{k,j}^N) - \varphi_k^{n-s}(\tau_{1,n-s})| \mathbf{1}_{B_{k,j}^N}(\tau_{1,n-s}). \quad (24) $$

Notice that there exists at most one $(k, j) \in [1, m(n, s)] \times J_N(k)$ such that $\tau_{1,n-s} \in B_{k,j}^N$. As $\varphi_k^{n,s} \in C_b^1(\mathbb{R}^{n-s})$, the mean value inequality further implies that

$$ \forall k \le m(n, s), \forall j \in J_N(k), \forall z \in B_{k,j}^N, \quad |\varphi_k^{n-s}(z_{k,j}^N) - \varphi_k^{n-s}(z)| \le \|D\varphi_k^{n,s}\|_{\infty} d_N, $$
where $d_N$ denotes the maximum of the diameters of $d$-dimensional hypercubes of side length $2^{-N}$, for $d \le n_{\max}$. Letting $M = \max_{n,s,k} \|D\varphi_k^{n,s}\|_{\infty}$, it follows that:

$$ \forall x \in E, \quad |\phi_N(x) - \phi(x)| \le \|\phi\|_{\infty} \mathbf{1}_{\{\|\tau\|_{\infty} > N\}} + Md_N $$

Hence $\phi_N$ converges point-wise to $\phi$. Thus, by dominated convergence,

$$ \langle \tilde{\eta}_T^X, \phi_N \rangle \underset{N \to \infty}{\longrightarrow} \langle \tilde{\eta}_T^X, \phi \rangle. \qquad (25) $$

Furthermore, it follows from Equation (4.1.3) that for any $K \ge 1$,

$$ \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \phi_N \rangle - \langle \tilde{\zeta}_T^X|K, \phi \rangle|] \le \|\phi\|_{\infty} \sup_{K \ge 1} \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \mathbf{1}_{\{\|\tau(\cdot)\|_{\infty} > N\}} \rangle|] + Md_N. $$

Reasoning as in the proof of Proposition 4.7, and using Assumption 3.1, we obtain that

$$ \lim_{N \to \infty} \sup_{K \ge 1} \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \mathbf{1}_{\{\|\tau(\cdot)\|_{\infty} > N\}} \rangle|] = 0, $$

and as a consequence,

$$ \lim_{N \to \infty} \sup_{K \ge 1} \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \phi_N \rangle - \langle \tilde{\zeta}_T^X|K, \phi \rangle|] = 0. \qquad (26) $$

Noticing that

$$ \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \phi \rangle - \langle \tilde{\eta}_T^X, \phi \rangle|] \le \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \phi - \phi_N \rangle|] + \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \phi_N \rangle - \langle \tilde{\eta}_T^X, \phi_N \rangle|] + |\langle \tilde{\eta}_T^X, \phi_N - \phi \rangle| $$

together with Equations (26), (23) and (25) finally yields that

$$ \lim_{K \to \infty} \mathbb{E}[|\langle \tilde{\zeta}_T^X|K, \phi \rangle - \langle \tilde{\eta}_T^X, \phi \rangle|] = 0. $$

This concludes the proof. $\square$

### 4.1.4 Proof of Theorem 3.2

The previous results are sufficient to establish Theorem 3.2. Indeed, it follows from Propositions 4.7 and 4.8 that from every subsequence of $(\zeta^K)_{K \ge 1}$, one may extract a subsubsequence converging in $D(\mathbb{R}_+, (M_F(E), w))^2$ to a solution of Equation (6) which is continuous with respect to the total variation norm. As by assumption, $\zeta_0^K$ converges in law to $\eta_0 \in \mathfrak{M}_1$, Proposition 4.2 implies that all of these subsubsequences converge to the same limit $\eta$, which is the unique solution of Equation (6) with initial condition $\eta_0$. As $\eta_0 \in \mathfrak{M}_1$, Proposition 4.2 further ensures that $\eta \in D(\mathbb{R}_+, \mathfrak{M}_1)$. This establishes the convergence of $(\zeta^K)_{K \ge 1}$ to $\eta$ in $D(\mathbb{R}_+, M_1(E))^2$.

## 4.2 Proof of Theorem 3.3

This section is devoted to extracting dynamical system (8a-c) from the measure-valued integral equation (6), under the assumption that $\nu$ is the exponential distribution of parameter $\gamma$.

### 4.2.1 Preliminary study of the dynamical system

Before establishing Theorem 3.3 itself, let us start by showing that dynamical system (8a-c) endowed with initial condition (9) admits at most a unique solution. Existence will follow from the proofs of the forthcoming subsections, since they construct a solution to the dynamical system.
For this section, let us rewrite dynamical system (8a-c) 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 C^1(\mathbb{R}_+, \mathbb{R}^d)$ and $f : \mathbb{R}^d \to \mathbb{R}^d$ be defined such that dynamical system (8a-c) amounts to

$$y'(t) = f(y(t)) \quad \forall t \ge 0. \tag{27}$$

The components of $y$ (and resp. $f$) will be called $s$, $i$ and $n_{S,I}^X$ (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_{S,I}^X \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_{S,I}^X : 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}$

$$f_{X,S,I}(y) = -\left[ \left( \lambda_X I + \frac{\tau_X(y)}{s} + \tau_G(y) \right) S - \gamma I \right] n_{S,I}^X + \gamma(I+1)n_{S,I+1}^X \mathbf{1}_{\{S+I<n_{\max}\}} \\ + \left( \lambda_X(I-1) + \frac{\tau_X(y)}{s} + \tau_G(y) \right) (S+1)n_{S+1,I-1}^X \mathbf{1}_{\{I \ge 1\}}.$$

Also, notice that there are some natural constraints that we expect the solution of dynamical system (8a-c) to satisfy. Clearly, $s$, $i$ and $n_{S,I}^X$ should belong to $[0, 1]$. Also, as the population is partitioned into susceptible, infected and removed individuals, it follows that $s + i \le 1$. Similarly, as all individuals belong to exactly one household and one workplace, and as $n_{S,I}^X$ 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_{S,I}^X \le 1, \text{ and } \sum_{(S,I) \in \mathbb{S}} Sn_{S,I}^X \le m_X s. \tag{28}$$

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

$$V = \left\{ y \in [0, 1]^d : s + i \le 1, \sum_{(S,I) \in \mathbb{S}} n_{S,I}^X \le 1 \text{ and } m_X s - \sum_{(S,I) \in \mathbb{S}} Sn_{S,I}^X \ge 0 \quad \forall X \in \{H, W\} \right\}.$$

**Proposition 4.10.** Let $y^* \in V$. Then the following assertions hold:

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

(ii) For any $T \ge 0$, this problem admits at most a unique solution $y$ on $[0, T]$.

(iii) In particular, for any $\varepsilon > 0$, the dynamical system (8a-c) endowed with initial condition (9) admits at most a unique solution.

The proof of this proposition is available in Appendix C. It relies on establishing the Lipschitz continuity of $f$ on $V$, from which uniqueness is deduced using Gronwall's lemma.
### 4.2.2 Some properties of the limiting measure η

The results of this section focus on the limiting measure η given by Theorem 3.2, and will be useful for establishing Theorem 3.3.

Let us introduce the following notations. For $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$ and $T \ge t \ge 0$, define

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

We further define, for $X \in \{H, W\}$ and $t \ge 0$, the following quantity which relates to the infectious pressure exerted on susceptibles outside of their structure of type X:

$$\Lambda_X(t) = \frac{\lambda_{\bar{X}}}{s_{\bar{X}}(t)} \langle \eta_t^{\bar{X}}, \mathbf{si} \rangle + \beta_G \frac{i_H(t)}{n_H}.$$

We can now state a result which is similar in spirit to Proposition 4.4. Notice that it holds under the same Assumptions as Theorem 3.2, and is not restricted to the Markovian case.

**Proposition 4.11.** Let $\eta$ be the unique solution in $C(\mathbb{R}_+, \mathfrak{M}_1)$ of Equation (6). Then for any $T \ge 0$ and $t \in [0, T]$, for any measurable bounded function $f: \mathbb{R}_+ \times E \to \mathbb{R}$, it holds for $X \in \{H, W\}$ that

$$\langle \eta_t^X, f_{T,t} \rangle = \langle \eta_0^X, f_{T,0} \rangle + \lambda_X \int_0^t \langle \eta_u^X, \mathbf{si}(f_{T,u}^I - f_{T,u}) \rangle dt + \int_0^t \Lambda_X(u) \langle \eta_u^X, \mathbf{s}(f_{T,u}^I - f_{T,u}) \rangle du. \quad (30)$$

In particular, the application $t \mapsto \langle \eta_t^X, f_{T,t} \rangle$ is continuous.

*Proof.* Start by noticing that for any $f \in C_b^1(\mathbb{R}_+ \times E, \mathbb{R})$, for $X \in \{H, W\}$,

$$\langle \eta_t^X, f_t \rangle = \langle \eta_0^X, f_{t,0} \rangle + \lambda_X \int_0^t \langle \eta_u^X, \mathbf{si}(f_{t,u}^I - f_{t,u}) \rangle dt + \int_0^t \Lambda_X(u) \langle \eta_u^X, \mathbf{s}(f_{t,u}^I - f_{t,u}) \rangle du. \quad (31)$$

Indeed, the proof of Equation (31) follows the exact same lines as the proof of Proposition 4.4, showing that for any $f \in C_b^1(\mathbb{R}_+ \times E)$, Equation (31) leads to Equation (6) using Lemma 4.3.

Consider now a measurable bounded function $f: \mathbb{R}_+ \times E \to \mathbb{R}$. Proceeding as in the proof of Lemma 4.1, we consider a mollifier $\psi$ on $\mathbb{R}^{1+n_{max}}$ and let $\psi_k(t, \tau) = k^{n_{max}+1}\psi(kt, k\tau)$ for any $(t, \tau) \in \mathbb{R} \times \mathbb{R}^{n_{max}}$. Letting $f_k(t, (n, s, \tau)) = f(\cdot, (n, s, \cdot)) * \psi_k(t, \tau)$, we obtain by convolution a sequence of smooth functions $(f_k)_{k \ge 1}$ which converges point-wise to $f$.

Then on the one hand, for any $t \in [0, T]$ and $x \in E$, as $k$ tends to infinity, $(f_k)_{T,t}(x)$ converges to $f_{T,t}(x)$, and further $(f_k)_{T,t}^I(x)$ converges to $f_{T,t}^I(x)$ by dominated convergence. Define $g_k \in C_b^1(E, \mathbb{R})$ by $g_k(t, x) = (f_k)_{T,t}(x)$, it then holds that $(g_k)_{t,u} = (f_k)_{T,u}$ and $(g_k)_{t,u}^I = (f_k)_{T,u}^I$. Thus applying Equation (31) to $g_k$ and using dominated convergence as $k$ goes to infinity yields the desired result.

Finally, the continuity of $t \mapsto \langle \eta_t^X, f_{T,t} \rangle$ on $[0, T]$ is a consequence of Equation (30), as the integrands of the right-hand-side are bounded. $\square$

Proposition 4.11 allows us to establish the following result under the assumption that $\nu$ is the exponential distribution. In particular, it implies that within each structure, at any time, the remaining infectious periods of currently infectious individuals are independent and identically distributed, of common law $\nu$.

**Proposition 4.12.** Assume that $\nu$ is the exponential law of parameter $\gamma$, and consider $\eta$ as defined in Theorem 3.2 with initial condition $\eta_0 = \eta_{0,\varepsilon}$. Let $(X_n)_{n \ge 0}$ be a sequence of independent identically distributed random variables of common law $\nu$. Let $T \ge 0$, $n \in [1, n_{max}]$ and $s \in [0, n-1]$. For any $m \in [0, n-s]$, any functions $f \in B_b(\mathbb{R}_+ \times \mathbb{R}^m)$, $g_1, \dots, g_{n-s-m} \in B_b(\mathbb{R}_+ \times \mathbb{R})$
and any $j_1 < \dots < j_m$ and $k_1 < \dots < k_{n-s-m}$ such that $\{k_1, \dots, k_{n-s-m}\} \cup \{j_1, \dots, j_m\} = [1, n-s]$, define

$$F(t, x) = \mathbf{1}_{\left\{ \begin{subarray}{l} \mathbf{n}(x)=n, \mathbf{s}(x)=s, \\ \tau_{j_\ell}(x)>0 \quad \forall 1 \le \ell \le m \end{subarray} \right\}} \left( f(t, \tau_{j_1}(x), \dots, \tau_{j_m}(x)) - \mathbb{E}[f(t, X_1, \dots, X_m)] \right) \prod_{\ell=1}^{n-s-m} g_\ell(t, \tau_{k_\ell}(x)).$$

Then

$$\forall t \in [0, T], \quad \langle \eta_t^X, F_{T,t} \rangle = 0. \qquad (32)$$

*Proof.* Let $n \in [1, n_{\max}]$. Consider any $m, f, g, \{j_1, \dots, j_m\}$ and $\{k_1, \dots, k_m\}$ satisfying the constraints given in the proposition, and define $F$ as above. Throughout the proof, we let $E_m(t) = \mathbb{E}[f(t, X_1, \dots, X_m)]$.

Notice that $F$ satisfies the assumptions of Proposition 4.11. Letting $C = (\lambda_X n_{\max} + \lambda_{\bar{X}} n_{\max} + \beta_G)n_{\max}$, it follows from Equation (30) that:

$$|\langle \eta_t^X, F_{T,t} \rangle| \le |\langle \eta_0^X, F_{T,0} \rangle| + C \int_0^t \left( |\langle \eta_u^X, F_{T,u} \rangle| + |\langle \eta_u^X, F_{T,u}^I \rangle| \right) du. \quad (33)$$

Let $(Y_k)_{k \ge 1}$ be a sequence of independent identically distributed random variables of common law $\nu$ which is independent from $(X_k)_{k \ge 1}$, and write $p_{n,s}^X = \pi_n^X \binom{n}{s} (1-\varepsilon)^s \varepsilon^{n-s}$. Then on the one hand, Equation (7) ensures:

$$\langle \eta_0^X, F_{T,0} \rangle = p_{n,s}^X \left( \begin{aligned} & \mathbb{E}[f(T, Y_{j_1} - T, \dots, Y_{j_m} - T) \mathbf{1}_{\{Y_{j_1} - T > 0, \dots, Y_{j_m} - T > 0\}}] \\ & - \mathbb{E}[f(T, X_1, \dots, X_m)] \mathbb{E}[\mathbf{1}_{\{Y_{j_1} - T > 0, \dots, Y_{j_m} - T > 0\}}] \end{aligned} \right) \prod_{\ell=1}^{n-s-m} \mathbb{E}[g_\ell(T, Y_{k_\ell} - T)].$$

Usual properties of the exponential distribution thus lead to $\langle \eta_0^X, F_{T,0} \rangle = 0$.

On the other hand, let us compute $F_{T,u}^I(x)$ for any $x \in E$. Distinguishing the cases $n-s \in \mathbb{J}_m = \{j_1, \dots, j_m\}$ and $n-s \notin \mathbb{J}_m$ leads to:

$$F_{T,u}^I(x) = \mathbf{1}_{\left\{ \begin{subarray}{l} \mathbf{n}(x)=n, \mathbf{s}(x)=s+1, \\ \tau_j(x)>T-u \quad \forall j \in \mathbb{J}_m \setminus \{n-s\} \end{subarray} \right\}} \left( \mathbf{1}_{\{n-s \in \mathbb{J}_m\}} e^{-\gamma(T-u)} a_{T,u}(x) \right. \\ \left. + \mathbf{1}_{\{n-s \notin \mathbb{J}_m\}} \mathbb{E}[g_{n-s}(X_{n-s} - (T-u))] b_{T,u}(x) \right)$$

where

$$\begin{align*} a(t, x) &= \left( \mathbb{E}[f(t, \tau_{j_1}(x), \dots, \tau_{j_{m-1}}(x), X_{n-s})] - E_m(t) \right) \prod_{j \in [1, n-s] \setminus \mathbb{J}_m} g_j(t, \tau_{k_j}(x)), \\ b(t, x) &= \left( f(t, \tau_{j_1}(x), \dots, \tau_{j_m}(x)) - E_m(t) \right) \prod_{j \in [1, n-s-1] \setminus \mathbb{J}_m} g_j(t, \tau_{k_j}(x)). \end{align*}$$

We are now ready to proceed by induction on $s \in [0, n-1]$. First, consider the case $s=n-1$. Then either $m=0$ and the result is immediate, or $m=1$ and $\mathbb{J}_1 = \{1\}$. Hence necessarily $n-s \in \mathbb{J}_1$ and $a(x) = 0$ as $\mathbb{E}[f(X_0-(T-t))\mathbf{1}_{\{X_0>T-t\}}] = \mathbb{E}[f(X_1)]$. Thus $F_{T,t}^I(x) = 0$ for any $x \in E$, and Equation (33) reduces to

$$|\langle \eta_T^X, F \rangle| = |\langle \eta_T^X, F_{T,T} \rangle| \le C \int_0^T |\langle \eta_t^X, F_{T,t} \rangle| dt.$$

Recalling that $t \mapsto \langle \eta_t^X, F_{T,t} \rangle$ is continuous on $[0, T]$ according to Proposition 4.11, Gronwall's inequality ensures that $|\langle \eta_T^X, F_{T,t} \rangle| = 0$ for every $t \in [0, T]$.
Suppose now that $s < n - 1$, and that the result holds for $s + 1$. It then follows from the induction hypothesis that for any $u \in [0, T]$,

$$ \langle \eta_u^X, \mathbf{1}_{\left\{ \begin{subarray}{l} \mathbf{n}(x)=n, \mathbf{s}(x)=s+1, \\ \tau_j(x)>T-u \quad \forall j \in \mathbb{J}_m \setminus \{n-s\} \end{subarray} \right\}} a_{T,u}(\cdot) \rangle = \langle \eta_u^X, \mathbf{1}_{\left\{ \begin{subarray}{l} \mathbf{n}(x)=n, \mathbf{s}(x)=s+1, \\ \tau_j(x)>T-u \quad \forall j \in \mathbb{J}_m \setminus \{n-s\} \end{subarray} \right\}} b_{T,u}(\cdot) \rangle = 0. $$

Hence $\langle \eta_t^X, F_{T,u}^I \rangle = 0$ for any $u \in [0, T]$, and Gronwall's inequality allows to conclude as previously. $\square$

### 4.2.3 Proof of Theorem 3.3

From now on, we assume that $\eta_0 = \eta_{0,\varepsilon}$ as given by (7), and that $\nu$ is the exponential distribution with parameter $\gamma$. Throughout this section, let $h$ be the Heaviside step function, i.e. $h(z) = \mathbf{1}_{\{z>0\}}$ for any real number $z$.

Let us establish the following proposition, which serves as a starting point of the proof of Theorem 3.3.

**Proposition 4.13.** Under the assumptions of Theorem 3.3, it holds that for any $X \in \{H, W\}$, $(s_X(t))_{t \ge 0}$ and $(i_X(t))_{t \ge 0}$ satisfy

$$ \begin{aligned} \frac{d}{dt}s_X(t) &= -n_X \left( \frac{\lambda_X}{n_X} \langle \eta_t^X, \mathbf{si} \rangle + \frac{\lambda_{\bar{X}}}{n_{\bar{X}}} \langle \eta_t^{\bar{X}}, \mathbf{si} \rangle + \beta_G \frac{i_H(t)}{n_H} \frac{s_X(t)}{n_X} \right), \\ \frac{d}{dt}i_X(t) &= -\frac{d}{dt}s_X(t) + \gamma i_X(t). \end{aligned} \quad (34) $$

Further,

$$ s_X(0) = (1 - \varepsilon)n_X \quad \text{and} \quad i_X(0) = \varepsilon n^X. \quad (35) $$

*Proof.* Notice that $\mathbf{s} \in C_b^1(E, \mathbb{R})$ is such that $\mathcal{A}\mathbf{s} = 0$ and $\mathbf{s}^T(x) - \mathbf{s}(x) = -1$ for all $x \in E$. It thus follows immediately from Equation (6) that, for any $X \in \{H, W\}$,

$$ s_X(T) = s_X(0) - \int_0^T \left( \lambda_X \langle \eta_t^X, \mathbf{si} \rangle + \lambda_{\bar{X}} \frac{s_X(t)}{s_{\bar{X}}(t)} \langle \eta_t^{\bar{X}}, \mathbf{si} \rangle + \beta_G \frac{i_H(t)}{n_H} s_X(t) \right) dt. $$

Further, since $\eta \in C(\mathbb{R}_+, (\mathfrak{M}_1, \|\cdot\|_{TV}))$, it follows that for any $t \ge 0$, $(\zeta_t^K)_{K \ge 1}$ converges in law to $\eta_t$. As $\mathbf{s}$ is continuous and bounded on $E$, this implies that $\langle \zeta_t^{X|K}, \mathbf{s} \rangle \to \langle \eta_t^X, \mathbf{s} \rangle$ when $K$ tends to infinity. The analogous result holds for $\mathbf{n}$. Hence Lemma 2.2 ensures that, for any $t \ge 0$, $s_X(t)/s_{\bar{X}}(t) = n^X/n_{\bar{X}}$. In other words,

$$ s_X(T) = s_X(0) - n_X \int_0^T \left( \frac{\lambda_X}{n_X} \langle \eta_t^X, \mathbf{si} \rangle + \frac{\lambda_{\bar{X}}}{n_{\bar{X}}} \langle \eta_t^{\bar{X}}, \mathbf{si} \rangle + \beta_G \frac{i_H(t)}{n_H} \frac{s_X(t)}{n_X} \right) dt. $$

As $\eta \in C(\mathbb{R}_+, (\mathfrak{M}_1, \|\cdot\|_{TV}))$, 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 (34) comes from the fundamental theorem of calculus.

Recall that $\eta_0 = \eta_{0,\varepsilon}$ as defined in Equation (7). In particular, we now have $n_X = m_X$, hence

$$ s_X(0) = \langle \eta_{0,\varepsilon}^X, \mathbf{s} \rangle = \sum_{n=1}^{n_{\max}} \pi_n^X \sum_{s=0}^{n} s \binom{n}{s} (1-\varepsilon)^s \varepsilon^{n-s} = (1-\varepsilon) \sum_{n=1}^{n_{\max}} n \pi_n^X = (1-\varepsilon)n^X. \quad (36) $$

This yields the first part of Equation (35).
It remains to take an interest in $i_X(t)$. As $\mathbf{i}$ does not belong to $C_b^1(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$. Thus, we may apply Proposition 4.11 and obtain that

$$ \begin{aligned} \langle \eta_T^X, \mathbf{i} \rangle = \langle \eta_0^X, \mathbf{i}_{T,0} \rangle &+ \lambda_X \int_0^T \langle \eta_t^X, \mathbf{si}(\mathbf{i}_{T,t}^T - \mathbf{i}_{T,t}) \rangle dt \\ &+ \lambda_{\bar{X}} \int_0^T \frac{1}{s_{\bar{X}}(t)} \langle \eta_t^{\bar{X}}, \mathbf{si} \rangle \langle \eta_t^X, \mathbf{s}(\mathbf{i}_{T,t}^T - \mathbf{i}_{T,t}) \rangle dt + \beta_G \int_0^T \frac{i_H(t)}{n_H} \langle \eta_t^X, \mathbf{s}(\mathbf{i}_{T,t}^T - \mathbf{i}_{T,t}) \rangle dt. \end{aligned} \quad (37) $$

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

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

hence

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

Injecting this into Equation (37) and using as before that $s_X(t)/s_{\bar{X}}(t) = n_X/n_{\bar{X}}$ yields

$$ i_X(T) = \langle \eta_0^X, \mathbf{i}_{T,0} \rangle + n_X e^{-\gamma T} \int_0^T e^{\gamma t} \left( \frac{\lambda_X}{n_X} \langle \eta_t^X, \mathbf{si} \rangle + \frac{\lambda_{\bar{X}}}{n_{\bar{X}}} \langle \eta_t^{\bar{X}}, \mathbf{si} \rangle + \beta_G \frac{i_H(t)}{n_H} \frac{s_X(t)}{n_X} \right) dt. \quad (38) $$

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

$$ \langle \eta_0^X, \mathbf{i}_{T,0} \rangle = \langle \eta_{0,\varepsilon}^X, \mathbf{i}_{T,0} \rangle = e^{-\gamma T} \sum_{n=1}^{n_{\max}} \pi_n^X \sum_{s=0}^{n} \binom{n}{s} (1-\varepsilon)^s \varepsilon^{n-s} (n-s) = e^{-\gamma T} \varepsilon n^X. \quad (39) $$

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

$$ \frac{d}{dT} i_X(T) = n_X \left( \frac{\lambda_X}{n_X} \langle \eta_T^X, \mathbf{si} \rangle + \frac{\lambda_{\bar{X}}}{n_{\bar{X}}} \langle \eta_T^{\bar{X}}, \mathbf{si} \rangle + \beta_G \frac{i_H^H}{n_H} \frac{s_X(T)}{n_X} \right) - \gamma i_X(T) = - \frac{d}{dT} s_X(T) - \gamma i_X(T). $$

We thus have recovered the second line of Equation (34). Finally, the second half of Equation (35) is obtained by a computation analogous to (36). This concludes the proof. $\square$

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

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

For $t \ge 0$, let $n_{S,I}^X(t) = \langle \eta_t^X, f^{S,I} \rangle$, which defines a continuous function on $\mathbb{R}_+$ as $f^{S,I}$ is bounded and measurable and $\eta \in C(\mathbb{R}_+, (\mathfrak{M}_1, \| \cdot \|_{TV}))$. 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\} = \{x \in E : \mathbf{s}(x)\mathbf{i}(x) > 0, (\mathbf{s}(x), \mathbf{i}(x)) \in \mathcal{S}\}. $$

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

$$ \frac{d}{dt} s_X(t) = -n_X \left( \frac{\lambda_X}{n_X} \sum_{(S,I) \in \mathcal{S}} S \mathrm{In}_{S,I}^X(t) + \frac{\lambda_{\bar{X}}}{n_{\bar{X}}} \sum_{(S,I) \in \mathcal{S}} S \mathrm{In}_{(S,I)}^{\bar{X}}(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_{\bar{X}}}{s_{\bar{X}}(t)} \langle \eta_t^{\bar{X}}, \mathbf{si} \rangle + \beta_G \frac{i_H(t)}{n_H} = \frac{\lambda_{\bar{X}}}{s_{\bar{X}}(t)} \sum_{(S,I) \in \mathcal{S}} S \mathrm{In}_{(S,I)}^{\bar{X}}(t) + \beta_G \frac{i_H(t)}{n_H}, $$
which may also be written in terms of the notations of Equation (8a-c):

$$ \Lambda_X(t) = \left( \frac{s\bar{X}(t)}{n_X} \right)^{-1} \tau_{\bar{X}}(t) + \beta_G \frac{i_H(t)}{n_H}. $$

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

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

$$ \begin{aligned} \frac{d}{dt} n_{S,I}^X(t) ={}& \gamma \left( (I+1)n_{S,I+1}^X(t) \mathbf{1}_{\{S+I<n_{\max}\}} - In_{S,I}^X(t) \right) \\ &+ \lambda_X \left( (S+1)(I-1)n_{S+1,I-1}^X(t) \mathbf{1}_{\{I\ge 1\}} - SIn_{S,I}^X(t) \right) \\ &+ \Lambda_X(t) \left( (S+1)n_{S+1,I-1}^X(t) \mathbf{1}_{\{I\ge 1\}} - Sn_{S,I}^X(t) \right). \end{aligned} \quad (40) $$

Further

$$ n_{S,I}^X(0) = \binom{S+I}{I} \pi_{S+I}^X (1-\varepsilon)^S \varepsilon^I. \quad (41) $$

*Proof.* Let us start by establishing the initial condition of Equation (41). Let $(S, I) \in \mathbb{S}$. We will make use of another expression of $f^{S,I}$, which will actually be of use throughout the proof. For two integers $j \le 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 $x \in E$,

$$ f^{S,I}(x) = \sum_{n=S+I}^{n_{\max}} \mathbf{1}_{\{\mathbf{n}(x)=n\}} \sum_{\mathbf{v} \in \mathbb{B}(n-S,I)} H^{n,S,\mathbf{v}}(\tau(x)), $$

where, for $\mathbf{v} \in \mathbb{B}(n-S, I)$ and any $\tau \in \mathbb{R}^{n_{\max}}$,

$$ H^{n,S,\mathbf{v}}(\tau) = \prod_{j \in \mathbf{v}} h(\tau_j) \prod_{j \in \llbracket 1, n-S \rrbracket \setminus \mathbf{v}} (1-h(\tau_j)). $$

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$ 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.

Let $\varepsilon > 0$. Throughout the following, for $X \in \{H, W\}$, $n \in \llbracket 1, n_{\max} \rrbracket$ and $s \in \llbracket 0, n \rrbracket$, let $p_{n,s}^X = \pi_n^X(\binom{n}{s})(1-\varepsilon)^s \varepsilon^{n-s}$. Notice that by definition of $\eta_{0,\varepsilon}$,

$$ n_{S,I}^X(0) = \sum_{n=S+I}^{n_{\max}} p_{n,S}^X \sum_{\mathbf{v} \in \mathbb{B}(n-S,I)} \langle \nu, h \rangle^{\#\mathbf{v}} \langle \nu, 1-h \rangle^{(n-S)-\#\mathbf{v}}. $$

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$ and leads to Equation (41).

Next, let us apply Proposition 4.11 to $f^{S,I}$. A brief computation, based on distinguishing the cases where $n-s \in \mathbf{v}$ or $n-s \notin \mathbf{v}$ for any $\mathbf{v} \in \mathbb{B}(n-S, I)$, yields

$$ (f^{S,I})_{T,t}^{\mathcal{I}} = \mathbf{1}_{\{I \ge 1\}} e^{-\gamma(T-t)} f_{T,t}^{S+1,I-1} + \mathbf{1}_{\{S+1 < n_{\max}\}} (1 - e^{-\gamma(T-t)}) f_{T,t}^{S+1,I}. $$

As a consequence, Proposition 4.11 ensures that

$$ \begin{aligned} n_{S,I}^X(T) ={}& \langle \eta_0^X, f_{T,0}^{S,I} \rangle - \int_0^T \left( \lambda_X \langle \eta_T^X, \mathbf{s}\mathbf{i} f_{T,t}^{S,I} \rangle + \Lambda_X(t) \langle \eta_t^X, \mathbf{s}f_{T,t}^{S,I} \rangle \right) dt \\ &+ \mathbf{1}_{\{I \ge 1\}} \int_0^T e^{-\gamma(T-t)} \left( \lambda_X \langle \eta_T^X, \mathbf{s}\mathbf{i} f_{T,t}^{S+1,I-1} \rangle + \Lambda_X(t) \langle \eta_t^X, \mathbf{s}f_{T,t}^{S+1,I-1} \rangle \right) dt \\ &+ \mathbf{1}_{\{S+I < n_{\max}\}} \int_0^T (1 - e^{-\gamma(T-t)}) \left( \lambda_X \langle \eta_T^X, \mathbf{s}\mathbf{i} f_{T,t}^{S+1,I} \rangle + \Lambda_X(t) \langle \eta_t^X, \mathbf{s}f_{T,t}^{S+1,I} \rangle \right) dt. \end{aligned} \quad (42) $$
We will thus focus on differentiating with respect to T the different expressions composing the right-hand-side.

Let us start with the term $\langle \eta_0^X, f_{T,0}^{S,I} \rangle$. By definition of $f_{T,0}^{S,I}$ and $\eta_{0,\epsilon}^X$, it holds that

$$ \langle \eta_0^X, f_{T,0}^{S,I} \rangle = \sum_{n=S+I}^{n_{\max}} p_{n,S}^X \binom{I}{n-S} e^{-\gamma I T} (1-e^{-\gamma T})^{n-S-I}. $$

Hence, using the Equality $(n-S-I)\binom{I}{n-S} = (I+1)\binom{n-S}{I}$, we obtain that

$$ \frac{d}{dT} \langle \eta_0^X, f_{T,0}^{S,I} \rangle = -\gamma I \langle \eta_0^X, f_{T,0}^{S,I} \rangle + \gamma (I+1) \sum_{n=S+I+1}^{n_{\max}} p_{n,S}^X \binom{n-S}{I} e^{-\gamma(I+1)T} (1-e^{-\gamma T})^{n-S-I-1} $$

and we thus recognise that

$$ \frac{d}{dT} \langle \eta_0^X, f_{T,0}^{S,I} \rangle = -\gamma I \langle \eta_0^X, f_{T,0}^{S,I} \rangle + \gamma (I+1) \langle \eta_0^X, f_{T,0}^{S,I+1} \rangle. \quad (43) $$

Let us now focus on the remaining terms of the right-hand side of Equation (42). This motivates a closer study of $f_{T,t}^{S,I}$ for any $(S, I) \in \mathbb{S}$. By definition, for $x \in E$,

$$ f_{T,t}^{S,I}(x) = \sum_{n=S+I}^{n_{\max}} \mathbf{1}_{\{\mathbf{n}(x)=n\}} \sum_{\mathbf{s}(x)=S} \prod_{\mathbf{v} \in \mathbb{B}(n-S,I)} h(\tau_j(x) - (T-t)) \prod_{\mathbf{v} \in [1,n-S]\setminus\mathbf{v}} (1-h(\tau_j(x) - (T-t))). $$

In particular, $f_{T,t}^{S,I}(x) > 0$ requires that at time $t$, in a structure of type $x$, there are $J \ge I$ infected individuals, out of which exactly $I$ must have a remaining infectious period exceeding $(T-t)$. Hence, consider $(S, I)$ and $T$ to be fixed and define, for any $J \in [I, n_{\max} - S]$ and $x \in E$,

$$ g_J(t, x) = \sum_{n=S+J}^{n_{\max}} \mathbf{1}_{\{\mathbf{n}(x)=n\}} \sum_{\mathbf{s}(x)=S} \sum_{\mathbf{v}_0 \in \mathbb{B}(n-S,J)} \sum_{\substack{\mathbf{v} \subseteq \mathbf{v}_0 \\ \#\mathbf{v}=I}} G^{n, \mathbf{v}_0, \mathbf{v}}(t, \tau(x)), $$

where, for $t \ge 0$ and $\tau \in \mathbb{R}^{n_{\max}}$,

$$ G^{n, \mathbf{v}_0, \mathbf{v}}(t, \tau) = \prod_{j \in \mathbf{v}} h(\tau_j - (T-t)) \prod_{j \in \mathbf{v}_0 \setminus \mathbf{v}} \mathbf{1}_{\{0 < \tau_j \le T-t\}} \prod_{j \in [1, n-S] \setminus \mathbf{v}_0} (1-h(\tau_j)). $$

Dependence of $g_J$ and $G^{n, \mathbf{v}_0, \mathbf{v}}$ on $(S, I)$ and $T$ is omitted in these notations for readability. It then holds that

$$ \forall x \in E, \quad f_{T,t}^{S,I}(x) = \sum_{J=I}^{n_{\max}-S} g_J(t, x). $$

Let $J \in [I, n_{\max} - S]$, $\mathbf{v}_0 = \{j_1, \dots, j_J\} \in \mathbb{B}(n-s, J)$ and $\mathbf{v} \subseteq \mathbf{v}_0$ such that $\#\mathbf{v} = I$. Proposition 4.12 applied to the functions $f: \mathbb{R} \times \mathbb{R}^J \to \mathbb{R}$ and $g_k: \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ defined by

$$ f(t, \tau_{j_1}, \dots, \tau_{j_J}) = \prod_{j \in \mathbf{v}} h(\tau_j - (T-t)) \prod_{j \in \mathbf{v}_0 \setminus \mathbf{v}} \mathbf{1}_{\{0 \le \tau_j \le T-t\}}, \\ g_k(t, \tau) = 1 - h(\tau) \quad \forall k \in [1, n-S-J], $$

leads to the following equality, for all $t \in [0, T]$:

$$ \begin{aligned} \langle \eta_t^X, g_J(t, \cdot) \rangle &= \sum_{n=S+J}^{n_{\max}} \sum_{\substack{\mathbf{v}_0 \in \mathbb{B}(n-S,J) \\ \mathbf{v} \subseteq \mathbf{v}_0 : \#\mathbf{v}=I}} e^{-\gamma(T-t)I} (1-e^{-\gamma(T-t)})^{J-I} \langle \eta_t^X, \mathbf{1}_{\{\mathbf{n}(\cdot)=n\}} H^{n,S, \mathbf{v}_0}(\tau(\cdot)) \rangle \\ &= \binom{J}{I} e^{-\gamma(T-t)I} (1-e^{-\gamma(T-t)})^{J-I} n_{S,J}^X(t). \end{aligned} $$
As a consequence, we obtain in particular that

$$ \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I} \rangle = \sum_{J=I}^{n_{\max}-S} \binom{J}{I} e^{-\gamma(T-t)I} (1 - e^{-\gamma(T-t)})^{J-I} S J n_{S,J}^X(t). $$

Differentiating with respect to $T$ yields

$$ \partial_T \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I} \rangle = -\gamma I \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I} \rangle + \mathbf{1}_{\{S+I < n_{\max}\}} \gamma(I+1) \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I+1} \rangle. $$

In particular, $(T, t) \mapsto \partial_T \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I} \rangle$ is continuous, as for any $(S, J) \in \mathbb{S}$, $n_{S,J}^X$ is continuous thanks to the continuity of $t \mapsto \eta_t^X$ with regard to the total variation norm. Let $g: \mathbb{R} \to \mathbb{R}_+$ be such that $t \mapsto g(t) \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I} \rangle$ is continuous on $\mathbb{R}_+$ for any $(S, I) \in \mathbb{S}$. It then holds that

$$ \frac{d}{dT} \int_0^T g(t) \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I} \rangle dt = g(T) S I n_{S,I}^X(T) - \gamma I \int_0^T g(t) \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I} \rangle dt \tag{44} \\ + \mathbf{1}_{\{S+I < n_{\max}\}} \gamma(I+1) \int_0^T g(t) \langle \eta_t^X, \mathbf{s} \mathbf{i} f_{T,t}^{S,I+1} \rangle dt. $$

Similarly, for $g: \mathbb{R} \to \mathbb{R}_+$ such that $t \mapsto g(t) \langle \eta_t^X, \mathbf{s} f_{T,t}^{S,I} \rangle$ is continuous on $\mathbb{R}_+$ for any $(S, I) \in \mathbb{S}$,

$$ \frac{d}{dT} \int_0^T g(t) \langle \eta_t^X, \mathbf{s} f_{T,t}^{S,I} \rangle dt = g(T) S n_{S,I}^X(T) - \gamma I \int_0^T g(t) \langle \eta_t^X, \mathbf{s} f_{T,t}^{S,I} \rangle dt \tag{45} \\ + \mathbf{1}_{\{S+I < n_{\max}\}} \gamma(I+1) \int_0^T g(t) \langle \eta_t^X, \mathbf{s} f_{T,t}^{S,I+1} \rangle dt. $$

In particular, we may apply Equations (44) and (45) to $g(t) = 1$, as well as $g(t) = \Lambda_X(t)$ and $g(t) = e^{-\gamma t} \Lambda_X(t)$. Indeed, $t \mapsto \Lambda_X(t) S n_{S,I}^X(t)$ is continuous and well defined for any $t \ge 0$, thanks to the inequality $\langle \eta_t^X, \mathbf{s} \mathbf{i} \rangle S n_{S,I}^X(t) \le (n_{\max} n_X / n_X) s_X(t)^2$ which ensures that despite the division by $s_X$ in the definition of $\Lambda_X$, there are no singularities.

In conclusion, Equation (43), together with Equations (44) and (45), allows to differentiate the right-hand-side of Equation (42). In particular, regrouping the terms factorised by $-\gamma I$ and $\gamma(I+1)$ allows to distinguish $-\gamma I n_{S,I}^X(T)$ and $\gamma(I+1) n_{S,I+1}^X(T)$, using Equation (42). This computation finally leads to Equation (40). $\square$

We may finally focus on the main result of this section, namely Theorem 3.3.

*Proof of Theorem 3.3.* Before concluding, we need to emphasize that it would have been possible to choose $X = W$ when replacing $S$ and $I$ by $K_X S_X$ and $K_X I_X$ for $\mathcal{I}_G$ in Proposition 4.5. All of the subsequent results still hold, simply replacing the household-related quantities in the definition of the rate for mean-field infections by their workplace-related counterparts.

As a consequence, Propositions 4.13 and 4.14 show that for any $X \in \{H, W\}$,

$$ y_X = \left( \frac{s_X}{m_X}, \frac{i_X}{m_X}, n_{S,I}^X : (S, I) \in \mathbb{S}, n_{(S,I)}^X : (S, I) \in \mathbb{S} \right) $$

satisfies the Cauchy problem (27) with initial condition (9). However, Proposition 4.10 ensures uniqueness of the solutions to this Cauchy problem. It hence is sensible to define, for $t \ge 0$,

$$ s(t) = \frac{s_H(t)}{m_H} = \frac{s_W(t)}{m_W} \quad \text{and} \quad i(t) = \frac{i_H(t)}{m_H} = \frac{i_W(t)}{m_W}. $$

This leads to dynamical system (8a-c) with initial conditions (9), and concludes the proof. $\square$
<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of the stochastic model starting from a single infected, with its large population approximation for two different choices of initial condition. The first initial condition is given by Equation (9) for $\varepsilon = 0.01$. The second initial condition is obtained by simulating a large number (> 2000) of stochastic epidemic trajectories, starting from a single infected until the proportion of infected reaches one percent. The initial condition corresponds to the average of the initial conditions observed in each simulation. Regarding the stochastic model, for this figure, 100 epidemics starting from a single infected were simulated. Similarly to Figure 2, only those reaching a threshold of 3% of infected are represented, and a time shift is applied to ease comparison between model outputs. Structure size distributions are those of Figure 1. Epidemic parameters: ($\beta_G, \lambda_H, \lambda_W, \gamma$) = (0.085, 0.1, 0.001, 0.125), $R_0$ = 1.7.</figcaption></figure>

Discussion

This paper has focused on proposing a new reduction for an SIR model with two levels of mixing, which explicitly includes households and workplaces. This reduced model was obtained as its large population limit, and the associated convergence of the stochastic model was established.

A possible model extension would be to consider a local level of mixing containing an arbitrary, yet finite, number of layers. As long as within each layer, each node is part of exactly one clique, and as long as cliques within each layer are constituted independently from one another as in the case for households and workplaces, the adaptation of the aforementioned results is expected to be straightforward. In particular, for exponentially distributed infectious period lengths, the dimension of the corresponding dynamical system should still be of order $O(n_{max}^2)$, implying that the model should remain tractable.

Furthermore, we have compared the reduced model obtained in this work with the corresponding EBCM in the line of [42]. In the case of our household-workplace model with two levels of mixing, the EBCM seems the less appropriate choice, as it is less parsimonious and only approaches the epidemic well if the initial proportion of infected is very small. However, this may change if a more general contact structure within layers is considered, such as a configuration model for the global level, in which case it seems sensible to assume that EBCM-like equations will appear.

Finally, let us emphasize that by essence, the large population limit obtained here corresponds to a situation where the number of infected individuals is of the same order as the population size. In a realistic scenario, however, an epidemic is initiated by very few infected individuals. In the case of a large epidemic outbreak, the number of infected subsequently grows until it no longer is negligible when compared to the population size, at which point the large population limit correctly captures the dynamics of the outbreak. This raises the question: which initial
condition is pertinent for the large population approximation? For uniformly mixing population, this is rather straightforward, for two main reasons. On the one hand, at each time, infected individuals are interchangeable in terms of infectious pressure exerted on susceptibles. On the other hand, the presence of recovered individuals at time t can be neglected in the study of the epidemic dynamics over the time interval $[t, \infty)$, simply by restricting the study to all other individuals, which still constitute a uniformly mixing population. As a consequence, in such a setting, it makes sense to suppose that at time zero, there are only infected and susceptible individuals, and that infected individuals are chosen uniformly at random in the population, with independent and identically distributed infectious period lengths.

This idea can of course be extended to our setting, and corresponds to the initial condition proposed in Theorem 3.3, while similar initial conditions have also been used in the literature in related settings [42, 11]. In our model however, one actually needs to know how infected and recovered individuals are distributed among households and workplaces, meaning that neither can recovered be ignored, nor is there any reason to believe that infected individuals are distributed uniformly at random in the population. Indeed, Figure 3 illustrates that when compared to stochastic simulations starting from a single infected, the large population approximation with initial condition given by Equation (9) fails to reproduce the epidemic dynamics, while they are correctly captured when using an initial condition which is inferred from stochastic simulations. As a consequence, it seems of interest to get a better understanding of this realistic initial condition, which arises from an epidemic started by a single infected. This may be achieved using a branching process approximation of the epidemic, which is designed to approach the initial, stochastic phase of the epidemic, and hence would represent a reduced model which complements the large population limit obtained in the present work.

# A Implementation of the large population limit

## A.1 Automatic implementation of the dynamical system

It is possible to implement dynamical system (8a-c) in an automated way, in the sense that equations do not need to be written individually. The key lies in the fact that the set $\mathbb{S}$ can be constructed automatically, with an intrinsic organization of the states $(S, I)$ it contains. For example, one may arrange them by growing number $n$ of susceptible and infected members of the structure, and for each $n$, by growing number $i$ of infected, leading to

$$ \mathbb{S} = \{(2,0), (1,1), \dots, (n_{\max}, 0), (n_{\max}-1, 1), \dots, (1, n_{\max}-1)\}. $$

This in turn allows to make an explicit correspondence between any state $(S, I) \in \mathbb{S}$ and e.g. some position in a vector containing all functions of our dynamical system of interest. A similar idea was already employed in [34], for another purpose. With the previous structure of $\mathbb{S}$, one may for instance notice that for any $n \in [2, n_{\max}]$ and $i \in [0, n-1]$, the state $(n - i, i)$ is the $c(n - i, i)$-th state enumerated in $\mathbb{S}$, where $c(n - i, i) = (n - 1)n/2 + i$. As a consequence, the general expression of Equation (8c) may be used to handle all the dynamics of the functions $n_{S,I}^X$, for $(S, I) \in \mathbb{S}$ and $X \in \{H, W\}$.

Also, notice that in practice, household sizes tend not to be as big as workplace sizes. It thus makes sense to distinguish explicitly a maximal size for each type of structure. This allows to avoid implementing unnecessary equations corresponding e.g. to household sizes that are not actually observed, and which thus artificially increase the dimension of the system.

## A.2 Computational performance

The aim of this section is to numerically assess the computational cost associated to solving the large dimensional dynamic system (8a-c) in comparison to stochastic simulations using
<table><caption>Table 1: Considered values of the contact rates and final times, grouped by value of R₀ and proportions of infections per layer (pG, pH, pW) characterizing the scenarios.</caption><thead><tr><th></th><th>R<sub>0</sub> = 1.2</th><th>R<sub>0</sub> = 1.4</th><th>R<sub>0</sub> = 1.7</th><th>R<sub>0</sub> = 2.0</th><th>R<sub>0</sub> = 2.5</th><th>(p<sub>G</sub>, p<sub>H</sub>, p<sub>W</sub>)</th></tr></thead><tbody><tr><td>β<sub>G</sub></td><td>0.03</td><td>0.035</td><td>0.045</td><td>0.05</td><td>0.06</td><td></td></tr><tr><td>λ<sub>H</sub></td><td>0.05</td><td>0.07</td><td>0.09</td><td>0.15</td><td>0.2</td><td>(0.2, 0.4, 0.4)</td></tr><tr><td>λ<sub>W</sub></td><td>0.0015</td><td>0.0016</td><td>0.0018</td><td>0.002</td><td>0.0022</td><td></td></tr><tr><td>T</td><td>130</td><td>130</td><td>105</td><td>85</td><td>75</td><td></td></tr><tr><td>β<sub>G</sub></td><td>0.06</td><td>0.07</td><td>0.085</td><td>0.1</td><td>0.125</td><td></td></tr><tr><td>λ<sub>H</sub></td><td>0.06</td><td>0.07</td><td>0.1</td><td>0.15</td><td>1.5</td><td>(0.4, 0.4, 0.2)</td></tr><tr><td>λ<sub>W</sub></td><td>0.00075</td><td>0.0008</td><td>0.001</td><td>0.0011</td><td>0.00115</td><td></td></tr><tr><td>T</td><td>145</td><td>130</td><td>95</td><td>80</td><td>55</td><td></td></tr></tbody></table>

Gillespie's algorithm, also referred to as SSA (stochastic simulation algorithm). In order to do so, the average execution times of one stochastic simulation (SSA) and of one resolution of the associated dynamic system using the ODE solver odeint from the scipy.integrate library are compared.

Let us start by describing the general procedure. Each of the two scripts (stochastic simulation or reduced model) is executed one hundred times, all runs being independent from one another. For each run and each script, the computation time of the script of interest is measured, as well as the computation time of a reference function (summing all integers up to one billion with a simple for-loop). The ratios of the runtimes of both the script of interest and the reference function are computed. Comparison of the computation times for the stochastic and the reduced model is then based on the comparison of the averages of those normalised runtimes.

It remains to take an interest in the choice of the model parameters, namely the structure size distributions, the epidemic parameters i.e. the contact rates $β_G$, $λ_H$, $λ_W$ and the removal rate $γ$, as well as the initial proportion of infected $ε$ and the time interval $[0, T]$ on which the epidemic is simulated. For the stochastic model, the population size $K$ will be fixed to ten thousand individuals. For all scenarios considered here, the structure size distributions will be those of Figure 1, and the initial proportion of individuals will be set to $ε = 0.005$. Different values of the epidemic parameters will be considered, as to obtain scenarios that differ both in terms of $R_0$ and in terms of the proportions of infections occurring within the general population, within households or within workplaces, respectively referred to as $p_G$, $p_H$ and $p_W$. The removal rate $γ$ will be fixed at 0.125, and only the contact rates will effectively vary. In total, ten different scenarios will be used, characterized by their values of $R_0 \in \{1.2, 1.4, 1.7, 2.0, 2.5\}$ and $(p_G, p_H, p_W) \in \{(0.2, 0.4, 0.4), (0.4, 0.4, 0.2)\}$.

Finally, parameter $T$ will be chosen as follows. For each set of epidemic parameters detailed above, the reduced model is used to compute the time $T_*$ at which the epidemic falls below one percent of infected individuals in the population, after the epidemic peak. $T$ then is determined by rounding down $T_*$ to the closest multiple of five.

Figure 4 uses the reduced model to plot the trajectories of the proportion of susceptible and infected individuals in the population, for each scenario. The corresponding parameters are summarized in Table 1. Notice that in particular, this includes the parameters of Figure 2.

Let us now turn to the results. For each scenario of Table 1, measurement of average normalised computation times was repeated three times. The results are shown in Figure 5, which indicates for each scenario the ratio of the average normalised runtime for one resolution of dynamical system (8a-c) over the average normalised runtime of one stochastic simulation. Let us first take an interest in the datasets labeled $(p_G, p_H, p_W) = (0.2, 0.4, 0.4)$ and $(p_G, p_H, p_W) =$
<figure><img src="image_6.png" /><figcaption>Figure 4: Proportion of susceptible (S) and infected (I) in the population, for each scenario detailed in Table 1, as given by dynamical system (8a-c). Scenarios are separated by values of ($p_G, p_H, p_W$) of infections per layer, namely (0.2, 0.4, 0.4) and (0.4, 0.4, 0.2) for the top and bottom panels, respectively. The corresponding values of $R_0$ are indicated by the color shades, as shown in the legend. The black crosses indicate for each curve that the proportion of infected falls below the threshold of one percent.</figcaption></figure>

(0.4, 0.4, 0.2). One may notice first that for each scenario, the results of all three repeats are close to one another, indicating that the results are reproducible. Further, for both possible values of ($p_G, p_H, p_W$), the results indicate a shared general trend. Indeed, for values of $R_0$ close to the critical case $R_0 = 1$, the ratio exceeds one, and diminishes subsequently, falling below one between $R_0 = 1.4$ and $R_0 = 1.7$ and attaining values of order $10^{-1}$. This behavior suggests that solving dynamic system (8a-c) is advantageous in terms of computation time for intermediate or high values of $R_0$, being up to one order of magnitude faster than one stochastic simulation. As the time interval $[0, T]$ on which the epidemic is studied originally depends on the scenario and is significantly shorter for larger values of $R_0$, one may wonder whether this difference influences the results. As a consequence, we have repeated the same procedure for all of the scenarios characterised by ($p_G, p_H, p_W$) = (0.2, 0.4, 0.4), with fixed $T = 75$. Figure 5 shows that the associated results are very similar to those obtained previously, pleading against this hypothesis.

Of course, this comparison could be pushed further. For instance, the most basic version of the SSA algorithm was used, and more advanced methods such as τ-leaping are expected to accelerate stochastic simulations. Also, a more thorough exploration of the parameter space would be pertinent, assessing for instance the influence of the structure size distributions.

## B Edge-based compartmental model
<figure><img src="image_7.png" /><figcaption>Figure 5: Ratio of the average normalised computation time for solving once dynamical system (8a-c) over the average normalised computation time for one stochastic simulation (SSA). This ratio was computed three times for each scenario of Table 1. The results are presented as a function of $R_0$, while colors indicate the value of $(p_G, p_H, p_W)$. Unless stated otherwise, the parameter $T$ from Table 1 was used. The dotted line indicates the threshold of one.</figcaption></figure>

## B.1 Presentation of the EBCM

Let us start by describing how to obtain the population structure of the local level of mixture described in Section 2.1 using a clique configuration model (CCM). In our case, each node belongs to exactly one clique within each layer (one household and one workplace, respectively). Let us briefly notice that whenever a node is picked uniformly at random, the probability of it belonging to a structure of type $X$ and size $n$ is given by $\hat{\pi}_n^X = n\pi_n^X/m_X$, for any $n \in [1, n_{max}]$ and $X \in \{H, W\}$. As a consequence, the layer corresponding to structures of type $X \in \{H, W\}$ is obtained by the following two steps. First, associate to each node a structure size distributed according to the size-biased law $\hat{\pi}^X$. This is done independently for each node. Second, for $k \in [1, n_{max}]$, form cliques of size $k$ by drawing uniformly without replacement $k$-tuples in the set of nodes of associated structure size $k$. This step stops when all nodes of associated clique size $k$ belong to a clique. This procedure is repeated independently for each layer, allowing to assemble households and workplaces.

Let us now turn to deriving the EBCM. Consider $s$ and $i$ the proportions of susceptible and infected individuals in the population, respectively. Let $\theta_n^X(t)$ for $X \in \{H, W\}$ and $n \in [1, n_{max}]$ be the chance of a susceptible belonging to a structure of type $X$ and size $n$ to escape infection within this structure, and $\theta^G(t)$ the chance of escaping infection through the mean-field level, up to time $t \ge 0$. The key idea is that a node is susceptible at time $t$ if and only if it has escaped infection up to time $t$, and the risks of infection within each layer are independent from one another. This makes use of properties of the CCM, which heuristically explain the decoupling of the risk of infection in the two local layers from one another. Further, the fact that in an infinite population, each individual structure has a negligible impact on the proportion of infected yields the intuition behind the decorrelation of the risks of infection at the local and global level. This
<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison of the stochastic model with the large population approximation given by dynamical system (8a-c) and the corresponding EBCM. Household and workplace distributions are those of Figure 1. Epidemic parameters are set to ($\beta_G, \lambda_H, \lambda_W, \gamma$) = (0.125, 1.5, 0.00115, 0.125). Initial conditions correspond to $\varepsilon \in \{0.001, 0.01, 0.05\}$ as indicated for each panel. For each of these scenarios, Gillespie's algorithm is used to simulate 50 trajectories of the stochastic model defined in Proposition 2.1 in a population of $K = 10000$ individuals (faint lines). For Panel (a), only trajectories reaching a threshold proportion of 0.005 infected are kept, and time is shifted so that time 0 corresponds to the moment when this threshold is reached. Finally, the deterministic solution $(s, i)$ of both dynamical system (8a-c) (thick lines) and the EBCM (dashed lines) are represented for each scenario. For Panel (a), the same time shifting procedure as for simulations is applied.</figcaption></figure>
leads to

$$s = \theta^G \prod_{X \in \{H, W\}} \left( \sum_{n=1}^{n_{\max}} \hat{\pi}_n^X \theta_n^X \right).$$

As we are considering an SIR model, it follows that $i'(t) = -s'(t) + \gamma i(t)$, so that the difficulty resides in understanding the dynamics of $\theta^G$ and $\theta_n^X(t)$, for $X \in \{H, W\}$ and $n \in [1, n_{\max}]$.

Define for $X \in \{H, W\}$ and $n \in [2, n_{\max}]$:

$$m_n^X = \theta^G \hat{\pi}_n^X \theta_n^X \sum_{k=1}^{n_{\max}} \hat{\pi}_k^X \theta_k^X,$$

which corresponds to the proportion of individuals who are susceptible and belong to a structure of type X and size n. Also, let $n_{(S,I,R)}^X$ be the proportion of susceptibles belonging to a structure of type X containing exactly S susceptibles, I infected and R removed individuals. This allows us to introduce the following quantities, which participate in the rates at which a member of a structure of type X and size n is infected, either within the considered structure or outside of it, respectively:

$$T_n^X = \lambda_X \sum_{\substack{(S,I,R) \in \mathbb{N}^3 \\ S+I+R=n}} SI n_{(S,I,R)}^X \quad \text{and} \quad \tau_n^X = \left( \beta_G i + \sum_{k=1}^{n_{\max}} T_k^{\overline{X}} \right) \frac{\hat{\pi}_n^X \theta_n^X}{\sum_{k=1}^{n_{\max}} \hat{\pi}_k^X \theta_k^X}.$$

One obtains the following dynamics:

$$\frac{d}{dt} \theta^G = -\beta_G i \theta^G, \quad \text{and} \quad \forall X \in \{H, W\}, \forall n \in [2, n_{\max}], \quad \frac{d}{dt} \theta_n^X = -\frac{T_n^X}{m_n^X} \theta_n^X.$$

Further, for any $X \in \{H, W\}$, $n \in [2, n_{\max}]$ and $(S, I, R) \in \mathbb{N}^3$ such that $S+I+R=n$ and either $S \ge 2$ or $SI \ge 1$:

$$\frac{d}{dt} n_{(S,I,R)}^X = - \left( \lambda_X SI + \frac{\tau_n^X}{m_n^X} S + \gamma I \right) n_{(S,I,R)}^X \\ + \gamma(I+1)n_{(S,I+1,R-1)}^X \mathbf{1}_{\{R \ge 1\}} \\ + \left( \lambda_X(S+1)(I-1) + \frac{\tau_n^X}{m_n^X}(S+1) \right) n_{(S+1,I-1,R)}^X \mathbf{1}_{\{I > 1\}}.$$

Additionally, as in a structure of size one, no infection may occur within the structure itself, $\theta_1^H$ and $\theta_1^W$ are constant over time. Finally, it remains to define the initial conditions. Following Volz et al. [42], we consider the case $\varepsilon \ll 1$. Then the only quantities which are not null at time zero are: for any $X \in \{H, W\}$, $n \in [1, n_{\max}]$ and $I \in [1, n-1]$,

$$\begin{align}
i(0) &= \varepsilon, \nonumber \\
\theta^G(0) &= \theta_n^X(0) = 1 - \varepsilon, \nonumber \\
n_{(n,0,0)}^X &= \frac{1}{n} \hat{\pi}_n^X (1 - \varepsilon)^n, \tag{46} \\
n_{(n-I,I,0)}^X &= \hat{\pi}_n^X (1 - \varepsilon)^{n-I} \varepsilon^I. \nonumber
\end{align}$$

The proportions of susceptible and infected as predicted by both the EBCM and dynamical system (8a-c) are shown in Figure 6, for different values of $\varepsilon$. Let us first notice that in the case $\varepsilon = 0.001$, corresponding to Panel (a) of Figure 6, the solutions $(s, i)$ of both the EBCM and dynamical system (8a-c) are in perfect accordance, emphasizing the fact that for very small
<table><caption>Table 2: Numerical assessment of the computation time needed to solve either dynamical system (8a-c) or the EBCM introduced in Appendix B. Model parameters: household and workplace size distribution from Figure 1; ($\beta_G$, $\lambda_H$, $\lambda_W$, $\gamma$) = (0.125, 1.5, 0.00115, 0.125); initial proportion of infected $\varepsilon = 0.005$; resolution of the numerical system over the time interval [0, 30].</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Runs</th><th colspan="3">Normalised runtimes</th></tr><tr><th>Average</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>Dynamical system (8a-c)</td><td>50</td><td>0.15</td><td>0.14</td><td>0.17</td></tr><tr><td>EBCM</td><td>10</td><td>2076</td><td>1887</td><td>2254</td></tr></tbody></table>

values of $\varepsilon$, the EBCM seems to yield the correct asymptotic population dynamics. However, for larger values of $\varepsilon$, the EBCM struggles to reproduce these dynamics. The problem for capturing the epidemic dynamics for higher values of $\varepsilon$ lies in the fact that defining the proper initial condition for the EBCM is not straightforward, leading to initial conditions consisting in an approximation which is only sensible whenever $\varepsilon$ is very small.

## B.2 Computational performance

In order to compare the computation times needed to solve either dynamical system (8a-c) or the dynamical system associated to the EBCM which has been introduced above, we will proceed similarly as in Appendix A.2, making use of the ODE solver odeint from the scipy.integrate library in both cases. However, this time, only one parameter set will be used, corresponding to the parameters chosen for Panel (a) of Figure 2. Further, the average normalised computation time is only computed once, instead of having three repeats as in Appendix A.2. Considering the relatively small fluctuations between repeats for all scenarios in Figure 5, this is not expected to significantly affect the qualitative result.

The model parameters and the associated average runtimes are shown in Table 2. Due to the excessive computation needed to solve the EBCM, only 10 runs of this script were performed. However, considering that the average normalised runtime for solving the EBCM is several orders of magnitude higher than the average normalised runtime for solving dynamical system (8a-c), this again is not expected to significantly alter the results. Finally, the computation times necessary for solving the EBCM are relatively homogeneous over all runs, indicating that the average computation time is not biased by an outlier.

## C Proof of Proposition 4.10

Let us start with the following lemma, which will be needed afterwards.

**Lemma C.1.** Consider a solution $y$ of dynamical system (8a-c) and let $\Delta(t) = m_X s(t) - \sum_{(S,I) \in \mathcal{S}} \mathrm{Sn}_{S,I}^X(t)$. Then

$$ \frac{d}{dt} \Delta(t) = \gamma n_{1,1}^X(t) - \left( \tau_G(t) + \frac{1}{s(t)} \tau_{\bar{X}}(t) \right) \Delta(t). $$

*Proof of Lemma C.1.* Let $X \in \{H, W\}$. First, notice that

$$ \{(S + 1, I - 1) : (S, I) \in \mathcal{S}, I \ge 1\} = \{(S, I) \in \mathcal{S} : S > 1\}. $$
As $S - 1 = 0$ whenever $S = 1$, we thus obtain that

$$ \begin{aligned} & \sum_{(S,I) \in \mathbb{S}} S^2 I n_{S,I}^X - \sum_{(S,I) \in \mathbb{S}} S(S+1)(I-1) n_{S+1,I-1}^X \mathbf{1}_{\{I \ge 1\}} \\ & = \sum_{(S,I) \in \mathbb{S}} S^2 I n_{S,I}^X - \sum_{(S,I) \in \mathbb{S}} (S-1) S I n_{S,I}^X = \sum_{(S,I) \in \mathbb{S}} S I n_{S,I}^X. \end{aligned} \quad (47) $$

Similarly, $\{(S, I + 1) : (S, I) \in \mathbb{S}, S + I < n_{\max}\} = \mathbb{S} \setminus \{(S, I) \in \mathbb{S} : I = 0\} \cup \{(1, 1)\}$. As $SI = 0$ whenever $I = 0$, it follows that

$$ \sum_{(S,I) \in \mathbb{S}} \gamma S I n_{S,I}^X - \sum_{(S,I) \in \mathbb{S}} \gamma S (I+1) n_{S,I+1}^X \mathbf{1}_{\{S+I < n_{\max}\}} = \gamma n_{1,1}^X. \quad (48) $$

The desired conclusion then results directly from Equations (8a-c), regrouping the terms of the form of Equations (47) and (48) in order to simplify the expression. $\square$

We are now ready to focus on the desired result.

*Proof of Proposition 4.10.* (i) By assumption, $y(0) \in V$. Let us start by checking that all components of $y$, as well as $\Delta$, stay non-negative over time.

Let $t_0 \ge 0$ be such that $y(t_0) \in V$. If $i(t_0) = 0$, then $i'(t_0) = -s'(t_0) \ge 0$ by assumption, which ensures that $i$ will not become negative on a neighbourhood of $t_0$. Similar arguments hold for the lower bounds of $\Delta$ and $s$, using Lemma C.1 and Inequality (28), respectively.

Let us now turn our attention to $n_{S,I}^X$ for $X \in \{H, W\}$ and $(S, I) \in \mathbb{S}$. Recall from Equation (8c) that its derivative may be ill defined, due to the division by $s(t)$. However, inequality (28) ensures that the ratios $\tau_{\bar{X}} n_{S,I}^X/s$ are well defined at all time, for any $(S, I) \in \mathbb{S}$. As a consequence, we may now notice as previously that, if $n_{S,I}^X(t_0) = 0$, Equation (8c) ensures that $\frac{d}{dt} n_{S,I}^X(t_0) \ge 0$.

The desired conclusion follows: whenever either of the quantities of interest reach zero, their derivatives are non-negative which ensures that they do not become negative shortly thereafter.

Next, let us have a look at the upper bounds. For $X \in \{H, W\}$, a brief computation yields

$$ \frac{d}{dt} \sum_{(S,I) \in \mathbb{S}} n_{S,I}^X(t) = -\gamma n_{1,1}^X(t) \le 0. $$

This assures that starting from $y^* \in V$, the inequality $\sum_{(S,I) \in \mathbb{S}} n_{S,I}^X(t) \le 1$ holds. For $X \in \{H, W\}$ and $(S, I) \in \mathbb{S}$, it follows that if $n_{S,I}^X(t_0) = 1$, then for any $(S', I') \in \mathbb{S} \setminus (S, I)$, $n_{S',I'}^X(t_0) = 0$. Thus

$$ \frac{d}{dt} n_{S,I}^X(t_0) = - \left( \lambda_X S I + \tau_{\bar{X}}(t_0) \frac{S}{s(t_0)} + \tau_G(t_0) S + \gamma I \right) \le 0, $$

from which one may deduce that $n_{S,I}^X$ remains less than or equal to one. The remaining upper bounds on $s$, $i$ and $s + i$ may be obtained using similar arguments.

We thus have established that if $y^* \in V$, then $y(t) \in V$ for all $t$ such that $y$ is well defined.

(ii) Consider any $T \ge 0$. Let $y = (s, i, n_{S,I}^X : X \in \{H, W\}, (S, I) \in \mathbb{S})$ be a solution to the Cauchy problem (27) with initial condition $y^* \in V$. Then it follows from Inequality (28) that

$$ s'(t) \ge -[(\lambda_H + \lambda_W)n_{\max} + \beta_G]s(t), $$

and as further $s \in C^1(\mathbb{R}_+)$ according to the Cauchy problem, we obtain by comparison that $s(t) \ge s(0) \exp(-[(\lambda_H + \lambda_W)n_{\max} + \beta_G]t)$ for any $t \in [0, T]$. As a consequence, on $[0, T]$, $s$ is bounded from below by

$$ \epsilon_T = s(0) \exp(-[(\lambda_H + \lambda_W)n_{\max} + \beta_G]T). $$
In order to prove that there exists at most a unique solution $y$ for any initial condition $y^* \in V$, we will distinguish two cases.

First, if $s(0) = 0$, then it follows that $s(t) = 0$ for any $t \ge 0$, and hence $n_{S,I}^X(t) = 0$ for any $t \ge 0$ and any $X \in \{H, W\}, (S, I) \in \mathbb{S}$. Subsequently, the equation for $i$ reduces to $i'(t) = -\gamma i(t)$ on $\mathbb{R}_+$, ensuring uniqueness of $y$ on $\mathbb{R}_+$.

Second, if $s(0) > 0$, it follows from (i) and from our lower bound on $s$ over $[0, T]$ that $y(t) \in V_T = V \cap \{s \ge \epsilon_T\}$ for any $t \in [0, T]$. Our aim is to show that $f$ is Lipschitz continuous on $V_T$. Let $y = (s, i, n_{S,I}^X : X \in \{H, W\}, (S, I) \in \mathbb{S})$ and $\hat{y} = (\hat{s}, \hat{i}, \hat{n}_{S,I}^X : X \in \{H, W\}, (S, I) \in \mathbb{S})$ be two elements of $V$. First, consider $f_s$. Let $X \in \{H, W\}$, and define $c_X = \lambda_X \#\mathbb{S}(n_{\max})^2/m_X$, then

$$|\tau_X(y) - \tau_X(\hat{y})| \le \frac{\lambda_X}{m_X} \sum_{(S,I) \in \mathbb{S}} SI |n_{S,I}^X - \hat{n}_{S,I}^X| \le c_X \|y - \hat{y}\|_{\infty}.$$

It further holds that

$$|\tau_G(y)s - \tau_G(\hat{y})\hat{s}| \le \beta_G(|i||s - \hat{s}| + |i - \hat{i}||\hat{s}|) \le 2\beta_G \|y - \hat{y}\|_{\infty}.$$

Thus, letting $c_s = c_H + c_W + 2\beta_G$, it follows that

$$|f_s(y) - f_s(\hat{y})| \le c_s \|y - \hat{y}\|_{\infty}.$$

Similarly, letting $c_i = c_s + \gamma$,

$$|f_i(y) - f_i(\hat{y})| \le |f_s(y) - f_s(\hat{y})| + \gamma |i - \hat{i}| \le c_i \|y - \hat{y}\|_{\infty}.$$

It remains to focus on $f_{X,S,I}$. Proceeding as above, letting $c'_X = (2\lambda_X n_{\max} + 4\beta_G + \gamma)n_{\max}$, we obtain that

$$|f_{X,S,I}(y) - f_{X,S,I}(\hat{y})| \le c'_{X} \|y - \hat{y}\|_{\infty} + \left| \frac{\tau_{\bar{X}}(y)}{s} Sn_{S,I}^X - \frac{\tau_{\bar{X}}(\hat{y})}{\hat{s}} S\hat{n}_{S,I}^X \right| + \left| \frac{\tau_{\bar{X}}(y)}{s} (S+1)n_{S+1,I-1}^X - \frac{\tau_{\bar{X}}(\hat{y})}{\hat{s}} (S+1)\hat{n}_{S+1,I-1}^X \right| \mathbf{1}_{\{I \ge 1\}}.$$

Notice that as $y$ belongs to $V_T$, it follows that for $k_X = \lambda_X n_{\max} m_X$,

$$(\tau_{\bar{X}}(y)Sn_{S,I}^X, s) \in D_T = \{(x, y) : \epsilon_T \le y \le 1, 0 \le x \le k_X y^2\}.$$

Let $(x, y)$ and $(u, v)$ be two elements of $D_T$. It then holds that

$$\left|\frac{x}{y} - \frac{u}{v}\right| \le \frac{1}{v}\left(\frac{x}{y}|v-y| + |x-u|\right) \le \epsilon_T^{-1}(1 \vee k_X)(|v-y| + |x-u|).$$

Letting $k_{X,T} = \epsilon_T^{-1}(1 \vee k_X)(\lambda_X n_{\max}^2/m_X + 1)n_{\max}$, we obtain for any $X \in \{H, W\}$ and $(S, I) \in \mathbb{S}$,

$$\left| \frac{\tau_{\bar{X}}(y)}{s} Sn_{S,I}^X - \frac{\tau_{\bar{X}}(\hat{y})}{\hat{s}} S\hat{n}_{S,I}^X \right| \le \epsilon_T^{-1}(1 \vee k_X)(|\tau_{\bar{X}}(y)Sn_{S,I}^X - \tau_{\bar{X}}(\hat{y})S\hat{n}_{S,I}^X| + |s - \hat{s}|) \\ \le k_{X,T} \|y - \hat{y}\|_{\infty}.$$

As a consequence, we conclude that

$$|f_{X,S,I}(y) - f_{X,S,I}(\hat{y})| \le (c'_{X} + k_{X,T}) \|y - \hat{y}\|_{\infty}.$$

This establishes the desired Lipschitz continuity of $f$ on $V_T$, with associated Lipschitz constant $c_T = \max(c_s, c_i, c'_H + k_{H,T}, c'_W + k_{W,T})$.
Suppose now that there are two solutions $y$ and $\hat{y}$ of Equations (8a-c) such that $y(0) = \hat{y}(0)$. It then holds that

$$ \|y(T) - \hat{y}(T)\|_{\infty} \le \int_{0}^{T} \|f(y(t)) - f(\hat{y}(t))\|_{\infty} dt \le c_{T} \int_{0}^{T} \|y(t) - \hat{y}(t)\|_{\infty}. $$

Thus Gronwall's lemma ensures that $\|y(t) - \hat{y}(t)\|_{\infty} = 0$ for any $t \le T$. The desired conclusion on uniqueness follows.

(iii) In order to establish (iii), it remains to show that the initial condition $y^*$ defined by Equation (9) belongs to $V$. Let us start by noticing that, following Equation (9), for any $X \in \{H, W\}$,

$$ \sum_{(S,I) \in \mathcal{S}} n_{S,I}^{X}(0) = \sum_{n=2}^{n_{\max}} \pi_{n}^{X} \sum_{I=0}^{n-1} \binom{n}{I} \varepsilon^{I} (1-\varepsilon)^{n-I} = \sum_{n=2}^{n_{\max}} \pi_{n}^{X} (1-\varepsilon^{n}) \le 1. $$

Similarly,

$$ \begin{aligned} \sum_{(S,I) \in \mathcal{S}} S n_{S,I}^{X}(0) &= \sum_{n=2}^{n_{\max}} \pi_{n}^{X} \sum_{I=0}^{n-1} (n-I) \binom{n}{I} \varepsilon^{I} (1-\varepsilon)^{n-I} = \sum_{n=2}^{n_{\max}} \pi_{n}^{X} n(1-\varepsilon) \\ &= (m_{X} - \pi_{1}^{X})(1-\varepsilon) \le m_{X}s(0), \end{aligned} $$

where we have used the fact that $s(0) = 1 - \varepsilon$, and recognizing $\mathbb{E}[n - B]$ for $B \sim \mathcal{B}(n, \varepsilon)$ to deduce the second equality. The other conditions following immediately from (9), we conclude that $y^* \in V$. $\square$

**Acknowledgment.** I would like to thank my thesis supervisor Vincent Bansaye for his continuous support during the elaboration of this work and detailed feedback on the manuscript, Viet Chi Tran for fruitful discussions on tightness criteria, as well as Frank Ball and Marie Doumic for pertinent discussions on related topics. Finally, I am very grateful to the late Elisabeta Vergu for her guidance and stimulating remarks. She will be missed.

## References

[1] David Aldous. Stopping Times and Tightness. *The Annals of Probability*, 6(2):335-340, April 1978.

[2] Håkan Andersson and Tom Britton. Density dependent jump Markov processes. In *Stochastic Epidemic Models and Their Statistical Analysis*, Lecture Notes in Statistics. Springer, New York, NY, 2000.

[3] Frank G. Ball, Lorenzo Pellis, and Pieter Trapman. Reproduction numbers for epidemic models with households and other social structures II: Comparisons and implications for vaccination. *Mathematical Biosciences*, 274:108-139, April 2016.

[4] Frank G. Ball, David J. Sirl, and Pieter Trapman. Epidemics on random intersection graphs. *The Annals of Applied Probability*, 24(3):1081-1128, June 2014.

[5] Vincent Bansaye, François Deslandes, Madeleine Kubasch, and Elisabeta Vergu. The epidemiological footprint of contact structures in models with two levels of mixing. Available on arXiv: https://arxiv.org/abs/2303.05287, 2023.
[6] Rosanna C. Barnard, Istvan Z. Kiss, Luc Berthouze, and Joel C. Miller. Edge-Based Compartmental Modelling of an SIR Epidemic on a Dual-Layer Static-Dynamic Multiplex Network with Tunable Clustering. *Bulletin of Mathematical Biology*, 80(10):2698-2733, October 2018.

[7] Patrick Billingsley. *Convergence of Probability Measures*. Wiley Series in Probability and Statistics. Probability and Statistics Section. Wiley, New York, 2nd ed edition, 1999.

[8] Tom Britton and Etienne Pardoux. Chapter 3 A General Two-Level Mixing Model. In Tom Britton and Etienne Pardoux, editors, *Stochastic Epidemic Models with Inference*, Lecture Notes in Mathematics, pages 159-213. Springer International Publishing, Cham, 2019.

[9] Laurent Decreusefond, Jean-Stéphane Dhersin, Pascal Moyal, and Viet Chi Tran. Large graph limit for an SIR process in random network with heterogeneous connectivity. *The Annals of Applied Probability*, 22(2), April 2012.

[10] María del Valle Rafo, Juan Pablo Di Mauro, and Juan Pablo Aparicio. Disease dynamics and mean field models for clustered networks. *Journal of Theoretical Biology*, 526:110554, October 2021.

[11] Francesco Di Lauro, Luc Berthouze, Matthew D. Dorey, Joel C. Miller, and István Z. Kiss. The Impact of Contact Structure and Mixing on Control Measures and Disease-Induced Herd Immunity in Epidemic Models: A Mean-Field Model Perspective. *Bulletin of Mathematical Biology*, 83(11):117, November 2021.

[12] Zhilan Feng, Dashun Xu, and Haiyun Zhao. Epidemiological Models with Non-Exponentially Distributed Disease Stages and Applications to Disease Control. *Bulletin of Mathematical Biology*, 69(5):1511-1536, July 2007.

[13] Raphael Forien, Guodong Pang, and Etienne Pardoux. Recent Advances in Epidemic Modeling: Non-Markov Stochastic Models and their Scaling Limits. *The Graduate Journal of Mathematics*, 7(2):19-75, September 2022.

[14] Nicolas Fournier and Sylvie Méléard. A microscopic probabilistic description of a locally regulated population and macroscopic approximations. *The Annals of Applied Probability*, 14(4):1880-1919, November 2004.

[15] Laurent Hébert-Dufresne, Pierre-André Noël, Vincent Marceau, Antoine Allard, and Louis J. Dubé. Propagation dynamics on networks featuring complex topologies. *Physical Review E*, 82(3):036115, September 2010.

[16] Thomas House and Matt J. Keeling. Deterministic epidemic models with explicit household structure. *Mathematical Biosciences*, 213(1):29-39, May 2008.

[17] Thomas House and Matt J. Keeling. Insights from unifying modern approximations to infections on networks. *Journal of the Royal Society Interface*, 8(54):67-73, January 2011.

[18] Karly A. Jacobsen, Mark G. Burch, Joseph H. Tien, and Grzegorz A. Rempała. The large graph limit of a stochastic epidemic model on a dynamic multilayer network. *Journal of Biological Dynamics*, 12(1):746-788, January 2018.

[19] Svante Janson, Malwina Luczak, and Peter Windridge. Law of large numbers for the SIR epidemic on a random graph with given degrees. *Random Structures & Algorithms*, 45(4):726-763, 2014.
[20] Anatole Joffe and Michel Metivier. Weak convergence of sequences of semimartingales with applications to multitype branching processes. *Advances in Applied Probability*, 18(1):20–65, March 1986.

[21] Benjamin Jourdain, Sylvie Méléard, and Wojbor A. Woyczynski. Lévy flights in evolutionary ecology. *Journal of Mathematical Biology*, 65(4):677–707, October 2012.

[22] Istvan Z. Kiss, Eben Kenah, and Grzegorz A. Rempala. Necessary and sufficient conditions for exact closures of epidemic equations on configuration model networks, August 2022.

[23] István Z. Kiss, Joel C. Miller, and Péter L. Simon. Mathematics of Epidemics on Networks: *From Exact to Approximate Models*, volume 46 of *Interdisciplinary Applied Mathematics*. Springer International Publishing, Cham, 2017.

[24] Thomas G. Kurtz. Epidemic models. In *Approximation of Population Processes*, number 36 in CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, Pa, 1981.

[25] Alun L. Lloyd. Realistic Distributions of Infectious Periods in Epidemic Models: Changing Patterns of Persistence and Dynamics. *Theoretical Population Biology*, 60(1):59–71, August 2001.

[26] Giancarlo De Luca, Kim Van Kerckhove, Pietro Coletti, Chiara Poletto, Nathalie Bossuyt, Niel Hens, and Vittoria Colizza. The impact of regular school closure on seasonal influenza epidemics: A data-driven spatial transmission model for Belgium. *BMC Infectious Diseases*, 18(1):29, January 2018.

[27] Junling Ma, Pauline van den Driessche, and Frederick H. Willeboordse. Effective degree household network disease model. *Journal of Mathematical Biology*, 66(1):75–94, January 2013.

[28] Sylvie Meleard and Vincent Bansaye. *Stochastic Models for Structured Populations: Scaling Limits and Long Time Behavior*. Springer International Publishing, Cham, 2015.

[29] Alba Mendez-Brito, Charbel El Bcheraoui, and Francisco Pozo-Martin. Systematic review of empirical studies comparing the effectiveness of non-pharmaceutical interventions against COVID-19. *Journal of Infection*, 83(3):281–293, September 2021.

[30] Joel C. Miller. Percolation and epidemics in random clustered networks. *Physical Review E*, 80(2):020901, August 2009.

[31] Joel C Miller. A note on a paper by Erik Volz: SIR dynamics in random networks. *Journal of Mathematical Biology*, 62:349–358, 2011.

[32] Mark E. J. Newman. Random graphs with clustering. *Physical Review Letters*, 103(5):058701, July 2009.

[33] Lorenzo Pellis, Neil M. Ferguson, and Christophe Fraser. Threshold parameters for a model of epidemic spread among households and workplaces. *Journal of The Royal Society Interface*, 6(40):979–987, November 2009.

[34] Lorenzo Pellis, Neil M. Ferguson, and Christophe Fraser. Epidemic growth rate and household reproduction number in communities of households, schools and workplaces. *Journal of Mathematical Biology*, 63(4):691–734, October 2011.
[35] Sylvie Roelly-Coppoletta. A criterion of convergence of measure-valued processes: Application to measure branching processes. *Stochastics: An International Journal of Probability and Stochastic Processes*, 17(1-2):43–65, 1986.

[36] Neil Sherborne, Joel C. Miller, Konstantin B. Blyuss, and Istvan Z. Kiss. Mean-field models for non-Markovian epidemics on networks. *Journal of Mathematical Biology*, 76(3):755–778, February 2018.

[37] Mario Ignacio Simoy and Juan Pablo Aparicio. Socially structured model for COVID-19 pandemic: Design and evaluation of control measures. *Computational and Applied Mathematics*, 41(1):14, December 2021.

[38] Guillaume St-Onge, Laurent Hébert-Dufresne, and Antoine Allard. Heterogeneous transmission in groups induces a superlinear force of infection, February 2023.

[39] Viet Chi Tran. *Modèles particulaires stochastiques pour des problèmes d'évolution adaptative et pour l'approximation de solutions statistiques*. PhD thesis, Université de Nanterre - Paris X, December 2006.

[40] Viet Chi Tran. *Une ballade en forêts aléatoires*. Habilitation a Diriger les Recherches, Université Lille 1, November 2014.

[41] Erik M. Volz. SIR dynamics in random networks with heterogeneous connectivity. *Journal of Mathematical Biology*, 56(3):293–310, March 2008.

[42] Erik M. Volz, Joel C. Miller, Alison Galvani, and Lauren Ancel Meyers. Effects of Heterogeneous and Clustered Contact Patterns on Infectious Disease Dynamics. *PLoS Computational Biology*, 7(6):e1002042, June 2011.

[43] Robert R. Wilkinson, Frank G. Ball, and Kieran J. Sharkey. The relationships between message passing, pairwise, Kermack–McKendrick and stochastic SIR epidemic models. *Journal of Mathematical Biology*, 75(6-7), April 2017.