## Abstract

The Covid-19 pandemic outbreak was followed by an huge amount of modeling studies in order to rapidly gain insights to implement the best public health policies. However, most of those compartmental models used a classical ordinary differential equations (ODEs) system based formalism that came with the tacit assumption the time spent in each compartment does not depend of the time already spent in it. To overcome this “memoryless” issue, a widely used workaround is to artificially increase and chain the number of compartments of an unique reality (*e*.*g*. many compartments for infected individuals). It allows for a greater heterogeneity and thus be closer to the observed situation, at the cost of rendering the whole model more difficult to apprehend and parametrize. We propose here an alternative formalism based on a partial differential equations (PDEs) system instead of ordinary differential equations, which provides naturally a memory structure for each compartment, and thus allows to keep a restrained number of compartments. We use such a model applied to the French situation, accounting for vaccinal and natural immunity. The results seem to indicate that the vaccination rate is not enough to ensure the end of the epidemic, but, above all, highlight a huge uncertainty attributable to the age-structured contact matrix.

## 1 Introduction

Shortly after the Covid-19 outbreak late 2019, many efforts were put in diverse research areas to understand both the disease and its aetiological agent, SARS-CoV-2, and produce the tools needed to deal with what quickly became a pandemic. Among those areas, mathematical modeling studies proliferated to better grasp the epidemics’ dynamics on a —at first— short and medium-term scale. Stochastic models were more appropriate early on to take into account the randomness of the transmission events [Kucharski et al., 2020; Hellewell et al., 2020; Beneteau et al., 2021], but they were rapidly complemented by deterministic models since many epidemics were settled within a couple of months in many countries. These models’ results provided valuable insights to guide public health policies, often on a nationwide scale [RSTB, 2021].

Many of the deterministic models developed were variations of the classical Susceptible–Infected–Recovered (SIR) compartmental model, usually involving a system of ordinary differential equations (ODEs) [Keeling and Rohani, 2008]. Such a simple formalism was adapted at first given the unknowns regarding the natural history of the disease. However, the knowledge accumulated in a matter of months made it clear that several assumptions were biologically unrealistic. In particular, the residence times in various compartments were not distributed exponentially, and the related “lack of memory” (also named Markov property) —meaning that the time spent in a compartment does not depend on the time already spent in the compartment, as implied by the ODE formalism— became particularly detrimental to short-term forecasting (see Supplementary Figure S1). For example, for the time elapsed between infection and potential hospital admission, this memoryless hypothesis does not hold [Salje et al., 2020; Sofonea et al., 2021]. A workaround to this issue consisted in adding new compartments, *e*.*g*. for exposed people but not yet infectious, thereby increasing the heterogeneity in the infected period. Earlier studies indeed show that the addition of intermediate compartment transform the sum of exponential distributed waiting times into an hypoexponential distribution [A. L. Lloyd, 2001]. Note that taking memory into account can also be achieved by completely different formalism such as discrete-time modelling, and thus be tailored to epidemiological data the time resolution of which is almost always that of the day [Sofonea et al., 2021].

Nevertheless, depending on the issue dealt with by the model, sources of heterogeneity might emerge and increase the number of host categories and thus the total number of equations and parameters. With the onset of vaccination campaigns, this phenomenon became even more pronounced [Kiem et al., 2021; Moore et al., 2021]. Even if this approach remained a useful approximation, the initial gain in simplicity progressively vanished, making the models increasingly difficult to analyse and parameterize.

Virus evolution and the emergence of variants of concern (VOC) [Davies, Abbott, et al., 2021], coupled with some pre-existing unknowns regarding the behaviour of natural and vaccine-induced immune responses [Zellweger et al., 2020; Alizon and Sofonea, 2021], further increased the need for modelling approaches. However, even for medium or long-term projections, ODE-based approaches remain far from ideal since immunity waning may occur rapidly (at least from a prospective point of view) and might not be memoryless.

To address these issues, we introduce an alternate formalism relying on partial differential equations (PDEs) instead of ODEs. Similarly to Richard et al. [2021], our model is structured in terms of the age of the host (in years) and the infection age (in days). The originality is that we incorporated vaccination status and clearance memory. More precisely, we also record the time since vaccination as well as the time since clearance (both in days). These daily time structures allow us to keep track of the time spent in each compartment, thereby providing a convenient way to correct the memory problem while limiting the number of epidemiological compartments. Despite sharing similarities and simplicity in its formalism, PDEs system based model needs a supplementary initial effort for parameterization, but in exchange it allows for more flexibility without adding more compartments.

We first present the details of the model and derive some of its main properties, such as the basic reproduction number. Then, we perform a sensitivity analysis to identify which parameters affect epidemiological dynamics the most. We also use the model to analyse the past French Covid-19 epidemic and formulate projections regarding future trends. Finally, we discuss perspectives to extend the model and tackle additional questions regarding Covid-19 epidemic.

## 2 Model

### 2.1 Model overview

The density of susceptible individuals of age *a* ∈ [0, *a*_{max}] at time *t* is denoted by *S*(*t, a*). Susceptible individuals leave the compartment either by being infected, at a rate *λ*(*t, a*) corresponding to the force of infection, or by becoming vaccinated, at a rate *ρ*(*t, a*).

Infected individuals are denoted *I*^{ℓ} ∈ (*t, a, i*). The exponent *ℓ* ∈ {*m, s, d*} corresponds to the three types of infections: *m* for mild and asymptomatic cases, *s* for severe cases that will require hospitalization at some point before recovering, and *d* for severe cases that always result in the patient’s death. Each of these categories of infected individuals are further stratified according to the time since infection, which is indexed by *i* ∈ [0, *i*_{max}]. Practically, this affects the recovery rates (*γ*^{m}(*a, i*) and *γ*^{s}(*a, i*)) and the death rate (*µ*(*a, i*)), as well as the different transmission rates *β*^{ℓ}(*a, i*), *ℓ* ∈ {*m, s, d*}; all of which are functions of *i*. The number of new mildly infected individuals at a given time *t* is given by the boundary condition
where *p*_{a} is the proportion of infections that lead to severe cases for individuals of age *a*.

We add a similar time structure *j* to record time since clearance for the density of recovered individuals, *R*(*t, a, j*), in order to account for a possible post-infection immunity waning at a rate *σ*(*a, j*). Recovered individuals are assumed to be vaccinated at the same rate as susceptible individuals, *ρ*(*t, a*). The number of newly recovered individuals of age *a* at time *t* is given by the boundary condition

The density of vaccinated individuals, *V* (*t, a, k*), also has its own time-structure *k* to capture the time since vaccination. This allows to take into account the immunity waning, *σ*^{v}(*a, k*), or any temporal variation in vaccine efficacy. The number of newly vaccinated individuals is given by the boundary condition

Since vaccine efficacy may be imperfect, we assume that vaccinated individuals can still be infected by the virus, but at a rate reduced by 1 − *ε*(*a, k*) compared to susceptible unvaccinated individuals. If the infection is mild, infected vaccinated hosts move to the *I*^{mv}(*t, a, i, k*) compartment, which is separated from mildinfected former susceptible individuals to allow for reduced transmission at a rate 1 − *ξ*(*a, k*). Vaccinated individuals can also develop a severe form of Covid-19 following infection but at a rate reduced by (1 − *ν*(*a, k*)). And therefore reduced by (1 − *ε*(*a, k*))(1 − *ν*(*a, k*)) compared to susceptible individuals. Hence, the number of severely infected individuals of age *a* at time *t* is given by the boundary condition
and
where ifr_{a} denotes the infection fatality rate (IFR), that is the fraction of individuals of age *a* who die from the infection. It is worth to note that due to VOC emergence inducing increase in virulence, both *p*_{a} and ifr_{a} will be scaled by *κ* accounting for this increase.

Regarding the infected vaccinated individuals who develop mild symptoms, the boundary conditions are and

The model flowchart is displayed in Figure 1. Notice that our model only has 8 compartments. Note also that vaccines can act in three non mutually exclusive ways by decreasing the risk of being infected (*ε*(*a, k*)), the probability to develop severe symptoms if infected (*ν*(*a, k*)), and the transmission rate if infected (*ξ*(*a, k*)).

The change over time, including the leaving of each compartment, is provided by the following PDE system, coupled with the boundary conditions (1)–(7):
with
for any (*t, a, i, j, k*) ∈ ℝ^{+} × [0, *a*_{max}] × [0, *i*_{max}] × [0, *j*_{max}] × [0, *k*_{max}]. Here, *K*(*a, a′*) is the kernel giving the mean contact rate between two individuals belonging respectively to the age classes *a* and *a′*. We also introduce efficacy, denoted *c*, of non-pharmaceutical interventions (NPIs) in reducing the contact rates between individuals independently of their age. We assume NPIs affect all individuals indifferently, no matter the compartment they belong to. Therefore, since both susceptibles and infected individuals are targeted, the reduction of the contact rate is a squared term.

The above system is associated with Assumption S1 in Supplementary Methods and the following initial conditions

Notice that the well-posedness of System (1)–(16) is analysed in Supplementary Methods A.2.

### 2.2 Model parametrization

In this study, we focus on the French Covid-19 epidemic in 2021. The values used are shown in Table 1 along with the (French) data we use for parameter inference. Additional details about these can be found in Supplementary Methods B.

The basic reproduction number is fixed but varies in time due to the emergence of the *α* and *δ* VOCs. The *α* VOC was first detected in France in early January and rapidly became dominant. Therefore the ℛ_{0} retained starting in January was 4.5 [Haim-Boukobza et al., 2021; Davies, Abbott, et al., 2021]. By the month of July, the *α* VOC was supplanted by the *δ* VOC, increasing the ℛ_{0} up to 6 [Alizon, Haim-Boukobza, et al., 2021].

Regarding the modelling of vaccine efficacy, for simplicity, we neglect immune waning, *i*.*e*. the decrease of immunity with time, meaning that *σ* ≡ 0 and *σ*^{v} ≡ 0. This assumption is motivated by the fact that we consider a medium-term scenario and it could readily be modified. We also assume that the three types of vaccine efficacies (against reinfection, severe symptoms, and transmission) are not maximal upon entry in the vaccinated compartment. More precisely, we assume a double sigmoid curve to capture two vaccine injections (Figure S2). The different levels of efficacy are based on the Public Health England [2021] report, and additional details are provided in the Supplementary Methods B.5. The vaccination rate *ρ*(*t, a*) is based on the observed French data (see Supplementary Methods B.6 for details about this implementation).

The different transmission rates *β*^{ℓ}(*a, i*), *ℓ* ∈ {*m, s, d*}, are simply the generation time, weighted to correct for the possibility for individuals to leave the infected compartments before the generation time becomes null.

Concerning some age-stratified parametrization functions, we assume no differences between age-groups. This assumption is either made for parsimony reasons (*i*.*e*. for *γ*^{m}(*a, i*), *γ*^{s}(*a, i*), and *µ*(*a, i*)) or because of lack of information (*e*.*g*. for *β*^{m}(*a, i*) and *β*^{s}(*a, i*)).

### 2.3 Contact matrices

Due to the available data and following the parametrization relative to the severity disease, the kernel *K*(*a, a′*) is also given for a finite number of age classes, thus providing a contact matrix. And this contact matrix *K*(*a, a′*) is also an important part that needs to be define as it will define the age-structure of the population regarding an age-severity differentiated infectious disease [Valle, Hyman, and Chitnis, 2013; Jacco Wallinga, Teunis, and Kretzschmar, 2006]. In that regard, we decide to present two competing choices. The first one from Béraud et al. [2015] was estimated to better apprehend the spread of infectious diseases. The second source of contact matrices comes from the French health agency (Santé Publique France) and the French national health insurance (CNAM). They provide Covid-19 specific contact matrices for 38 weeks ranging from August, 2020 to April, 2021.

The latter reveal pronounced changes across weeks. These are most likely due to a variety of reasons such as control restrictions policies or school holidays. Interestingly, these changes do not affect all age classes in the same way (Figure 2).

The two sources of contact matrix also exhibit qualitative pattern differences, as illustrated in Figure 3. Indeed, that from Béraud et al. [2015] gives more weight to relatively young people who tend to have contact with people close in age, such as colleagues or friends, which could be representing the active population. Furthermore, in this matrix, older people have few contacts. Conversely, SPF matrices seem to have more extra-generational contacts, which could increase the role of transmission within households.

Therefore, we included all (normalized) contact matrices in the sensitivity analysis.

### 2.4 Model outputs and fitting procedures

The main model outputs are population sizes of the compartments over time for the year 2021 in France. The French publicly available hospital admission data is not stratified by age so we only use the global incidence data for parameterization and comparison purposes. In our model, incidence data in hospital admissions dynamics corresponds to the entry in the severe infection compartments (*I*^{s}(*t, a, i*) and *I*^{d}(*t, a, i*)) with a twelve days lag [Salje et al., 2020].

For each compartment dynamic, we build a 95% confidence interval using the 0.025 and 0.975 quantiles at each time step of all model runs used for the sensitivity analysis (see below).

Regarding parameter inference, we consider daily minimal sum of squares between the data and simulations. We first fit the vaccination rate *ρ*(*t, a*) on French data as detailed in Supplementary Methods B.6. Due to high computational cost, we fit the NPI policies efficacy only on the median trajectory (defined as the trajectory obtained using the median parameter set but with the Béraud et al. [2015] contact matrix).

### 2.5 Sensitivity analysis

We perform a variance-based sensitivity analysis to assess the robustness of the model given its inputs. We compute the Sobol main sensitivity indices for each model parameter and for each time step [Saltelli et al., 2008]. For an input parameter *X*_{i} and a given day, this index reflects the fraction of the variance in the output *Y* (here the daily hospital admissions) and is defined by

The difference between the sum of the main indices and 1 corresponds to the variance originating from the interactions between all the parameters. The analysis was performed on 30, 400 model runs with different parameters sets chosen using a Latin Hypercube Sampling within the ranges detailed in Table S1.

Assessing the sensitivity of model outputs depending on the contact matrix is more delicate since drawing each matrix coefficient would be numerically too costly and drawing an entire matrix would cause a loss of information regarding the role of the different age classes. However, we possess 38 weekly contact matrices from SPF and another contact matrix from Béraud et al. [2015]. Therefore, for each age class, we randomly draw the corresponding age class column (*i*.*e*. the rate of being infected for the given age class by all the age classes) among the 39 available matrices. As discussed in Section 2.3, the two sources of contact matrices exhibit qualitatively different patterns, suggesting potential difference in terms of within-household transmission or active population transmission patterns. To avoid giving more weight to a specific pattern, the Béraud et al. [2015] matrix was weighted 38 times more than the SPF matrix.

Additional details regarding the sensitivity analysis can be found in Supplementary Results C.

## 3 Results

### 3.1 Basic reproduction number and NPIs

The basic reproduction number, denoted ℛ_{0}, is a widely used metrics in epidemiology because it corresponds to the average number of secondary infections caused by an infected host in an otherwise fully susceptible population [Anderson and Robert M May, 1992]. Calculating it for our PDE system is not trivial and for this we use the next generation operator approach [Diekmann, Heesterbeek, and Metz, 1990; Inaba, 2012]. More precisely, we show that the number of new infections in individuals of age *a* at time *t* in a fully susceptible population, denoted *I*_{N} (*t, a*), satisfies the renewal equation (see Supplementary Methods A.3 for details)
where Ω(*a, i*) can be interpreted as the infectiousness expectation of an individual of age *a* infected since time *i* and is defined by
where *π*_{ℓ} is the survival probability (*i*.*e*. remaining in the compartment) of infected individuals of the *I*^{ℓ} compartment. Mathematically, , with *f*_{ℓ} = *γ*^{m}, *γ*^{s}, *µ*, respectively, for *ℓ* = *m, s, d*.

Following the Next Generation Theorem, the basic reproduction number ℛ_{0} is calculated as the spectral radius, noted *r*(*U*), of the next generation operator *U* defined from into itself by

For parametrization purpose, we assume that the contact matrix *K*(·, *·*) is given up to a positive constant *β* (to be determined), such that , and satisfies . Consequently, we find that *β* is given by
where *Ū* is the operator defined from into itself by
with

In the following, the ℛ_{0} is set to correspond to that of the *α* and then the *δ* VOC, which are both higher than that of the initial lineages. Note that, within this paper, we scale *K*(·, ·) by
rather than *β* because we estimate the level of NPI efficacy (*c*) beforehand on real data, and for prospective scenarios after the current date, we arbitrarily set it to the desired value.

### 3.2 Sensitivity analysis

Performing per time-point sensitivity analyses on the daily hospital admissions for all the model parameters (Figure 4), we noticed that most of the variance originated from the contact matrix (and its 9 parameters), especially the younger age classes. This effect was even more striking when considering the raw variance originating from each parameter (Figure S3).

We also observed important time-variations of some parameters, such as the generation time Weibull’ scale parameter. The time period where this is the most predominant also corresponds to the period with few newly hospital admissions (Figure 6), and therefore a lesser variance. The sharp decrease of this parameter sensitivity coincides with the epidemic’s growth reprisal.

Furthermore others parameters explained variance, such as the VOC-increased virulence (in red) at first or more notably the interactions between parameters progressively increasing over time reaching half of the variance explained at the end of the year.

### 3.3 Inferred dynamics

By parameterizing our model with existing data and inferring additional parameters, we could estimate past epidemic dynamics and investigate scenarios for future trends (Figure 5). The vaccinal coverage modeling did follow quite well real data, even though dissimilarities emerged in the summer (corresponding to the French summer holidays).

We may also observe a slight rebound in infected individuals mid-March, which follows by two weeks the end of the winter holidays. We also see this phenomenon on the new hospital admissions (Figure 6), with a supplementary delay corresponding to the delay between infection and hospital admissions.

It seems that the current vaccination rate is not high enough to avoid a new epidemic wave. However, we observe that the uncertainty is huge and does not allow to have a precise idea of what might happen.

## 4 Discussion

Mathematical modelling has emerged as a central tool to control and anticipate the SARS-Cov-2 pandemic. This importance is likely to increase now that vaccination has become the cornerstone of the public health response. However, limitation of current vaccination models lies in either neglecting memory effects or compensating by highly dimensional models with dozens of ordinary differential equations. In this study, we used partial differential equations to develop a model than can capture medium and long term hospital admission dynamics in a population with natural and vaccine-induced immunity only with 8 general compartments.

Regarding the sensitivity analysis, the contact matrix unexpectedly contributed more variance in daily hospital admissions than the VOC related increase of virulence itself. This predominant role is somehow surprising because although there are transmission and sensitivity differences based on age (e.g. [Davies, Klepac, et al., 2020]), the strongest age differences appear in the IFR. And more precisely, contact to younger age groups appeared to be the most important contributor to outcome variance, although they were, and by far, the less likely to be hospitalized.

An important limitation of the model is that the contact matrix is assumed not to vary over the course of a simulated epidemic. As suggested by the temporal variance in the SPF matrix data (Figure 2), this may be oversimplistic. For instance, we observed a difference of patterns in simulations whether they assumed an high or low contact rates among younger age-classes (as shown on Supplementary Figure S4). A baseline for the different contacts rates, if such a concept can even exist biologically, would most likely be impossible to determine because of the variety of events over a year inducing changes in social interactions such as calendar events (*e*.*g*. school holidays), implementation of control policies (*e*.*g*. lockdown, curfew), or even mediatization of the epidemic (resulting in spontaneous behavioural change with respect to the perception of the epidemic).

The importance of the age-structure of the host population in shaping Covid-19 epidemics is widely acknowledged. However, this effect is usually studied in the clinical context of disease severity and less so for transmission dynamics [Salje et al., 2020; Sofonea et al., 2021]. However, using a PDE formalism, Richard et al. [2021] find the population structure to be the parameter that brings the relatively the most variance to their model’s output. Both Richard et al. [2021] and Keeling, Hill, et al. [2021] use a constant contact matrix, but they explore the impact of age-differentiated NPI policies. That being said, none of those studies (including ours), are able to fully assess the role of the age-structure since there is potentially additional unknown patterns impacting medium-term forecasting. For example, in absence of external data it seems impossible to distinguish between the “true contact matrix” and age-differentiated NPI policies. The problem is that the two would likely yield different outputs in NPI-lifting scenarios.

The variance brought by the other parameters is quite low, even if we can mention some such as the increase in virulence associated with infections caused by VOCs. The originality of our approach is that it shows that this effect represents a large fraction (nearly 20%) of the observed variance in hospital admissions for a wide range of parameters. Note that the relative importance of virulence is less in the prospective part of the model (*i*.*e*. after August 2021) but this is potentially because other parameters such as the ones related to vaccination and interactions became more important over time. We also did not consider an increase in virulence in the *δ* VOC compared to the *α* VOC (but recent data shows this might very well be the case [Sheikh et al., 2021]).

The generation time Weibull’s scale parameter, which has an impact on the mean generation time, also affects hospital admission dynamics, especially in June/July and in November/December. However, we shall put it in perspective since the impact of this parameter arise only when there is few variance (Figure S3) and a decreasing epidemic (Figure 6). This can be explained by the fact that a shorter mean generation time for a given reproduction number is known to increase the epidemic’s growth rate [Nishiura, 2010; Wallinga and Lipsitch, 2007]. On this aspect, it might be relevant to note that the *δ* VOC seems to have a shorter generation time than the wildtype strain which was not taken into account [Zhang et al., 2021].

By applying our model to the context of the French epidemic, we show that vaccination is unlikely to be able to prevent a new epidemic on its own without the use of non-pharmaceutical interventions (NPIs). A strong caveat is that anticipating variations in the vaccination rate is extremely difficult as it relies on sociological factors. Increasing the uncertainty range for this parameter would most likely increase the variance in the sensitivity analysis. However, given past vaccination dynamics, we do not expect this to qualitatively affect the results. Regarding the medium-term forecasting, we did not include potential weather-related variations in behaviour or infectivity, which have been reported for SARS-CoV-2 [Ma et al., 2021].

To analyse our model, we had to make a number of simplifying assumptions, which are common to differential equation-based models. The two major ones are the lack of spatial heterogeneity and the contact homogeneity among a given age-class. The lack of spatial heterogeneity implies an identical contact rate across the whole country. This is not problematic at the start of an epidemic but is not adapted for long-term modelling as it affects persistence of the disease [Alun L. Lloyd and Robert M. May, 1996; Hagenaars, Donnelly, and Ferguson, 2004]. Furthermore, age contact patterns allow us to capture some of the heterogeneity in the population but there could be other social heterogeneities that could, for instance, correlate with vaccination status. As shown in the case of influenza virus, these could affect epidemiological dynamics [Barclay et al., 2014].

One advantage of this PDE model is the restrained number of compartments. An alternative to the problem for ODE system based models would be to chain and multiply compartments. This would also require rewriting the formalism, depending on whether we consider short, medium or long-term temporal scale. Here, this problem can be alleviated by PDE models, such as the one presented in our study due to the presence of memory structures. The same model could be used to monitor new hospital admissions or the need of a new vaccination campaign years later if it happen to have time-induced loss of immunity.

Another advantage of the PDE formalism is that our model can be extended in several ways without adding any compartments thanks to the structuring in time since an event (infection, recovery, or vaccination). For instance, we could account for a difference in the building up of vaccine immunity in susceptible versus recovered individuals. This would be consistent with the fact that the latter enter the vaccinated compartment at a later vaccination age (*i*.*e. k >* 0) and a single vaccine dose appears to be sufficient to build strong immunity [Mazzoni et al., 2021]. Also, in the context of waning immune, our model can be used to investigate the long-term benefits or costs of implementing vaccine boosters depending on assumptions regarding vaccine efficacy or duration of natural immunity. More generally, we can investigate the effect implementing age-stratified vaccination policies.

One downside for PDE system based model we can think of is an higher computation time.

On a more prospective side, such a model could be used to investigate virus evolution with explicit interplay between change in susceptibility, contagiousness, virulence and immune escapes (post-infection and vaccine) and trade-offs between them.

## Data Availability

Data and scripts are available in the linked git repository.

## Conflict of interest

We declare no conflict of interest.

## Funding

BR is funded by the Ministère de l’Enseignement Supérieur et de la Recherche.

## Data, script and code availability

The model was implemented in `R` (v4.1.1) [R Core Team, 2021]. The data and scripts used within this study are available in the repository https://gitlab.in2p3.fr/ete/covid19_vaccination_model.

## Supplementary information

### Supplementary figures and table

### A Model

System (1)–(16) is considered under the following general assumption

*p*_{a}∈ [0, 1],*and*ifr_{a}∈ [0,*p*_{a}]*for all a*∈ ℝ_{+};,

*with*0 ≤*ρ*(*t, a*) ≤ 1,*for all*(*t, a*) ∈ ℝ_{+}× ℝ_{+};*with*0 ≤*ε*(·, ·),*ξ*(·, ·),*ν*(·, ·) ≤ 1;;

;

*Transmission rates satisfy**for each ℓ*∈ {*m, s, mv*}.

#### A.1 Implementation

The model was implemented in `R` [R Core Team, 2021], using `Rcpp` [Eddelbuettel and François, 2011] to maximize computational efficiency.

The PDE system was implemented using an Euler explicit scheme.

#### A.2 Well-posedness

Let us introduce the Banach space as well as its positive cone

Now, we define the subspaces by:

It follows that there exists a unique linear operator for each *ℓ* ∈ {1, 2} such that Π_{1}*ψ* = *ψ*(·, 0, *·*) and Π_{2}*ψ* = *ψ*(·, *·*, 0) for all *ψ* ∈ *Y*_{2}. Next, let *A* : *D*(*A*) ⊂ *X* →*X* be the linear operator on the domain
and defined by

Finally, let us introduce the following nonlinear map :
wherein *ϕ*(*t*) is the function:

From here, System (1)–(16) rewrites as the following nondensely defined Cauchy problem:

Therefore, under Assumption S1, we have the well-posedness of System (S1); that is, the Cauchy problem (S1) generates a unique globally defined, positive and bounded non-autonomous semiflow.

The proof of this result is based on a rather standard methodology combining an integrated semigroup approach and Volterra integral formulation in the context of multiple structured variables (*e*.*g*., Richard, Choisy, et al. [2022] and the references therein) and existence of the semiflow for non-autonomous systems (*e*.*g*., [Pazy, 2012; Magal, 2001]).

#### A.3 Basic reproduction number and NPI policices

We consider that there are no vaccinated, *i*.*e ρ*(*t, a*) ≡ 0 and *V* (*t, a, k*) ≡ 0, nor recovered people, *i*.*e. R*(*t, a, j*) ≡ 0, and that the number of susceptible individuals is very close to the total population size. For simplicity, we first introduce the “survival” probability (*i*.*e*. the probability to remain in the compartment) of infected individuals of, respectively, the *I*^{m}, *I*^{s}, and *I*^{d} compartments:

By linearizing System (1)–(16), we obtain the following Volterra formulation for *I*^{m}, *I*^{s}, and *I*^{d} compartments:
with *λ*_{0}(*t, a*) defined as *λ*(*t, a*) with no control policies (*c* = 0),

Let *I*_{N} (*t, a*) = *λ*_{0}(*t, a*)*S*(0, *a*) be the density of newly infected of age *a* at time *t*. Then, by (S5)–(S7) it comes
where
and *f* (*t, a*) is accounting for the initial population.

The basic reproduction number ℛ_{0} is then the spectral radius, denoted by *r*(*U*), of the next generation operator *U* defined from into itself by

For parametrization purpose, we assume that the contact matrix *K*(·, *·*) is given up to a positive constant *β* (to be determined), such that , and satisfies . Consequently, we find that *β* is given by
where *Ū* is the operator defined from into itself by
with

Note that, within this paper, we scale *K*(·, ·) by
rather than *β* since the NPI level efficacy was fitted beforehand on real data.

To go further steps in the computation of *r*(*Ū*), in addition to the general Assumption S1, we also assume that

*Functions S*_{0}, *K*, *are positive almost everywhere*.

Then, we can show that *r*(*Ū*) is given by the spectral radius of the following linear operator, defined from into itself:

The spectral radius of this later operator is computed easily since the age *a* is numerically divided into *n* ∈ ℕ^{∗} classes, so that the term inside the integral of the latter equation is a *n* × *n* matrix. Finally, the scaling parameter *β* is obtain from (S9).

Importantly, the symmetric property of the contact matrix *K* is not strictly necessary for the the computation of *r*(*Ū*). However, in addition to Assumptions S1 and S2, if *K* is a symmetric function, then the Rayleigh quotient formulation leads to (see Proposition F.2 in Richard, Alizon, et al., 2021)

### B Model parametrization

In this section, we describe the parametrization and the assumptions made in the main text. The uncertainty ranges retained for each parameter are displayed in the Table S1.

#### B.1 Proportion of severe cases, IFR and increase in virulence

The proportion of severe cases corresponds here to the fraction of the population who will be hospitalized following a SARS-CoV-2 infection. This parameter is age-dependent and follows the infection fatality rate (IFR) by Verity et al. [2020].

However, studies show that the virulence of the infection increased (taken into account by the *κ* parameter) by more than 60% with the *α* VOC [Davies et al., 2021; Challen et al., 2021]. Estimation for the *δ* VOC are still at an early stage, so we used the virulence from the *α* VOC, even if the former seems more virulent [Sheikh et al., 2021].

#### B.2 Generation time and transmission rates

For the underlying generation time, we use that provided by Ferretti et al. [2020], which follows a Weibull distribution, with a shape parameter of 2.826 (95% CI [1.75 − 4.7]) and scale parameter of 5.665 (95% CI [4.7 − 6.9]). We assume that the transmission rates from mild infections *β*^{m}(*a, i*) and severe cases *β*^{s+d}(*a, i*) are equal to the generation time corrected by the probability for individuals to leave their infected compartment.

#### B.3 Recovery rates

From the data shown in Salje et al. [2020], we retrieve the time severely-infected individuals spend as infected, whether they required ICU admission or not, by adding up the different exponential distributions of the different infected compartments of their model (which will be denoted *E*_{1}, *E*_{2}, *I*^{hosp}, *I*^{non ICU}, *I*^{ICU}, *H*_{1}, *H*_{2}, *H*_{ICU}, *ICU* _{1}, *ICU* _{2} in a later study by Kiem et al. [2021]). This gives us the probability of remaining in the infected compartment over time, thereby allowing us to infer the recovery rates. We also know from Lefrancq et al. [2021] the probability for hospitalized individuals to require ICU admission, which provides us with appropriate recovery rates weighting for severe cases. We apply the same method for mildly-infected individuals.

#### B.4 Initial conditions

For this PDE model, the initialization is not as straightforward as for ODE models since within a compartment individuals do not have the same age of infection (for infected individuals) or time since clearance (for recovered individuals). Initializing over all the domain of definition of each compartments is difficult since a uniform initialization would almost immediately be counterbalanced by a higher probability to leave the compartment for higher ages (*i*.*e. i >* 0 and *j >* 0). Put differently, this would produce distributions different from what we might expect with a constant inflow in the compartments. To overcome this issue, we start the different runs with a 45 days delay (not shown) in order to let the different compartments stabilize around some distribution.

Hozé et al. [2021] estimate that the proportion of recovered adults in France was of 0.149 (95% CI [0.132 − 0.169]) on January 15th, 2021. For simplicity, and in absence of more detailed data, we assumed this proportion to be constant across age classes, including the younger age groups.

For infected individuals, we initialize the density with a qualitative value. The *β* coefficient associated to the ℛ_{0} and the NPI was fitted such that the number of daily hospital admissions was close to the real data. The overall number of infected individuals in the model at January 1st, was around 267, 000 [170, 000 − 310, 000].

#### B.5 Vaccine properties

We assume that the three types of vaccine efficacies (against infection, severe forms, and transmission) follow a double-sigmoid temporal pattern starting from the day of injections (Figure S2). The reduction in transmission rate corresponds to function *ξ*(*a, k*) in our model, the infection immunity corresponds to function *ε*(*a, k*). The total reduction of virulence corresponds to the cumulative effects of *ε*(*a, k*) and *ν*(*a, k*).

The order of magnitude of the final (full) efficacy levels are based on that from the Pfizer-BioNTech vaccine after two doses provided by Public Health England [2021]. However, note that the outcomes referenced in the report were used as overall proxy in this study (“symptomatic disease” and “infection” for *ε*(*a, k*), “hospitalisation” and “mortality” for *ν*(*a, k*) and “transmission” for *ξ*(*a, k*)) as they do not exactly match our implementation. The reduction of transmission estimation of 14 days after the second dose was not available, so we applied a rule of thumb so that the increase between the two injections was similar to the other vaccine properties. Hence, we assume a 75% efficacy for transmission reduction 14 days after the second dose.

We also assume there was no difference in efficacy between age classes, and that the different efficacy levels remains constant 14 days after the second dose until the end of our projections (*i*.*e*. no immune waning).

#### B.6 Vaccination rate

We model the vaccination rate using a sigmoid function,
where *v*_{tot} denotes the total number of vaccinated individuals at the end of the year (which is a model input), and *θ*_{1} and *θ*_{2} are the sigmoid curve parameters fitted to the observed data. This gives us the number of newly vaccinated people at each time step.

The number of doses attributed to each age-group at each time step depends on initial weights (*ω*_{a}(*t*_{0})), which can be interpreted as the age-based strategy vaccination prioritization, the proportion of the age group targeted (*t*_{a}), assuming that the total number of vaccinated in each age class may vary and is lower than 100%, and the proportion of individuals already vaccinated within each age-group at time *t* (*p*_{v}(*t, a*)).

Therefore, at time *t* + Δ*t*, the splitting of the number of doses is given by *ω*_{a}(*t* + Δ*t*) which is defined by

Finally, we have

The initial weights, *ω*_{a}(*t*_{0}), and *t*_{a} is fitted on a ordinary least squares metric to reproduce at best the real vaccination rate.

Vaccination is assumed to start on January, 1st, 2021.

### C Sensitivity analysis

To perform the sensitivity analysis, we use the `lhs` package [Carnell, 2020] to generate a Latin Hypercube Sample (LHS). The parameters were drawn in a uniform distribution within the confidence interval specified for each parameters and shown in Table S1.

For each parameter combination, a model run was computed. In total, 30, 400 model runs were performed. Then, for each time step, we used the `multisensi` package [Bidot, Lamboni, and Monod, 2018] to compute the Sobol main indices, given by
as implemented in the `sensitivity` package [Iooss et al., 2021]. The difference between the sum of all the main indices and 1 corresponds to the effect of interactions between parameters. More explanations are available in Saltelli et al. [2008].

Due to numerical approximations, some indices may sometimes be negative (the lowest was − 0.004). These were rounded to 0.

## Acknowledgments

The authors acknowledge the ISO 9001 certified IRD i-Trop HPC (South Green Platform) at IRD montpellier for providing HPC resources that have contributed to the research results reported within this paper.

We thank all the ETE modeling team for thoughtful discussions.

We thank the French public health agency (Santé Publique France) and the French national health insurance (Caisse Nationale d’Assurance Maladie) for providing Covid-19 contact matrices, and especially Alexandra Mailles and Jonathan Bastard.