# An integro-ODE-PDE model for Covid-19 pandemic with vaccination and loss of immunity: well-posedness, RBF-FD numerical scheme and simulation

Robin Vaudry, Candy Sonveaux, Didier Georges, Clémentine Prieur

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

## To cite this version:

Robin Vaudry, Candy Sonveaux, Didier Georges, Clémentine Prieur. An integro-ODE-PDE model for Covid-19 pandemic with vaccination and loss of immunity: well-posedness, RBF-FD numerical scheme and simulation. IMA Journal of Mathematical Control and Information, 2025, 42 (3), pp.dnaf020. 10.1093/imamci/dnaf020. hal-05155149

HAL Id: hal-05155149

https://inria.hal.science/hal-05155149v1

Submitted on 9 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# An integro-ODE-PDE model for Covid-19 pandemic with vaccination and loss of immunity: well-posedness, RBF-FD numerical scheme and simulation

Robin Vaudry*, Candy Sonveaux†, Didier Georges‡and Clémentine Prieur§

## Abstract

During an epidemic, the way individuals spread the infection can vary depending on the number of days they have been infected. Vaccination is a key factor in the control of the transmission of the pathogen. We propose a SIR-like model to simulate the propagation of a virus in a closed population. The model presents a continuous structure in age of infectivity and in age of immunity. By using a fixed point theorem, we then prove the well-posedness of the model and properties of its solution. Then we present the radial basis function-generated finite difference (RBF-FD) method we use to numerically solve our problem. Finally, numerical simulations are explored to witness the impact of the parameter values in the model.

**Key words:** SIR-like models, RBF-FD, well-posedness.

## 1 Introduction

The health crisis caused by the COVID-19 pandemic showed once again the need to model the propagation of epidemics. Compartmental models remain the dominant tool for meeting this need. The best-known compartmental model to date is the SIR model introduced by Kermack and McKendrick in 1927 in Kermack and McKendrick (1927). At each time, the population is divided into three categories: S for Susceptible, I for Infectious individuals, and R for Removed individuals who cannot infect nor get infected anymore.

The equations governing the model are a system of ODE-PDE equations:

$$ 
\left\{ 
\begin{aligned}
\frac{dS}{dt}(t) &= -S(t) \int_{0}^{\infty} \beta(\alpha')I(t, \alpha')d\alpha' \\
\frac{\partial I}{\partial t}(t, \alpha) + \frac{\partial I}{\partial \alpha}(t, \alpha) &= -\gamma(\alpha)I(t, \alpha) \\
\frac{dR}{dt}(t) &= \int_{0}^{\infty} \gamma(\alpha')I(t, \alpha')d\alpha' \\
I(t, 0) &= S(t) \int_{0}^{\infty} \beta(\alpha')I(t, \alpha')d\alpha'
\end{aligned}
\right. 
$$

<footnote>*Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000, Grenoble, France</footnote>

<footnote>†Univ. of Namur, Department of Mathematics and naXys, Belgium</footnote>

<footnote>‡Univ. Grenoble Alpes, Grenoble INP, GIPSA-lab, 38000 Grenoble, France</footnote>

<footnote>§Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000, Grenoble, France</footnote>
The variable $\alpha$ here refers to the age of infection of the individuals. In other words, $I(t, \alpha)$ is the proportion of individuals that has been infected for a time $\alpha$. The function $\beta$ is called the transmission rate and depicts how fast the epidemic spreads in the population. The function $\gamma$ is called recovery rate and describes how fast individuals of compartment I recover from the infection.

A variety of epidemiological models has been derived since 1927. In particular, models with an age of infection structure have been studied in Chen et al. (2013), Liu et al. (2015), Chekroun and Kuniya (2020), Magal and McCluskey (2013), McCluskey (2012), Dietz and Schenzle (1985) and more recently in Richard et al. (2021). Other studies focused on the immunity loss process, as in Bhattacharya and Adler (2012), Li et al. (2010), Duan et al. (2014), Barbarossa et al. (2018b). Epidemiologists have also been interested in the impact of vaccination on classical compartmental models, as in Shulgin et al. (1998), Zaman et al. (2008), Chauhan et al. (2014), Okuwa et al. (2021). Lastly, models that combine vaccination and waning immunity have also been studied, under the name of immuno-epidemiological models, as in Heffernan and Keeling (2009), Barbarossa and Röst (2015), Barbarossa et al. (2018a), Ehrhardt et al. (2019).

In this paper, we propose a comprehensive model adapted to a regional scale to track individuals who are susceptible, exposed, infected, hospitalized, vaccinated, recovered or died, taking into account the age of infection for infected individuals and the loss of immunity for the recovered and vaccinated ones, in a vaccination strategy context. After presenting our SEIRHVD model, we will focus on two key aspects: the analysis of its well-posedness, and a cost-effective numerical method based on the approximation of the differential operators in ages using the radial basis function-generated finite difference (RBF-FD) method Fornberg and Flyer (2015) for generating COVID-19 pandemic evolution scenarios. Some simulation results from the model will complete the presentation, showing the impact of the parameter selection. At this stage, it is important to emphasise that the scientific approach we are proposing is generic and could be interesting for epidemics other than Covid-19, as well as other fields of application involving complex, potentially non-linear dynamics in a heterogeneous environment.

The paper is organized as follows: in Section 2, we present the SEIRHVD model. Section 3 is devoted to the extensive statement and proof of the well-posedness of the model. In Section 4, the numerical scheme is exposed. Some numerical simulations are presented and analyzed in Section 5. Finally, we draw some conclusions and perspectives in Section 6.

## 2 SEIRHVD model

Our main aim is to include in our model vaccination and loss of immunity for both recovered or vaccinated individuals. The transmission rate and the recovery rate depends on the age of infectivity $\alpha$ as in Kermack and McKendrick (1927), whereas the loss of immunity depends on the age of immunity $\tilde{\alpha}$. We also include a specific compartment for hospitalized individuals, as the hospital maximal capacity is one key factor in pandemic control. The scheme of the compartmental model built in this paper is presented in Figure 1 below. It is composed with 7 compartments described hereafter:
- *S*: Susceptible. *S* represents individuals that are not infectious but can get infected by contact with an infectious individual.

- *E*: Exposed. These individuals are infected but not infectious yet. This is a transitory compartment between *S* and *I*.

- *I*: Infectious. Infectious individuals are infected individuals that can transmit the pathogen to the *S* individuals. They can present symptoms or not.

- *R*: Recovered. This compartment receives individuals that got cured from the disease, are still immune and have not been vaccinated since they last recovered.

- *H*: Hospitalized. This depicts hospitalized individuals that have mild or severe symptoms.

- *D*: Deceased. This compartment represents the individuals that died due to the infection to the pathogen.

- *V*: Vaccinated. In this compartment, individuals have been vaccinated and still hold an immunity against infection.

<figure><img src="image_3.png" /><figcaption>Figure 1: Compartmental model developed</figcaption></figure>

Our model is still driven by an ODE-PDE system of coupled equations: for
$T > 0, t \in ]0, T], \alpha > 0, \tilde{\alpha} > 0,$

$$
\left\{
\begin{array}{l}
\displaystyle \frac{dS}{dt}(t) = -S(t) \left( \int_{0}^{\infty} \beta(t, \alpha') I(t, \alpha') d\alpha' \right) + \left( \int_{0}^{\infty} (\sigma_{R}(\tilde{\alpha}') R(t, \tilde{\alpha}') + \sigma_{V}(\tilde{\alpha}') V(t, \tilde{\alpha}')) d\tilde{\alpha}' \right) - \theta(t) S(t) \\
\\
\displaystyle \frac{dE}{dt}(t) = S(t) \left( \int_{0}^{\infty} \beta(t, \alpha') I(t, \alpha') d\alpha' \right) - \delta E(t) \\
\\
\displaystyle \frac{dH}{dt}(t) = \int_{0}^{\infty} \gamma_{I}(\alpha') \tau_{H} I(t, \alpha') d\alpha' - \gamma_{H} H(t) \\
\\
\displaystyle \frac{\partial I}{\partial t}(t, \alpha) + \frac{\partial I}{\partial \alpha}(t, \alpha) = -\gamma_{I}(\alpha) I(t, \alpha) - \theta(t) \tau_{A} I(t, \alpha) \\
\\
\displaystyle \frac{\partial R}{\partial t}(t, \tilde{\alpha}) + \frac{\partial R}{\partial \tilde{\alpha}}(t, \tilde{\alpha}) = -\theta_{R}(t, \tilde{\alpha}) R(t, \tilde{\alpha}) - \sigma_{R}(\tilde{\alpha}) R(t, \tilde{\alpha}) \\
\\
\displaystyle \frac{\partial V}{\partial t}(t, \tilde{\alpha}) + \frac{\partial V}{\partial \tilde{\alpha}}(t, \tilde{\alpha}) = -\sigma_{V}(\tilde{\alpha}) V(t, \tilde{\alpha}) \\
\\
\displaystyle \frac{dD}{dt}(t) = \gamma_{H} \tau_{D} H(t) \\
\\
I(t, 0) = \delta E(t) \\
\\
R(t, 0) = \int_{0}^{\infty} (1 - \tau_{H}) \gamma_{I}(\alpha') I(t, \alpha') d\alpha' + (1 - \tau_{D}) \gamma_{H} H(t) \\
\\
V(t, 0) = \int_{0}^{\infty} \theta(t) \tau_{A} I(t, \alpha') d\alpha' + \int_{0}^{\infty} \theta_{R}(t, \alpha') R(t, \tilde{\alpha}') d\tilde{\alpha}' + \theta(t) S(t)
\end{array}
\right.
\qquad (1)
$$

with the constraint that

$$
S(0)+E(0)+H(0)+\int_{0}^{\infty} I(0, \alpha') d \alpha'+\int_{0}^{\infty} R(0, \tilde{\alpha}') d \tilde{\alpha}'+\int_{0}^{\infty} V(0, \tilde{\alpha}') d \tilde{\alpha}'+D(0)=1.
\qquad(2)
$$

We will state later in Section 3.3 (see Theorem 4) that the system described by
equations in 1 preserves mass at all times, that is (2) is not only true at time $t = 0$
but for any $t \ge 0$. For simplicity, we assume in this work that the population
is closed. With this assumption we neglect births and deaths unrelated to the
pathogen, as far as the fact that presymptomatic and asymptomatic cases can
spread the disease between populations through travel. Such an assumption may
be reasonable if the evolution of the system is studied within a short observation
time interval. Let us now provide more details on the parameters involved in
(1):

* $\beta(t, \alpha)$ is the transmission rate of the model, which may depend on time to model, e.g., evolution of sanitary measures;

* $\theta(t)$ represents the vaccination rate of individuals in $S$ and asymptomatic individuals in $I$, $\theta_R(t, \tilde{\alpha})$ represents the vaccination rate of individuals in $R$;

* $\sigma_R(\tilde{\alpha})$ and $\sigma_V(\tilde{\alpha})$ depict the progressive immunity loss in compartments $R$ and $V$;

• $\delta$ is the incubation rate;

* $\gamma_I$ and $\gamma_H$ are the transition rates in $I$ and $H$;

* $\tau_A$ is the proportion of infectious individuals that are asymptomatic;
• $\tau_H$ is the proportion of infected individuals that require hospitalization;

• $\tau_D$ is the proportion of hospitalized individuals that die due to their symptoms.

In the following, we assume that all these parameters are non-negative and bounded by 1. We plot in Figure 1 the state transition diagram completed with parameters involved in each transition between two states in the model.

<figure><img src="image_4.png" /><figcaption>Figure 2: Compartmental model with its parameters</figcaption></figure>

To conclude this section, note that the number of compartments in our model has been chosen as a compromise between a very sophisticated model, probably closer to the reality, and a sparse model that can still be used for inference and control of the epidemics. The model we propose could benefit from a structure in age of individuals and a geographical structure. E.g., a spatial model with an age structure (age of individuals) would lead to an integro-ODE-PDE model in dimension 6 ($t, \alpha, \tilde{\alpha}, a, x, y$, with $a$ the age of individuals and ($x, y$) the spatial coordinates).

# 3 Well-posedness and bounds of the system

This section is devoted to the statement and proof of the well-posedness of the system described by (1). We first state in Theorem 1 existence and unicity of the solution. Then we state in Theorem 3 that the solution depends continuously on the initial and boundary conditions. Finally, we state in Theorem 4 that, under boundedness and regularity assumptions on the initial and boundary conditions, the solution of (1) is bounded between 0 and 1. The proof of these results is a generalization of the proof of analogous results for the coupled PDE-ODE model in Batool and Bajcinca (2021) to coupled ODE-PIDE (partial integro-differential equation) models. The main ingredient of the proof consists in decoupling the equations of System (1) by replacing the variables intervening in the integral terms by any continuous functions. It leads to a new system of independent non-autonomous ordinary differential equations described in (3), whose solutions are explicitly known. This defines an operator (see Definition 4) whose fixed points are the solutions of the initial coupled model. Note that the study of
well-posedness relies on a good choice of solution space (see Definition 2) for which existence and uniqueness can be obtained. For convenience, the chosen functional space on which the operator is defined requires little regularity ($L^1$), but with explicit bounds. Before stating the results, we introduce some useful notation and definitions.

**Definition 1** (Space of initial and boundary conditions). The space $\Xi$ of initial and boundary conditions is defined by:

$$ \Xi = \left\{ \begin{array}{l} \xi = (S_0, E_0, H_0, f_I, f_R, f_V, h_I, h_R, h_V) \text{ such that } S_0, E_0, H_0 \in [0, 1]; \\ f_I, f_R, f_V \in C^1(\mathbb{R}_+, \mathbb{R}_+); h_I, h_R, h_V \in C^1([0, T], \mathbb{R}_+) \text{ with} \\ h_I(0) = f_I(0), h_R(0) = f_R(0), h_V(0) = f_V(0); \\ S_0 + E_0 + H_0 + \int_0^\infty f_I(\alpha')d\alpha' + \int_0^\infty f_R(\tilde{\alpha}')d\tilde{\alpha}' + f_V(\tilde{\alpha}')d\tilde{\alpha}' \in [0, 1] \end{array} \right\} $$

endowed with the following norm:

$$ \begin{array}{l} \|\xi\|_{\Xi} = |S_0| + |E_0| + |H_0| + \|f_I\|_1 + \|f_R\|_1 + \|f_V\|_1 + \|h_I\|_{\infty} + \|h_R\|_{\infty} + \|h_V\|_{\infty} \\ \text{with } |\cdot| \text{ the absolute value, } \|\cdot\|_1 \text{ the } L^1\text{-norm and } \|\cdot\|_{\infty} \text{ the } L^{\infty}\text{-norm.} \end{array} $$

The notation used for the components of $\xi$ are to be linked with the compartments of the model. $S_0, E_0, H_0$ refer to the initial conditions for $S, E$ and $H$ while $f_I, f_R, f_V$ refer to the initial conditions for $I, R, V$ and $h_I, h_R, h_V$ refer to the boundary conditions on $I, R, V$.

In Definition 2 below, we define the state space for any time $T' \in (0, T]$ as we then prove the well-posedness on $[0, nT']$ for all $n = 1, 2, ..., \frac{T}{T'}$ by induction.

**Definition 2** (State space). For any $T' \in (0, T]$, let $X_{T'} = (L^1([0, T]))^3 \times (L^1([0, T'] \times \mathbb{R}_+))^3$. For any $\xi \in \Xi$, let $Y_{T'}^\xi \subseteq X_{T'}$ defined as the set of functions $(y_1, y_2, y_3, y_4, y_5, y_6) \in X_{T'}$ such that:

- $y_1, y_2, y_3 \in L^1([0, T'])$ with $0 \le y_1(t) \le \bar{S}(t, \xi)$, $0 \le y_2(t) \le \bar{E}(t, \xi)$, $0 \le y_3(t) \le \bar{H}(t, \xi)$ almost everywhere, $\bar{S}$ and $\bar{E}$ real-valued functions defined on $[0, T] \times \Xi$ by

$$ \bar{S}(t, \xi) = S_0 + (\|f_R\|_1 + \|f_V\|_1)t + \frac{\|h_R\|_{\infty} + \|h_V\|_{\infty}}{2}t^2, $$

$$ \begin{aligned} \bar{E}(t, \xi) = & E_0 + S_0 \|f_I\|_1 t + \left( \frac{S_0 \|h_I\|_{\infty}}{2} + \frac{\|f_R\|_1 + \|f_V\|_1}{2} \|f_I\|_1 \right) t^2 \\ & + \left( \frac{\|h_R\|_{\infty} + \|h_V\|_{\infty}}{6} \|f_I\|_1 + \frac{\|f_R\|_1 + \|f_V\|_1}{3} \|h_I\|_{\infty} \right) t^3 \\ & + \frac{\|h_R\|_{\infty} + \|h_V\|_{\infty}}{8} \|h_I\|_{\infty} t^4, \end{aligned} $$

$$ \bar{H}(t, \xi) = H_0 + \|f_I\|_1 t + \frac{\|h_I\|_{\infty}}{2} t^2; $$
- $y_4, y_5, y_6 \in L^1([0, T'] \times \mathbb{R}_+)$ with $0 \le \int_0^\infty y_4(t, \alpha')d\alpha' \le \bar{I}(t, \xi)$, $0 \le \int_0^\infty y_5(t, \tilde{\alpha}')d\tilde{\alpha}' \le \bar{R}(t, \xi)$, $0 \le \int_0^\infty y_6(t, \tilde{\alpha}')d\tilde{\alpha}' \le \bar{V}(t, \xi)$ almost everywhere, $\bar{I}, \bar{R}$ and $\bar{V}$ real-valued functions defined on $[0, T] \times \Xi$ by

$$ \bar{I}(t, \xi) = \|f_I\|_1 + \|h_I\|_{\infty} t, \quad \bar{R}(t, \xi) = \|f_R\|_1 + \|h_R\|_{\infty} t, \quad \bar{V}(t, \xi) = \|f_V\|_1 + \|h_V\|_{\infty} t. $$

For $x = (x_1, x_2, x_3, x_4, x_5, x_6) \in X_{T'}$, the norm of $x$ is defined as $\|x\|_{X_{T'}} = \int_0^{T'} \|x\|_t dt$ with for any $t \in [0, T']$:

$$ \|x\|_t := |x_1(t)| + |x_2(t)| + |x_3(t)| + \int_0^\infty |x_4(t, \alpha')| d\alpha' + \int_0^\infty |x_5(t, \tilde{\alpha}')| d\tilde{\alpha}' + \int_0^\infty |x_6(t, \tilde{\alpha}')| d\tilde{\alpha}'. $$

Then we have the following result whose proof is postponed to Appendix A.

**Lemma 1.** For any $\xi \in \Xi$, for any $T' \in (0, T]$, $(Y_{T'}^\xi, ||\cdot||_{X_{T'}})$ is a non-empty complete metric space.

## 3.1 Uniqueness of solution

Theorem 1 below states the existence and the uniqueness of the solution to (1).

**Theorem 1.** Let $\xi \in \Xi$. There exists a unique solution, $y \in Y_T^\xi$ of the coupled system described by (1) on $[0, T]$.

The rest of the section is devoted to the proof of Theorem 1. Following the proof in Batool and Bajcinca (2021), we introduce a decoupled ODE-PIDE system built from the original coupled system defined in (1).

**Definition 3** (Decoupled system). Let $\bar{y}_1, \bar{y}_2, \bar{y}_3 \in C^0([0, T])$ and $\bar{y}_4, \bar{y}_5, \bar{y}_6 \in C^0([0, T] \times \mathbb{R}_+)$. We define a decoupled ODE-PIDE system as:

$$
\left\{
\begin{aligned}
\frac{dS}{dt}(t) &= -S(t) \left( \int_0^\infty \beta(t, \alpha') \bar{y}_4(t, \alpha') d\alpha' + \theta(t) \right) + \int_0^\infty (\sigma_R(\tilde{\alpha}') \bar{y}_5(t, \tilde{\alpha}') + \sigma_V(\tilde{\alpha}') \bar{y}_6(t, \tilde{\alpha}')) d\tilde{\alpha}' \\
\frac{dE}{dt}(t) &= \bar{y}_1(t) \left( \int_0^\infty \beta(t, \alpha') \bar{y}_4(t, \alpha') d\alpha' \right) - \delta E(t) \\
\frac{dH}{dt}(t) &= \int_0^\infty \gamma_I(\alpha') \tau_H \bar{y}_4(t, \alpha') d\alpha' - \gamma_H H(t) \\
\frac{\partial I}{\partial t}(t, \alpha) + \frac{\partial I}{\partial \alpha}(t, \alpha) &= -\gamma_I(\alpha) I(t, \alpha) - \theta(t) \tau_A I(t, \alpha) \\
\frac{\partial R}{\partial t}(t, \tilde{\alpha}) + \frac{\partial R}{\partial \tilde{\alpha}}(t, \tilde{\alpha}) &= -\theta_R(t, \tilde{\alpha}) R(t, \tilde{\alpha}) - \sigma_R(\tilde{\alpha}) R(t, \tilde{\alpha}) \\
\frac{\partial V}{\partial t}(t, \tilde{\alpha}) + \frac{\partial V}{\partial \tilde{\alpha}}(t, \tilde{\alpha}) &= -\sigma_V(\tilde{\alpha}) V(t, \tilde{\alpha}) \\
I(t, 0) &= \delta \bar{y}_2(t) := k_I(t) \\
R(t, 0) &= \int_0^\infty \gamma_I(\alpha') (1 - \tau_H) \bar{y}_4(t, \alpha') d\alpha' + \gamma_H (1 - \tau_D) \bar{y}_3(t) := k_R(t) \\
V(t, 0) &= \int_0^\infty \theta(t) \tau_A \bar{y}_4(t, \alpha') d\alpha' + \int_0^\infty \theta_R(t, \tilde{\alpha}') \bar{y}_5(t, \tilde{\alpha}') d\tilde{\alpha}' + \theta(t) \bar{y}_1(t) := k_V(t)
\end{aligned}
\right.
\quad (3)
$$
and initial conditions:

$$ egin{cases} S(0) = S_0, E(0) = E_0, H(0) = H_0 \\ I(0, \alpha) = f_I(\alpha), R(0, \tilde{\alpha}) = f_R(\tilde{\alpha}), V(0, \tilde{\alpha}) = f_V(\tilde{\alpha}) \end{cases} $$

where $S_0, E_0, H_0 \in [0, 1]$ and $f_I, f_R, f_V \in C^0(\mathbb{R}_+)$ non-negative functions so that $S_0 + E_0 + H_0 + \int_0^\infty f_I(\alpha')d\alpha' + \int_0^\infty f_R(\tilde{\alpha}') + f_V(\tilde{\alpha}')d\tilde{\alpha}' \in [0, 1]$, and with $f_I(0) = k_I(0), f_R(0) = k_R(0), f_V(0) = k_V(0)$.

Note that the structure of this system is similar to the one of the original coupled system described in (1), except that each equation can be solved independently from the others. The decoupled system of equations defined in (3) can be solved easily. The first three equations are linear ODEs of order 1 whose solution is straightforward. The last three ones are linear PDEs of order 1 in time and order 1 in $\alpha$ (or $\tilde{\alpha}$) and can be solved individually using the method of characteristics. Finally, solving these equations we get that the unique continuous solution of the decoupled system described in (3) is given by $Sol^\xi(\tilde{y})$, where $Sol^\xi$ is the solution operator on the state space $Y_T^\xi$ (see Definition 4 below).

**Definition 4 (Solution operator).** For $\xi = (S_0, E_0, H_0, f_I, f_R, f_V, h_I, h_R, h_V) \in \Xi$ and any $T' \in (0, T]$, we define the solution operator $Sol^\xi: Y_{T'}^\xi \to X_{T'}$ as $Sol^\xi(y) = (S^\xi[y], E^\xi[y], H^\xi[y], I^\xi[y], R^\xi[y], V^\xi[y])$ with, for almost every $t \in [0, T']$, $\alpha \in \mathbb{R}_+, \tilde{\alpha} \in \mathbb{R}_+$,

$$ S^\xi[y](t) = \exp\left(-\int_0^t F(\tau)d\tau\right)\left(S_0 + \int_0^t \exp\left(\int_0^\tau F(\tau')d\tau'\right)G(\tau)d\tau\right), $$

$$ E^\xi[y](t) = e^{-\delta t}\left(\int_0^t e^{\delta\tau}y_1(\tau)\left(\int_0^\infty \beta(\tau, \alpha')y_4(\tau, \alpha')d\alpha'\right)d\tau + E_0\right), $$

$$ H^\xi[y](t) = e^{-\gamma_H t}\left(\int_0^t e^{\gamma_H t}\left(\int_0^\infty \gamma_I(\alpha')\tau_H y_4(\tau, \alpha')d\alpha'\right)d\tau + H_0\right), $$

$$ I^\xi[y](t, \alpha) = \begin{cases} f_I(\alpha - t) \exp\left(-\int_0^t \gamma_I(\tilde{u} + \alpha - t) + \tau_A \theta(\tilde{u}) d\tilde{u}\right) & \text{if } t \le \alpha, \\ h_I(t - \alpha) \exp\left(-\int_0^\alpha \gamma_I(\tilde{u}) + \tau_A \theta(\tilde{u} + t - \alpha) d\tilde{u}\right) & \text{if } t > \alpha, \end{cases} $$

$$ R^\xi[y](t, \tilde{\alpha}) = \begin{cases} f_R(\tilde{\alpha} - t) \exp\left(-\int_0^t \theta_R(\tilde{u} + \tilde{\alpha} - t, \tilde{u}) + \sigma_R(\tilde{u} + \tilde{\alpha} - t) d\tilde{u}\right) & \text{if } t \le \tilde{\alpha}, \\ h_R(t - \tilde{\alpha}) \exp\left(-\int_0^{\tilde{\alpha}} \theta_R(\tilde{u}, \tilde{u} + t - \tilde{\alpha}) + \sigma_R(\tilde{u}) d\tilde{u}\right) & \text{if } t > \tilde{\alpha}, \end{cases} $$

$$ V^\xi[y](t, \tilde{\alpha}) = \begin{cases} f_V(\tilde{\alpha} - t) \exp\left(-\int_0^t \sigma_V(\tilde{u} + \tilde{\alpha} - t) d\tilde{u}\right) & \text{if } t \le \tilde{\alpha}, \\ h_V(t - \tilde{\alpha}) \exp\left(-\int_0^{\tilde{\alpha}} \sigma_V(\tilde{u}) d\tilde{u}\right) & \text{if } t > \tilde{\alpha}, \end{cases} $$

$$ \text{where} \quad \left\{ \begin{array}{l} F: \mathbb{R}^+ \to \mathbb{R}^+ \\ t \mapsto \theta(t) + \int_0^\infty \beta(t, \alpha')y_4(t, \alpha')d\alpha' \quad \text{and} \quad \left\{ \begin{array}{l} G: \mathbb{R}^+ \to \mathbb{R}^+ \\ t \mapsto \int_0^\infty (\sigma_R(\tilde{\alpha}')y_5(t, \tilde{\alpha}') + \sigma_V(\tilde{\alpha}')y_6(t, \tilde{\alpha}'))d\tilde{\alpha}'. \end{array} \right. \end{array} \right. $$

Let us now state Proposition 2 which provides the relation between the unique continuous solution to the coupled system described in (1) with the existence of a fixed point for $Sol^\xi$. We first recall that for $y = (y_1, y_2, y_3, y_4, y_5, y_6) \in$
$Y_{T'}^{\xi}$, if there exists $z \in Y_{T'}^{\xi}$ such that $z - y = 0$ almost everywhere and $z$ is a continuous function, then $z$ is called a continuous representative of $y$ and is denoted by $\bar{y}$.

**Proposition 2.** Let $T' \in (0, T]$ and $y = (y_1, y_2, y_3, y_4, y_5, y_6) \in Y_{T'}^{\xi}$ with a continuous representative $\bar{y} = (\bar{y}_1, \bar{y}_2, \bar{y}_3, \bar{y}_4, \bar{y}_5, \bar{y}_6)$. Then $y$ is a fixed point of $\text{Sol}^{\xi}$ on $Y_{T'}^{\xi}$ if and only if $\bar{y}$ solves the coupled system (1) on $[0, T']$.

*Proof of Proposition 2.* The proof is decomposed in two steps.

- First, let $y$ a fixed point of $\text{Sol}^{\xi}$ then $\text{Sol}^{\xi}[y] = y$ almost everywhere. Therefore, the continuous representative of $\text{Sol}^{\xi}[y]$, denoted by $\text{Sol}^{\xi}[y]$, is equal to $\bar{y}$. Now, from the expression of $\text{Sol}^{\xi}[y]$, we have that $\text{Sol}^{\xi}[\bar{y}] = \overline{\text{Sol}^{\xi}[y]}$. Therefore $\bar{y} = \text{Sol}^{\xi}[\bar{y}]$, thus $\bar{y}$ is the unique continuous solution of the decoupled system defined by (3). Therefore, $\bar{y}$ is also solution to the coupled system defined by (1).

- Let now $\bar{y}$ solves the coupled system, then in particular $\bar{y}$ solves the decoupled system. This solution is the unique solution of the decoupled system and is thus equal to $\text{Sol}^{\xi}[\bar{y}]$ which is, as already mentioned in the first step, equal to the continuous representative of $\text{Sol}^{\xi}[y]$. As $\bar{y}$ is also the continuous representative of $y$, we get $\text{Sol}^{\xi}[y] = y$ almost everywhere, that is $y$ is a fixed point of $\text{Sol}^{\xi}$. $\square$

Let $T' \in (0, T]$ and $\xi \in \Xi$. In order to prove Theorem 3, we first state in Lemmas 2 and 3 that the solution operator $\text{Sol}^{\xi}: Y_{T'}^{\xi} \to Y_{T'}^{\xi}$ is Lipschitz continuous, with a Lipschitz constant tending to zero with $T'$.

**Lemma 2.** Let $T' \in (0, T]$ and $\xi \in \Xi$. Then, for any $y \in Y_{T'}^{\xi}$, $\text{Sol}^{\xi}(Y_{T'}^{\xi}) \subseteq Y_{T'}^{\xi}$.

**Lemma 3.** Let $\xi \in \Xi$ and $T' \in [0, T]$. There exists a non-negative constant $C_{\text{Sol}}(\xi)$ such that for any $y_1, y_2 \in Y_{T'}^{\xi}$ and for almost every $t \in [0, T']$,

$$||\text{Sol}^{\xi}[y_1] - \text{Sol}^{\xi}[y_2]||_{X_{T'}} \le T' C_{\text{Sol}}(\xi) ||y_1 - y_2||_{X_{T'}}.$$

The proof of Lemma 2 is postponed to Appendix C, while the one of Lemma 3 is postponed to Appendix B.

*Sketch of proof for Theorem 1.* We only provide here a sketch of the proof as the complete proof is postponed to Appendix B. First, for any $0 < T' \le T$, we consider the solution operator $\text{Sol}^{\xi}: Y_{T'}^{\xi} \to Y_{T'}^{\xi}$. We know from Proposition 2 that a fixed-point of this operator is solution of the coupled system defined by (1) on $[0, T']$. From Lemma 3, there exists $0 < T' \le T$ such that the solution operator $\text{Sol}^{\xi}: Y_{T'}^{\xi} \to Y_{T'}^{\xi}$ is a contraction. Therefore, by Banach fixed-point theorem, it has a unique fixed point that is the solution of the coupled system defined by (1) on $[0, T']$. Note that it is always possible to choose $T'$ such that $\frac{T}{T'} = N \in \mathbb{N}^*$. Then we prove by induction the existence and unicity of a solution on each interval $[0, nT']$, for $n \le N$. Finally, for $n = N$, Theorem 1 is proved. $\square$
## 3.2 Continuity of the solution

This section is devoted to the statement of the continuity of the solution to (1) with respect to the initial and boundary conditions.

**Theorem 3.** Let $\xi^a \in \Xi$ and $y^a$ the corresponding solution to (1) on $[0, T]$. Then, for any $\eta > 0$, there exists a constant $C(\xi^a, \eta)$ such that for any $\xi^b$ with $y^b$ the corresponding solution to (1), $||\xi^a - \xi^b||_\Xi < \eta$ entails:

$$ ||y^a - y^b||_{X_T} \leq C(\xi^a, \eta) ||\xi^a - \xi^b||_\Xi. $$

The proof of Theorem 3 is postponed to Appendix D. It is obtained as a byproduct of the proof of Theorem 1.

## 3.3 Biological properties of the solution

This section is devoted to the statement that the solution to (1) has a biological meaning, in the sense that its components are non-negative and bounded by 1 (since they represent proportions of a population).

**Theorem 4.** Assume that $S_0, E_0, H_0, D_0 \in [0, 1]$, $f_I, f_R, f_V \in C^1(\mathbb{R}_+, \mathbb{R}_+)$, $h_I, h_R, h_V \in C^1([0, T], \mathbb{R}_+)$ with

$$ S_0 + E_0 + H_0 + D_0 + \int_0^\infty f_I(\alpha')d\alpha' + \int_0^\infty f_R(\tilde{\alpha}') + f_V(\tilde{\alpha}')d\tilde{\alpha}' = 1 \quad (4) $$

$$ \text{and} \quad 0 \le f_I, f_R, f_V, h_I, h_R, h_V \le 1. \quad (5) $$

Assume that there exists $\Psi_1, \Psi_2$ two integrable functions on $\mathbb{R}_+$ such that, for any $t \in [0, T]$, $\alpha, \tilde{\alpha} \ge T$,

$$ |f_I(\alpha - t)|, |f_I'(\alpha - t)| \le \Psi_1(\alpha), \quad (6) $$

$$ |f_R(\tilde{\alpha} - t)|, |f_R'(\tilde{\alpha} - t)|, |f_V(\tilde{\alpha} - t)|, |f_V'(\tilde{\alpha} - t)| \le \Psi_2(\tilde{\alpha}). \quad (7) $$

Let $(S, E, H, I, R, V, D)$ be the unique solution to the coupled system described by (1) with initial and boundary conditions $\xi = (S_0, E_0, H_0, f_I, f_R, f_V, h_I, h_R, h_V)$. Then, for any $t \in [0, T]$,

$$ S(t)+E(t)+H(t)+\int_{0}^{\infty} I(t, \alpha')d\alpha'+\int_{0}^{\infty} R(t, \tilde{\alpha}')d\tilde{\alpha}'+\int_{0}^{\infty} V(t, \tilde{\alpha}')d\tilde{\alpha}'+D(t) = 1, \quad (8) $$

$$ 0 \le S(t), E(t), H(t), \int_{0}^{\infty} I(t, \alpha')d\alpha', \int_{0}^{\infty} R(t, \tilde{\alpha}')d\tilde{\alpha}', \int_{0}^{\infty} V(t, \tilde{\alpha}')d\tilde{\alpha}', D(t) \le 1, \quad (9) $$

$$ \text{and for any } \alpha, \tilde{\alpha} \ge 0, \ 0 \le I(t, \alpha), R(t, \tilde{\alpha}), V(t, \tilde{\alpha}) \le 1. \quad (10) $$

We will now discuss the initial and boundary conditions required for Theorem 4. Conditions (4) and (5) ensure that, at the initial time, the sum of the proportions of all health states is equal to one. Conditions (6) and (7) are more technical. They ensure mass conservation over time. The complete proof of Theorem 4 is given in Appendix E. We only provide here a sketch of the proof.
Sketch of the proof for Theorem 4. We note first that for $t = 0$, Equation (8) reduces to Equation (4). Then, we prove that the derivative with respect to time of the sum of the proportions in all the compartments is equal to 0. This is obtained through the equations in (1) and the ability to permute the derivative with respect to time and the integral over the age operators. This second point is proved in Lemma 6 stated and proven in Appendix E. □

# 4 Numerical scheme

In this section, we present the numerical scheme we propose to solve the coupled system described by (1). It consists in a method of lines based on a radial basis function-generated finite difference (RBF-FD) method to compute the differential operators in age $\alpha$ and $\tilde{\alpha}$. The RBF-FD method is an original approach to solve PDEs numerically. It is presented and developed in Fornberg and Flyer (2015), Flyer et al. (2016). This approach is particularly interesting because it enables high precision to be achieved even at low resolution. It is a meshless finite difference method that generates sparse differentiation matrices. In contrast to other traditional methods, such as the finite difference method, the finite element method or the method of characteristics, the RBF-FD method is easy to generalize to dimensions significantly greater than 1, making it easy to introduce additional spatial variables and operators in our model (such as age or geographical location of individuals in a domain with complex boundaries). The obtained system of nonlinear differential equations is then discretized in time using a partially implicit Euler scheme of order 1. Section 4.1 is devoted to the presentation of the RBF-FD method. Section 4.2 focuses on the computation of the integrals that appear in the equations of the model. Finally we provide in Section 4.3 the equations driving the numerical scheme we propose for our model.

## 4.1 RBF-FD method

In this section, we explain how we discretize the PDEs in (1). The main reason we use RBF-FD approximation here is that we have in mind future generalization to spatialized epidemiological models. Indeed, RBF-FD approximations generalize classical grid-based finite differences, replacing mesh generation by scattered node generation, increasing by the way the geometric flexibility of the discretizations and making it easier to carry out local refinement in critical areas. In order to illustrate the approach, we focus on the hyperbolic equation on $I$ given by:

$$ \left\{ \begin{array}{l} \frac{\partial I}{\partial t}(t, \alpha) + \frac{\partial I}{\partial \alpha}(t, \alpha) = -\gamma_I(\alpha)I(t, \alpha) - \theta(t)\tau_A I(t, \alpha), \\ I(t, 0) = \delta E(t) \end{array} \right. $$

for $t \in [0, T]$ and $\alpha \ge 0$. The same method can be applied for the equations on $R$ and $V$. In our study, the differentiation operator in time $\partial/\partial t$ is approximated by an Euler scheme. So we only discuss in this section the approximation of the operator $\partial/\partial \alpha$ as well as the Dirichlet boundary condition on the edge of the domain.
First, we introduce a truncation age $\alpha_{max} > 0$. We then propose a discretization of the age of disease $\alpha$ on the bounded interval $J_b = [0, \alpha_{max}]$. We will discuss later the choice of the sequence of $N_\alpha + 1$ discretized ages $\{\alpha_i\}_{i=0}^{N_\alpha}$. The method consists in approximating the linear operator $\partial/\partial\alpha$ at a point $\alpha_*$ in the domain $J_b$ by a linear combination of the function values, $\{f_k\}_{k=1}^n$ at the nearest $n << N_\alpha$ node locations $\{\alpha_\ell\}_{\ell \in \Lambda(\alpha_*)}$; this group of nodes composes the stencil, denoted by $S_n(\alpha_*) = \{\alpha_\ell\}_{\ell \in \Lambda(\alpha_*)}$. In other words, we seek weights $\{w_\ell\}_{\ell \in \Lambda(\alpha_*)}$ such that:

$$L f|_{\alpha=\alpha_*} := \frac{\partial f}{\partial \alpha}|_{\alpha=\alpha_*} \approx \sum_{\ell \in \Lambda(\alpha_*)} w_\ell f_\ell,$$

with $n$ and $w_\ell$ being known as the stencil size and differentiation weights, respectively. To calculate $w_\ell$, we enforce that this linear combination of function values is exact for the interpolant function $\tilde{f}$ of the form

$$\tilde{f}(\alpha) = \sum_{\ell \in \Lambda(\alpha_*)} \lambda_{\ell} \psi(||\alpha - \alpha_{\ell}||) + \sum_{j=0}^{m} \gamma_j \alpha^j, \quad (11)$$

with $\psi: r \in \mathbb{R}_+ \mapsto \phi(r) = r^3$ the polyharmonic spline and with the constraints $\sum_{\ell \in \Lambda(\alpha_*)} \lambda_\ell \alpha_\ell^j = 0$, for $j = 0, \dots, m$. Note that polynomials help in improving accuracy Flyer et al. (2016). By interpolant we mean that $\tilde{f}(\alpha_\ell) = f(\alpha_\ell) = f_\ell$ for $\ell \in \Lambda(\alpha_*)$. For sake of clarity, we assume (without loss of generality) in the following that $\Lambda(\alpha_*) = \{1, \dots, n\}$. The weights are obtained by solving the following linear system (see, e.g., Fornberg and Flyer (2015)):

$$H^T \begin{pmatrix} w_1 \\ \vdots \\ w_n \\ w_{n+1} \\ w_{n+2} \\ \vdots \\ w_{n+m} \end{pmatrix} := \begin{pmatrix} A & \begin{pmatrix} 1 & \alpha_1 & \cdots & \alpha_1^m \\ \vdots & \vdots & \ddots & \vdots \\ 1 & \alpha_n & \cdots & \alpha_n^m \end{pmatrix} \\ \begin{pmatrix} 1 & \cdots & 1 \\ \alpha_1 & \cdots & \alpha_n \\ \vdots & & \vdots \\ \alpha_1^m & \cdots & \alpha_n^m \end{pmatrix} & \begin{pmatrix} w_1 \\ \vdots \\ w_n \\ w_{n+1} \\ w_{n+2} \\ \vdots \\ w_{n+m} \end{pmatrix} \end{pmatrix} \begin{pmatrix} L\psi(||\alpha - \alpha_1||)|_{\alpha=\alpha_*} \\ \vdots \\ L\psi(||\alpha - \alpha_n||)|_{\alpha=\alpha_*} \\ L1|_{\alpha=\alpha_*} \\ L\alpha|_{\alpha=\alpha_*} \\ \vdots \\ L\alpha^m|_{\alpha=\alpha_*} \end{pmatrix}$$

with

$$A = \begin{pmatrix} \psi(||\alpha_1 - \alpha_1||) & \psi(||\alpha_1 - \alpha_2||) & \cdots & \psi(||\alpha_1 - \alpha_n||) \\ \psi(||\alpha_2 - \alpha_1||) & \psi(||\alpha_2 - \alpha_2||) & \cdots & \psi(||\alpha_2 - \alpha_n||) \\ \vdots & \vdots & \ddots & \vdots \\ \psi(||\alpha_n - \alpha_1||) & \psi(||\alpha_n - \alpha_2||) & \cdots & \psi(||\alpha_n - \alpha_n||) \end{pmatrix}.$$
This system has been proven to be invertible, that is matrix H invertible, as soon as $m + 1 \le n$ (see, e.g., De Marchi (2013)). In the following, we choose $f(\cdot) := I(t, \cdot)$ and $m = 0$. We then get the following approximation:

$$ \begin{aligned} \frac{\partial I}{\partial \alpha}(t, \alpha_*) &= (L\psi(\|\alpha - \alpha_1\|)|_{\alpha=\alpha_*}, \dots, L\psi(\|\alpha - \alpha_n\|)|_{\alpha=\alpha_*}, L1|_{\alpha=\alpha_*}) H^{-1}(f_1, \dots, f_n, 0)^T \\ &= (w_1, \dots, w_n) (I(t, \alpha_1), \dots, I(t, \alpha_n))^T. \end{aligned} $$

By concatenating the above system for $\alpha_* = \alpha_i$, $0 \le i \le N_\alpha$, we obtain a $(N_\alpha + 1) \times (N_\alpha + 1)$ matrix $M_\alpha$ such that:

$$ \frac{\partial I}{\partial \alpha}(t) = M_{\alpha}I(t) $$

with $I(t) := (I(t, 0), I(t, \alpha_1), \dots, I(t, \alpha_{N_\alpha}))^T :=: (I^{bound}(t), (I^{int}(t))^T)$. Then we introduce matrices $M_{1,\alpha}$ and $M_{2,\alpha}$ defined by:

$$ M_{1,\alpha} = (M_{\alpha}[i,j])_{1\le i,j \le N_{\alpha}}, \quad M_{2,\alpha} = (M_{\alpha}[i,0])_{1\le i \le N_{\alpha}}. $$

The system discretized in age corresponding to the PDE on I can then be written as an ODE:

$$ \left\{ \begin{array}{l} \frac{dI^{int}}{dt}(t) = -M_{1,\alpha}I^{int}(t) - \delta E(t)M_{2,\alpha} + F_I(t, I^{int}(t)) \\ I^{bound}(t) = \delta E(t) \end{array} \right. $$

$$ \text{wih } F_I(t, I^{int}(t)) = \begin{pmatrix} -(\gamma(\alpha_1) + \tau_A \theta(t))I(t, \alpha_1) \\ \vdots \\ -(\gamma(\alpha_{N_\alpha}) + \tau_A \theta(t))I(t, \alpha_{N_\alpha}) \end{pmatrix}. $$

Similarly, we build the matrix $M_{\tilde{\alpha}}$ for the approximation of the differential operator $\frac{\partial}{\partial \tilde{\alpha}}$ on the nodes $(\tilde{\alpha}_i)_{i=0}^{N_{\tilde{\alpha}}}$ and the corresponding submatrices $M_{1,\tilde{\alpha}}$ and $M_{2,\tilde{\alpha}}$.

## 4.2 QMC computation of the integrals

In this section, we explain how the integrals in (1) are computed using a Quasi-Monte Carlo sampling (QMC). In order to illustrate the approach, we focus on the computation of the integral

$$ \int_{0}^{\infty} \beta(t, \alpha') I(t, \alpha') d\alpha'. $$

We first truncate the integral on the bounded interval $J_b = [0, \alpha_{max}]$ as for the discretization of the age operators. For an integral defined on the $d$-dimensional unit hypercube $[0, 1]^d$, QMC numerical integration consists in evaluating the intergand on the $n$ first terms of a low discrepancy sequence of points over
$[0, 1]^d$ and then in approximating the integral by computing the empirical mean of integrand evaluations. For the $n$ first terms of a low discrepancy sequence over $[0, 1]^d$, $\{x_i\}_{i=1}^n$, the well known Koksma-Hlawka inequality provides:

$$ \left| \frac{1}{n} \sum_{i=1}^{n} f(x_i) - \int_{[0,1]^d} f(x) dx \right| \le D_n^* \times V_{HK}(f), $$

where $f: [0, 1]^d \to \mathbb{R}$, $D_n^*$ is the so-called star discrepancy, $V_{HK}(f)$ is the total variation of $f$ in the sense of Hardy and Krause (see, e.g., Kuipers and Niederreiter (2012)). In the following, we only compute one-dimensional integrals (with respect to $\alpha$ or $\tilde{\alpha}$). Integrals are computed with a Sobol' sequence $\{x_i\}_{i=0}^{N_\alpha}$ on $[0, 1]$, for which we know that $D_n^* = \mathcal{O}(\log N_\alpha / N_\alpha)$ (see, e.g., Lemieux (2009)). We then transform this sequence into $\{\alpha_i\}_{i=0}^{N_\alpha} = \{\alpha_{max}x_i\}_{i=0}^{N_\alpha}$ on $[0, \alpha_{max}]$. Finally, QMC approximation of $\int_0^\infty \beta(t, \alpha')I(t, \alpha')d\alpha'$ is given by:

$$ \frac{\alpha_{max}}{N_{\alpha} + 1} \sum_{i=0}^{N_{\alpha}} \beta(t, \alpha_i) I(t, \alpha_i). $$

It is worth noticing that the same sequence $\{\alpha_i\}_{i=0}^{N_\alpha} = \{\alpha_{max}x_i\}_{i=0}^{N_\alpha}$ is used for the RBF-FD approximation of the age operators, which avoids additional changes of coordinates between the variables used for the age operators and the integrals.

## 4.3 Resulting recursive equations

Based on the approximations of the age operators and the integrals, we will now consider two different treatments of the compartimental equations, distinguishing between variables in age-distributed equations ($I$, $R$ and $V$) and those in ordinary differential equations ($S$, $E$, $H$ and $D$).

For the variables structured in age, we choose to implement a semi-explicit Euler scheme. We start from the following system obtained in Section 4.1:

$$ \left\{ \begin{array}{l} \frac{dI^{int}}{dt}(t) = -M_{1,\alpha}I^{int}(t) - \delta E(t)M_{2,\alpha} + F_I(t, I^{int}(t)) \\ I^{bound}(t) = \delta E(t). \end{array} \right. $$

We then discretize time interval $[0, T]$. More precisely, for $N \in \mathbb{N}^*$, we introduce $\Delta t = \frac{T}{N}$ and define, for $0 \le n \le N-1$, $t_n = n \times \Delta t$. We then implement the following scheme:

$$ \left\{ \begin{array}{l} \frac{\hat{I}_{n+1}^{int} - \hat{I}_{n}^{int}}{\Delta t} = -M_{1,\alpha}\hat{I}_{n+1}^{int} - \delta E_n M_{2,\alpha} + F_I(t_n, \hat{I}_{n}^{int}) \\ I_{n+1}^{bound} = \delta E_{n+1} \end{array} \right. \quad (12) $$

where $\hat{I}_{n}^{int}$ is the approximation of $I^{int}(t_n)$, $E_n$ the one of $E(t_n)$ and $I_{n}^{bound}$ the one of $I^{bound}(t_n)$. This method is semi-implicit in the sense that the implication only concerns the approximate age operator, i.e., the term $M_{1,\alpha}\hat{I}_{n+1}^{int}$ is
$S_n$ and to express $S_n$ as $S_n = 1 - E_n - H_n - D_n - \text{Int}_n(I) - \text{Int}_n(R) - \text{Int}_n(V)$. Such an elimination of a compartimental equation in SIR-type system is classical in the literature; see, for instance, Kendal (1956) for an early reference.

It follows that the time discretization of system (1) is given by the following system of recursive equations:

$$
\left\{
\begin{aligned}
E_{n+1} &= E_n + \Delta t \left[ S_n \text{Int}_n(\beta, I) - \delta E_n \right] \\
H_{n+1} &= H_n + \Delta t \left[ \tau_H \text{Int}_n(\gamma_I, I) - \gamma_H H_n \right] \\
D_{n+1} &= D_n + \Delta t \gamma_H \tau_D H_n \\
\hat{I}_{n+1}^{\text{int}} &= [\text{Id} + \Delta t M_{1,\alpha}]^{-1} [\hat{I}_n - \Delta t \delta E_n M_{2,\alpha} + \Delta t F_I(t_n, \hat{I}_n^{\text{int}})] \\
\hat{R}_{n+1}^{\text{int}} &= [\text{Id} + \Delta t M_{1,\hat{\alpha}}]^{-1} \left[ \hat{R}_n^{\text{int}} - \Delta t \left( \left( (1-\tau_H)\text{Int}_n(\gamma_I, I) + \gamma_H(1-\tau_D)H_n \right) M_{2,\hat{\alpha}} + F_R(t_n, \hat{R}_n^{\text{int}}) \right) \right] \\
\hat{V}_{n+1}^{\text{int}} &= [\text{Id} + \Delta t M_{1,\hat{\alpha}}]^{-1} \left[ \hat{V}_n^{\text{int}} - \Delta t \left( \left( \tau_A \text{Int}_n(\theta, I) + \text{Int}_n(\theta_R, R) + \theta(t_n)S_n \right) M_{2,\hat{\alpha}} + F_R(t_n, \hat{R}_n^{\text{int}}) \right) \right] \\
S_{n+1} &= 1 - E_{n+1} - H_{n+1} - D_{n+1} - \text{Int}_{n+1}(I) - \text{Int}_{n+1}(R) - \text{Int}_{n+1}(V)
\end{aligned}
\right.
\quad (13)
$$

where

$$
F_I(t_n, \hat{I}_n^{\text{int}}) = \begin{pmatrix} -(\gamma(\alpha_0) + \tau_A \theta(t)) \hat{I}_n^{\text{int}}(\alpha_0) \\ \vdots \\ -(\gamma(\alpha_{N_\alpha}) + \tau_A \theta(t_n)) \hat{I}_n^{\text{int}}(\alpha_{N_\alpha}) \end{pmatrix},
$$

$$
F_R(t_n, \hat{R}_n^{\text{int}}) = \begin{pmatrix} -(\sigma_R(\tilde{\alpha}_0) + \theta_R(t_n, \tilde{\alpha}_0))R(t_n, \tilde{\alpha}_0) \\ \vdots \\ -(\sigma_R(\tilde{\alpha}_{N_{\tilde{\alpha}}}) + \theta_R(t_n, \tilde{\alpha}_{N_{\tilde{\alpha}}}))R(t_n, \tilde{\alpha}_{N_{\tilde{\alpha}}}) \end{pmatrix},
$$

and

$$
F_R(t_n, \hat{R}_n^{\text{int}}) = \begin{pmatrix} -\sigma_V(\tilde{\alpha}_0)V(t_n, \tilde{\alpha}_0) \\ \vdots \\ -\sigma_V(\tilde{\alpha}_{N_{\tilde{\alpha}}})V(t_n, \tilde{\alpha}_{N_{\tilde{\alpha}}}) \end{pmatrix}.
$$

# 5 Numerical simulations

In this section, we illustrate the effectiveness of the proposed model and numerical scheme with a simulation based on a set of realistic parameters and the population of Isère department in France. Table 1 presents the scalar parameters. Their values are found in the literature. Table 2 displays the functional parameters involved in the description of the model.

$\gamma_I$, $\sigma_R$ and $\sigma_V$ are defined as sigmoid functions with an inflection age corresponding to the average residence time in the compartment (17 days for infection and 180 days for loss of immunity). In our simulations, $\alpha \mapsto \beta(\alpha)$ does not depend on time and corresponds to $\beta_m$ in (Richard et al., 2021, Eq. (9) and Fig.
<table><caption>Table 1: Scalar parameters of the model</caption><thead><tr><td>Variable</td><td>Value</td><td>Reference</td></tr></thead><tbody><tr><td>δ</td><td><sup>1</sup>&frasl;<sub>4.5</sub></td><td>Calba et al. (2021)</td></tr><tr><td>γ<sub>H</sub></td><td><sup>1</sup>&frasl;<sub>14.5</sub></td><td>Boëlle et al. (2020)</td></tr><tr><td>τ<sub>A</sub></td><td><sup>20</sup>&frasl;<sub>100</sub></td><td>Byambasuren et al. (2020)</td></tr><tr><td>τ<sub>H</sub></td><td><sup>9</sup>&frasl;<sub>100</sub></td><td>Oliver et al. (2022)</td></tr><tr><td>τ<sub>D</sub></td><td><sup>8</sup>&frasl;<sub>1000</sub></td><td>Estimated value in accordance with data.gouv.fr</td></tr></tbody></table>

<table><caption>Table 2: Function parameters of the model</caption><thead><tr><td>Variable</td><td>Expression (values are in days)</td></tr></thead><tbody><tr><td>β(t, α)</td><td>see (Richard et al., 2021, Eq. (9) and Fig. 2.A) and Equation (14)</td></tr><tr><td>γ<sub>I</sub>(α)</td><td><sup>1</sup>&frasl;<sub>1 + e<sup>-(α-17)</sup></sub></td></tr><tr><td>σ<sub>R</sub>(α̃)</td><td><sup>1</sup>&frasl;<sub>1 + e<sup>-(α̃-180)</sup></sub></td></tr><tr><td>σ<sub>V</sub>(α̃)</td><td><sup>1</sup>&frasl;<sub>1 + e<sup>-(α̃-180)</sup></sub></td></tr><tr><td>θ(t)</td><td>from DataVaccin-COVID (Isère department, France)</td></tr><tr><td>θ<sub>R</sub>(t, α̃)</td><td>θ(t) × <sup>1</sup>&frasl;<sub>1 + e<sup>-(α̃-21)</sup></sub></td></tr></tbody></table>

2.A):

$$ \beta : \alpha \mapsto \begin{cases} 0.2 \times \frac{3}{5.65} \times \left(\frac{\alpha}{5.65}\right)^2 \times e^{-\left(\frac{\alpha}{6.5}\right)^3} & \text{for } \alpha \le 5.2, \\ 0.2 \times \left(\frac{1}{2}\right)^{(\alpha-5.2)} \times \frac{3}{5.65} \times \left(\frac{\alpha}{5.65}\right)^2 \times e^{-\left(\frac{\alpha}{6.5}\right)^3} & \text{for } \alpha > 5.2. \end{cases} \quad (14) $$

The expression of $\beta$ makes use of the Weibull distribution as proposed in Richard et al. (2021), since it appears to best fit the infection transmission timing. $\beta$ also reflects that individuals do not isolate before symptoms (here when $\alpha \le 5.2$), and that isolation is stronger when symptoms are more severe.

The vaccination inputs $\theta(t)$ and $\theta_R(t, \tilde{\alpha})$ are obtained with vaccination data available on the website DataVaccin-COVID for the Isère Department in France. Figures 3 and 4 show the transmission rate $\beta$ and the vaccination rate $\theta$.

We now present our choice of initial conditions at time $t = t_0 = 0$. Let
<figure><img src="image_5.png" /><figcaption>Figure 3: Disease transmission rate $\beta$ as a function of disease age $\alpha$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Vaccination rate of $S$ and $I$ asymptomatic individuals as a function of time.</figcaption></figure>

$S_0, E_0, H_0, I(t = 0, \tilde{\alpha}), R(t = 0, \tilde{\alpha}), V(t = 0, \tilde{\alpha}), D_0$ be the initial conditions of our model. We wish these quantities:

i) to correspond to an epidemiological reality: for all $\alpha, \tilde{\alpha} > 0$,

$$0 \le S_0, E_0, H_0, I(t = 0, \alpha), R(t = 0, \tilde{\alpha}), V(t = 0, \tilde{\alpha}), D_0 \le 1;$$
ii) to be compatible with boundary conditions:

$$I(t = 0, \alpha = 0) = \delta E_0,$$

$$R(t = 0, \tilde{\alpha} = 0) = \int_{0}^{\infty} (1 - \tau_{H})\gamma_{I}(\alpha')I(t = 0, \alpha')d\alpha' + (1 - \tau_{D})\gamma_{H}H_{0},$$

$$V(t = 0, \tilde{\alpha} = 0) = \int_{0}^{\infty} \theta(0)\tau_{A}I(t = 0, \alpha')d\alpha' + \int_{0}^{\infty} \theta_{R}(0, \alpha')R(t = 0, \tilde{\alpha}')d\tilde{\alpha}' + \theta(0)S_{0};$$

iii) to satisfy the population conservation constraint at $t = 0$.

In the simulation, we choose to work with zero vaccination terms in the initial state, which corresponds to the start of the epidemic. Thus, we have

$$\theta(0) = 0 \quad \text{and} \quad \forall \tilde{\alpha} > 0, \theta_{R}(0, \tilde{\alpha}) = 0.$$

We also choose an initial population such that $H_0 = D_0 = 0$ et $0 \le S_0, E_0 \le 1$. Hence, constraint ii) gives us for $V$:

$$V(t = 0, \tilde{\alpha} = 0) = 0.$$

Next, we select the initial condition of the infectious individuals:

$$I(t = 0, \alpha) = \delta E_0 e^{-\alpha}.$$

This expression satisfies conditions i) and ii). In addition, an exponential decrease in the number of infected individuals with age $\alpha$ is consistent with epidemiological observations.

Condition ii) on $R$ and the sigmoidal expression of $\gamma_I$ then give us:

$$\begin{aligned} R(t = 0, \tilde{\alpha} = 0) &= \delta E_0 (1 - \tau_H) \int_0^\infty \frac{0.5 e^{-\alpha'}}{1 + e^{-0.5(\alpha' - 10)}} d\alpha' \\ &= 0.5 \times \delta E_0 (1 - \tau_H) e^{-10} \int_{-10}^\infty \frac{e^{-u}}{1 + e^{-\varepsilon u}} du, \end{aligned}$$

with

$$e^{-10} \int_{-10}^{\infty} \frac{e^{-u}}{1 + e^{-\varepsilon u}} du \le e^{-10} \int_{-10}^{\infty} e^{-u} du \le e^{-10} e^{10} \le 1.$$

It follows that $0 \le R(t = 0, \tilde{\alpha} = 0) \le 1$.

We then choose for the expressions of $R(t = 0, \tilde{\alpha})$ and $V(t = 0, \tilde{\alpha})$ for all $\tilde{\alpha} \ge 0$:

$$R(t = 0, \tilde{\alpha}) = R(t = 0, \tilde{\alpha} = 0)e^{-\kappa \tilde{\alpha}} \text{ with } \kappa > 0 \text{ and } V(t = 0, \tilde{\alpha}) = 0.$$

These expressions satisfy conditions i) and ii). We now seek to determine whether there exists $\kappa > 0$ such that condition iii) is also verified. We have:

$$\int_{0}^{\infty} I(t = 0, \alpha')d\alpha' = \delta E_0,$$

$$\int_{0}^{\infty} R(t = 0, \tilde{\alpha}')d\tilde{\alpha}' = \frac{1}{\kappa}\delta E_0(1 - \tau_H)\int_{0}^{\infty} \gamma_I(\alpha')e^{-\alpha'}d\alpha',$$

$$\int_{0}^{\infty} V(t = 0, \tilde{\alpha}')d\tilde{\alpha}' = 0.$$
<table><caption>Table 3: Simulation parameters</caption><thead><tr><th>Variables</th><th>Values (time and ages are in days)</th></tr></thead><tbody><tr><td>(N<sub>α</sub>, N<sub>̃α</sub>)</td><td>(512, 512)</td></tr><tr><td>(α<sub>max</sub>, ̃α<sub>max</sub>)</td><td>(50, 400)</td></tr><tr><td>(N, Δt)</td><td>(540, 1)</td></tr><tr><td>(S<sub>0</sub>, E<sub>0</sub>, H<sub>0</sub>, D<sub>0</sub>)</td><td>(0.95, 0.05, 0, 0)</td></tr></tbody></table>

We then have:

$$ 
\begin{aligned}
S_0 + E_0 + H_0 + \int_0^\infty I(t=0, \alpha')d\alpha' + \int_0^\infty R(t=0, \tilde{\alpha}')d\tilde{\alpha}' + \int_0^\infty V(t=0, \tilde{\alpha}')d\tilde{\alpha}' + D_0 \\
&= S_0 + E_0 + \int_0^\infty I(t=0, \alpha')d\alpha' + \int_0^\infty R(t=0, \tilde{\alpha}')d\tilde{\alpha}' \\
&= S_0 + E_0 + \delta E_0 + \frac{1}{\kappa}\delta E_0(1 - \tau_H) \int_0^\infty \gamma_I(\alpha')e^{-\alpha'}d\alpha' \\
&= S_0 + (1 + \delta)E_0 + \frac{1}{\kappa}\delta E_0(1 - \tau_H) \int_0^\infty \gamma_I(\alpha')e^{-\alpha'}d\alpha'.
\end{aligned}
$$

We therefore choose $S_0 \in ]0, 1[$ and $E_0 \in ]0, 1[$ such that

$$S_0 + (\delta + 1)E_0 < 1.$$

To ensure condition iii), we then select $\kappa > 0$ as

$$\kappa = \frac{\delta E_0 (1 - \tau_H) \int_0^\infty \gamma_I(\alpha') e^{-\alpha'} d\alpha'}{1 - S_0 - (\delta + 1)E_0},$$

so that

$$ 
\begin{aligned}
S_0 + E_0 + H_0 + \int_0^\infty I(t=0, \alpha')d\alpha' + \int_0^\infty R(t=0, \tilde{\alpha}')d\tilde{\alpha}' + \int_0^\infty V(t=0, \tilde{\alpha}')d\tilde{\alpha}' + D_0 \\
= S_0 + (1 + \delta)E_0 + \frac{1}{\kappa}\delta E_0(1 - \tau_H) \int_0^\infty \gamma_I(\alpha')e^{-\alpha'}d\alpha' = 1.
\end{aligned}
$$

Table 3 provides parameter values used for the simulation. The simulations are computed on a time horizon of 540 days (18 months). Simulation results are plotted in Figures 5 to 9. In Figures 5 and 6, we have plotted $I : t \mapsto \int_0^\infty I(t, \alpha')d\alpha'$, $R : t \mapsto \int_0^\infty R(t, \tilde{\alpha}')d\tilde{\alpha}'$ and $V : t \mapsto \int_0^\infty V(t, \tilde{\alpha}')d\tilde{\alpha}'$. Figure 6 zooms on the initial transient of state variables $H$, $I$, $R$ and $D$. Figures 8, 7, and 9 focus on the initial transient in time and age of the variables $I(t, \alpha)$, $R(t, \tilde{\alpha})$ and $V(t, \tilde{\alpha})$, respectively.

With our choice of vaccination rates and model parameter values, we observe one and only one infection peak, followed by a non-zero proportion of individuals conserving their immunity over time. The vaccination of individuals in $R$ contributes to the slow decrease of the number of individuals in $R$ from day 40
to day 220, before this number reaches 0. We can check that the simulation satisfies the conservation of the population mass at any time step, as required.

To demonstrate the value of the SEHIRVD model distributed in ages of infectivity and immunity, we simulated its non-distributed version, in the same simulation conditions and with $\beta = 0.073$, $\gamma_I = \frac{1}{17}$, $\sigma_R = \sigma_V = \frac{1}{180}$ and $\theta_R(t) = \frac{1}{21}\theta(t)$, while the other parameters were kept unchanged. The value of $\beta$ has been chosen so that the maximum of the peak of infected individuals is the same in the non-distributed model as in the age-structured model. Here, the values of $\gamma_I$, $\sigma_R$ and $\sigma_V$ have been chosen to ensure consistency between asymptotic solution of the non-distributed model and the integral of the solution of the age-distributed system. Figure 10 shows the dynamics of the non-distributed model. Compared to Figure 5, some notable differences can be observed, in particular a slower flattening of the trajectories of $I(t)$ and $R(t)$ in the non-distributed case.

<figure><img src="image_7.png" /><figcaption>Figure 5: Simulation of the SEHIRDV model distributed in ages of infectivity and immunity. Here, $I, R, V$ denote the integral of the $I(t, \alpha), R(t, \tilde{\alpha}), V(t, \tilde{\alpha})$ in ages, respectively.</figcaption></figure>

It is interesting to note that the simulation in Python (including the computations of RBF-FD operators) takes less than 2 seconds on an Intel core-i7 computer without use of GPU parallelization, demonstrating the cost-effectiveness of the numerical scheme we propose in this paper, despite the complexity of the model.

# 6 Conclusions and perspectives

In this paper, we propose a compartmental model to simulate the spread of the SARS-CoV-2. The model is illustrated in Isère department, France. Our model
<figure><img src="image_8.png" /><figcaption>Figure 6: Simulation of the SEHIRDV model. Only *H*, the integral of *I* in *α*, the integral of *R* in *$\tilde{\alpha}$* et *D* are shown.</figcaption></figure>

is structured in age of infection $\alpha$ and age of immunity $\tilde{\alpha}$. The equations that drive the model from a system of PDE-ODEs with integrals in source terms and boundary conditions. We proved that, under realistic assumptions, the system is well-posed: there is a unique solution to the system of equations and this solution depends continuously on the initial conditions of the problem. To solve numerically these equations, we propose a simulation scheme based on the RBF-FD method, that is an effective alternative to the classical finite difference schemes to discretize linear differential operators. In order to compute the integral terms, we use a Quasi-Monte Carlo quadrature scheme. At this stage, we do not pretend to be able to predict the spread of SARS-CoV-2 in Isère. In fact, the main goal of this work is to explore and propose a methodology to model, analyze and simulate structured compartmental models in epidemiology. It is important to emphasise that the scientific approach we have proposed in this paper is generic and could be interesting for epidemics other than Covid-19, by introducing information about the distribution of the infected and immunized individuals. We are convinced that this approach can be generalized to any epidemic where the pathogen is spread by contact between individuals. In this paper, the model presented is structured in age of disease $\alpha$ and age of immunity $\tilde{\alpha}$. Other structures can also be considered. The age of the individuals $a$ is an interesting factor that can be added to the model in order to display the difference in the interactions between individuals and in the response of the immune system with respect to the age. Moreover, the addition of the spatial variable $x$ depicting the geographical location of individuals in a given area would allow the model to emphasize on the movement of individuals and on where the contacts leading to infection occur, thanks to the introduction of diffusion operators as proposed in Oliver et al. (2022). This work can also be continued by considering the observation of system parameters and states. Hospitalization and mortality data are available at the scale of Isère for this purpose. Finally, this model
<figure><img src="image_9.png" /><figcaption>Figure 7: $I(t, \alpha)$ in 2D. Warmer colors correspond to higher individual densities.</figcaption></figure>

can be used for optimal control of vaccination. The vaccination laws $\theta(t)$ and $\theta_R(t, \tilde{\alpha})$ could be computed so that the $H$ or $I$ compartments respect health resource constraints. For instance, the $H$ compartment should not exceed the maximal capacity of the hospitals in a given area. An important conclusion of our work is that the prediction of the spread of an epidemic differs depending on whether the model used incorporates age structure. Therefore, it seems important to take this structure into account when designing a control strategy.

## A Completeness of ($Y_{T'}^{\xi}, ||.||_{X_{T'}}$)

*Proof of Lemma 1.* We first note that $Y_{T'}^{\xi}$ is a non-empty space as $0 \in Y_{T'}^{\xi}$, and that $(X_{T'}, ||.||_{X_{T'}})$ is a Banach space as a product of Banach spaces. Thus, to prove that $(Y_{T'}^{\xi}, ||.||_{X_{T'}})$ is a Banach space, it suffices to prove that $Y_{T'}^{\xi}$ is a closed subset of $X_{T'}$. Let $L_{+}^{1}$ the subset of functions in $(L^{1}([0, T]))^{6}$ almost everywhere positive. We will prove that $Y_{T'}^{\xi} = \Gamma_{1}^{-1}(L_{+}^{1}) \cap \Gamma_{2}^{-1}(L_{+}^{1})$ with $\Gamma_{1}$ and $\Gamma_{2}$ continuous functions on $X_{T'}$. Then, as $L_{+}^{1}$ is a closed subset of $(L^{1}([0, T]))^{6}$, it entails that $Y_{T'}^{\xi}$ is a closed subset in $X_{T'}$. Define $\Gamma_{1}$ as:

$$ \Gamma_1: \left\{ \begin{array}{ll} X_{T'} & \rightarrow (L^1([0, T]))^6 \\ (x_1, x_2, x_3, x_4, x_5, x_6) & \mapsto \left( x_1, x_2, x_3, \int_0^\infty x_4(., \alpha')d\alpha', \int_0^\infty x_5(., \tilde{\alpha}')d\tilde{\alpha}', \int_0^\infty x_6(., \tilde{\alpha}')d\tilde{\alpha}' \right). \end{array} \right. $$
<figure><img src="image_10.png" /><figcaption>Figure 8: R(t, α) in 2D. Warmer colors correspond to higher individual densities.</figcaption></figure>

For $x = (x_1, x_2, x_3, x_4, x_5, x_6)$, $y = (y_1, y_2, y_3, y_4, y_5, y_6) \in X_{T'}$:

$$ \begin{aligned} ||\Gamma_1(x) - \Gamma_1(y)||_{L^1} &= \int_0^{T'} (|x_1 - y_1|(t) + |x_2 - y_2|(t) + |x_3 - y_3|) (t)dt \\ &\quad + \int_0^{T'} \left( \left| \int_0^\infty (x_4(., \alpha') - y_4(., \alpha')) d\alpha' \right| + \left| \int_0^\infty (x_5(., \tilde{\alpha}') - y_5(., \tilde{\alpha}')) d\tilde{\alpha}' \right| \right) (t)dt \\ &\quad + \int_0^{T'} \left( \left| \int_0^\infty (x_6(., \tilde{\alpha}') - y_6(., \tilde{\alpha}')) d\tilde{\alpha}' \right| \right) (t)dt \\ &\le ||x - y||_{X_{T'}} \end{aligned} $$

so that $\Gamma_1$ is Lipschitz continuous. The same holds true for $\Gamma_2$ defined as:

$$ \Gamma_2: \left\{ \begin{array}{ll} X_{T'} & \to (L^1([0, T]))^6 \\ x & \mapsto \left( \bar{S}(., \xi) - x_1, \bar{E}(., \xi) - x_2, \bar{H}(., \xi) - x_3, \right. \\ & \qquad \left. \bar{I}(., \xi) - \int_0^\infty x_4(., \alpha') d\alpha', \bar{R}(., \xi) - \int_0^\infty x_5(., \tilde{\alpha}') d\tilde{\alpha}', \bar{V}(., \xi) - \int_0^\infty x_6(., \tilde{\alpha}') d\tilde{\alpha}' \right). \end{array} \right. $$

It concludes the proof of Lemma 1.

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

## B Proof of Theorem 1

The aim of this section is to complete the proof of Theorem 1. The first step consists in proving Lemma 3 stated in Section 3.1. The proof of Lemma 3 relies on two intermediate lemmas stated hereafter.

**Lemma 4.** Let $\xi^a \in \Xi$ and $\eta > 0$. There exists a non-negative constant $K(\xi^a, \eta)$ such that for all $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_\Xi < \eta$, for all $t \in [0, T]$, and for $i \in \{a, b\}$:

$$ \bar{S}(\xi^i, t), \bar{E}(\xi^i, t), \bar{H}(\xi^i, t), \bar{I}(\xi^i, t), \bar{R}(\xi^i, t), \bar{V}(\xi^i, t) \le K(\xi^a, \eta). $$
<figure><img src="image_12.png" /><figcaption>Figure 9: $V(t, \alpha)$ in 2D. Warmer colors correspond to higher individual densities.</figcaption></figure>

The proof of Lemma 4 is postponed to Appendix C.

**Lemma 5** (Proposition 3.16 in Blandfort (2015)). Let $T' \in (0, T]$. Let $\xi^a \in \Xi, \eta > 0$, $u$ and $v$ non-negative functions on $[0, T'] \times X_{T'} \times \Xi$ and $z_0$ a function on $\Xi$. For $t \in [0, T']$, $x \in X_{T'}, \xi \in \Xi$, we define:

$$Z(t, x, \xi) := \exp\left(-\int_{0}^{t} u(\tau, x, \xi)d\tau\right) \left(z_{0}(\xi) + \int_{0}^{t} \exp\left(\int_{0}^{\tau} u(\tau', x, \xi)d\tau'\right) v(\tau, x, \xi)d\tau\right).$$

We assume that for any $y^a \in Y_{T'}^{\xi^a}$, there exist two non-negative constants $C(\xi^a, \eta)$ and $C'(\xi^a, \eta)$ such that we have, for any $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_{\Xi} < \eta$, any $y^b \in Y_{T'}^{\xi^b}$ and any $t \in [0, T']$:

$$\begin{aligned} & u(t, y^a, \xi^a), u(t, y^b, \xi^b), v(t, y^a, \xi^a), v(t, y^b, \xi^b) \le C(\xi^a, \eta), \\ & |u(t, y^a, \xi^a) - u(t, y^b, \xi^b)| \le C'(\xi^a, \eta)||y^a - y^b||_t, \\ & |v(t, y^a, \xi^a) - v(t, y^b, \xi^a)| \le C'(\xi^a, \eta)||y^a - y^b||_t. \end{aligned}$$

Then, for any $y^a \in Y_{T'}^{\xi^a}$, there exists a non-negative constant $\tilde{C}(\xi^a, \eta)$ such that we have, for any $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_{\Xi} < \eta$, any $y^b \in Y_{T'}^{\xi^b}$ and any $t \in [0, T']$:

$$|Z(t, y^a, \xi^a) - Z(t, y^b, \xi^b)| \le \tilde{C}(\xi^a, \eta) \left( |z_0(\xi^a) - z_0(\xi^b)| + \int_0^t ||y^a - y^b||_\tau d\tau \right).$$

We now use Lemmas 4 and 5 to prove Propositions 5 and 6 below.

**Proposition 5.** Let $T' \in (0, T]$. Let $\xi^a \in \Xi, \eta > 0$. There exist non-negative constants $C_S(\xi^a, \eta)$, $C_E(\xi^a, \eta)$, $C_H(\xi^a, \eta)$ such that we have, for any $\xi^b \in \Xi$
<figure><img src="image_13.png" /><figcaption>Figure 10: Simulation of the SEHIRVD model without ages.</figcaption></figure>

with $||\xi^a - \xi^b||_\Xi < \eta$, any $y^a \in Y_{T'}^{\xi^a}$, any $y^b \in Y_{T'}^{\xi^b}$ and for almost all $t \in [0, T']$:

$$|S^{\xi^a}[y^a](t) - S^{\xi^b}[y^b](t)| \le C_S(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right),$$

$$|E^{\xi^a}[y^a](t) - E^{\xi^b}[y^b](t)| \le C_E(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right),$$

$$|H^{\xi^a}[y^a](t) - H^{\xi^b}[y^b](t)| \le C_H(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right).$$

**Proposition 6.** Let $T' \in (0, T]$. Let $\xi^a \in \Xi, \eta > 0$. There exist non-negative constants $C_I(\xi^a, \eta)$, $C_R(\xi^a, \eta)$, $C_V(\xi^a, \eta)$ such that we have, for any $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_\Xi < \eta$, any $y^a \in Y_{T'}^{\xi^a}$, any $y^b \in Y_{T'}^{\xi^b}$ and for almost all $t \in [0, T']$:

$$\int_0^\infty |I^{\xi^a}[y^a](t, \alpha') - I^{\xi^b}[y^b](t, \alpha')| d\alpha' \le C_I(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right),$$

$$\int_0^\infty |R^{\xi^a}[y^a](t, \tilde{\alpha}') - R^{\xi^b}[y^b](t, \tilde{\alpha}')| d\tilde{\alpha}' \le C_R(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right),$$

$$\int_0^\infty |V^{\xi^a}[y^a](t, \tilde{\alpha}') - V^{\xi^b}[y^b](t, \tilde{\alpha}')| d\tilde{\alpha}' \le C_V(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right).$$

*Proof of Proposition 5.* We only prove the inequality on $S$, as the inequalities on $H$ and $E$ can be proven similarly. For $t \in [0, T']$, $x \in X_{T'}$, $\xi \in \Xi$, we define:

$$u^S(t, x, \xi) = \theta(t) + \int_0^\infty \beta(t, \alpha') x_4(t, \alpha') d\alpha'$$

$$v^S(t, x, \xi) = \int_0^\infty (\sigma_R(\tilde{\alpha}') x_5(t, \tilde{\alpha}') + \sigma_V(\tilde{\alpha}') x_6(t, \tilde{\alpha}')) d\tilde{\alpha}'$$

$$z_0^S(\xi) = S_0.$$
It is obvious that $u^S$ and $v^S$ are non-negative functions. Let now consider $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_\Xi < \eta$ and $y^a \in Y_{T'}^{\xi^a}, y^b \in Y_{T'}^{\xi^b}$. Then, for $i \in \{a, b\}$:

$$ 
\begin{aligned}
u^S(t, y^i, \xi^i) &= \theta(t) + \int_0^\infty \beta(t, \alpha') y_4^i(t, \alpha') d\alpha' \\
&\le 1 + \int_0^\infty y_4^i(t, \alpha') d\alpha' \le 1 + \bar{I}(\xi^i, t) \le 1 + K(\xi^a, \eta) \text{ by Lemma 4,}
\end{aligned}
$$

$$ 
\begin{aligned}
v^S(t, y^i, \xi^i) &= \int_0^\infty (\sigma_R(\tilde{\alpha}') y_5^i(t, \tilde{\alpha}') + \sigma_V(\tilde{\alpha}') y_6^i(t, \tilde{\alpha}')) d\tilde{\alpha}' \\
&\le \int_0^\infty (y_5^i(t, \tilde{\alpha}') + y_6^i(t, \tilde{\alpha}')) d\tilde{\alpha}' \le \bar{R}(\xi^i, t) + \bar{V}(\xi^i, t) \le 2K(\xi^a, \eta) \text{ by Lemma 4.}
\end{aligned}
$$

Therefore, with $C(\xi^a, \eta) = 2 \max\{1, K(\xi^a, \eta)\}$, it follows that $u^S(t, y^i, \xi^i), v^S(t, y^i, \xi^i) \le C(\xi^a, \eta)$. We also have:

$$ 
\begin{aligned}
|u^S(t, y^a, \xi^a) - u^S(t, y^b, \xi^b)| &= \left| \int_0^\infty \beta(t, \alpha')(y_4^a(t, \alpha') - y_4^b(t, \alpha')) d\alpha' \right| \\
&\le \int_0^\infty |y_4^a(t, \alpha') - y_4^b(t, \alpha')| d\alpha' \le ||y^a - y^b||_t, \\
|v^S(t, y^a, \xi^a) - v^S(t, y^b, \xi^b)| &= \left| \int_0^\infty \sigma_R(\tilde{\alpha}') (y_5^a(t, \tilde{\alpha}') - y_5^b(t, \tilde{\alpha}')) d\tilde{\alpha}' + \int_0^\infty \sigma_V(\tilde{\alpha}') (y_6^a(t, \tilde{\alpha}') - y_6^b(t, \tilde{\alpha}')) d\tilde{\alpha}' \right| \\
&\le \int_0^\infty |y_5^a(t, \tilde{\alpha}') - y_5^b(t, \tilde{\alpha}')| d\tilde{\alpha}' + \int_0^\infty |y_6^a(t, \tilde{\alpha}') - y_6^b(t, \tilde{\alpha}')| d\tilde{\alpha}' \le 2 ||y^a - y^b||_t.
\end{aligned}
$$

Therefore we have, with $C'(\xi^a, \eta) = 2$,

$$ |u^S(t, y^a, \xi^a) - u^S(t, y^b, \xi^b)|, |v^S(t, y^a, \xi^a) - v^S(t, y^b, \xi^b)| \le C'(\xi^a, \eta) ||y^a - y^b||_t. $$

We finally obtain the inequality stated on $S$ by applying Lemma 5. $\square$

**Proof of Proposition 6.** We only prove the inequality on $I$, as the inequalities on $R$ and $V$ can be proven similarly. Let $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_\Xi < \eta$ and $y^a \in Y_{T'}^{\xi^a}, y^b \in Y_{T'}^{\xi^b}$. Then, for almost all $t \in [0, T]$:

$$ 
\begin{aligned}
\int_0^\infty |I^{\xi^a}[y^a] - I^{\xi^b}[y^b]| (t, \alpha') d\alpha' &= \int_0^t |I^{\xi^a}[y^a] - I^{\xi^b}[y^b]| (t, \alpha') d\alpha' + \int_t^\infty |I^{\xi^a}[y^a] - I^{\xi^b}[y^b]| (t, \alpha') d\alpha' \\
&=: I_1 + I_2.
\end{aligned}
$$

* For $I_1$, Definition 4 implies that

$$ 
\begin{aligned}
I_1 &= \int_0^t |h_I^a - h_I^b|(t - \alpha') \exp\left(-\int_0^{\alpha'} \gamma_I(\tilde{u}) + \tau_A \theta(\tilde{u} + t - \alpha') d\tilde{u}\right) d\alpha' \\
&\le \int_0^t |h_I^a - h_I^b|(t - \alpha') d\alpha' \le T ||h_I^a - h_I^b||_\infty \le T ||\xi^a - \xi^b||_\Xi.
\end{aligned}
$$

* For $I_2$, Definition 4 implies that

$$ 
\begin{aligned}
I_2 &= \int_t^\infty |f_I^a - f_I^b|(t - \alpha') \exp\left(-\int_0^{\alpha'} \gamma_I(\tilde{u} + \alpha' - t) + \tau_A \theta(\tilde{u}) d\tilde{u}\right) d\alpha' \\
&\le \int_0^t |f_I^a - f_I^b|(t - \alpha') d\alpha' \le ||f_I^a - f_I^b||_1 \le ||\xi^a - \xi^b||_\Xi.
\end{aligned}
$$
Therefore, by choosing $C_I(\xi^a, \eta) = T' + 1$, we prove the inequality on $I$ stated in Proposition 6. $\square$

Then, from Propositions 5 and 6 we prove Proposition 7 and Corollary 8 below.

**Proposition 7.** Let $T' \in (0, T]$. Let $\xi^a \in \Xi, \eta > 0$. There exists a non-negative constant $C_{Sol}(\xi^a, \eta)$ such that for any $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_\Xi < \eta$, any $y^a \in Y_{T'}^{\xi^a}$, any $y^b \in Y_{T'}^{\xi^b}$ and for almost all $t \in [0, T']$:

$$ ||Sol^{\xi^a}[y^a] - Sol^{\xi^b}[y^b]||_t \le C_{Sol}(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right). $$

*Proof of Proposition 7.* This result directly follows from the definition of $||\cdot||_t$ (see Definition 2), Proposition 5 and Proposition 6. $\square$

**Corollary 8.** Let $T' \in (0, T]$. Let $\xi^a \in \Xi, \eta > 0$. There exists a non-negative constant $C_{Sol}(\xi^a, \eta)$ such that for any $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_\Xi < \eta$, any $y^a \in Y_{T'}^{\xi^a}$, any $y^b \in Y_{T'}^{\xi^b}$ and for almost all $t \in [0, T']$:

$$ ||Sol^{\xi^a}[y^a] - Sol^{\xi^b}[y^b]||_{X_T'} \le T'C_{Sol}(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + ||y^a - y^b||_{X_T'} \right). $$

*Proof of Corollary 8.* From Proposition 7 and by definition of $||\cdot||_{X_{T'}}$ we have, for any $\xi^b \in \Xi$ with $||\xi^a - \xi^b||_\Xi < \eta$, any $y^a \in Y_{T'}^{\xi^a}$, any $y^b \in Y_{T'}^{\xi^b}$ and for almost all $t \in [0, T']$:

$$ \begin{aligned} ||Sol^{\xi^a}[y^a] - Sol^{\xi^b}[y^b]||_{X_T'} &= \int_0^{T'} ||Sol^{\xi^a}[y^a] - Sol^{\xi^b}[y^b]||_t dt \\ &\le \int_0^{T'} C_{Sol}(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^t ||y^a - y^b||_\tau d\tau \right) dt \\ &\le \int_0^{T'} C_{Sol}(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + \int_0^{T'} ||y^a - y^b||_\tau d\tau \right) dt \\ &\le T' C_{Sol}(\xi^a, \eta) \left( ||\xi^a - \xi^b||_\Xi + ||y^a - y^b||_{X_T'} \right). \end{aligned} $$

$\square$

*Proof of Lemma 3.* Lemma 3 directly follows from Corollary 8 with $\xi^b = \xi^a$. $\square$

*Proof of Theorem 1.* As already mentioned in Section 3, we first consider, for any $0 < T' \le T$, the solution operator $Sol^\xi : Y_{T'}^\xi \to Y_{T'}^\xi$. We know from Proposition 2 that a fixed-point of this operator is solution of the coupled system defined by (1) on $[0, T']$. From Lemma 3, there exists $0 < T' \le T$ such that the solution operator $Sol^\xi : Y_{T'}^\xi \to Y_{T'}^\xi$, is a contraction. Therefore, by Banach fixed-point theorem, it has a unique fixed point that is the solution of the coupled system defined by (1) on $[0, T']$. Note that it is always possible to choose $T'$ such that $\frac{T}{T'} = N \in \mathbb{N}^*$. It remains now to prove by induction the existence and unicity of a solution on each interval $[0, nT']$, for $n \le N$ so that finally, for $n = N$, Theorem 1 is proved.
* Initialization step (n = 1): Already proven, using Proposition 2.

* Induction step: Assume that the unique solution $y^n = (y_1^n, y_2^n, y_3^n, y_4^n, y_5^n, y_6^n) \in Y_{nT'}^\xi$ of the coupled system on $[0, nT']$ is built for $n \in \mathbb{N}$ with $1 \le n \le N-1$. Now, the goal is to build the unique solution $y^{n+1} \in Y_{(n+1)T'}^\xi$ of the coupled system on $[0, (n+1)T']$. In order to do so, we define $\tilde{Y}_{(n+1)T'}^\xi \subseteq Y_{(n+1)T'}^\xi$ as the set of functions $y \in Y_{(n+1)T'}^\xi$ such that for almost all $t \in [0, nT']$, $(y_1^\xi, y_2^\xi, y_3^\xi, y_4^\xi, y_5^\xi, y_6^\xi)(t) = (y_1^n, y_2^n, y_3^n, y_4^n, y_5^n, y_6^n)(t)$. Note that $(Y_{(n+1)T'}^\xi, ||\cdot||_{X_{(n+1)T'}})$ is a non-empty complete metric space as $(n+1)T' \le T$. Moreover, $\tilde{Y}_{(n+1)T'}^\xi$ is a closed subset of $(Y_{(n+1)T'}^\xi, ||\cdot||_{X_{(n+1)T'}})$. Therefore, $(\tilde{Y}_{(n+1)T'}^\xi, ||\cdot||_{X_{(n+1)T'}})$ is also a non-empty complete metric space. Moreover, $\text{Sol}^\xi : \tilde{Y}_{(n+1)T'}^\xi \to \tilde{Y}_{(n+1)T'}^\xi$. Indeed, if $y \in \tilde{Y}_{(n+1)T'}^\xi$, then $\text{Sol}^\xi[y] \in Y_{(n+1)T'}^\xi$, thanks to the definition of the bounds in Definition 2, and as for almost all $t \in [0, nT']$, $\text{Sol}^\xi[y](t) = \text{Sol}^\xi[y^n](t) = y^n(t)$. It remains now to prove that $\text{Sol}^\xi$ is a contraction on $(\tilde{Y}_{(n+1)T'}^\xi, ||\cdot||_{X_{(n+1)T'}})$. Consider $f, g \in \tilde{Y}_{(n+1)T'}^\xi$, by definition of $\tilde{Y}_{(n+1)T'}^\xi$ and by Proposition 7 we get:

$$ \begin{aligned} ||\text{Sol}^\xi[f] - \text{Sol}^\xi[g]||_{X_{(n+1)T'}} &= \int_0^{(n+1)T'} ||\text{Sol}^\xi[f] - \text{Sol}^\xi[g]||_t dt \\ &= \int_{nT'}^{(n+1)T'} ||\text{Sol}^\xi[f] - \text{Sol}^\xi[g]||_t dt \\ &\le \int_{nT'}^{(n+1)T'} C_{\text{Sol}}(\xi, \eta) \int_0^t ||f - g||_\tau d\tau dt \\ &\le \int_{nT'}^{(n+1)T'} C_{\text{Sol}}(\xi, \eta) \int_0^{(n+1)T'} ||f - g||_\tau d\tau dt \\ &\le T' C_{\text{Sol}}(\xi, \eta) \int_0^{(n+1)T'} ||f - g||_\tau d\tau. \end{aligned} $$

Finally, $||\text{Sol}^\xi[f] - \text{Sol}^\xi[g]||_{X_{(n+1)T'}} \le T' C_{\text{Sol}}(\xi, \eta) ||f - g||_{X_{(n+1)T'}}$. Since $T' C_{\text{Sol}}(\xi, \eta) < 1$ it follows that $\text{Sol}^\xi$ is a contraction on $(\tilde{Y}_{(n+1)T'}^\xi, ||\cdot||_{X_{(n+1)T'}})$. Using Banach fixed-point theorem, this implies that $\text{Sol}^\xi$ has a unique fixed point $y^{n+1}$ on $\tilde{Y}_{(n+1)T'}^\xi$. Remark that $y^{n+1}$ is also a fixed-point of $\text{Sol}^\xi$ on $Y_{(n+1)T'}^\xi$ as $\tilde{Y}_{(n+1)T'}^\xi \subseteq Y_{(n+1)T'}^\xi$. Thus, from Proposition 2, the continuous representative of $y^{n+1}$ is a solution to the coupled system described by (1) on $[0, (n+1)T']$. It remains to prove that it is the unique solution to the coupled system described by (1) on $[0, (n+1)T']$. Suppose $x^{n+1} \in Y_{(n+1)T'}^\xi$ is another solution of the coupled system on $[0, (n+1)T']$. Then, by Proposition 2, $x^{n+1}$ is a fixed point of $\text{Sol}^\xi$ on $[0, (n+1)T']$ thus also on $[0, nT']$. Thus, by induction hypothesis, $x^{n+1} = y^n$ almost everywhere in $[0, nT']$. Therefore, $x^{n+1} \in \tilde{Y}_{(n+1)T'}^\xi$ and is a fixed-point of $\text{Sol}^\xi$. By uniqueness of the fixed point of $\text{Sol}^\xi$ on $\tilde{Y}_{(n+1)T'}^\xi$, we have $x^{n+1} = y^{n+1}$. Thus the coupled system described by (1) has a unique
solution $y^{n+1} \in Y_{(n+1)T'}^{\xi}$ on $[0, (n+1)T']$. This concludes the induction step. $\square$

### C Proof of intermediate lemmas

*Proof of Lemma 2.* Let $T' \in (0, T]$ and $\xi \in \Xi$. We have to prove that for any $y \in Y_{T'}^{\xi}$,

$$S^{\xi}[y](t) \in [0, \bar{S}(\xi, t)], \quad E^{\xi}[y](t) \in [0, \bar{E}(\xi, t)], \quad H^{\xi}[y](t) \in [0, \bar{H}(\xi, t)], \quad (15)$$

$$\int_{0}^{\infty} I^{\xi}[y](t, \alpha')d\alpha' \in [0, \bar{I}(\xi, t)], \quad \int_{0}^{\infty} R^{\xi}[y](t, \tilde{\alpha}')d\tilde{\alpha}' \in [0, \bar{R}(\xi, t)] \quad (16)$$

$$\text{and} \quad \int_{0}^{\infty} V^{\xi}[y](t, \tilde{\alpha}')d\tilde{\alpha}' \in [0, \bar{V}(\xi, t)], \quad (17)$$

with $\bar{S}(\xi, t)$, $\bar{E}(\xi, t)$, $\bar{H}(\xi, t)$, $\bar{I}(\xi, t)$, $\bar{R}(\xi, t)$ and $\bar{V}(\xi, t)$ given in Definition 2. From Definition 4 and from the definition of $\Xi$, all the components of the solution operator are non-negative functions. For the upper-bounds, we first prove (16) and (17). We have, for almost all $t \in [0, T']$:

$$\begin{aligned} \int_{0}^{\infty} I^{\xi}[y](t, \alpha')d\alpha' &= \int_{0}^{t} I^{\xi}[y](t, \alpha')d\alpha' + \int_{t}^{\infty} I^{\xi}[y](t, \alpha')d\alpha' \\ &= \int_{0}^{t} h_{I}(t - \alpha') \exp\left(-\int_{0}^{\alpha'} \gamma_{I}(\tilde{u}) + \tau_{A}\theta(\tilde{u} + t - \alpha')d\tilde{u}\right) d\alpha' \\ &\quad + \int_{t}^{\infty} f_{I}(\alpha' - t) \exp\left(-\int_{0}^{t} \gamma_{I}(\tilde{u} + \alpha' - t) + \tau_{A}\theta(\tilde{u})d\tilde{u}\right) d\alpha' \\ &\le \int_{0}^{t} h_{I}(t - \alpha')d\alpha' + \int_{t}^{\infty} f_{I}(\alpha' - t)d\alpha' \\ &\le ||h_{I}||_{\infty}t + ||f_{I}||_{1} = \bar{I}(\xi, t). \end{aligned}$$

Similarly we can prove:

$$\int_{0}^{\infty} R^{\xi}[y](t, \tilde{\alpha}')d\tilde{\alpha}' \le ||h_{R}||_{\infty}t + ||f_{R}||_{1} = \bar{R}(\xi, t),$$

$$\int_{0}^{\infty} V^{\xi}[y](t, \tilde{\alpha}')d\tilde{\alpha}' \le ||h_{V}||_{\infty}t + ||f_{V}||_{1} = \bar{V}(\xi, t).$$

Let us now prove the upper-bound in (15). We have:

$$S^{\xi}[y](t) = \exp\left(-\int_{0}^{t} F(\tau)d\tau\right)\left(S_{0} + \int_{0}^{t} \exp\left(\int_{0}^{\tau} F(\tau')d\tau'\right)G(\tau)d\tau\right)$$

with

$$\begin{aligned} G(\tau) &= \int_{0}^{\infty} (\sigma_{R}(\tilde{\alpha}')y_{5}(\tau, \tilde{\alpha}') + \sigma_{V}(\tilde{\alpha}')y_{6}(\tau, \tilde{\alpha}'))d\tilde{\alpha}' \\ &\le \int_{0}^{\infty} (y_{5}(\tau, \tilde{\alpha}') + y_{6}(\tau, \tilde{\alpha}'))d\tilde{\alpha}' \le \bar{R}(\xi, \tau) + \bar{V}(\xi, \tau) \end{aligned}$$
as $0 \le \sigma_R, \sigma_V \le 1$. Also, $F(\tau) = \theta(\tau) + \int_0^\infty \beta(\tau, \alpha')y_4(\tau, \alpha')d\alpha' \ge 0$. We therefore have

$$
\begin{aligned}
S^\xi[y](t) &\le \exp\left(-\int_0^t F(\tau)d\tau\right) \left(S_0 + \int_0^t \exp\left(\int_0^\tau F(\tau')d\tau'\right) \left(\bar{R}(\xi, \tau) + \bar{V}(\xi, \tau)\right) d\tau\right) \\
&\le S_0 \exp\left(-\int_0^t F(\tau)d\tau\right) + \int_0^t \exp\left(-\int_\tau^t F(\tau')d\tau'\right) \left(\bar{R}(\xi, \tau) + \bar{V}(\xi, \tau)\right) d\tau \\
&\le S_0 + \int_0^t \left(\bar{R}(\xi, \tau) + \bar{V}(\xi, \tau)\right) d\tau \\
&\le S_0 + \int_0^t \left(\left(||f_R||_1 + ||f_V||_1\right) + \left(||h_R||_\infty + ||h_V||_\infty\right)\tau\right) d\tau \\
&\le S_0 + \left(||f_R||_1 + ||f_V||_1\right)t + \frac{||h_R||_\infty + ||h_V||_\infty}{2} t^2 = \bar{S}(\xi, t)
\end{aligned}
$$

Similar computations lead to:

$$
E^{\xi}[y](t) = e^{-\delta t} \left( E_0 + \int_0^t e^{\delta \tau} y_1(\tau) \left( \int_0^{\infty} \beta(\tau, \alpha') y_4(\tau, \alpha') d\alpha' \right) d\tau \right).
$$

Thus, as $0 \le \beta \le 1$:

$$
\begin{aligned}
E^{\xi}[y](t) &\le e^{-\delta t} \left( E_0 + \int_0^t e^{\delta \tau} y_1(\tau) \left( \int_0^{\infty} y_4(\tau, \alpha') d\alpha' \right) d\tau \right) \\
&\le e^{-\delta t} \left( E_0 + \int_0^t e^{\delta \tau} \bar{S}(\xi, \tau) \bar{I}(\xi, \tau) d\tau \right) \\
&\le e^{-\delta t} E_0 + \left( \int_0^t e^{-\delta(t-\tau)} \bar{S}(\xi, \tau) \bar{I}(\xi, \tau) d\tau \right) \\
&\le E_0 + \left( \int_0^t \bar{S}(\xi, \tau) \bar{I}(\xi, \tau) d\tau \right) \\
&\le E_0 + \int_0^t \left( S_0 + (||f_R||_1 + ||f_V||_1)\tau + \frac{||h_R||_\infty + ||h_V||_\infty}{2} \tau^2 \right) (||f_I||_1 + ||h_I||_\infty \tau) d\tau \\
&\le E_0 + \int_0^t \left( S_0 ||f_I||_1 + \left(S_0 ||h_I||_\infty + (||f_R||_1 + ||f_V||_1)||f_I||_1\right)\tau \right. \\
&\qquad \left. + \frac{||h_R||_\infty + ||h_V||_\infty}{2} ||f_I||_1 + (||f_R||_1 + ||f_V||_1)||h_I||_\infty \right) \tau^2 + \left( \frac{||h_R||_\infty + ||h_V||_\infty}{2} ||h_I||_\infty \right) \tau^3 d\tau \\
&\le E_0 + S_0 ||f_I||_1 t + \left( \frac{S_0 ||h_I||_\infty}{2} + \frac{||f_R||_1 + ||f_V||_1}{2} ||f_I||_1 \right) t^2 \\
&\qquad + \left( \frac{||h_R||_\infty + ||h_V||_\infty}{6} ||f_I||_1 + \frac{||f_R||_1 + ||f_V||_1}{3} ||h_I||_\infty \right) t^3 + \left( \frac{||h_R||_\infty + ||h_V||_\infty}{8} ||h_I||_\infty \right) t^4 \\
&= \bar{E}(\xi, t).
\end{aligned}
$$

Finally, we write $H^\xi[y](t)$ as:

$$
H^{\xi}[y](t) = e^{-\gamma_H t} \left( H_0 + \int_0^t e^{\gamma_H \tau} \left( \int_0^{\infty} \gamma_I(\alpha') \tau_H y_4(\tau, \alpha') d\alpha' \right) d\tau \right).
$$
Thus, as $0 \le \gamma_I, \tau_H \le 1$:

$$
\begin{align*}
H^{\xi}[y](t) &\le e^{-\gamma_H t} \left( H_0 + \int_0^t e^{\gamma_H \tau} \left( \int_0^\infty y_4(\tau, \alpha') d\alpha' \right) d\tau \right) \\
&\le e^{-\gamma_H t} H_0 + \left( \int_0^t e^{-\gamma_H (t-\tau)} \left( \int_0^\infty y_4(\tau, \alpha') d\alpha' \right) d\tau \right) \\
&\le H_0 + \left( \int_0^t \left( \int_0^\infty y_4(\tau, \alpha') d\alpha' \right) d\tau \right) \\
&\le H_0 + \int_0^t \bar{I}(\xi, \tau) d\tau \le H_0 + \int_0^t \left( \|f_I\|_1 + \|h_I\|_\infty \tau \right) d\tau \\
&\le H_0 + \|f_I\|_1 t + \frac{\|h_I\|_\infty}{2} t^2 = \bar{H}(\xi, t).
\end{align*}
$$

□

*Proof of Lemma 4.* We only make the proof for $\bar{S}$ and $\bar{I}$, as the proof for $\bar{E}$ and $\bar{H}$ is similar to the one for $\bar{S}$ and the proof for $\bar{R}$ and $\bar{V}$ is similar to the one for $\bar{I}$. Let $\xi^b \in \Xi$ such that $\|\xi^a - \xi^b\|_{\Xi} < \eta$ and let $t \in [0, T]$. We have:

$$
\begin{align*}
\bar{S}(\xi^a, t) &= S_0^a + (\|f_R^a\|_1 + \|f_V^a\|_1)t + \frac{\|h_R^a\|_\infty + \|h_V^a\|_\infty}{2}t^2 \\
&\le S_0^a + (\|f_R^a\|_1 + \|f_V^a\|_1)T + \frac{\|h_R^a\|_\infty + \|h_V^a\|_\infty}{2}T^2 := K_S^a(\xi^a, \eta).
\end{align*}
$$

Also, using the triangle inequality:

$$
\begin{align*}
\bar{S}(\xi^b, t) &= S_0^b + (\|f_R^b\|_1 + \|f_V^b\|_1)t + \frac{\|h_R^b\|_\infty + \|h_V^b\|_\infty}{2}t^2 \\
&\le S_0^b - S_0^a + (\|f_R^b - f_R^a\|_1 + \|f_V^b - f_V^a\|_1)t + \frac{\|h_R^b - h_R^a\|_\infty + \|h_V^b - h_R^a\|_\infty}{2}t^2 + S_0^a \\
&\quad + (\|f_R^a\|_1 + \|f_V^a\|_1)t + \frac{\|h_R^a\|_\infty + \|h_V^a\|_\infty}{2}t^2 \\
&\le \eta + \eta T + \eta T^2 + K_S^a(\xi^a, \eta) := K_S^b(\xi^a, \eta).
\end{align*}
$$

This constant $K_S^b(\xi^a, \eta)$ does not depend on the choice of $\xi^b$. Similarly:

$$
\begin{align*}
\bar{I}(\xi^a, t) &= \|f_I^a\|_1 + \|h_I^a\|_\infty t \\
&\le \|f_I^a\|_1 + \|h_I^a\|_\infty T := K_I^a(\xi^a, \eta).
\end{align*}
$$

Using the triangle inequality, we have:

$$
\begin{align*}
\bar{I}(\xi^b, t) &= \|f_I^b\|_1 + \|h_I^b\|_\infty t \\
&\le \|f_I^b - f_I^a\|_1 + \|h_I^b - h_I^a\|_\infty t + \|f_I^a\|_1 + \|h_I^a\|_\infty t \\
&\le \eta + \eta T + K_I^a(\xi^a, \eta) := K_I^b(\xi^a, \eta).
\end{align*}
$$

Therefore, by defining $K(\xi^a, \eta)$ as the maximum of the $K_Z^i(\xi^a, \eta)$ over all $Z \in \{S, E, H, I, R, V\}$ and $i \in \{a, b\}$, the result follows. $\square$
## D Proof of the continuity of the solution

*Proof of Theorem 3.* Proposition 2 implies that $y^a$ is a fixed point of $\text{Sol}^{\xi^a}$ and $y^b$ is a fixed point of $\text{Sol}^{\xi^b}$. Therefore, for $t \in [0, T]$,

$$ ||y^a - y^b||_t = ||\text{Sol}^{\xi^a}[y^a] - \text{Sol}^{\xi^b}[y^b]||_t. $$

Then, by Proposition 7,

$$ \begin{aligned} ||y^a - y^b||_t &\le C_{\text{Sol}}(\xi^a, \eta) \left( ||\xi^a - \xi^b||_{\Xi} + \int_0^t ||y^a - y^b||_{\tau} d\tau \right) \\ &\le C_{\text{Sol}}(\xi^a, \eta) ||\xi^a - \xi^b||_{\Xi} + \int_0^t C_{\text{Sol}}(\xi^a, \eta) ||y^a - y^b||_{\tau} d\tau. \end{aligned} $$

Then, by using Gronwall's inequality with $f(t) = C_{\text{Sol}}(\xi^a, \eta)||\xi^a - \xi^b||_{\Xi}$ and $g(t) = C_{\text{Sol}}(\xi^a, \eta)$, we obtain:

$$ \begin{aligned} ||y^a - y^b||_t &\le C_{\text{Sol}}(\xi^a, \eta)||\xi^a - \xi^b||_{\Xi} \\ &\quad + \int_0^t C_{\text{Sol}}(\xi^a, \eta)||\xi^a - \xi^b||_{\Xi} C_{\text{Sol}}(\xi^a, \eta) \exp\left(\int_{\tau}^t C_{\text{Sol}}(\xi^a, \eta) d\tau'\right) d\tau \\ &\le C_{\text{Sol}}(\xi^a, \eta)||\xi^a - \xi^b||_{\Xi} + \int_0^t (C_{\text{Sol}}(\xi^a, \eta))^2 ||\xi^a - \xi^b||_{\Xi} \exp(TC_{\text{Sol}}(\xi^a, \eta)) d\tau \\ &\le (C_{\text{Sol}}(\xi^a, \eta) + T(C_{\text{Sol}}(\xi^a, \eta))^2) \exp(TC_{\text{Sol}}(\xi^a, \eta)) ||\xi^a - \xi^b||_{\Xi}. \end{aligned} $$

Finally,

$$ ||y^a - y^b||_{X_T} = \int_0^T ||y^a - y^b||_t dt \le C(\xi^a, \eta)||\xi^a - \xi^b||_{\Xi} $$

with $C(\xi^a, \eta)$ given by

$$ T(C_{\text{Sol}}(\xi^a, \eta) + T(C_{\text{Sol}}(\xi^a, \eta))^2) \exp(TC_{\text{Sol}}(\xi^a, \eta)). $$

This implies that the solution of the coupled system described by (1) on $[0, T]$ depends continuously on the initial and boundary conditions. □

## E Proof of biological bounds for the solution

In order to prove Theorem 4, we need to permute the derivative and the integral operators. As stated in Lemma 6 below, this is the case under some appropriate assumptions.

**Lemma 6.** Suppose that Assumptions (4) to (7) of Theorem 4 hold true. Then,

$$ \begin{aligned} \frac{d}{dt} \left( \int_{0}^{\infty} I(., \alpha') d\alpha' \right) (t) &= \left( \int_{0}^{\infty} \frac{\partial I}{\partial t}(t, \alpha') d\alpha' \right), \\ \frac{d}{dt} \left( \int_{0}^{\infty} R(., \tilde{\alpha}') d\tilde{\alpha}' \right) (t) &= \left( \int_{0}^{\infty} \frac{\partial R}{\partial t}(t, \tilde{\alpha}') d\tilde{\alpha}' \right), \\ \frac{d}{dt} \left( \int_{0}^{\infty} V(., \tilde{\alpha}') d\tilde{\alpha}' \right) (t) &= \left( \int_{0}^{\infty} \frac{\partial V}{\partial t}(t, \tilde{\alpha}') d\tilde{\alpha}' \right). \end{aligned} $$
Let us now give the proof of Theorem 4.

*Proof of Theorem 4.* For $t = 0$, the result is directly given by (4) and the other assumptions. As mentioned in the proof of Lemma 6, $I(t, \cdot)$ is integrable on $[t, \infty]$. In particular, this means that $I(t, +\infty) := \lim_{\alpha \to \infty} I(t, \alpha)$ exists and is equal to 0. Using the same methodology, one can prove that for all $t \in [0, T]$, $R(t, +\infty) := \lim_{\tilde{\alpha} \to \infty} R(t, \tilde{\alpha})$ and $V(t, +\infty) := \lim_{\tilde{\alpha} \to \infty} V(t, \tilde{\alpha})$ exist and are both equal to 0. Let now $t \in [0, T]$ and compute:

$$ 
\begin{aligned}
\mathcal{D} &=: \frac{d}{dt} \left( S + E + H + \int_{0}^{\infty} I(\cdot, \alpha') d\alpha' + \int_{0}^{\infty} R(\cdot, \tilde{\alpha}') d\tilde{\alpha}' + \int_{0}^{\infty} V(\cdot, \tilde{\alpha}') d\tilde{\alpha}' + D \right) (t) \\
&= \frac{dS}{dt}(t) + \frac{dE}{dt}(t) + \frac{dH}{dt}(t) + \frac{d}{dt} \left( \int_{0}^{\infty} I(\cdot, \alpha') d\alpha' \right) (t) \\
&\quad + \frac{d}{dt} \left( \int_{0}^{\infty} R(\cdot, \tilde{\alpha}') d\tilde{\alpha}' \right) (t) + \frac{d}{dt} \left( \int_{0}^{\infty} V(\cdot, \tilde{\alpha}') d\tilde{\alpha}' \right) (t) + \frac{dD}{dt}(t).
\end{aligned}
$$

By Lemma 6, this is equivalent to:

$$ \mathcal{D} = \frac{dS}{dt}(t) + \frac{dE}{dt}(t) + \frac{dH}{dt}(t) + \int_{0}^{\infty} \frac{\partial I}{\partial t}(t, \alpha') d\alpha' + \int_{0}^{\infty} \frac{\partial R}{\partial t}(t, \tilde{\alpha}') d\tilde{\alpha}' + \int_{0}^{\infty} \frac{\partial V}{\partial t}(t, \tilde{\alpha}') d\tilde{\alpha}' \\ + \frac{dD}{dt}(t). \qquad (18) $$

Yet, using the PDE controlling $I$ in the system described by (1), we get:

$$ 
\begin{aligned}
\int_{0}^{\infty} \frac{\partial I}{\partial t}(t, \alpha') d\alpha' &= \int_{0}^{\infty} -\frac{\partial I}{\partial \alpha}(t, \alpha') - (\gamma_I(\alpha') + \theta(t)\tau_A)I(t, \alpha') d\alpha' \\
&= -\int_{0}^{\infty} \frac{\partial I}{\partial \alpha}(t, \alpha') - \int_{0}^{\infty} (\gamma_I(\alpha') + \theta(t)\tau_A)I(t, \alpha') d\alpha' \\
&= I(t, 0) - I(t, +\infty) - \int_{0}^{\infty} (\gamma_I(\alpha') + \theta(t)\tau_A)I(t, \alpha') d\alpha'.
\end{aligned}
$$

Finally

$$ \int_{0}^{\infty} \frac{\partial I}{\partial t}(t, \alpha') d\alpha' = \delta E(t) - \int_{0}^{\infty} (\gamma_I(\alpha') + \theta(t)\tau_A) I(t, \alpha') d\alpha'. \quad (19) $$

Similarly, we can write:

$$ \int_{0}^{\infty} \frac{\partial R}{\partial t}(t, \tilde{\alpha}') d\tilde{\alpha}' = \int_{0}^{\infty} \gamma_I(\alpha')(1 - \tau_H)I(t, \alpha') d\alpha' + \gamma_H(1 - \tau_D)H(t) \\ - \int_{0}^{\infty} (\theta_R(t, \tilde{\alpha}') + \sigma_R(\tilde{\alpha}))R(t, \tilde{\alpha}') d\tilde{\alpha}', \quad (20) $$

$$ \int_{0}^{\infty} \frac{\partial V}{\partial t}(t, \tilde{\alpha}') d\tilde{\alpha}' = \int_{0}^{\infty} \theta(t) \tau_A I(t, \alpha') d\alpha' + \int_{0}^{\infty} \theta_R(t, \tilde{\alpha}') R(t, \tilde{\alpha}') d\tilde{\alpha}' + \theta(t) S(t) \\ - \int_{0}^{\infty} \sigma_V(\tilde{\alpha}') V(t, \tilde{\alpha}') d\tilde{\alpha}'. \quad (21) $$

We now use (19), (20) and (21) to rewrite the right hand side in (18) and obtain $\mathcal{D} = 0$. Therefore, as (8) is satisfied for $t = 0$, it is also satisfied for all $t \in [0, T]$.
The quantities $S, E, H, I, R$ and $V$ are non-negative from Definition 4. Moreover, from the equation on $D$ in (1), it follows that $D(t) = D_0 + \int_0^t \gamma_H \tau_D H(\tilde{u}) d\tilde{u}$, so it is also non-negative. Finally, the upper-bound in Inequality (9) is deduced from (8). It remains to prove (10). From Definition 4, we have that for all $t \in [0, T]$, for all $\alpha, \tilde{\alpha} \ge 0$

$$I(t, \alpha) \le \max \{f_I(\alpha - t), h_I(t - \alpha)\} \le 1.$$

Similarly

$$R(t, \tilde{\alpha}) \le \max \{f_R(\tilde{\alpha} - t), h_R(t - \tilde{\alpha})\} \le 1 \text{ and } V(t, \tilde{\alpha}) \le \max \{f_V(\tilde{\alpha} - t), h_V(t - \tilde{\alpha})\} \le 1.$$

This completes the proof. □

*Proof of Lemma 6.* We only prove the first equality since the other two can be proven similarly. The proof uses a corollary of the Lebesgue's dominated convergence theorem.

* For all $t \in [0, T]$, $\alpha \mapsto I(t, \alpha)$ is Lebesgue-integrable on $\mathbb{R}_+$. Indeed, for $t \in [0, T]$, we have from Definition 4:

$$I(t, \alpha) = \begin{cases} f_I(\alpha - t) \exp\left(-\int_0^t \gamma_I(\tilde{u} + \alpha - t) + \tau_A \theta(\tilde{u}) d\tilde{u}\right) & \text{if } t \le \alpha, \\ h_I(t - \alpha) \exp\left(-\int_0^\alpha \gamma_I(\tilde{u}) + \tau_A \theta(\tilde{u} + t - \alpha) d\tilde{u}\right) & \text{if } t > \alpha. \end{cases}$$

Therefore if $t \in [0, T]$ and $\alpha > t$, the previous expression gives

$$|I(t, \alpha)| = \left| f_I(\alpha - t) \exp\left(-\int_0^t \gamma_I(\tilde{u} + \alpha - t) + \tau_A \theta(\tilde{u}) d\tilde{u}\right) \right| \\ \le f_I(\alpha - t).$$

Then

$$\int_t^\infty |I(t, \alpha')| d\alpha' = \int_t^\infty I(t, \alpha') d\alpha' \le \int_t^\infty f_I(\alpha' - t) d\alpha' \le 1.$$

by Assumption (4).

* For all $\alpha \in \mathbb{R}_+$, $t \mapsto I(t, \alpha)$ is differentiable on $[0, T]$ as $I \in C^1([0, T] \times \mathbb{R}_+)$.

* Now, we have to define a Lebesgue-integrable function $\phi : \mathbb{R}_+ \to \mathbb{R}_+$ such that for all $t \in [0, T]$, $\alpha \in \mathbb{R}_+$,

$$\left| \frac{\partial I}{\partial t}(t, \alpha) \right| \le \phi(\alpha).$$

As $I \in C^1([0, T] \times \mathbb{R}_+)$, $\frac{\partial I}{\partial t}$ is continuous on the compact set $[0, T]^2 \subseteq [0, T] \times \mathbb{R}_+$. Therefore there exists $M \in \mathbb{R}_+$ so that for all $(t, \alpha) \in [0, T]^2$

$$\left| \frac{\partial I}{\partial t}(t, \alpha) \right| \le M.$$
Now, consider $t \in [0, T], \alpha \in ]T, +\infty[$. Then

$$I(t, \alpha) = f_I(\alpha - t) \exp\left(-\int_0^t \gamma_I(\tilde{u} + \alpha - t) + \tau_A \theta(\tilde{u}) d\tilde{u}\right)$$

and therefore

$$\frac{\partial I}{\partial t}(t, \alpha) = -f'_I(\alpha - t) \exp\left(-\int_0^t \gamma_I(\tilde{u} + \alpha - t) + \tau_A \theta(\tilde{u}) d\tilde{u}\right) + f_I(\alpha - t)(-\tau_A \theta(t)) - \gamma_I(\alpha - t) \times \exp\left(-\int_0^t \gamma_I(\tilde{u} + \alpha - t) + \tau_A \theta(\tilde{u}) d\tilde{u}\right).$$

Thus, for $t \in [0, T], \alpha \in ]T, +\infty[$,

$$\left|\frac{\partial I}{\partial t}(t, \alpha)\right| \le \left|f'_I(\alpha - t)\right| + \left|f_I(\alpha - t)\right|(\tau_A \theta(t) + \gamma_I(\alpha - t)) \\ \le \left|f'_I(\alpha - t)\right| + 2\left|f_I(\alpha - t)\right| \le 3\Psi_1(\alpha).$$

We define $\phi: \mathbb{R}_+ \to \mathbb{R}_+$ as:

$$\phi(\alpha) = \begin{cases} M & \text{if } \alpha \le T, \\ 3\Psi_1(\alpha) & \text{if } \alpha > T. \end{cases}$$

The first equality in Lemma 6 thus follows from a corollary of the Lebesgue's dominated convergence theorem. □

## Data Availability Statement

No new data were generated or analysed in support of this research.

## Acknowledgments

The authors would like to thank Elisabeta Vegu, who sadly passed away, for her valuable contribution to the modeling.

## Funding

This work was partially funded by the MODCOV19 platform https://modcov19.math.cnrs.fr/ implemented in March 2020 by the National Institute for Mathematical Sciences and their Interactions of the CNRS to coordinate modelling activities in France in connection with the multiples facets of the pandemic crisis.

## References

* Barbarossa, M.V., Polner, M., Röst, G., 2018a. Temporal evolution of immunity distributions in a population with waning and boosting. bioRxiv URL: https://api.semanticscholar.org/CorpusID:52883594.
Barbarossa, M.V., Polner, M., Röst, G., et al., 2018b. Temporal evolution of immunity distributions in a population with waning and boosting. Complexity 2018.

Barbarossa, M.V., Röst, G., 2015. Immuno-epidemiology of a population struc-
tured by immune status: a mathematical study of waning immunity and
immune system boosting. Journal of mathematical biology 71, 1737–1770.

Batool, I., Bajcinca, N., 2021. Well-posedness of a coupled pde–ode model of stem cell lineage involving homeostatic regulation. Results in Applied Mathematics 9, 100135. doi:10.1016/j.rinam.2020.100135.

Bhattacharya, S., Adler, F.R., 2012. A time since recovery model with varying rates of loss of immunity. Bulletin of mathematical biology 74, 2810–2819.

Blandfort, P., 2015. Mathematical modeling of synaptic transmission at chemical synapses [master thesis]. Technische Universität Kaiserslautern.

Boëlle, P.Y., Delory, T., Maynadier, X., Janssen, C., Piarroux, R., Pichenot, M., Lemaire, X., Baclet, N., Weyrich, P., Melliez, H., et al., 2020. Trajectories of hospitalization in covid-19 patients: an observational study in france. Journal of clinical medicine 9, 3148.

Byambasuren, O., Cardona, M., Bell, K., Clark, J., McLaws, M.L., Glasziou, P., 2020. Estimating the extent of asymptomatic covid-19 and its potential for community transmission: systematic review and meta-analysis. Official Journal of the Association of Medical Microbiology and Infectious Disease Canada 5, 223–234.

Calba, C., Zhu-Soubise, A., Mahdjoub-Assaad, S., Villegas-Ramirez, B., Tarantola, A., Barataud, D., King, L., Morvan, E., et al., 2021. Linked transmission chains of imported sars-cov-2 variant b. 1.351 across mainland france, january 2021. Eurosurveillance 26, 2100333.

Chauhan, S., Misra, O.P., Dhar, J., 2014. Stability analysis of sir model with vaccination. American journal of computational and applied mathematics 4, 17–23.

Chekroun, A., Kuniya, T., 2020. An infection age-space structured sir epidemic model with neumann boundary condition. Applicable Analysis 99, 1972–1985.

Chen, Y., Yang, J., Zhang, F., 2013. The global stability of an sirs model with infection age. Mathematical Biosciences & Engineering 11, 449–469.

De Marchi, S., 2013. Four lectures on radial basis functions. Department of Mathematics, University of Padua.

Dietz, K., Schenzle, D., 1985. Proportionate mixing models for age-dependent infection transmission. Journal of mathematical biology 22, 117–120.

Duan, X., Yuan, S., Li, X., 2014. Global stability of an svir model with age of vaccination. Applied Mathematics and Computation 226, 528–540.
Ehrhardt, M., Gašper, J., Kilianová, S., 2019. Sir-based mathematical modeling of infectious diseases with vaccination and waning immunity. Journal of Computational Science 37, 101027.

Flyer, N., Fornberg, B., Bayona, V., Barnett, G.A., 2016. On the role of polynomials in rbf-fd approximations: I. interpolation and accuracy. Journal of Computational Physics 321, 21–38.

Fornberg, B., Flyer, N., 2015. A primer on radial basis functions with applications to the geosciences. SIAM, CBMS-NSF Regional Conference Series in Applied Mathematics.

Fornberg, B., Lehto, E., 2011. Stabilization of rbf-generated finite difference methods for convective pdes. Journal of Computational Physics, 2270–2285.

Heffernan, J., Keeling, M.J., 2009. Implications of vaccination and waning immunity. Proceedings of the Royal Society B: Biological Sciences 276, 2071–2080.

Kendal, D., 1956. Deterministic and stochastic epidemics in closed populations. Proc. 3rd Berkeley Symp. on Math. Statist. and Prob., 149–165.

Kermack, W.O., McKendrick, A.G., 1927. A contribution to the mathematical theory of epidemics. Proceedings of the royal society of London. Series A, Containing papers of a mathematical and physical character 115, 700–721.

Kuipers, L., Niederreiter, H., 2012. Uniform distribution of sequences. Courier Corporation.

Lemieux, C., 2009. Monte Carlo and Quasi-Monte Carlo Sampling. Springer New York. URL: http://dx.doi.org/10.1007/978-0-387-78165-5, doi:10.1007/978-0-387-78165-5.

Li, X.Z., Wang, J., Ghosh, M., 2010. Stability and bifurcation of an sirs epidemic model with treatment and age of vaccination. Applied Mathematical Modelling 34, 437–450.

Liu, L., Wang, J., Liu, X., 2015. Global stability of an seir epidemic model with age-dependent latency and relapse. Nonlinear Analysis: real world applications 24, 18–35.

Magal, P., McCluskey, C., 2013. Two-group infection age model including an application to nosocomial infection. SIAM Journal on Applied Mathematics 73, 1058–1095.

McCluskey, C.C., 2012. Global stability for an sei epidemiological model with continuous age-structure in the exposed and infectious classes. Mathematical Biosciences & Engineering 9, 819–841.

Okuwa, K., Inaba, H., Kuniya, T., 2021. An age-structured epidemic model with boosting and waning of immune status. Math. Biosci. Eng 18, 5707–5736.

Oliver, M., Georges, D., Prieur, C., 2022. Spatialized epidemiological forecasting applied to covid-19 pandemic at departmental scale in france. Systems & Control Letters 164, 105240.
Richard, Q., Alizon, S., Choisy, M., Sofonea, M.T., Djidjou-Demasse, R., 2021. Age-structured non-pharmaceutical interventions for optimal control of covid-19 epidemic. PLoS computational biology 17, e1008776.

Shulgin, B., Stone, L., Agur, Z., 1998. Pulse vaccination strategy in the sir epidemic model. Bulletin of mathematical biology 60, 1123–1148.

Zaman, G., Kang, Y.H., Jung, I.H., 2008. Stability analysis and optimal vaccination of an sir epidemic model. BioSystems 93, 240–249.