Disclaimer

This text is for research and educational purposes only not a decision-making advice. The following calculations are made using voluntarily simplified assumptions.

For any question, please refer to the French public health authority, Santé Publique France, where one can find all official instructions: https://www.gouvernement.fr/info-coronavirus

The World Health Organisation (WHO) provides as well a comprehensive website: https://www.who.int/fr/emergencies/diseases/novel-coronavirus-2019

Context

On March 11th, 2020, the WHO announced that the COVID-19 outbreak had reached pandemic stage. On March 12th, the French authorities announced the closure of school structures starting from March 16th. On March 14th, the Prime Minister announced the transition to the third stage of the fight plan against the epidemic and the closure of most of the country’s infrastructure.

Several research groups published more or less detailed versions of mathematical models of the spread of the epidemic online. The structure of the deterministic mathematical model described below includes aspects of some of these models. The internet sites (most of them in English) are referenced at the end of the note.

This model is more complex than the classical (SIR) model, but its structure remains very simple. Indeed, all individuals are assumed to react on average in the same way to infection (there are no differences in age, sex, contacts). Moreover, there is no spatial structure in the model where everyone is potentially in contact with everyone. In summary, at best this model captures the dynamics of the epidemic in a city in France with limited differences across ages.

In this note, we present the model and its hypotheses. Then we calculate the expression of the basic reproduction number (\(R_0\)) as a function of the model parameters. Finally, we perform intervention simulations using default parameter values and compare an epidemic suppression strategy (intense but short control) with a mitigation strategy (lighter but sustainable control).

Life cycle of the infection

Simple model

COVID-19 is a respiratory disease that is transmitted both directly (from person to person) and indirectly (by viruses present in the environment shed by infected people). The latter route of transmission is neglected here. Susceptible people (noted \(S\)) who are infected go through a stage where they are infected but they are neither infectious nor symptomatic (noted \(E\) for exposed). Then they become infectious while remaining asymptomatic (noted \(A\) for asymptomatic). This stage appears to be particularly important in the spread of COVID-19. Symptoms then develop, with the individuals remaining infectious (denoted \(I\) for infectious). Finally, the infection ends when people are either immune or dead (denoted \(R\) for removed). This life cycle can be represented using the following flow chart:

Simple COVID-19 life cycle.

In order to simulate epidemics, it is necessary to know the rate at which exposed people become infectious (\(\epsilon\)), the rate at which infectious people become symptomatic (\(\sigma\)), the rate at which symptomatic people cease to be infectious (\(\gamma\)), the rate of transmission from asymptomatic people (\(\beta_A\)), and the rate of transmission from symptomatic people (\(\beta_I\)). Finally, parameter \(c\) corresponds to the intensity of control (percentage of infected people who are confined). On this diagram, we have not represented the arrival of new hosts of type E by migration (at a rate \(\nu\)).

Model with differential disease severity

Using existing data, this model can be made more complex without adding too many unknown parameters. In practice, we will assume that a fraction \(p\) of infections are mild and a fraction \(1-p\) are severe and require hospitalization. This translates into the life cycle below where mild cases are indexed with a “1” and severe cases with a “2”.

COVID-19 life cycle with different levels of case severity.

This model requires some additional parameters: \(\alpha\) is the mortality rate for severe cases, \(\gamma_2\) is the recovery rate for severe cases (which may differ from that for mild cases, \(\gamma_1\)) and potentially different transmission rates from asymptomatic or mild (\(\beta_1\)) and from severe (\(\beta_2\)) cases, since here we assume that asymptomatic and symptomatic cases of a given severity degree have identical transmission rates. Finally, on the diagram, \(\lambda\) is the strength of infection, which is the rate at which susceptible individuals become infected. It is defined by the formula \[\begin{align} \lambda = \left[ \beta_1~(A_1+I_1+A_2) + \beta_2~I_2 \right] (1-c). \end{align}\] For the sake of clarity we do not show transmission and migration events on the diagram.

(Rather) technical point : some models incorporate the bifurcation between mild and severe cases after either stage \(E\), or stage \(A\) or even stage \(I\). This is problematic because, for very high values of p (the proportion of mild cases), one finds very low transition rates to severe cases. The problem is that this means that these events occur on average later than events with high rates. However, the fact that an event is rare (e.g. dying from COVID-19), does not mean that it occurs very late in the course of the infection. To avoid this issue we made the waiting time depend on an infection event (passage from \(S\) to \(E\)). Indeed, it is possible that “by chance” a susceptible person is infected very quickly (after one day) or very late (in two months), whereas for a duration of infection this is biologically unrealistic (you don’t heal in 1 day or 2 months). End of the technical parenthesis.

Calculation of \(R_0\)

One of the strengths of the \(R_0\) concept is that it can be calculated from the life cycle of the infectious agent. For our model with differential severity of infections, we can thus show using the Next Generation theorem (Diekmann et alii 1990 J. Math. Biol., Hurford et alii 2010 J R Soc Interface) that, if natural mortality is neglected, then \[\begin{align} R_0 = \frac{\beta_1~S_0}{\sigma} + p~\frac{\beta_1~S_0}{\gamma_1} + (1 - p)~\frac{\beta_2~S_0}{\alpha + \gamma_2} . \end{align}\]

The three terms of this expression have an intuitive interpretation because they correspond to a fraction of the secondary infections generated by an infected person during the course of his or her infection (which is the definition of \(R_0\)). The first term corresponds to infections caused by asymptomatic people \(A_1\) and \(A_2\) (in fact, \(1/\sigma\) is the average time spent in the asymptomatic state). The second term corresponds to infections caused by symptomatic mild infections \(I_1\) (the duration of this phase can be found through \(1/\gamma_1\) as well as the proportion of mild cases \(p\)). Finally, the last term corresponds to secondary infections caused by severe and symptomatic cases \(I_2\).

With this expression, we can see that there are several ways to lower \(R_0\) and thus control the epidemic. For example, the number of susceptible people and the rates of transmission can be reduced (by confining the population, reducing contacts and wearing masks). One can also increase the recovery rates by isolating symptomatic people (from the point of view of the virus, a person who no longer has contacts is equivalent to a recovery, since he or she no longer transmits the virus).

Parameter ranges

The following table lists the model parameters, the default values chosen and the biological ranges listed. These ranges are often wide due to large variations from one infection to another. For more details, please refer to the MIDAS consortium compilations.

Notation Description Default Range Reference
\(\epsilon\) rate of end of latency \(1/4.2\ \text{jours}^{-1}\) \([0.21\ ;\ 0.27]\) Li et alii (2020, NEJM)
\(\sigma\) rate of symptoms onset \(1\ \text{jour}^{-1}\) \([0.9\ ;\ 1.1]\) Ferguson et alii (2020, Rapport 8)
\(\gamma_1\) recovery rate of mild cases \(1/17\ \text{jours}^{-1}\) \([0.025\ ;\ 0.1]\) Zhou et alii (2020, Lancet)
\(\gamma_2\) recovery rate of severe cases \(1/20\ \text{jours}^{-1}\) \([0.039\ ;\ 0.13]\) Zhou et alii (2020, Lancet)
\(R_0\) basic reproduction number \(2.5\) \([2\ ;\ 3]\) our report
\(p\) proportion of infections that do not require hospitalization \(0.9\) \([0.85\ ;\ 0.95]\) Santé Publique France
\(\theta\) case fatality ratio of hospital patients \(0.15\) \([0.135\ ;\ 0.165]\) Santé Publique France
\(\alpha\) mortality rate of severe infections calculé — see the formula
\(\nu\) migration rate \(10^{-6}\ \text{jour}^{-1}\) — ad hoc
\(b\) decrease in transmission due to hospitalization \(0.2\) — ad hoc
\(c\) fraction of the \(R_0\) decreased by public health policies — — variable

Note that the mortality rate \(\alpha\) can be expressed as a function of the case-fatality rate of severe cases through the formula \[\begin{align} \alpha & =\gamma_2~\frac{\theta}{1-\theta}. \end{align}\]

In addition, the transmission rate \(\beta\) is parameterized indirectly using the values of \(R_0\) (between 2 and 3 depending on the country, with values close to 2.5 for France). Reworking the previous equation, we find that \[\begin{align} \beta & = \frac{R_0}{S_0} \gamma_1 \sigma \frac{\alpha + \gamma_2}{(\gamma_1+p \sigma)(\alpha + \gamma_2) + b \gamma_1 \sigma (1- p)} . \end{align}\]

Group immunity and epidemic size

As a reminder, as mentioned in a previous report, the percentage \(P^\star\) of the population that needs to be immunized (i.e. in compartment \(R\)) in order to prevent a future epidemic is: \[\begin{align*} P^\star & = 1 - \frac{1}{R_0}. \end{align*}\]

Furthermore, in the absence of any control measures, the total fraction of the population \(Q^\star\) that will be affected by the epidemic is found by solving the equation \[\begin{align*} Q^\star~R_0 + \log(1-Q^\star) & = 0. \end{align*}\]

Numerical simulations

In these simulations, we normalized the total population size and represent the percentage of the population infected. This is voluntary because, as noted above, this model has no structure and therefore cannot be applied directly to an entire country but rather to a small town or village.

Epidemic wave without control policy

In the absence of any control, the epidemic is growing exponentially over time. On the following graph, each colour corresponds to one of the life cycle categories. The grey dashed line corresponds to the fraction of the population that can be accommodated in resuscitation or intensive care (approximately 12,000 beds in total, which means \(1.79E-4\) per inhabitant).

In this figure, we can see that the majority of hosts are in the \(I_1\) state. This makes sense because the majority of infected people do not have a severe infection requiring hospitalization and because this is also the life-cycle stage with the longest duration. Finally, we also note that the proportion of infected people requiring hospitalization (\(I_2\), the line in turquoise blue) quickly exceeds the proportion of hospital beds (dashed line).

If we look further back in time, we see (see figure below) that in the absence of any intervention, the epidemic reaches its peak 150 days after the beginning. This corresponds approximately to the point where the curves for those still susceptible (\(S\) in yellow) and those immunized following recovery (\(R\) in blue) cross. Finally, we see, as indicated in our previous report, that the epidemic does not stop once the threshold for herd immunity (materialized by the black dashed line) is reached.

Because of the simplicity of the assumptions made in the model, we prefer not to discuss the mortality associated with this scenario. But the fact that mortality is visible on this graph (in red) which represents the French population suggests that it would be not negligible (in addition we have neglected here the additional indirect mortality due to the saturation of the health system).

Epidemic wave with strong control but of limited duration

We then modelled the implementation of a public health policy aimed at strongly cutting off the transmission of the virus (here \(c=90\%\)). This type of policy is called “suppression” or “containment”. Here we have started this policy on the 60th day of the outbreak and for a duration of 30 days. This is shown by the shaded rectangle in the figure below.

As expected, the epidemic is suppressed during the implementation period of the control policy. It can be seen that the curve for people requiring intensive care (\(I_2\) in turquoise) remains well below the threshold value for the number of available hospital beds (grey dashed lines). Unfortunately, once the containment measures cease, the epidemic starts again very strongly and the number of available hospital beds is quickly exceeded.

In the end, the policy results in shifting the epidemic peak (the curves intersect at 200 days instead of 150). However, the total number of people affected remains the same (the blue curve for people who recovered reaches a similar value).

Epidemic wave with moderate but persistent control

Another approach, championed by some countries, is to slow the growth of the epidemic but not to stop it. This is called mitigation or attenuation. This approach can be modelled by e.g. limiting only \(c=40 \%\) of virus transmission for 365 days. As shown in the figure below, this strategy does not prevent the spread of the epidemic. It can also be seen that the number of severe \(I_2\) cases is rapidly rising above the threshold of the proportion of available intensive care and resuscitation beds.

In the end, the course of the epidemic is very similar to that without control, but we see that the epidemic peak occurs much later (after 300 days). Moreover, the total proportion of the population infected is much lower than in the absence of control. We can also see that we are close to the threshold of herd immunity (black dashed line). However, we are below the threshold at the time the measures cease and a new small epidemic (the number of susceptible people is decreasing) is observed.

Limitations and Future Studies

The purpose of this report is to illustrate the structure of a compartmentalized epidemiological model. We also wanted to illustrate the effects of different policies rather than trying to identify an optimal policy (which a priori would combine suppression and mitigation). Finally, we made several simplistic assumptions that will be released in models in preparation. In particular:

  • we have chosen a value for each of the parameters of our model but these vary a priori in ranges that should be explored,

  • this model has not been calibrated to real data (notably the dates of the first cases and the records of the number of deaths); thus temporal aspects such as the date of the epidemic peak have only relative value,

  • we have neglected the indirect transmission through fomites (door handles, fabrics, money,…), which could have an important role,

  • the population is not age structured though mortality from COVID-19 infection is known to increase with age,

  • we have assumed that only people with a severe infection (\(I_2\)) die, whereas in reality, the saturation of the hospital system risks increasing the mortality of these people and also increasing the mortality of uninfected people (\(S\) and \(R\)),

  • control strategies are applied in a uniform manner here, whereas in principle they vary over time,

  • the delays in models with continuous time differential equations are average and this can greatly affect the temporality of the dynamics.

Technical appendix

Model equations with differential severity

We present here a model without demographics (so births and natural deaths are neglected): \[\begin{align} \renewcommand{\d}{\mbox{d} } \newcommand\dt[1]{\frac{\d{#1}}{\d t}} \dt{S} & = - \lambda~S \\ \dt{E_1} & = p~\lambda~S -\epsilon~E_1 + p~\nu \\ \dt{A_1} & = \epsilon~E_1 - \sigma~A_1 \\ \dt{I_1} & = \sigma~A_1 - \gamma_1~I_1 \\ \dt{R_1} & = \gamma_1~I_1 \\ \dt{E_2} & = (1-p)~\lambda~S -\epsilon~E_2 + (1-p)~\nu \\ \dt{A_2} & = \epsilon~E_2 - \sigma~A_2 \\ \dt{I_2} & = \sigma~A_2 - \left(\gamma_2 +\alpha \right)~I_2 \\ \dt{R_2} & = \gamma_2~I_2 \\ \dt{M} & = \alpha~I_2 \nonumber \\ \text{with } N & = S+ \sum_j \left( E_j+A_j+I_j+R_j \right) \nonumber \\ \text{and } \lambda & = (1-c) \left(\beta_1~A_1 + \beta_1~I_1 + \beta_1~A_2 + \beta_2~I_2 \right) \nonumber \end{align}\]

Calculation of \(R_0\)

Our goal is to perturb a system that is initially free of infections, thus where the variables \((E_1,A_1,I_1,E_2,A_2,I_2)\) are equal to 0.

The Jacobian matrix of the system of differential equations listed above and evaluated in its initial state is: \[\begin{align} \mathcal{J} & = \begin{bmatrix} - \epsilon & p \beta_1 S_0 & p \beta_1 S_0 & 0 & p \beta_1 S_0 & p \beta_2 S_0 \\ \epsilon & -\sigma & 0 & 0 & 0 & 0 \\ 0 & \sigma & -\gamma_1 & 0 & 0 & 0 \\ 0 & (1-p) \beta_1 S_0 & (1-p) \beta_1 S_0 & - \epsilon & (1-p) \beta_1 S_0 & (1-p) \beta_2 S_0 \\ 0 & 0 & 0 & \epsilon & -\sigma & 0 \\ 0 & 0 & 0 & 0 & \sigma & -\gamma_2 - \alpha \end{bmatrix} \end{align}\]

If the real part of the dominant eigenvalue of \(\mathcal{J}\) is strictly greater than 0, the initial state (i.e. no epidemic) is unstable, implying that the epidemic can spread. However, the associated expressions are inelegant.

The Next Generation Theorem helps to overcome this (Diekmann et alii 1990, Hurford et alii 2010). To apply it, we divide the Jacobian matrix into a matrix of births \(\mathcal{F}\) and a matrix of deaths \(\mathcal{V}\), such as \(\mathcal{J}=\mathcal{F}-\mathcal{V}\) with \(\mathcal{F} \geq 0\), \(\mathcal{V}^{-1} \geq 0\); the modulus of the dominant eigenvalue of \(\mathcal{V}\) strictly positive.

Here we define: \[\begin{align} \mathcal{F} & = \begin{bmatrix} 0 & p \beta_1 S_0 & p \beta_1 S_0 & 0 & p \beta_1 S_0 & p \beta_2 S_0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & (1-p) \beta_1 S_0 & (1-p) \beta_1 S_0 & 0 & (1-p) \beta_1 S_0 & (1-p) \beta_2 S_0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \end{align}\]

and

\[\begin{align} \mathcal{V} & = \begin{bmatrix} \epsilon & 0 & 0 & 0 & 0 & 0 \\ - \epsilon & \sigma & 0 & 0 & 0 & 0 \\ 0 & - \sigma & \gamma_1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \epsilon & 0 & 0 \\ 0 & 0 & 0 & - \epsilon & \sigma & 0 \\ 0 & 0 & 0 & 0 & - \sigma & \gamma_2 + \alpha \end{bmatrix} \end{align}\]

If the modulus of the dominant eigenvalue of \(\mathcal{F}.\mathcal{V}^{-1}\) is greater than 1, then the initial state is unstable. This dominant eigenvalue corresponds to \(R_0\). The resolution (calculation of \(\mathcal{V}^{-1}\) and then of \(\mathcal{F}.\mathcal{V}^{-1}\) and then of its eigenvalues) leads to the expression given in the main text.

Other modeling sites

Sources and Acknowledgements

  • The ETE modeling team is composed of Samuel Alizon, Thomas Bénéteau, Marc Choisy, Gonché Danesh, Ramsès Djidjou-Demasse, Baptiste Elie, Yannis Michalakis, Bastien Reyné, Quentin Richard, Christian Selinger, Mircea T. Sofonea.

  • The calculation of \(R_0\) was facilitated by the use of the Next Generation theorem. For more details see:

    • Diekmann O, Heesterbeek JA, Metz JA. 1990. On the definition and the computation of the basic reproduction ratio \(R_0\) in models for infectious diseases in heterogeneous populations. J. Math. Biol. 28:365–82 https://doi.org/10.1007/BF00178324

    • Hurford A, Cownden D, Day T. 2010. Next-generation tools for evolutionary invasion analyses. J. R. Soc. Interface 7:561–71 https://doi.org/10.1098/rsif.2009.0448

  • Contribution to this work :

    • design of work: the whole team

    • carrying out analyses: BR, CS and SA

    • report writing: SA (translation YM)

    • contact :

    • approval: the whole team

  • Creative Commons License
    This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.