Mathematical modeling of the dynamics of COVID-19 variants of concern: Asymptotic and finite-time perspectives.

The COVID-19 pandemic has seen multiple waves, in part due to the implementation and relaxation of social distancing measures by the public health authorities around the world, and also caused by the emergence of new variants of concern (VOCs) of the SARS-Cov-2 virus. As the COVID-19 pandemic is expected to transition into an endemic state, how to manage outbreaks caused by newly emerging VOCs has become one of the primary public health issues. Using mathematical modeling tools, we investigated the dynamics of VOCs, both in a general theoretical framework and based on observations from public health data of past COVID-19 waves, with the objective of understanding key factors that determine the dominance and coexistence of VOCs. Our results show that the transmissibility advantage of a new VOC is a main factor for it to become dominant. Additionally, our modeling study indicates that the initial number of people infected with the new VOC plays an important role in determining the size of the epidemic. Our results also support the evidence that public health measures targeting the newly emerging VOC taken in the early phase of its spread can limit the size of the epidemic caused by the new VOC (Wu et al., 2139Wu, Scarabel, Majeed, Bragazzi, & Orbinski, ; Wu et al., 2021).

Introduction

The COVID-19 pandemic, caused by the infection of SARS-Cov-II virus, has become one of the most severe and deadly pandemics in recent history. By May 2022, more than two years after its first known outbreak in December 2019, the WHO reported over 6.28 million COVID-19 deaths and over half a billion confirmed COVID-19 cases (World Health Organization, 2021), while the total number of people infected with COVID-19 is believed to be much greater. It may take many years from now to fully ascertain the health burden and socioeconomic impact of the pandemic.

The COVID-19 pandemic has included multiple waves. These waves are mainly due to the implementation and relaxation of non-pharmaceutical interventions and the emergence of new variants of concern (VOCs). These VOCs have been observed worldwide. A COVID-19 variant contains one or more mutations in its viral genome. Certain COVID-19 variants have higher transmissibility and severity in populations than other COVID-19 variants. Emerging VOCs are those variants that are considered to have a distinguishable and significant health impact. Global travel and the timing of non-pharmaceutical interventions has made it difficult to determine when the importation of new VOCs may enter a given geographic region. The Alpha, Beta, and Gamma variants were first detected between October and November 2020 in the United Kingdom, South Africa, and Brazil, respectively (National Collaborating Centre for Infectious Diseases, 2022). These three variants then spread to other countries through global travel, and countries around the world experienced dissimilar transmission dynamics of these variants (Chen et al., 2021).

In Canada, these VOCs were first detected in December 2020 (Alpha), January 2021 (Beta), and March 2021 (Gamma). Each of these variants contributed differently to the transmission dynamics observed in the third wave. More recently, the emergence and spread of the Omicron variant (BA.1 and BA.2) was a major driver of the fifth and sixth waves in Canada. An in depth understanding of the COVID-19 VOCs transmission dynamics during the previous waves can provide valuable new insights on how to effectively prevent and control future waves of VOCs.

Mathematical modelling has been widely used as a research and policy tool. During the COVID-19 pandemic, academic researchers and government agencies worldwide have used mathematical models to help understand the spatial and temporal dynamics of the COVID-19 disease. These models have incorporated dynamical features motivated by disease dynamics, emerging variants, vaccine dynamics, and public health policies (Wu et al., 2139; Wu et al., 2021; Layton & Sadria, 2022; Yuan et al., 1101; Auger & Moussaoui, 2021; Callaway, 2020; Korber et al., 2020). The public health policies included the following strategies: travel restrictions, social distancing, isolation, testing for cases, and contact tracing. COVID-19 models have provided valuable insights and evidence that helped inform policies in a continuously changing pandemic.

A common feature of VOCs of COVID-19 has been their increased transmissibility, and they are expected to have a higher basic reproduction number. By the standard theory of the multi-strain competition (Andreasen, Lin, & Levin, 1997), under general assumptions of strain competition, the emerging VOC with the highest basic reproduction number will be able to invade a population and replaces the wild-type or existing variants according to the competitive exclusion principle in ecology (Andreasen et al., 1997; Bremermann & Thieme, 1989; Fath & Jørgensen, 2018; Levin, 1970). Several mechanisms for co-existence of strains of the pathogen in a population have been established in the literature, including super-infection and co-infection, mutation of one strain to another, cross immunity among strains, and population age heterogeneity in which different strains preferentially infect different age groups (see (Martcheva, 2015) for reviews). Time-periodic infection rates caused by seasonality and environmental influences have also been shown to lead to strain co-existence (Martcheva, 2009). These mechanisms for co-existence are generally not applicable to COVID-19 epidemics, and we assume the competitive exclusion principle holds for COVID-19 VOCs. The mathematical theories of multi-strain competition and competitive exclusion are based on the asymptotic behaviours of solutions to mathematical models when time is infinitely large. For finite-time horizon real-world epidemics such as epidemic waves of COVID-19, variant and strain dominance and coexistence often are not as clearly defined as in the theory of competition. As were shown in (Wu et al., 2139; Wu et al., 2021; Layton & Sadria, 2022), public health interventions can play a key role in mitigating, and possibly preventing, an emerging VOC from becoming dominant in the population.

The main objective of our study was to understand the dynamics of variant competition using a two-variant mathematical model and to interpret these behaviours in the context of emerging COVID-19 VOCs. The overall analysis will include both asymptotic behaviours based on the theory of mathematical epidemiology, and finite-time dynamics during a single epidemic wave using public health COVID-19 data from the provinces of Alberta and British Columbia, Canada. Our investigation focused on the following questions:

How to interpret dominance and coexistence of variants within the finite-time horizon of a COVID-19 wave, in comparison to the infinite-time horizon of the asymptotic limits?

What are the different characteristics of variant dominance or coexistence in finite-time and infinite-time horizons?

How to distinguish between dominance and coexistence of variants during a COVID-19 wave?

What public health measures can be implemented to prevent emerging variants from becoming dominant and/or mitigate the spread and size of the resulting epidemic?

Multiple variant dynamics of COVID-19 have not been widely discussed in the modeling literature, especially in the context of the finite-time dynamics during an epidemic or a single epidemic wave. With both an asymptotic and finite-time perspective, our study enriches the theory of variant/strain competition and it provides actionable insights to public health interventions related to preventing, mitigating, and managing emerging COVID-19 VOCs.

In the next section, we illustrate the derivation of a general two-variant model for infectious diseases that includes COVID-19 as a special case. Section 3 provides a detailed mathematical analysis of the model and Section 4 presents numerical simulation results using public health data from Alberta, Canada.

Derivation of the model

We developed an SIRS type of compartmental model for the transmission of two viral variants in the population. Historically, COVID-19 in Canada has presented with a clearly identified dominant variant and a newly emerging VOC. This has motivated our consideration of using a two variant rather than a multi-variant model. Furthermore, the mathematical analysis was simplified using a two-variant model. To further reduce the technicality in the mathematical analysis, we considered a model of SIRS type. More complex models such as SEIR and SEIAR types that include latent and asymptomatic compartments have been used for COVID-19 dynamics. In (Roda, Varughese, Han, & Li, 2020), the authors have shown that SIR models perform better than SEIR or more complex compartmental models to represent the public health data of COVID-19.

Our mathematical model has four compartments: number of individuals susceptible to the viral infection at time t, S(t); number of individuals infected with variants 1 and 2 at time t and not detected by the public health surveillance, I 1(t) and I 2(t); and number of individuals recovered (detected or undetected from testing) and remain protected against infection at time t, R(t). Individuals detected from testing transition to R and these individuals do not contribute to further transmission. Undetected individuals transition to R based on an average infectious period and these individuals also do not contribute to further transmission.

For COVID-19, the evidence of an individual having co-infection from two or more variants has been rare, with only a number of case reports of co-infection (see (da Silva Francisco et al., 2021; BBC News, 2022; Rockett et al., 2745; 1016Vatteroni, Capria, Spezia, Frateschi, & Pistello, ; Bolze et al., 1101)). Considering the evidence of co-infection during the COVID-19 pandemic, we assumed that co-infection from two or more variants is negligible and we have a single recovered compartment, R, in the model. The transmission coefficients for variants 1 and 2 are given by β 1 and β 2, respectively. Each transmission coefficient is a product of two factors: the average contact rate between the susceptible and infected individuals, and the probability of transmission per contact. During COVID-19, social-distancing measures and face masking were aimed at reducing the contact rate among individuals and the probability of transmission, respectively. The parameters γ 1 and γ 2 are the natural recovery rates from variant 1 infection and variant 2 infection, respectively, for infected people who are not detected by the public health surveillance. Parameters ρ 1 and ρ 2 are case-infection ratios for variants 1 and 2, respectively. A COVID-19 case is an individual diagnosed as COVID-19 positive by a PCR test as is recorded in the public health surveillance system. An undetected infection is an individual infected with COVID-19 but not recorded in the public health surveillance system. The case-infection ratio is defined as the number of new cases divided by the number of people living with undetected COVID-19 infections, or the so-called hidden infections. The ratio 1/ρ , i = 1, 2, measures the number of undetected COVID-19 infections by variant i in the community for each newly diagnosed variant i case, and it is a measure of the effectiveness of public health surveillance. Daily positive cases identified through testing is denoted by ρ 1 I 1(t) and ρ 2 I 2(t) for variants 1 and 2, respectively. Once detected through testing, individuals are assumed to not infect others until recovered, which is analogous to isolation requirements and/or reduced social interactions while ill. The rate at which individuals lose immunity is given by the parameter δ. The parameter Λ is the influx of susceptible individuals from birth and migration. When modeling a short duration epidemic, such as a COVID-19 wave, births and baseline deaths are often negligible, and, in this case, the parameter Λ is set to 0. Since the COVID-19 pandemic is transitioning to an endemic state, the assumption of Λ > 0 is helpful to assess the long-term effects of the infection. Also, with the consideration of the COVID-19 endemic state, we assumed the death rates (baseline or infection-related) are positive.

The model is depicted in the transfer diagram in Fig. 1 . Model parameters are listed in Table 1 , together with their biological meanings. Based on our assumptions and the transfer diagram in Fig. 1, the following system of differential equations can be derived for the model:

Fig. 1

Transfer diagram for the SIRS model for two variants.

Table 1

List of model parameters and their biological interpretations.

Parameter	Description
β1	transmission rate of variant 1
β2	transmission rate of variant 2
γ1	recovery rate of variant 1
γ2	recovery rate of variant 2
ρi	case-infection ratio for variant i, i = 1, 2
δ	rate of immunity loss
dS	background death rate of susceptible people
d_I_1	death rate of infected people (variant 1)
d_I_2	death rate of infected people (variant 2)
dR	death rate of recovered people
S0	initial susceptible population size
I01	initial number of infections for variant 1
I02	initial number of infections for variant 2
R0	initial people immune

Model analysis

By examining the direction of the vector field of system (1) on coordinate subspaces of , we can verify that the nonnegative orthant is positively invariant under the flow of model (1), namely, solutions with nonnegative initial conditions will remain nonnegative, and the model is well-posed.

Adding all equations in system (1) leads towhere . This implies that lim sup(S(t) + I 1(t) + I 2(t) + R(t)) ≤ Λ/d. We study system (1) in the following feasible region:which is positively invariant and contains the global attractor of model (1) in .

Equilibria and stability analysis

Model (1) always has the disease-free equilibrium . Since waning immunity is included, when the infection is not present, previously acquired immunity will be lost at the equilibrium P 0 and the entire population will be susceptible. There are two possible single-variant equilibria: and , where

The single-variant equilibrium P 1 exists in Γ ifit is outside of if , and it coincides with P 0 if . Similarly, P 2 exists ifit is outside of if , and it coincides with P 0 if .

We note that, for i = 1, 2, threshold parameter is the basic reproduction number for the variant i when it is the only variant present in the population. Let

Then is the basic reproduction number for the two-variant model (1), namely, it measures the average number of the secondary infections caused by a single infective with either variants during its entire infectious period.

Can a coexistence equilibrium exist when ? Assuming , then S* needs to satisfy

simultaneously. This is only possible if , namely, when . Furthermore, in such a case, infinitely choices of are possible as solutions of a linear equationand accordingly, a line segment defined by this equation connects P 1 and P 2 and consists entirely of positive equilibria (see Appendix for a proof).

Stability analysis

Local stability analysis of equilibria can be carried out using the method of linearization and the Routh-Hurwitz criteria. We state the following result that summarizes the existence and stability of the equilibria. Technical proofs are presented in the Appendix. It can be shown that similar results hold for SEIR or more complex models of this type.

Let, andbe defined in (3) - (6). The following statements hold.

If , then the disease-free equilibrium P 0 is the only equilibrium in the feasible region Γ and it is asymptotically stable.

If , then P 0 is unstable. Furthermore,

if , then the single-variant equilibrium exists and is asymptotically stable, while P 2 does not exist in Γ.

if , then both single-variant equilibria P 1 and P 2 exists in Γ. P 1 is asymptotically stable while P 2 is unstable in the direction pointing to the interior of Γ.

if , then the single-variant equilibrium exists and is asymptotically stable, while P 1 does not exist in Γ.

if , then both single-variant equilibria P 1 and P 2 exists in Γ. P 2 is asymptotically stable while P 1 is unstable in the direction pointing to the interior of Γ.

if , then both P 1 and P 2 exist, and there exists a line segment in Γ that consists entirely of positive equilibria and connects P 1 and P 2. Each positive equilibrium is neutrally stable in the direction of the line segment, and is asymptotically stable in directions transversal to the line segment.

Biologically, results in Theorem 1 on the existence and stability of equilibria infer that outcomes of the variants in the population are determined by the variant-specific reproduction numbers and as defined in (4), (5). As illustrated in Fig. 2 , in region I, relation holds, both variants are able to establish in the population. Variant 1 has the larger basic reproduction number and will become dominant and eventually drive the variant 2 to extinction; in region II, the reverse relation holds, and variant 2 will be dominant and drive variant 1 to extinction; and in region III, both and are less than 1, and neither variant can establish itself in the population and the disease dies out. This is consistent with .

Fig. 2

A diagram illustrating the results in Theorem 1. In region I, where holds, variant 1 will become dominant and eventually drive the variant 2 to extinction; in region II, where relation holds, variant 2 is dominant and drives variant 1 to extinction; and in region III, both and are less than 1, and neither variant can establish itself in the population and the disease dies out. On the half line defined by , both variants are able to coexist in the population.

We pay a special attention to the case in the diagram in Fig. 2 when . This is the only scenario under which both variants can coexist in the population under our model assumptions. Mathematically, a line of equilibria (case (e) in Theorem 1) is non-generic and not all equilibria on the line will survive under small perturbations. Furthermore, the half line given by has measure 0 in the 2-dimensional parameter region . These facts suggest that coexistence of the two variants in the sense of positive asymptotic limits is unlikely.

In real-world epidemics such as the COVID-19, dominance and coexistence of variants are often discussed within a finite time horizon (e.g. a single epidemic wave or pandemic) rather than an infinitely long time (or within asymptotic limits). In this paper, we use practical and theoretical dominance and coexistence to distinguish between finite and infinite time analysis and observations.

Public health data during the COVID-19 pandemic provides real-world examples of practical dominance and coexistence between variants and the original strain. Fig. 3 and Fig. 4 describe incident cases and percentage contributions of cases by VOCs for COVID-19 in British Columbia and Alberta (Jan 1 to May 30, 2021), respectively. In Fig. 3, the Alpha variant and Gamma variant appear to demonstrate practical coexistence during the second COVID-19 wave, having comparable levels of variant-specific incident cases and percent contributions. In contrast, the Beta variant was not able to establish itself in the population during the same period, and both Alpha and Gamma variants showed practical dominance over the Beta variant. Based on our model analysis and results in Theorem 1, the VOC data from British Columbia in Fig. 3 suggests that the reproduction numbers of Alpha and Gamma variants are similar, and both are larger than the reproduction number of the Beta variant.

Fig. 3

Weekly reported public health data in the Province of British Columbia, Canada, shows that variants Alpha and Gamma are able to coexist at comparable levels, while variant Beta was not able to establish itself in the population. The data covers the period from January 3 to June 20, 2021. Source of data: http://www.bccdc.ca/health-info/diseases-conditions/covid-19/about-covid-19/variants, accessed on June 25, 2021.

Fig. 4

Daily reports of public health data in the Province of Alberta, Canada, show that the Alpha variant dominated the Gamma variant in both (a) case numbers and (b) case percentages. This is in clear contrast to the variants situation in British Columbia as shown in Fig. 3. The data covers the period from December 15 2020 to June 15, 2021. Source of data: https://www.alberta.ca/stats/covid-19-alberta-statistics.htm#variants-of-concern, accessed on June 25, 2021. Note: the drop in the top curve in (b) is artificial and was due to the temporary stoppage of typing of variants during that period.

In Fig. 4, the VOC data for Alberta during the same time period as Fig. 3 describes a different situation. The Alpha variant showed a practical dominance over both Beta and Gamma variants in terms of incident cases and percentage case contributions. Although the VOC data from British Columbia showed practical coexistence of the Alpha and Gamma variant suggesting similar reproduction numbers, this was not observed in Alberta during the same time period. Alberta experienced practical dominance of the Alpha variant following the original strain. How can the mathematical theory of competition of variants be used to explain the apparently different situations of VOCs shown in Fig. 3, Fig. 4? Can variant-specific public health interventions explain the different situations in the two provinces? We carried out numerical investigations to address these questions in Section 4.

Numerical investigations and implications for endemic states of the COVID-19 pandemic

Numerical simulations were carried out using model (1) to investigate questions related to the dominance and coexistence of variants raised in the previous two sections. Subsection 4.1 illustrates analytical results from numerical simulations of theoretical (long-term) dominance and coexistence as described in Theorem 1. Parameter values in these simulations were independent of time and they are described in figure captions. Subsection 4.2 focuses on simulations to demonstrate practical (short-term) dominance and coexistence of variants. In these simulations, parameter values are fitted to public health data from Alberta, Canada, and include time-dependency to account for phased implementation and relaxation of public health measures (including testing policy changes).

Theoretical dominance and coexistence of variants

Our analytic results in Theorem 1 show that the variant with a larger basic reproduction number will become dominant, and will drive the main circulating variant to extinction as time tends to infinity. Furthermore, when the two variants have the same basic reproduction number, they can coexist as time tends to infinity. These long-term behaviours are termed theoretical dominance and coexistence dynamics. To illustrate the concepts, numerical simulations of number of the number (right) and percentage (left) of hidden infections by each of the variants are shown in Fig. 5, Fig. 6 .

Fig. 5

Simulations of model (1) that demonstrate theoretical dominance of variant 1 by variant 2 when . Parameter values used for simulations are β2 = 1.5β1, γ2 = γ2 = 0.1, ρ2 = ρ1, I01 = 9000 and I02 = 100. We note that even when the initial number of infected of variant 2 is smaller than that of variant 1, variant 2 still become dominant in the long term because of its basic reproduction number is larger.

Fig. 6

Simulations of model (1) demonstrating the coexistence of variants 1 and 2 when . Parameter values used in the simulations are β2 = β2, γ2 = γ2 = 0.1, ρ2 = ρ1, I01 = 9000 and I02 = 100. Simulation results show that the variant 2 with smaller initial condition I02 has a smaller limit in both numbers (a) and percentages (b).

Impact of on the theoretical dominance of variants

We assume that γ 1 = γ 2 and ρ 1 = ρ 2, and then select β 2 = 1.5β 1 so that . By Theorem 1, we expect that variant 2 will dominate variant 1, irrespective of their initial values I 10 > 0 and I 20 > 0, and that , and I 1(t) → 0 as t → ∞ (see Fig. 5). From the simulation results, we observe the following:

The number of infection I 2(t) is much larger than I 1(t) during the course of the epidemic, suggesting that variant 2 dominates variant 1 (Fig. 5 (a)).

The asymptotic limit in Fig. 5 (a) appears to be very close to 0 in comparison to the peak value of I 2(t) because of the scale, the infected populations I 1(t) and I 2(t) both appear to converge to 0 at the end of the epidemic. This leaves doubt about the dominance of variant 2 based on simulations of I 1(t) and I 2(t).

The percentage contributions of I 1(t) and I 2(t) in Fig. 5 (b) demonstrate a clear dominance of variant 2 over variant 1.

In conclusion, using percentage contributions of I (t) instead of their numbers can better identify the dominant variant for analyzing both numerical simulations (Fig. 5) and public-heath data (see Fig. 3, Fig. 4).

Theoretical coexistence of VOCs

In Fig. 6, the parameter values included the following assumptions: γ 1 = γ 2, ρ 1 = ρ 2, and β 1 = β 2 so that . Based on Theorem 1, we expect that both variants will coexist, and I 1(t), I 2(t) both converge to positive limits as t → ∞. Furthermore, the limits are determined by the initial conditions. From our simulation results in Fig. 6 we observed the following:

In Fig. 6 (a), I 1(t) is much larger than I 2(t) during the course of the epidemic, and the variant 1 appear to dominate variant 2.

In Fig. 6 (a), variant 1 has a larger initial condition than variant 2, and variant 1 produces a much larger epidemic than variant 2. This illustrates the conclusion (e) of Theorem 1 that the limits of number of infected for coexisting variants are dependent on the initial conditions.

In Fig. 6 (b), the percentage contributions of I 1 and I 2 illustrate that both variants take a positive percentage throughout the epidemic. In such a case, we say that variant 2 has a low-level coexistence with variant 1.

Overall, two variants coexist when . In this case, the variant with a greater number of initially infected will have a higher epidemic curve and percentage contribution of total infections. The significance of this theoretical result for real-world epidemics is that when there are multiple variants having similar reproductive numbers (e.g. Omicron BA.4 and BA.5) at an early stage within a similar time period, the variant with the largest initial infections will have the greatest percentage contribution of total infections. In situations where VOCs differ in their severity in terms of symptoms and case-fatality, control measures during the early stage of an epidemic or wave can be designed to limit the importation and spread of a more severe variant. If control measures cannot be implemented that target all circulating variants, suppressing the most severe variant would allow for a milder variant to take a higher percentage contribution and help reduce severe outcomes.

Practical dominance and coexistence of VOCs during the COVID-19 pandemic

In this subsection, we will examine the concepts of dominance and coexistence of new variants during real-world epidemics, or epidemic waves that typically only last for only a finite time. The COVID-19 pandemic provided an ideal context for such a study, since each epidemic wave was caused by a new variant that became dominant during the epidemic wave. We will use public health data in the provinces of Alberta and British Columbia, Canada, to inform the model simulations.

The main objective of the following numerical investigations is to provide plausible explanations for differences in transmission dynamics in finite time (a few months) between Alpha and Gamma variants in Alberta and British Columbia, as shown in Figs. 3 and 4. In these simulations, it is important to determine conditions that allow for variant dominance or the coexistence of two variants during an epidemic wave. Since these dynamics are challenging to observe through infections alone (see Fig. 5, Fig. 6), percentage contributions of variants were also provided. This is observed in Alberta's variants data in Fig. 4, where, the Alpha variant appears to dominate the Gamma variant in the number of cases (Fig. 4 (a)). However, it is clear from the percentage contribution by variant (Fig. 4 (b)) that the Gamma variant showed low-level coexistence with the Alpha variant. This shows the percentage contribution by variant is a better gauge for variant coexistence than incident case numbers.

To realistically model the COVID-19 dynamic, the incorporation of public health interventions such as social distancing, lock-down measures, testing, quarantine, isolation, and contact tracing are important. Some key parameters require time dependent characteristics to reflect the policy changes at different phases of the pandemic. For example, time-dependent transmission coefficients β (t) will reflect impacts on transmission such as social distancing and lock-down measures. Time-dependent ρ (t) will account for the impacts of COVID-19 testing, contact tracing, quarantine, and isolation measures. Accordingly, we will carry out our simulations using the following modified model:

We note that the influx of susceptibles Λ and non-COVID-19 death are set to zero since the epidemic only lasted three months and the impact of birth and background death on the size of the susceptible population were negligible during epidemic.

Determination of β(t) and ρ(t) in the model

A step-wise function used in time-dependent transmission parameters, β (t), i = 1, 2, represented the easing lock-down measures in Alberta between January and June 2021. The baseline transmission value was obtained from the endpoints of prior modeling results based on Alberta public health data in January 2021. One step-wise increase of 30% was introduced 18 days after the simulation start date of January 25, 2021. The simulation end date was June 30, 2021.

The time-dependent case-infection ratio ρ (t), i = 1, 2, represented the effects of population health-seeking behaviours and behavioural change during the COVID-19 pandemic. They were generated and scaled within the simulations. The minimum and maximum values for ρ (t) was obtained from prior fitting results between March and May 2020, informed by the testing and health link call data from Alberta Health Services. During each simulation, a burn-in run would help generate a case detection curve, which would be scaled such that the minimum and maximum values would be limited to previous fitting results. The minimum value was limited between 2.5 × 10−2 and 3.8 × 10−2. The maximum value was limited between 6.6 × 10−2 and 7.9 × 10−2.

The Affine Invariant Ensemble Markov Chain Monte Carlo (MCMC) algorithm (Roda, 2020; Roda et al., 2020) was the calibration procedure used to estimate time dependent terms ρ 1(t)I 1(t) and ρ 2(t)I 2(t), baseline transmission rates, the infectious period, and initial conditions (I 1, I 2, S and R) using confirmed cases for variant 1 and 2 on day t, respectively. The Matlab 2020a software was used to run the calibration procedures. Prior distributions of parameters at the start of simulations (January 25, 2021) were informed from epidemiological information and previous calibration results conducted before January 2022.

Dominance of alpha variant over the beta variant during the third wave in alberta and British Columbia

From Fig. 3, Fig. 4, it is apparent that the Alpha variant dominated the Beta variant during the third COVID-19 waves in Alberta and British Columbia, in both case number and case percentage. This suggested that the basic reproduction number of the Alpha variant was sufficiently larger than that of the Beta variant in both provinces.

Based on the numerical simulations using model (7) (see Fig. 7 ), results showed that when the basic reproduction number of the Alpha variant (variant 2) is 50% larger than the Beta variant (variant 1), namely, , the Alpha variant will dominate the Beta variant even when the Alpha variant had a smaller initial number of infected people.

Fig. 7

Simulations of model (7) demonstrating the dominance of the Beta variant (variant 1) by the Alpha variant (variant 2) in Alberta as observed in Fig. 4, with a relation . Parameter values used for simulations are β2 = 1.5β2, γ2 = γ2 = 0.1, ρ2 = ρ1, I01 = 9000 and I02 = 100. Even when the initial number of the infected is much lower for the Alpha variant, its sufficiently larger basic reproduction number allows the Alpha variant to become dominant.

Practical coexistence of the Alpha and Gamma variants during the third wave of COVID-19 epidemic in British Columbia

Fig. 3 describes VOC trends for COVID-19 from British Columbia during the third wave. Case numbers and percentage contributions by variant type show that the Alpha and Gamma variants coexisted during the third wave and dominated the Beta variant.

From our theoretical results in Theorem 1, two variants can coexist long-term if and only if they have the same basic reproduction number. The observed coexistence of the Alpha and Gamma variants during the third COVID-19 wave in British Columbia suggests that both the Alpha and Gamma variants may have similar basic reproduction numbers. We performed numerical simulations on the model (7) to verify this possibility. From our simulation results in Fig. 8 , we observed the following:

Fig. 8

Simulation results of model (7) demonstrating the coexistence of the Alpha variant (variant 1) and Gamma variant (variant 2) when they have similar basic reproduction numbers , as can be observed in both (a) case numbers and (b) case percentages. Parameter values used in the simulations are β2 = 1.06β2, γ2 = γ2 = 0.1, ρ2 = ρ1, I01 = 9000 and I02 = 2000.

When variant 1 (Alpha) and variant 2 (Gamma) have similar basic reproduction numbers, , they can coexist, in both the infection numbers and percentage contributions (see Fig. 8).

In Fig. 8, variant 2 had a slightly larger basic reproduction number but a smaller number of initially infected people compared to variant 1. Based on this simulation, variant 1 had greater number of infections and percentage contributions compared to Variant 2. This suggested that when the basic reproduction numbers are comparable between two variants, the variant with the more infected people early on will produce a larger epidemic wave. This agrees with our observations of theoretical coexistence in Sub-section 4.1.2. The use of real data asserts the importance of early public health responses during an epidemic given an emerging variant.

Early public health responses given an emerging variant can help control replacement dynamics in a population: Alberta's experience with COVID-19 variants, Alpha and Gamma

During the third wave, the Gamma and Alpha variant coexisted in British Columbia, while in Alberta, coexistence of the Gamma variant occurred at low-levels with the Alpha variant dominating the epidemic wave (see Fig. 3, Fig. 4). The practical coexistence of the Alpha and Gamma variants in both provinces suggest that the two variants have similar basic reproduction numbers. In subsection 4.1.2, simulations showed that the variant with more initially infected people can lead to a larger epidemic size (see Fig. 6). The following simulations aim to address differences in public health responses that may have resulted in variant dynamics observed in both provinces.

As outbreaks of the Gamma variant emerged in Alberta, enhanced public health responses were implemented to contain its spread through increased testing, contact tracing, and isolation measures ((Hinshaw, 2021), (Lachacz, 2022)). Can variant-specific public health measures explain the low-level coexistence of the Gamma variant in Alberta? We carried out numerical simulations to investigate the impacts of containment measures (i.e. testing, contact tracing, and isolation) on variant dynamics.

These simulations aimed to investigate the impacts of increasing the case-infection ratio ρ (t), which incorporates the effects of case detection through testing, contact tracing, and isolation. Variant 1 (Alpha) and variant 2 (Gamma) had similar reproduction numbers since they were able to coexist in both provinces. To reflect enhanced surveillance of the emerging variant 2, the case-infection ratio ρ 1(t) was increased to two times of ρ 1(t). In Fig. 9 , simulations show low-level coexistence of the variant 1 and variant 2, with the latter was at much lower levels for both infection numbers and percentage contributions. An important result in Fig. 9 was that while variant 1 and variant 2 had a similar basic reproduction number, the targeted increase in testing, contact tracing, and isolation of the emerging variant was effective to mitigate the total number of infections.

Fig. 9

Simulation results of model (7) demonstrating that variant-specific public health interventions can prevent the Gamma variant from taking hold in Alberta. We have assumed that the Gamma variant (variant 2) has a slightly higher transmission rate β2 than the transmission rate β1 of the Alpha variant (variant 1), β2 = 1.17β1. We also assumed that the Gamma variant has a higher ρ2(t) than the Alpha variant, ρ2(t) = 2ρ1(t). The initial number of infected are chosen as I01 = 9000 and I02 = 100. The choices of ρ(t) and I0 reflect the additional effort in testing, DNA typing, contact tracing, and isolation directed at cases of the Gamma variants implemented in Alberta.

These simulations results highlighted that when an emerging variant has a similar reproduction number as the existing variant, early public health measures that increase case-infection ratios (i.e. testing, contact tracing, and isolation) of the emerging variant can prevent its replacement and/or allow for low level coexistence in the population. In addition, if the emerging variant poses a greater population-level risk of severity, this targeted approach can be effective to reduce the level of severe disease outcomes in the community. This was particularly evident in Alberta, where targeted public health measures that increased case-infection-ratios for the emerging Gamma variant was effective at mitigating its total infections and percent contributions (Fig. 4).

Summary and discussions

We investigated disease transmission dynamics of two variants (strains) using deterministic modelling and focused on coexistence of variants and dominance of one variant over the other. Analysis and simulations considered asymptotic limits (termed theoretical coexistence) and finite-time horizon for e.g. within a single epidemic wave (termed practical coexistence). While the simulations focused on COVID-19 specific data, the lessons learned can offer insights for general infectious disease epidemics and public health interventions.

The theoretical concept of coexistence of variants and dominance of a variant in mathematical epidemiology is based on the asymptotic limit in an infinite time horizon. The outcomes are determined by the variant-specific basic reproduction numbers, irrespective of the respective number of initially infected individuals. In the case of simple competition, the principle of competitive exclusion prevails: a variant with the largest basic reproduction number dominates and drive the other variants to extinction. The theory focuses on the asymptotic limits without providing much information about finite-time relations among the variants, and it is less useful when applied to actual epidemics such as seasonal influenza and epidemic waves of COVID-19. One useful insight from this theory based on infinite-time horizon is the possibility that variants can coexist when their reproduction numbers are similar and when they emerge concurrently. When applied to real-world epidemics, this can help us to infer that two variants with comparable number of confirmed cases should have similar basic reproduction numbers. In Fig. 8, we similar reproductive number and different initial conditions for infectious people.

Real-world epidemics typically last for only a short period of time and the number of infected caused by each variant will rise at the start, peak, then decline to low levels. Assessing patterns of coexistence or dominance in simulations was challenging using the number of infections alone. Percentage contributions was a better outcome to assess coexistence and dominance relations. Unlike the asymptotic theory of variant dominance that depends solely on the variant-specific basic reproductive numbers, in finite-time variant dynamics, the number of initially infected individuals can also play an important role to determine which variant may cause an higher level of epidemic or a greater final size.

The dynamics of the Alpha and Gamma variant in Alberta and British Columbia during the third wave was of particular interest in this analysis since coexistence of both variants were observed in both provinces. However, the percent contributions of these variants were shared more equally in British Columbia (Hogan et al., 2021) compared to Alberta (Figs. 3 and 4). Despite magnitude differences, the overall coexistence observed indicated that both variants likely had similar basic reproduction numbers based on theoretical insights from mathematical epidemiology. Despite comparable reproductive numbers for both variants, early and enhanced containment measures (e.g. testing, contact tracing, and/or isolation) targeted on the emerging Gamma variant may have also been more effective in keeping infectious cases at lower levels in Alberta. Given variants with comparable reproductive numbers, there may be some interaction between initial size of infectious people (i.e. initial number and size of outbreaks) (Fig. 8) and the effectiveness of containment measures (including asymptomatic testing (Fraser, Riley, Anderson, & Ferguson, 2004)) (Fig. 9) that could also explain differences in variant dynamics between Alberta and British Columbia.

Overall, the mathematical analysis and simulations of real-world examples such as COVID-19 provided valuable insights for future events involving variants emerging at similar times with comparable reproductive numbers. Earlier containment measures that effectively target emerging variants that pose a greater risk for severity can impact finite-time variant replacement dynamics by forcing the more severe variant to coexist at a much lower level in a population compared to less severe variants. While differences in percent contributions of the Alpha and Gamma variants between Alberta and British Columbia were described using modeling, there may be other factors not captured in the mathematical model that could play a role such as vaccine coverage and importation through travel. In Canada, as we transition from the pandemic to endemic phase of COVID-19, the expectation of emerging VOCs is a reality. These valuable insights offered from the analysis of past waves can help with the future management of COVID-19.

Declaration of competing interest

The authors claim no conflicts of interest.