Analysis of multi-strain infection of vaccinated and recovered population through epidemic model: Application to COVID-19.

In this work, an innovative multi-strain SV EAIR epidemic model is developed for the study of the spread of a multi-strain infectious disease in a population infected by mutations of the disease. The population is assumed to be completely susceptible to n different variants of the disease, and those who are vaccinated and recovered from a specific strain k (k ≤ n) are immune to previous and present strains j = 1, 2, ⋯, k, but can still be infected by newer emerging strains j = k + 1, k + 2, ⋯, n. The model is designed to simulate the emergence and dissemination of viral strains. All the equilibrium points of the system are calculated and the conditions for existence and global stability of these points are investigated and used to answer the question as to whether it is possible for the population to have an endemic with more than one strain. An interesting result that shows that a strain with a reproduction number greater than one can still die out on the long run if a newer emerging strain has a greater reproduction number is verified numerically. The effect of vaccines on the population is also analyzed and a bound for the herd immunity threshold is calculated. The validity of the work done is verified through numerical simulations by applying the proposed model and strategy to analyze the multi-strains of the COVID-19 virus, in particular, the Delta and the Omicron variants, in the United State.

1 Introduction

The growing threat of infectious diseases with resistance to drugs and vaccinations, causing large number of deaths worldwide, is a cause for concern to the medical community and the general population. Scientists around the world are working to learn more about such diseases in order to study how likely an emerging variant of the disease can spread more easily than existing original variants. More data and analyses are needed for such study. These analyses will shed more light on the possibility of reinfections in people who already recovered from original strain, and infections in people who are fully vaccinated against original or previous strains. An example of such disease is an emerging virus called the corona virus 2019 (COVID-19) virus that has infected and killed millions around the world within a period of two years. The virus was caused by the virus species ‘severe acute respiratory syndrome corona virus’, named SARS-CoV-2. The airborne transmission occurs by inhaling droplets loaded with SARS-CoV-2 particles that are expelled by infectious people. Symptoms of the virus appear 2–14 days of exposure to certain strains of the virus. Several SARS-CoV-2 variants have emerged around the world, with each new variants having different characteristics. The COVID-19 virus evolves as changes in its genetic code occur during replication of the genome. The United States Centers for Disease Control and Prevention (CDC) and the US government SARS-COV-2 Inter-agency Group (SIG) [1] have confirmed the emergence of at least twelve new variants of the virus. The SIG group evaluates the risk of the circulating SARS-CoV-2 variants in the United States and make recommendations about the impact, severity, and how spread the virus is. Lineages of these variants that are classified as variant being monitored (VBM) and designated as Variant of Concerned (VOC) in the United States may lead to more severe cases of the virus. By lineages, we mean a group of related virus variant from a common ancestor. The variant, called Omicron-B.1.1.529 [2] was first detected in specimens collected on November 11, 2021 in Botswana and on November 14, 2021 in South Africa. It was first detected in the United States on December 1, 2021 and classified as a VOC. There are increased attention given to the Omicron variant. It is now the main circulating variant in the United States as of early January, 2022, accounting for about 95% of the US reported cases because of its high infection rate. Studies show that even vaccinated individuals can still be infected with this new variant if proper care is not taken. Scientist around the world are working to learn more about this new variant, and to gather more data for the purpose of studying if the variant causes more illness, hospitalizations, and death than infection with other variants. As of early December 2021, one of the new variants, called the Delta-B.1.617.2 variant was said to be the main circulating variant in the United States, and also classified as variant of concern by the SIG group. It was first identified in India and was reported to spread much faster and causes more severe cases than other early variants in the United States, probably causing twice [3] as many infections. Other variants, called the Alpha-B.1.1.7, Beta-B.1.351, Gamma-P.1, Epsilon B.1.427 & B.1.429, Eta B.1.525, Iota B.1.526, Kappa B.1.617.1, Zeta P.2, Mu B.1.621, are classified as VBM. The Alpha-B.1.1.7, Beta-B.1.351, and Gamma-P.1 variants were first discovered in the United Kingdom, South Africa, and Japan/Brazil, respectively. We aim to study how these variants are being transmitted and the impact that vaccines are having on mitigating the number of infection cases in a particular population.

Several mathematical models [4-36] have been developed to study the transmission of infectious diseases. Some of these works [5, 11, 20, 24, 31, 32, 37–43] discussed the transmission of infectious disease caused by the variants and lineages of the COVID-19 virus. SIR models including a modified SIR model with two strains and vaccinated group [10], a generalized SIR model with n- strains [39], a SIR model with complete cross-protection and nonlinear force of infection [44], a coupled multi-strain SIR epidemic model [24] have been developed to describe the transmission of the COVID-19 strains. Other works such as a multi-strain SEIR models with optimal control [5, 6], multi-strain SEIR models with saturated and general incidence rates [5, 15], SIRD [45] and SEIPAHRF model [46] with Caputo fractional derivative, SCIRP model incorporating media influence [47], and statistical analysis [48] have also been considered in describing the transmission of the virus and its strains. We direct the readers to the work of Hattaf et al. [49, 50] for more recent information about the fractional differential equations and its generalization.

Viruses undergo changes and these changes are cause for concern for people who have recovered from the virus, and also to those who are vaccinated against certain strain of the virus. As discussed in Fudolig et al. [10], a highly infectious emergent strain can infect the susceptible population before the original strain, thereby impeding the spread of the original strain or causing the two strains to coexist in an endemic equilibrium. For this reason, the need to determine conditions in which a newly emerged strain and an existing strains that have a means of immunity will coexist in a population is of utmost important.

Most of the papers mentioned above often utilize simpler versions of the multi-strain network, and provides less mathematical details on the asymptotic behavior of the model. In this work, we develop a multi-strain compartmental model by assuming a population is completely susceptible to n- different variants of a particular virus at the beginning of an epidemic, with the population of the region partitioned into compartments consisting of susceptible population, population vaccinated against strain k (1 ≤ k ≤ n) of the virus, population exposed to strain k of the virus, infected asymptomatic population with strain k of the virus, infected symptomatic population with strain k of the virus, and population that recovered from strain k of the infection. By denoting P = {1, 2, ⋯, n} and as a subset of the power set 2 with r number of strains, r = 0, 1, ⋯, n, with r = 0 representing disease-free case, we study the existence and stability conditions for the equilibrium point corresponding to the scenario where only strains in survive. This result is used to estimate the secondary number of infections produced by strain-k infected individual when introduced into a completely susceptible population. This number helps public health expert control the spread of the virus as new variants emerge. Conditions under which an endemic with more than one strain of the virus exist are calculated. The question as to whether this is possible is answered by studying the stability analysis of endemic strains. This model also helps us to understand the correlation between the daily number of administered vaccines and the number of infected, exposed, and recovered population in the region better. To understand the disease dynamics better, we study a powerful quantitative concept that can be used to characterize the contagiousness of each strain of the infectious disease and how transmissible they are. The basic reproduction number, which is the expected number of secondary cases produced by a typical infectious individual in a completely susceptible population in the presence and absence of vaccination are calculated. This study also helps to shed more light on the possibility of a disease being eliminated from a population if enough individuals are immune due to either vaccination or recovery from prior exposure to the disease. A bound for the herd immunity threshold, which is the minimum proportion of the population that must be vaccinated in order to stop the spreading of the disease in the population is calculated and analyzed for each variants.

This paper is organized as follows: In Section 2, a model is developed for the transmission of multi strain infectious diseases for the case where individuals vaccinated against specific strains are immune to that strain and its predecessors but can still be infected by newer emerging strains.The validity of the model, together with the existence and uniqueness of its solution is proved in this section. The reproduction numbers for the cases where the population is vaccinated and when it is not vaccinated are calculated. With these, the effect of vaccination in mitigating infection is studied. Existence of equilibrium points for the case where certain strains persist in the population is examined. In Section 3, the local and global stability of all equilibrium points for model (1) is investigated. In addition to the assumption made in Section 2, an extension of model (1) for the case where those that recovered from a particular strain can get infected by emerging strain is derived and analyzed in Section 4. Existence and stability of equilibrium points for the model is also discussed. Numerical simulations are carried out in Section 5 by applying models (1) and (35) to analyze COVID-19 data. Summary of the work done in this work is discussed in Section 6.

2 Materials and methods

2.1 Model formulation

By assuming a population is completely susceptible to n- different variants of a particular virus at the beginning of an epidemic, we partitioned the population of the region into compartments consisting of susceptible (denoted S(t)) population, population V(t) vaccinated against strain k of the virus, population E(t) of exposed individuals infected with strain/variant k of the virus, infected asymptomatic population A(t) with strain/variant k of the virus, infected symptomatic population I(t) with strain/variant k of the virus, and population R(t) that recovered from strain k of the virus at a given time t, for k = 1, 2, ⋯, n. These state variables are described in Table 1. We assume that the order of existence of newer strains follow k = 1, 2, ⋯, n. That is, the population is first infected by strain 1, followed by strain 2, ⋯, n. We assume those who are vaccinated and recovered from strain k are immune to previous and present strains j = 1, 2, ⋯, k, but can be infected by newer strains j = k + 1, k + 2, ⋯, n. We can see, for example, that this assumption is satisfied in the case of COVID-19 and supported by the United States Centers for Disease Control and Prevention (CDC), who claimed that current vaccines approved in the United States are expected to protect against severe illness, hospitalizations, and deaths caused by variants of the virus, although people who are fully vaccinated can still be infected. We assume that the migration rate into the susceptible population not receiving vaccination against any of the strains of the virus is (1 − q)μ, while the migration rate into the population vaccinated against strain k is μq, where . The parameter β denotes the per capita contact rate of susceptible individuals with symptomatic infected strain j; γ denotes the per capita contact rate of susceptible individuals with asymptomatic infected strain j; and denote the reduced per capita contact rates of symptomatic and asymptomatic strain j infected individual, respectively, that are vaccinated against strain k, j = k + 1, ⋯, n; h and ϵ denote the per capita contact rates of symptomatic and asymptomatic strain j infected individual, respectively, that recovered from strain k, j = k + 1, ⋯, n. These parameters are described in detail in Table 2.

Table 1

Description of variables for the epidemic model.

Variable	Description
S	Population of susceptible individuals
V k	Population of vaccinated individuals immune to strain k of the infection
E k	Population of individuals exposed to strain k of the infection
A k	Population of asymptomatic individuals infected with strain k of the infection
I k	Population of symptomatic individuals infected with strain k of the infection
R k	Population of individuals who recovered from the k-th strain

Table 2

Description of parameters for the epidemic model.

Parameter	Description
β k	Transmission rate of symptomatic infected individuals with strain k of infection interacting with susceptible population
γ k	Transmission rate of asymptomatic infected individuals with strain k of infection interacting with susceptible population
β¯_k,j	Transmission rate of symptomatic infected individuals with strain k interacting with population vaccinated against strain j of infection, k = j + 1, ⋯, n
γ¯_k,j	Transmission rate of asymptomatic infected individuals with strain k interacting with population vaccinated against strain j of infection, k = j + 1, ⋯, n
h k, j	Transmission rate of symptomatic infected individuals with strain k interacting with population that recovered from strain j of infection, k = j + 1, ⋯, n
ϵ k, j	Transmission rate of asymptomatic infected individuals with strain k interacting with population that recovered from strain j of infection, k = j + 1, ⋯, n
μ	Natural birth rate
p	Fraction of infection cases that are asymptomatic
q k	Fraction of population vaccinated against strain k of infection
λk	Transition rate of individuals with strain k of infection from exposed to infected class
r k	Asymptomatic recovery rate of those with strain k of infection
θ k	Symptomatic recovery rate of those with strain k of infection
q	∑j=1nq_j

In this section, we study a case where individuals vaccinated against specific strains are immune to that strain and its predecessors but can still be infected by newer emerging strains. Without loss of generality, we assume the population is normalized so that the sizes S, V, E, A, I, and R are in percentages, for k = 1, 2, ⋯, n. The model governing the transmission of the infectious disease for this case follows the system of deterministic differential equation where and the states and parameters in the model are described in Tables 1 and 2, respectively.. An extension of model (1) to include a case where those that recovered from a particular strain can get infected by emerging strain is discussed in Section 4.

In order to better understand the transmission dynamics described in (1), we give a schematic diagram of the model in Fig 1. The circle compartments represent group of individuals. An arrow pointing out of a compartment represents migration of individuals out of the compartment.

Fig 1

Schematic diagram for the epidemic model (1).

The circle compartments represent group of individuals.

Remark 1. Since individuals in the vaccinated group V are only assumed to be immune to strain j and its predecessors, and not to future strains k > j, we assume they are also susceptible to future strains k > j. For this reason, we expect that the rate at which an infectious individual with strain k make contact with susceptible individuals and individuals in group V, j < k, should be the same. That is, for j = 1, 2, ⋯, k − 1. Likewise, we expect β = h, , for j = 1, 2, ⋯, k − 1. Vaccination against past and current strains do not provide any protection against new emerging strains. In the case where the infectivity of strain k (k > j) is different for the vaccinated group V due to certain circumstances so that and , we present a comparison of the reproduction number for the case where , and the case where , for j = 1, 2, ⋯, k − 1 in Remark 3. Also, based on model (1), an individual infected with certain strain k cannot be infected with another strain j ≠ k at a given time t.

2.2 Validity of the epidemic model (1)

In this section, we discuss the validity of the proposed epidemic model (1), that is, we discuss the existence and uniqueness of the solution of (1). Let . Since the population is assumed to be normalized, we have denoting the entire population. It follows that dN/dt = μ − μN, N(t0) = 1, so that the population size is constant over time. This suggests that the birth and death rates (with the recruitment rate into the susceptible S and vaccinated V classes in the absence and presence of vaccines being (1 − q)μ and qμ, respectively, k = 1, 2, ⋯, n) are assumed to be the same so that the population is constant over certain period of time. The most appropriate epidemiological feasible region where solution exist is the set The existence, uniqueness, and positiveness of the solution of (1) is shown for the case where S0 > 0, V > 0, E > 0, A > 0, I > 0, R > 0 using results from Kelley and Peterson [51].

Theorem 1. If S0 > 0, V > 0, E > 0, A > 0, I > 0, R > 0, k = 1, 2, ⋯, n, then there exist a positive unique solution of (1) in the feasible region for all t ≥ 0.

Proof. Let S0 > 0, V > 0, E > 0, A > 0, I > 0, R > 0 for model (1). Define so that Eq (1) can be written as for some vector function h(t). It can be shown that is continuous with continuous first-order partial derivatives with respect to . It follows from Theorem 3.1 of Kelley and Peterson [51] that there exist a unique solution S(t), V(t), E(t), A(t), I(t), R(t), k = 1, 2, ⋯, n of (1) for all t ≥ 0. Let . It follows from (1) that S(t) satisfies for all t ≥ 0. Likewise, we can show in a similar manner that V(t) > 0, E > 0, A(t) > 0, I(t) > 0, and R(t) > 0 for all k = 1, 2, ⋯, n, t ≥ 0. The result follows since dN/dt = μ − μN, N(t0) = 1.

The long term behavior of the solutions of model (1) depends on certain thresholds of the reproduction number. The reproduction number and the thresholds are calculated in sections to come.

2.3 Reproduction number

Define

The disease-free equilibrium, denoted , of the system (1) is given by The dynamic model (1) can be written in the form where denotes the rate of appearance of new infections in compartment j, with and denoting the rate of transfer of individuals in and out of compartment j, respectively [52]. For any vector u = (u1⋯, u), define the n × n matrix (u) by the diagonal matrix Define where q0 = 0. Let and be two vectors with entries respectively, where and . Define F and V such that where (i, j) with 1 ≤ i, j ≤ 3n corresponds to the index of the infected compartments. Based on our model (1), there are three infected compartments E, A, and I, each assumed to have n-different strains so that the infected compartments are the first 3n entries of x in (8). Let 0 represents the zero square matrix of order n. The corresponding matrices F and V are calculated as where (u) is as defined in (9) and , , , , , and are vectors. From the matrix where ()−1 is the diagonal matrix , we see that the average length of time an individual spent being exposed, asymptomatically infected, symptomatically infected with strain k is 1/e, 1/a, and 1/w, respectively. Also, the average length of time an individual exposed to strain k spent being asymptomatic with the same strain during its life time is expected to be . The average length of time an individual exposed to strain k spent being symptomatic with strain k during its life time is expected to be . Since the transmission rates of symptomatic and asymptomatic individuals with strain k in a susceptible population are β and γ, respectively, it follows that at the emergence of a new strain k (k ≥ 2), the expected number of new infections produced by individual with such strain in a completely susceptible non-vaccinated population, and a population completely vaccinated against strain j, are and , respectively, j = 1, 2, ⋯, n, where we set q0 = 0.

The expected number, , of new infections produced by individual with such strain k in a population containing susceptible and vaccinated individuals is obtained as where 1 − q is the proportion of those that are susceptible but not vaccinated. Using the next generation matrix [52], since F is non-negative and V is a non-singular M-matrix [53], the number obtained in (11) is the (k, k)-th entry of the next generation matrix FV−1, for k = 1, 2, ⋯, n. It is the expected number of new infections in compartment k (compartment exposed to strain k) produced by infected individual originally introduced into the same compartment.

Remark 2. Effect of Vaccination

We remark here that the expected number of new infections in the compartment exposed to strain k depends on the vaccination rates q, l = k, k + 1, ⋯, n. If no one is vaccinated in the population (that is, if q = 0 for all l = 1, 2, ⋯, n), then the expected number of infections caused by strain k is obtained to be It follows from Remark 1 that this number is clearly more than the reproduction number obtained in (11) where some individuals are receiving vaccination in the population. The reproduction number can be written in terms as showing that the ratio of numbers of infection caused by strain k in a population with vaccination to a population without vaccination is 1 − q + Q : 1. If individuals are not vaccinated against strain k (that is, if q = 0 for fixed k), then the expected number of infections caused by strain k increases to which is the same as . This shows that a non vaccinated strain k infected individual produces more infections than a vaccinated strain k infected individual. If on the other hand, everyone in the population is vaccinated (that is, q = 1), then the expected number of infections caused by strain k significantly reduced to That is, a completely vaccinated population produces lesser infections than population not completely vaccinated. These analyses show the importance of being vaccinated in the population.

The reproduction number for the system, which is defined as the expected number of secondary infections produced by a typical infected individuals over the course of its infectious period, is calculated using the next generation matrix approach by van den Diessche et al. [52] as We shall later show that if , then we expect a typical infected individual to produce less than one new infected individual, meaning all strains of the disease will eventually die out in the population (in the presence of vaccination). Although it is well known that a value of the reproduction number greater than one means that epidemic will persist in the population [8, 52], we shall later show that a strain with a reproduction number greater than 1 can still die out on the long run if a newer emerging strain has a greater reproduction number. The value on the other hand, can be interpreted as the reproduction number for typical individual infected with strain k of the disease, for k = 1, 2, ⋯, n.

Remark 3. In the case where and such that and for all 1 ≤ k, j ≤ n due to some form of partial immunity as a result of vaccines, then the strain k reproduction number and the reproduction number in (11) and (14), respectively, are obtained as where and we set ∀ k = 1, 2, ⋯, n. It is easy to show that in this case, so that the estimates of the reproduction numbers given in (15) are greater than those given in (11) and (14), respectively. This shows that the reproduction number reduces as vaccines reduce the infectivity of the viruses (as expected).

Remark 4. The reproduction number is for the case where some members of the population are vaccinated against certain strain of the disease, that is, case with compartments V1, V2, ⋯, V in the population. As explained in Remark 2, the reproduction number, denoted , for a completely susceptible non-vaccinated population is given by

2.4 Existence of equilibrium points

We discuss conditions under which equilibrium points of (1) exist. System (1) has many equilibrium points. Let P = {1, 2, ⋯, n} denotes set of indices representing order of existence of new strain in the population, with 2 denoting the power set of P. Let denotes a subset of 2 with r number of strains, r = 0, 1, ⋯, n, with r = 0 representing disease-free case. We study the existence conditions for the equilibrium point corresponding to the scenario where only strains in survives, r = 0, 1, ⋯, n. The case r = 0 exists for the disease-free equilibrium case. The disease-free equilibrium is given in (6). For the case r = 1 representing case where only one strain survives, we shall denote the strain by strain m, m = 1, 2, ⋯, n, and the equilibrium referred to as the strain m equilibrium point and denoted . We study the conditions under which such equilibrium point exists, and in general, we also study conditions under which strains with indices in survives, for r = 1, 2, ⋯, n.

Theorem 2. The strain m unique equilibrium point, denoted , for the epidemic model (1) exists in the feasible region provided .

Proof. The strain m equilibrium point is obtained by solving the equation h(t) = 0 and setting E = A = I = R = 0 for 1 ≤ k ≠ m ≤ n, where h(t) is a vector in (4) containing the right hand side of (1). It follows from the equation governing A and I that and . Substituting these into the equations governing S, V, k = 1, 2, ⋯, n, E, and R, we have where c and are defined in (5). If , the strain m equilibrium point is obtained as where and for 1 ≤ k ≠ m ≤ n. We can show that by using (5) to show that

Remark 5. Theorem 2 shows that strain for 1 ≤ k ≠ m ≤ n (guaranteeing for 1 ≤ k ≠ m ≤ n) and . Theorem 2 is also valid for the case where and , j = 1, 2, ⋯, k − 1. The proof of this is shown in Theorems 17 and 18.

Suppose with strain 1 ≤ τ1 < τ2 ≤ n, , j = 1, ⋯, 2. We denote the equilibrium point corresponding to the case where strains τ1, τ2 survive by . We give theorem under which the equilibrium exists.

Theorem 3. The epidemic model (1) has an equilibrium point

Proof. Condition (20) implies that . The equilibrium point is obtained as where and if . If , then

The result follows.

Remark 6 Since (20) implies that τ1 before the emergence of strain τ2 caused an endemic. This is also confirmed in the work of Fudolig et al [10] where an (20) is equivalent to where and . We analyze the condition geometrically in Appendix B in S1 Appendix. It shows that for only strains τ1 and τ2 to remain in the system on the long run, the number of infection produced by strain τ2-infected individuals must be more than but not up to the number produced by strain τ1-infected individuals. Once falls outside this region, then only strain τ1 or τ2 remains in the system on the long run. For instance, if , it follows from (22) that the compartmental equilibrium value and , so that the values , , and are positive but , , and are negative, showing that only strain τ1 remains on the long run. In the same sense, we can show that only one strain remains in the system on the long run if . This endemic region is shown graphically in Appendix B in S1 Appendix.

Also, it can be shown that

Suppose with strain 1 ≤ τ < τ ≤ n, j = 2, ⋯, r, and . We denote the equilibrium point corresponding to the case where strains τ1, τ2, ⋯, τ survive by . We give theorem under which the equilibrium , r = 2, ⋯, n, exists.

Theorem 4. The epidemic model (1) has an equilibrium point

Proof. Condition (24) implies that . The equilibrium point is obtained as where and for 1 ≤ k ≠ τ1, τ2, ⋯, τ ≤ n if condition (24) is satisfied. If condition (24) is satisfied, then and for k = 2, 3, ⋯, r − 1, , and .

Remark 7. Condition (24) implies That is, the system is already in endemic state with strain τ before the emergence of strain τ, i = 1, 2, ⋯, r − 1, caused an endemic. Theorem 4 shows that if only strains τ1, τ2, ⋯, τ survive, then the system (without these strains) must converge to the disease-free equilibrium state.

2.4.1 Endemic equilibrium

The result for the endemic equilibrium can be calculated from Theorem 4 by extending the set to . We state the result without proof in the next theorem.

Theorem 5. The epidemic model (1) has an endemic equilibrium point in the feasible region provided

In the next section, we discuss the convergence of the system in (4). We study conditions under which all strain infections are eradicated in the population on the long run. We also discuss the condition under which certain strain of the disease persists.

3 Stability analysis

In this section, we discuss the convergence of system (1) under certain conditions. That is, we study condition(s) under which the system converges to disease-free or equilibriums , r = 1, 2, ⋯, n.

3.1 Stability analysis of the disease-free equilibrium point

Theorem 6. The disease-free equilibrium if .

Given the initial condition , Theorem 6 shows that if the system satisfying (1) starts near the initial point , then the system converges to the equilibrium point (that is, all the strain of infections are eradicated) provided the threshold . We use the idea presented in [27] to prove Theorem 6.

Proof. Let be the n × n identity matrix, , = {β, ⋯, βn}, γ = {γ, ⋯, γn}, the matrices and defined by so that where is defined in (9). The linearization of the model (1) at the equilibrium point is derived as where In order to show that the equilibrium point is locally stable, we need to show that the maximum real part, s(A), of A is negative. From the structure of the matrix A, this reduces to showing that the maximum real part of the eigenvalues of matrix is negative (or equivalently, ). The determinant of the matrix is If , then . Matrix is a non-singular Z matrix that can be written in the form where and are lower and upper diagonal matrices, respectively, with positive diagonals obtained as where , and for j = 1, 2, …, 3n, with The diagonals and if . Therefore, the matrix is a non-singular M matrix, and hence, a P-matrix. It follows from Berman [54] and Plemmons [53] that the maximum real part of the eigenvalues of is negative.

Remark 8. Theorem 6 shows that if the population starts sufficiently close to − q, the proportion of the size of the vaccinated class immune to strain

We prove in the next theorem that the population, irrespective of where it starts from, converges to the point if .

Theorem 7. The disease-free equilibrium is globally stable in the feasible region if .

Proof. Define the Lyapunov function by where where . Let s = S/S0, . If , it follows that where and using (6). Using the fact that the arithmetic mean of a list of non-negative real numbers is greater than or equal to the geometric mean of the same list [55], it follows that if , then . Also, if S = S0 = 1 − q, I = E = R = 0, (or equivalently, if s = 1, I = E = R = 0, v = 1) for all k = 1, 2, ⋯, n. Since is the largest invariant set in the subset of where , its global stability follows by the LaSalle’s Invariance Principle [56]

Remark 9. Theorem 7 shows that the susceptible class converges to a fraction 1 − q of the entire population size (this is simply the population size without the vaccinated group), the vaccinated class immune to strain k converges to a fraction qk of population size, and no infected class, hence no need of recovery on the long run if . This suggests the threshold for disease eradication is the number .

3.2 Bound for the critical vaccination threshold

Herd immunity is a state where significant proportion of the population is immune to an infection so that only few susceptible individuals can be infected and transmit the infection [57]. Classical vaccine-induced herd-immunity threshold suggests that the spread of a disease can be stopped by vaccinating certain fraction of the population. This might be invalid and biased due to many factors such as emergence of multi variant strains of the virus/disease, the dynamic nature of virus transmission, presence of immunity due to infection, changes in implementation and adherence to public health measures, and uncertainties in vaccine effectiveness and duration of immunity [57, 41]. These can cause the estimation of the Herd-immunity threshold to be imprecise. In this section, we aim to estimate a bound for the minimum proportion of the population that must be vaccinated in order for certain infection to die out in the population.

Let E ∈ (0, 1] denotes the proportion of vaccinated individuals against strain j who are protected by vaccines, for j = 1, 2, ⋯, n. It follows from Remark 4 that we can write in (14) in terms of in (17) as follows: for some constant c ∈ (0, 1). If the vaccine is perfect and provides 100% immunity, then the vaccine effectiveness E = 1, otherwise, E < 1 if it provides only partial immunity. According to Theorem 7, the population reaches herd immunity with respect to strain j, with incidence of the infection declining, if . This yield The number is the minimum proportion of a population infected that must be vaccinated to stop strain j from spreading. This number is referred to as the herd immunity threshold [41]. In the presence of multiple vaccines for strains j = 1, 2, ⋯, n with effectiveness E1, E2, ⋯, E, respectively, let H and H denote the minimum and maximum herd immunity thresholds for the strains j = 1, 2, ⋯, n, respectively. The interval contains the maximum proportion of the population expected to be vaccinated to stop the spread of the disease.

Remark 10. Caveats to the estimate

We note here that the bound in (33) is calculated for a population satisfying the dynamics given in (1). The bound is expected to change in the emergence of a new strain of the virus. Also, the bound is calculated without taking into consideration the proportion of those who have immunity from the virus. These and more are some of the caveats to this estimates.

Sometimes new variants emerge and disappear. Other times, new variants persist. We study conditions under which new variant, say strain m, persists in the population. In the next section, we study the behavior of the system in the case where disease is not completely eradicated in the system.

3.3 Stability analysis of strain m equilibrium

In the next theorem, we study how the population behaves on the long run if the number, , of new infections produced by infected individual with strain m is greater than one, while for 1 ≤ k ≠ m ≤ n.

Theorem 8. The strain m equilibrium for the epidemic model (1) is globally stable in the feasible region if for all 1 ≤ k ≠ m ≤ n and .

Proof. According to Theorem 2, the strain m equilibrium exists and non-negative if . Define the Lyapunov function by where with φ0 = 0. Let s = S/S*, , , , . It follows that if and for 1 ≤ k ≠ m ≤ n, where and using (19). Using the fact that the arithmetic mean of a list of non-negative real numbers is greater than or equal to the geometric mean of the same list [56], it follows that . Equality holds if S = S*, I = E = R = 0 for all k ≠ m, for all k = 1, 2, ⋯, n, . Since is the largest invariant set in the subset of where , its global stability follows by the LaSalle’s Invariance Principle [56].

Remark 11. Theorem 8 shows that the system converges to on the long run if for 1 ≤ k ≠ m ≤ n and irrespective of the starting point. That is, if the number, , of new infections is greater than one while for 1 ≤ k ≠ m ≤ n, then on the long run, the susceptible class converges to a fraction of the entire population size, no exposure to any strain other than strain m in the population (hence no infection other than those caused by strain m, and no need for recovery), and before the existence of strain m (that is, k < m), the vaccinated class immune to strain k < m converges to a fraction of the population size. We see from Theorem 7 that if there is no endemic in the population, the susceptible population S converges on the long run to (1 − q) × 100% of the population. This number reduces by in the presence of an emerging strain m. Also, from Theorem 7, the vaccinated class Vk converges to q × 100% of the population size if there is no endemic in the population. This number also reduces by in the presence of an emerging strain m. If k > m on the other hand, since for 1 ≤ k ≠ m ≤ n, we have the vaccinated class immune to strain k (k > m) converging to a fraction of the population, no exposure to any strain other than strain m in the population (hence no infection and no need for recovery).

3.4 Stability analysis of the equilibrium point

We give the proof of the stability of the equilibrium point in the next theorem.

Theorem 9. The equilibrium point is globally stable in the feasible region if for 1 ≤ k ≠ τ, τ ≤ n and condition (20) is satisfied.

Proof. Assume for all k ∉ τ, τ and condition (20) is satisfied. The existence of the equilibrium point follows from Theorem 3. Define the Lyapunov function by where with φ0 = 0. The derivative of computed along the solution of (1) is It follows from (34) and Remark 6 that and Let s = S/S+, , , , . Define We have where using (22). It follows that . Equality holds if S = S+, I = E = R = 0 for all , for all k = 1, 2, ⋯, n, . Since is the largest invariant set in the subset of where , its global stability follows by the LaSalle’s Invariance Principle [56].

3.5 Stability analysis of the equilibrium point

Theorem 10. For r = 3, 4, ⋯, n, the equilibrium point is globally stable in the feasible region if for all and condition (24) is satisfied.

Proof. The proof of Theorem 10 is similar to that of Theorem 9 by extending to , r = 3, 4, ⋯, n.

Remark 12. Theorem 10 can be extended to a case where .

4 Model for re-infected recovered and vaccinated population

There have been confirmed cases of the COVID-19 reinfections around the world [38, 40, 42]. In this section, in addition to assuming that individuals vaccinated against strain k can gets infected with emerging strains j > k, we also discuss the case where individuals who recovered from strain k can be infected with emerging strains j > k. For this additional assumption, we extend (1) to the form The schematic diagram of model (35) is given in Fig 2.

Fig 2

Schematic diagram for the epidemic model (35).

The circle compartments represent group of individuals.

4.1 Validity of the epidemic model (35)

In this section, we discuss the validity of the proposed epidemic model (35).

Theorem 11. If S0 > 0, V > 0, E > 0, A > 0, I > 0, R > 0, then there exist a positive unique solution of (35) in the feasible region for all t ≥ 0.

Proof. The proof is similar to Theorem 1.

4.2 Existence of equilibrium points for model (35)

The disease-free equilibrium of the system (35) is the same as that of the system (1), and given by

4.3 Reproduction number for model (35)

It can also be shown in a similar manner that model (35) has the same reproduction numbers and as that of model (1). This shows that the number of infection in a population where every individual who recovered from strain k is immune to all possible strains is the same for the population where individuals who recovered from strain k can still be infected with emerging strains j > k. Hence, for strain k infected individual, the expected number, of new infections produced by individual with strain k in a susceptible population satisfying (35) is the same as (11). Likewise, the reproduction number for the system (35) is obtained in a similar manner to (14).

4.4 Existence of endemic equilibrium points for model (35)

Theorem 12. The strain (35) exists in the feasible region (1) and (35) and obtained in Theorem (2).

Remark 13. Theorem 12 shows that the existence of a single strain endemic (strain m in this case) does not depend on whether the recovered class is immune to the strain or not. Regardless of the immunity status of the recovered strain k-class to the strain, the equilibrium point of the system will be if .

Define It can be shown that . Conditions for existence of other equilibrium points are given in Theorem 13.

Theorem 13. The epidemic model (35) has an equilibrium point satisfying for k ≠ τ1, τ2, in the feasible region provided is satisfied.

Proof. The proof is similar to that of Theorem 4.

Remark 14 Unlike Remark 13, we see here that the existence of two endemic equilibrium points depends on the immunity status of the recovered class. The exposed and infectious equilibrium points , , and for model (35) are greater compared to that of model (1). This is because unlike model (1), model (35) suggests that the recovered population R is not fully immune to emerging strains j > k, increasing the possibility of populating the exposed and infected classes. In addition, condition (39) implies that . The condition is equivalent to , where and .

4.5 Global analysis of equilibrium points for model (35)

Theorem 14 The disease-free equilibrium (35) is locally asymptotically stable in the feasible region

Proof. The proof is similar to that given in Theorems 6 and 7.

Theorem 15. The strain m equilibrium for the epidemic model (35) is globally stable in the feasible region if for all 1 ≤ k ≠ m ≤ n and .

Proof. The proof is similar to the proof of Theorem 8.

Theorem 16. The equilibrium point for model (35) is globally stable in the feasible region if for 1 ≤ k ≠ τ1, τ2 ≤ n and condition (39) is satisfied.

Proof. The proof is similar to the proof of Theorem 9.

5 Results: Covid-19 data analysis

We apply models (1) and (35) to analyze the United States daily COVID-19 cases (number of infection cases, recovery cases, and vaccination). The daily COVID-19 cases data are available on the CDC website [58] for the COVID-19 periods 04/01/2020 till present. We analyze the confirmed COVID-19 infection cases, and vaccination cases as reported by U.S. states, U.S. territories, New York City, and the District of Columbia from the previous day. Since the two recent variant of concerns (VOC) in the United States are the Delta and Omicron variants, we consider the case where n = 2 using model (35), with the Delta variant as Strain τ1 = 1 and the Omicron variant as strain τ2 = 2. Two analyses are performed in this section. The first analysis is shown in Section 5.1 to confirm the validity of the results derived in this work. Using published and estimated COVID-19 parameters, we confirm the stability results for the disease-free equilibrium, strain 1 equilibrium, strain 2 equilibrium, and the endemic equilibrium . In section 5.2, the real COVID-19 cases for the United States is analyzed using model (35).

5.1 Simulation results using published and estimated parameters

Using model (35), we confirm the existence, and stability of the disease-free equilibrium, strain 1 equilibrium, strain 2 equilibrium, and the endemic equilibrium for the case where strains 1 and 2 represent the Delta and Omicron variants, respectively. The CDC data for vaccination shows that about 63% of the population of the United States are fully vaccinated (either taken the two dozes of Pfizer or Moderna, or the single dose of Jannsen) as at January 25, 2022. For this reason, we chose q1 and q2 to be in the interval [0, 0.63]. In their paper, Saldana [59] shows that about 80% of infection was asymptomatic with average incubation and recovery time of 5 and 10 days, respectively. Also, Hay et al. [37] shows that the mean duration of Delta and Omicron’s infections is 10.9 days and 9.87 days, with 95% confidence intervals (8.83, 10.9) and (9.41, 12.4), respectively. For this reason, we set r1, θ1 ∈ [1/10.9, 1/8.83] and r2, θ2 ∈ [1/12.4, 1/9.41]. Based on the five COVID-19 Pandemic Planning Scenarios estimated by CDC, the number of infections that are asymptomatic is uncertain and in the interval [0.15, 0.7], with the best estimate of 30%. We set μ = 0.0124. Bernal et al. [60] shows in their studies that the effectiveness of two doses of BNT162b2 (Pfizer) vaccines was 88.0% (95% CI, 85.3 to 90.1) among those with the delta variant. On November 16, 2020, the company Moderna announced that their vaccine is more than 94% effective at preventing COVID-19, based on an analysis of 95 cases [61]. We use these estimates for the Herd immunity plot. Following results from Jing et al, we select values for the incubation period to be between 2 days and 7 days, so that λ1, λ2 ∈ [1/7, 1/2]. The range of parameters used in the simulation is shown in Table 3.

Table 3

Model parameter values.

Parameter	Units	Range	References
β 1	day−1	[0.1, 1.4007]	[37, 62–66]
β 2	day−1	[0.1, 1.6761]	Assumed
γ 1	day−1	[0.07, 0.9567]	[62, 63, 65, 66]
γ 2	day−1	[0.07, 0.9567]	Assumed
q 1	day−1	[0, 0.63]	CDC
q 2	day−1	[0, 0.63]	CDC
p	day−1	[0.4, 0.6]	CDC
λ1	day−1	[1/7, 1/2]	[67]
λ2	day−1	[1/7, 1/2]	[67]
r 1	day−1	[1/10.9, 1/8.83]	[37]
r 2	day−1	[1/12.4, 1/9.41]	[37]
θ 1	day−1	[1/10.9, 1/8.83]	[37]
θ 2	day−1	[1/12.4, 1/9.41]	[37]

Simulation result for the case where the population is free of disease on the long run is shown in Fig 3.

Fig 3

Stability analysis for the disease-free equilibrium.

Case where and . Here, we set β1 = 0.07, β2 = 0.2, γ1 = 0.1, γ2 = 0.12, μ = 1/80.3, q1 = 0.5, q2 = 0.4, λ1 = 1/5, λ2 = 1/4, r1 = 1/10.5, r2 = 1/9, θ1 = 0.09, θ2 = 0.1. The reproduction numbers and so that . On the long run, the susceptible population reduces to 10% of the population size while the population receiving vaccination against strains 1 and 2 increases to 50% and 40% of the population size, respectively, in this scenario. The exposed and infected population size converges to zero on the long run. The disease-free equilibrium is calculated as .

Fig 4 shows simulation result for the case where only strain 1 endemic exists in the population. Similarly, Fig 5 shows simulation result where only strain 2 endemic exists in the population. In Fig 6, we answer the question as to whether it is possible to have an endemic with more than one strain of the virus. The figure shows that this is possible if condition (39) is satisfied. Fig 7 shows the importance of condition (39) by confirming that no two strains remain in the population on the long run even if but . Fig 8 shows the simulation result for the case where .

Fig 4

Stability analysis for the strain 1 equilibrium.

Case where and . This is a case where β1 = 1.2, β2 = 0.4, γ1 = 0.5, γ2 = 0.08, μ = 1/80.3, q1 = 0.1, q2 = 0.5, λ1 = 1/5, λ2 = 1/4, r1 = 1/10.5, r2 = 1/9, θ1 = 0.09, θ2 = 0.1. A scenario where the population of those receiving vaccination against strain 1 (with high transmission rate) dropped, leading to strain 1 endemic. In this case, and so that . The strain 1 equilibrium is obtained as

Fig 5

Stability analysis for strain 2 equilibrium.

Case where and . In this case, we set β1 = 0.18, β2 = 0.9, γ1 = 0.1, γ2 = 0.4, μ = 1/80.3, q1 = 0.5, q2 = 0.1, λ1 = 1/5, λ2 = 1/4, r1 = 1/10.5, r2 = 1/9, θ1 = 0.09, θ2 = 0.1. In this case, and so that . The strain 2 equilibrium is obtained as . This vector shows proportions that each of the population sizes S, V1, V2, E1, E2, A1, A2, I1, I2, R1, R2 converge to on the long run.

Fig 6

Stability analysis for the equilibrium .

Case where . In this case, we set β1 = 1.2, β2 = 0.9, γ1 = 0.5, γ2 = 0.4, μ = 1/80.3, q1 = 0.1, q2 = 0.5, λ1 = 1/5, λ2 = 1/4, r1 = 1/10.5, r2 = 1/9, θ1 = 0.09, θ2 = 0.1. In this case, and so that . We see in this case that and , implying the existence of endemic with more than one strain. We see an endemic with both strains 1 and 2 because the population is already in an endemic state with strain 1 before strain 2 caused an endemic, and the number of secondary infection caused by strain 2 is more than but not up to that caused by strain 1. The endemic equilibrium in this case is

Fig 7

Stability analysis for the case where but condition (39) is not satisfied.

In this case, we set β1 = 1, β2 = 0.05, γ1 = 0.5, γ2 = 0.4, μ = 1/80.3, q1 = 0.1, q2 = 0.5, p = 0.6, λ1 = 0.2, λ2 = 0.25, r1 = 1/10.5, r2 = 1/9, θ1 = 0.09, θ2 = 0.1. The endemic equilibrium in this case is the strain 1 equilibrium . Here, , , and condition (39) is not satisfied since . This condition implies, from (38), that and , so that the values , , and are positive but , , and are negative. This shows that only strain 1 endemic exists on the long run. A similar analysis is presented in Appendix B in S1 Appendix geometrically.

Fig 8

Stability analysis for the case where .

In this case, we set β1 = 0.9, β2 = 1.0, γ1 = 0.4, γ2 = 0.45, μ = 1/80.3, q1 = 0.12, q2 = 0.1, p = 0.8, λ1 = 0.21, λ2 = 0.2, r1 = 1/9, r2 = 1/10.5, θ1 = 0.1, θ2 = 0.09. Here, and so that . The endemic equilibrium in this case is {0.1739, 0.0267, 0.1000, 0, 0.0410, 0, 0.0609, 0, 0.0160, 0, 0.5815}.

5.1.1 What happens if but condition (39) is not satisfied?

We study a case where but condition (39) is not satisfied. Although the condition is similar to that in Fig 6, we see here that the system converges to the strain 1 equilibrium point. That is, even though the reproduction numbers and are more than one, strain 2 still gets eradicated from the system on the long run while strain 1 caused an endemic. This study shows that the second inequality in condition (39) is necessary for the existence of endemic with more than one strain.

5.1.2 What happens if

We study an interesting case where although the reproduction number for strain 1 is greater than 1, the strain still gets eradicated from the system on the long run. This happens because a newer strain 2 caused an endemic with a higher reproduction number . This case is analyzed geometrically in Appendix B in S1 Appendix.

5.2 Real data: Analysis using the Delta variant cases for the period 10.09.2021 to 12.18.2021

The work done here is applied to real data using the overall SARS-CoV-2 weekly Variant Proportions [68] for the United States collected from the Centers for Disease Control and Prevention (CDC) for the period 10.09.2021 to 12.18.2021. The proportion of the Delta (B.1.617.2) variant is collected and used together with the number of daily cases (collected from CDC [58]) of all variants in the United States. Within this period, data shows that on the average, the proportion of the Delta variant is about 95% of the total cases in the United State, suggesting that the Delta variant is the dominating variant during the period. For this reason, we set n = 1 in (35) and estimated the parameters in Table 4 using the Nelder-Mead simplex algorithm as described in Lagarias et al. [69]. These parameters are estimated by minimizing the sum of square error between the real Delta infection cases and the simulated cases, where the real Delta infection case is generated by multiplying the weekly proportion of the Delta variant by the weekly COVID-19 cases (generated from the daily cases [58]). We use N = 331893745 for the total population of the United States, with initial point , and endemic equilibrium point . The real and estimated weekly COVID-19 cases for the United States for period 10.09.2021 to 12.18.2021 are given in Fig 9 and generated using model (35).

Table 4

Parameter estimates for model (35) with n = 1.

Parameter	β^_1	γ^_1	μ^	q^_1	p^	λ^_1	r^_1	θ^_1
Estimate	1.3800	0.6000	0.2045	0.01037	0.3	0.1500	0.1500	0.0300

Fig 9

Real (blue) and estimated (red) COVID-19 weekly cases for the Delta variant in the United States.

Fig 10 shows the estimated results for the population sizes S, V1, E1, A1, I1, and R1 for the Delta variant.

Fig 10

Estimated trajectory path for S, V1, E1, A1, I1, and R1 for the Delta variant.

The analysis suggests that the trajectories of the population of those receiving vaccination against the delta variant and those recovering from the variant are decreasing, while the susceptible, exposed and asymptomatic populations are increasing in the time period 10.09.2021 to 12.18.2021. The symptomatic population decreases from 10.09.2021 to 10.23.2021, after which it started increasing until the end of the analysis in 12.18.2021. The reproduction number for the variant, using (14), was calculated to be 1.94.

Using 88% and 70% for the vaccine effectiveness of two doses of Pfizer vaccine among those with Delta and Omicron variants, we calculate the Herd Immunity threshold bound to be [97%, 1). We give a plot of the Herd immunity as a function of the measure of the effectiveness of available vaccines in Fig 11.

Fig 11

Herd immunity as a function of the measure of the effectiveness of vaccines for the Delta variant in the United States.

5.3 Real data: Analysis using the Delta and Omicron variant cases for the period 12.11.2021 to 01.15.2022

The World Health Organization (WHO) and CDC classified a new variant, B.1.1.529, as a VOC and named it Omicron. The work done here is applied to real data using the SARS-CoV-2 weekly Variant Proportions [68] for the United States collected by CDC. The proportion of the Delta (B.1.617.2) and the Omicron (B.1.1.529) variants shown in Fig 12 are collected and used together with the number of weekly cases (collected from CDC [58]) of each variants in the United States. Data suggests the existence of the two variants as VOC and that the infectivity of the Delta variant is slowing down while the infectivity of the Omicron variant is increasing during the period 12.11.2021 to 01.15.2022. For this reason, we set n = 2 in (35) and estimated the parameters , , , , , , , , , , , , , and in Table 5 using the Nelder-Mead simplex algorithm [69]. The Delta and Omicron infected data are generated by multiplying the weekly proportion of the Delta variant by the weekly COVID-19 cases. We use N = 331893745 for the total population of the United States, with initial point . The real and estimated weekly COVID-19 cases for the infected cases I1 and I2 for period 12.11.2021 to 01.15.2022 are given in Fig 13 and generated using model (35). The reproduction number for the delta variant is less than one. On the other hand, the infectivity of the Omicron variant in this period was so high that the reproduction number obtained was . The strain 2 endemic equilibrium point suggests that the Delta variant’s infectivity and exposure will die out on the long run, while there is an Omicron variant endemic. The proportion of the population vaccinated against the delta variant converges to 0.18% on the long run, and that of the Omicron variant converges to 15% on the long run. The proportion of the susceptible, exposed, asymptomatic infected, symptomatic infected, and recovered Omicron population converge to 14.1%, 28.64%, 8.30%, 12.74%, and 21.04% of the population size, respectively, on the long run.

Fig 12

Proportion of the Delta (a) and the Omicron (b) variant for the period 12.11.2021 to 01.15.2022 in the United States collected from CDC2.

Table 5

Parameter estimates for model (35) for the Delta and Omicron variants with n = 2.

Parameter	β^_1	β^_2	γ^_1	γ^_2	μ^	q^_1	q^_2	p^	λ^_1	λ^_2	r^_1	r^_2	θ^_1	θ^_2
Estimate	0.0397	2.1614	0.2157	5.6282	0.1500	0.0105	0.1500	0.3946	0.2550	0.2204	0.9987	0.1500	0.2972	0.1500

Fig 13

Real (dotted blue) and estimated (red) COVID-19 weekly cases for the Delta variant and the Omicron variant for the period 11.27.2021 to 01.15.2022 in the United States.

The estimated results for all compartments are shown in Fig 14.

Fig 14

Estimated results for the Delta and Omicron variants.

The result shows that the number of Delta exposure and infection cases is decreasing (slowing down) while the number of Omicron exposure and infection cases is increasing (speeding up) between December 11, 2021 and January 15, 2022. These bring about a decrease in the population of those who are vaccinated against the Delta variant and an increase (at a slowing pace) in the number of those vaccinated against the Omicron variant. The number of those susceptible to the variants rises between the period 12.11.2021 to 01.01.2022, after which it started falling until 01.15.2022. The population of those who recovered from the Delta variant is speedingly descreasing while the count of those who recovered from Omicron variant is increasing at a speeding pace. The result of the trajectory of the estimated Delta and Omicron cases is similar to the plot shown in the work of Nyberg et al. [70].

5.3.1 Sensitivity analysis

In this section, we study the impact of each epidemiological parameters on the disease transmission and prevalence. We are interested in discovering the parameters that have high impact on the basic reproduction number , k = 1, 2, ⋯, n. This can be achieved using sensitivity analysis by calculating the sensitivity index of each parameters on . The normalized forward sensitivity index of a variable F that depends differentiably on a parameter u is defined as [71, 72] The analytic expression for the sensitivity index of , k = 1, 2, ⋯, n, with respect to the parameters in (35) is obtained as The positive sensitivity index with respect to β, γ, and λ shows that an increase in the value of the strain k’s symptomatic, asymptomatic transmission rates, and the latency rate leads to an increase in the basic reproduction number, , as expected. Likewise, the negative sensitivity index for μ, q, l = k, k + 1, ⋯, n, r, and θ shows that an increase in the natural death rate, vaccination rate, symptomatic and asymptomatic recovery rates, respectively, leads to a decrease in , as expected. Also, the zero value for the sensitivity index for q, l = 1, 2, ⋯, k − 1, shows that being vaccinated against predecessor strains l = 1, 2, ⋯, k − 1 does not have any impact on the current and future strains j ≥ k infection count. The sensitivity index is positive if . That is, an increase in the fraction of infection cases that are asymptomatic will lead to an increase in the basic reproduction number, , if the total number of infection caused by an asymptomatic strain k infectious individual is more than that caused by a symptomatic strain k individual. The magnitude of the sensitivity of to changes in these parameters can be calculated by evaluating for each parameter estimates u in Tables 4 and 5 for case where n = 1 and n = 2, respectively. For a particular parameter u, a higher value of suggests that is more sensitive to u. The sensitivity index of to each parameters estimated in Table 5 for model (35) is given in Table 6 for k = 1, 2.

Table 6

Sensitivity index of to parameters generated in Table 5 for model (35).

	Parameters u
	ΔuR_k	β k	γ k	μ	q k	p k	λk	r k	θ k
k = 1	ΔuR_1	0.4204	0.5796	−3.9138	−1.1912	0.7744	1.4524	−0.5046	−0.9401
k = 2	ΔuR_2	0.3707	0.6293	−5.0312	−1.1765	0.9823	1.8374	−2.0975	−1.2358

Analysis suggests that changes in the asymptomatic transmission rate contribute more to the changes in the reproduction number than the symptomatic transmission rate for the Delta (k = 1) and Omicron (k = 2) variants. Also, is more sensitive to changes in the infection rate λ than the transmission rates and the fraction p of cases that are asymptomatic. For the Delta variant, the reproduction number is more sensitive to the recovery rate θ1 of the symptomatic class than r1. Likewise for the Omicron variant, the asymptomatic recovery rate r2 has much influence over than the symptomatic recovery rate θ2.

6 Summary and discussion

In this work, a multi-variant epidemic model for analyzing the emergence and dissemination of viral multi-strains of an infectious disease in a population that is assumed to be completely susceptible to n different strains of the disease is developed and analyzed. A major assumption on the viral strains is that those who are vaccinated and recovered from a specific strain k ≤ n are immune to that strain and its predecessors but can still be infected by newer emerging strains. A study of how well-poised the model is is carried out by showing the existence of non-negative solution of the derived model. The model compares the cases where the force of infection on the susceptible vaccinated and unvaccinated populations are different and the same. The reproduction number for each specific strains j = 1, 2, ⋯, n is obtained and analyzed. We show that the reproduction number for the system where individuals who recovered from certain strain can be infected by emerging strain is the same for the system where such individuals cannot be re-infected by emerging strain. In order to shed more light on the possibility of an endemic with more than one strain of the virus, global stability analysis is obtained for different equilibrium points , r = 0, 1, 2, ⋯, n, of the system, with denoting the disease-free equilibrium. It was shown with respect to model (35) that for an endemic with strains τ1 ≤ n and τ2 ≤ n (where τ1 < τ2) to occur, it must be that the reproduction number . This condition shows that for an endemic with strains τ1 ≤ n and τ2 ≤ n to occur, the population must have been in an endemic state with the first emerged strain τ1 and the number of secondary infection caused by strain τ1 must be greater than that of strain τ2, with a necessary condition that . The importance of this necessary condition is emphasized numerically by showing that a system that only satisfies the condition does not converge to the endemic equilibrium with strains τ1 and τ2, but instead converges to strain τ1 equilibrium . The later condition reduces to for the case where individuals who recovered from certain strain cannot be infected by emerging strain. A general condition is obtained for the case where endemic with strains τ1, τ2, ⋯, τ (with τ1 < τ2 < ⋯ < τ) exists in the population. Our analysis shows that for an endemic with strains τ1, τ2, ⋯, τ to exist, the population must first be in endemic with strain τ1, followed by τ2, ⋯, τ, so that condition (24) satisfied. Also, we showed numerically that a strain with a reproduction number greater than 1 can still die out on the long run if a newer emerging strain has a greater reproduction number. The effect of vaccines on the population is also analyzed. The herd immunity for each strain j = 1, 2, ⋯, n is computed as a function of the effectiveness of vaccines against the strain. A threshold for the herd immunity for the case where there is endemic of more than one strain of the virus is also calculated. The result is applied to analyze real COVID-19 data for the Delta and Omicron variants in the United States. The reproduction numbers for the Delta and Omicron variants cases for the period 12.11.2021 to 01.15.2022 suggest that the Delta variant cases is dying down, while there is an endemic of the Omicron variant. Using model (35), possible trajectories for the susceptible, vaccinated, exposed, infected and recovered population are plotted by first estimating the parameters in the model. The parameters are estimated by minimizing the sum of square error between the real and estimated infection cases using the Nelder-Mead simplex algorithm [69]. Further research on the COVID-19 cases is ongoing and the results of the research will be made known once it is available. In the future, we plan to modify the model to fit emerging characteristics of the virus and extends limitations in the model to a more general case.

(PDF)

Click here for additional data file.