key: cord-0300822-tm70gwe6
authors: Reyne, B.; Richard, Q.; Nous, C.; Selinger, C.; Sofonea, M. T.; Djidjou-Demasse, R.; Alizon, S.
title: The importance of the population age-structure: insights from Covid-19 dynamics model structured by age, time since infection and acquired immunity
date: 2021-10-03
journal: nan
DOI: 10.1101/2021.09.30.21264339
sha: efd11adeeddafb02feb1f16b240792fdb01fcb91
doc_id: 300822
cord_uid: tm70gwe6

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.

Shortly after the Covid-outbreak late , many e orts were put in diverse research areas to understand both the disease and its aetiological agent, SARS-CoV-, 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 rst-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., ; 5 Hellewell et al., ; Beneteau et al., ], 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, ]. Many of the deterministic models developed were variations of the classical Susceptible-Infected-Recovered (SIR) compartmental model, usually involving a system of ordinary di erential equations (ODEs) 10 [Keeling and Rohani, ] . Such a simple formalism was adapted at rst 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 compart-15 ment, as implied by the ODE formalism-became particularly detrimental to short-term forecasting (see Supplementary Figure S ). For example, for the time elapsed between infection and potential hospital admission, this memoryless hypothesis does not hold [Salje et al., ; Sofonea et al., ] . 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 inter-20 mediate compartment transform the sum of exponential distributed waiting times into an hypoexponential distribution [A. L. Lloyd, ] . Note that taking memory into account can also be achieved by completely di erent 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., ]. Nevertheless, depending on the issue dealt with by the model, sources of heterogeneity might emerge 25 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., ; Moore et al., ]. Even if this approach remained a useful approximation, the initial gain in simplicity progressively vanished, making the models increasingly di cult to analyse and parameterize.

Virus evolution and the emergence of variants of concern (VOC) [Davies, Abbott, et al., ] , coupled 30 with some pre-existing unknowns regarding the behaviour of natural and vaccine-induced immune responses [Zellweger et al., ; Alizon and Sofonea, ] , 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 di erential equations

. Model overview 50 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 55 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 strati ed according to the time since infection, which is indexed by i ∈ [0, i max ]. Practically, this a ects the recovery rates (γ m (a, i) and γ s (a, i)) and the death rate (µ(a, i)), as well as the di erent 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 60 I m (t, a, 0) = (1 − p a )λ(t, a)S(t, a), ( ) 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 65 R(t, a, 0) = imax 0 γ s (a, i)I s (t, a, i) + γ m (a, i)I m (t, a, i) di.

( )

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 e cacy. The number of newly vaccinated individuals is given by the boundary condition V (t, a, 0) = ρ(t, a)S(t, a) + ρ(t, a) jmax 0 R(t, a, j)dj.

( )

Since vaccine e cacy 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 70 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-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 75 I s (t, a, 0) = p a 1 − ifr a p a λ(t, a) S(t, a) + kmax 0 (1 − ε(a, k))(1 − ν(a, k))V (t, a, k)dk ( ) 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 Figure : Model owchart. On this owchart, subscripts denote additional structure beside time t for each compartment: a stands for the host age, i for time since infection, j for time since clearance and k for time since vaccination. Exponent denotes for di erent types of infections: m for mild or asymptomatic cases, s for severe cases who will survive, d for cases who will die and mv for mild infection though vaccinated host.

and

The model owchart is displayed in Figure . 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 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 system, coupled with the boundary conditions ( )-( ):

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 e cacy, denoted c, of non-pharmaceutical interventions (NPIs) in reducing the contact rates between individuals independently of their age. We assume NPIs a ect all individuals indi erently, no matter the compartment they belong to. Therefore, since both susceptibles and infected individuals are targeted, the 90 reduction of the contact rate is a squared term.

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

, for ∈ {m, s, d}, I mv (t = 0, ·, ·, ·) = I mv 0 (·, ·, ·) ∈ L ∞ is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; In this study, we focus on the French Covid-epidemic in . The values used are shown in Table  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 xed but varies in time due to the emergence of the α and δ VOCs. The α VOC was rst detected in France in early January and rapidly became dominant. Therefore the R 0 100 retained starting in January was 4.5 [Haim-Boukobza et al., ; Davies, Abbott, et al., ] . By the month of July, the α VOC was supplanted by the δ VOC, increasing the R 0 up to 6 [Alizon, Haim-Boukobza, et al., ].

Regarding the modelling of vaccine e cacy, 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 105 we consider a medium-term scenario and it could readily be modi ed. We also assume that the three types of vaccine e cacies (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 The di erent 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-strati ed parametrization functions, we assume no di erences 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)).

Due to the available data and following the parametrization relative to the severity disease, the kernel K(a, a ) is also given for a nite number of age classes, thus providing a contact matrix. And this contact matrix 120 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint Figure : Coe cients of variation of each element of the SPF/CNAM contact matrix over the weeks available. The higher the value, the greater the variability in contacts between age-classes over the di erent weeks.

K(a, a ) is also an important part that needs to be de ne as it will de ne the age-structure of the population regarding an age-severity di erentiated infectious disease [Valle, Hyman, and Chitnis, ; Jacco Wallinga, Teunis, and Kretzschmar, ] . In that regard, we decide to present two competing choices. The rst one from Béraud et al. [ ] 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 125 national health insurance (CNAM). They provide Covid-speci c contact matrices for weeks ranging from August, to April, . 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 a ect all age classes in the same way ( Figure ) .

The two sources of contact matrix also exhibit qualitative pattern di erences, as illustrated in Figure ] 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.

.

The main model outputs are population sizes of the compartments over time for the year in France. The French publicly available hospital admission data is not strati ed by age so we only use the global incidence data for parameterization and comparison purposes. In our model, incidence data in hospital admissions 140 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., ]. For each compartment dynamic, we build a 95% con dence interval using the 0.025 and 0.975 quantiles at each time step of all model runs used for the sensitivity analysis (see below). is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint ] (on the right). All matrices were normalized in order to be compared.

Regarding parameter inference, we consider daily minimal sum of squares between the data and simu- 

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

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

Assessing the sensitivity of model outputs depending on the contact matrix is more delicate since drawing each matrix coe cient would be numerically too costly and drawing an entire matrix would cause a loss of information regarding the role of the di erent age classes. However, we possess 38 weekly contact matrices from SPF and another contact matrix from Béraud et al. [ ]. 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 . , the two sources of contact matrices exhibit qualitatively di erent patterns, suggesting potential di erence in terms of within-household transmission or active population transmission patterns. To avoid giving more weight to a speci c pattern, the Béraud et al. [ ] matrix was weighted 38 times more than the SPF matrix.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10. 1101 /2021 Additional details regarding the sensitivity analysis can be found in Supplementary Results C.

The basic reproduction number, denoted R 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 sus- ]. Calculating it for our PDE system is not trivial and for this we use the next generation operator approach [Diekmann, Heesterbeek, and Metz, ; Inaba, ] . 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), satis es the renewal equation (see Supplementary Methods A. for details)

where Ω(a, i) can be interpreted as the infectiousness expectation of an individual of age a infected since time i and is de ned by

where π is the survival probability (i.e. remaining in the compartment) of infected individuals of the I compartment. Mathematically, π (a, i) = e − i 0 f (a,r)dr , with f = γ m , γ s , µ, respectively, for = m, s, d. Following the Next Generation Theorem, the basic reproduction number R 0 is calculated as the spectral 180 radius, noted r(U ), of the next generation operator U de ned from L 1

For parametrization purpose, we assume that the contact matrix K(·, ·) is given up to a positive constant β (to be determined), such that K(·, ·) = βK(·, ·), and K satis es i j K(a i , a j ) = 1. Consequently, we nd that β is given by

In the following, the R 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 e cacy (c) beforehand on real data, and for prospective scenarios after the current date, we arbitrarily set it to the desired value.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint Figure : Temporal sensitivity analysis. We represent the main Sobol indices for each time step. These indices give the relative variance explained by each parameter. There are 9 parameters associated to the contact matrix corresponding to the rate of being infected for each age class from younger to older (bottom to top). 'Virulence' corresponds to the (increased) virulence of the VOC.

.

Performing per time-point sensitivity analyses on the daily hospital admissions for all the model parameters ( Figure ) , we noticed that most of the variance originated from the contact matrix (and its 9 parameters), especially the younger age classes. This e ect was even more striking when considering the raw variance originating from each parameter ( Figure S ). We also observed important time-variations of some parameters, such as the generation time Weibull' 195 scale parameter. The time period where this is the most predominant also corresponds to the period with few newly hospital admissions ( Figure ) , 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 rst or more notably the interactions between parameters progressively increasing over time reaching half of the 200 variance explained at the end of the year.

By parameterizing our model with existing data and inferring additional parameters, we could estimate past epidemic dynamics and investigate scenarios for future trends (Figure ) . The vaccinal coverage modeling did follow quite well real data, even though dissimilarities emerged in the summer (corresponding to the French 205 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 ) , 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, 210 we observe that the uncertainty is huge and does not allow to have a precise idea of what might happen.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint Figure : Epidemic dynamics of the French SARS-CoV-epidemic in . For vaccinated hosts, we show the real data (green dots). Vertical lines indicate imposed changes in the basic reproduction number R 0 (in red) or estimated changes in the e cacy of NPIs (in blue). Plain lines and shaded areas respectively represent the median and 95% con dence interval computed from the simulations used for the sensitivity analysis. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10.1101/2021.09.30.21264339 doi: medRxiv preprint

Mathematical modelling has emerged as a central tool to control and anticipate the SARS-Cov-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 e ects or com-215 pensating by highly dimensional models with dozens of ordinary di erential equations. In this study, we used partial di erential 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 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 220 surprising because although there are transmission and sensitivity di erences based on age (e.g. [Davies, Klepac, et al., ]), the strongest age di erences 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 225 of a simulated epidemic. As suggested by the temporal variance in the SPF matrix data ( Figure ) , this may be oversimplistic. For instance, we observed a di erence of patterns in simulations whether they assumed an high or low contact rates among younger age-classes (as shown on Supplementary Figure S ) . A baseline for the di erent 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 230 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-epidemics is widely acknowledged. However, this e ect is usually studied in the clinical context of disease severity and less so for ] use a constant contact matrix, but they explore the impact of age-di erentiated 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 240 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-di erentiated NPI policies. The problem is that the two would likely yield di erent 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 245 shows that this e ect represents a large fraction (nearly %) 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 ) 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 250 [Sheikh et al., ] ). The generation time Weibull's scale parameter, which has an impact on the mean generation time, also a ects 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 S ) and a decreasing epidemic (Figure ) . This can be explained by the fact that a shorter mean generation time for 255 a given reproduction number is known to increase the epidemic's growth rate [Nishiura, ; Wallinga and Lipsitch, ] . 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., ] . 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 260 strong caveat is that anticipating variations in the vaccination rate is extremely di cult as it relies on sociolog-. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10.1101/2021.09.30.21264339 doi: medRxiv preprint ical 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 a ect 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- [Ma et al., ].

To analyse our model, we had to make a number of simplifying assumptions, which are common to di erential 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 a ects persistence of the disease [Alun L. Lloyd and Robert M. May, ; Hagenaars, Don-270 nelly, and Ferguson, ]. 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 in uenza virus, these could a ect epidemiological dynamics [Barclay et al., ]. One advantage of this PDE model is the restrained number of compartments. An alternative to the 275 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 di erence 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 su cient to build strong 285 immunity [Mazzoni et al., ] . Also, in the context of waning immune, our model can be used to investigate the long-term bene ts or costs of implementing vaccine boosters depending on assumptions regarding vaccine e cacy or duration of natural immunity. More generally, we can investigate the e ect implementing age-strati ed 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-o s between them.

Data, script and code availability 305 The model was implemented in R (v . . ) [R Core Team, ] . The data and scripts used within this study are available in the repository https://gitlab.in2p3.fr/ete/covid19_vaccination_model. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint ). "Age-dependent e ects in the transmission and control of COVID-epidemics". en. In: Nature Medicine . , pp.

-.

: -X. : 10. 1038/s41591-020-0962-9.

: https://www.nature.com/articles/s41591-020-0962- is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint ). "Feasibility of controlling COVID-outbreaks by isolation of cases and contacts". English. In: The Lancet Global Health . . Publisher: Elsevier, e -e .

: -X. ). "On a new perspective of the basic reproduction number in heterogeneous environments". en. In: Journal of Mathematical Biology . , pp.

-. : -.

: 10 . 1007/s00285-011-0463-z.

: . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint ). "Vaccination and non-pharmaceutical interventions for COVID-: a mathematical modelling study". en. In: The Lancet Infectious Diseases, S . : .

: 10 . 1016 / S1473 -3099(21 ) 00143 -2.

: https : / / linkinghub . elsevier . com / retrieve/pii/S1473309921001432 (visited on / / ). Nishiura, Hiroshi ( ). "Time variations in the generation time of an infectious disease: Implications for sampling to appropriately quantify transmission potential". en. In: Mathematical Biosciences and Engineering . , pp.

- is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 

Without memory Providing memory to any compartment might yield di erent results in terms of epidemic dynamics. As illustrated here, in both cases half of the people leave the compartment at approximately twelve days. However, the number of people still in the compartment at or days might be radically di erent, given the initial population, whether memory is provided or not. By assumption . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; Figure S : Younger age classes contact rates e ect on hospital admission dynamics. We show randomly selected trajectories among runs made with younger age-classes having the lowest (resp. highest) contact rates. We selected 50 trajectories among runs made with the 3 contacts matrices having the lowest (resp. highest) contact rates among [0 − 9] y.o., plus 50 for the [10 − 19] y.o.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10. 1101 /2021 System ( )-( ) is considered under the following general assumption

. p a ∈ [0, 1], and ifr a ∈ [0, p a ] for all a ∈ R + ;

. ρ ∈ L ∞ + (R + × R + ), with 0 ≤ ρ(t, a) ≤ 1, for all (t, a) ∈ R + × R + ;

. ε(·, ·), ξ(·, ·), ν(·, ·) ∈ L ∞ + (R + × R + ) with 0 ≤ ε(·, ·), ξ(·, ·), ν(·, ·) ≤ 1;

. σ(·, ·), µ(·, ·), γ m (·, ·), γ s (·, ·), γ mv (·, ·), ω(·, ·) ∈ L ∞ + (R + × R + );

. K ∈ L ∞ + (R + × R + );

. Transmission rates satisfy β ∈ L ∞ + (R + × R + ) for each ∈ {m, s, mv}.

The model was implemented in R [R Core Team, ], using Rcpp [Eddelbuettel and François, ] to maximize computational e ciency.

The PDE system was implemented using an Euler explicit scheme.

Let us introduce the Banach space

as well as its positive cone

It follows that there exists a unique linear operator Π ∈ L(Y 2 , L 1 (R 2 + )) 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

CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10.1101/2021.09.30.21264339 doi: medRxiv preprint and de ned by

Finally, let us introduce the following nonlinear map F : D(A) → X:

wherein φ(t) is the function: φ(t) = S(t, ·), 0, I m (t, ·, ·), I s (t, ·, ·), I d (t, ·, ·), 0, I mv (t, ·, ·, ·), 0, R(t, ·, ·), 0, V (t, ·, ·) ∈ D(A).

From here, System ( )-( ) rewrites as the following nondensely de ned Cauchy problem:

Therefore, under Assumption S , we have the well-posedness of System (S ); that is, the Cauchy problem (S ) generates a unique globally de ned, positive and bounded non-autonomous semi ow.

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. [ ] and the references therein) and existence of the semi ow for non-autonomous systems 505 (e.g., [Pazy, ; Magal, ] ).

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 rst introduce the "survival" probability (i.e. the probability to remain in the compartment)

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 of infected individuals of, respectively, the I m , I s , and I d compartments:

By linearizing System ( )-( ), we obtain the following Volterra formulation for I m , I s , and I d compartments:

with λ 0 (t, a) de ned 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 (S )-(S ) it 

where Ω(a, i) = β m (a, i)(1 − p a )π m (i) + β s (a, i)p a 1 − ifr a p a π s (i) + β d (a, i) ifr a π d (i), and f (t, a) is accounting for the initial population. The basic reproduction number R 0 is then the spectral radius, denoted by r(U ), of the next generation operator U de ned from L 1

For parametrization purpose, we assume that the contact matrix K(·, ·) is given up to a positive constant β (to be determined), such that K(·, ·) = βK(·, ·), and K satis es i j K(a i , a j ) = 1. Consequently, we nd that β is given by

where U is the operator de ned from L 1 + ([0, a max ], R) into itself by is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 with

Note that, within this paper, we scale K(·, ·) by 520 β * := (1 − c) 2 β = (1 − c) 2 R 0 r U rather than β since the NPI level e cacy was tted beforehand on real data.

To go further steps in the computation of r(U ), in addition to the general Assumption S , we also assume that Assumption S Functions S 0 , K, Ω are positive almost everywhere.

Then, we can show that r(U ) is given by the spectral radius of the following linear operator, de ned from

The spectral radius of this later operator is computed easily since the age a is numerically divided into n ∈ 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 (S ). Importantly, the symmetric property of the contact matrix K is not strictly necessary for the the com- 

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

The proportion of severe cases corresponds here to the fraction of the population who will be hospitalized following a SARS-CoV-infection. . 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.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10. 1101 /2021 From the data shown in Salje et al. [ ], we retrieve the time severely-infected individuals spend as infected, whether they required ICU admission or not, by adding up the di erent exponential distributions of the di erent 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. [ ]) . 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 555 Lefrancq et al. [ ] 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. . 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 β coe cient associated to the R 0 and the NPI was tted 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 st, was around 267, 000 [170, 000 − 310, 000]. 575 We assume that the three types of vaccine e cacies (against infection, severe forms, and transmission) follow a double-sigmoid temporal pattern starting from the day of injections (Figure S ). 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 e ects of ε(a, k) and ν(a, k).

The order of magnitude of the nal (full) e cacy levels are based on that from the P zer-BioNTech 580 vaccine after two doses provided by Public Health England [ ]. 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 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

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 tted 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 )), 595 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 %, 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 de ned by 600 ω a (t + ∆t) = a (t + ∆t) a a (t + ∆t) a) ).

Finally, we have ρ(t, a) · [S(t, a) + R(t, a)] = ω a (t) · f (t; θ 1 , θ 2 , v tot ).

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

Vaccination is assumed to start on January, st, .

To perform the sensitivity analysis, we use the lhs package [Carnell, ] to generate a Latin Hypercube Sample (LHS). The parameters were drawn in a uniform distribution within the con dence interval speci ed for each parameters and shown in Table S . 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, ] 610 to compute the Sobol main indices, given by ]. Due to numerical approximations, some indices may sometimes be negative (the lowest was −0.004).

These were rounded to 0.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10. 1101 /2021 project.org/package=lhs. Challen, Robert et al. (Mar. ) . "Risk of mortality in patients infected with SARS-CoV-variant of concern / : matched cohort study". en. In: BMJ : https://www.R-project.org/. Richard, Quentin, Samuel Alizon, et al. (Mar. ) . "Age-structured non-pharmaceutical interventions for optimal control of COVID-epidemic". en. In: PLOS Computational Biology . . Publisher: Public Library of Science, e . : -.

: 10 . 1371 / journal . pcbi . 1008776. :

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted October 3, 2021. ; https://doi.org/10. 1101 /2021 

Rapid spread of the SARS-CoV-Delta variant in some French regions

The authors acknowledge the ISO certi ed IRD i-Trop HPC (South Green Platform) at IRD mont-295 pellier 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-contact matrices, and especially Alexandra Mailles and Jonathan Bastard.