Impact of Lockdown on the Epidemic Dynamics of COVID-19 in France.

The COVID-19 epidemic was reported in the Hubei province in China in December 2019 and then spread around the world reaching the pandemic stage at the beginning of March 2020. Since then, several countries went into lockdown. Using a mechanistic-statistical formalism, we estimate the effect of the lockdown in France on the contact rate and the effective reproduction number R e of the COVID-19. We obtain a reduction by a factor 7 (R e = 0.47, 95%-CI: 0.45-0.50), compared to the estimates carried out in France at the early stage of the epidemic. We also estimate the fraction of the population that would be infected by the beginning of May, at the official date at which the lockdown should be relaxed. We find a fraction of 3.7% (95%-CI: 3.0-4.8%) of the total French population, without taking into account the number of recovered individuals before April 1st, which is not known. This proportion is seemingly too low to reach herd immunity. Thus, even if the lockdown strongly mitigated the first epidemic wave, keeping a low value of R e is crucial to avoid an uncontrolled second wave (initiated with much more infectious cases than the first wave) and to hence avoid the saturation of hospital facilities.

1. Introduction

COVID-19 epidemic was reported in the Hubei province in China in December 2019 and then spread around the world reaching the pandemic stage at the beginning of March 2020 (1). To slow down the epidemic, several countries went into lockdown with different levels of restrictions. In the Hubei province, where the lockdown has been set long before the other countries (on January 23), the epidemic has reached a plateau, with only sporadic new cases by April 15 [from the data of Johns Hopkins University Center for Systems Science and Engineering (2)]. In France, the first cases of COVID-19 were detected on January 24, and the lockdown has been set on March 17. This national lockdown means important restrictions on movement, with a mandatory home confinement except for essential journeys including food shopping, care, 1 h individual sporting activity and work when teleworking is not possible, and closing of the borders of the Schengen area. It also includes closures of schools and universities as well as all non-essential public places, including shops (except for food shopping), restaurants, cafés, cinemas, and nightclubs.

The basic reproduction number R0 corresponds to the expected number of new cases generated by a single infectious case in a fully susceptible population (3). Several studies, mostly based on Chinese data, aimed at estimating the R0 associated with the COVID-19 epidemic, leading to values from 1.4 to 6.49, with an average of 3.28 (4). As the value of R0 can be interpreted as the product of the contact rate and of the duration of the infectious period, and since the objective of the lockdown and associated restriction strategies are precisely to decrease the contact rate, an important effect on the number R of secondary cases generated by an infectious individual is to be expected. This value R is often referred to as “effective reproduction number,” and corresponds to the counterpart of R0 in a population that is not fully susceptible (5). If R > 1, the number of infectious cases in the population follows an increasing trend, and the larger R, the faster this trend. On the contrary, if R < 1, the epidemic will gradually die out. The control measures in China have been shown to have a significant effect on the COVID-19 epidemic, with growth rates that shifted from positive to negative values (corresponding to R < 1) within 2 weeks (6). The study (7) showed that containment policies in Hubei province also led to a subexponential growth in the number of cases, consistent with a decrease in the effective reproduction number R. Fitting a SEIR epidemic model to time series of reported cases from 31 provinces in China, Tian et al. (8) found a basic reproductive number R0 = 3.15 before the implementation of the emergency response in China, a value that was divided by more than 20 once the control measures were fully effective. Using contact surveys data for Wuhan and Shanghai it was estimated in Zhang et al. (9) that the effective reproduction number was divided by a factor 7 in Wuhan and 11.5 in Shanghai.

Standard epidemiological models generally rely on SIR (Susceptible-Infected-Removed) systems of ordinary differential equations and their extensions [for examples of application to the COVID-19 epidemic, see (10, 11)]. With these models, and more generally for most deterministic models based on differential equations, when the loss of information due to the observation process is heavy, specific approaches have to be used to bridge the gap between the models and the data. One of these approaches is based on the mechanistic-statistical formalism, which uses a probabilistic model to connect the data collection process and the latent variable described by the ODE model. Milestone articles and textbook have been written about this approach or related approaches (12), which is becoming standard in ecology (13, 14). The application of this approach to human epidemiological data is still rare.

In a previous study (15), we applied this framework to the data corresponding to the beginning of the epidemic in France (from February 29 to March 17), with a SIR model. Our primary objective was to assess the infection fatality ratio (IFR), defined as the number of deaths divided by the number of infected cases. As the number of people that have been infected is not known, this quantity cannot be directly measured, even now (on April 15). The mechanistic-statistical framework allowed us to compute an IFR of 0.8% (95%-CI: 0.45–1.25%), which was consistent with previous findings in China (0.66%) and in the UK (0.9%) (16) and lower than the value previously computed on the Diamond Princess cruse ship data (1.3%) (17). In this previous study, we also computed the R0 in France, and we found a value of 3.2 (95%-CI: 3.1–3.3). Although the number of tests at that stage was low, an advantage of working with the data from the beginning of the epidemic was that the initial state of the epidemic was known.

Here, we develop a new mechanistic-statistical approach, based on a SIRD model (D being the dead cases compartment), in the aim of

estimating the effect of the lockdown in France on the contact rate and the effective reproduction number R;

estimating the number of infectious individuals and the fraction of the population that has been infected by the beginning of May (at the official date at which the lockdown should be relaxed).

2. Materials and Methods

2.1. Data

We obtained the number of positive cases and deaths in France, day by day from Santé Publique France (18), from March 31 to April 14. We obtained weekly data on the number of individuals tested (in private laboratories and hospitals) from the same source. We assumed that during each of these weeks the number of tests per day was constant. This assumption is consistent with the small variations between the number of tests during the first week (111,690) and the second week of observation (132,392). As the data on the number of positive cases are not fully reliable (fewer cases during weekends with a rebound on Monday), we smoothed the data with a moving average over 5 days. Official data on the number of deaths by COVID-19 since the beginning of the epidemic in France only take into account hospitalized people. About 728, 000 people in France live in nursing homes [EHPAD, source: DREES (19)]. The number of deaths in these structures has only been reported recently, and cannot be obtained day by day. Latest data from Santé Publique France indicate a total number of 10, 643 deaths at hospital and 6, 524 deaths in nursing homes by April 15. The total number of deaths therefore corresponds to about 1.6 times the number of deaths at hospital. The same factor had been estimated in Roques et al. (15) based on local dataset in the French Grand Est region.

2.2. Mechanistic-Statistical Framework

The mechanistic-statistical framework consists in the combination of a mechanistic model that describes the epidemiological process, a probabilistic observation model and an inference procedure.

2.2.1. Mechanistic Model

The dynamics of the epidemic are described by the following SIRD compartmental model:

with S the susceptible population, I the infectious population, R the recovered population, D the number of deaths due to the epidemic and N the total population. For simplicity, we assume that N is constant, equal to the current French population, thereby neglecting the effect of the small variations of the population on the coefficient α/N. The parameter α is the contact rate (to be estimated) and 1/β is the mean time until an infectious becomes recovered. Based on the results in Zhou et al. (20), the median period of viral shedding is 20 days, but the infectiousness tends to decay before the end of this period: the results in He et al. (21) indicate that infectiousness starts 2–3 days before symptom onset and declines significantly 8 days after symptom onset. Based on these observations we assume here that the mean duration of the infectiousness period is 1/β = 10 days. In Li et al. (22), the duration of the incubation period was estimated to have a mean of 5.2 days. Thus, the mean duration of the non-infectious exposed period is relatively short (about 2–3 days), and can be neglected without much differences on the results, as shown in Liu et al. (23). Inclusion of an exposed compartment (as in SEIR models) is particularly relevant when exposed individuals can indirectly transmit the disease e.g., through insect vectors [e.g., (24)], which is seemingly not the case for coronaviruses. The parameter γ corresponds to the death rate of the infectious (to be estimated).

2.2.1.1. Initial conditions

The model is started at a date t0 corresponding to April 1st. The initial number of infectious I(t0) = I0 is not known and will be estimated. The total number of recovered at time t0 is also not known. However, as the compartment R has no feedback on the other compartments, we may assume without loss of generality that R(t0) = 0, thereby considering only the new recovered individuals, starting from the date t0. We fixed D(t0) = 3523, the number of deaths at hospital by March 31. The initial S population at the beginning of the period, should still be close to the total French population: by March 31 only 52,128 cases had been observed in France, corresponding to 0.08% of the total population. A factor 8 had been estimated in Roques et al. (15) between the cumulated number of observed cases and the actual number of cases at the beginning of the epidemic. Even though this factor may have changed, this means that the proportion of the total population that has been infected by March 31 is still small. We can get an upper bound for the cumulated number of cases by March 31 by dividing the number of hospital deaths at the end of the observation period (10,129 by April 14) by the hospital IFR [0.5%, as estimated in (15)] leading to about 2 million cases. This means that the value of S(t0) is between 65 and 67 million cases. For our computation, we assumed that , corresponding to about 98.5% of the French population. As shown in Figure S3, our results are not much sensitive to the value of S(t0) (at least when S/N remains close to 1).

2.2.1.2. Numerical method

The ODE system (1) was solved thanks to a standard numerical algorithm, using Matlab® ode45 solver.

2.2.2. Observation Model

The number of cases tested positive on day t, denoted by , is modeled by independent binomial laws, conditionally on the number of tests n carried out on day t, and on p the probability of being tested positive in this sample:

The tested population consists of a fraction of the infectious cases and a fraction of the susceptibles: n = τ1(t)I(t)+τ2(t)S(t). Thus,

with κ: = τ2(t)/τ1(t), the relative probability of undergoing a screening test for an individual of type S vs an individual of type I. We assumed that the ratio κ was independent of t over the observation period. The coefficient σ corresponds to the sensitivity of the test. In most cases, RT-PCR tests have been used and existing data indicate that the sensitivity of this test using pharyngeal and nasal swabs is about 63−72% (25). We assumed here σ = 0.7 (70% sensitivity).

Each day, the number of new observed deaths (excluding nursing homes), denoted by , is modeled by independent Poisson distributions conditionally on the process D(t), with mean value D(t) − D(t − 1) (which measures the daily increment in the number of deaths):

Note that the time t in (1) is a continuous variable, while the observations and are reported at discrete times. For the sake of simplicity, we used the same notation t for the days in both the discrete and continuous cases. In the formulas (2) and (3) I(t), S(t), and D(t) are computed at the end of day t.

2.2.3. Statistical Inference

The unknown parameters are α, γ, κ, and I0. We used a Bayesian method (26) to estimate the posterior distribution of these parameters.

2.2.3.1. Computation of the likelihood function

The likelihood is defined as the probability of the observations (here, the increments ) conditionally on the parameters. Using the observation models (2) and (3), and using the assumption that the increments and are independent conditionally on the underlying SIRD process and that the number of tests n is known, we get:

with t the date of the first observation and t the date of the last observation. In this expression depends on α, γ, κ, I0 through p and D(t).

2.2.3.2. Posterior distribution

The posterior distribution corresponds to the distribution of the parameters conditionally on the observations:

where π(α, γ, κ, I0) corresponds to the prior distribution of the parameters (detailed below) and C is a normalization constant independent of the parameters.

2.2.3.3. Prior distribution

Regarding the contact rate α, the initial number of infectious cases I0 and the probability κ, we used independent non-informative uniform prior distributions in the intervals α ∈ (0, 1), and κ ∈ (0, 1). To overcome identifiability issues, we used an informative prior distribution for γ. This distribution, say f, was obtained in Roques et al. (15) during the early stage of the epidemic (f is depicted in Figure S1). In Roques et al. (15), the number of infectious cases I0 at the beginning of the epidemic was known (equal to 1), and did not need to be estimated. Thus, we estimated in Roques et al. (15) the distribution of the parameter γ by computing the distribution of the infectious class and using the formula D′(t) = γI(t) together with mortality data (which were not used for the estimation of the other parameters, unlike in the present study). Finally, the prior distribution is defined as follows:

The numerical computation of the posterior distribution is performed with a Metropolis-Hastings (MCMC) algorithm, using 5 independent chains, each of which with 106 iterations, starting from the posterior mode. To find the posterior mode we used the BFGS constrained minimization algorithm, applied to , via the Matlab® function fmincon. In order to find a global minimum, we applied this method starting from 4,000 random initial values. The Matlab® codes are available as Supplementary Material.

3. Results

3.1. Model Fit

Denote by the posterior mode, and S*(t), I*(t), R*(t), D*(t) the solutions of the system (1) associated with these parameter values. The observation model (2) implies that the associated expected number of cases tested positive on day t is (expectation of a binomial) with

The observation model (3) implies that the expected cumulated number of deaths on day t is D*(t).

To assess model fit, we compared these expectations and the observations, i.e., the cumulated number of cases tested positive, with C0 the number of cases tested positive by March 31 (C0 = 52, 128) and the cumulated number of deaths , with M0 the number of reported deaths (at hospital) by March 31 (M0 = 3123). The results are presented in Figure 1. We observe a good match with the data.

Figure 1

Expected number of observed cases and deaths associated with the posterior mode vs. number of cases actually detected (total cases). The blue curve corresponds to the expected number of cases tested positive given by the model, the red curve corresponds to the expected cumulated number of deaths D*(t) (excluding nursing homes). The crosses correspond to the observations (blue crosses: cumulated number of positive cases, red crosses: cumulated number of deaths). C0 is the number of cases tested positive on March 31 (C0 = 52 128).

The pairwise posterior distributions of the parameters (α, I0), (α, γ), (α, κ), (γ, I0), (γ, κ), (κ, I0) are depicted as Figure S2. With the exception of the parameter γ (Figure S1), for which we chose an informative prior, the posterior distribution is clearly different from the prior distribution, showing that new information was indeed contained in the data.

3.2. Contact Rate and Effective Reproduction Number

The effective reproduction number can be simply derived from the relation R = α/(β + γ) when S is close to N (3). The distribution of R is therefore easily derived from the marginal posterior distribution of the contact rate α (since we assumed β = 1/10; see section 2.2). It is depicted in Figure 2. We observe a mean value of R of 0.47 (95%-CI: 0.45–0.50).

Figure 2

Posterior distribution of the effective reproduction number R in France.

3.3. Dynamics of the Infectious Class

The marginal posterior distribution of I0 indicates that the number of infectious individuals at the beginning of the considered period (i.e., April 1st) is 1.4·106 (95%-CI: 1.1·106 − 1.8·106). The computation of the solution of (1) with the posterior distribution of the parameters leads to a number of infectious and a total number of infected cases (including recovered) at the end of the observation period (April 14). By May 10, if the restriction policies remain unchanged, we get a forecast of I(T) = 1.6·105 infectious cases (95%-CI: 1.3·105 − 2.1·105) and (I + R)(T) = 2.5·106 infected cases including recovered (95%-CI: 2.0·106 − 3.2·106). The dynamics of the distributions of I and I + R are depicted in Figure 3. By May 10, the total number of infected cases (including recovered) therefore corresponds to a fraction of 3.7% of the total French population. This value does not include the recovered cases before April 1st.

Figure 3

Distribution of the number of infectious cases I(t) and cumulated number of infected cases I(t) + R(t) across time. Solid lines: average value obtained from the posterior distribution of the parameters. Shaded areas: 0.025–0.975 interquantile ranges.

4. Discussion

Many studies focused on the estimation of the basic reproductive number R0 of the COVID-19 epidemic, based on data-driven methods and mathematical models [e.g., (4, 27)] describing the epidemic from its beginning. In average, the estimated value of R0 was about 3.3. We focused here on an observation period that began after the lockdown was set in France.

We obtained an effective reproduction number that was divided by a factor 7, compared to the estimate of the R0 carried out in France at the early stage of the epidemic, before the country went into lockdown [a value R0 = 3.2 was obtained in (15)]. This indicates that the restriction policies were very efficient in decreasing the contact rate and therefore the number of infectious cases. In particular, the value R = 0.47 is significantly below the threshold value 1 were the epidemic starts dying out.

The decay in the number of infectious cases can also be observed from our simulations. It has to be noted that, although the number of infectious cases is a latent, or “unobserved” process, the mechanistic-statistical framework allowed us to estimate its value (Figure 3). The cumulated number of infected cases that we obtained by May 10 (I + R) corresponds to a fraction of 3.7% (95%-CI: 3.0–4.8%) of the total French population, without taking into account the number of recovered individuals before April 1st, which is not known. Based on a value R0 = 3.2, the herd immunity threshold, corresponding to the minimum fraction of the population that must have immunity to stop the epidemic, would be 1 − 1/R0≈69% [a threshold of 80% was proposed in (28)]. This proportion will probably not be reached by May 10. As emphasized by Angot (29), a too fast relaxation of the lockdown-related restrictions before herd immunity is reached or efficient prophylaxis is developed), would expose the population to an uncontrolled second wave of infection. In the worst-case scenario, the effective reproduction number R would approach the initially estimated value of R0, and the second wave would start with about 1.6·105 infectious individuals (in comparison with the few cases that initiated the first wave in France) and about 64·106 susceptible individuals. Keeping a low value of R is therefore crucial to avoid the saturation of hospital facilities.

We deliberately chose a parsimonious mechanistic model with a few parameters to avoid identifiability issues. Possible extensions include stage-structured models, where the infectious class I and the contact rate α would depend on another variable: I = I(t, τ) and α = α(t, τ) with τ the time since infection, to take into account the dynamics of the viral load on the infectiousness. See e.g., Murray (3) (chapter 19.6) for an introduction to such modeling approaches. Another insightful extension would consist in using spatially-explicit models, e.g. reaction-diffusion models (30) to describe the spatial spread of the epidemic, and to be able to estimate local values for the parameter R and the number of susceptible cases. Although herd immunity is far from being reached at the country scale, it is likely that the fraction of immune individuals strongly varies over the territory, with possible local immunity effects [e.g., by April 4 the proportion of people with confirmed SARS-CoV-2 infection based on antibody detection was of 41% in a high-school located in Northern France (31)].

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://www.gouvernement.fr/info-coronavirus/carte-et-donnees https://geodes.santepubliquefrance.fr and https://ourworldindata.org/coronavirus-testing.

Author Contributions

LR, EK, JP, AS, and SS conceived the model and designed the statistical analysis. LR and SS wrote the paper. LR carried out the numerical computations. All authors reviewed the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.