$\\{(S,I+1):(S,I)\in\mathbb{S},S+I<n_{\max}\\}=\\{(S,I)\in\mathbb{S}:S+I\geq
3,I\geq 1\\}.$

Further, both integrands are bounded by a constant which does not depend on
$t$. The continuity of $t\mapsto\langle\eta^{X}_{t},G^{(S,I)}_{T,t}\rangle$ on
$[0,T]$ follows.

We are now ready to reason by induction. Notice that

$\\{(S,I)\in\mathbb{S}:S+I\geq 3,I\geq
1\\}=\bigcup_{k=3}^{n_{\max}}\bigcup_{\ell=1}^{k-1}\\{(k-\ell,\ell)\\}.$

For any $k\in\llbracket 3,n_{\max}\rrbracket$, for any $\ell\in\llbracket
1,k-1\rrbracket$, Equation (34) becomes

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

Let us start with the case $k=n_{\max}$. If further $\ell=1$, Equation (36)
reduces to

$\lvert\langle\eta^{X}_{T},G^{(n_{\max}-1,1)}\rangle\rvert\leq
C\int_{0}^{T}\lvert\langle\eta^{X}_{t},G^{(n_{\max}-1,1)}_{T,t}\rangle\rvert
dt.$

As $G^{(n_{\max}-1,1)}=G^{(n_{\max}-1,1)}_{T,T}$, the desired conclusion
follows from Gronwall’s lemma. Suppose now that Equation (35) holds for
$(S,I)=(n_{\max}-\ell+1,\ell-1)$ with $\ell\in\llbracket
1,n_{\max}-1\rrbracket$. As further $k=n_{\max}$, we may notice that Equation
(36) simplifies to

$\lvert\langle\eta^{X}_{T},G^{(n_{\max}-\ell,\ell)}\rangle\rvert\leq
C\int_{0}^{T}\lvert\langle\eta^{X}_{t},G^{(n_{\max}-\ell,\ell)}_{T,t}\rangle\rvert
dt.$

Again, Gronwall’s lemma implies that Equation (35) holds for
$(S,I)=(n_{\max}-\ell,\ell)$. We thus have established that Equation (35) is
true for $(S,I)=(n_{\max}-\ell,\ell)$ for all $\ell\in\llbracket
0,n_{\max}-1\rrbracket$.

Les us now assume that Equation (35) is satisfied for
$(S,I)=(k+1-\ell^{\prime},\ell^{\prime})$ for some $k\in\llbracket
3,n_{\max}-1\rrbracket$ and for all $\ell^{\prime}\in\llbracket
1,k-1\rrbracket$. Then in particular,
$\langle\eta^{X}_{t},G^{(k,1)}_{T,t}\rangle=0$ for all $t\in[0,T]$. Hence
proceeding like before leads to the conclusion that
$\langle\eta^{X}_{t},G^{(k-1,1)}_{T,t}\rangle=0$ for all $t\in[0,T]$ as the
first two terms of the integrand of the associated inequality derived from
Equation (36) are null. The same procedure allows to show that if Equation
(35) holds for $(S,I)=(k-\ell+1,\ell-1)$ for some $\ell\in\llbracket
2,k-1\rrbracket$, then using the induction hypothesis, it also holds for
$(S,I)=(k-\ell,\ell)$.

This concludes the proof. ∎

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

###### Proof of Theorem 2.2..

Before concluding, we need to emphasize that it would have been possible to
chose $X=W$ when replacing $S$ and $I$ by $K_{X}S^{X}$ and $K_{X}I^{X}$ for
$\mathcal{I}_{G}$ in Proposition 3.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 3.12 and 3.13 show that for any
$X\in\\{H,W\\}$,

$\left(\frac{s^{X}}{m_{X}},\frac{i^{X}}{m_{X}},n^{X}_{(S,I)}:(S,I)\in\mathbb{S},n^{\overline{X}}_{(S,I)}:(S,I)\in
C_{\overline{X}}\right)$

satisfies the Cauchy problem (15) with initial condition (7). However,
Proposition 3.8 ensures uniqueness of the solutions to this Cauchy problem. It
hence is sensible to define, for $t\geq 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 eqs. 6b to 6f with initial conditions (7), and
concludes the proof. ∎

## 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 approximation, and the
associated convergence of the stochastic model was established. In the general
case of an arbitrary distribution $\nu$ of the length of infectious periods,
the stochastic model converges under the large graph limit to the unique
solution of a measure-valued equation. In the particular case where $\nu$ is
an exponential distribution, we have shown that the epidemic dynamics can
actually be reduced to a closed, finite-dimensional dynamical system. These
results further are of interest since, to our knowledge, none of the
previously proposed model reductions for epidemic models on networks with
cliques, let alone those which in addition consider two levels of mixing, have
been proven to coincide with the large graph limit of their stochastic
counterparts.

Regarding the implementation of the reduced model, we have focused on the case
where $\nu$ is the exponential distribution. Comparing the solution to the
dynamical system with stochastic simulations of epidemic trajectories shows
good accordance in practice, as expected. We also have emphasized that the
dynamical system can be implemented in an automatic way, which implies that
despite its high number of equations, it can be used even for larger structure
sizes. Finally, some numerical evidence points to the reduced model being
significantly less demanding than stochastic simulations in terms of
computational cost for sufficiently large values of the reproduction number,
which means that it is of interest for numerical explorations. A broader study
of computational cost would be pertinent to refine these results.

Notice here that 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 is
the case for households and workplaces, the adaptation of the aforementioned
results is expected to be straightforward. In particular, the dimension of the
corresponding dynamical system should still be of order $O(n_{\max}^{2})$,
implying that the model may remain tractable.

Finally, we have compared the reduced model obtained by the large population
approximation with the corresponding EBCM in the line of [29]. 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.

Acknowledgment. I would like to thank Vincent Bansaye for his continuous
support during the elaboration of this work and detailed feedback on the
manuscript, as well as Viet Chi Tran and Frank Ball for discussions on related
topics. Finally, I am very grateful to 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] 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.
  * [3] 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.
  * [4] Vincent Bansaye, François Deslandes, Madeleine Kubasch, and Elisabeta Vergu.  The epidemiological footprint of contact structures in models with two levels of mixing, March 2023.
  * [5] 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\.
  * [6] Patrick Billingsley.  Convergence of Probability Measures.  Wiley Series in Probability and Statistics. Probability and Statistics Section. Wiley, New York, 2nd ed edition, 1999.
  * [7] 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\.
  * [8] 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.
  * [9] 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.
  * [10] 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\.
  * [11] Thomas House and Matt J. Keeling.  Deterministic epidemic models with explicit household structure.  Mathematical Biosciences, 213(1):29–39, May 2008.
  * [12] 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.
  * [13] 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.
  * [14] 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.
  * [15] 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.
  * [16] 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.
  * [17] Sylvie Meleard and Vincent Bansaye.  Stochastic Models for Structured Populations: Scaling Limits and Long Time Behavior.  Springer International Publishing, Cham, 2015.
  * [18] 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.
  * [19] Joel C. Miller.  Percolation and epidemics in random clustered networks.  Physical Review E, 80(2):020901, August 2009.
  * [20] Joel C Miller.  A note on a paper by Erik Volz: SIR dynamics in random networks.  Journal of Mathematical Biology, 62:349–358, 2011.
  * [21] Mark E. J. Newman.  Random graphs with clustering.  Physical Review Letters, 103(5):058701, July 2009.
  * [22] 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.
  * [23] 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.
  * [24] 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, May 2007.
  * [25] 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.
  * [26] Guillaume St-Onge, Laurent Hébert-Dufresne, and Antoine Allard.  Heterogeneous transmission in groups induces a superlinear force of infection, February 2023.
  * [27] 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.
  * [28] Erik M. Volz.  SIR dynamics in random networks with heterogeneous connectivity.  Journal of Mathematical Biology, 56(3):293–310, March 2008.
  * [29] 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.

##  Appendix A Implementation of the large population approximation

###  A.1 Automatic implementation of the dynamical system

It is possible to implement dynamical system eqs. 6b to 6f 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 [23], for another purpose. With the previous structure of
$\mathbb{S}$, one may for instance notice that for any $n\in\llbracket
2,n_{\max}\rrbracket$ and $i\in\llbracket 0,n-1\rrbracket$, the state $(n,i)$
is the $c(n,i)$-th state enumerated in $\mathbb{S}$, where
$c(n,i)=(n-1)n/2+i$. As a consequence, the general expression of Equation (6f)
may be used to handle all the dynamics of the functions $n^{X}_{(S,I)}$, 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 eqs. 6b to 6f in
comparison to stochastic simulations using 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 form 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 average runtimes over all hundred runs of
both the script of interest and the reference function are computed. As
computation times tend to fluctuate between executions of a given script, the
average runtime of the script of interest is subsequently normalised by the
average runtime of the reference function. Comparison of the computation times
for the stochastic and the reduced model is then based on the comparison of
those normalised average 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 $\beta_{G}$, $\lambda_{H}$, $\lambda_{W}$ and the removal rate $\gamma$,
as well as the initial proportion of infected $\varepsilon$ and the time
interval $[0,T]$ on which the epidemic is simulated. For the stochastic model,
the population size $N$ will be fixed to ten thousand individuals. For all
scenarios considered here, the structure size distributions will be those of
Figure 1. Further, the initial proportion of individuals will be set to
$\varepsilon=0.005$, as Figure 2 indicates that for the stochastic model, the
risk of extinction at the beginning of the epidemic due to stochastic
fluctuations is small when starting from $50$ infected at time $0$, even for
low values of $R_{0}$. However, 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 $\gamma$ 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 A1: Proportion of susceptible (S) and infected (I) in the population,
for each scenario detailed in Table A1, as given by dynamical system eqs. 6b
to 6f. 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.

Figure A1 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 A1. Notice that in
particular, this includes the parameters of Figure 2.

| $R_{0}=1.2$ | $R_{0}=1.4$ | $R_{0}=1.7$ | $R_{0}=2.0$ | $R_{0}=2.5$ | $(p_{G},p_{H},p_{W})$
---|---|---|---|---|---|---
$\beta_{G}$ | 0.03 | 0.035 | 0.045 | 0.05 | 0.06 | (0.2, 0.4, 0.4)
$\lambda_{H}$ | 0.05 | 0.07 | 0.09 | 0.15 | 0.2
$\lambda_{W}$ | 0.0015 | 0.0016 | 0.0018 | 0.002 | 0.0022
$T$ | 130 | 130 | 105 | 85 | 75
$\beta_{G}$ | 0.06 | 0.07 | 0.085 | 0.1 | 0.125 | (0.4, 0.4, 0.2)
$\lambda_{H}$ | 0.06 | 0.07 | 0.1 | 0.15 | 1.5
$\lambda_{W}$ | 0.00075 | 0.0008 | 0.001 | 0.0011 | 0.00115
$T$ | 145 | 130 | 95 | 80 | 55
Table A1: Values of the contact rates $\beta_{G}$, $\lambda_{H}$ and
$\lambda_{W}$ considered, grouped by value of $R_{0}$ and proportions of
infections per layer $(p_{G},p_{H},p_{W})$ characterizing the scenarios.

Let us now turn to the results. For each scenario of Table A1, measurement of
normalised average computation times was repeated three times. The results are
shown in Figure A2, which indicates for each scenario the ratio of the
normalised average runtime for one resolution of dynamical system eqs. 6b to
6f over the normalised average 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})=(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 eqs. 6b to
6f 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 A2 shows that the associated results are very similar to those
obtained previously, pleading against this hypothesis.

Figure A2: Ratio of the normalised average computation time for solving once
dynamic system eqs. 6b to 6f over the normalised average computation time for
one stochastic simulation (SSA). This ratio was computed three times for each
scenario of Table A1. The results are presented in function of $R_{0}$, while
colors indicate the value of $(p_{G},p_{H},p_{W})$. Unless stated otherwise,
the parameter $T$ from Table A1 was used. The dotted line indicates the
threshold of one.

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

##  Appendix B Edge-based compartmental model

###  B.1 Presentation of the EBCM

Figure B1: Comparison of the stochastic model with the large population
approximation given by dynamical system eqs. 6b to 6f and the corresponding
EBCM. Household and workplace distributions are those of Figure 1. Epidemic
parameters were 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
was used to simulate $50$ trajectories of the stochastic model defined in
Proposition 1.1 in a population of $K=10000$ individuals (faint lines; the
same 50 realizations of the population structure as in Figure 2 were used for
each subplot). For Panel (a), only trajectories reaching a threshold of
$0.005$ infected were kept, and time was shifted so that time $0$ corresponds
to the moment when this threshold is reached. Finally, the deterministic
solution $(s,i)$ of both dynamical system eqs. 6b to 6f (thick lines) and the
EBCM (dotted lines) are represented for each scenario. For Panel (a), the same
time shifting procedure as for simulations is applied.

Let us start by describing how to obtain the population structure of the local
level of mixture described in Section 1.1 using a 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}^{X}_{n}=n\pi^{X}_{n}/n^{X}$, for any
$n\in\llbracket 1,n_{\max}\rrbracket$ 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\llbracket
1,n_{\max}\rrbracket$, 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^{X}_{n}(t)$ for $X\in\\{H,W\\}$ and $n\in\llbracket
1,n_{\max}\rrbracket$ 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\geq 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
leads to

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

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

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

$m^{X}_{n}=\theta^{G}\hat{\pi}^{X}_{n}\theta^{X}_{n}\sum_{k=1}^{n_{\max}}\hat{\pi}^{\overline{X}}_{k}\theta^{\overline{X}}_{k},$

which corresponds to the proportion of individuals who are susceptible and
belong to a structure of type $X$ and size $n$. Also, let $n^{X}_{(S,I,R)}$ 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^{X}_{n}=\lambda_{X}\sum_{\begin{subarray}{c}(S,I,R)\in\mathbb{N}^{3}\\\
S+I+R=n\end{subarray}}SI\,n^{X}_{(S,I,R)}\quad\text{and}\quad\tau^{X}_{n}=\left(\beta_{G}i+\sum_{k=1}^{n_{\max}}T^{\overline{X}}_{k}\right)\frac{\hat{\pi}^{X}_{n}\theta^{X}_{n}}{\sum_{k=1}^{n_{\max}}\hat{\pi}^{X}_{k}\theta^{X}_{k}}.$

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\llbracket
2,n_{\max}\rrbracket,\;\frac{d}{dt}\theta^{X}_{n}=-\frac{T^{X}_{n}}{m^{X}_{n}}\theta^{X}_{n}.$

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

$\displaystyle\frac{d}{dt}n^{X}_{(S,I,R)}$
$\displaystyle=-\left(\lambda_{X}SI+\frac{\tau^{X}_{n}}{m^{X}_{n}}S+\gamma
I\right)n^{X}_{(S,I,R)}$
$\displaystyle+\gamma(I+1)n^{X}_{(S,I+1,R-1)}\mathbf{1}_{\left\\{R\geq
1\right\\}}$
$\displaystyle+\left(\lambda_{X}(S+1)(I-1)++\frac{\tau^{X}_{n}}{m^{X}_{n}}(S+1)\right)n^{X}_{(S+1,I-1,R)}\mathbf{1}_{\left\\{I>1\right\\}}.$

Additionally, as in a structure of size one, no infection may occur within the
structure itself, $\theta^{H}_{1}$ and $\theta^{W}_{1}$ are constant over
time. Finally, it remains to define the initial conditions. Following Volz _et
al._ [29], 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\llbracket
1,n_{\max}\rrbracket$ and $I\in\llbracket 1,n-1\rrbracket$,

$\displaystyle i(0)$ $\displaystyle=\varepsilon,$ (B1)
$\displaystyle\theta^{G}(0)=\theta^{X}_{n}(0)$ $\displaystyle=1-\varepsilon,$
$\displaystyle n^{X}_{(n,0,0)}$
$\displaystyle=\frac{1}{n}\hat{\pi}^{X}_{n}(1-\varepsilon)^{n},$
$\displaystyle n^{X}_{(n-I,I,0)}$
$\displaystyle=\hat{\pi}^{X}_{n}(1-\varepsilon)^{n-I}\varepsilon^{I}.$

The proportions of susceptible and infected as predicted by both the EBCM and
dynamical system eqs. 6b to 6f are shown in Figure B1, for different values of
$\varepsilon$. Let us first notice that in the case $\varepsilon=0.001$,
corresponding to Panel (a) of Figure B1, the solutions $(s,i)$ of both the
EBCM and dynamical system eqs. 6b to 6f are in perfect accordance, emphasizing
the fact that for very small 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. At first, as
for $\varepsilon=0.01$ in Panel (b) of Figure B1, the difference is mainly
visible at the beginning of the epidemic and seems to dampen afterwards. For
even larger values of $\varepsilon$, the solution to the EBCM does not succeed
in reproducing the dynamics over time, as shown in Panel (c) of Figure B1. The
problem for capturing the epidemic dynamics for higher values of $\varepsilon$
likely 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. Of
course, it would be possible to try to correct for this behaviour by adapting
the initial condition, but it is not clear which initial condition would be a
sensible choice.

###  B.2 Computational performance

In order to compare the computation times needed to solve either dynamical
system eqs. 6b to 6f 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 A2, this is not expected to significantly affect
the qualitative result.

Structure size distribution |  Figure 1
---|---
Epidemic parameters | $(\beta_{G},\lambda_{H},\lambda_{W},\gamma)=(0.125,1.5,0.00115,0.125)$
Initial proportion of infected | $\varepsilon=0.001$
Time interval $[0,T]$ | $T=30$
Average runtime (100 runs) |  Dynamical system eqs. 6b to 6f | 6 seconds
Reference function | 39 seconds
Average runtime (45 runs) | EBCM | 33266 seconds (9.2 hours)
reference function | 36 seconds
Table B2: Numerical assessment of the computation time needed to solve either
dynamical system eqs. 6b to 6f or the EBCM introduced in Appendix B.1. The
model parameters used are indicated together with the average runtimes.

The model parameters and the associated average runtimes are shown in Table
B2. Due to the excessive computation needed to solve the EBCM, only 45 runs of
this script were performed. However, considering that the normalised average
runtime for solving the EBCM is several thousand times higher than the
normalised average runtime for solving dynamical system eqs. 6b to 6f, this
again is not expected to significantly alter the results. Finally, as shown in
Figure B2, the computation times necessary for solving the EBCM are relatively
homogeneous over all 45 runs, always exceeding eight hours, indicating that
the average computation time is not biased by an outlier.

Figure B2: Evaluation of the time needed to numerically solve the EBCM
introduced in Section 2.2.2, using the ODE solver odeint from the
scipy.integrate library. The script has been executed independently 45 times,
and both the function of interest and a common reference function have been
timed for each execution. The model parameters are the same as in Panel (a) of
Figure 2.

##  Appendix C Proofs

###  C.1 Proof of Proposition 3.8

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

######  Lemma C.1.

Consider a solution $y$ of the dynamical system eqs. 6b to 6f and let
$\Delta(t)=m_{X}s(t)-\sum_{(S,I)\in\mathbb{S}}Sn^{X}_{(S,I)}(t)$. Then

$\small\frac{d}{dt}\Delta(t)=\gamma
n^{X}_{(1,1)}(t)-\Big{(}\tau_{G}(t)+\frac{1}{s(t)}\tau_{\overline{X}}(t)\Big{)}\Delta(t).$

###### Proof of Lemma C.1..

Let $X\in\\{H,W\\}$. First, notice that

$\\{(S+1,I-1):(S,I)\in\mathbb{S},I\geq 1\\}=\\{(S,I)\in\mathbb{S}:S>1\\}.$

As $S-1=0$ whenever $S=1$, we thus obtain that

$\displaystyle\sum_{(S,I)\in\mathbb{S}}$ $\displaystyle
S^{2}In^{X}_{(S,I)}-\sum_{(S,I)\in\mathbb{S}}S(S+1)(I-1)n^{X}_{(S+1,I-1)}\mathbf{1}_{\left\\{I\geq
1\right\\}}$ (C1)
$\displaystyle=\sum_{(S,I)\in\mathbb{S}}S^{2}In^{X}_{(S,I)}-\sum_{(S,I)\in\mathbb{S}}(S-1)SIn^{X}_{(S,I)}=\sum_{(S,I)\in\mathbb{S}}SIn^{X}_{(S,I)}.$

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

$\sum_{(S,I)\in\mathbb{S}}\gamma
SIn^{X}_{(S,I)}-\sum_{(S,I)\in\mathbb{S}}\gamma
S(I+1)n^{X}_{(S,I+1)}\mathbf{1}_{\left\\{S+I<n_{\max}\right\\}}=\gamma
n^{X}_{(1,1)}.$ (C2)

The desired conclusion then results directly from Equations eqs. 6b to 6f,
regrouping the terms of the form of Equations (C1) and (C2) in order to
simplify the expression. ∎

We are now ready to focus on the desired result.

###### Proof of proposition 3.8..

(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}\geq 0$ be such that $y(t_{0})\in V$. If $i(t_{0})=0$, then
$i^{\prime}(t_{0})=-s^{\prime}(t_{0})\geq 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
(16), respectively.

Regarding $n^{X}_{(S,I)}$ for $X\in\\{H,W\\}$ and $(S,I)\in\mathbb{S}$, let us
first emphasize that, by inequality (16), the derivative of $n^{X}_{(S,I)}$ is
well defined at $t_{0}$. Indeed, recall from Equation (6f) that the terms
$n^{X}_{(S,I)}(t)/s(t)$ and $n^{X}_{(S+1,I-1)}(t)/s(t)$ may be ill defined,
due to the division by $s(t)$. However, inequality (16) ensures that these
ratios stay bounded at all time. As a consequence, we may now notice as
previously that, if $n^{X}_{(S,I)}(t_{0})=0$, Equation (6f) ensures that
$\frac{d}{dt}n^{X}_{(S,I)}(t_{0})\geq 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^{X}_{(S,I)}(t)=-\gamma
n^{X}_{(1,1)}(t)\leq 0.$

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

$\frac{d}{dt}n^{X}_{(S,I)}(t_{0})=-\left(\lambda_{X}SI+\tau_{\overline{X}}(t_{0})\frac{S}{s(t_{0})}+\tau_{G}(t_{0})S+\gamma
I\right)\leq 0,$

from which one may deduce that $n^{X}_{(S,I)}$ 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) In order to prove the existence of a unique global solution $y$ for any
initial condition $y^{*}\in V$, it is enough to show that $f$ is Lipschitz
continuous on $V$, as global existence and uniqueness then follow from Cauchy-
Lipschitz Theorem.

Let us first consider $f_{s}$. Let
$y=(s,i,n^{X}_{(S,I)}:X\in\\{H,W\\},(S,I)\in\mathbb{S})$ and
$\hat{y}=(\hat{s},\hat{i},\hat{n}^{X}_{(S,I)}:X\in\\{H,W\\},(S,I)\in\mathbb{S})$
be two elements of $V$, and let $X\in\\{H,W\\}$. Letting
$c_{X}=\lambda_{X}\\#\mathbb{S}(n_{\max})^{2}/m_{X}$, we obtain that

$\lvert\tau_{X}(y)-\tau_{X}(\hat{y})\rvert\leq\frac{\lambda_{X}}{m_{X}}\sum_{(S,I)\in\mathbb{S}}SI\lvert
n^{X}_{(S,I)}-\hat{n}^{X}_{(S,I)}\rvert\leq c_{X}\left\lVert
y-\hat{y}\right\rVert_{\infty}.$

Further, recall that $0\leq\tau_{G}(y)\leq\beta_{G}$ for any $y\in V$.

Letting $c_{s}=c_{H}+c_{W}+\beta_{G}$, it follows that

$\lvert
f_{s}(y)-f_{s}(\hat{y})\rvert\leq\lvert\tau_{H}(y)-\tau_{H}(\hat{y})\rvert+\lvert\tau_{W}(y)-\tau_{W}(\hat{y})\rvert+\lvert\tau_{G}(y)s-\tau_{G}(\hat{y})\hat{s}\rvert\leq
c_{s}\left\lVert y-\hat{y}\right\rVert_{\infty}.$

Similarly, letting $c_{i}=c_{s}+\gamma$,

$\lvert f_{i}(y)-f_{i}(\hat{y})\rvert\leq\lvert
f_{s}(y)-f_{s}(\hat{y})\rvert+\gamma\lvert i-\hat{i}\rvert\leq
c_{i}\left\lVert y-\hat{y}\right\rVert_{\infty}.$

Finally, notice that inequality (16) implies that
$\max(Sn^{X}_{(S,I)}/s,S\hat{n}^{X}_{(S,I)}/\hat{s})\leq m_{X}$ for any
$X\in\\{H,W\\}$ and $(S,I)\in\mathbb{S}$. As a consequence, we obtain that

$\lvert\tau_{\overline{X}}(y)\frac{Sn^{X}_{(S,I)}}{s}-\tau_{\overline{X}}(y)\frac{S\hat{n}^{X}_{(S,I)}}{\hat{s}}\rvert\leq
m_{X}\lvert\tau_{\overline{X}}(y)-\tau_{\overline{X}}(\hat{y})\rvert\leq
m_{X}c_{\overline{X}}\left\lVert y-\hat{y}\right\rVert_{\infty},$

from which we may deduce, letting
$c=2(\lambda_{X}(n_{\max})^{2}+m_{X}c_{\overline{X}}+\beta_{G}+\gamma)$, that

$\lvert f_{X;(S,I)}(y)-f_{X;(S,I)}(\hat{y})\rvert\leq c\left\lVert
y-\hat{y}\right\rVert_{\infty}.$

This establishes the desired Lipschitz continuity of $f$ on $V$.

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

$\sum_{(S,I)\in\mathbb{S}}n^{X}_{(S,I)}(0)=\sum_{n=2}^{n_{\max}}\pi^{X}_{n}\sum_{I=0}^{n-1}\binom{n}{I}\varepsilon^{I}(1-\varepsilon)^{n-I}=\sum_{n=2}^{n_{\max}}\pi^{X}_{n}(1-\varepsilon^{n})\leq
1.$

Similarly,

$\displaystyle\sum_{(S,I)\in\mathbb{S}}Sn^{X}_{(S,I)}(0)$
$\displaystyle=\sum_{n=2}^{n_{\max}}\pi^{X}_{n}\sum_{I=0}^{n-1}(n-I)\binom{n}{I}\varepsilon^{I}(1-\varepsilon)^{n-I}=\sum_{n=2}^{n_{\max}}\pi^{X}_{n}n(1-\varepsilon)$
$\displaystyle=(m_{X}-\pi^{X}_{1})(1-\varepsilon)\leq m_{X}s(0),$

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 (7), we conclude
that $y^{*}\in V$. ∎

###  C.2 Proof of Lemma 3.9

Let $X\in\\{H,W\\}$. To begin with, let us introduce the family of sets

$E_{(n,s)}=\\{x\in
E:\mathbf{n}(x)=n,\mathbf{s}(x)=s\\},\;\;\forall(n,s)\in\llbracket
1,n_{\max}\rrbracket\times\llbracket 0,n\rrbracket$

which constitutes a partition of $E$. In the following, let
$(n_{t},s_{t},\tau_{t})$ be distributed according to $\eta^{X}_{t}$, for any
$t\geq 0$. Let $f$ be a non-negative measurable function on $E$. Then

$\displaystyle\langle\eta^{X}_{t},f\rangle$
$\displaystyle=\mathbb{E}\left[f(n_{t},s_{t},\tau_{t})\right]$
$\displaystyle=\sum_{n=1}^{n_{\max}}\sum_{s=0}^{n}\mathbb{E}[f(n_{t},s_{t},\tau_{t})|n_{t}=n,s_{t}=s]\mathbb{P}(n_{t}=n,s_{t}=s).$

Conditionally on $\\{n_{t}=n,s_{t}=s\\}$, almost surely
$(n_{t},s_{t},\tau_{t})\in E_{(n,s)}$ and
$(\tau_{t})_{n-s+1:n_{\max}}=0_{n_{\max}-(n-s)}$. Thus

$\displaystyle\langle\eta^{X}_{t},f\rangle$
$\displaystyle=\sum_{n=1}^{n_{\max}}f(n,n,0_{n_{\max}})\mathbb{P}(n_{t}=n,s_{t}=s)$
(C3)
$\displaystyle+\sum_{n=1}^{n_{\max}}\sum_{s=0}^{n-1}\mathbb{E}\left[f\left(n,s,\sum_{k=1}^{n-s}(\tau_{t})_{k}e_{k}\right)\Big{|}n_{t}=n,s_{t}=s\right]\mathbb{P}(n_{t}=n,s_{t}=s).$

Let $p^{X}_{t}(n,s)=\mathbb{P}(n_{t}=n,s_{t}=s)$. Then on the one hand, define
for any $n\in\llbracket 1,n_{\max}\rrbracket$,

$\rho^{X,n,n}_{t}=p^{X}_{t}(n,n).$

On the other hand, for $s<n$, let $\tilde{\eta}^{X,(n,s)}_{t}$ be the
distribution of $((\tau_{t})_{1},\dots,(\tau_{t})_{n-s})$ conditionally on
$\\{n_{t}=n,s_{t}=s\\}$. If we manage to show that
$\tilde{\eta}^{X,(n,s)}_{t}$ is absolutely continuous with regard to the
Lebesgue measure on $\mathbb{R}^{n-s}$, then letting
$\tilde{\rho}^{X,(n,s)}_{t}$ be its density, the result follows with

$\forall\tau\in\mathbb{R}^{n-s},\;\rho^{X,n,s}_{t}(\tau)=p^{X}_{t}(n,s)\tilde{\rho}^{X,(n,s)}_{t}(\tau).$

In order to show that for any $n\in\llbracket 1,n_{\max}\rrbracket$ and
$s\in\llbracket 1,n-1\rrbracket$, $\tilde{\eta}^{X,(n,s)}_{t}$ is absolutely
continuous with regard to the Lebesgue measure on $\mathbb{R}^{n-s}$, we will
proceed by induction.

Let us first consider the case $s=n-1$. We want to apply Equation (4) to test
functions $\varphi^{*}\in\mathcal{C}^{1}_{b}(E,\mathbb{R}_{+})$ which are
equal to zero outside of the set $\\{x\in
E:\mathbf{n}(x)=n,\mathbf{s}(x)=n-1\\}$. Hence, let
$\varphi\in\mathcal{C}^{1}_{b}(\mathbb{R})$ be non-negative, and define
$\varphi^{*}\in\mathcal{C}^{1}_{b}(E,\mathbb{R}_{+})$ by

$\forall x\in
E,\quad\varphi^{*}(x)=\varphi(\tau_{1}(x))\mathbf{1}_{\left\\{(\mathbf{n}(x),\mathbf{s}(x))=(n,s)\right\\}}.$

Let $T>0$ and for $t\in[0,T]$, define $f_{t}(x)=\varphi^{*}(\Psi(x,T,t))$ for
$x\in E$. Using Equation (C3) and applying Equation (4) to $f_{t}$ after
noticing that $\mathcal{A}f_{t}(x)=0$ for any $(t,x)\in\mathbb{R}_{+}\times
E$, we obtain that

$\displaystyle p^{X}_{t}(n,n-1)$
$\displaystyle\langle\tilde{\eta}^{X,(n,n-1)}_{T},\varphi\rangle=\langle\eta^{X}_{T},\varphi^{*}\rangle=\langle\eta^{X}_{T},f_{T}\rangle$
$\displaystyle=\langle\eta^{X}_{0},f_{0}\rangle+\lambda_{H}\int_{0}^{T}\langle\eta^{X}_{t},\mathbf{s}\mathbf{i}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt+\int_{0}^{T}\Lambda^{X}_{t}\langle\eta^{X}_{t},\mathbf{s}(f^{\mathcal{I}}_{t}-f_{t})\rangle
dt.$

As $\varphi\geq 0$ and letting
$C=n_{\max}(\lambda_{H}n_{\max}+\lambda_{W}n_{\max}+\beta_{G})$, it follows
that

$p^{X}_{t}(n,n-1)\langle\tilde{\eta}^{X,(n,n-1)}_{T},\varphi\rangle\leq\langle\eta^{X}_{0},f_{0}\rangle+C\int_{0}^{T}\langle\eta^{X}_{t},f^{\mathcal{I}}_{t}\rangle
dt.$ (C4)

Notice that Equation (5) with a change of variables $z=\sigma-T$ implies that

$\langle\eta^{X}_{0},f_{0}\rangle=n\pi^{X}_{n}(1-\varepsilon)^{n-1}\varepsilon\int_{\mathbb{R}}\varphi(z)\mathbf{1}_{\left\\{z>-T\right\\}}\gamma\mathrm{e}^{-\gamma(z+T)}dz.$

Furthermore, by definition, for any $x\in E$,

$f^{\mathcal{I}}_{t}(x)=\mathbf{1}_{\left\\{\mathbf{n}(x)=n,\mathbf{s}(x)=n\right\\}}\int_{0}^{+\infty}\varphi(\sigma-(T-t))\gamma\mathrm{e}^{-\gamma\sigma}d\sigma.$

Thus Equation (C3), Fubini’s theorem and a change of variables
$z=\sigma-(T-t)$ lead to

$\int_{0}^{T}\langle\eta^{X}_{t},f^{\mathcal{I}}_{t}\rangle
dt=\int_{\mathbb{R}}\varphi(z)\int_{0}^{T}p^{X}_{t}(n,n)\gamma\mathrm{e}^{-\gamma(z+(T-t))}dtdx.$

Let

$\displaystyle H^{X,(n,n-1)}_{T}:z\in\mathbb{R}\mapsto$
$\displaystyle\left(n\pi^{X}_{n}(1-\varepsilon)^{n-1}\varepsilon\mathbf{1}_{\left\\{z>-T\right\\}}\gamma\mathrm{e}^{-\gamma(z+T)}\right)$
$\displaystyle+\int_{0}^{T}p^{X}_{t}(n,n)\gamma\mathrm{e}^{-\gamma(z+(T-t))}dt.$

Then $H^{X,(n,n-1)}_{T}$ is integrable on $\mathbb{R}$ with respect to the
Lebesgue measure, as its integral is equal to the right-hand side of Equation
(C4) which is finite. It further satisfies

$p^{X}_{t}(n,s)\langle\tilde{\eta}^{X,(n,n-1)}_{T},\varphi\rangle\leq\int_{\mathbb{R}}H^{X,(n,n-1)}_{T}(z)\varphi(z)dz.$

In particular, let $B$ be a Borel set which is null for the Lebesgue measure.
Consider a mollifier $\psi$ on $\mathbb{R}$ in the same sense as for the Proof
of Lemma 3.1. For $\alpha>0$, let
$\psi_{\alpha}:z\in\mathbb{R}\mapsto\alpha^{-1}\psi(\alpha z)$. Then for any
$\alpha>0$, we may define $\varphi_{\alpha}=\mathbf{1}_{B}*\psi_{\alpha}$,
which is an element of $\mathcal{C}^{1}(\mathbb{R})$ with compact support.
Thus for every $\alpha>0$,

$p^{X}_{T}(n,s)\langle\tilde{\eta}^{X,(n,n-1)}_{T},\varphi_{\alpha}\rangle\leq\int_{\mathbb{R}}H^{X,(n,n-1)}_{T}(z)\varphi_{\alpha}(z)dz.$
(C5)

As $\varphi_{\alpha}$ further is bounded by $1$ for every $\alpha>0$, and as
$\tilde{\eta}^{X,(n,n-1)}_{T}$ is a finite measure and $H^{X,(n,n-1)}_{T}$ is
integrable on $\mathbb{R}$ with respect to the Lebesgue measure, dominated
convergence allows to let $\alpha$ go to $0$ on both sides of inequality (C5).
As $B$ is a null set for the Lebesgue measure, the right-hand side goes to
$0$. This demonstrates that $\tilde{\eta}^{X,(n,n-1)}_{T}$ is absolutely
continuous with respect to the Lebesgue measure on $\mathbb{R}$.

Let us now suppose that $\tilde{\eta}^{X,(n,s+1)}_{T}$ is absolutely
continuous with respect to the Lebesgue measure on $\mathbb{R}^{n-s-1}$ for
some $n^{\prime}\in\llbracket 1,n_{\max}\rrbracket$ and some $s\in\llbracket
0,n-2\rrbracket$. This time, let
$\varphi\in\mathcal{C}^{1}_{b}(\mathbb{R}^{n-s})$ be non-negative and define
$\varphi^{*}\in\mathcal{C}^{1}_{b}(E,\mathbb{R}_{+})$ by

$\forall x\in
E,\quad\varphi^{*}(x)=\mathbf{1}_{\left\\{\mathbf{n}(x)=n,\mathbf{s}(x)=s\right\\}}\varphi(\tau_{1}(x),\dots,\tau_{n-s}(x)).$

For any $T\geq 0$ and $t\in[0,T]$, let $f_{t}:x\in
E\to\varphi^{*}(\Psi(x,T,t))$. Proceeding like before, we obtain that

$p^{X}_{T}(n,s)\langle\tilde{\eta}^{X,(n,s)}_{T},\varphi\rangle\leq\langle\eta^{X}_{0},f_{0}\rangle+C\int_{0}^{T}\langle\eta^{X}_{t},f^{\mathcal{I}}_{t}\rangle
dt.$ (C6)

In order to simplify notations in the following, define for $T\geq t\geq 0$
and $(\tau_{1},\dots,\tau_{n-s})\in\mathbb{R}^{n-s}$ the function

$\varphi_{T,t}(\tau_{1},\dots,\tau_{n-s})=\varphi(\tau_{1}-(T-t),\dots,\tau_{n-s}-(T,t)).$

Then similarly to the case $(n,n-1)$, Equation (5) leads to

$\langle\eta^{X}_{0},f_{0}\rangle=\pi^{X}_{n}\binom{n}{s}(1-\varepsilon)^{s}\varepsilon^{n-s}\int_{\mathbb{R}^{n-s}}\varphi_{T,0}(\sigma_{1},\dots,\sigma_{n-s})\prod_{j=1}^{n-s}(\mathbf{1}_{\left\\{\sigma_{j}>0\right\\}}\gamma\mathrm{e}^{-\gamma\sigma_{j}})d\sigma_{1}\dots
d\sigma_{n-s}.$ (C7)

Notice here that for any $T>0$ and $t\in[0,T]$, the application
$z\in\mathbb{R}^{n-s}\mapsto z-\sum_{k=1}^{n-s}(T-t)e_{k}$, where $e_{k}$ is
the $k$-th vector of the canonical basis of $\mathbb{R}^{n-s}$, defines a
$\mathcal{C}^{1}$-diffeomorphism from $\mathbb{R}^{n-s}$ into itself, whose
Jacobian matrix is the identity matrix of dimension $n-s$. Using the case
$t=0$, Equation (C7) becomes by change of variables:

$\small\langle\eta^{X}_{0},f_{0}\rangle=\pi^{X}_{n}\binom{n}{s}(1-\varepsilon)^{s}\varepsilon^{n-s}\int_{\mathbb{R}^{n-s}}\varphi(z_{1},\dots,z_{n-s})\prod_{j=1}^{n-s}(\mathbf{1}_{\left\\{z_{j}>-T\right\\}}\gamma\mathrm{e}^{-\gamma(z_{j}+T)})dz_{1}\dots
dz_{n-s}.$

Moreover, notice that this time, for any $x\in E$ and $\sigma>0$,

$f^{\mathcal{I}}_{t}(x)=\mathbf{1}_{\left\\{\mathbf{n}(x)=n,\mathbf{s}(x)-1=s\right\\}}\varphi_{T,t}(\tau_{1}(x),\dots,\tau_{n-s-1}(x),\sigma).$

Thus, for $t\in[0,T]$,

$\displaystyle\langle\eta^{X}_{t},f^{\mathcal{I}}_{t}\rangle=\int_{\mathbb{R}^{n-s-1}}\Bigg{(}\int_{0}^{+\infty}\varphi_{T,t}(\tau_{1},\dots,\tau_{n-s-1},\sigma)\gamma\mathrm{e}^{-\gamma\sigma}d\sigma\Bigg{)}p^{X}_{t}(n,s+1)\tilde{\eta}^{X,(n^{\prime},s^{\prime}+1)}_{t}(d\tau_{1},\dots,d\tau_{n^{\prime}-s^{\prime}-1}).$

Fubini’s theorem together with the induction hypothesis yield

$\displaystyle\langle\eta^{X}_{t},f^{\mathcal{I}}_{t}\rangle=\int_{\mathbb{R}_{+}}\int_{\mathbb{R}^{n-s-1}}\varphi_{T,t}(\tau_{1},\dots,\tau_{n-s-1},\sigma)p^{X}_{t}(n,s+1)\tilde{\rho}^{X,(n,s+1)}_{t}(\tau_{1},\dots,\tau_{n-s-1})d\tau_{1}\dots
d\tau_{n-s-1}\gamma\mathrm{e}^{-\gamma\sigma}d\sigma.$

The previously introduced family of $\mathcal{C}^{1}$-diffeomorphisms may
serve again for a change of variables, allowing to obtain that

$\displaystyle\langle\eta^{X}_{t},f^{\mathcal{I}}_{t}\rangle=p^{X}_{t}(n,s+1)\int_{\mathbb{R}^{n-s}}\Bigg{(}\varphi(z_{1},\dots,z_{n-s-1},z_{n-s})\mathbf{1}_{\left\\{z_{n-s}>-(T-t)\right\\}}\gamma\mathrm{e}^{-\gamma(z_{n-s}+(T-t))}$
$\displaystyle\times\;\;\rho^{X,n,s+1}_{t}(z_{1}+(T-t),\dots,z_{n-s-1}+(T-t))\Bigg{)}dz_{1}\dots
dz_{n-s}.$

For $z=(z_{1},\dots,z_{n-s})\in\mathbb{R}^{n-s}$, define

$\displaystyle
H^{X,(n,s)}_{T}(z)=\pi^{X}_{n}\binom{n}{s}(1-\varepsilon)^{s}\varepsilon^{n-s}\prod_{j=1}^{n-s}(\mathbf{1}_{\left\\{z_{j}>-T\right\\}}\gamma\mathrm{e}^{-\gamma(z_{j}+T)})$
$\displaystyle+\int_{0}^{T}\rho^{X,n,s+1}_{t}(z_{1}+(T-t),\dots,z_{n-s-1}+(T-t))\mathbf{1}_{\left\\{z_{n-s}>-(T-t)\right\\}}\gamma\mathrm{e}^{-\gamma(z_{n-s}+(T-t))}dt.$

It then follows like before from Equation (C6) that $H^{X,(n,s)}_{T}$ is
integrable on $\mathbb{R}^{n-s}$ with respect to the Lebesgue measure, and

$p^{X}_{T}(n,s)\langle\tilde{\eta}^{X,(n,s)}_{T},\varphi\rangle\leq\int_{\mathbb{R}^{n-s}}H^{X,(n,s)}_{T}(z_{1},\dots,z_{n-s})dz_{1}\dots
dz_{n-s}.$

The absolute continuity of $\tilde{\eta}^{X,(n,s)}_{T}$ with respect to the
Lebesgue measure on $\mathbb{R}^{n-s}$ follows using the same arguments as
previously. This concludes the proof. ∎