Mechanistic-statistical SIR modelling for early estimation of the actual number of cases and mortality rate from COVID-19

Data

The COVID-19 screening data in France and South Korea used for preparing this note are those available from January 22, 2020 to March 17, 2020. These data describe the number of positive cases and deaths, day by day (source: Johns Hopkins University Center for Systems Science and Engineering, https://github.com/CSSEGISandData/COVID-19). The number of tests carried out is only available from February 22 (Sources: Santé publique France and Korean Center for Disease Control). As some data (positive cases, deaths) are not fully reliable (example: 0 new cases detected in France on March 12, 2020), we smoothed the data with a moving average over 5 days.

Mechanistic model

SIR models are the most standard ODE-based epidemiological models. These models are compartmental: they divide the population into subclasses: the susceptibles (S), the infected (I) and the recovered (R, immune individuals in our case). The simplest SIR compartmental model does not take into account the demography of the S compartment: with N = S + I + R the total population, which is constant. The impact of the compartment D (number of dead people) on the dynamics of the SIR system is therefore neglected here. The compartment D satisfies:

We use this equation later to compute the death rate due to the disease. The initial conditions are S(t₀) = N − 1, I(t₀) = 1 and R(t₀) = 0, where N corresponds to the population size in France or South Korea (67·10⁶ and 52·10⁶ people). The SIR model is started at t = t₀, which should approach the date of introduction of the virus in the considered country (this point is shortly discussed at the end of this note).

Note that I′(t) = β I (R₀ S/N − 1), with R₀ = α/β the basic reproduction rate (Murray, 2002). When R₀ < 1, we observe that I′ < 0, and so the epidemic cannot spread in the population. When R₀ > 1, the infected compartment I increases as long as R₀ S > N = S + I + R.

The system (1) can be analytically solved, using a change of time variable which requires a numerical integration. Here, we rather solve the ODE system thanks to a standard numerical algorithm, using Matlab® ode23s solver.

Observation model

We denote by the number of cases tested positive on day t. We suppose that these increments follow independent binomial laws, conditionally on the number of tests, on I(t) and on S(t): where n_t corresponds to the number of tests carried out on day t and p_t the probability of being tested positive in this sample. The tested population consists of a fraction of the infecteds and a fraction of the susceptibles: n_t = τ₁(t) I(t) + τ₂(t) S(t). Thus, with κ := τ₂(t)/τ₁(t), the relative probability of undergoing a screening test for an individual of type S vs an individual of type I (probability of being tested conditionally on being S / probability of being tested conditionally on being I). We assume that the ratio κ does not depend on t at the beginning of the epidemic (i.e., over the period that we use to estimate the parameters of the model). The daily number of deaths caused by the virus is assumed to be known exactly.

In order to compute the maximum likelihood estimator (the MLE, i.e., the parameters that maximize ℒ), we use the BFGS constrained minimization method, applied to − ln(ℒ), via the Matlab ® function fmincon. In order to find a global maximum of ℒ, we apply this method starting from random initial values for α, t₀, κ drawn uniformly in the following intervals:

For each country, the minimization algorithm is applied to 2000 random initial values of the parameters.

The posterior distribution of the parameters (α, t₀, κ) is computed with a Bayesian method, using uniform prior distributions in the intervals given by (4). This posterior distribution corresponds to the distribution of the parameters conditionally on the observations: where π(α, t₀, κ) corresponds to the prior distribution of the parameters (therefore uniform) and C is a normalization constant independent of the parameters. The numerical computation of the posterior distribution (which is only carried out for French data) is performed with a Metropolis-Hastings (MCMC) algorithm, using 5 independent chains, each of which with 10⁶ iterations, starting from random values close to the MLE.

Unless otherwise mentioned, the data used to compute the MLE and the posterior distribution are those corresponding to the period from February 29 to March 17.

Parameter distribution

The joint posterior distributions of the three pairs of parameters (α, κ), (t₀, α) and (t₀, κ) for France are presented in Appendix A (Fig. 6). Note that the distributions are very different from the uniform prior distribution. However, the distributions of t₀ and κ are quite dispersed. The joint distribution of (t₀, κ), presented in Appendix A shows a correlation between t₀ and κ. Based on this distribution, in order to reduce the uncertainty on κ, we assume that t₀ is between January 13 and 30.

• Download figure
• Open in new tab

Figure 6:

Joint posterior distributions of (α, κ), (t₀, α) and (t₀, κ) in France.

Fig. 3 depicts the marginal posterior distribution of the basic reproduction rate R₀. The value of R₀ corresponding to the MLE in France is . A similar computation in South Korea, based on the data used in Fig. 2 gives .

• Download figure
• Open in new tab

Figure 3:

Posterior distribution of the basic reproduction rate R₀ in France.

Actual number of infected cases

Using the posterior distribution of the model parameters, with the constraint that t₀ is between January 13 and 30, we can compute the daily distribution of the actual number of infected peoples. This distribution is represented across time in Fig. 4. We deduce the following ratios between the actual number of infected people and the observations, (with the sum of the observed infected cases at time t). Thus, in France, the estimated ratio between the actual number of infected and observed cases is 15 (95%-CI: 4-33).

• Download figure
• Open in new tab

Figure 4: Distribution of the actual number of infected cases in France across time.

Solid line: average value obtained from the posterior distribution of the parameters. Dotted curves: 0.025 and 0.975 pointwise posterior quantiles.

Actual mortality rate

The mortality rate corresponds to the fraction of the infected who die, that is γ(t)/(γ(t) + β). The term γ(t) is computed using the formula (2) and the mortality data. We thus obtain, on March 17, a mortality rate in France of 5.2/1000 (95%-CI: 1.5/1000 − 11.7/1000). The temporal dynamics of the mortality rate are depicted in Fig. 5.

• Download figure
• Open in new tab

Figure 5: Dynamics of the mortality rate in France.

Solid line: average value obtained from the posterior distribution of the parameters. Dotted curves: 0.025 and 0.975 pointwise quantiles.

On the number of infecteds and the mortality rate

The actual number of infected individuals in France is probably much higher than the observations (we find here a factor × 15), which leads at a lower mortality rate than that calculated on the basis of the observed cases. However, if the virus led to contaminate 80% of the French population (Ferguson et al., 2020), the total number of deaths to deplore in the absence of variation in the mortality rate (increase induced for example by the saturation of hospital structures, or decrease linked to better patient care) would be 277, 000 (95%-CI : 81, 000 − 629, 000)). This estimate could be corroborated or invalidated when 80% of the population will be infected, eventually over several years, assuming that an infected individual is definitively immunized. It has to be noted that measures of confinement or social distancing can decrease both the percentage of infected individuals in the population and the degree of saturation of hospital structures.

On the differences between France and South Korea

The mechanistic-statistical SIR model achieves a satisfactory goodness-of-fit for French data, but does not capture the decline in the number of cases observed in South Korea. Grounding the analysis of the South Korean situation on the first phase of the epidemic improves the goodness-of-fit and yields an estimated R₀ twice as low as the French R₀, resulting on a slower epidemic dynamic. The difference between the dynamics predicted by the SIR model and the South Korean data is probably linked to a different management of the epidemic in Korea, having a strong impact on the epidemic dynamics (more important screening, tracing, social distancing in South Korea).

On the value of R₀

The value of R₀ obtained in South Korea is consistent with recent estimates for COVID-19 (2.0,2.6); see Ferguson et al. (2020). The estimated distribution in France is therefore surprisingly high. This difference could be due to a different definition of R₀ depending on the type of model used to calculate it. A direct estimate, by a non-mechanistic method, of the parameters (ρ, t₀) of a model of the form gives t₀ = 30 (January 30) and ρ = 0.19. With the SIR model, I′(t) ≈ I (α − β) for small times (S ≈ N), which leads to a growth rate equal to ρ ≈ α − β, and a value of α ≈ 0.24, that is to say R₀ = 4.8, which is consistent with the distribution presented in Fig. 3. Note that β = 1/20 corresponds to the median period of viral shedding of 20 days described by Zhou et al. (2020). A shorter period would lead to a lower value of R₀.

On the uncertainty linked to the data

The uncertainty on the actual number of infected and therefore the mortality rate are very high. We must therefore interpret with caution the inferences that can be made based on the data we currently have in France. In addition, we do not draw forecasts here: the future dynamics will be strongly influenced by the containment measures that will be taken and should be modeled accordingly.

On the hypotheses underlying the model

The data used here contain a limited amount of information, especially since the observation period considered is short and corresponds to the initial phase of the epidemic dynamics, which can be strongly influenced by discrete events. This limit led us to use a particularly parsimonious model in order to avoid problems of identifiability for the parameters. The assumptions underlying the model are therefore relatively simple and the results must be interpreted with regard to these assumptions. For instance, the date of the introduction t₀ must be seen as an efficient date of introduction for a dynamics where a single introduction would be decisive for the outbreak and the other (anterior and posterior) introductions would have an insignificant effect on the dynamics.