key: cord-0859012-ef0b0035
authors: Sofonea, M. T.; Reyne, B.; Elie, B.; Djidjou-Demasse, R.; Selinger, C.; Michalakis, Y.; Alizon, S.
title: Epidemiological monitoring and control perspectives: application of a parsimonious modelling framework to the COVID-19 dynamics in France
date: 2020-05-24
journal: nan
DOI: 10.1101/2020.05.22.20110593
sha: 0a816b1d7788bc3bfc266f235c06c15f68fd99b4
doc_id: 859012
cord_uid: ef0b0035

SARS-Cov-2 virus has spread over the world creating one of the fastest pandemics ever. The absence of immunity, asymptomatic transmission, and the relatively high level of virulence of the COVID-19 infection it causes led to a massive flow of patients in intensive care units (ICU). This unprecedented situation calls for rapid and accurate mathematical models to best inform public health policies. We develop an original parsimonious model that accounts for the effect of the age of infection on the natural history of the disease. Analysing the ongoing COVID-19 in France, we estimate the value of the key epidemiological parameters, such as the basic reproduction number (R0), and the efficiency of the national control strategy. We then use our deterministic model to explore several scenarios posterior to lock-down lifting and compare the efficiency of non pharmaceutical interventions (NPI) described in the literature.

In Dec 2019, a rapidly increasing number of 'pneumonia of unknown etiology' cases in Wuhan, 25 China, triggered a specific surveillance mechanism implemented by China's Center for Disease Control (CDC) in the wake of the 2003 SARS-CoV pandemic (Li et al., 2020a) . In Jan 2020, China CDC identified the new inter-human transmitted pathogen as a coronavirus, which was later named SARS-CoV-2 after its phylogenetically close predecessor (Coronaviridae Study Group of the ICTV, 2020). Meanwhile, the respiratory disease for which it is responsible, called 30 COVID-19, was better described from the clinical point of view, especially by identifying age and co-morbidities as risk factors of symptomatic aggravation and fatality (Chen et al., 2020; Huang et al., 2020) . Because of its relatively high basic reproduction number (R 0 ≈ 2.2) in the Wuhan outbreak (Li et al., 2020a) and the high transimissibility of asymptomatic/undocumented cases (Nishiura et al., 2020a; Li et al., 2020b) , as well as the intensity of international travel, the 35 virus rapidly spread in mainland China and then all over the world.

On Mar 11, 2020, the WHO announced that the COVID-19 outbreak had reached pandemic stage. To date, no pharmaceutical treatment has proven to be efficient, neither prophylactically nor therapeutically. 40 The documented importation of SARS-CoV-2 on the French metropolitan territory dates back to Jan 24, 2020 with the detection of three cases with a travel history to Wuhan (Stoecklin et al., 2020) . These were also the first COVID-19 cases detected in Europe.

From Feb 27, the daily incidence of detected cases started to increase exponentially, followed on Mar 8 by the daily hospital mortality. On Feb 29, the government announced that France from Mar 17. 55 days later, on May 11, the lock-down was lifted.

Mathematical modelling was involved early on to estimate the magnitude of the COVID-19 epidemic in Wuhan. Estimates of the serial interval were obtained rapidly allowing to estimate the basic reproduction number (Li et al., 2020a) . Later on, analyses on the number of travelers 55 from Wuhan to Europe and infected by COVID-19 allowed to estimate the magnitude of the epidemics in Wuhan (Imai et al., 2020) or the incubation time of the infection (Backer et al., 2020) . Furthermore, stochastic models allowed to better assess the probability to control initial outbreaks (Hellewell et al., 2020) . Early models also stressed the importance of detecting infections early on to control the epidemics (Ferretti et al., 2020) . 60 It quickly became clear that the epidemic had reached a large enough size to escape stochastic forces in Wuhan and that, in spite of a full lock down and travel restrictions, it had spread all over the world, making deterministic models more appropriate. This also led to an increase in modelling efforts because deterministic compartmental models based on ordinary differential equations are more commonly used in epidemiological modelling. However, although these 65 models have useful analytical properties when analysed over a long time period, they are perform poorly on short time scales. One reason for this is that they are essentially Markovian or 'memoryless'. This means for instance that an individual that has been infected for 10 days has the same probability to clear the infection that someone infected for less than a day. Over a long time period this effect averages out but on a shorter time scale, even if the dynamics are 70 deterministic, this Markovian assumption makes it difficult to reconcile the model with data.

Focusing on the French epidemic, Di Domenico et al. (2020) developed a variant of an SEIR model with multiple levels of infection severity and estimated the R 0 to 3.0 (with a 95% CI of 2.8-3.2). Using an SIR model, Magal and Webb (2020) estimated it to 4.5 (no confidence interval computed). Using a meta-population model, Roux et al. (2020) estimated it in various agent-based simulation model with 194 parameters, which allows them to estimate R 0 to 3.1 80 (no confidence interval computed) and investigate the effect of different types of interventions on the spread of the epidemics. Finally, Salje et al. (2020) , with priority access to restricted ICU patient data, developed a detailed SEIR model with explicit ICU admissions tailored to the French epidemics. Using a detailed statistical inference model, they estimated the R 0 to 2.9 (with a 95% CI of 2.8-2.99) and an R t after the lock-down of 0.67 (with a 95% CI of 0.65-0.68). 85 The model we introduce is inspired by an earlier optimal control model (Djidjou-Demasse et al., 2020) . Taking advantage from its two key qualities, namely being simultaneously deterministic and tailored for non-exponentially distributed transition times, accurately follow the main properties of COVID-19 epidemic, while limiting the number of compartments and parameters.

In this study, by using the French COVID-19 epidemics as an example, we introduce a more 90 mechanistic compartmental model, which explicitly features two categories of critically ill patients. This allows us not only to describe the past state of the epidemic but also to explore the effect of original and ea control strategies.

Unlike the majority of deterministic models used in epidemiology (Keeling and Rohani, 2008) , 95 ours runs in discrete time, with a time step equal to one day. This choice implies a formalism of numerical sequences that is closer to algorithmic syntax, less compact and less intuitive than that of ordinary differential equations. However, discrete-time models have the great advantage of implementing process memory. Indeed, the usual continuous-time models are said to be memory-less because the probability for an individual to leave a certain compartment, 100 for instance to recover, is independent of the time already spent in this compartment.

The absence of memory does not alter the qualitative properties of the model when an asymptotic behaviour is studied. It can however become extremely oversimplifying when studying phenomena over short time scales, that is in the transient regime of the system. Indeed, failing to account for the past may buffer sudden changes in the system. In the case of COVID-impossible to interpret the incidence data.

Although there are ways to correct ordinary differential equation models to deal with each of 110 the above-mentioned situations (e.g. delays, impulses, non-autonomous equations), the discretetime approach makes it possible to easily accumulate all of these options. It also approaches the very format of the data collected, monitored, and communicated during epidemics on a daily basis.

While accounting for non-Markovian transition times between clinical-epidemiological com-115 partments, our model remains deterministic. This means that contrarily to agent-based models, it can explore a variety of scenarios in a computationally efficient way. It thus combines the advantages from both stochastic modelling (namely simulating non-exponential-like waiting times) and law of large numbers (namely capturing the general trend of the system by averaging). Moreover, the deterministic nature of the model makes it quite parsimonious, allowing 120 unknown parameter values to be estimated, thereby preventing arbitrary choices. Likewise, its simple formulation opens ways to a variety of extensions and facilitates broad applications to health authorities.

By calibrating our model using French data, we show that epidemiological dynamics can be 125 accurately captured. More precisely, we estimate the R 0 of the French epidemics to be between 2.59 and 3.39, a number that decreased to between 21.3 and 27.1% of its value after the lockdown (95% likelihood interval).

By comparing our calibrated model to one without memory, we show the importance of allowing model parameters to depend on the age of the infection to accurately capture short-130 term dynamics. We also investigate the effect of earlier or later implementation of a national lock-down on epidemic peak intensity and timing.

Finally, we use our flexible and statistically robust framework to compare 3 main types of non-pharmaceutical interventions (NPI) that have been proposed in the literature: focusing the control on sub-populations who are more at risk, implementing periodic control, and deploying 2 Methods 140

Following classical epidemiological models (Kermack and McKendrick, 1927; Keeling and Rohani, 2008) , we group individuals with the same contribution to the dynamics into compartments whose densities are tracked over time. We consider age structured dynamics, by splitting each compartment into an arbitrary number of age groups, which are hereafter denoted by an 145 index i. This enhances the model with two key features. First, many COVID-19 clinical parameters are age-dependent, especially the infection fatality rate (Verity et al., 2020a) . With this age structure, we can adjust nationwide averages and capture demographic effects by matching demographic data to age-stratified medical data. Second, we can investigate agedifferentiated non-pharmaceutical interventions (NPI) . Note that this model can be extended 150 to formally take into account any kind of finer stratification, e.g. age, sex and comorbidities simultaneously. Furtermore, adding a discrete implicit spatial (also known as meta-population) structure is straightforward.

Since hospital admissions and critical cases are less frequent than non-severe infections, we assume that all transmission events occur in the community. A specific extension of the model 155 could be to focus on nosocomial transmission.

The structure of the system is shown in Figure 1 . Initially, all individuals in group i belong to the susceptible compartment, the density of which is denoted S i . These individuals can be infected with a probability Λ i called the force of infection (the i subscript indicates that not all ages are equally susceptible to infection). Each square represents a group of individuals who share the same clinical kinetics and who contribute equally to the epidemic dynamics. Contiguous squares form a compartment, in which each individual progresses day after day, therefore allowing to capture memory effects of the infection age. Pink boxes correspond to infected individuals in the community (depicted by the yellow area). Light blue boxes represent the hospitalized critical cases (the light blue area depicting the hospital). The purple-grey area corresponds to removed compartments that do not contribute to the epidemic. Arrows between boxes correspond to the daily flow of individuals from one compartment to the other. Dotted arrows depict transitions that occur with probability 1. The i subscript indicates the age group. For the sake of simplicity, only one group is depicted here and only one of the two probabilities is shown for each bifurcating transition (the other being its complementary to 1). Details about the compartments, flows and notations are provided in the main text. probability 1 − η k to move to the Y i,k+1 subgroup. No individual of a cohort of critical cases remains in the Y i,⋅ compartment after h days.

As detailed in Appendix S2.5, two groups of hospitalized critical patients are considered.

Those who have a substantial chance of recovering and who will benefit from a long stay (at least greater than one day) in an intensive care unit (ICU), are denoted by H i,⋅ . Those who will die, either after a short stay in ICU or in another ward are denoted by W i,⋅ . Upon hospitalization, a proportion of the incoming Y i,⋅ moves to H i,⋅ while the remaining part moves to W i,⋅ . Time

In the SIR-like modelling framework, the force of infection refers to the infection rate per 180 capita of susceptibles, often expressed as λ ∶= βI (Keeling and Rohani, 2008) . Equivalently, the instantaneous incidence is βIS = λS, which is the translation of the mass action law implied by the mean-field approximation made by such spatially unstructured models.

In our discrete time model, the force of infection Λ i is not a rate but a daily probability of infection (per capita of susceptibles from group i) that saturates with the prevalence. Indi-185 vidual contributions of infected individuals are not additive when prevalence is high because a susceptible host surrounded by infected individuals can be infected by several of them the same day. When prevalence is low, the probability of contamination by multiple infectors the same a day is low and the force of infection is well approximated by the sum of contributions of each infected individual. Λ i is therefore a monotonically increasing function of prevalence, 190 bounded by 1 and with a positive initial slope.

As we show in Appendix S2.2, deriving an expression for the force of infection is far from trivial for two reasons. First, we do not have a single class of infected individuals I as in most simple models. We therefore need to define and calculate an effective infectious density at time t (denoted I(t)), which can be seen as the number of individuals in the J and Y classes, 195 weighted by both the level of per-capita contact ratio c (t) and the generation time distribution, here approximated by the serial interval.

As detailed in Appendix S2.2, we find that the force of infection can be written as

where c i denotes the relative contact rate per capita of i-individuals (i.e. the current contact rate divided by the pre-epidemic baseline), I the effective infectious density, S 0 the initial

This model can be applied to any population where COVID-19 incidence data can be collected.

However, it also requires additional details, such as the time spent in ICU or the mortality, 205 which is why we apply it to France where at least part of the data is available.

All the model parameters are detailed in Appendix S1. We attempted to calibrate the models to the best of our ability based on the data available at the time of writing this report.

For instance, the calculation of µ i is based on reports from Santé Publique France. 

For several key processes in the model, the probability for an event to occur depends on the elapsed time. This is the case for the interval between contamination and hospitalization, through the probability distributions that underlie the η k sets of parameters. This is motivated 220 by the fact that, as early noticed by clinicians (Bouadma et al., 2020) , respiratory complications of COVID-19 arise in a quite narrow time window approximately a week after symptom onset.

Memory is also involved in the transmission term. capturing the empirical and effective non-homogeneous infectiosity period without recourse to 10 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020 . . https://doi.org/10.1101 /2020 additional compartments (such as a latent/exposed E, and convalescent densities) and their corresponding parameters.

We use two parameters to describe non-exponential distribution, which we assume to be Weibull distributions with a shape parameter greater than 1 (which captures the 'ageing' prop-235 erty) and a scale parameter. Other components of the model could also, in theory, accommodate memory processes, namely the time distribution from hospitalization to ICU discharge or death (ρ k and υ k ). However, preliminary fitting attempts reveal that an exponential (memory-less) distribution is more parsimonious. In Appendix S2.6, we explain in further details how the distributions for the various waiting times were defined and estimated from the data.

For comparison purposes, we also considered the continuous-time memory-less translation of this model (detailed in Appendix S2.9), which proved to be less adequate to capture the dynamics of the French epidemic in France, as shown by Fig.3. 

Parameter estimation is required prior to analyzing the outcome of the model when i) strong 245 simplifying assumptions about the phenomenon under study are made, and ii) some parameter values are not known with sufficient certainty.

First, the strongest assumption of the present model is the mean-field approximation: the population is supposed to be well-mixed, which is obviously not the case at the country scale.

Nonetheless, earlier works have shown that non spatially structured models produce conserva-250 tive estimates from a public health viewpoint, while their parsimony and tractability outweighs the greater precision provided by finer models (Keeling, 1999; Trapman et al., 2016) .

Second, the value of the basic reproduction number, of the day of onset of the epidemic wave, or of the lock-down effect have only been estimated in recent modelling works (ETE Modelling Team, 2020b; Danesh et al., 2020; Salje et al., 2020) . Those values might not directly apply 255 to our model and corrections might be needed. Likewise, known COVID-19 generation time distributions (Nishiura et al., 2020b) or IFR (Verity et al., 2020b) originate from field studies outside France. Besides, the estimated distribution of the ICU admission to death interval does not capture reporting delays (note that deaths and ICU admission are not coming up through the same channels). Parameter fitting is thus used to account for the uncertainty in initial parameter values, thereby improving predictions by re-calibrating their value.

Parameter inference was performed using nationwide daily ICU admissions, current ICU bed occupancy, as well as the cumulative number of deaths, all published by Santé Publique France and available on the French government data repository (Santé Publique France, 2020b).

We first located the region of highest likelihood using initial parameter values estimated 265 from data or compatible with the literature. Then, the maximum likelihood estimates (MLE) and associated 95%-intervals were calculated stepwise with respect to daily ICU admissions, ICU discharges, and finally daily mortality time series.

Confidence intervals for the simulation outputs are based on a collection of parameter sets assumed to be equally likely. These parameter sets originate from random draws according 270 to a multivariate Gaussian distribution (centered around the maximum likelihood parameter set and variances based on the confidence interval of each parameter) and only the resulting parameter sets whose likelihoods are not significantly different from that of the MLE are kept.

The 95% confidence intervals of the model's output are then calculated as the daily sample quantiles across all runs.

Further details about the procedure used for parameter estimation can be found in Appendix S2.7. The methods used to obtain predictions of the dynamics are shown in Appendix S2.8.

By analysing nationwide hospital data, publicly provided by Santé Publique France (Santé 280 Publique France, 2020b), we obtain maximum likelihood estimates for all parameters except those determining the generation time distribution, that was kept fixed following Nishiura et al. (2020a) .

The estimates and their likelihood intervals are summarised in Table S The fitted non-markovian discrete time model accurately captures the dynamics of both the daily hospital mortality and the daily number of ICU admissions since most of the data points fall into the 95% confidence intervals. As can be observed in Figure 2 (top), the model correctly approaches the number of daily admissions in the vicinity of the peak, which is crucial for hospital management at the local but also national level.

13 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . The blue and pink curves respectively represent the median daily ICU admissions and the median daily (hospital) mortality as generated by the fitted model.Turquoise triangles and red circles are the (rolling 7-day average) data counterparts. The black curve shows the median daily temporal reproduction number calculated from the simulated epidemic. The dotted horizontal line shows the reproduction number threshold value, i.e. 1. Bottom panel. The blue and pink curves respectively represent the median number of occupied beds in ICU nationwide and the median cumulative (hospital) mortality as generated by the fitted model. The turquoise triangles and red circles are the (rolling 7-day average) data counterparts. The purple dotted horizontal line shows the initial French ICU capacity, ca. 5,000 beds. The green curve shows the median proportion of the population that has recovered (and is assumed to be immune). The green dotted horizontal line corresponds to the median herd immunity threshold. The two vertical lines show respectively (from left to right) the beginning and the end of the French national lock-down. Shaded areas correspond to 95% confidence intervals.

14 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020 . . https://doi.org/10.1101 /2020 We also present the temporal reproduction number (R t ), which rapidly drops below unity following the onset of the national lock-down and was equal to 0.71 [0.69, 0.74] by May 11 according to the model. Notice that R t started to decrease before the lock-down onset, due to density dependent effects. In a model with strong host heterogeneity and so-called 'superspreaders', this effect would be even more pronounced. Other explanations include local satu-300 ration, staggered implementation of pre-containment measures, such as health communication campaigns, or improved patient management as diagnosis and therapy became more effective.

However, neither the structure of the model nor the level of detail in the available data makes it possible to identify the isolated impact of each measure on epidemiological dynamics. dynamics (though with a slight tendency to overestimate the declining ICU bed occupancy), which is essential in assessing the risk of a saturation of such hospital units, which would lead to an excess-mortality. The fitting of these data points could be improved with access to non-aggregated patient data or distributions of ICU residency time. The cumulative mortality curve is fitted with great accuracy, which allows us to use it as a comparison criteria between 310 control strategies in further analyses. Finally, the figure also shows that the level of population immunisation, the median of which we estimate at 2.37% ([2.27, 2 .48] % 95%-CI) by May 11, which is far below the classical group immunity threshold.

We also performed the parameter value inference by censoring the data to the right in order to assess the relevance of estimates obtained earlier in the epidemics. As shown in Appendix 315 S3.2, estimates with a censoring on Apr 15, that is one month before the final data point shown in Figure 2 , were already accurate. Estimates with an earlier censoring are qualitatively correct but the confidence intervals larger.

To illustrate the ability of our discrete model to capture COVID-19 short-term dynamics, we show in Figure 3 the best outcome of parameter inference for a Markovian (memory-less) 320 model. In spite of one additional degree of freedom compared to the focal model (see Appendix S2.9 for more details), the best fitted curves with the memory-less model (respectively in blue and pink) fail to capture the timing and amplitude of the peaks, while their decline is slower than the data. 15 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 

Using our estimated parameters, we then explored the effects of implementing the national lock-down a week earlier or a week later. As shown in Figure 4 , the peak was reached on Apr 8, with 7019 ICU beds occupied (the model estimates it to Apr 12 and 6920 beds). Enforcing the lock-down a week earlier (in green) would have led to an earlier and smaller epidemic peak with less than 1,500 ICU beds occupied on Apr 5. Conversely, another week of delay (in red) 330 would have led to a peak above 32,000 beds occupied in ICU on Apr 18, which is largely above the ICU capacity at the time (approximately 5,000 beds). Overall, in the range studied, each elapsed week multiplies ICU occupancy peak by more than 4.5, while delaying it by 3 weeks.

These differences also translate in terms of mortality. Implementing the lock-down a week earlier could have led, according to the model, to 13,300 [12,900-13,700] less deaths, while 16 All rights reserved. No reuse allowed without permission.

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

This model cannot estimate the effects of the end of the national lock-down in France on May 11, as it relies on data notifying events that occur on average two (ICU admission) to four (death) 340 weeks after contamination. However, assuming our estimated value hold after the lock-down is ended, we can explore the future dynamics as a function of the post-lock-down reproduction number. The aim is not to make predictions, since there is no data to test them, but rather to highlight the interplay between the post-lock-down reproduction number and NPI-enforcement timing with respect to a potential epidemic rebound. For illustration purposes, the chronicles 345 of four scenarios are provided in Appendix S3.4.

We also provide a graphical online interface (provided as a supplementary data) to maximise the reach of the model.

Following a more systematic approach of end of the national lock-down, we investigated 17 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . the effect of the NPI reinforcement delay on the ICU peak height and timing. The results are 350 shown in Figure 5 . As expected, a reproduction number less or equal to 1 (in green) does not require any control reinforcement. Conversely, higher values of R t trigger an epidemic rebound, that can saturate France's ICU capacity (ca. 5,000 regular beds) with values as low as 1.2 [1.1, 1.3] if an appropriate control response is not implemented before mid-July. For R t ≥ 1.3, we find that a reinforcement as early as mid-June is necessary to preserve the national health 355 system ( Figure 5 , top panel). Additionally, Figure 5 (bottom panel) shows that if a massive peak is to occur (R t ≥ 1.7), it will likely be in the second part of July, even if an appropriate response in mounted in mid-June or later. For lower peaks (R t < 1.5), the height of which can be substantially reduced by timely reinforcement of control, occur within one month (early responses) to two weeks (late responses) after NPI implementation.

18 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. .

Colors indicate post-lock-down reproduction number, which ranges from 1 (green) to 2.1 (brown) (see the legend with the 95%-CI). The abscissa indicate the date of implementation of renewed NPI that bring the reproduction number to R t = 0.91 [0.85, 0.98]. Each dot represents the median highest number of ICU occupied beds (top panel) and the median date at which the peak is reached (bottom) (the bars indicate the 95% CI). Dots at the bottom of each panel correspond to an absence of epidemic rebound (therefore the peak is artificially considered to be on May 12). The vicinity with the threshold value 1 for R t explains the large CI of the peak dates for the green and turquoise scenarios.

In this section we use the model previously fitted on the data available as of May 12 to explore control strategies that belong to three main classes of non-pharmaceutical interventions (NPI):

adaptive (indicator-triggered) lock-down, periodic alternation of lock-down and release phases (not indicator-triggered), and age-differential control.

19 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . https://doi.org/10. 1101 /2020 We compare scenarios using the predicted cumulative hospital mortality on Dec 31 2020.

The choice of this criterion is motivated by three reasons. First, the mortality time series is the one captured with the greatest accuracy by the model. Second, in the absence of pharmaceutical solutions, deaths are a more stable indicator of the state of the epidemic (e.g. the timing of hospitalizations may change as knowledge of the disease accumulates). Third, peak ICU activity 370 and concern is more likely to vary among countries, which makes it less generic.

Importantly, the following simulations are in no way intended to be statistical forecasts of the COVID-19-related death toll at the end of 2020 in France or elsewhere, but rather a numerical illustration of the non-trivial interplay between the degrees of freedom in each of the NPI strategies considered. For the sake of parsimony, we set the implementation of all 375 further considered strategies by May 12, thus avoiding to consider multiple scenarios between lock-down lifting and new NPI reinforcement (this aspect having been addressed above).

To facilitate the comparison in terms of control burden of the following strategy, we hereafter reason in terms of per capita contact ratio (PCCR), denoted by c. This dimensionless number aims to quantify the average potentially infectious contacts an individual has per unit of time, 380 relative to the pre-epidemic baseline (see Appendixes S2.2 and S2.3 for more details about how c is formally introduced in the model). This definition has the following consequences:

• the average PCCR before lock-down is c pre-lock = 1,

• under the proportionate mixing assumption, the encounter rate at date t of two individuals belonging to age groups i and j respectively is proportional to

• the threshold PCCR value corresponding to a reduction of the reproduction number to

By standardizing the control parameters of the model on a relative scale, the PCCR allows us 390 to perform comparisons of control strategies that are independent from the estimated (past) or arbitrary (scenarios) values of both basic and temporal reproduction numbers. In addition, it provides an easy way to picture the control burden of each strategy. Indeed, the average PCCR c (t) ranges from c lock (the strongest control implemented so far) to 1 (the fully relaxed, 20 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . 'pre-lock-down', situation). Between these two extremes, the remarkable value c t demarcates 395 the threshold above which an epidemic rebound is made possible.

Notice that this formalism caputres both reduction in contact rate (e.g. working from home) and probability of transmission per contact (e.g. wearing a mask).

The adaptive lock-down popularised by Ferguson et al. (2020) consists in triggering a lock-400 down (or possibly other broad restrictive measures) if an epidemic rebound is detected by a relevant indicator. This requires to carefully identify the threshold above which the lock-down is triggered. Following Ferguson et al. (2020), we use ICU admissions to quantify this threshold because they achieve a good balance between reflecting the epidemiological dynamics in the general population (the sampling of newly infected individuals is more homogeneous than with 405 testing) and limiting the delay between the data and the state of the epidemics (we estimate this delay to be 2 weeks, while we estimate 4 weeks for mortality data). Figure 6 shows that cumulative mortality is exponentially correlated with the daily ICU admission threshold, over the investigated range. As a consequence, there is no remarkable inflexion point that would justify a particular value, which would leave decision makers to 410 coldly balance socio-economical costs (implied by low thresholds) with potentially saved lives (endangered by higher thresholds). Such a question is once again out of the scope of the present work.

21 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 22 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. For the sake of clarity, only three parameter sets, chosen at random among the pool, illustrate here the dynamics generated by an adaptive lock-down strategy. In this scenario, the on and off thresholds are set to 15 daily ICU admissions nationwide. Symbols as in Fig.2 .

By definition, an epidemic remains under control if R t is kept below 1. If control measures alternate according to some periodic pattern, the instantaneous value of R t is of little interest and we consider instead its average value over one cycle of interventions, denoted by R t . For the sake of simplicity, we consider a NPI that alternates between a hard phase (typically, a lock-down) and a relaxed (or release) phase, independently from any field indicator. Let us 425 denote by c r and c h the average per capita contact ratio during the relaxed and hard phase respectively and p r the time proportion of each NPI cycle allowed to the relaxed phase.

The epidemic is under control if

The NPI cycle can contain more than two phases, but one can always pool them into those favorable to social and economic life, represented by the couple (c r , p r ), and those concentrating the necessary public health restrictions, namely the couple (c h , 1 − p r ).

23 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. .

Assuming that hard phases have the same effect as the first lock-down, we can set c 2 h = 1 − κ (see Appendix S2.3). Then, the maximum proportion of time allowed to the relaxed phase is given by

which is a decreasing function of c r . This expression captures the intuitive trade-off, shown in Figure S -8, which is that the lower the control in the relaxed phases, the shorter they can be.

If there is no control in the relaxed phase (c r = 1), then, based on our estimates for R 0 and κ (Table S -5), controlling the epidemics requires that hard control be implemented 88% of the time.

Here, we maximize the product between intensity and duration (p r,max c r ) because it partly captures the absolute degree of freedom for a fixed period of time (i.e. try to be the less constrained as possible for the longest time period as possible). Noticing that d dcr (p r,max c r ) < 0 leads to the conclusion that the optimum is to seek for the lowest per capita contact ratio, which is obtained by setting equation 3 equal to 1. This yields c ⋆ r = c t = R Importantly, the simple product p r,max c r does not need to be the objective quantity to maximize and parametrizations of p α r,max c β r , with α, β > 0, based on socio-economical arguments, should be investigated. However, this question is out of the scope of the present work.

24 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. We show the cumulative mortality by the end of the year 2020 if a periodic control is implemented from May 12. For each of the four NPI cycle durations, the model is run for four values for c r shown in different colors, with a proportion of time satisfying eq. 3 (truncated to an integer number of days). For example, in the red scenario, p r,max (0.59) = 88%, so for the weekly cycle the relaxed phase was set to floor (0.88 ⋅ 7) = 6 days.

Up to now we only investigated the proportion of the time spent in a relaxed or hard control 445 phase. However, a 0.5 proportion can correspond to a 1 week / 2 week periodicity or to a 1 month / 2 month periodicity. In fact, the trade-off relation shown in Figure S -8 theoretically works for any time split where a proportion p r,max is spent in a relaxed phase with per capita contact ratio c r and a proportion 1 − p r,max ) corresponds to a strong control phase where per capita contact ratio is restricted to κ. To be implemented, such periodic control should consider 450 phases lasting at least several days. Furthermore, a periodicity greater than a month exposes populations to a greater risk during the relaxed phase. Figure 8 suggests that the duration of the overall cycle has an effect up to 5,000 deaths, over the explored range. A weekly cycle yields the lowest death tolls, both in terms of median and CI. Furthermore, if the period lasts a month or less, a lower control of the epidemic in 455 the relaxed phase yields the lowest cumulative mortality. Morecisely, these results shows that allowing 1 moderately relaxed (the average PCCR being equal to 80%) day per week within a prolonged lock-down has a better mortality outcome than the adaptive lock-down set with 15 25 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . daily ICU admission as a threshold (approx. 8,000 deaths less by the end of the year). 460 Age-group targeted control measures are motivated by the age dependent risk profile (Verity et al., 2020a; Collaborative et al., 2020) . However, this model could easily include other risk factors such as heart disease, diabetes, or obesity. One of the simplest strategies is to incite people at risk to stay at home. On the other side of the age pyramid, a finer age-specific strategy could be to close schools and universities. More generally, political and economical 465 incentives to telework could also be considered as age-specific control measures.

We divided the population into three groups according to the cut-off ages of 25 and 65 years old. The first group is motivated by the control leverage represented by school and university closure, while the third is motivated by preventing viral circulation within the age group having the highest IFR. Each age group was assigned a fixed PCCR until the end of the year: 50, 58, 470 62, 65 and 75%. This led to 125 scenarios that are shown in Figure 9 . The ratio of daily interaction between two individuals relative to the pre-epidemic baseline is the product of their individual contact ratios.

The results of the numerical investigation of 3-age group-differential controls are shown in Figure 9 . This scatter plot shows that cumulative mortality is well explained by demographic-475 averaged PCCR. Therefore, the variance in contact ratio among age groups (not shown on the plot) has little impact on the long-term epidemic dynamics. This result suggests that there is no easy way of taking advantage of the strong correlation between COVID-19 complication risk and age, unless implementing an unrealistically restrictive quarantine for at-risk groups.

26 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Figure 9 : Effect of age-differential control interventions on cumulative mortality by the end of 2020 Each dot represents a simulation of the model with age-differential restrictions implemented after May 11. The abscissa corresponds to the demographic-weighted average in PCCR and the ordinates represent the median final death toll by the end of the year if no other solution is implemented (bars indicante 95%-CI). The color of the dots shows the temporal reproduction number by May 18: in green the epidemic is under control and the 97.5% quantile is less that 1, in purple the median is below 1 but the 97.5% quantile is above 1, in pink the median is above 1 but the 2.5% quantile is below 1, and in grey the epidemic is out of control, the 2.5% quantile is greater than 1. Vertical bars show the 95%-likelihood interval of the expected per capita contact ratio threshold above which the reproduction number is above 1, in a homogeneously-controlled population.

27 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. However, given the magnitude of the risks at stake, the diversity of fields involved, and the number of unknowns (such as climate effects, immunity duration), mathematical modelling 490 predictions should be handled with caution in the decision process.

Informing decision making represents a challenge because most mathematical modelling in epidemiology relies on continuous time deterministic models. These offer a wide palette of analytical tools, but they also become less accurate on short time scales. Conversely, stochastic models, whether agent based or not, offer a much more precise picture of the early stages of the 495 epidemics. However, they become less necessary once the outbreak threshold has been crossed and epidemic dynamics are essentially deterministic.

We therefore developed an original framework at the crossroads of individual-centered and compartmental modelling approaches, tailored to the COVID-19 natural history, which has two great advantages. First, its discrete time structure, shared with that of common epidemiological 500 data, allows us to assume any distribution for model processes, therefore introducing what is known in the literature as memory effects (or 'non-Markovian' processes), related to the age of the infection. As a result, we obtain a much better fit than classical memory-less (or 'Markovian') models on intermediate timescales (weeks or months). Second, the deterministic nature of the model allows us to perform extremely fast simulations, especially compared to 505 agent-based modelling that requires the drawing of millions of random numbers for a single simulation. The computational performance combined to the great parsimony, that still allows it to accurately fit the observed epidemiological dynamics, makes this model a relevant tool that can be easily transposed and deployed to other settings, countries or scales -even with data limited to those publicly available, as it is the case in the present work.

First, we used our model to infer three key parameters of the epidemics. Our R 0 is comparable to that already computed in France based on their 95% confidence interval (Di Domenico et al., 2020; Roques et al., 2020; Salje et al., 2020; Hoertel et al., 2020) . We also estimate the temporal reproduction number R t after the lock-down to 0.71, which is higher than that estimated by earlier studies (Salje et al., 2020; Hoertel et al., 2020) , and the number of lives 515 potentially saved (more than 50,000) compared to a lock-down implemented a week later.

Once the parameters estimated, we investigated potential scenarios for control, based on the assumption that the major characteristics of the epidemics will not change, meaning for instance that a seasonal effect would have a limited impact on transmission dynamics. Our results suggest that, in the case of a rise of the reproduction number above the unit threshold 520 value, a reinforcement needs to be implemented within 5 weeks after lock-down lifting to prevent a massive second epidemic peak.

Then, we explore two types of cycling strategies. The first strategy, known as adaptive lockdown and popularised as 'stop-and-go', consists in alternating strict lock-down and absence of control (Ferguson et al., 2020) , the shift between the two being based on the incidence of daily 525 patient admission in ICU. We show that this requires a low threshold without a particular value standing out. The second strategy explored is the fixed-period alternation between lockdown and relaxed phases. Our work highlights that, for the adaptive lock-down to perform better than the latter, both in terms of casualties and time spent under relaxed control, part of the restrictions have to be maintained between lock-downs. Besides, we investigated age-530 differential control strategies and show that there is little public health benefit to be gained from the variance of control restrictions between age groups, unless implementing measures even more restrictive than the past lock-down. However, synergy between age differential and periodic control strategies is left to be explored.

The model makes several strong assumptions. First, there is no spatial structure. This (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 24, 2020. . https://doi.org/10.1101/2020.05.22.20110593 doi: medRxiv preprint contact rate. We also neglected fomite transmission (see Ferretti et al. (2020) for an example) 540 and assumed perfect and lifelong immunity against reinfection due to currently insufficient data on immunity. One simplifying assumption we made is that mortality probabilities do not vary over time, whereas in practice hospital saturation could affect mortality, whether related to COVID-19 or not.

In terms of outlook, this work lays the foundation for an online application soon to be 545 released, as an update of the earlier version COVIDSIM-FR (ETE Modelling Team, 2020a).

Next challenges include taking to account for possible changes in parameter values with time, mainly detecting and estimating seasonal effects. If these are of significant impact, model fitting could be adapted, either by estimating parameters within defined time windows or by left-censoring the data as time goes by. Lastly, the NPI analysis exposed in this work sets 550 the ground for the exploration of finely tuned public health measures accounting for spatial heterogeneity and combining advantages of adaptive, periodic and group-differential modalities, with the aim of avoiding an epidemic rebound while minimizing its population burden.

30 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. We acknowledge support from the University of Montpellier, the CNRS and the IRD.

We have no competing interests. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

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

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

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

The copyright holder for this preprint this version posted May 24, 2020. 37 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020 . . https://doi.org/10.1101 /2020 (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 24, 2020 . . https://doi.org/10.1101 /2020 

Summary table with all the notations used in the study. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted May 24, 2020. 3 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 4 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. Instead of classical ordinary differential equations (ODE), the dynamics of the model satisfy a system of recurrence relations (one could as well write as finite difference equations (FDE)), which is detailed below. For the sake of simplicity, we omit time dependence in the notations.

Instead, X i denotes a density at a given time t and X ′ i the density at time t + 1. The equations are formally identical for all age groups i. Between-group dynamics are coupled through the 35 forces of infection Λ ⋅ defined in the next subsection.

(S-1)

As indicated in the main text, our goal is to have the force of infection function follow the well-known Michaelis-Menten (or Holling type II physiological response) function of the form

x ↦ x a+x . However, we cannot use the prevalence as the argument of Λ i here because all 40 infected individuals, whether they are critically ill (Y ) or not (J), do not contribute equally to transmission events. This heterogeneity in contagiosity originates from differences in infection ages (individuals contaminated 6 days earlier are more contagious than those 10 days earlier) and in contact rates.

To address this issue, we introduce the effective infectious density I (t),

All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . https://doi.org/10. 1101 /2020 which is the sum over all infected community compartments weighted by both the generation time ζ k , which is the time between the infection of an 'infector' and the infection of his or her 'infectee', and the per-capita contact ratio c i (t). The latter is defined as the current contact rate per-capita of individuals of age group i (k i (t)) relative to their pre-epidemic baseline contact rate (k i,0 ),

With I kept constant, Λ i is expected to display a Michaelis-Menten behavior (i.e. positive initial slope, increasing and upper-bounded) with respect to the current per-capita rate as well.

Consequently, the force of infection should satisfy

, (S-4) as shown in Fig.S-1 .

We then derive the expression of a using known parameters. To do so, we consider the 55 probability for a given susceptible individual to be part of the first generation of cases, that is, to have been infected by the index case. Let us denote Λ i,index (t) the probability of being infected by the index case on day t. By definition of the basic reproduction number, the index case infects R 0 secondary cases, on average, by the end of its contagious period. Under the mean-field approximation, the probability of being part of these secondary cases is simply R 0 S 0 60 (which is an extremely rare event). Summing over all possible days, we therefore have

For simplicity, we make the following three assumptions:

• the index case is not critically ill (less than 5% of cases are),

• the index case belongs to the same focal age group i (More generally, one should average the per-capita contact rate of the index case weighted by the (unknown) probability

• the index case has infected all his or her secondary cases before public health measures are implemented.

It follows from these assumptions that if we set the contamination day of the index case to t = 0, the effective infectious density I restricted to the index case can simply be written as

(note that k i (t) = k i,0 over the considered period of time, hence c i (t) = 1 and all densities are equal to 0 except J i,t = 1).

Applying these results to equation S-4, the daily probability of infection by the index case 75 therefore becomes

Now, from the magnitude comparison (S-5), we have a ≫ k i (t) ζ t , and hence

In the limit of low prevalence (or low contact rates), the mass action law is recovered. Using equation (S-5) and the fact that the sequence (ζ t ) t≥1 being a mass function, it sums up to 1, we get

Combining equations (S-4), (S-3) and (S-6) finally leads to the generic expression of the force of infection:

.

This is the equation indicated in the main text.

In the special case where strong public health control measures such as lock-down are being implemented, all individuals may exhibit similar per capita contact rates, k i (t) = k lock . From 7 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . Figure S -1: Daily force of infection as a function of the product of daily contact rate and effective infectious density Λ i is the probability for one susceptible individual from group i to be infected a given day. The initial non-zero slope comes from the law of mass action implied by the mean-field approximation. The derivation of Λ i is based on the remarkable coordinates shown near the origin of the graph (not to scale). When the effective infectious density is equal to 1 (i.e. as if the generation time distribution were concentrated in a single day) and the contact rate is that in absence of any health measure (denoted by k 0,i ), then, by definition of the basic reproduction number, Λ i equals R 0 S 0 .

(S-4), it is then possible to express the force of infection in the following way:

where I + is the effective infectious density as if it were calculated in absence of health measures 85 (i.e. with c i = 1). In equation (S-8), one can interpret the quantity c 2 i,lock ∶= k lock k i,0 2 < 1 as a factor lowering the basic reproduction number R 0 . The lock-down effect, defined as the reduction of R 0 due to this measure, can be calculated by averaging c i,lock over age groups (according to demography). Hence,

The largest and most reliable nationwide data for the COVID-19 epidemic in France is that of daily COVID-related death toll in hospitals, communicated daily since Feb 16 2020 by the national public health agency (Santé Publique France). This is why our model neglects COVID-related deaths occurring outside hospitals. In particular, we removed nursing homes (or EHPAD) from calculations.

Starting from Mar 18 2020, two additional time series are communicated: daily admissions in intensive care units (ICU) and current ICU occupied beds. While the former capture the dynamics from contamination to ICU admission, the latter captures moreover the kinetics of ICU stay.

These time series are altered by week-ends and bank days: e.g. death tolls are notably lower 100 on Sundays than previous days, while it increases the next Mondays. It has even been suggested that reporting delays propagate also to Tuesdays. In order to smooth these artifactual weekly oscillations, a right-shifted 7-days moving average was performed over all time series prior to analysis. We will refer to these smoothed datasets asM for daily hospital mortality,Ã for daily ICU admissions andH for current ICU occupied beds. SinceH contains part of the cumulative 105 information ofÃ, we also consideredD, the cumulative counterpart ofM to equilibrate the first step of the fitting procedure (see below).

Symptom severity of COVID-19 is increasingly classified into mild, moderate, severe and critical. Because we rely on hospital mortality and ICU flow data, we focus on critical cases, i.e. the 110 ones concerned by intensive care and COVID fatality. In the absence of detailed large-scale hospitalization data, we made the following assumptions:

• non-critical cases are not admitted into ICU (even though some do need hospitalization),

• all critical cases need hospitalization, and only survive if they go through intensive care.

For medical reasons not addressed here, in France not all critical cases are admitted soon 115 enough into ICU. Part of them die in non-intensive care wards, while others die shortly after entering the ICU, therefore not contributing to ICU bed occupancy. We therefore need to 9 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . https://doi.org/10.1101/2020.05.22.20110593 doi: medRxiv preprint estimate two key criticality-related probabilities, namely θ i , the proportion of critical cases within age group i, and ψ i , the proportion of critical cases in age group i that contribute to ICU bed occupancy (i.e. their stay in the ward exceeds one day). Because these probabilities 120 differ among age classes and improper averaging could lead to substantial bias (see below), we first need to make calculations focused on the smallest age stratification unit (usually a decade). In the following, θ (a) and ψ (a) are the age-specific critical case frequency and longstay admission given critical illness respectively. In addition, P a [X 1 X 2 ] reads as the probability of event X 1 given X 2 has occurred for an individual of age a. 125 Let us consider the four events needed to derive θ (a) and ψ (a) from available data:

• I, being infected by SARS-CoV-2,

• U, being hospitalized,

• B, occupying an ICU bed for more than a day,

• D, dying at the hospital from The interplay between these events is formally depicted by the tree diagram in Fig We need to find two independent equations involving probabilities related to these events for which age-stratified data is known in order to solve θ (a) and ψ (a). First, we need to identify such data:

• P a [D I] is the proportion of deaths among COVID-19 infected individuals, better known as the Infection Fatality Rate (IFR), which has been calculated using the Diamond Princess data by Verity et al. and corrected for non-uniform attack rate by Ferguson et al. and will be denoted by IFR (a) hereafter,

• P a [D B ∩ U ∩ I] is the proportion of deaths among COVID+ ICU hospitalized patients; this age stratified data has been communicated by Santé Publique France as a weekly 140 epidemic report on May 7 2020 (Santé Publique France, 2020a), and will be denoted by (Tables S1 and S2) .

A first equation comes by noticing that the probability for an infected individual to die from COVID-19 is the probability of developing critical illness if infected (namely θ (a)) times the proportion of deaths among critical cases. The latter is the sum of critical cases that die in ICU after a stay longer than one day (µ (a) ψ (a)) and the critical cases that are not lengthily 

All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020 . . https://doi.org/10.1101 /2020 The unknown probability can be calculated as follows

After elementary calculations, we finally get

, hence both θ (a) and ψ (a) can be calculated from age-stratified available data.

The mean proportion of critical cases among a specific age group, θ i , is simply the demographicweighted average of θ (a) over the considered ages, as infection samples uniformly the susceptible compartment (NB: the IFR stratified data used here from (Ferguson et al., 2020) is already 160 corrected for non-uniform attack rate). However, as the probability of the next events related to critical illness are not homogeneous with respect to age, the age group averages ψ i and µ i cannot be weighted directly with relative demographic age frequencies. Instead, they must be calculated by taking into account that ages with higher ψ (a) and then µ (a) will be over-sampled in Y i → H i and H i → D i transitions respectively.

To account for this bias and as well to allow adjusting the parameters to both the model (for ψ (a) and µ (a)) and the French epidemic (for the IFR), we introduce a series of corrections detailed hereafter in the calculation of the parameters used to run the model.

First, we account for the fact that the ICU fatality rate might mix both short and longstay patients, while our model splits these two flows. The corrected age-specific long-stay ICU 170 fatality rate will be denoted byμ (a) and calculated as the product of the corresponding data and a correcting factor denoted by C M , i.e.μ (a) ∶= C M µ (a). Likewise, the age-specific IFR (which was not estimated from French data) will be corrected asÎFR (a) ∶= C F IFR (a). Now, equations S2.5 and S2.5 rewrite aŝ

1 − (1 −μ (a))ψ (a) .

12 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. .

Note that there is no need for a fourth correction factor as C M , C F and C Ψ already capture possible corrections forθ. These corrections are yet of unknown values but will be estimated by the fitting procedure detailed below. As the corrected values are still probabilities, the correction factors are necessarily upper bounded, respectively by (max a µ (a)) −1 , (max a IFR (a)) −1 , and (max a ψ (a)) −1 .

The last step consists in age-group averaging as mentioned above. Let us denote by f (a) the 180 frequency of individuals of age a in the French metropolitan population (after having removed the ca 730,000 individuals living in nursing homes). We call the i-group relative frequency of age a as the standardized age frequency

where A i is the set of ages belonging to age group i. As previously implied, the frequency of critical cases in age group i is the straightforward demographic weighted average

The frequency of long-stay ICU patients among hospitalized critical cases in age group i is then weighted by both the relative age frequencies and the ratio of critical illness probability to the group average, i.e.

Finally, as the last event to occur, the average fatality ratio for long-stay ICU patients belonging to group i must be corrected by the ratio of the product of the two previous frequencies 190 relative to the group average φ i , i.e.

Four time distributions underlie the dynamics: the generation time (or index-contaminationto-secondary-contamination interval), the contamination-to-hospitalization interval of critical 195 cases, the ICU length of stay and the hospitalization-to-outside-ICU death interval of critical 13 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020 . . https://doi.org/10.1101 /2020 cases. Each of these events can be seen as a random waiting time variable, denoted by Z ○ (capital zeta), H ○ (capital eta), P ○ (capital rho) and Υ ○ (capital upsilon) respectively.

Initially, all four random variables are assumed to follow Weibull distributions, with shape parameters greater than one. Such distributions are widely used in the biomedical literature, 200 along with Gamma and Lognormal distributions, for fitting the probability density function (PDF) of ageing processes, for which the probability for the focal event to occur increases with lapsed time (Bolker, 2008) . Preliminary exploratory fittings of the model to daily mortality data indicated that maximum likelihood estimates of the shape parameter of P ○ and Υ ○ were close to unity. Because the computational procedures used for maximum likelihood estimation 205 misbehave in the vicinity of parameter range boundaries, the shape parameter of these two distributions was set to 1, turning them into exponential distributions by fitting (the exponential distribution being a special case of Weibull distributions), though not by assumption.

Weibull distributions have a right-unbounded support [0, ∞), which means that true distributions require truncation for obvious computational reasons. Let us introduce the generic 210 notations Ξ ≡ Z, H, P, Υ and x ≡ g, h, r, u. We construct the right-truncated analogous distributions by setting the finite upper boundary of their support x ∶= min {n ∈ N ∶ F Ξ ○ (n) ≥ 0.99}, i.e. the upper-integer-rounded 99%-quantile of the original distribution Ξ ○ , where F ⋅ denotes the cumulative distribution function (CDF). The truncated distributions Ξ are therefore such that their CDF satisfy F Ξ = F Ξ ○ F Ξ ○ (x), and having as their supports F Ξ (Ω) = [0; x]. 215 Dynamics unfold in discrete time in our model, which means these continuous distributions need to be discretized into sequences to be implemetend into the framework. The generation time sequence is straigtforwardly defined as ζ k ∶= F Z (k) − F Z (k − 1) for 1 ≤ k ≤ g. Indeed, transmission events do not affect the progression of the infector within its compartment. However, the three other sequences of parameters, ξ k ≡ η k , ρ k , υ k , represent the proportion of individuals 220 that leave the compartment k days after having entered it. These parameters need to capture the probability that the corresponding event occurs on day k but they also need to be standardized by the probability of not having left the compartment by that day. They are therefore calculated as

The fitting procedure was performed using the mle2 routine from the bbmle package (Bolker et al., 2020) implemented in R (R Core Team, 2020). Starting from the initial parameter values v 0 , an ordinary least square optimum v 1 was found by minimizing the euclidean distance tõ we calculated their likelihood with respect toM and kept only those whose likelihood was not significantly different from that of v 3 , according to Wilk's theorem (Wilks, 1938) . Sampling stopped once 10 3 draws have satisfied the condition. Importantly, we considered all the 15 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . retained parameters equivalent from the likelihood point of view, i.e. they are assumed to rep-255 resent equally likely versions of adjusted parameter sets. Their diversity thus allows to account for uncertainty in the real parameter values.

For any further analysis of the model, system S-1 was run independently with each of the 10 3 parameter sets. The confidence intervals of simulated tracked densities (S (t), J (t)...) as well as any derived quantity (R (t), ι (t)...) were then simply calculated for each time point as 260 the unweighted 2.5% and 97.5% sample quantiles of the 10 3 outputs at the given time point.

The central estimations correspond to the median value of these distributions.

Tracked densities are the number of individuals in each clinical-epidemiological compartment X i,k , the dynamics of which satisfy S-1 and are thus directly provided by numerical iteration of 265 the recurrence relations. However, several quantities of interest require additional calculations, hereafter exposed.

Let us first introduce two notations for the sake of concision. Tracked densities without indices will refer to sum over all groups, and, for multiple-days compartments, over all possible days of progression as well: X ≡ ∑ i ∑ k X i,k (e.g. J (t) represents all individuals belonging to J i,k 270 for all groups i and all ages of infection k). The daily difference ∆X (t) ≡ X (t) − X (t − 1) is straightforwardly used to extract the instantaneous dynamics from a cumulative time series.

The following three time series are crucial as they are used for likelihood calculations with respect to their data counterpart:

The next times series are intermediate calculations required for further key quantities.

• I (t) ∶= J + Y is the community infectious density and represents all not hospitalized infected individuals, which can be used to estimate the expected proportion of COVID 280 PCR+ in the general population, 16 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . https://doi.org/10. 1101 /2020 • C (t) ∶= ∑ τ ≤t ∑ i (J i,1 (τ ) + Y i,1 (τ )) = S 0 − S (t) and ∆C (t) are the cumulative and instantaneous incidence respectively. The latter is used for the calculation of the most scrutinized indicator in epidemic monitoring, namely the temporal (or effective) reproduction number, R (t), which we calculate here 285 following Wallinga & Lipsitch, at the time of the infectees' contamination:

In this work, (population) immunization ι refers to the proportion of individuals that have been infected by SARS-CoV-2. We assume waning immunity is negligible at this timescale. Its calculation is simply

In absence of waning immunity from the host and antigenic drift from the virus, and assum-290 ing public health measures are fully relaxed, further epidemic can only be passively prevented if the herd immunity threshold is reached, i.e. ι ≥ ι h where

A classical result from Kermack & McKendrick is that if the epidemic has started spreading, the final cumulative relative incidence will not stop at the herd immunity threshold, but continue to a greater value, known as the final size proportion, that has no close form solution, but can 295 implicitly be defined as

Finally, current prevalence in the community is simply given by π (t) = I S + I + R .

17 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . https://doi.org/10.1101/2020.05.22.20110593 doi: medRxiv preprint

The Markovian continuous-time model analogous to S-1 corresponds to the following set of ordinary differential equations

where the force of infection is given by

Note that a latent compartment E has been added to account for a delay between contamination time and the beginning of the infectious period. The average latency and infectious period are equal to ω −1 and γ −1 . By construction, all transition times are exponentially distributed. 18 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 19 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 20 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 21 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 22 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. 23 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. Here are shown four runs of the model with distinct reproduction number after lock-down lifting and control response timings. Importantly, these plots are not statistical predictions of the future but only illustration-purpose qualitative explorations. The blue and pink curves respectively represent the median number of occupied beds in ICU nationwide and the median cumulative (hospital) mortality as generated by the fitted model. The turquoise triangles and red circles are the (rolling 7-day average) data counterparts. The black curve shows the median daily temporal reproduction number calculated from the simulated epidemic. The dotted horizontal line shows the reproduction number threshold value, i.e. 1. The purple dotted horizontal line shows the initial French ICU capacity, ca. 5,000 beds. The vertical corresponds to the end of the French national lock-down (May 11). Shaded areas correspond to 95% confidence intervals Note: the abscissa scales differ across panels. Top left panel. After lockdown lifting, the reproduction number R t = 0.91 ([0.85, 0.98] ), the epidemic is under control and vanishes spontaneously. Top right panel. After lock-down lifting, the reproduction number R t = 0.91 ([0.85, 0.98]), a second wave arises in absence of any control response. Bottom left panel. As previously, but a reinforced set of NPIs is implemented by Jul 15, reducing the reproduction number below 1. Bottom right panel. As previously, but a reinforced set of NPIs is implemented earlier, by Jun 15.

24 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. The maximum proportion of time possible to be spent in the relaxed phase of a NPI cycle without triggering a new epidemic, denoted by p r,max , is governed by three quantities, following eq.3: the basic reproduction number, R 0 , the per capita contact ratio during the hard phase of the NPI cycle (here set to the same value as during lock-down, κ), and the per capita contact ratio during the relaxed phase, c r . This figure illustrates the trade-off relation between p r,max and c r : the central curve represents the median relation while the surrounding shaded areas correspond to the 95% confidence interval. Vertical dashed bars indicate the 95% CI of the average per capita contact ratio threshold associated to R t = 1 (therefore applicable the whole time, without need of any harder phase), this interval is [56, 60] %. Horizontal dashed bars indicate the 95% of the proportion of the relaxed phase in each NPI cycle that guarantees associated to R t = 1, this interval is [6.2, 12] %. 25 All rights reserved. No reuse allowed without permission.

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

The copyright holder for this preprint this version posted May 24, 2020. . https://doi.org/10. 1101 /2020 

Incubation period of 2019 novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China

Ecological models and data in R

bbmle: Tools for general maximum likelihood estimation

Severe Salje, H

Estimating the burden of SARS-CoV-2 in France

COVID-19 : Point épidémiologique hebdomadaire du 7 mai 2020

This work is based on a report and online shiny application published online on Apr 6, 2020