A multi-strain epidemic model for COVID-19 with infected and asymptomatic cases: Application to French data.

Many SARS-CoV-2 variants have appeared over the last months, and many more will continue to appear. Understanding the competition between these different variants could help make future predictions on the evolution of epidemics. In this study we use a mathematical model to investigate the impact of three different SARS-CoV-2 variants on the spread of COVID-19 across France, between January-May 2021 (before vaccination was extended to the full population). To this end, we use the data from Geodes (produced by Public Health France) and a particle swarm optimisation algorithm, to estimate the model parameters and further calculate a value for the basic reproduction number R0. Sensitivity and uncertainty analysis is then used to better understand the impact of estimated parameter values on the number of infections leading to both symptomatic and asymptomatic individuals. The results confirmed that, as expected, the alpha, beta and gamma variants are more transmissible than the original viral strain. In addition, the sensitivity results showed that the beta/gamma variants could have lead to a larger number of infections in France (of both symptomatic and asymptomatic people).

Introduction

Since the beginning of 2020, a pandemic of a novel coronavirus disease (COVID-19) has spread throughout the world and caused millions of deaths. COVID-19 is an infectious respiratory disease caused by the SARS-CoV-2 virus. The World Health Organization (WHO) learned of the existence of the SARS-CoV-2 on December 31, 2019 when an outbreak of “viral pneumonia” cases was notified in Wuhan, People’s Republic of China (Zhou et al., 2020, Tang et al., 2020). This disease has then spread to the whole world and very quickly led to saturated hospitals.

The fast spread of SARS-CoV-2 also led to the emergence of different variants (i.e., viral genomes that may contain one or more mutations) for this virus, who have been circulating around the world since the beginning of this pandemic. While most of the genetic mutations observed in the circulating SARS-CoV-2 variants do not significantly change virus biology and its properties, some fitness-enhancing mutations (e.g., mutations leading to an increase in transmissibility, or in ability to evade the immune response) have been observed after the first months of virus spread (Harvey et al., 2021). The circulating variants are classified as: Variants Being Monitored (VBM) – those for which data indicate a potential impact on approved treatments but might not pose yet a significant and imminent risk to public health, Variants of Interest (VOI) – for which there is predicted increase in transmissibility of disease severity, and reduced efficacy of treatments, Variants of Concern (VOC) – for which there is evidence of an increase in transmissibility, increase in the severity of disease, and reduced effectiveness of treatments (Centers for Disease Control and Prevention, 2021, European Centre for Disease Prevention and Control, 2021). Variants can be re-classified based in their attributes and prevalence. As of October 2021, the European Centre for Disease Prevention and Control (ECDC) has listed three VOC (Beta, Gamma and Delta) and two VOI (Mu and Lambda) (European Centre for Disease Prevention and Control, 2021).

To investigate the spread of COVID-19, including the role of the different SARS-CoV-2 variants on this spread, and to propose measures to slow-down this spread, researchers have focused their attention on various mathematical compartmental models, either stochastic (He et al., 2020, Faranda and Alberti, 2020) or deterministic (Calafiore and Fracastoro, 2021, He et al., 2020, Tang et al., 2020, Bentout et al., 2021, Yavuz et al., 2021, Arruda et al., 2021, Gonzalez-Parra et al., 2021, Khyar and Allali, 2020, Sridhar et al., 2021, Yagan et al., 2021); see also references therein. The great majority of these models focus on one single variant and investigate, for example, the impact of hospitalisation and quarantine (He et al., 2020, Tang et al., 2020), the impact of vaccination (Yavuz et al., 2021), or the impact of different age classes (Bentout et al., 2021, Calafiore and Fracastoro, 2021). Over the last few months, when it became clear the importance of different VOC on the fast increase in the number of COVID-19 cases in different countries, a series of mathematical models have been derived to investigate the role of multiple variants on the spread of SARS-CoV-2 (Arruda et al., 2021, Gonzalez-Parra et al., 2021, Khyar and Allali, 2020, Sridhar et al., 2021, Yagan et al., 2021). For example, Sridhar et al. (2021) and Yagan et al. (2021) used network models to investigate the effectiveness of mask-wearing in limiting the spread of COVID-19 in the context of viral mutations. Arruda et al. (2021) considered a SEIR model that incorporated several viral strains and also reinfection due to vaning immunity. They studied time-varying control strategies in the context of lockdown measures, and so they focused on the cost of infection control (i.e., lockdown) vs. the cost of elevated infection levels to the healthcare system, over a two-year time interval. Khyar and Allali (2020) focused on the global stability analysis of a two-strain SEIR epidemic model with two general incidence rates. They also calculated the basic reproduction number for their epidemic model. Finally, Gonzalez-Parra et al. (2021) extended a two-strain SEIR model to include also asymptomatic, hospitalized and dead individuals, to investigate the transmission of COVID-19 in Columbia. They consider variations in the contagiousness of the two strains, to see their impacts on the number of infections, hospitalisations and deaths.

In this study we also focus on the role of different SARS-CoV-2 variants on the transmission of this virus across France in the absence of vaccination. However, unlike the previous studies, here we consider a simpler compartmental model (i.e., a generalisation of a SIR model) that includes also asymptomatic cases and dead individuals infected with the original variant, as well as cases infected with other variants. To calculate the basic reproductive number , we first focus on a general model with variants, and derive a formula for . Then, we apply this formula for the case with two variants (), corresponding to the alpha and beta/gamma variants (see discussion on French data in Section 2 ). This model is parametrised (with the help of a particle swarm optimisation algorithm) using French data selected over the longest possible time interval that can take into account the emergence of different SARS-CoV-2 variants without including also the vaccination. Sensitivity analysis is performed to understand the impact of uncertainty in model parameters on the overall outcome (i.e., symptomatic and asymptomatic infections). This sensitivity analysis allows us to investigate the severity of infections with different variants in France, in the absence of any vaccination.

Fig. 2

The evolution of and , as given by data on Geodes (see also Table 4, Table 5, Table 6). This is the rescaled data that includes the assumption of 13 days in the infected compartments.

Table 4

Percentage of initial strain.

Date	Value	Date	Value	Date	Value
2021-02-12-2021-02-18	42,9	2021-02-13-2021-02-19	40,7	2021-02-14-2021-02-20	39,7
2021-02-15-2021-02-21	39,5	2021-02-16-2021-02-22	37,3	2021-02-17-2021-02-23	35,3
2021-02-18-2021-02-24	33,6	2021-02-19-2021-02-25	32,0	2021-02-20-2021-02-26	30,3
2021-02-21-2021-02-27	29,5	2021-02-22-2021-02-28	29,4	2021-02-23-2021-03-01	27,2
2021-02-24-2021-03-02	25,7	2021-02-25-2021-03-03	24,3	2021-02-26-2021-03-04	22,9
2021-02-27-2021-03-05	21,6	2021-02-28-2021-03-06	20,9	2021-03-01-2021-03-07	20,7
2021-03-02-2021-03-08	19,1	2021-03-03-2021-03-09	18,1	2021-03-04-2021-03-10	17,0
2021-03-05-2021-03-11	15,7	2021-03-06-2021-03-12	14,7	2021-03-07-2021-03-13	14,2
2021-03-08-2021-03-14	14,1	2021-03-09-2021-03-15	13,0	2021-03-10-2021-03-16	12,2
2021-03-11-2021-03-17	11,5	2021-03-12-2021-03-18	11,0	2021-03-13-2021-03-19	10,4
2021-03-14-2021-03-20	10,2	2021-03-15-2021-03-21	10,1	2021-03-16-2021-03-22	9,4
2021-03-17-2021-03-23	8,9	2021-03-18-2021-03-24	8,4	2021-03-19-2021-03-25	8,0
2021-03-20-2021-03-26	7,6	2021-03-21-2021-03-27	7,4	2021-03-22-2021-03-28	7,3
2021-03-23-2021-03-29	6,9	2021-03-24-2021-03-30	6,7	2021-03-25-2021-03-31	6,5
2021-03-26-2021-04-01	6,3	2021-03-27-2021-04-02	6,1	2021-03-28-2021-04-03	6,0
2021-03-29-2021-04-04	6,0	2021-03-30-2021-04-05	5,9	2021-03-31-2021-04-06	5,6
2021-04-01-2021-04-07	5,4	2021-04-02-2021-04-08	5,4	2021-04-03-2021-04-09	5,3
2021-04-04-2021-04-10	5,2	2021-04-05-2021-04-11	5,2	2021-04-06-2021-04-12	5,1
2021-04-07-2021-04-13	4,9	2021-04-08-2021-04-14	4,8	2021-04-09-2021-04-15	4,6
2021-04-10-2021-04-16	4,3	2021-04-11-2021-04-17	4,3	2021-04-12-2021-04-18	4,2
2021-04-13-2021-04-19	4,0	2021-04-14-2021-04-20	3,9	2021-04-15-2021-04-21	3,8
2021-04-16-2021-04-22	3,8	2021-04-17-2021-04-23	3,7	2021-04-18-2021-04-24	3,7
2021-04-19-2021-04-25	3,7	2021-04-20-2021-04-26	3,6	2021-04-21-2021-04-27	3,5
2021-04-22-2021-04-28	3,4	2021-04-23-2021-04-29	3,4	2021-04-24-2021-04-30	3,4
2021-04-25-2021-05-01	3,4	2021-04-26-2021-05-02	3,4	2021-04-27-2021-05-03	3,3
2021-04-28-2021-05-04	3,4	2021-04-29-2021-05-05	3,3	2021-04-30-2021-05-06	3,4
2021-05-01-2021-05-07	3,4	2021-05-02-2021-05-08	3,4	2021-05-03-2021-05-09	3,4
2021-05-04-2021-05-10	3,3	2021-05-05-2021-05-11	3,2	2021-05-06-2021-05-12	3,3
2021-05-07-2021-05-13	3,2

Table 5

Percentage of alpha variant.

Date	Value	Date	Value	Date	Value
2021-02-12-2021-02-18	44,3	2021-02-13-2021-02-19	46,0	2021-02-14-2021-02-20	46,9
2021-02-15-2021-02-21	47,1	2021-02-16-2021-02-22	49,2	2021-02-17-2021-02-23	51,0
2021-02-18-2021-02-24	52,8	2021-02-19-2021-02-25	54,3	2021-02-20-2021-02-26	55,9
2021-02-21-2021-02-27	56,6	2021-02-22-2021-02-28	56,8	2021-02-23-2021-03-01	59,0
2021-02-24-2021-03-02	60,6	2021-02-25-2021-03-03	62,1	2021-02-26-2021-03-04	63,5
2021-02-27-2021-03-05	64,8	2021-02-28-2021-03-06	65,6	2021-03-01-2021-03-07	65,8
2021-03-02-2021-03-08	67,2	2021-03-03-2021-03-09	68,1	2021-03-04-2021-03-10	69,1
2021-03-05-2021-03-11	70,1	2021-03-06-2021-03-12	71,2	2021-03-07-2021-03-13	71,6
2021-03-08-2021-03-14	71,7	2021-03-09-2021-03-15	73,0	2021-03-10-2021-03-16	73,9
2021-03-11-2021-03-17	74,6	2021-03-12-2021-03-18	75,2	2021-03-13-2021-03-19	75,6
2021-03-14-2021-03-20	75,9	2021-03-15-2021-03-21	75,9	2021-03-16-2021-03-22	76,5
2021-03-17-2021-03-23	76,9	2021-03-18-2021-03-24	77,7	2021-03-19-2021-03-25	78,1
2021-03-20-2021-03-26	78,8	2021-03-21-2021-03-27	79,2	2021-03-22-2021-03-28	79,3
2021-03-23-2021-03-29	80,1	2021-03-24-2021-03-30	80,6	2021-03-25-2021-03-31	80,8
2021-03-26-2021-04-01	81,2	2021-03-27-2021-04-02	81,4	2021-03-28-2021-04-03	81,5
2021-03-29-2021-04-04	81,5	2021-03-30-2021-04-05	81,4	2021-03-31-2021-04-06	81,8
2021-04-01-2021-04-07	82,1	2021-04-02-2021-04-08	82,2	2021-04-03-2021-04-09	82,3
2021-04-04-2021-04-10	82,6	2021-04-05-2021-04-11	82,6	2021-04-06-2021-04-12	82,7
2021-04-07-2021-04-13	82,5	2021-04-08-2021-04-14	82,3	2021-04-09-2021-04-15	82,4
2021-04-10-2021-04-16	82,4	2021-04-11-2021-04-17	82,3	2021-04-12-2021-04-18	82,3
2021-04-13-2021-04-19	82,2	2021-04-14-2021-04-20	82,5	2021-04-15-2021-04-21	82,6
2021-04-16-2021-04-22	82,5	2021-04-17-2021-04-23	82,4	2021-04-18-2021-04-24	82,3
2021-04-19-2021-04-25	82,3	2021-04-20-2021-04-26	82,2	2021-04-21-2021-04-27	81,9
2021-04-22-2021-04-28	81,6	2021-04-23-2021-04-29	81,3	2021-04-24-2021-04-30	80,8
2021-04-25-2021-05-01	80,5	2021-04-26-2021-05-02	80,5	2021-04-27-2021-05-03	80,2
2021-04-28-2021-05-04	79,9	2021-04-29-2021-05-05	79,7	2021-04-30-2021-05-06	79,4
2021-05-01-2021-05-07	79,1	2021-05-02-2021-05-08	79,1	2021-05-03-2021-05-09	79,1
2021-05-04-2021-05-10	78,2	2021-05-05-2021-05-11	77,9	2021-05-06-2021-05-12	77,6
2021-05-07-2021-05-13	77,4

Table 6

Percentage of beta variant.

Date	Value	Date	Value	Date	Value
2021-02-12-2021-02-18	4,5	2021-02-13-2021-02-19	4,6	2021-02-14-2021-02-20	4,6
2021-02-15-2021-02-21	4,6	2021-02-16-2021-02-22	4,6	2021-02-17-2021-02-23	4,8
2021-02-18-2021-02-24	4,9	2021-02-19-2021-02-25	5,0	2021-02-20-2021-02-26	5,2
2021-02-21-2021-02-27	5,4	2021-02-22-2021-02-28	5,4	2021-02-23-2021-03-01	5,5
2021-02-24-2021-03-02	5,4	2021-02-25-2021-03-03	5,3	2021-02-26-2021-03-04	5,2
2021-02-27-2021-03-05	5,0	2021-02-28-2021-03-06	4,9	2021-03-01-2021-03-07	4,8
2021-03-02-2021-03-08	4,7	2021-03-03-2021-03-09	4,7	2021-03-04-2021-03-10	4,9
2021-03-05-2021-03-11	5,0	2021-03-06-2021-03-12	5,0	2021-03-07-2021-03-13	5,0
2021-03-08-2021-03-14	5,0	2021-03-09-2021-03-15	5,0	2021-03-10-2021-03-16	5,0
2021-03-11-2021-03-17	4,9	2021-03-12-2021-03-18	4,8	2021-03-13-2021-03-19	4,7
2021-03-14-2021-03-20	4,6	2021-03-15-2021-03-21	4,6	2021-03-16-2021-03-22	4,6
2021-03-17-2021-03-23	4,5	2021-03-18-2021-03-24	4,4	2021-03-19-2021-03-25	4,4
2021-03-20-2021-03-26	4,4	2021-03-21-2021-03-27	4,3	2021-03-22-2021-03-28	4,3
2021-03-23-2021-03-29	4,3	2021-03-24-2021-03-30	4,2	2021-03-25-2021-03-31	4,2
2021-03-26-2021-04-01	4,1	2021-03-27-2021-04-02	4,1	2021-03-28-2021-04-03	4,0
2021-03-29-2021-04-04	4,1	2021-03-30-2021-04-05	4,1	2021-03-31-2021-04-06	4,0
2021-04-01-2021-04-07	3,9	2021-04-02-2021-04-08	3,9	2021-04-03-2021-04-09	3,9
2021-04-04-2021-04-10	3,8	2021-04-05-2021-04-11	3,8	2021-04-06-2021-04-12	3,8
2021-04-07-2021-04-13	3,9	2021-04-08-2021-04-14	4,0	2021-04-09-2021-04-15	4,0
2021-04-10-2021-04-16	4,1	2021-04-11-2021-04-17	4,2	2021-04-12-2021-04-18	4,2
2021-04-13-2021-04-19	4,4	2021-04-14-2021-04-20	4,6	2021-04-15-2021-04-21	4,7
2021-04-16-2021-04-22	4,9	2021-04-17-2021-04-23	5,0	2021-04-18-2021-04-24	5,1
2021-04-19-2021-04-25	5,1	2021-04-20-2021-04-26	5,1	2021-04-21-2021-04-27	5,1
2021-04-22-2021-04-28	5,2	2021-04-23-2021-04-29	5,3	2021-04-24-2021-04-30	5,4
2021-04-25-2021-05-01	5,5	2021-04-26-2021-05-02	5,5	2021-04-27-2021-05-03	5,6
2021-04-28-2021-05-04	5,7	2021-04-29-2021-05-05	5,7	2021-04-30-2021-05-06	5,7
2021-05-01-2021-05-07	5,7	2021-05-02-2021-05-08	5,7	2021-05-03-2021-05-09	5,7
2021-05-04-2021-05-10	5,9	2021-05-05-2021-05-11	6,0	2021-05-06-2021-05-12	5,9
2021-05-07-2021-05-13	6,0

The paper is structured as follows. In Section 2 we introduce a general compartment model that considers a generic number of N variants, that can lead to symptomatic as well as asymptomatic cases. We use this model to calculate a general formula for in terms of the parameters associated with different variants. In Section 3 we discuss SARS-CoV-2 data for France, and based on this data we focus on a simpler model with variants: alpha and beta. Using available data, we parametrise this simpler model with the help of a particle swarm optimisation algorithm. In Section 4 we estimate the basic reproductive number for the three SARS-CoV-2 variants considered here, and we perform a sensitivity analysis to investigate the impact of uncertainties in the previously-identified model parameters on the model outcome (i.e., the number of infectious individuals, as well as ).

A general model with N variants

We start by presenting a general compartmental model that takes into account a generic number N of variants, with . We choose N variants because it allows us to include all current and possible future VOC variants. Symptomatic and asymptomatic cases are taken into consideration with the assumption that asymptomatic people do not die from the disease. Natural births and deaths can be overlooked, due to the large size of the French population. Loss of immunity can also be ignored because the number of individuals infected twice is still low (Brouqui et al., 2021), and this loss of immunity is not present for all variants. For example, variants with the E484K mutation in the spike protein may be responsible for immune escape (Nonaka et al., 2021, Wise, 2021).

Model description

For the model considered in this study, we denote the initial strain by , the first variant by , …, the N-th variant by . The population is divided into different classes: susceptible individuals (S) who are healthy and can contract the disease; individuals infected with the original virus and having symptoms (); individuals infected with the original virus and having no symptoms (); individuals infected with the variants , …, and having symptoms (, …, ); individuals infected with the variants , …, and having no symptoms (, …, ); dead individuals (D); recovered individuals (R) who are immune to the disease.

The dynamics of the population is described by the system below (see also Fig. 1 for a schematic description of the interactions between different model variables).

Fig. 1

Diagram describing the interactions incorporated into model (1).

In the above system (1), the equation for describes the evolution of individuals infected with the initial strain and having symptoms (). Susceptible persons S come into contact with persons infected by the initial strain with symptoms () or without symptoms (). There is a certain probability of being infected by the initial virus which results in the terms and . A proportion of these individuals will have symptoms. Finally, infected individuals can die at a rate or can recover at a rate . Similar terms can be found in the equations describing the evolution of the other symptomatic infected compartments ,…,; see also Table 1 .

Table 1

Summary of parameters involved in model (1). Parameters and are respectively equal to and , with N the size of the population considered (hence the units for these parameters are ). Here, and are the average duration times between two contaminating contacts (with the strain k) by a symptomatic person and by an asymptomatic person.

Symbol	Biological interpretation	Unit
μ^V_0,μ^V_1,..,μ^V_N	Mortality rate due to the initial virus and to other N variant	day^-1
1δ^V_0,1δ^V_1,…,1δ^V_N	Average time of infection for symptomatic individuals	day
1b^V_0,1b^V_1,…,1b^V_N	Average time of infection for asymptomatic individuals	day
αiV_0,αiV_1,…,αiV_N	Transmission rate of variants and initial virus from symptomatic people	size^-1×day^-1
αaV_0,αaV_1,…,αaV_N	Transmission rates of variants and initial virus from asymptomatic people	size^-1×day^-1
p^V_0,p^V_1,…,p^V_N	Probability that an individual infected with the initial virus, variant 1, .. variant N, is symptomatic	–

The equation for in (1) describes the evolution of the individuals infected by the initial strain and having no symptoms (). A proportion of susceptibles can become infected and have no symptoms, and as before there is a probability of infection following contact with other symptomatic individuals or asymptomatic individuals . In addition, individuals can leave this asymptomatic class following recovery at a rate . We remind the reader that the asymptomatic persons cannot die. Similar terms can be found in the equations describing the evolution of the other compartments ,…,; see also Table 1.

The equation for in (1) describes the evolution of recovered individuals. The terms that appear in this equation describe the recovery of infected and asymptomatic individuals, as discussed above.

The equation for in (1) describes the time-evolution of COVID-19-dead individuals. The terms in this equation are the death terms that we discussed above for the symptomatic infected compartments.

The equation for in (1) describes the time-evolution of the susceptible people. These individuals leave the compartment S when becoming infected (after contact with symptomatic or asymptomatic people).

All parameters involved in this general model (1) with N viral variants are summarized in the Table 1. In the following we discuss briefly the calculation of the basic reproduction number for this general model.

calculation

A fundamental tool in epidemiology is the basic reproduction number , which represents the number of secondary infections resulting from a single infectious individual introduced into a fully susceptible population (Diekmann et al., 1990, Van den Driessche and Watmough, 2002, Perasso, 2018). Note that in general a population is not fully susceptible, with some individuals already having immunity/cross-immunity to the virus. Therefore, even if throughout this study we refer to the basic reproductive number , we actually understand the effective reproduction number (i.e., the number of secondary infections resulting from a single infectious individual introduced into a population formed of susceptible and non-susceptible individuals); this will be more clear in Section 3 in the context of alpha and beta variants for COVID-19.

One method that can be used to calculate the basic reproduction number for the finite-dimensional system (1) makes use of the next generation matrix (Van den Driessche and Watmough, 2002). This method is based on the definition of as the dominant eigenvalue of the “next generation matrix” (Diekmann et al., 1990), i.e., a matrix that relates the numbers of newly infected individuals in various categories in consecutive generations (Diekmann et al., 2010). Applying this method to model (1) we obtain the next generation matrix M below:

Thus, we obtain the following formula for :

Each eigenvalue consists only of the parameters corresponding to a single strain. It is like there are several epidemics at the same time, and they don’t interfere with each other. By obtaining a value of from the available data (see Section 3), we will be able to conclude which variant takes over the others.

Parameters estimation for a 2 variant-model (): alpha and beta/gamma variants

Throughout the rest of this paper we focus on a simpler model with only two variants, (i.e., ), and we use the available data to obtain an estimate for (or rather the effective reproduction number, as discussed above). We estimate the 18 parameters involved in model (1) and in Eq. (2) using a Particle Swarm Optimization algorithm.

Data

The data we use in this study is from Geodes (i.e., the cartographic observatory of epidemiological indicators produced by Public Health France1 ), and covers the period between 12th February and 7th May 2021, because before the 12 February the different strains of the virus were not recorded (only virus presence was recorded). Moreover, we consider data only until 7th May to avoid the effects of vaccines (since in France the vaccination really started off in May). We focus on the percentages of positive RT-PCR tests identifying the initial strain, the alpha variant (which appeared in France since mid-December 2020) and the beta & gamma variants (present in France since end of December 2020 & beginning of February 2021, respectively), calculated over a 7-day period. This period allows to remove the effect of reduced testing on Sundays. The data can be seen in Table 4, Table 5, Table 6, in Appendix A. Note that on Geodes, the data for the beta and gamma variants are presented together, and for this reason in this study we combine them into a single variant (called variant 2 in the next section).

On Geodes, the data included also percentages of a so-called “indeterminate variant”. To remove this indeterminate variant (so we can have only the initial strain, the alpha variant and the beta/gamma variant), we considered the following approach. Denote the alpha strain by , the beta strain by and the initial strain by . Then we calculate:with

Similar calculations were performed to obtain the new percentages for the beta/gamma variants, as well as the initial strain.

The number of new confirmed cases per day in France is available on the official french public data site 2 . These data are summarised in Fig. 7, in Appendix. To remove the effects of Sunday testing, when fewer RT-PCR tests are carried out, we take an average over 7 days at each time point, to correspond to the values of the percentages of the variants. These data are shown in Fig. 8 in Appendix.

Table 7

Number of new confirmed cases in France per day.

Date	Value	Value	Value	Value	Value
2021-02-12-2021-02-16	20701	21231	16546	4376	19590
2021-02-17-2021-02-21	25018	22501	24116	22371	22046
2021-02-22-2021-02-26	4646	20064	31519	25403	25207
2021-02-27-2021-03-03	23996	19952	4703	22857	26788
2021-03-04-2021-03-08	25279	23507	23306	21825	5327
2021-03-09-2021-03-13	23302	30303	27166	25229	29759
2021-03-14-2021-03-18	26343	6471	29975	38501	34998
2021-03-19-2021-03-23	35088	35327	30581	15792	14678
2021-03-24-2021-03-28	65373	45641	41869	42619	37014
2021-03-29-2021-04-02	9094	30702	59038	50659	46677
2021-04-03-2021-04-07	13917	66794	10793	8045	12951
2021-04-08-2021-04-12	84999	41243	43284	34895	8536
2021-04-13-2021-04-17	39113	43505	38045	36442	35861
2021-04-18-2021-04-22	29344	6696	43098	34968	34318
2021-04-23-2021-04-27	32340	32633	24465	5952	30317
2021-04-28-2021-05-02	31539	26538	24299	25670	9888
2021-05-03-2021-05-07	3760	24371	26000	21712	19124
2021-05-08-2021-05-12	20745	9128	3292	19791	21498
2021-05-13	19461

Table 8

Number of new confirmed cases in France per day, averaged over 7 days.

Date	Value	Date	Value	Date	Value
12/02-18/02	18566,14286	13/02-19/02	19054,00000	14/02-20/02	19216,85714
15/02-21/02	20002,57143	16/02-22/02	20041,14286	17/02-23/02	20108,85714
18/02-24/02	21037,57143	19/02-25/02	21452,14286	20/02-26/02	21608,00000
21/02-27/02	21840,14286	22/02-28/02	21541,00000	23/02-01/03	21549,14286
24/02-02/03	21948,14286	25/02-03/03	21272,28571	26/02-04/03	21254,57143
27/02-05/03	21011,71429	28/02-06/03	20913,14286	01/03-07/03	21180,71429
02/03-08/03	21269,85714	03/03-09/03	21333,42857	04/03-10/03	21835,57143
05/03-11/03	22105,14286	06/03-12/03	22351,14286	07/03-13/03	23273,00000
08/03-14/03	23918,42857	09/03-15/03	24081,85714	10/03-16/03	25035,14286
11/03-17/03	26206,28571	12/03-18/03	27325,14286	13/03-19/03	28733,57143
14/03-20/03	29529,00000	15/03-21/03	30134,42857	16/03-22/03	31466,00000
17/03-23/03	29280,71429	18/03-24/03	33119,57143	19/03-25/03	34640,00000
20/03-26/03	35608,71429	21/03-27/03	36650,42857	22/03-28/03	37569,42857
23/03-29/03	36612,57143	24/03-30/03	38901,71429	25/03-31/03	37996,71429
26/03-01/04	38713,57143	27/03-02/04	39400,42857	28/03-03/04	35300,14286
29/03-04/04	39554,42857	30/03-05/04	39797,14286	31/03-06/04	36560,42857
01/04-07/04	29976,57143	02/04-08/04	34882,28571	03/04-09/04	34106,00000
04/04-10/04	38301,28571	05/04-11/04	33744,28571	06/04-12/04	33421,85714
07/04-13/04	37860,14286	08/04-14/04	42225,00000	09/04-15/04	35517,28571
10/04-16/04	34831,42857	11/04-17/04	33771,00000	12/04-18/04	32978,00000
13/04-19/04	32715,14286	14/04-20/04	33284,42857	15/04-21/04	32064,85714
16/04-22/04	31532,42857	17/04-23/04	30946,42857	18/04-24/04	30485,28571
19/04-25/04	29788,28571	20/04-26/04	29682,00000	21/04-27/04	27856,14286
22/04-28/04	27366,28571	23/04-29/04	26254,85714	24/04-30/04	25106,14286
25/04-01/05	24111,42857	26/04-02/05	22029,00000	27/04-03/05	21715,85714
28/04-04/05	20866,42857	29/04-05/05	20075,14286	30/04-06/05	19385,71429
01/05-07/05	18646,42857	02/05-08/05	17942,85714	03/05-09/05	17834,28571
04/05-10/05	17767,42857	05/05-11/05	17113,14286	06/05-12/05	16470,00000
07/05-13/05	16148,42857

From the percentages of each variant as well as the number of people infected at each time point t (i.e., day), we obtain three datasets corresponding to the number of new cases per day for each of the three strains. These three datasets are of size 85 corresponding to the period from February 12 to May 7. However, the different symptomatic infected compartments (I) in system (1) (with ) count the current numbers of infected individuals and not the newly infected per day. To transform the data so that we can use it to parametrise our model, we assume that individuals stay an average of 13 days in a symptomatic compartment before recovering (George et al., 2021). By including this recovery delay, we reduced the size of the three datasets from 85 to 73 days:with denoting the rescaled cumulative data that consider the assumption of 13 days infections (i.e., we sum up all daily infections that took place over the past 13 days), which is then used to parametrise the mathematical model.

The evolution of the three infected compartments can be seen in Fig. 2. Note that the curve for the infected population with the alpha variant (panel (b)) starts decreasing around day (marked by the vertical dashed line). This corresponds to 2 April 2021, when the French government imposed a 4-weeks national lockdown (which lasted until 3 May). At the same time, the beta/gamma variant starts to increase slightly until , after which it decreases as well.

Data for initial population sizes. For the numerical simulations performed throughout the rest of the paper, we consider the following initial population sizes:

Using the latest census data from mainland France (Décembre, 2020), in this study we assume that the total population of France is . We consider this initial condition because on 12 February 2021 (when we start our investigation into the dynamics of these different variants) there were 3,267,386 people in France who recovered from the disease (i.e., ) and 81,448 people who died (i.e., ) (Santé publique france, 2021).

For the infected sub-populations, we use the datasets in Table 4, Table 5, Table 6: (the first value of the first dataset), (the first value of the second dataset), (the first value of the third dataset). Regarding the infected people, we assume that are symptomatic and are asymptomatic (Xiao et al., 2020), and that the asymptomatic are not tested the same way as the symptomatic infected individuals and thus we take , .

For the Susceptible population, we need to subtract from the total population the number of individuals infected with the different variants:

Particle Swarm Optimization

The Particle Swarm Optimization (PSO) algorithm was created in 1995 by Russel Eberhart, electrical engineer, and James Kennedy, social psychologist (Kennedy and Eberhart, 1995), and it is inspired by the collective behaviour of flocks of birds or schools of fish. For the PSO, the particle (i.e., potential solution of the model) moves by making allowances between getting closer to the optimal solution visited and getting closer to the solutions found in its neighborhood. This optimization algorithm performs well on parameter optimization for ordinary differential equation models (Akman et al., 2018), but convergence to the overall optimal solution is not always guaranteed. In the following we describe how we apply this PSO algorithm to identify the parameters of our model (1) with two SARS-CoV-2 variants, using the data described above.

Objective function and search bounds for PSO

For each viral strain we calculate the residual sum of squareswhere

.

describe the empirical observations at time corresponding to the number of individuals infected by strain .

describe the numerical predictions at time for the individuals infected by strain (as given by system (1) with ).

The objective function to be minimized is therefore

Based on reasonable parameter bounds (see discussion below), we constrain our search space with lower and upper bounds for :

This means that infected people take at most 25 days and at least 1 day to infect a person (Amiri Mehra et al., 2020). Indeed, the transmission rate has been estimated at for the whole world (Marinov et al., 2020) which corresponds to days to transmit the virus (because N = 7,870,000,000 is the world population). The symptomatic proportion is between 0.55 and 0.9 (Xiao et al., 2020), a symptomatically-infected person takes between 3 days (Pottier, 2020) and 13 days (George et al., 2021) to recover, and between 5 days (Fernández-Villaverde and Jones, 2020) and 13 days to die (also to match our infection hypothesis (George et al., 2021)).

Parameter estimation

In Fig. 3 we show the best fit between the data collected from Geodes (Geodes, 2021) on infected individuals (black dots), and the numerically-predicted number (red curve) of individuals infected with (a) initial strain, (b) alpha strain, (c) beta/gamma strain. The parameters for which the numerical solutions were obtained are summarised in Table 2 , while Table 3 gives the biological interpretation of these parameters. The parameters and are respectively equal to and , with N the size of the population considered, and the average duration between two contaminating contacts by a symptomatic person and by an asymptomatic person by the strain k. For the biological interpretation of the other parameters see Table 1.

Fig. 3

Fitting the numerically-predicted infection curves (red curves) for the three strains, to the data collected from Geodes (black dots). The fit was performed using the Particle Swarm Optimisation (PSO) algorithm.

Table 2

Parameter values estimated from Geodes data using the PSO algorithm.

Symbol	Value	Symbol	Value	Symbol	Value
αiV_0	2.198814e−09	αiV_1	6.748454e−09	αiV_2	1.970476e−09
αaV_0	2.580832e−09	αaV_1	1.099686e−08	αaV_2	1.353279e−08
p^V_0	7.594355e−01	p^V_1	8.531318e−01	p^V_2	8.229163e−01
b^V_0	1.490951e−01	b^V_1	2.949318e−01	b^V_2	2.110977e−01
δ^V_0	7.926397e−02	δ^V_1	3.122869e−01	δ^V_2	1.194538e−01
μ^V_0	7.727916e−02	μ^V_1	1.184327e−01	μ^V_2	9.728384e−02

Table 3

Biological interpretation of parameters estimated here.

Viral strain	Parameter interpretation	Value
V_0	Number of days for a contaminating contact with a symptomatic	7.11
	Number of days for a contaminating contact with a asymptomatic	6.05
	Percentage of symptomatic individuals	0.76
	Asymptomatics: number of days to recover	6.71
	Symptomatics: number of days to recover	12.62
	Number of days before dying	12.94
V_1 (alpha)	Number of days for a contaminating contact with a symptomatic	2.31
	Number of days for a contaminating contact with a asymptomatic	1.42
	Percentage of symptomatic individuals	0.85
	Asymptomatics: number of days to recover	3.39
	Symptomatics: number of days to recover	3.20
	Number of days before dying	8.44
V_2 (beta/ gamma)	Number of days for a contaminating contact with a symptomatic	7.93
	Number of days for a contaminating contact with a asymptomatic	1.15
	Percentage of symptomatic individuals	0.82
	Asymptomatics: number of days to recover	4.73
	Symptomatics: number of days to recover	8.37
	Number of days before dying	10.28

The fit of the individuals infected with the initial and alpha strain is relatively good. In contrast, the fit of individuals infected with the beta strain is quite poor for time days and days.

The evolution of asymptomatic infected people is presented in the Fig. 4 . We see that the number of asymptomatic individuals infected with the initial strain increases at the beginning and then decreases (following the trend of symptomatic infected individuals). The asymptomatic individuals infected with the other two variants (alpha and beta) follow the same general pattern as the symptomatic individuals but with fewer cases.

Fig. 4

Time-evolution of the three asymptomatic compartments of model (1).

Looking at the numbers in Table 3, Table 2 we see that the people who die the fastest are those affected by the alpha variant. Indeed, . (Note that this result is also supported by the values of the case fatality rate (CFR) for each of the three variants: .) In addition, people who take the least time to transmit the virus are those affected by the alpha variant which is consistent with the fact that it is said to be more transmissible. There are a higher number of days for the transmission of the disease to symptomatic people than for the asymptomatic people, and this is true for all three cases.

Results

Estimation of the basic reproduction number

To calculate the value of , we have computed three non-zero eigenvalues for the next generation matrix: (corresponding to the initial strain), (corresponding to the alpha strain), (corresponding to the beta/gamma strain). It is clear that the two variants (with ) are more contaminating than the initial strain (with ). Both alpha and beta variants have very close values, so we can conclude at first that they have a similar intensity, at least for France, in the absence of vaccination. Note that this value is consistent with the value calculated by other French studies that investigated COVID-19 transmission during the same period of time (Gaymard et al., 2020).

Local sensitivity analysis

We now investigate the impact on when we vary by one day the biological parameter estimates given in Table 3. For example, if we estimated at the number of days required for the transmission of the disease by a symptomatic person, we now consider the interval days. Each time we vary one parameter, while fixing all other parameters. The results are presented in Fig. 5 .

Fig. 5

Local sensitivity of when we vary: (a) parameters associated with the original strain ; (b) parameters associated with the alpha strain ; (c) parameters associated with the beta/gamma strain . The right-hand sub-panel in (c) shows, on a much larger vertical axis, as we vary .

First, we see in Fig. 5(a) that the parameters of the initial strain do not impact because the corresponding value is never the maximum among the three eigenvalues of the next generation matrix. For the alpha variant (Fig. 5(b)), we see that is especially sensitive to and , parameters describing the transmission of infection following contacts with symptomatic and asymptomatic infected individuals. Thus, reducing the contacts (e.g., via lockdown) could lead to a reduction in . Also notice that none of the parameters related to the alpha variant increase above 2. For the beta/gamma variants (Fig. 5(c)), we see that is very sensitive to . Changes in this parameter can increase up to (unlike the case for the alpha variant, where was always below 2). This shows us the importance of controlling the parameter to be able to control the outbreak of the beta/gamma variant. Another parameter that impacts is , the recovery rate of asymptomatic individuals. This shows the key role of the asymptomatic people in the rapid spread of the beta/gamma strain.

Global sensitivity analysis for the infected compartments

Since knowledge about model parameters is incomplete (especially knowledge about the role of asymptomatic people on disease transmission), we are now performing a global sensitivity and uncertainty analysis of the different infected compartments, where we vary all parameters at the same time. To this end, we consider a classical approach that combines the Latin Hypercube Sampling (LHS) with the Partial Rank Correlation Coefficient (PRCC) (Marino et al., 2008). The parameters are varied in the same way as before (i.e., by one day for each of the rates), and the sample size for the LHS is 100.

In Fig. 6 we show the variation in the number of individuals infected with the original strain, as we vary all model parameters. The black curves show the median time evolution of and , while the gray regions show the range between quantiles and . We can see that both and are very sensitive for days. For example, a 1-day change in the various transition rates incorporated into the model, leads on day (i.e. February 22) to a variation between 50,000 to 80,000 in the symptomatic infected people (), and a variation between 10,000 to 25,000 asymptomatic infected people approximately. In Fig. 7 we show the variation in the number of individuals infected with the alpha strain, as we vary all model parameters. Note here the change in the maximum peak of both and : not only the amplitude but also the shift in time: from day for the median black curve, to day for the quantiles. If we compare these results with the results in Fig. 8 for the beta/gamma strain, we observe that while the maximum for is reached around days, for the maximum is reached for days. Therefore, the beta/gamma strains are more aggressive. Moreover, by comparing Fig. 7, Fig. 8, we see that and could reach much higher amplitudes (i.e., 10-fold higher) compared to and . We will return to this aspect in the Discussion Section 5.

Fig. 6

Uncertainty in the and , as we vary by a maximum of day all model parameters (according to the LHS scheme).

Fig. 7

Uncertainty in the and as we vary by a maximum of day all model parameters (according to the LHS scheme).

Fig. 8

Uncertainty in the and as we vary by a maximum of day all model parameters (according to the LHS scheme).

Finally, for completeness, in Fig. 9 we also show the time-variation in the recovered and dead people. Again, we see a large variability in the outcome between days and which, for lasts also for days.

Fig. 9

Uncertainty in the and as all model parameters are varied by one day.

To quantify the sensitivity of model outcomes (i.e., number of symptomatic and asymptomatic infected individuals) to changes in model parameters, and identify the critical parameters, in Fig. 10, Fig. 11, Fig. 12 we show the Partial Rank Correlation Coefficient (PRCC) – a sampling-based measure for nonlinear but monotonic relationships between inputs and outputs (Marino et al., 2008). We see here that the parameters with the greatest impact on the outcome (i.e., parameters with PRCC indexes above and below ) are:

Fig. 10

PRCC of and .

Fig. 11

PRCC of and .

Fig. 12

PRCC of and .

For : parameters and ; while for : parameters , and ;

For : parameters and ; while for : parameters and ;

For : parameter and ; while for : parameters and .

First we observe that due to the interactions between the model variables, parameters associated with the second variant () impact the infections with the original () strain and the infections with the first variant (). For example, impacts not only and , but also and , which suggest that the people infected by beta/gamma strain could have a strong impact on the overall evolution of epidemics. This is despite the fact that for , the impact of beta/gamma variant was decoupled from the impact of the original variant and alpha variant (see Eqn. (2)).

Discussion & conclusion

In this paper we investigated the evolution of COVID-19 epidemic in France between February-May 2021 in the presence of multiple viral strains, while estimating the parameters that characterise this epidemic. The time interval for the data was chosen to reduce the impact of vaccination on the results (although other control measures such as masks and social distancing were still in place in France during this period). We focused on three viral variants identified in the data: the original (i.e., March 2020) SARS-CoV-2 variant, the alpha variant (first detected in the UK) and the beta/gamma variant (first detected in South Africa and in Brazil). The last two were considered together here, because in France the tests did not discriminate between them.

First, we showed that the value for was given by the maximum of three basic reproduction values corresponding to the three different variants. This suggested a sort of decoupling between the evolution of different variants, the most transmissible one having the greatest influence on the evolution of the outbreak. Using French data on Geodes (2021) and a PSO algorithm, we estimated , consistent with the value estimated by other French studies (Gaymard et al., 2020). Moreover, this value was similar for the alpha and beta/gamma variants, but larger than the value corresponding to the original strain (thus suggesting that between February-May the alpha and beta/gamma variants were more transmissible than the original strain). A local sensitivity analysis for showed that the parameters with the largest impact were the transmission rates of symptomatic and asymptomatic individuals; see Fig. 5.

In regard to parameter estimation (see Table 3) we also found a higher number of days necessary to create contaminating contacts for the symptomatic people than for the asymptomatic people (for all three variants). This unexpected result could be explained by the fact that asymptomatic people come into contact with more peers ignoring distancing rules. The results in Table 3 show also that the number of days to recover when symptomatic is slightly lower than for asymptomatic people for the alpha variant, while it is almost twice more for the other strains. It would be interesting to investigate whether there is a specificity of the alpha variant in contrast to the other variants. At this point we could not find any epidemiological data to test the validity of this theoretical result.

Second, a global sensitivity and uncertainty analysis identified the variations in the amplitude of symptomatic and asymptomatic infections for all three strains, as well as the day when these infections peak (see Fig. 6, Fig. 7, Fig. 8). Moreover, the PRCC analysis identified some critical parameters for the evolution of the epidemics (see Fig. 10, Fig. 11, Fig. 12). The interesting result was that parameters associated with the variant () seemed to impact also the evolution of infections with the and variants, thus suggesting that the infections with the different viral strains are not really decoupled (as suggested by the formula (2)).

The uncertainty analysis also showed that the evolution of those infected with the beta/gamma strain i.e. has more uncertainty compared to the other two strains (when we vary the interaction rates by one day).

Also in regard to the uncertainty of the results in this study, we need to mention the fact that since we could not find mortality data for each of the different strains discussed here, the exact identification of parameters actually depends on all other model parameters listed in Table 2.

Throughout this study we chose the average infection period to be 13 days, since this was similar to other published studies (George et al., 2021). Even reducing this period by 1–2 days did not seem to lead to significant changes. However, in the future it would be interesting to investigate model dynamics when we change (increase/decrease) this infection period in a more significant way.

This study focused on the dynamics of COVID-19 epidemics in the presence of alpha and beta/gamma strains, with distancing measures (and masks worn indoor and outdoor) but no vaccinations (as the vaccination program in France really took off in May 2021). Our results about the potential large amplitude in the infections with beta and gamma strains, can be understood in this context: these strains would have caused a disastrous outbreak situation if the situation had not changed with the arrival of the vaccine. The subsequent emergence of the delta variant also changed the dynamics of the epidemics, and if we include it into our model (1) the simulation results and predictions will inevitably change.

The study can be further extended to the new SARS-CoV-2 variants that are appearing, which might avoid the anti-viral immune responses generated by the vaccination. In regard to this, a very recent review (Malik et al., 2022) discussed vaccine effectiveness in the context of the multiple SARS-CoV-2 variants that emerged over the last months, and concluded that since a large number of variants have mutations mainly associated with the spike protein, which is also a key component of most of the vaccines on the market, vaccine efficacy needs to be assessed for each variant. While some vaccine efficacy studies have been published for earlier variants (Bian et al., 2021), many more such studies are ongoing.

CRediT authorship contribution statement

Mathilde Massard: Conceptualization, Writing – original draft, Writing – review & editing. Raluca Eftimie: Supervision. Antoine Perasso: Supervision. Bruno Saussereau: Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.