How Does Vaccine-Induced Immunity Compare to Infection-Acquired Immunity in the Dynamics of COVID-19?

TL;DR

This study compares how immunity from vaccines and infections affects the spread of COVID-19, finding that infection-acquired immunity has a stronger impact on disease dynamics.

Contribution

The study introduces a compartmental model showing that waning infection-acquired immunity significantly influences disease spread more than vaccine-induced immunity.

Findings

A backward bifurcation may occur when the vaccinated reproduction number (Rv) is equal to unity.

Global sensitivity analysis shows that the waning of infection-acquired immunity (ω2) has a higher influence on cumulative infections than vaccine-induced immunity (ω1).

Simulations indicate that changes in ω2 lead to more pronounced disease dynamics compared to changes in ω1.

Abstract

Five years into the COVID-19 pandemic, the availability of effective vaccines has substantially reduced new cases, hospitalizations, and mortality. However, the waning of immunity has been a topic of particular interest in relation to disease control. The objective of this study is to investigate the impact of the decline in vaccine-induced immunity (ω1) and infection-acquired immunity (ω2) on disease dynamics. For this purpose, we use a compartmental model with seven compartments that accounts for differential morbidity, vaccination, and waning immunity. A compartmental model divides a population into distinct groups depending on their disease status. The temporal changes in the compartments are represented through ordinary differential equations (ODEs). The model is mathematically analyzed to show that a backward bifurcation (i.e., a perverse outcome) may occur when the vaccinated reproduction number (Rv) is equal to unity. Both local and global sensitivity analysis on the reproduction number reveal that the vaccine efficacy, waning of vaccine-induced immunity, vaccine coverage rate, coefficients of transmissibility, and the recovery rate for mild infections are the most sensitive parameters. The global sensitivity analysis on the cumulative number of infections shows that ω1 and ω2 are both pivotal parameters, while ω2 has a higher influence. Simulations on infections and mortality suggest that the changes in ω2 result in dynamics that are more pronounced compared to the dynamics resulting from the changes in ω1, thus indicating the importance of the duration of infection-acquired immunity in disease spread.

Figures

Full text

1. Introduction

COVID-19 is a highly contagious viral disease caused by a novel type of coronavirus known as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1]. COVID-19 is mainly transmitted by the breathing of air carrying aerosol particles containing the SARS-CoV-2 virus or by the eyes, nose, or mouth having contact with such particles [1]. Symptoms usually appear within two weeks after exposure to the virus, and can include a variety of conditions ranging from mild symptoms to severe illness and Long COVID (or Post-COVID conditions) [2]. After its first detection in Wuhan, China in December 2019, COVID-19 developed into a pandemic within the span of a year, and has currently become endemic in many countries [3]. As of 1 December 2024, the cumulative numbers of COVID-19 cases and disease-related deaths stand at 776,973,432 and 7,077,725, respectively [3].

Initial disease control and prevention efforts focused on the identification of infected individuals through aggressive testing and extensive contact tracing as well as on the reduction in exposure through the isolation of infected individuals, quarantine of exposed individuals, lockdowns, travel bans, social distancing, and the wearing of masks [1]. Later efforts focused on vaccination as a better strategy for disease control. There are currently several effective vaccines, with the Pfizer/BioNTech Comirnaty vaccine, Moderna COVID-19 vaccine, Oxford–AstraZeneca COVID-19 vaccine, Johnson & Johnson COVID-19 vaccine, Sinopharm COVID-19 vaccine, and the Sinovac-CoronaVac COVID-19 vaccine being some of the most administered vaccines in the world [3,4]. By the beginning of 2024, more than 13 billion COVID-19 vaccine doses have been administered worldwide [3].

While all existing COVID-19 vaccines are leaky (i.e., not completely effective at preventing infection), most of them are considerably efficacious against severe morbidity and mortality [5]. An independent analysis by the FDA based on the initial phase 3 clinical data for the Pfizer COVID-19 vaccine suggests that the vaccine can be 95% efficacious in preventing severe morbidity and mortality [6]. Some vaccines may only block infections, whereas others are protective against severe symptoms and death, but not infection or transmission [7,8,9]. For example, one study examines which of these two types of COVID-19 vaccines is more effective at a population level [9]. Despite the widespread accessibility of vaccines that are continually being enhanced, the control of SARS-CoV-2 infections has been particularly challenging due to reinfections and the emergence of new strains or variants [10,11]. Therefore, the waning of vaccine-induced immunity and infection-acquired immunity (i.e., natural immunity) are key factors that merit attention in a study focused on the epidemiological impact of vaccines.

COVID-19 is perhaps the most heavily studied disease of the recent past. The impressive corpus of mathematical modeling on COVID-19 is based on a variety of approaches, including continuous-time deterministic (ordinary differential equations (ODEs) and partial differential equations (PDEs)), discrete-time deterministic, statistical, and stochastic models [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. Many of these studies have used compartmental models, where a population is categorized based on each individual’s disease status at a given moment [27]. While most studies employ a binary classification of symptomatic/asymptomatic infections [12,13,17,19,20,23], only a few account for differential morbidity based on the severity of symptoms [9,24,28,29]. A broader clinical spectrum of SARS-CoV-2 symptoms includes asymptomatic, presymptomatic, mild, moderate, severe, and critical infections [30].

In the mathematical modeling of epidemiology, the phenomenon of backward bifurcation entails dire epidemiological consequences. This undesirable outcome results from the disease-free equilibrium (i.e., steady state of the model free from infections) co-existing with two endemic equilibria (i.e., steady states of the model with persistent infections), one stable and the other unstable, when the reproduction number is less than unity [31,32,33]. The disease can become endemic even when $[eqn]$ (i.e., classical threshold), but is greater than a certain value $[eqn]$ , which lies in the open interval (0,1). Thus, additional measures may be required to prevent the disease from becoming endemic or to bring $[eqn]$ below $[eqn]$ [32].

Several COVID-19 models exhibit the phenomenon of backward bifurcation [34,35,36,37,38]. For example, one ODE-based compartmental model accounts for the effects of vaccination and waning immunity on three sub-classes of infected individuals (i.e., symptomatic, asymptomatic, and hospitalized) [34]. A backward bifurcation in the model occurs in the context of low vaccination coverage and infection of vaccinated individuals. Moreover, an extended SVIR model of COVID-19 that categorizes the infected individuals as symptomatic or asymptomatic shows that a backward bifurcation can result from the coupling of the imperfect nature of the vaccine with the waning of natural immunity [37]. Another study uses distributed delay equations to examine a compartmental model for an endemic scenario of COVID-19, with decline in immunity from both vaccination and natural infection [35]. The direction of the bifurcation at the unity of the reproduction number of this model depends on the immunity waning rates.

Another COVID-19 model follows an extended version of the SVEIR structure to represent disease endemicity as well as waning immunity, and categorizes the infected individuals as presymptomatic infectious and symptomatic infectious [38]. This model also undergoes the phenomenon of backward bifurcation, which can be eliminated by increasing vaccine efficacy or decreasing disease transmissibility. The study shows that the waning of vaccine-induced immunity has a more pronounced effect on disease dynamics than does the waning of post-recovery immunity [38].

All the above studies exemplify the importance of waning immunity on disease dynamics and control. However, none of those studies seems to investigate the dynamics of subsequent disease waves. While at least one latter wave appears in the epidemic trajectories depicted in two works [35,37], these studies do not discuss the implications of the successive infection waves. Different studies have explored COVID-19 disease waves, but using tools such as time-varying transmission rates [39], stochastic models [40], and statistical data analysis [41].

An epidemiological model for an emerging infection by Hewage et al. [28] comprises seven compartments that account for vaccination and differential morbidity. The distinct compartments encompass susceptible, vaccinated, exposed, mildly symptomatic, moderately symptomatic, severely symptomatic, and recovered individuals. The model examines the impact of vaccine hesitancy on disease burden by means of a functional form for vaccine coverage, a feature not used in other studies. This model, however, does not account for the waning of natural or vaccine-induced immunity, an important factor in disease dynamics [22,35,37,38].

In this study, we build upon the between-host model in a previous study [28] to understand the dynamics of COVID-19 in the presence of imperfect vaccination and the waning of both vaccine-induced immunity and infection-acquired immunity. The incorporation of waning immunity allows us to examine the effects of recurring infections. We also include three infected classes, namely mild, moderate, and severe cases, to understand how differential morbidity affects the epidemic trajectories. The primary aims of our study are presented in terms of the following research questions:

1. What are the implications of the bifurcation at the unity of the reproduction number?
2. How sensitive are the model parameters to the vaccinated reproduction number; cumulative infections; first peak of infections; and the time until the first peak?
3. How do the waning of vaccine-induced immunity and the waning of infection-acquired immunity affect the disease dynamics within the population? Which factor has the higher impact on disease dynamics?
4. How efficacious of a vaccine is required to balance the effects of waning immunity?

We perform a comprehensive mathematical investigation that includes some preliminary results on the model’s qualitative nature, the derivation of equilibria and vaccinated reproduction number, and a bifurcation analysis. Then, we conduct both a local and global sensitivity analysis to understand how the input parameters influence the uncertainty of some model outputs. We also investigate the herd immunity threshold for the disease and how the two types of waning immunity affect the temporal variations of the time-dependent effective reproduction number. We next perform computer simulations on the model to explore how the two types of waning immunity affect the course of an epidemic and subsequent infection waves. Lastly, we quantify the level of vaccine effectiveness required to compensate for the two types of waning immunity.

The organization of this paper follows: Section 1 provides a review of the literature with a background to the current study. Section 2 elaborates the development of the model. A mathematical analysis of the model is provided in Appendix A. We present the results of the sensitivity analyses and numerical investigations in Section 3. The conclusions of the study are outlined in Section 4.

2. Model Development

The current study examines the dynamics of COVID-19 infections through an ODE-based compartmental model. Our model classifies a population into seven mutually exclusive compartments (with no overlaps), depending on each individual’s disease status at time t: unvaccinated susceptible (S); vaccinated susceptible (V); exposed (E); infected and experiencing mild symptoms ( $[eqn]$ ); infected and experiencing moderate symptoms ( $[eqn]$ ); infected and experiencing severe symptoms ( $[eqn]$ ); and recovered (R). This model is based on supplementing the epidemiological model in a previous study [28] to include the waning of vaccine-induced immunity and infection-acquired immunity. We also exclude disease mortality from the compartment of moderately symptomatic individuals in light of COVID-19 disease dynamics [1,42]. Additionally, we assume that the vaccine coverage rate remains constant over time.

Figure 1 schematically illustrates the transmission of COVID-19 in a population. The set of ordinary differential equations that mathematically represents the disease dynamics within the community is given in System (1). These ODEs are formulated based on the flow of individuals into and out of the seven compartments and the overall community over time:

$$[eqn]$$

The current model consists of 16 parameters as defined in Table 1. For the convenience of subsequent mathematical analysis, we use the bulk parameters $[eqn]$ , $[eqn]$ , and $[eqn]$ to represent the total loss of the infected compartments $[eqn]$ and $[eqn]$ , respectively. The main assumptions used in formulating the model are provided below.

• Demographics: The per capita natural community mortality rate is $[eqn]$ . New individuals are recruited into the community at the rate $[eqn]$ , where $[eqn]$ . However, the population is not closed because of disease-induced deaths.
• Vaccination: Susceptible individuals are vaccinated at the rate p. The vaccine is both leaky and waning, but not all-or-nothing (i.e., providing no protection to a certain fraction of the vaccinated individuals and perfect lifelong protection to others [43]). The vaccine offers an individual some protection against infection, but not against severe disease or death in the event of an infection. The level of protection that the vaccine provides against infection is $[eqn]$ , which is defined as the probability of a vaccinated individual being successfully protected from infection. The vaccine-induced immunity wears off at the rate $[eqn]$ .
• Force of infection and incidence rate: The time-varying and density-dependent (per capita) force of infection of the model is given by $[eqn]$ . Unvaccinated susceptible and vaccinated susceptible individuals become exposed at the rates $[eqn]$ and $[eqn]$ , respectively. Thus, the disease incidence rate is given by $[eqn]$ .
• Disease transmission and progression: There is no vertical transmission of disease in the community. Disease transmission occurs by mass action that results from interactions between susceptible and infected individuals. The unvaccinated and the vaccinated (with a leaky protection) may acquire an infection. The disease incubation period is $[eqn]$ . Individuals sequentially progress from mild symptoms through to severe symptoms (i.e., $[eqn]$ ). Individuals progress from $[eqn]$ to $[eqn]$ at the rate $[eqn]$ , and from $[eqn]$ to $[eqn]$ at the rate $[eqn]$ .
• Differential morbidity: The more severe the symptoms, the less likely it is that an individual will mingle with the community and will be able to transmit the disease (i.e., $[eqn]$ ). It takes longer to recover from more severe infections (i.e., $[eqn]$ ). Only the individuals with severe symptoms succumb to the disease (at the rate $[eqn]$ ).
• Recovery and natural immunity: Individuals in $[eqn]$ , $[eqn]$ , and $[eqn]$ recover at the rates $[eqn]$ , $[eqn]$ , and $[eqn]$ , respectively. Recovered individuals transition to the same recovered compartment (R), irrespective of the severity of previous morbidity. Recovery from an infection does not provide permanent immunity. Infection-acquired immunity wanes at the rate $[eqn]$ ; the waning of natural immunity causes the recovered individuals to become susceptible again. This immunity does not depend on the severity of previous morbidity. Moreover, since the vaccine is effective merely against blocking an infection, only the infection-acquired immunity takes effect upon recovery, regardless of the vaccination status of individuals.

The units and values of the parameters used in the numerical investigations are presented in Table 2.

We refer the reader to Appendix A for an analytical investigation of Model (1). Appendix A.1 explores the qualitative nature of Model (1), including the existence and uniqueness of solutions, non-negativity and boundedness of solutions, and the positive invariance of the solution space. Appendix A.2 derives the disease-free equilibrium of Model (1) as well as the two important epidemiological metrics, namely the basic reproduction number ( $[eqn]$ , the reproduction number of the model in the absence of vaccination) and the vaccinated reproduction number ( $[eqn]$ , the reproduction number of the model in the presence of vaccination), along with some related results. Appendix A.3 analyzes the endemic equilibria of Model (1) and the direction of the bifurcation at $[eqn]$ (i.e., whether it is forward or backward). The model undergoes a backward bifurcation at $[eqn]$ , given that the parameters satisfy the inequality in Equation (A9). This condition is further examined with respect to some key factors, including the vaccination parameters and immunity waning rates.

3. Numerical Results

In Section 3.1, we investigate both the local sensitivity and global sensitivity of some model outputs with respect to the model parameters. Numerical simulations on the vaccinated reproduction number, effective reproduction number, prevalence of infections, and disease mortality are provided in Section 3.2. Section 3.3 outlines the quantification of results in relation to vaccine efficacy and immunity waning rates.

3.1. Sensitivity Analysis of the Model

Here, we use both local and global sensitivity analysis methods to examine the sensitivity of Model (1) to changes in its parameters. Local sensitivity analysis is helpful in assessing the sensitivity of a model output in the neighborhood of a specific parameter space [51]. Conversely, global sensitivity analysis investigates the sensitivity of a model output to variations in parameters in their entire ranges, thus providing more comprehensive insights [52].

Section 3.1.1 explores the sensitivity of $[eqn]$ with respect to the model parameters using the direct sensitivity analysis method that is based on partial derivatives. In Section 3.1.2, we use the Partial Rank Correlation Coefficient (PRCC) method in tandem with Latin Hypercube Sampling (LHS) to examine the global sensitivity of $[eqn]$ , cumulative infections over 1500 days, first peak of $[eqn]$ , and the time until the first peak of $[eqn]$ to changes in the model parameters.

3.1.1. Differential Sensitivity Analysis

Local sensitivity analysis aims to analyze the impact of small variations in the input parameters on a specific model output [51]. Here, we use the direct method (i.e., differential sensitivity analysis) to examine the local sensitivity of $[eqn]$ to changes in its parameter space over a small neighborhood.

Given a model parameter $[eqn]$ , the formula for the normalized forward sensitivity index of $[eqn]$ with respect to $[eqn]$ is given by $[eqn]$ [53]. The sensitivity index $[eqn]$ provides valuable insights into the percentage change in $[eqn]$ , which results from a small percentage change in the parameter $[eqn]$ . Specifically, $[eqn]$ has a biological interpretation: a 1% change in the parameter $[eqn]$ causes $[eqn]$ to change by $[eqn]$ %. If $[eqn]$ , then decreasing the value of $[eqn]$ lowers $[eqn]$ . If $[eqn]$ , then increasing the value of $[eqn]$ lowers $[eqn]$ .

The bar plot in Figure 2 illustrates the local sensitivity indices of the model parameters with respect to $[eqn]$ . For example, $[eqn]$ means that a 1% increase (or decrease) in the value of $[eqn]$ causes a $[eqn]$ % increase (or decrease) in $[eqn]$ . Since $[eqn]$ is independent of $[eqn]$ , we have $[eqn]$ . The most negatively sensitive parameters are $[eqn]$ and $[eqn]$ , whereas the most positively sensitive parameters are $[eqn]$ , and $[eqn]$ . At a moderate level, the parameters $[eqn]$ and p are negatively sensitive, while $[eqn]$ and $[eqn]$ are positively sensitive.

$[eqn]$ is much less sensitive to the variations in the parameters $[eqn]$ and $[eqn]$ . Changes in the incubation rate ( $[eqn]$ ) usually have a very mild impact on $[eqn]$ . However, it is somewhat surprising that changes in the disease progression rates ( $[eqn]$ and $[eqn]$ ) have a minor impact on $[eqn]$ . Even though $[eqn]$ and $[eqn]$ are the most sensitive, as demographic parameters, these are less likely to be receptive to interventions. Thus, the local sensitivity analysis shows that the vaccinated reproduction number can greatly be reduced with control measures that prioritize $[eqn]$ , and $[eqn]$ . Therefore, some effective interventions include enhancing the vaccine efficacy, increasing the recovery rates for mildly symptomatic individuals, and decreasing disease transmissibility.

3.1.2. Partial Rank Correlation Coefficient (PRCC) Sensitivity Analysis

Local sensitivity analysis provides only limited insights into the sensitivity of model outputs because of its inability to explore the variations in the parameters over all admissible ranges. As such, PRCC sensitivity analysis is an efficient approach to global sensitivity analysis based on rank correlations. We use the Latin hypercube sampling (LHS) method to sample the parameters for the PRCC approach. The LHS method is a stratified sampling procedure without replacement that allows us to survey the whole parameter space. A total of 1000 samples are drawn from the ranges given in Table 2 as uniform probability distributions between the minimum value of the range (min) and the maximum value of the range (max).

We next calculate various model outputs (such as the reproduction number and the peak of infections) for each of the generated sample sets of parameters. We assess the strength and direction of the relationships between each model output and the parameters using PRCC values as sensitivity indices. The PRCC method works by comparing the ranks of model parameters with the ranks of the model outputs [54]. PRCC values vary between −1 and 1; a positive PRCC value indicates a positive correlation between the relevant parameter and the model output, whereas a negative PRCC value indicates a negative correlation [54]. The closer the PRCC value is to 1 or −1, the more that parameter contributes monotonically to changes in the output.

Here, we use the PRCC method to perform a global sensitivity analysis of the model parameters to $[eqn]$ ; the sum of the cumulative number of infections across mild ( $[eqn]$ ), moderate ( $[eqn]$ ), and severe ( $[eqn]$ ) cases over the first 1500 days; the number of infections at the first peak, summed across the three infected classes $[eqn]$ , and $[eqn]$ ; and the sum of the times to $[eqn]$ , and $[eqn]$ . The monotonicity of the relationships was verified, as required for the PRCC analysis.

As depicted in Figure 3a, the most influential parameters to the vaccinated reproduction number are the vaccine efficacy ( $[eqn]$ ), vaccine coverage rate (p), and the vaccine-induced immunity waning rate ( $[eqn]$ ). The first two transmission parameters ( $[eqn]$ ) and the first two recovery parameters ( $[eqn]$ ) are also influential. It is expected that any reproduction number will be sensitive to transmission and recovery parameters. Intuitively, the vaccine efficacy and coverage rate are strongly negatively correlated with $[eqn]$ ; a more efficacious vaccine and a higher coverage rate should slow the spread of disease. The rate of waning of vaccine-induced immunity is strongly positively correlated with $[eqn]$ . This indicates that the faster the vaccine-induced immunity wanes, the faster the disease spreads.

The equation for $[eqn]$ does not include the disease-acquired immunity waning rate ( $[eqn]$ ), and hence $[eqn]$ does not appear in the PRCC plot for $[eqn]$ . This, however, does not mean that $[eqn]$ has no effect on the dynamics of the disease. In order to discover its effects, we next examine the overall disease burden, measured as the sum of the cumulative number of infections across $[eqn]$ , $[eqn]$ , and $[eqn]$ over the first 1500 days. In Figure 3b, both waning rates ( $[eqn]$ and $[eqn]$ ) show a nearly equal and high effect on this disease burden metric. Numerical simulations suggest that the infection curves undergo at least two waves over the course of the first 1500 days due to the effects of waning immunity. Thus, the high negative sensitivity indices for $[eqn]$ and $[eqn]$ here are reasonable. Moreover, the PRCC value of $[eqn]$ is slightly higher than that of $[eqn]$ , indicating that the infection-acquired immunity waning has a higher influence on the cumulative infections than does the vaccine-induced immunity waning.

In the trajectories of infections across most parameter samples, the infected classes peak twice, and in some cases, there are three or four peaks. However, the first peak is always the highest. In Figure 3c,d, we show that the vaccine efficacy and coverage rate ( $[eqn]$ and p) are very influential to the first (and largest) peak of infections as well as to the time to this peak. The respective recovery and transmission parameters are sensitive as well. The disease incubation rate ( $[eqn]$ ) has the highest PRCC value. Neither waning rate is strongly influential to the peaks across the parameter ranges. Interestingly, the bars representing the PRCC values for the peak (Figure 3c) for the time to this peak (Figure 3d), indicating how the parameters that drive the peak higher also cause the trajectory to reach that peak sooner.

Overall, the global sensitivity analysis reveals that the vaccine-induced immunity waning rate, vaccine efficacy, and the vaccine coverage rate ( $[eqn]$ ) drive the changes in the vaccinated reproduction number, whereas both waning rates, vaccine efficacy, and the vaccine coverage rate ( $[eqn]$ ) drive the changes in the cumulative number of infections. However, the disease incubation rate, transmission coefficients, and the recovery rates ( $[eqn]$ ) have a large influence on the peak of infections and the time to this peak. On a comparison of the impact of vaccine-induced immunity waning ( $[eqn]$ ) and infection-acquired immunity waning ( $[eqn]$ ) on the model outputs, $[eqn]$ is more influential for the changes in the cumulative infections, whereas $[eqn]$ is more influential for the peak of infections. Thus, some efficient disease control measures include producing vaccines with higher efficacy and longer protection, increasing the coverage of vaccination programs, reducing disease transmission, and expediting the recovery of infected individuals.

3.2. Computer Simulations

We now conduct some numerical simulations on Model (1) to understand the disease dynamics of COVID-19 under the effects of vaccination and waning immunity in order to propose the key factors that should be taken into account when planning control measures. The software packages used for the numerical simulations of the model include MATLAB R2023b and Mathematica 14.1. The initial values of the state variables are given below. SVE $[eqn]$ $[eqn]$ $[eqn]$ R99999001000

The vaccinated reproduction number, evaluated at the baseline parameter values (Table 2), is $[eqn]$ . This number indicates that an outbreak will occur. In addition, we obtain $[eqn]$ ; this inequality implies that vaccination will reduce the reproduction number.

The concept of herd immunity refers to the herd effect or the indirect protection against a communicable disease that susceptible individuals receive when a sufficiently large number of individuals in the population is immune to the infection [55]. Thus, the herd immunity threshold is defined as the minimum proportion of the population that must be inoculated to reach the herd immunity level [55,56]. A narrative review of research on herd immunity in fighting COVID-19 is provided in the study by Suryawanshi et al. [57]. Moreover, in-depth, data-driven analyses of herd immunity in the context of COVID-19 can be found in several studies [58,59]. Specifically, Randolph et al. [58] discuss pertinent metrics for evaluating the societal cost of achieving global herd immunity against COVID-19. In addition, some studies emphasize the challenges in achieving herd immunity against COVID-19 even with mass vaccination [60,61].

The herd immunity threshold of a disease, denoted H, is computed by means of the formula $[eqn]$ [56]. A modified version, known as the critical vaccination level, is

$$[eqn]$$

where $[eqn]$ is the effectiveness of the vaccine against transmission [55]. We have $[eqn]$ when the vaccine induces complete protection against infection (i.e., $[eqn]$ ). For Model (1), we obtain $[eqn]$ and $[eqn]$ . Thus, to attain herd immunity, a vaccination coverage of 52% is required for a flawless vaccine; a vaccination coverage of 74% is required for a vaccine with 70% efficacy. A drawback of this approach is that Equation (2) is a general formula that does not account for the specific dynamics of a model at hand, except the numerical value of its basic reproduction number. In particular, this formula fails to account for other important factors such as waning immunity.

We now use numerical simulations to show how $[eqn]$ changes against the three vaccination parameters, vaccine coverage rate (p), vaccine efficacy ( $[eqn]$ ), and waning rate for vaccine-induced immunity ( $[eqn]$ ). Figure 4a depicts that $[eqn]$ varies non-linearly with a substantial rational decay (resembling the decline observed in a hyperbolic rational function) when p increases over the given interval. This relationship indicates that $[eqn]$ drops sharply as p is increased up to 0.01. Specifically, the vaccination coverage rate should be greater than 0.01 in order to bring $[eqn]$ down from the critical level. Moreover, $[eqn]$ seems to level off at about 0.7 when $[eqn]$ .

$[eqn]$ decreases linearly with $[eqn]$ at a considerably steep slope (Figure 4b). However, when the remaining parameters are fixed (Table 2), $[eqn]$ cannot be reduced below unity for any admissible value of $[eqn]$ . This means that even with a perfectly effective vaccine, an outbreak cannot be avoided. Figure 4c indicates that $[eqn]$ is an increasing function of $[eqn]$ , but the slope of the curve decreases as $[eqn]$ increases. Thus, the reproduction number can be substantially decreased with vaccines that remain effective for at least 11 months (i.e., $[eqn]$ , the point where the increase in $[eqn]$ seems to slow down). In particular, we need a vaccine with a protection that lasts longer than 2.7 years (i.e., $[eqn]$ , the point that corresponds to $[eqn]$ ) to prevent an outbreak in the community. Further, $[eqn]$ appears to level off at around 1.8 for $[eqn]$ .

The variations in $[eqn]$ with respect to pairs of vaccination parameters are depicted through the three-dimensional plots in Figure 5. Figure 5a shows that $[eqn]$ is significantly lower when the vaccine coverage rate takes larger values ( $[eqn]$ ) and the waning rate for vaccine-induced immunity takes smaller values on the range $[eqn]$ (that corresponds to a vaccine lasting up to 1.37 years). Given a certain level for p, $[eqn]$ increases steadily as $[eqn]$ increases. If $[eqn]$ , it is difficult to lower $[eqn]$ below 1, even if the vaccine-induced immunity does not wear off. Figure 5b indicates that given some level for $[eqn]$ , $[eqn]$ increases rapidly as $[eqn]$ increases. However, $[eqn]$ can be substantially reduced with higher values for the vaccine efficacy ( $[eqn]$ ) and lower values for the waning rate ( $[eqn]$ ). Specifically, the requirements for decreasing $[eqn]$ below unity are approximately $[eqn]$ and $[eqn]$ .

The effective reproduction number ( $[eqn]$ ) of a disease is a more pragmatic alternative to $[eqn]$ (or $[eqn]$ ) [62]. While $[eqn]$ gives the average number of secondary infections caused by an initial case in a completely susceptible population, $[eqn]$ estimates the expected number of secondary infections produced by an infectious case at any time point t during the course of an epidemic [63]. Thus, $[eqn]$ provides more insightful information than $[eqn]$ does in decision making while an epidemic is underway (when some individuals in the population are no longer susceptible). Several studies have estimated the effective reproduction number for COVID-19 using epidemic data [63,64,65]. $[eqn]$ is also defined as a data-driven measure of disease transmission that provides an estimate of the average number of new cases produced by each infectious individual on day t of an epidemic [66]. An interactive, continuously updated illustration of the monthly variations in $[eqn]$ and its 95% credible interval for the US and each state is provided on the webpage [66].

The effective reproduction number is formulated as $[eqn]$ , where $[eqn]$ is the basic reproduction number of the disease and $[eqn]$ is the proportion of susceptible individuals in the population at time t [62]. This formula can be revised based on the current model by changing $[eqn]$ to $[eqn]$ to account for the vaccinated individuals who become susceptible due to vaccine failure. Hence, we obtain $[eqn]$ , which yields $[eqn]$ ; this gives the relation $[eqn]$ .

Even though $[eqn]$ does not contain the waning rate for infection-acquired immunity ( $[eqn]$ ), $[eqn]$ does, making $[eqn]$ an apt tool for exploring the effects of $[eqn]$ on the reproduction number. Figure 6 illustrates how $[eqn]$ changes over time for different cases based on the presence of vaccination as well as the waning of vaccine-induced and disease-acquired immunity.

A comparison of Figure 6a,c as well as Figure 6b,f illustrates how vaccination helps decrease the reproduction number. Figure 6b shows that $[eqn]$ undergoes some oscillations leading to several outbreaks over a period of 1500 days, when the disease-acquired immunity wanes in the absence of vaccination. However, there are fewer oscillations and no subsequent epidemics when the disease-acquired immunity wears off under the presence of vaccination, despite the decay of vaccine-induced immunity (Figure 6f).

An epidemic can be prevented and the oscillations can be greatly suppressed with a vaccine that provides long-lasting immunity (i.e., $[eqn]$ ), as seen in Figure 6c,d. Figure 6a,c,e show a disease spread that is greatly controlled, when the recovery from an infection provides lasting immunity (i.e., $[eqn]$ ). The consistent downward trend in Figure 6c shows that vaccination is highly effective in disease control when there is no waning of either type of immunity. Overall, the waning of infection-acquired immunity ( $[eqn]$ ) has a higher influence on the oscillatory behavior of $[eqn]$ than does the waning of vaccine-induced immunity ( $[eqn]$ ). In the inevitable event of natural immunity waning, vaccination (even with decaying protection) plays a pivotal role in mitigating the effective reproduction number (and hence disease burden).

Figure 7 tracks the changes in the number of infected individuals ( $[eqn]$ , and $[eqn]$ ) over time. The disease propagation exhibits several waves over the course of 5000 days. The first peak of each curve is reached within 145–155 days. The peak for less intense infections is slightly earlier. At the first peak, the number of mild cases is about twice that of severe cases; there are approximately 67% more moderate cases than severe cases. These relative differences seem to remain consistent over the subsequent peaks. However, the latter waves are substantially flattened, and case numbers are also much lower during those waves. Eventually, the disease moves towards becoming endemic at lower case levels. The graph also confirms the absence of a periodic orbit (a solution curve that repeats itself over time) in the model.

Figure 8 shows how the number of individuals in each infected compartment ( $[eqn]$ and $[eqn]$ ) evolves over time for different values of the waning of vaccine-induced immunity ( $[eqn]$ ). The number of infections across all classes increases as $[eqn]$ increases. When $[eqn]$ , there are no subsequent waves after the first peak. Thus, a non-waning vaccine resolves the infection after the first peak, but when the vaccine protection wears off, reinfections cause subsequent peaks. The number of cases at the first peak of each $[eqn]$ (for $[eqn]$ ) increases by approximately one-third as $[eqn]$ varies from 0 to 0.02.

For each level of morbidity, the time until the first peak is approximately the same for all displayed values of $[eqn]$ , but the peak is delayed to a minuscule extent as $[eqn]$ decreases. However, the peaks for the subsequent waves are substantially delayed as $[eqn]$ decreases. This reveals that the faster waning of vaccine-induced protection causes the disease to peak earlier in the population. The relative drops in the peak number of infections (as $[eqn]$ decreases from 0.02 to 0) lessen as the symptoms become more lethal (Figure 8a–c), meaning that the peaks of mild (or moderate) infections drop more sharply than do the peaks of moderate (or severe) infections as the vaccine-induced immunity takes longer to wear off. The differences in the curves are negligible between the peak of the first wave and the beginning of the second wave, but the trajectories differ substantially thereafter.

Figure 9 shows the variations in $[eqn]$ and $[eqn]$ over time for different values of the waning rate for disease-acquired immunity ( $[eqn]$ ). If $[eqn]$ , the disease completely dies out after the first peak. Thus, no subsequent infection peaks occur in the absence of disease-acquired immunity waning. There are three discernible waves over the course of 1500 days, but for most values of $[eqn]$ , the third peak is not reached within the displayed time period. For each level of morbidity, the curves for various values of $[eqn]$ fall on top of each other until the first peak, but the curves become considerably different afterwards. This observation for the first wave in Figure 9 contrasts with the corresponding observation for the first wave in Figure 8.

Given each $[eqn]$ , the peak of the first wave delays very slightly as $[eqn]$ increases (contrary to the observations on the changes of $[eqn]$ ). However, the larger the value of $[eqn]$ , the earlier the peaks of the subsequent waves are (in agreement with the changes of $[eqn]$ ). Larger values of $[eqn]$ correspond to higher peaks in each wave. Thus, the disease burden increases greatly as the infection-acquired immunity wanes faster. The graphs in Figure 8 and Figure 9 reiterate that the infection trajectories undergo several waves, with each wave being less severe than the one before.

A comparison of Figure 8 and Figure 9 suggests that $[eqn]$ has a greater impact on the continuation of infections than does $[eqn]$ . The curves for $[eqn]$ in Figure 8 drop to almost 0 at the end of the first wave, whereas the curves for $[eqn]$ in Figure 9 drop only to about two-thirds of the peaks (for higher values of $[eqn]$ ) before rising again. The second peaks of the curves for $[eqn]$ (Figure 9) occur earlier than those of the curves for $[eqn]$ (Figure 8). The number of cases at the first peak of the curves for $[eqn]$ (Figure 9) vary in smaller ranges for each $[eqn]$ compared to the number of cases at the first peak of the curves for $[eqn]$ (Figure 8). Moreover, for larger waning rates, the peaks of the curves for $[eqn]$ (Figure 8) are higher than the peaks of the curves for $[eqn]$ (Figure 9) across all three infections, meaning that the vaccine-induced immunity waning ( $[eqn]$ ) has a higher influence on the first wave of infections than does the infection-acquired immunity waning ( $[eqn]$ ).

There are no subsequent waves after the first wave in either of the cases $[eqn]$ and $[eqn]$ or $[eqn]$ and $[eqn]$ . A separate numerical simulation suggests that when $[eqn]$ , there are no subsequent waves for any value of $[eqn]$ ; but, when $[eqn]$ , subsequent waves do occur for $[eqn]$ . Thus, the main driver of the infection waves is the waning of disease-acquired immunity ( $[eqn]$ ). As Figure 8 and Figure 9 indicate, infection-acquired immunity waning ( $[eqn]$ ) can lead to substantially large subsequent waves, when coupled with the vaccine-induced immunity waning ( $[eqn]$ ). Therefore, the waning of infection-acquired immunity has a higher negative impact on the infection dynamics than does the waning of vaccine-induced immunity.

We next examine the variation in the disease-induced deaths per day over time with respect to $[eqn]$ and $[eqn]$ , using the term $[eqn]$ . Figure 10a illustrates at least three waves for the surface plot for higher values of vaccine-induced immunity waning ( $[eqn]$ ). However, only one peak occurs when $[eqn]$ , suggesting that a vaccine with a non-waning immunity can largely control case fatalities. At the first peak of the surface plot that occurs after approximately 140 days, the number of deaths per day stands at about 110 for all displayed values of $[eqn]$ . The value of the peak increases almost imperceptibly as $[eqn]$ increases. Moreover, the number of deaths per day at the subsequent peaks remains well below 40.

In Figure 10b, lower values of infection-acquired immunity waning ( $[eqn]$ ) give rise to more peaks over the time period displayed. There is only one peak for the surface when $[eqn]$ , implying that the disease dies out after the first wave in the event that the infection-acquired immunity is permanent. The surface plot forms an erratic shape when $[eqn]$ and $[eqn]$ . The first peak occurs at approximately 150 days, and the peak seems to slightly decrease as $[eqn]$ increases (contrary to the observations on $[eqn]$ ). Moreover, the dynamics of Figure 8 and Figure 10a as well as the dynamics of Figure 9 and Figure 10b seem comparable.

There are several waves of infections and deaths in the population in the presence of the waning of both types of immunity. The primary driving force of these waves is the waning of infection-acquired immunity. Sufficiently large levels of infection-acquired immunity waning can lead to the continuation of infections, even when the vaccine-induced protection is non-waning. However, the waning of vaccine-induced immunity has a larger impact on the peak of the first wave of infections than does the waning of infection-acquired immunity.

3.3. Quantification of Results on Waning Immunity

This section provides a quantification of the effects on disease burden that result from the waning of vaccine-induced immunity ( $[eqn]$ ) and infection-acquired immunity ( $[eqn]$ ). Here, we intend to answer the research question of how efficacious of a vaccine is required to reduce the population level impact of waning immunity.

To this end, we approximate the efficacy of a vaccine ( $[eqn]$ ) required to reduce the peaks of the first wave of infections across all three classes ( $[eqn]$ , $[eqn]$ , $[eqn]$ ) by 40% (as compared to the reductions in those peaks with no immunity waning). First, we calibrate the value of vaccine efficacy ( $[eqn]$ ) that yields a 40% reduction in $[eqn]$ under different levels of $[eqn]$ , where $[eqn]$ denotes the peak of the infected compartment $[eqn]$ (see Table 3).

For a given value of $[eqn]$ , the vaccine efficacy required to decrease the peak of $[eqn]$ increases as $[eqn]$ increases (i.e., symptoms become more lethal). Given each $[eqn]$ , larger efficacy levels are required to lower the peak as $[eqn]$ increases. Specifically, if the vaccine-induced immunity wanes after 6.67 months, then $[eqn]$ is required to reduce the peak of mild infections by 40%. However, the efficacy level needed to achieve the same goal is much lower ( $[eqn]$ ), when the vaccine-induced protection lasts for at least 33.33 months.

We next calibrate the vaccine efficacy ( $[eqn]$ ) required to realize a 40% drop in $[eqn]$ for each $[eqn]$ under three different values of infection-acquired immunity waning ( $[eqn]$ ), namely 0.001, 0.003, and 0.005. The resulting quantifications are presented in Table 4. Given some value for $[eqn]$ , the vaccine efficacy required to attain the above goal slightly increases as the symptoms become more severe (i.e., i increases over $[eqn]$ for $[eqn]$ ). Moreover, higher levels of waning necessitate higher efficacies to reduce the peak of infections.

If the infection-acquired immunity takes about 33.33 months to wane, then a vaccine with $[eqn]$ is sufficient to reduce the peak of each infected class by 40%. However, if the immunity wanes faster (within 11.11 months), a vaccine with $[eqn]$ is required to cut off the peak of all infections by 40%. When the time for waning is even shorter (within 6.67 months), the efficacy of the vaccine must be increased up to $[eqn]$ to achieve the same goal. Interestingly, the relative increases in $[eqn]$ required in the case of increasing $[eqn]$ are lower than those required in the case of increasing $[eqn]$ .

Figure 11a depicts the percentage reduction in the peak of total infections with respect to the vaccine efficacy ( $[eqn]$ ) and vaccine-induced immunity waning ( $[eqn]$ ). The peak infections can only be reduced by an optimal level of 60% even in the extreme case: $[eqn]$ and $[eqn]$ . Overall, $[eqn]$ and $[eqn]$ are required to cut the peak of total infections by 50%. Figure 11b illustrates how the percentage reduction in the peak of total cases changes with respect to the vaccine efficacy ( $[eqn]$ ) and infection-acquired immunity waning ( $[eqn]$ ). When the vaccine efficacy is less than 0.6, the reduction is no more than 20%, regardless of the value of $[eqn]$ . However, $[eqn]$ and $[eqn]$ can yield very high levels of reduction (by 50%) in the peak of infections. The variations in Figure 11a are starker than the variations in Figure 11b. For a given level of vaccine efficacy, the level of $[eqn]$ required to achieve a certain reduction in $[eqn]$ is lower than the required level of $[eqn]$ . Thus, longer vaccine-induced immunity is more important in reducing the peak of infections than longer infection-acquired immunity.

A comparison of Figure 11a,b as well as Table 3 and Table 4 indicates that the vaccine efficacy necessitated to reduce the first peak of infections by a certain percentage is higher for a given level of vaccine-induced immunity waning ( $[eqn]$ ) than for the same level of disease-acquired immunity waning ( $[eqn]$ ). Thus, the peaks of infections are more sensitive to the variations in $[eqn]$ (given a sufficiently large value for $[eqn]$ ) than to the variations in $[eqn]$ (given a sufficiently large value for $[eqn]$ ). Moreover, these results are consistent with the global sensitivity analysis (Figure 3c) as well as with the simulations on the first peaks of infection curves (Figure 8 and Figure 9). Thus, the disease burden, especially the overwhelming of healthcare facilities, can be reduced substantially with vaccines that provide long-lasting protection.

4. Conclusions

In this study, we investigated an extended SVEIR model that incorporates differential morbidity as well as the waning of both vaccine-induced immunity and infection-acquired immunity. The model was examined, both mathematically and numerically, to determine the direction of the bifurcation at the unity of the reproduction number, the sensitivity of the model parameters to some model outputs, how the waning of both types of immunity affects the disease dynamics within the community, and the efficacy of a vaccine that is required to eliminate the negative impact of waning immunity.

We first showed some preliminary mathematical results on the qualitative nature of Model (1). Then, we computed the disease-free equilibrium of the model and derived a quadratic equation that yields its endemic equilibria. The vaccinated reproduction number ( $[eqn]$ ) of the model was computed by means of the NGM method. We performed a bifurcation analysis of the model using the sign of derivative method and the center manifold theory. Thus, we proved that the model exhibits the phenomenon of backward bifurcation at $[eqn]$ when it satisfies a certain inequality that consists of some model parameters. This perverse outcome disappears when the vaccine is perfect ( $[eqn]$ ) or when the disease-acquired immunity is permanent ( $[eqn]$ ). Moreover, the waning of infection-acquired immunity plays a more prominent role in the occurrence of a backward bifurcation than does the waning of vaccine-induced immunity.

A local sensitivity analysis of the model through the differential sensitivity method revealed that the parameters that are most sensitive to $[eqn]$ are the vaccine efficacy ( $[eqn]$ ), transmission coefficients for infected individuals with mild and moderate symptoms ( $[eqn]$ and $[eqn]$ ), and the recovery rate for infected individuals with mild symptoms ( $[eqn]$ ). The PRCC/LHS global sensitivity analysis showed that the vaccine efficacy ( $[eqn]$ ), vaccine coverage rate (p), and the waning rate for vaccine-induced immunity ( $[eqn]$ ) are the most sensitive parameters to $[eqn]$ . The results provided by the two methods are somewhat different because the differential sensitivity method examines the parameters one at a time in a local neighborhood of their point values, whereas the PRCC method examines all parameters simultaneously over their global ranges.

Moreover, the PRCC sensitivity analysis showed that both vaccine-induced immunity waning ( $[eqn]$ ) and infection-acquired immunity waning ( $[eqn]$ ) are highly sensitive to the cumulative number of infections; $[eqn]$ is more sensitive to this model output than $[eqn]$ . Therefore, some effective ways to curtail an epidemic of COVID-19 include producing vaccines with high efficacy and long-lasting protection, increasing the coverage rate of vaccination programs, quarantining infected individuals, and promptly treating infected individuals with therapy that hastens their recovery.

Numerical simulations were performed to assess the impact of vaccination and waning immunity on $[eqn]$ , the effective reproduction number ( $[eqn]$ ), the levels of infections with mild, moderate, and severe symptoms, and the disease-induced deaths. We calculated the herd immunity threshold of the disease and showed that an overall vaccination coverage of 74% is required to achieve the herd effect. Simulations on the effective reproduction number over time showed that $[eqn]$ largely controls the temporal variations in $[eqn]$ as opposed to $[eqn]$ . Specifically, the waning of immunity acquired from natural infection has a larger impact on maintaining $[eqn]$ below 1 than does the waning of immunity induced by a vaccine.

Figure 8, Figure 9, Figure 10 and Figure 11 demonstrated how infections and mortalities over time are much lower with smaller waning rates for vaccine-induced immunity and disease-acquired immunity. These simulations delineated that the waning of either immunity can lead to waves of disease spread over a long period of time. The faster the waning of either type of immunity, the larger the infection peaks become and the earlier the subsequent peaks occur. The quantification of results revealed the efficacy level that a vaccine must offer to reduce the peak of infections by 40% under different levels of vaccine-induced immunity waning ( $[eqn]$ ) and infection-acquired immunity waning ( $[eqn]$ ).

Overall, $[eqn]$ has a greater impact on the continuation of infections than does $[eqn]$ . However, the variations in $[eqn]$ affect the first peak of infections more potently than do the variations in $[eqn]$ . This result is consistent with the conclusions of a previous study [38], but that study does not examine subsequent infection waves. Hence, our study provides novel insights on the effects of the two types of waning immunity on overall disease dynamics.

There are several limitations to our study. The assumption that the progression of an infection takes the path $[eqn]$ may not reflect the clinical progression of symptoms in general for all individuals in a population. Alternate paths may include $[eqn]$ , $[eqn]$ , or $[eqn]$ . In addition, we did not examine the role of vaccines in reducing the severity of symptoms because of the complexity that such a phenomenon would bring forth, especially in terms of the mathematical analysis. Furthermore, we did not discuss virus mutation and the generation of mutant variants in the current study. Future research may also benefit from incorporating a time-varying functional form for the vaccine coverage rate.

In conclusion, we used a COVID-19 epidemiological model to show that the waning of disease-induced immunity impacts the overall disease dynamics more prominently than does the waning of vaccine-induced immunity. The study also revealed that sufficiently efficacious vaccines, high vaccine coverage rates, and the slow waning of vaccine-induced immunity can substantially reduce infections and disease mortality. These results suggest some important directions that could guide a future study, especially by examining whether a program encouraging the uptake of booster shots and the vaccination of recovered individuals can help keep infections at bay.