Mathematical assessment of the role of waning and boosting immunity against the BA.1 Omicron variant in the United States

1.   Introduction

Since December 2019, the world has been experiencing a devastating pandemic of a novel coronavirus (COVID-19), caused by SARS-CoV-2, on a scale never before seen since the 1918/1919 influenza pandemic ^[1]. As of mid July 2022, the SARS-CoV-2 pandemic has caused over $555$ million confirmed cases and over $6.35$ million deaths globally ^[2,3] (with the United States bearing the brunt of the burden, with over $88.6$ million confirmed cases and over $1$ million COVID-19 deaths) ^[3]. For most parts of the year 2020, the control and mitigation efforts against SARS-CoV-2 in the United States were restricted to the use of nonpharmaceutical interventions, such as social-distancing, quarantine of suspected cases, isolation of those with symptoms of SARS-CoV-2, use face coverings, community lockdowns, contact-tracing, etc. ^[4,5,6,7,8], until the Food and Drug Administration (FDA) provided Emergency Use Authorization (EUA) to two safe and highly-efficacious vaccines (developed by Pfizer Inc. and Moderna Inc.) in December of 2020 ^[9,10]. Both of the approved vaccines were primarily administered in two-dose regiments with three to four weeks apart, and each offer an estimated protective efficacy against symptomatic COVID-19 infection of about 95% ^[11,12]. Another vaccine, developed by Johnson & Johnson (administered as a single dose), received FDA-EUA in late February 2021 ^[13] (this vaccine has an estimated 75% efficacy in preventing severe/critical illness caused by COVID-19 ^[14]). The rapid development and administrative deployment of effective vaccines has played an extremely vital role in minimizing and mitigating the global burden of the pandemic ^[15,16]. Our study is focused on these three vaccines being used and administered in the United States.

Despite the rapid development and deployment of the effective vaccines, COVID-19 cases and mortality continued to rise in the United States for most part of 2021 (and even early 2022). This is largely due to emergence of deadly and highly-contagious SARS-CoV-2 variants of concern (notably the Alpha, Beta, Gamma, Delta and Omicron variants) ^{[17,18,19,20,21]}. Specifically, the emergence of the Omicron variant (B.1.1.529), in November of 2021, has dramatically changed the trajectory of the pandemic ^[22]. It is believed to be at least three times more contagious than Delta ^[22,23]. A subvariant of Omicron, BA.2 was first identified in the United States from a sample collected on December 14, 2021, in New Jersey ^[24]. It is believed to be more contagious than Omicron i.e. BA.1 ^[25,26].

Numerous clinical studies have shown that the efficacy of the SARS-CoV-2 vaccines wane over time (with estimated waning time of about 9 months) ^[8,27,28]. Consequently, the FDA approved the administration of booster shots, for all three vaccines, during August-November of 2021. Primarily, a booster dose (for persons aged 18 years and above), were approved because of waning vaccine effectiveness over time ^[29]. In late March 2022, the FDA authorized a second booster shot of COVID-19 vaccines for vulnerable populations in the U.S. (i.e., for people 50 years of age and older, and for individuals with certain immuno-compromising conditions who are at higher risk of severe disease, hospitalization and death). A second booster shot is equivalent to a fourth dose for people who received a Pfizer-BioNTech or Moderna mRNA series or a third dose for those who received the single-shot Johnson & Johnson vaccine. The goal of this study is to use mathematical modeling and analysis to assess the population-level impact of vaccination and booster shots (keeping in mind the waning efficacies of all the approved vaccinations) programs, based on using the three FDA-approved vaccines, on the dynamics of the Omicron SARS-CoV-2 variant in the United States. The focus of the study is on determining the minimum vaccination coverage (i.e., vaccine-derived herd immunity) needed to effectively curtail the spread of the highly-contagious Omicron variant in the United States.

Numerous mathematical models, of various types, have been formulated and used to gain insight and understanding on the dynamics of the COVID-19 pandemic (with majority being of the form of deterministic systems of nonlinear differential equations ^{[4,5,7,27,30,31,32,33]}). Here, too, a deterministic model will be developed and used to study the dynamics of the disease. A notable feature of the model to be developed is that it incorporates numerous pertinent aspects of the vaccination program and the current knowledge of the epidemiology of the COVID-19 pandemic, including the waning and boosting of both the vaccine-derived and natural immunity. The model will be parameterized using cumulative case data for the COVID-19 pandemic during the onset of the Omicron variant in the United States. The rest of the paper is organized as follows. The model is formulated in Section 2. The basic qualitative features of the model are also derived. Rigorous analysis of the model, with respect to the existence and asymptotic stability of its disease-free equilibrium, is carried out in Section 3. Expressions for vaccine-derived herd immunity thresholds are also derived in Section 3. The model is fitted with observed cumulative COVID-19 case data in Section 4. The methodology for implementing the data-fitting process and estimating the unknown parameters of the model is also described. Numerical simulations of the model are carried out in Section 5.

2.   Model formulation

To formulate the model for the transmission dynamics of SARS-CoV-2 in the presence of boosting and waning of both vaccine-derived and natural immunity, we split the total population at time $t$, denoted by $N(t)$, into mutually exclusive compartments of unvaccinated susceptible individuals ($S(t)$), fully-vaccinated susceptible individuals with high vaccine-derived immunity ($V_1(t)$), vaccinated susceptible individuals with moderate vaccine-derived immunity ($V_2(t)$) ^[34], vaccinated susceptible individuals with low vaccine-derived immunity ($V_3(t)$), exposed individuals (i.e., newly-infected individuals who are not yet infectious; $E(t)$), pre-symptomatic infectious individuals ($I_p(t)$), symptomatically-infectious individuals ($I_s(t)$), asymptomatically-infectious individuals ($I_a(t)$), hospitalized individuals ($I_h(t)$), recovered individuals with high infection-acquired natural immunity ($R_{n1}(t)$), recovered individuals with moderate infection-acquired natural immunity ($R_{n2}(t)$), recovered individuals with low infection acquired natural immunity ($R_{n3}(t)$), recovered individuals with high infection-acquired natural and vaccine-derived immunity ($R_{nv1}(t)$), recovered individuals with moderate infection-acquired natural and vaccine-derived immunity ($R_{nv2}(t)$) and recovered individuals with low infection-acquired natural and vaccine-derived immunity ($R_{nv3}(t)$). Thus,

Numerous clinical studies show that the vaccine-derived immunity against SARS-CoV-2 begin to wane after nine months of the receipt of the full vaccine doses ^[27,28,35]. Consequently, in our model formulation, individuals in the $V_1$ class (who enjoy high level of the protective efficacy of the vaccine) are those that are within nine months of receipt of full vaccine doses. Furthermore, individuals in the $V_2$ class are those who have received the full doses between 9 months to a year ago (hence, the vaccine efficacy is moderate). Finally, individuals in the $V_3$ class are assumed to have received the full vaccine doses from a year to two years ago (and the vaccine efficacy is very low). This study allows for the waning and boosting of vaccine-derived and natural immunity (boosting of natural immunity is assumed to occur due to treatment or the use of other immune-boosting supplements ^[36,37]).

The model is given by the following deterministic system of nonlinear differential equations, where a dot represents differentiation with respect to time $t$ (a streamlined/abbreviated flow diagram of the model is depicted in Figure 1, and the state variables and parameters of the model are described in Tables 1 and 2, respectively):

where,

is the infection rate. In (2.2), $\beta$ is the effective contact rate for individuals and $\eta_{j}$ (with $j = p, s, a, h$) is the modification parameter for the heterogeneity in the infectiousness of infected individuals in the presymptomatic ($I_p$), symptomatic ($I_s$), asymptomatic ($I_a$) and hospitalized ($I_h$) class, respectively.

In the model (2.1), $\Pi$ is the recruitment of individuals into the population, $\omega_{v_i}\, (i = 1, 2, 3)$ is the vaccine waning rate for vaccinated individuals in stage $V_i$, $\omega_{n_i}$ is the waning rate of natural immunity for recovered individuals in stage $R_{n_i}$, $\omega_{nv_i}$ is the waning rate of both vaccine-derived and natural immunity for individuals in stage $R_{nv_i}$, $\lambda$ is the infection rate (defined in (2.2)), $\xi_v$ is the per capita vaccination rate and $\mu$ is the natural death rate. Vaccinated individuals in $V_2$ and $V_3$ classes receive booster doses at the rate $\rho_{vi}$ ($i = 2, 3)$ and revert to the higher efficacy vaccination stage $V_1$ and $V_2$, respectively. Similarly, recovered individuals in the $R_{n2}$ and $R_{n3}$ classes receive immune booster at a rate $\rho_{n2}$ and $\rho_{n3}$, respectively (and revert, respectively, to stages $R_{n1}$ and $R_{n2}$). Individuals in $R_{nv2}$ and $R_{nv3}$ (that have both the vaccine-derived and natural immunity) receive a booster at a rate $\rho_{nv2}$ and $\rho_{nv3}$, respectively (and revert to $R_{nv1}$ and $R_{nv2}$, respectively).

The parameter $\varepsilon_{vi}$ is the average protective efficacy of the vaccine for vaccinated susceptible individuals in the $V_i\, (i = 1, 2, 3)$ compartment, while $\varepsilon_{n_i} \, (i = 1, 2, 3)$ is the efficacy of natural immunity to prevent recovered individuals (in the $R_{ni}$ class) from acquiring future SARS-CoV-2 infection and $\varepsilon_{nv_i} \, (i = 1, 2, 3)$ is the efficacy of natural and vaccine-derived immunity to prevent future SARS-CoV-2 infection of recovered individuals (in the $R_{nvi}$, $i = 1, 2, 3$, classes). Exposed individuals progress to the pre-symptomatic stage at the rate $\sigma_E$, and pre-symptomatic individuals progress to either become symptomatically-infectious, at a rate $r\sigma_p$ (where $0\leq r\leq 1$ is the proportion of these individuals that show clinical symptoms), or become asymptomatically-infectious, at the rate $(1-r)\sigma_p$. Symptomatic individuals are hospitalized at a rate $\phi_s$, and infectious individuals in stage $I_{j}$ (with $j = p, s, a, h)$ recover at a rate $\gamma_j\, (j = p, s, a, h)$. Finally, disease-induced mortality occur in the $I_p, I_s, I_a$ and $I_h$ classes at a rate $\delta_j\, (j = p, s, a, h)$.

Some of the main assumptions made in the formulation of the model (2.1) include:

(a) Homogeneous mixing: we assumed a well-mixed population, such that every member of the community is equally likely to mix with (and acquire infection from) every other member of the community.

(b) Vaccinated susceptible individuals (in the $V_1, \, V_2$ and $V_3$ classes) are assumed to have received the full required doses (i.e., two doses for Pfizer or Moderna vaccine, one dose for the Johnson & Johnson vaccine), and that enough time has elapsed for the body to develop immunity.

(c) The three SARS-CoV-2 vaccines that received FDA's Emergency Use Authorization (Pfizer, Moderna and Johnson & Johnson) are imperfect ^[10,13,38]. That is, the vaccines offer partial protective immunity (with average efficacy $0 < \varepsilon_{v_i} < 1$), which wanes over time (at a rate $\omega_{v_i}$), for $i = 1, \cdots, 3$ ^[28,35]. In other words, vaccinated individuals can experience breakthrough infection ^[39,40].

(d) We assumed gradual waning of both vaccine-derived and natural immunity over time, resulting, ultimately, in reverting to the wholly-susceptible class $S$ ^[28]. Moreover, the overall transitions from $V_1$ to $S$, $R_{n_1}$ to $S$ and $R_{nv_1}$ to $S$ hold gamma distribution ^[41].

(e) Vaccination is only offered to wholly-susceptible individuals or those who recovered naturally from COVID-19 infection but their natural immunity has waned completely or those recovered individuals who had acquired natural plus vaccine-derived immunity after recovering from COVID-19 infection but the immunity has completely waned over time. In other words, individuals who are currently infected are not vaccinated.

(f) Immunity level can be increased or strengthened, by using immunity boosters ^{[36,37,42,43]}, for the individuals in the $V_i, R_{n_i}$ and $R_{nv_i}$ ($i = 1, \cdots, 3)$ classes.

2.1.   Basic qualitative analysis

Before carrying out the asymptotic analysis and numerical simulations of the model (2.1), it is instructive to explore its basic qualitative features with respect to its well-posedness (i.e., with respect to the non-negativity, boundedness and invariance of its solutions). First of all, since the model (2.1) monitors the temporal dynamics of human populations, all its parameters are non-negative. It is convenient to define the following biologically-feasible region for the model (2.1):

where $N(t)$ is the total population. For the model (2.1) to be mathematically- and biologically-meaningful, it is necessary that the solutions of the model (2.1) remain non-negative for all non-negative initial conditions. That is, solutions that start in $\Omega$ remain in $\Omega$ for all time $t > 0$ (i.e., $\Omega$ is positively-invariant with respect to the model (2.1)). Furthermore, let

with $i = 1, 2, 3$, be the vector of initial solutions of the model (2.1). We claim the following result.

Theorem 2.1. Consider the model (2.1) with non-negative initial data $X(0)$. The region $\Omega$ is positively-invariant and bounded with respect to the model (2.1).

Proof. Adding all the equations of the model (2.1) gives

By the non-negativity of parameters for model (2.1), it follows from (2.3) that

Hence, if $N > \frac{\Pi}{\mu}$, then $\dot{N} < 0$. Thus, it follows, by applying a standard comparison theorem ^[44] on (2.4), that:

Hence, if $N(0) \leq \frac{\Pi}{\mu}$, then $N(t) \leq \frac{\Pi}{\mu}$. Thus, the solutions of the model (2.1) are bounded. Therefore, every solution of the model (2.1) with initial conditions in $\Omega$ remains in $\Omega$ for all time $t$. In other words, the region $\Omega$ is positively-invariant and attracts all initial solutions of the model (2.1).

The consequence of Theorem 2.1 is that it is sufficient to consider the dynamics of the flow generated by (2.1) in $\Omega$, since the model (2.1) is epidemiologically and mathematically well-posed ^[45] there.

3.   Asymptotic stability analysis of disease-free equilibrium

The disease-free equilibrium (DFE) of the model (2.1) is given by:

where,

with,

and all other components (for the infected and recovered compartments of the model) take the value zero.

The asymptotic stability property of the DFE (${\mathcal E}_0$) can be explored using the next generation operator method ^[46,47]. Using the notation in ^[46], it follows that the associated non-negative matrix of new infection terms ($F$) and the M-matrix of the linear transition terms ($V$) are given, respectively, by:

where $f_1 = \beta \eta_{p}\left(\dfrac{S^*+ A^*}{N^*}\right), f_2 = \beta \eta_{s}\left(\dfrac{S^*+A^*}{N^*}\right), f_3 = \beta \eta_{a}\left(\dfrac{S^*+A^*}{N^*}\right), f_4 = \beta \eta_{h}\left(\dfrac{S^*+A^*}{N^*}\right)$, with

It is convenient to define the quantity (where $\rho$ is the spectral radius):

The result below follows from Theorem 2 of ^[46].

Theorem 3.1. The disease-free equilibrium $({\mathcal E}_0)$ of the model (2.1) is locally-asymptotically stable (LAS) if ${\mathbb R}_{v} < 1$, and unstable if ${\mathbb R}_{v} > 1$.

The threshold quantity ${\mathbb R}_v$ is the vaccination reproduction number of the model (2.1), which measures the average number of new COVID-19 cases generated by a single infectious individual introduced into a population where a certain proportion is vaccinated. The epidemiological interpretation of Theorem 3.1 is that a small influx of COVID-19 cases will not generate a large outbreak in the community if the vaccination reproduction number (${\mathbb R}_v$) is brought to, and maintained at a, value less than unity.

In the absence of vaccination and other public health interventions, the vaccination reproduction number (${\mathbb R}_{v}$) reduces to the basic reproduction number (denoted by ${\mathbb R}_{0}$). That is,

3.1.   Global asymptotic stability of DFE: Special Case

Consider the special case of the model (2.1) with perfect vaccine protective efficacy against primary infection and re-infection and no waning of vaccine-derived and natural immunity (i.e., we consider the model (2.1) with $\varepsilon_{v_i} = \varepsilon_{n_i} = \varepsilon_{nv_i} = 1$, $\omega_{n_i} = \omega_{nv_i} = 0$, with $i = 1, 2, 3$). It is convenient to let:

Substituting $\varepsilon_{v_i} = \varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ and $\omega_{n_i} = \omega_{nv_i} = 0$ (with $i = 1, 2, 3$) into the model (2.1), it follows that the equation for the rate of change of the new compartment $\tilde{S}$ (defined in (3.5)) is given by (where the infection rate, $\lambda$, is as defined in (2.2)):

It can be shown that the special case of the model is positively-invariant and bounded in the region (as shown in Section 2.1)

Furthermore, the disease-free equilibrium of the special case of the model is given by:

For the aforementioned special case of the model, it can be seen that the associated next generation matrix of new infection terms, denoted by $\tilde{F}$, is given by (note that, for this special case, the next generation matrix of linear transition terms, $V$, remains the same, as given in (3.3). Further, $N^* = \Pi/\mu$):

Thus,

We claim the following result:

Theorem 3.2. Consider the special case of the model (2.1) with $\varepsilon_{v_i} = \varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ and $\omega_{n_i} = \omega_{nv_i} = 0$ (for $i = 1, 2, 3$). The disease-free equilibrium of the special case of the model $({\mathcal E}_{0R})$ is globally-asymptotically stable in $\Omega_*$ whenever ${\tilde{{\mathbb R}}}_v < 1$.

Proof. Consider the model (2.1) with $\varepsilon_{v_i} = \varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ and $\omega_{n_i} = \omega_{nv_i} = 0$ ($i = 1, 2, 3$). Further, let ${\tilde{{\mathbb R}}}_v < 1$. The proof is based on using a comparison theorem ^[48]. It can be shown, first of all, that the region $\Omega_*$ is positively-invariant and attracts all solutions of the aforementioned special case of the model (2.1) (as shown in Section 2.1). The equations for the infected compartments of the special case of the model (2.1) can be re-written in terms of the next generation matrices ($\tilde{F}$ and $V$) as below:

where,

and,

Since $S(t) \leq N(t)$ for all $t > 0$ in $\Omega_*$, it follows that the matrix $M$, defined in (3.10), is non-negative. Hence, Eq (3.9) can be re-written in terms of the following inequality:

It should be recalled from the local asymptotic stable result for the DFE (given in Theorem 3.1) that all eigenvalues of the next generation matrix ${\tilde F}V^{-1}$ are negative if ${\tilde{\mathbb R}}_v < 1$ (i.e., $\tilde{F}-V$ is a stable matrix). Thus, it can be concluded that the linearized differential inequality system (3.11) is stable whenever $\rho(\tilde{F}V^{-1}) < 1$. Hence, it follows that (for the linear system of ordinary differential equations (3.11)):

Substituting $E(t) = I_p(t) = I_s(t) = I_a(t) = I_h(t) = 0$ into the differential equations for the rate of change of the $R_{n_i}(t)$, $R_{nv_i}(t)$ and $\tilde{S}$ (with $\, i = 1, 2, 3$) compartments of the model (2.1) shows that (where $\tilde{S}^*$ is defined in (3.7)):

Thus, the DFE $({\mathcal E}_{0R}$) of the special case of the model (2.1), with $\varepsilon_{v_i} = \varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ and $\omega_{n_i} = \omega_{nv_i} = 0$ ($i = 1, 2, 3$), is globally-asymptotically stable in $\Omega_*$ whenever ${\tilde{\mathbb R}}_v < 1$.

Epidemiologically-speaking, Theorem 3.2 shows that, for the special case of the model (2.1) with $\varepsilon_{v_i} = \varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ and $\omega_{n_i} = \omega_{nv_i} = 0$ ($i = 1, 2, 3$), the disease can be eliminated from the community if the threshold quantity, ${\tilde{\mathbb R}}_v$, can be brought to (and maintained at) a value less than unity.

The global asymptotic stability of the disease-free equilibrium of the model (2.1) can also be established for another special case, as described below. Consider the special case of the model (2.1) in the absence of disease-induced mortality (i.e. $\delta_p = \delta_s = \delta_a = \delta_h = 0$) and no reinfection (i.e., $\varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ with $i = 1, 2, 3$). Setting $\delta_p = \delta_s = \delta_a = \delta_h = 0$ in the model (2.1), and adding all the equations of the model shows that $\dfrac{{dN}}{dt} = \Pi - \mu {N}$, from which it follows that $N(t) \rightarrow \dfrac{{\Pi}}{\mu}\, \, {\rm as}\, \, t \rightarrow \infty.$ From now on, we replace $N(t)$ with its limiting value, $N^* = \Pi/\mu$ (i.e., the standard incidence formulation for the infection rate is now replaced by a mass action incidence). Furthermore, it is convenient to define the following feasible region for the special case of the model (where $S^*$ and $V_i^*$, with $i = 1, 2, 3$, are as defined in Section 3):

It can be shown that the region $\Omega_{**}$ is positively-invariant with respect to the aforementioned special case of the model ^[49]. Further, for this special case of the model, it can be shown that the associated next generation matrices are given, respectively, by:

where (with $S^*$ and $A^*$ as defined in Section 3),

and,

It is convenient to define the following threshold quantity:

We claim the following result:

Theorem 3.3. Consider the special case of the model (2.1) in the absence of disease-induced mortality (i.e., $\delta_p = \delta_s = \delta_a = \delta_h = 0$) and no reinfection of recovered individuals (i.e., $\varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ with $i = 1, 2, 3$). The disease-free equilibrium of the special case of the model $({\mathcal E}_{0})$ is globally-asymptotically stable in $\Omega_{**}$ whenever ${\hat{{\mathbb R}}}_v < 1$.

The proof of Theorem 3.3, based on using a comparison theorem, is given in Appendix B.

3.2.   Derivation of vaccine-induced herd immunity threshold

Herd immunity, which is a measure of the minimum percentage of the number of susceptible individuals that need to be protected against the infection in order to eliminate community transmission of an infectious disease, can be attained through two main ways, namely natural immunity route (following natural recovery from infection with the disease) or by vaccination (which is widely considered to be the safest and the fastest way) ^[50,51]. For vaccine-preventable diseases, such as COVID-19, it is not practically possible to vaccinate every susceptible individual in the community due to various reasons, such as individuals with certain underlying medical conditions, infants, individuals who are pregnant, breastfeeding women or those who are unwilling to be vaccinated for COVID-19 due to some other reasons ^[30]. Therefore, it is crucial to determine the minimum proportion of the susceptible population that need to be vaccinated in order to protect those that cannot be vaccinated (so that vaccine-induced herd immunity is achieved in the population). Since we have three vaccination classes ($V_1$, $V_2$ and $V_3$), accounting for the three levels of vaccine-derived immunity (high, moderate and low), we will compute vaccine-derived herd immunity thresholds for the United States with respect to each of the vaccination classes. Specifically, we let

Here, $f_{v_1}$ represents the proportion of susceptible members of the population that have been fully-vaccinated (using any of the three approved vaccines), $f_{v_2}$ represents proportion of $V_2$ boosted and $f_{v_3}$ represents proportion of $V_3$ boosted, at the disease-free equilibrium $({\mathcal E}_0)$. Using the definition (3.13) in (3.4) gives:

which can be expressed in terms of the basic reproduction number (${\mathbb R}_0$) as:

Setting ${\mathbb R}_v = 1$ (i.e., the bifurcation point) in Eq (3.14), and simplifying, gives:

from which we can solve for the fraction fully-vaccinated (those in $V_1$ class), or received booster doses (i.e., those in $V_2$ and $V_3$ classes), at steady-state for each vaccinated class (denoted by $f_{v_i} {\rm with}\, \, i = 1, \cdots, 3$), in terms of the basic reproduction number, giving:

It follows from Eq (3.16) that ${\mathbb R}_v < (>)1$ if $f_{v_i} > (<)f_{v_i}^c$ $({\rm with}\, \, i = 1, 2, 3)$. Furthermore, ${\mathbb R}_v = 1$ whenever $f_{v_i} = f_{v_i}^c$ $({\rm with}\, \, i = 1, 2, 3)$. This result is summarized below:

Theorem 3.4. Vaccine-induced herd immunity (i.e., COVID-19 elimination) can be achieved in the United States, using any of the approved anti-COVID vaccine, if vaccination of susceptible individuals and boosting of vaccine-derived immunity resulted in $f_{vi} > f_{vi}^c$ (i.e., if ${\mathbb R}_v < 1$) for each corresponding $i$, for all $i = 1, 2, 3$. If $f_{vi} < f_{vi}^c$ (i.e., if ${\mathbb R}_v > 1$), then the vaccination program will fail to eliminate the COVID-19 pandemic.

Epidemiologically-speaking, Theorem 3.4 implies that the use of any of the approved COVID-19 vaccines can lead to the elimination of the pandemic in the United States if the proportion of susceptible individuals fully-vaccinated and with high level of vaccine-derived immunity (i.e., those in $V_1$ class) and boosted (i.e., those in $V_2$ and $V_3$ classes) at steady-state reached or exceeded the aforementioned critical threshold values. In other words, the SARS-CoV-2 pandemic is predicted to be eliminated in the United States if $f_{vi} > f_{vi}^c$ for each corresponding $i$, for all $i = 1, 2, 3$. On the contrary, the Vaccination program will fail to eliminate the pandemic if the proportion vaccinated (and boosted) at the disease-free equilibrium falls below the aforementioned critical herd immunity thresholds.

It should be mentioned that since the Pfizer and Moderna vaccines offer protective efficacy of about $95\%$ and $94\%$, respectively ^[52], and the Johnson & Johnson vaccine offers a protective efficacy of about $67\%$ ^[53], we set the average vaccine protective efficacy for individuals in the $V_1$ class to be

It is worth mentioning that the expression (3.17), for the average vaccine protective efficacy, is a dynamic quantity that depends on the actual combinations of SARS-CoV-2 vaccines adopted in the community at time $t$ (it should be mentioned that, in deriving the estimate for $\varepsilon_{v_1}$ in Eq (3.17), we used the values of the respective efficacies of the three vaccines estimated in ^[52,53] during the period around October, 2020 to January, 2021). Table 5 summarizes the assumed baseline efficacy levels for average vaccine-derived and natural immunity to be used in our numerical simulations. Using the baseline values of the fixed and fitted parameters in Tables 3 and 4, together with the baseline average vaccine-derived and natural immunity protective efficacy levels in Table 5, it follows from Eq (3.16) that the critical vaccine-derived herd immunity threshold for each of the vaccinated compartment is given, respectively, by $f_{v_1}^c = 0.59$, $f_{v_2}^c = 0.42$ and $f_{v_3}^c = 1.01$. In other words, based on the parameterization of the model (2.1) with the recent case data for Omicron BA.1 variant in the United States, population-level herd immunity can be achieved in the United States if the following conditions hold:

(a) at least 59% of the wholly-susceptible individuals are fully vaccinated (i.e., 59% of individuals in the $S$ class are fully-vaccinated and moved to the $V_1$ class);

(b) 42% of vaccinated individuals with moderate vaccine-derived immunity (i.e., those in $V_2$ class) are boosted;

(c) almost all of the vaccinated individuals whose level of vaccine-derived immunity is low (i.e., those in $V_3$) are boosted.

Hence, vaccine-derived herd immunity will be achieved via the aforementioned vaccination-boosting strategy that entails having at least 59% of the wholly-susceptible population to be fully-vaccinated followed by the boosting of an average of (42%+101%)/2 = 71.5% of the fully-vaccinated individuals with moderate and low vaccine-derived immunity.

Figure 2 depicts contour plots of the vaccination reproduction number (${\mathbb R}_v$), as a function of vaccination efficacy ($\varepsilon_{vi}$) and coverage of fully-vaccinated or boosted individuals at steady-state ($f_{vi}$), for $i = 1, 2, 3$. It follows from these plots that, for the overall average vaccine-protective efficacy set at 85% (as stated above), at least 59% of the wholly-susceptible population need to be vaccinated at steady-state to bring the vaccination reproduction number (${\mathbb R}_v$) below one (Figure 2(a)). For the case when the average vaccine-protective efficacy has waned to 50% (i.e., fully-vaccinated individuals now have moderate vaccine-derived immunity; here, $\varepsilon_{v2} = 0.5$, as given in Table 5), up to 42% of individuals in the $V_2$ class need to be boosted to bring the vaccination reproduction number to a value below one (Figure 2(b)). Finally, when vaccine-derived immunity has waned to the low level of 20%, the contour plot in Figure 2(c) shows that all of the fully-vaccinated individuals with low vaccine-derived immunity (i.e., individuals in the $V_3$ class) need to be boosted to bring the reproduction number to a value less than one. In summary, the results depicted in Figure 2 show that population-level herd immunity can be achieved in the United States via the implementation a vaccination program (based on using any of the three approved vaccines) that emphasizes the full vaccination of a sizable proportion of the susceptible pool (at least 59%) followed by the administration of booster doses to individuals in whom their vaccine-derived immunity has waned to moderate (at least 42%) or low (at least 100%) levels. Overall, our study shows that, for the case where the protective immunity offered by the vaccine for fully-vaccinated individuals in the $V_1$ class is 85%, vaccine-derived herd immunity can be achieved in the United States if at least 59% of the susceptible population is fully-vaccinated (with any of the three approved vaccines) followed by the boosting of at least 71.5% of the fully-vaccinated individuals whose vaccine-derived immunity has waned to moderate or low level.

It should be mentioned that for the case when the high level of the average vaccine-induced efficacy for individuals in the $V_1$ class is decreased to $55\%$, for instance (while the average vaccine protective efficacy for individuals in the $V_2$ and $V_3$ classes remain at the baseline level), our simulations showed that at least $91\%$ of the wholly-susceptible population need to be fully-vaccinated, followed by the (marginal) boosting (about $4\%$) of the vaccinated individuals whose vaccine-derived immunity has waned to moderate or low level. Thus, this study shows that lower average protective efficacy of the vaccines (for fully-vaccinated individuals) incurs higher requirement for the vaccination coverage of the susceptible population (followed by a correspondingly low level of boosting for the vaccinated individuals whose average vaccine protective efficacy has waned to moderate or low level) to achieve herd immunity. Vaccinating 91% of the wholly-susceptible population is, of course, not realistically feasible in large populations, such as the United States. Hence, it is imperative that highly efficacious vaccines are developed and used, and combined with boosting (at moderate to high levels) of vaccinated individuals whose immunity has waned to moderate or low level. In other words, using vaccines with higher protective efficacy (e.g., vaccines with 85% average protective efficacy, as computed in Eq (3.17)) incurs lower, and realistically attainable, requirement for the vaccination coverage (about 59%) and attainable (moderate to high) boosting level (about 71.5%) for the fully-vaccinated individuals whose vaccine-derived immunity has waned to moderate or low level. As of May 20, 2022, data from the CDC shows that about $66.5\%$ of the U.S. population is fully-vaccinated, and about $46.4\%$ of the population of fully-vaccinated individuals is boosted ^[54]. Thus, this study shows that, even for the scenario where the three vaccines offer such high average protective efficacy, a sizable proportion of the fully-vaccinated individuals need to be boosted in order to achieve vaccine-derived herd immunity in the United States (using the aforementioned combined vaccination-boosting strategy).

4.   Data fitting and parameter estimation

In this section, we fit the model (2.1) by using the available data for the observed cumulative COVID-19 cases for the U.S. (for the period November 28, 2021–February 23, 2022). The model (2.1) has several parameters, some of which are known from the literature (as tabulated in Table 3) and the remaining unknown parameters are obtained by fitting the model (2.1) with the cumulative case data obtained from the Johns Hopkins University COVID-19 repository ^[3]. The model was fitted using a standard nonlinear least squares approach, which involved using the inbuilt MATLAB minimization function "lsqcurvefit" to minimize the sum of the squared differences between each observed cumulative cases data points and the corresponding cumulative cases points obtained from the model (2.1) (i.e., $r \sigma_p I_p$). The unknown parameters which are estimated from the fitting are presented in Table 4.

The data fitting is done by splitting the available COVID-19 cumulative case data for the United States for the period from November 28, 2021 (when Omicron first emerged) to March 23, 2022 into two segments. The first segment of the data, from November 28, 2021 to February 23, 2022 (i.e., the region to the left of the dashed vertical cyan line), was used to fit the model (2.1) and to estimate the unknown parameters. The results obtained, depicted in Figure 3(a), show a very good fit for the model output (blue curve) and the observed cumulative case data (red dots). The data for the second segment of the data, for the period from February 24th, 2022 to March 23, 2022, was used for the cross validation of the fitted data. This also shows a very good fit for the model output (green curve) and for the remaining data points of the observed data (red dots) of Figure 3(a). This segment of the Figure 3(a) clearly shows that the model (2.1) cross validates the observed cumulative case data for the period from February 23, 2022 to March 23, 2022 perfectly (solid green curve). Furthermore, we show, in this figure, the prediction of the model for the cumulative COVID-19 cases for approximately a five-week period after March 24, 2022 (i.e., the region to the right of the dashed vertical black line), as illustrated by the solid magenta curve in Figure 3(a). The model was then simulated (using the fixed, fitted and assumed parameter baseline values in Tables 3, 4 and 5), and compared with the observed new daily cases data. The results obtained, depicted in Figure 3(b), show a very good fit for the daily new COVID-19 cases in the United States.

Furthermore, in fitting the model to observed cumulative case data for the United States, we set the average efficacies for individuals in the $V_i$ classes (with $i = 1, 2, 3)$ to their respective values in Table 5. Specifically, the value of the average vaccine efficacy for individuals in the $V_1$ class ($\varepsilon_{v_1}$) is set at 85% based on the expression given in Eq (3.17). The estimate for the average vaccine efficacy for individuals in the $V_2$ class ($\varepsilon_{v_2}$) was adapted from the empirical studies in ^[55,56]. First of all, these studies show variations in the timing of the attainment of the residual efficacy (after a few months of receipt of COVID-19 vaccination), which depend on the type of vaccine (Pfizer, Moderna or Johnson & Johnson). Further, the studies showed that the vaccine-derived efficacy of the Johnson & Johnson vaccine dropped from the initial $74.8\%$ to $59.4\%$ after five months of the receipt of the single-dose ^[55]. Similarly, the protective efficacy of the Pfizer vaccine decreased from $94.5\%$ to $75.7\%$ after 8 months of the receipt of the full doses of the Pfizer vaccine ^[56]. Finally, the efficacy of the Moderna vaccine decreased from $95.5\%$ to $84.3\%$ after 5 months of full vaccination ^[56]. Based on these estimates, we considered it plausible to set the average vaccine efficacy for individuals in the $V_2$ class to be $50\%$.

Additionally, numerous empirical studies have shown that the effectiveness of the Pfizer and Moderna vaccines against symptomatic COVID-19 is less than $20\%$ for the Omicron variant after the administration of the second dose of the vaccine ^[57,58]. Consequently, we deem it reasonable to set the average vaccine-derived efficacy for individuals in the $V_3$ class to be about $20\%$. Moreover, an empirical study showed that individuals who received the second booster dose of the Pfizer or Moderna vaccine have their vaccine effectiveness against the acquisition of breakthrough infection increased from $19\%$ to $49\%$ ^[59]. This is in line with the assumption in our model formulation regarding the transition of individuals from the $V_3$ class (where the average efficacy was set to be $20\%$) to the $V_2$ class (where the average vaccine efficacy is set at $50\%$).

5.   Numerical simulations

The model (2.1) will now be simulated to assess the population-level impact of waning and boosting of vaccine-derived and natural immunity on the dynamics of the Omicron variant in the United States. Unless otherwise stated, the simulations will be carried out using the baseline values of the parameters tabulated in Tables 3–5.

5.1.   Assessing the impact of waning of vaccine-derived immunity: with and without boosting

To assess the impact of waning of vaccine-derived immunity for this scenario, we simulate the model (2.1) using the following three (arbitrarily-chosen) levels of the parameters related to the waning of vaccine-derived immunity in the population:

(i) Low level of waning of vaccine-derived immunity: here, we consider vaccine-derived immunity to wane within 48 months (i.e., we set $\omega_{v_1} = \omega_{v_2} = \omega_{v_3} = 0.0007$ per day) but parameters related to natural immunity and combined natural and vaccine-derived immunity are kept at their respective baseline levels.

(ii) Baseline level of waning of vaccine-derived immunity: in this case, waning of vaccine-derived immunity is set to occur within 9 months (so that, $\omega_{v_1} = \omega_{v_2} = \omega_{v_3} = 0.0037$ per day) but parameters related to natural immunity and combined natural and vaccine-derived immunity are kept at their respective baseline levels.

(iii) High level of waning of vaccine-derived immunity: in this scenario, vaccine-derived immunity wanes within 3 months (i.e., $\omega_{v_1} = \omega_{v_2} = \omega_{v_3} = 0.0110$ per day), but parameters related to natural immunity and combined natural and vaccine-derived immunity are kept at their respective baseline levels.

For these simulations, all other parameters of the model (including those that involve the waning of natural immunity and combined waning of natural and vaccine-derived immunity) are maintained at their baseline values (given in Tables 3–5). Furthermore, these simulations are carried in the absence and presence of boosting of vaccine-derived immunity (recall that boosting of vaccine-derived immunity, maintained at baseline level, is achieved via the administration of the required doses of any of the approved SARS-CoV-2 booster vaccines used in the United States).

The simulation results obtained, depicted in Figure 4. First of all, these simulations also depict the fitting of the model's output for the daily new and cumulative cases with the observed data (used in Section 4) for the baseline scenario (as shown by the blue curves and the red dots in Figure 4). Furthermore, these simulations show that, in the absence of boosting of vaccine-derived immunity, waning of vaccine-derived immunity generally induces only a marginal impact on the average number of new daily COVID-19 cases in the United States, for each of the three waning levels considered in our simulations, in comparison to the baseline scenario. For example, under the fast waning scenario for vaccine-derived immunity (i.e., vaccine-derived immunity wanes within three months, but natural immunity is maintained at its baseline level) and no boosting of vaccine-derived immunity is implemented, the simulations show a marginal (about 2%) increase in the peak level of the daily new cases, in comparison to the peak baseline level (this is evident by comparing the blue and magenta curves in Figure 4$(a)$, and the zoomed-in version of the segments of the curves near the peaks shown in Figure 4$(b)$). For the slow waning scenario (i.e., if the vaccine-derived immunity wanes within 48 months, but natural immunity is maintained at baseline level), the increase in daily new cases at the peak (in comparison to the baseline) reduces to about 1.5% (compare the blue and green curves in Figures 4$(a)$ and $(b)$).

In the presence of boosting of vaccine-derived immunity (at baseline level), our simulations show a significant reduction in the average number of daily new cases at the peak recorded under the above waning scenarios without boosting of vaccine-derived immunity. For instance, for the case where vaccine-derived immunity wanes within three months (but natural immunity is maintained at baseline level), boosting of vaccine-derived immunity at the baseline level significantly reduces the increase in daily new cases at the peak (by about 90%), in comparison to the corresponding case without boosting of vaccine-derived immunity (compare the blue and magenta curves in Figure 4$(c)$, and the corresponding zoomed-in portions of the curves near the peaks shown in Figure 4$(d)$). Furthermore, under the slow waning scenario, boosting of vaccine-derived immunity at baseline level further increases the reduction in the peak daily new cases (compare the green and blue curves in Figures 4$(c)$ and $(d)$). Similar dynamics are observed (and illustrated) with respect to the cumulative number of new cases, for the three waning scenarios considered in these simulations without (Figures 4$(e)$ and $(f)$) and with (Figures 4$(g)$ and $(h)$) boosting of vaccine-derived immunity.

We further simulated the model to assess the impact of waning and boosting of vaccine-derived immunity (for the case where natural immunity is maintained at baseline) for the following two scenarios:

Scenario (a): Fast waning and slow boosting. Here, we assume that the waning of vaccine-derived immunity range between 3 to 6 months and the duration of boosting of vaccine-derived immunity range from 20 days to 180 days.

Scenario (b): Fast waning and boosting near the baseline level. Under this scenario, vaccine-derived immunity wanes within the same 3 to 6 months period (as in Scenario (a)), but boosting of vaccine-derived immunity is accelerated to be implemented within 10 to 20 days (i.e., near the baseline level of 14 days).

The results obtained are depicted in the form of heat maps for the vaccination reproduction number (${\mathbb R}_v$) of the model (2.1), as a function of the rates of waning ($\omega_v$) and boosting ($\rho_v$) of vaccine-derived immunity in Figure 5. This figure shows that, for the fast waning and slow boosting scenario (i.e., Scenario (a)), the values of the vaccination reproduction number lie in the range ${\mathbb R}_v\in [0.82, 1.21]$ (with a mean of ${\mathbb R}_v\approx 1.015)$, suggesting that the disease will persist in the population (this is in line with the theoretical result given in Theorem 3.1). In other words, this result shows that faster waning and slower boosting, in comparison to waning and boosting at baseline levels, increases the prospect for disease persistence in the population. For Scenario (b), our simulations (Figure 5$(b)$) show a marked decrease in the range of the reproduction number, with ${\mathbb R}\in [0.78, 0.89]$ (with a mean of ${\mathbb R}_v = 0.835$), suggesting possible elimination of the pandemic (in line with Theorems 3.1–3.3). Thus, boosting of vaccine-derived (near the baseline rate) enhances the prospect for pandemic elimination.

In summary, while the simulations in this section show that waning of vaccine-derived immunity generally induces only a marginal impact in the average number of new cases at the peak of the COVID-19 pandemic, boosting of vaccine-derived immunity (maintained at its baseline level) resulted in a dramatic reduction in the average number of new cases at the peak, in comparison to the case where boosting is not implemented. Furthermore, delay in boosting of vaccine-derived immunity, in comparison to the baseline level of boosting, could alter the trajectory of the disease from possible elimination (as measured by the vaccine reproduction number, ${\mathbb R}_v$, taking a value less than one) to persistence of the disease (as measured by the reproduction number being greater than one).

5.2.   Assessing the effect of waning of natural immunity: with and without boosting

Natural immunity can be boosted via treatment or the use of other immune-boosting supplements ^[36,37]. To assess the impact of waning of natural immunity, we simulated the model (2.1) using the following (arbitrarily-chosen) waning levels:

(i) Low level of waning of natural immunity: here, too, we consider natural immunity to wane within 48 months (i.e., we set $\omega_{n_1} = \omega_{n_2} = \omega_{n_3} = 0.0007$ per day), but vaccine-derived immunity and combined natural and vaccine-derived immunity are kept at their respective baseline levels.

(ii) Baseline waning of natural immunity: in this case, waning of natural immunity is set to occur within 9 months (so that, $\omega_{n_1} = \omega_{n_2} = \omega_{n_3} = 0.0037$ per day), but vaccine-derived immunity and combined natural and vaccine-derived immunity are maintained at their respective baseline levels.

(iii) High level of waning of natural immunity: here, too, natural immunity is assumed to wane within 3 months (i.e., $\omega_{n_1} = \omega_{n_2} = \omega_{n_3} = 0.0110$ per day), but vaccine-derived immunity and combined natural and vaccine-derived immunity are kept at their respective baseline levels.

For the simulations in this section, we set all other parameters (including those related to the waning of vaccine-derived immunity and combined waning of natural and vaccine-derived immunity) to their baseline values (given in Tables 3–5). The simulation results obtained, depicted in Figure 6, also showed that waning of the natural immunity in general only induces a marginal increase in the average number of new daily COVID-19 cases in the United States, in comparison to the baseline scenario (where the waning of natural immunity is assumed to occur within 9 months). In particular, if natural immunity wanes within three months and no boosting of natural immunity is implemented, the average number of new daily cases at the peak increases by about $4.2\%$, in comparison to the baseline scenario (compare the blue and magenta curves in Figure 6$(a)$, and the zoomed-in portions of the curves near the peaks, depicted in Figure 6$(b)$). An additional marginal increase in the average number of new daily cases at the peak is recorded under the slow waning scenario for the natural immunity, in comparison to the baseline scenario (compare the blue and green curves in Figures 6$(a)$ and $(b)$).

However, if natural immunity is boosted (at baseline level), our simulations show a marked reduction in the increase in the average daily new cases recorded at the peak, in comparison to the corresponding scenario without boosting of the natural immunity. Specifically, when natural immunity wanes within three months and boosting of natural immunity is implemented (maintained at its baseline level), the increase in the average number of new daily cases at the peak (in comparison to the baseline) reduces to about 1.6% (compare the blue and magenta curves in Figure 6$(c)$, and the zoomed-in portions near the peak depicted in Figure 6$(d)$). This represents an approximately 62.2% reduction in the average daily new cases at the peak, in comparison to the corresponding scenario where natural immunity is not boosted. It should be mentioned that boosting of vaccine-derived immunity (at baseline) plays a more significant role in reducing the average number of new daily cases, in comparison to the corresponding boosting of natural immunity (this can be seen by comparing the corresponding peaks in Figures 4 and 6). In particular, while boosting of vaccine-derived immunity (at baseline) will lead to about $90\%$ reduction in the number of new daily cases at the peak, boosting of natural immunity (at baseline) will lead to about $62\%$ reduction in the number of new daily cases at the peak).

Further significant reductions in the average number of new daily cases are recorded if the natural immunity wanes at a slower rate (compare the blue and green curves in Figures 6$(a)$ and $(c)$ or $(b)$ and $(d)$, without and with boosting of natural immunity). We illustrated similar dynamics with respect to the cumulative number of new cases without (Figures 6$(e)$ and $(f)$) and with (Figures 6$(g)$ and $(h)$) boosting of natural immunity.

In summary, like in the case of waning of vaccine-derived immunity discussed in Section 5.1, the simulations in this section show that while the waning natural immunity only causes a marginal increase in the average number of new cases at the peak, boosting natural immunity (at baseline) resulted in a significant reduction in the average number of new cases recorded at the peak, in comparison to the scenario where a strategy for boosting of natural immunity is not implemented in the community.

6.   Discussion and conclusions

The COVID-19 pandemic, caused by SARS-CoV-2, has made a significant impact on public health and the economy of almost every nation on earth since its emergence in December of 2019. The United States became the epicenter of the pandemic since late May, 2020 (recording the highest numbers of cumulative cases, hospitalizations and deaths). As of mid July, 2022, the virus had caused over $88.6$ million and 1 million deaths in the United States ^[2,3]. The rapid development, deployment, and administration of several safe and highly effective vaccines contributed significantly in curbing the spread of the virus worldwide. Three of these vaccines (the Pfizer-BioNTech, Moderna and Johnson & Johnson vaccines) have been approved by the FDA for use in the United States. The effectiveness of these vaccines in combating COVID-19 has been negatively affected by the emergence of various variants of SARS-CoV-2 (e.g., the Delta and Omicron variants). In particular, the Omicron (B.1.1.529) variant was declared a variant of concern by the World Health Organization in late November, 2021 ^[69], due to its exceptionally high transmissibility. Although all the available vaccines were developed for the original SARS-CoV-2 virus strain, they have been able to offer some level of cross-protection against other variants of concern. Furthermore, multiple studies have shown that the efficacy of vaccine-derived immunity wanes over time ^[28,35,70]. In order to overcome the waning effect of vaccine-derived immunity, booster vaccines were recommended by the CDC in November 2021 ^[42,71].

In this study, we developed a mathematical model to assess the population-level impact of the waning and boosting of vaccine-derived and natural immunity against the Omicron BA.1 variant of SARS-CoV-2 in the United States. The model was parameterized by fitting it to the observed cumulative COVID-19 case data for the United States for the period from November, 28, 2021 to February 23, 2022 ^[3]. We used the remaining segment of the available data (i.e., the segment from February 24, 2022 to March 23, 2022) to cross validate the model. This cross validation, together with simulations involving the new daily COVID-19 cases, showed a good match to the observed data.

The model was rigorously analyzed to gain qualitative insight into the dynamics and burden of the diseases. The analysis showed that the disease-free equilibrium (DFE) of the model is locally-asymptotically stable whenever the vaccination reproduction number (denoted by $\mathbb{R}_{v} < 1$) is below one. Using the baseline values of the fixed and estimated parameters of the model, we computed the numerical value of $\mathbb{R}_{v}$ during the period of the emergence and circulation of the Omicron variant (starting from late November of 2021). The computed value was $\mathbb{R}_{v} = 0.81$ (suggesting that Omicron was on a downward trajectory towards elimination in the United States). The numerical value of the basic reproduction number of the model (which is computed in the absence of any control measure implemented) was $\mathbb{R}_{0} = 2.051$. The computed value of $\mathbb{R}_{0}$ falls within the estimated range of $\mathbb{R}_{0}$ values for Omicron presented in several modeling studies (such as those in ^{[72,73,74,75]}). We showed that the disease-free equilibrium of the model is globally-asymptotically for two special cases ((a) when the vaccines offer 100% protection against acquisition of infection and no reinfection and waning of immunity occurs and (b) disease-induced mortality is negligible and reinfection does not occur) when the associated vaccination reproduction number is less than one. The epidemiological implication of this global asymptotic stability result is that the SARS-CoV-2 pandemic can be eliminated if the the associated vaccination reproduction number can be brought to (and maintained at) a value less than one (in other words, having the value of this reproduction threshold less than one is necessary and sufficient for the elimination of the pandemic in the United States).

Explicit expression for the vaccine-induced herd immunity threshold was derived, and we showed, using current data for COVID-19 cases in the United States, that, for the case where the three vaccines offer 85% average protective efficacy against the Omicron variant, vaccine-derived herd immunity will be achieved in the United States via a combined vaccination-boosting strategy that entails fully-vaccinating $59\%$ of the wholly-susceptible population combined with the boosting of at least $71.5\%$ of the population of the fully-vaccinated individuals whose vaccine-derived immunity has waned to moderate or low level. On the other hand, if the average protective efficacy offered by the three vaccines is reduced to a lower level, such as $55\%$ (as against $85\%$ above), at least $91\%$ of the wholly-susceptible population need to be vaccinated to achieve herd immunity. Furthermore, if the average vaccine efficacy offered by the aforementioned three vaccines drops to $50\%$ (or drops to a much lower value, as presented by Barnard et al. ^[16], or as stated in several studies on the effectiveness of these vaccines against Omicron ^{[57,72,76,77]}), then almost $100\%$ of the wholly-susceptible population need to be vaccinated in order to achieve vaccine-derived herd immunity in the United States. This very high level of vaccination coverage (i.e., vaccinating almost $100\%$ of the susceptible population) is not realistically attainable, especially in large populations such as the United States. Data related to COVID-19 from the United States Centers for Disease Control and Prevention show that, as of May 20, 2022, about $66.5\%$ of the U.S. population was fully-vaccinated and $46.4\%$ of this population received a booster ^[54]. Thus, our study suggests that, for the scenario that the three vaccines offer the reasonably high average protective efficacy of $85\%$ against acquisition of infection, herd immunity can realistically achieved in the United States by fully-vaccinating a moderate proportion (about $59\%$) of the wholly-susceptible and boosting about $71.5\%$ of this cohort in whom the vaccine-derived immunity has waned to moderate or low level.

We conducted extensive numerical simulations to assess the impact of waning and boosting of vaccine-derived and natural immunity for each three arbitrarily selected waning scenarios (slow, baseline, and fast). Our study showed, based on these simulations, that in the absence of boosting of vaccine-derived and natural immunity, waning of vaccine-derived and natural immunity only causes a marginal increase in the average number of daily cases (at the peak) and the number of cumulative COVID-19 cases, in comparison to the baseline scenario. In other words, we showed that waning of either vaccine-derived or natural immunity (or both) has only marginal impact, for each of the three waning scenarios we considered, on the dynamics of the SARS-CoV-2 pandemic (as measured in terms of increases in the average number of daily new cases recorded at the peak, in comparison to the case where baseline values of all the parameters of the model are used).

We also showed that if fully-vaccinated individuals with moderate or low level of vaccine-derived immunity are boosted (at baseline level), the effect of waning of immunity is a lot less pronounced, in comparison to the baseline scenario (in other words, dramatic reductions in the increase in the average number of daily new cases at the peak recorded (under the three waning scenarios) are achieved if both immunity types are boosted at baseline level, in comparison to the corresponding scenarios where the immunity wanes but no boosting is implemented. We further showed that boosting of vaccine-derived immunity is more beneficial (in reducing average number of new cases) than boosting of natural immunity. Specifically, for the fast waning scenario, boosting of vaccine-derived immunity (at baseline level) resulted in an approximate $90\%$ reduction in the average number of new daily cases at the peak, while boosting of natural immunity resulted in about $\approx 62\%$ reduction in the number of new daily cases at the peak (in comparison to the corresponding scenarios without boosting).

Furthermore, this study shows that boosting of vaccine-derived immunity (implemented near the baseline level) increased the prospects of altering the trajectory of the COVID-19 pandemic from persistence to possible elimination (even for the fast waning scenario of the vaccine derived-immunity) of the pandemic in the United States. Thus, the implementation of vaccination-boosting strategy greatly enhances the prospects of eliminating the COVID-19 pandemic in the United States.

In addition to the standard assumptions on which the model is built, some of the limitations of this study include the fact that we did not explicitly account for the impact of other control interventions (notably, the use of face masks, voluntary testing and detection of SARS-CoV-2 cases, isolation of confirmed cases, etc.), which also play important roles in the battle against the COVID-19 pandemic. Furthermore, this study assumes that the population is well-mixed and does not explicitly account for a number of heterogeneities, including age and risk structure, which may be relevant to gain insight into the dynamics of the disease. Furthermore, the current study did not account for the effects of other SARS-CoV-2 variants, including the BA.2 Omicron variant (which is more contiguous than the original BA.1 Omicron variant) ^[19,26]. We fitted our deterministic model to cumulative case data instead of raw (new daily case) data, which might lead to narrower confidence intervals and/or provide a misleading measure of uncertainty. Hence, the results of the study should be interpreted with some caution. Nonetheless, our study shows, overall, that the prospect for the effective control and mitigation (and, consequently, elimination) of the COVID-19 pandemic in the United States is very promising using a combined vaccination-boosting strategy, provide the vaccinate and boosting coverages are moderately high enough.

Acknowledgements

ABG acknowledges the support, in part, of the Simons Foundation (Award $\#$585022) and the National Science Foundation (Grant Number: DMS-2052363). CNN acknowledges the support of the Simons Foundation (Award #627346) and the National Science Foundation (Grant Number: DMS #2151870). SS acknowledges the support of the Fulbright Foreign Student Program. The authors are very grateful to the three anonymous reviewers for their very constructive comments.

Conflict of interest

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

Appendix A.   Full flow diagram of the model (2.1)

Appendix B.   Proof of Theorem 3.3

Proof. Consider the special case of the model (2.1) with $\delta_p = \delta_s = \delta_a = \delta_h = 0$ and $\varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ ($i = 1, 2, 3$). Further, let ${\hat{{\mathbb R}}}_v < 1$. The equations for the infected compartments of this special case of the model can be re-written in terms of the next generation matrices ($\hat{F}$ and $\hat{V}$) as follows:

where (with $S^*$ and $A^*$ as defined in Section 3),

and,

Since $S \leq S^*$, $V_1 \leq V_1^*$, $V_2 \leq V_2^*$ and $V_3 \leq V_3^*$ for all $t > 0$ in $\Omega_{**}$, it follows that the matrix $\hat{M}$, defined in (B.2), is non-negative. Hence, Eq (B.1) can be re-written in terms of the following inequality:

The proof is concluded the same way as in the proof of Theorem 3.2. Thus, the DFE $({\mathcal E}_{0}$) of the special case of the model (2.1) (with $\delta_p = \delta_s = \delta_a = \delta_h = 0$ and $\varepsilon_{n_i} = \varepsilon_{nv_i} = 1$ ($i = 1, 2, 3$)) is globally-asymptotically stable in $\Omega_{**}$ whenever ${\hat{\mathbb R}}_v < 1$.