Dynamic of a two-strain COVID-19 model with vaccination.

COVID-19 is a respiratory illness caused by an ribonucleic acid (RNA) virus prone to mutations. In December 2020, variants with different characteristics that could affect transmissibility emerged around the world. To address this new dynamic of the disease, we formulate and analyze a mathematical model of a two-strain COVID-19 transmission dynamics with strain 1 vaccination. The model is theoretically analyzed and sufficient conditions for the stability of its equilibria are derived. In addition to the disease-free and endemic equilibria, the model also has single-strain 1 and strain 2 endemic equilibria. Using the center manifold theory, it is shown that the model does not exhibit the phenomenon of backward bifurcation, and global stability of the model equilibria are proved using various approaches. Simulations to support the model theoretical results are provided. We calculate the basic reproductive number R 1 and R 2 for both strains independently. Results indicate that - both strains will persist when R 1 > 1 and R 2 > 1 - Stain 2 could establish itself as the dominant strain if R 1 < 1 and R 2 > 1 , or when R 2 > R 1 > 1 . However, because of de novo herd immunity due to strain 1 vaccine efficacy and provided the initial stain 2 transmission threshold parameter R 2 is controlled to remain below unity, strain 2 will not establish itself/persist in the community.

Introduction

The potential for COVID-19 virus circulating inside bats to mutate to humans was noted globally [1]. COVID-19 is a deadly respiratory disease caused by the Severe Acute Respiratory Syndrome (SARS) virus with sustained human-to-human transmission since December 2019 when the first case of the novel virus was detected in Wuhan, China [2], [3]. COVID-19 has a general mortality rate below 5%, with an average of 2.3% [4], but the older populations is the higher risk group with mortality rate of 8% for individuals between 70–79 years and 14.8% for people older than 80 years [5]. Despite the seemingly low mortality rate, the number of hospitalizations is quite high, presenting global health burden and a major challenge to health care systems worldwide [6]. COVID-19 is transmitted from human-to-human through direct contact with contaminated objects or surfaces and through inhalation of respiratory droplets from both symptomatic and infectious humans [7], [8]. The 2019 COVID-19 outbreak is still ongoing and represents a serious challenge for communities around the globe, endangering the health of millions of people, and resulting in severe socioeconomic consequences due to lock-down measures. In fact, Usaini et al. [9] noted that reducing the influx of immigrants could play a significant role in decreasing the number of infected individuals when the recruitment rate of immigrants is below a certain critical value.

COVID-19 transmission dynamics models are flourishing and abound in the literature [10], [11], [12], [13], to cite a few and the reference therein. Availability of COVID-19 vaccines brings hope to the potential end of the pandemic [13]. Vaccines provide a determining pharmacological measure in the struggle against the COVID-19 pandemic, as we now face a very different epidemiological landscape from the early pandemic [14], thereby opening the possibility to explore real scenarios that combine the effects of both non-pharmaceutical public health interventions (e.g., face mask, hand washing, social distancing) and therapeutic measures such as treatment and mass vaccination strategies [15]. Compartmental models have been crucial to study the evolution of several disease outbreaks. COVID-19 outbreak has provided a platform for several research activities based on compartmental-like epidemic models have been conducted to investigate different key aspects of the spread, control, and mitigation of the disease [15]. Deterministic compartmental disease transmission models are characterized by the subdivision of the population into compartments based on individuals’ health status. The history of mathematical epidemiologic models date as far back as Bernouilli [16], [17], [18]. Mathematical models of a two-strain disease are numerous in the literature [19]: malaria [20], influenza [21], [22], SARS-CoV-2 [23], dengue [24], disease with age structure and super-infection [25], influenza with a single vaccination [26], [27] to cite but a few and the references therein.

While studies on the dynamics of two viral infections have considered cross-immunity and co-infection [28], others described effects of two competing strains characterized by cross-immunity [29]. In fact, it has been reported that the number of strains increases logarithmically the maximum number of infected individuals and the mean mortality rate [30]. Reinfection (potentially due to waning immunity) and multiple viral strains are some of the major challenges in the current COVID-19 pandemic [31]. Because some COVID-19 strains are thought to be more transmissible than the original strain, Yagan et al. [32] discussed the importance of incorporating mutations and evolutionary adaptations into epidemic models. Other studies of multi-strain epidemic model for COVID-19 did not include vaccination [33], [34], [35].

Our proposed model is a mirror of a multi strain (two-strain) dynamics flu model with a single vaccination in [27] and modified in [26] to include the force of infection in both infected compartments and extending the incidence function to a more general form. With COVID-19 specificity, we included infections from vaccinated individuals against strain 1 (the resident strain), as well as strain 2 (the wild strain), since vaccination against strain 1 may not procure any or very limited protection against the second and more virulent strain 2. Because variant strains have the potential to substantially alter transmission dynamics and vaccine efficacy, Gonzalez-Parra et al. [23] investigated the impact of more infectious strain of the transmission dynamics of the COVID-19 pandemic, but they did not consider vaccination. They concluded that a new variant with higher transmissibility may cause more devastating outcomes in the population. While Puga et al. [36] investigated co-circulation of two infectious diseases and the impact of vaccination against one of them, our proposed model is seemingly new and to the best of our knowledge, no COVID-19 modeling study has accounted for strain 1 vaccination with possibility of infection with strain 1 when vaccinated against the original strain 1 as well as infection with strain 2 for which strain 1 vaccination may not provide much protection.

This paper is organized as follows. We formulate a deterministic compartmental epidemic model of the transmission dynamics of COVID-19 in a homogeneously mixed population in Section “Model Formulation” . Section “Model analysis” is devoted to well-posedness of the proposed mathematical model, derivation of its equilibria, the basic reproduction number, and analysis using dynamical systems theory of the COVID-19 transmission dynamics with strain 1 vaccination. Section “Numerical simulations” covers several numerical simulations of the disease dynamics in the presence of strain 1 vaccination in a community where treatment is administered to infected individuals. The graphical illustrations are based on various scenarios when the basic reproduction number is either greater or less than unity. The conclusion is provided in Section “Conclusion”.

Model formulation

It is assumed that the population is homogeneous-mixed and individuals have equal probability of acquiring the infection. Only human-to-human transmission of COVID-19 is considered. According to individuals diseases status, the human population at time denoted by is divided into sub-populations of susceptible individuals , vaccinated individuals , individuals infected with strain 1 , individuals infected with strain 2 and recovered . The total human population is given by

Homogeneous mixing of individuals in the population is assumed so that the standard incidence (rate of infection of strain 1 per unit time) is . Thus, at any given time , the probability that an individuals will carry strain 1 infection is . Infected individuals with strain 1 either die naturally at a constant rate or at a constant disease-induced death rate . The per capita life expectancy is given by while is the death adjusted average infectious life of individuals infected with strain 1.

While at the onset of the COVID-19 pandemic, the dynamics of the disease was much faster than that of birth (or recruitment) and deaths [37] because of the then short period of the pandemic [38], neglecting these demographic factors was well justified in the plethora of models in the literature. However, since the disease has been there for a while (since December 2019), at present, it is important to account for the model vital dynamics when describing the evolution of the COVID-19 pandemic. Therefore, disease-specific death rate respectively for strain 1 and strain 2 and natural death are accounted for. We incorporate both cohort vaccination (where a fraction of the newly recruited members of the community are vaccinated), and continuous vaccination program (where a fraction of susceptible individuals is vaccinated per unit time) [39]. We note however that although COVID-19 vaccination is not yet given at birth, its incorporation into our theoretical study does not impact of the conclusions as it has been shown that some analytical/numerical modeling results are independent of the type of vaccination program adopted [40]. This results is ascertain by graphically showing the potential impact of the model parameters on the initial disease transmission and . With the resurgence of COVID-19 after the multiple waves and various strains, and as COVID-19 vaccines are still evolving with lower age groups being considered, our proposed model is proactive in the even that the vaccine development includes vaccination of the birth cohort.

From the conceptual model flow diagram in Fig. 1, we derive the following deterministic system of nonlinear differential equations with initial conditions where The model system (1) involves both exogenous parameters such as the vaccination rate , and the recovery/treatment rates and - the latter represent the inverse of the length in days of the contagious period, and endogenous parameters such as the disease transmission rate and . The parameter is the vaccine effectiveness. Thus, implementation of a vaccination program causes the following transformation in the model: , that is the reduction in acquiring strain 1 infection for individuals already vaccinated against strain 1.

Fig. 1

Flowchart of the state progression of individuals in a population exposed to two strains of COVID-19. At time , susceptible individuals can become infected (primary infection) with strain 1 or strain 2 or vaccinated against strain 1. Vaccinated individuals can acquire COVID-19 strain 2. Infected individuals recover from both strains and move into class . Recovery is not permanent.

The model parameters, their description, values and sources are provided in the Table 1.

Table 1

Fundamental model parameters.

Parameter	Description	Value	Unity	Range	Reference
Λ	Recruitment or inflow into the population	100059×365	H × day^−1		[23], [41]
v	Continuous strain 1 vaccination rate		day^−1	[10^−5,8×10^−2]	[42], [43]
a	Effective contact rate	0.85	day^−1		[44], [45]
β_1	Transmission probability of strain 1		1	[0.127,0.527]	[23], [46]
β_2	Transmission probability of strain 2		1	[0.127,0.527]	[46]
ɛ	Strain 1 vaccine efficacy	0.87	day^−1		[47], [48]
μ	Natural death rate	159×365	day^−1		[23], [46]
δ_i	Strain i={1,2} disease-induced death rate	6.83×10^−5	day^−1		[23], [44]
σ	Rate of loss of immunity	190	day^−1		[49]
ρ	Cohort vaccination rate		day^−1	(0,0.99]	[42], [43]
τ_1	Recovery rate of strain 1 infected individuals		day^−1	[130,14]	[44], [46], [50]
τ_2	Recovery rate of strain 2 infected individuals		day^−1	[130,14]	[44], [46], [50]

Model analysis

Disease-free equilibrium and basic reproduction number

For system (1) with non-negative initial values, its solutions are non-negative and ultimately bounded. The proof is routine, see for example [38]. Positivity is important for biologically feasible solutions of the model while boundedness implies that solutions are finite. Next, we show that the solutions of model system (1) enter in a bounded region .

The closed set is positively invariant and attracting.

By adding all the equations of the model system (1), we obtain:

Using the comparison theorem as described in [51], [52], we have . If then . Thus, the region is positively invariant for the model, while if , then the solution enter in the region in finite time or when .

Thus, the region attracts all solutions in . So the system is positively invariant and attracting. ■

The disease free equilibrium of the system (1) is given by where

The basic reproduction number is with , and .

Using the next generation matrix method in [53] the associated next generation matrix is given by: and the rate of transfer of individual to the compartments is given by:

Hence, the new infection terms and the remaining transfer terms are respectively given by: and

The dominant eigenvalue or spectral radius of the next generation matrix which represents the basic reproductive number is given by:

The basic reproductive number of a disease, denoted is defined as the average number of secondary infections that a single infectious individual will give rise to over the duration of his infection, in an otherwise entirely susceptible population.

Let and

Then

Because it has been posited that modeling results are independent of the type of vaccination program adopted [40], we now investigate the potential impact of the cohort vaccination rate on long-term dynamics of the disease. To identify and rank pivotal model parameters whose uncertainties contribute significantly to the model prediction, we use the partial rank correlation coefficients (PRCC) to identify and quantify the statistical significance of each model parameter with respect to and (since ). For more details, see [54]. We note as expected from Fig. 2, Fig. 3 that the cohort vaccination rate does not have a significant impact on the initial disease transmission. The most sensitive parameters are the effective contact rate , the transmission probabilities of strains 1 and 2, respectively and , vaccination rate and vaccine efficacy , and the recovery/treatment rates of strain 1 and strain 2 . Thus, to mitigate the spread of the disease, implementation of control measures that target contact reduction, vaccination and treatment will be most effective in curbing the spread of the disease. We also show in Appendix “Potential impact of cohort vaccination rate on the long-term dynamics of COVID-19” that cohort vaccination has almost no impact on the long-term dynamics of the disease, Fig. 14, Fig. 15, Fig. 16, Fig. 17. Though this might seems counter intuitive, it is no surprise because birth cohort have little or minimal interaction with the larger population.

Fig. 2

PRCCs showing the effect of varying the model input parameters on .

Fig. 3

PRCCs showing the effect of varying the model input parameters on .

Fig. 14

Time series of strain 1 infected individuals for and .

Fig. 15

Time series of strain 2 infected individuals for and .

Fig. 16

Long-term dynamics of strains 1 and 2 infected individuals and for .

Fig. 17

Long-term dynamics of strains 1 and 2 infected individuals and for .

In Section “Numerical simulations”, we shall investigate four possible scenarios/combinations when or and or .

The disease-free equilibrium is unstable if while it is locally asymptotically stable if .

The Jacobian matrix associated with the model system (1) at the disease-free equilibrium is given by is given in Box I.

Thus the eigenvalues of are for and .

If , then and we obtain that the disease-free equilibrium of model (1) is locally asymptotically stable. If , then the disease-free equilibrium loses its stability. ■

The disease-free equilibrium is globally asymptotically stable if .

Consider the Lyapunov function Since , then and attains zero at .

Now, we need to show that .

Because is diminished at any time due to infection and is also decreased due to imperfect vaccination, let and , then

Furthermore, if and only if , so by using LaSalle’s invariant principle, this implies that is globally asymptotically stable in . ■

Using the standard comparison theorem as described in [51], [52] and rigorously applied in [41], [55], [56], this result can also be proved (see Appendix “GAS of DFE using comparison theorem”).

Because the impact of the vaccination on disease dynamics is key to our study, we write and as functions of the vaccination rate . That is, Thus,

From the first equation in (5), the strain 1 basic reproduction number decreases as vaccination rate increase. That is, as expected vaccination against strain 1 is always beneficial in controlling strain 1, however its impact on strain 2 depends on the effective contact rate and the transmission probability of strain 2 , as individuals vaccinated against strain 1 may be less likely to be infected by strain 2 than those who are not vaccinated [27]. Large vaccination rate could potentially lead to reducing to some value less than 1, thereby helping to mitigate the spread of strain 2. Since is a function of daily vaccination rate and vaccine efficacy , its variation for different values of these parameters is shown in Fig. 4. Dark blue color corresponds to high vaccine coverage and efficacy, indicating the possibility to decrease the value of . That, both high vaccine coverage and efficacy will contribute to the reduction of the value .

Fig. 4

Heat map of for different values of vaccination rate and vaccine efficacy . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Endemic equilibrium

The model (1) has

a unique single-strain 1 infection equilibrium if and only if .

a unique single-strain 2 infection equilibrium if and only if .

a 2-strain infection equilibrium when .

(1) Equilibrium is the solution of the system

From the first three equations in (6) above, we have the expression of , and in function of .

Replacing all this expression in the last equation, we obtain the following expression of in function of . with and

So we can therefore have and according to

Using , replacing and by their expressions, then after some algebraic calculations, we obtain that is the solution of the equation where

When , then the discriminant of the quadratic Eq. (7) is given by is positive (because and have distinct signs), so Eq. (7) admits two real solutions. In addition, the product of those two solutions is , implying that the two solutions have different signs. Hence, we can conclude that when , the system admit a unique single-strain 1 infection equilibrium.

(2) The proof follows the same approach and steps as in the above case (1).

(3) To find , we consider the system

where , with . After some algebraic manipulations we obtain , , and .

Replacing ; and with their values in the third and fourth equation yields the following system where and are monotone functions defined by If the system (10) admits a solution, then, the model system (1) will have an endemic equilibrium. Obtaining the explicitly expression for the exact solution of the non-linear autonomous system (10) is a daunting task. Also, it is not obvious if the system (10) admits multiple solutions, it is therefore important to explore the uniqueness and global stability of the 2-strain endemic equilibrium. To this end, we investigate if the model system (1) undergoes the phenomenon of backward bifurcation where a stable disease-free equilibrium co-exists with a stable endemic equilibrium when .  ■

Bifurcation analysis

In determining the possibility of bifurcation, we use the Centre Manifold theory approach [57]. For simplification of the notations and ease of algebraic manipulations, the following change of variables are made. Let , , , and , by using the vector notation (where denote the transpose), our model system (1) can be written in the form with as follows: where and .

The Jacobian of the system (11) at the DFE is given by

First, consider the case .

Consider the case when , which is the bifurcation point. Suppose, further that is chosen as a bifurcation parameter. Solving for from gives , and the Jacobian denoted by is given by

When , the Jacobian of (11) at has a right eigenvector given by , where, , , , and

, where , , and

.

Further, the Jacobian has a left eigenvector , where

To determine the direction of the bifurcation, we compute and determine the sign of the two bifurcation coefficients and . At the DFE, the bifurcation coefficient is given by

After some algebraic computations, we obtain

Since , and for our chosen model parameter values in Table 1, it is therefore evident that , since , and .

The second bifurcation coefficient is given by Because is negative and is positive, by Theorem 4.1 in [57], this precludes the model system (1) from exhibiting the phenomenon of backward bifurcation at , that is co-existence of both a stable disease-free equilibrium with a stable endemic equilibrium cannot hold. Consequently, at the threshold parameter , a forward or transcritical bifurcation occurs as depicted in Fig. 5. Thus, the endemic equilibrium is unique, and one could construct a suitable Lyapunov function to prove its global asymptotically stability [40], [55], [58], [59].

Fig. 5

Forward bifurcation diagram.

Numerical simulations

To illustrate the basic mechanisms underlying the model dynamics, several graphical representations depicting the dynamical behavior of the model system (1) when the fundamental threshold parameter is either greater or less than unity are presented to support the analytical results. The model parameter values used in our simulations are shown in Table 1. The unit of is person per day while all other parameters’ unit is per day. We consider the following initial conditions: . Because , four scenarios will be considered. Note that values of our proposed COVID-19 model’s basic reproduction number (with and denoting respectively the strain 1 and strain 2 basic reproduction numbers) are in agreement with previous COVID-19 modeling studies [60], [61], [62].

Simulations when with and

Time series solution of two infected classes and are plotted in Fig. 6, Fig. 7 when with . In this case, the solutions of the model system (1) approach the equilibria . These graphs show the occurrence of a second wave as predicted in [10], and possibly a third wave. However, we note that under very pessimistic conditions with availability of only strain 1 vaccine, strain 2 could become the dominant strain in the population if infections with strain 2 are more than double that of strain 1, see Fig. 7.

Fig. 6

Simulations showing individuals and for and respectively.

Fig. 7

Simulations showing individuals and for and respectively.

Fig. 8, Fig. 9 depict the case when with and . As depicted in Fig. 9, at the long run, strain 2 could establish itself as the dominant strain in the population. Note that in this case, the solution profiles approach the equilibrium

Fig. 8

Simulations showing individuals and for and respectively .

Fig. 9

Simulations showing individuals and for and respectively.

Simulations when with (a) and and (b)

Fig. 10, Fig. 11 display respectively the cases when (a) with and (left panel) and (b) (right panel). Again, the dynamical behavior of the graph of strain 1 depicts multiple waves. In both cases (a) and (b), strain 2 will not establish itself in the population as the solutions approach the equilibrium . It is surprising to notice that when both strain 1 and strain 2 reproduction numbers are equal and greater than unity, strain 2 could eventually die out. Several reasons could explain this. First, strain 2 emerged in the population when efforts to mitigate the strain 1 such as non-pharmaceutical interventions (including physical distancing, hand hygiene, and mask-wearing) as well as treatment were already well underway. Secondly, continuous vaccination against strain 1 could likely confer some protection to individuals against strain 2.

Fig. 10

Simulations showing individuals and for respectively.

Fig. 11

Simulations showing individuals and for and respectively.

Simulations when

When , that is both and are less than unity with , the evolutionary dynamics of the solutions approach the disease-free equilibrium , see Fig. 12. On the other hand, when with , it is striking to note that because is very closed to 1, Fig. 13, the strain 1 infection will not be quickly eradicated, and efforts to further reduce the average number of infections is warranted.

Fig. 12

Simulations showing individuals and for and respectively.

Fig. 13

Simulations showing individuals and for and respectively.

Conclusion

COVID-19 emerged in December 2019, has rapidly evolved as a pandemic with wide-ranging socio-economic consequences. The disease causes severe acute respiratory syndrome and results in substantial morbidity and mortality. While an effective vaccine is essential to containing the spread of COVID-19, the emergence of a second strain could complicate mitigation efforts. We developed a simple compartment model of the transmission dynamics of a 2-strain COVID-19 model to examine the impact of strain 2 in a population where vaccination against strain 1 is available. The proposed 2-strain COVID-19 model with strain 1 vaccination is derived as a deterministic system of nonlinear differential equations. The model is then theoretically analyzed, its basic reproduction number is derived as well as sufficient conditions for the stability of its equilibria. We calculate the basic reproductive numbers and for both strains independently. Using the center manifold theory, it is shown that the model does not exhibit bi-stability also known as backward bifurcation, but a forward bifurcation occurs at . Thus, global stability of the disease-free equilibrium when , and of the endemic equilibrium when is established using suitably constructed Lyapunov functions and other approaches such as the comparison method.

To gain insight into whether strain 2 will establish itself in the population as the dominant strain, several simulations to support the model theoretical results are provided. Results indicate that – both strains will persist when both and – Stain 2 could likely establish itself as the dominant strain if and , or when is at least two times . However, with the current knowledge of the epidemiology of the COVID-19 pandemic and the availability of treatment and effective vaccines against strain 1, strain 2 would eventually be eradicated in the population if the threshold parameter is controlled to remain below unity. That is, two co-circulating strains will not persist simultaneously but only one of the strains may persist in the long run. We note that if we ignore the model vital dynamics (recruitment and natural death) to mimic the early days of the epidemic, there is no noticeable impact on the dynamical behavior of the disease. There are however some contrasting findings with respect to the value of the basic reproduction number.

Under a very pessimistic condition, strain 2 could become the dominant strain in the population if infections with strain 2 are more than double that of strain 1, see Fig. 7.

When both strain 1 and strain 2 reproduction numbers are equal and greater than unity, strain 2 could eventually die out while strain 1 persists.

From observation (ii), while provides a good measure for disease dying out or persisting in a population, this threshold quantity when less than but close to unity might mislead the assessment of the transmission dynamics of the disease. Thus, ensuring that the value of is below unity may depend on how far from unity this value is in order to ascertain how quickly the disease eventually dies out. That is, even though the model does not exhibit the possibility of bi-stable behavior of its equilibria, strain 1 could persist for some time when , but close to 1. From this finding and in the face of waning adherence to physical distancing, the use of non-pharmaceutical interventions the word has relied upon (such as lock-downs, travel restrictions, contact tracing, mask wearing, and social/physical distancing), and the emergence of other COVID-19 variants, it is cautionary to ensure decision on relaxing/lifting these non-pharmaceutical prevention measures are not solely based on the value of the basic reproduction number being less than unity, but considerations should be made on how close to 1 this value actually is as well as other socio- and eco-epidemiological factors pertaining to the dynamics of COVID-19, and also account for regional heterogeneity in transmission and travel. Nevertheless, there is a glimpse of hope that if individuals concurrently continue to adhere to non-pharmaceutical interventions (including physical distancing, hand hygiene, and mask-wearing), and other pharmaceutical efforts to mitigate the strain 1 such as treatment and vaccination continue, strain 2 could be eradicated.

The proposed model has some limitations. While co-infection of the COVID-19 has not been a major issue, from a theoretical standpoint, the model could be extended to include the latent class, and individuals dually infected with both strain 1 and strain 2. In this case, one could compute the invasion reproductive number for strain 1 when strain 2 is at endemic equilibrium and vice-versa [63]. Due to the severity of the disease, explicitly incorporating the quarantine and hospitalized class is viable. While these suggestions will increase the complexity of the model analysis, by construction, there are often uncertainty around some parameter values, and a detailed uncertainty and sensitivity analyses to determine the parameters that have the highest effect on the model variables is important to guide policy and health decision makers on which parameter to prioritize or to target first [64].

CRediT authorship contribution statement

S.Y. Tchoumi: Conceptualization, Formal analysis, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. H. Rwezaura: Formal analysis, Methodology, Validation, Visualization, Writing – original draft. J.M. Tchuenche: Conceptualization, Methodology, Validation, Visualization, Writing – review & editing, 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.