key: cord-0737277-pt0pcwnj
authors: Mancuso, Marina; Eikenberry, Steffen E.; Gumel, Abba B.
title: Will vaccine-derived protective immunity curtail COVID-19 variants in the US?
date: 2021-09-09
journal: Infect Dis Model
DOI: 10.1016/j.idm.2021.08.008
sha: 086841a2013e68d0e9a6173363bcb60ac7a35617
doc_id: 737277
cord_uid: pt0pcwnj

Multiple effective vaccines are currently being deployed to combat the COVID-19 pandemic, and are viewed as the major factor in marked reductions of disease burden in regions with moderate to high vaccination coverage. The effectiveness of COVID-19 vaccination programs is, however, significantly threatened by the emergence of new SARS-COV-2 variants that, in addition to being more transmissible than the wild-type (original) strain, may at least partially evade existing vaccines. A two-strain (one wild-type, one variant) and two-group (vaccinated or otherwise) mechanistic mathematical model is designed and used to assess the impact of the vaccine-induced cross-protective efficacy on the spread the COVID-19 pandemic in the United States. Rigorous analysis of the model shows that, in the absence of any co-circulating SARS-CoV-2 variant, the vaccine-derived herd immunity threshold needed to eliminate the wild-type strain can be achieved if 59% of the US population is fully-vaccinated with either the Pfizer or Moderna vaccine. This threshold increases to 76% if the wild-type strain is co-circulating with the Alpha variant (a SARS-CoV-2 variant that is 56% more transmissible than the wild-type strain). If the wild-type strain is co-circulating with the Delta variant (which is estimated to be 100% more transmissible than the wild-type strain), up to 82% of the US population needs to be vaccinated with either of the aforementioned vaccines to achieve the vaccine-derived herd immunity. Global sensitivity analysis of the model reveal the following four parameters as the most influential in driving the value of the reproduction number of the variant strain (hence, COVID-19 dynamics) in the US: (a) the infectiousness of the co-circulating SARS-CoV-2 variant, (b) the proportion of individuals fully vaccinated (using Pfizer or Moderna vaccine) against the wild-type strain, (c) the cross-protective efficacy the vaccines offer against the variant strain and (d) the modification parameter accounting for the reduced infectiousness of fully-vaccinated individuals experiencing breakthrough infection. Specifically, numerical simulations of the model show that future waves or surges of the COVID-19 pandemic can be prevented in the US if the two vaccines offer moderate level of cross-protection against the variant (at least 67%). This study further suggests that a new SARS-CoV-2 variant can cause a significant disease surge in the US if (i) the vaccine coverage against the wild-type strain is low (roughly <66%) (ii) the variant is much more transmissible (e.g., 100% more transmissible), than the wild-type strain, or (iii) the level of cross-protection offered by the vaccine is relatively low (e.g., less than 50%). A new SARS-CoV-2 variant will not cause such surge in the US if it is only moderately more transmissible (e.g., the Alpha variant, which is 56% more transmissible) than the wild-type strain, at least 66% of the population of the US is fully vaccinated, and the three vaccines being deployed in the US (Pfizer, Moderna, and Johnson & Johnson) offer a moderate level of cross-protection against the variant.

Since its emergence late in December 2019, the novel coronavirus pandemic , caused by the RNA virus SARS-COV-2, has spread to every corner of the world. While other coronavirus pandemics have emerged in recent decades (including the 2002-2004 severe acute respiratory syndrome (SARS-CoV-1) and the 2012-present middle eastern respiratory syndrome (MERS-CoV) pandemics [60] ), COVID-19 represents the greatest global public health challenge since the 1918-1919 influenza pandemic. As of August 3, 2021, there have been over 198 million confirmed COVID-19 cases and over 4.2 million deaths worldwide, with the US being the most affected country (with over 35 million confirmed cases and now over 613,000 COVID-19 deaths) [21] .

The primary mode of transmission of COVID-19 is via respiratory droplets [16] , and the most common symptoms of the disease include dry cough, shortness of breath, fever, fatigue, muscle aches, and loss of taste or smell. The immunity, when compared to the wild-type SARS-CoV-2 strain, which does not contain any major mutations [15, 19, 83] . These concerns are compounded by the fact that, as more individuals became vaccinated, mask mandates and social-distancing measures were relaxed, or lifted altogether [3, 27, 78] . The SARS-CoV-2 strains sequenced in the US during the middle of May 2021 showed the B.1.1.7 variant accounting for nearly 70% of the variants sequenced, 8 .1% of the variants were identified as P.1, and 2.5% for the delta variant [20] . But, as of July 31, 2021, the Alpha variant was only 2.9% of the sequenced variants, while the Delta variant increased to over 83% [20] . The varying proportions of variant transmission, and how they respond to currently available vaccines, add more uncertainty to the effort to effectively suppress or eliminate the pandemic in the US.

Recent studies regarding the cross-protective immunity of existing vaccines on the SARS-CoV-2 VOCs have shown varying results [62, 83, 36, 84] . Some studies evaluating the cross-efficacy of the Pfizer-BioNtech vaccine on the B.1.1.7 and B.1.351 variants showed no significant differences of antibody neutralization when two doses of the vaccine were delivered 18-28 days apart [62, 84] . However, another study showed that less neutralizing antibodies were produced in individuals who acquired immunity through infection recovery (convalescent immunity) or only received a single dose of the vaccine [62, 4] . Another study compared the cross-efficacy of the B.1.1.7 and B.1.351 variants when two doses of the Moderna vaccine was administered 28 days apart [83] , which showed reduction in antibody neutralization for the B.1.351 strain, but not the B.1.1.7 strain. Further research is needed to validate and confirm these initial reports.

Despite the uncertainty regarding the level of cross-efficacy that current EUA vaccines provide against VOCs, some researchers showed that there is, at least, some degree of cross-protection [70, 15] . However, it is of major public health interest to understand if the current vaccination program will be enough to eliminate the pandemic in the United States, under a plausible range of variant transmissibility and vaccination cross-protection scenarios.

During the earliest days of the COVID-19 pandemic, before extensive experimental or case data was available, mathematical modeling was a useful tool to gauge the severity and potential impact that the disease would have on the community [42, 44, 51, 64, 66, 32, 85, 68, 31, 23, 48, 29, 46, 59, 38, 40] . With the increased availability and administration of the Pfizer, Moderna, and Johnson & Johnson vaccines in the US, as well as uncertainties surrounding the recently emerged SARS-CoV-2 variants, mathematical modeling can be a very useful tool to estimate the extent that cross-protective immunity from vaccines protect the community from the emerging COVID-19 variants. A number of recent modeling studies have used deterministic SEIR-type models to assess the potential impact of SARS-CoV-2 variants on the dynamics of COVID-19 in a population. For instance, Gonzalez-Parra et al. [34] used a two-strain SEIR modeling framework, which also includes disease transmission by asymptomatically-infectious individuals, to simulate COVID-19 dynamics in Colombia under varying levels of variant infectiousness (in comparison to the wildtype strain). A two-strain vaccination model was used by Betti et al. [8] to address the crucial question of whether a mutant SARS CoV-2 strain (mutant) could undermine vaccination efforts. Specifically, they estimated the time when a variant strain had overtaken the wild-type SARS CoV-2 strain in Ontario, Canada.

The current study uses a mathematical model to assess the impact of vaccination and vaccine-induced crossprotection against the B.1.1.7 and other SARS-CoV-2 variants circulating in the US as of June 2021. The two-strain and two-group model, which takes the form of a deterministic system of nonlinear differential equations, is formulated in Section 2.1. The model is fitted using COVID-19 mortality data, and a key unknown parameter of the model is estimated from the fitting, in Section 2.2. Basic qualitative features of the model, including the derivation of vaccinederived herd immunity threshold for the US and parameter sensitivity analysis, are discussed in Section 3. Numerical simulations of the model are also presented in this section.

In this section, a two-strain, two-group mechanistic model for the transmission dynamics of two SARS-CoV-2 strains (the wild-type and variant strains) in the US in the presence of vaccination will be formulated. Although multiple SARS-CoV-2 variants are co-circulating in the US, we consider only one variant strain in the model formulation. This is because almost all of the jurisdictions in the US that are currently experiencing major surges of SARS-CoV-2 have a dominant variant among the current co-circulating variant strains [20] . In formulating the model, the total population of individuals in the community at time , denoted by ( ), is stratified into two groups based on vaccination status, namely the total unvaccinated group, denoted by ( ) and the total vaccinated group denoted by ( ), so that ( ) = ( ) + ( ). The total unvaccinated population ( ( )) is further stratified into the mutually-exclusive compartments of individuals that are unvaccinated susceptible ( ), latent or "exposed" ( ), pre-symptomatically infectious ( ), symptomatically infectious ( ), asymptomatically infectious ( ), hospitalized ( ), and recovered ( ), with = 1, 2 (represents the SARS-CoV-2 strain), so that

Similarly, the total vaccinated population ( ( )) is sub-divided into the sub-populations of individuals that are vaccinated susceptible ( ( )), exposed ( ), pre-symptomatically infectious ( ), symptomatically infectious ( ), asymptomatically infectious ( ), hospitalized ( ), and recovered ( ), with = 1, 2. Thus,

Vaccinated individuals are assumed to have received the full required doses (i.e., they are fully vaccinated). The equations for the rate of change of the aforementioned state variables of the two-strain, two-group vaccination model for COVID-19 transmission dynamics and control in the US are given by the deterministic system of nonlinear differential equations given by Equation (A-8) of Appendix A. The flow diagram of the model is depicted in Figure 1 , and the state variables and parameters of the model are described in Tables A.7 and A.8, respectively.

In the formulation of the model (A-8), we assume a setting where two strains of SARS-CoV-2 are co-circulating in the population, namely:

(a) strain 1: this is the wild-type (original) strain. It is assumed to be the predominant strain circulating the US throughout 2020 [50] ;

(b) strain 2: this is a new dominant, variant strain assumed to be circulating in the US since the end of 2020 [22] ).

For mathematical tractability and convenience, we assume that susceptible individuals can be infected with either of the two strains, but not with both. Furthermore, we are not aware of real evidence for co-infection of multiple SARS-CoV-2 strains at the current moment. The three EUA vaccines in the US (i.e., the Pfizer, Moderna, and Johnson & Johnson vaccines) target the wild-type strain with efficacy , and are assumed to induce cross-protective efficacy ( ) against the dominant variant (strain 2) [62, 84, 36, 83] . Vaccination is offered to all eligible unvaccinated individuals, but for simplicity, the model does not consider vaccination of currently infected individuals (symptomatic or asymptomatic). In other words, the model only considers vaccinating those who are wholly-susceptible and those who recovered from the disease prior to being vaccinated. Although currently infectious individuals may be getting vaccinated, they comprise of extremely small fraction of the overall population at any given time, and are not included (for simplicity). Some of the other main assumptions made in the formulation of the model (A-8) include:

(i) Homogeneous mixing: the population is well-mixed, such that every member of the community is equally likely to mix with every other member of the community. Further, the waiting times in each epidemiological compartment are exponentially distributed [43] .

(ii) Because of the inherent age structure in the COVID-19 vaccine administration (where, in the context of the Pfizer vaccine, for instance, only people at age 12 or older are vaccinated [72] ), we include demographic (birth and death) parameters to account for the inflow of new susceptible individuals that are eligible for vaccination.

(iii) Vaccination is only offered to wholly-susceptible individuals or those who naturally recovered from COVID-19 infection (but not for currently-infected individuals).

(iv) Vaccinated individuals are assumed to have received the full recommended dosage (two doses for Pfizer or Moderna vaccine, one dose for the Johnson & Johnson vaccine), and that the vaccine administered was stored at the appropriate temperatures, and that enough time has elapsed for the body to develop immunity. (v) The vaccines are imperfect against infection (i.e, breakthrough infection can occur) [72, 71, 73] . The vaccines offer strong therapeutic benefits, vis a vis reducing severe disease, hospitalization, and death [72, 71, 73] .

(vi) Vaccine-derived and natural immunity may wane over time in individuals, in which they revert to the whollysusceptible class.

In this section, the model (A-8) is fitted using the cumulative mortality data from Johns Hopkins' Center for Systems Science & Engineering (from January 22, 2020 to March 6, 2021) [21] . Fitting is used to estimate the contact rate parameter for the wild-type strain (strain 1), 1 , during each wave of the COVID-19 pandemic in the US [21, 10] . The fitting is done in the absence of vaccination for the period from January 22, 2020 to January 21, 2021, and including vaccination for January 22, 2021 to March 6, 2021. We use a MATLAB routine to minimize the sum of squared differences between each observed cumulative mortality data point and the corresponding mortality point obtained from the model (A-8). COVID-19 mortality data, rather than case data, is used to fit the model since it is more reliable (the latter underestimates the true number of cases owing to the inability to track asymptomatic and pre-symptomatic cases, resulting from the absence of random wide-scale testing across the US) [18] . Daily vaccination data, for the period December 21, 2020 to March 6, 2021 [10] , was used in the fitting of the model to account for the transition of individuals from nonvaccinated compartments to attaining fully-vaccinated status.

The first case of the Alpha (B.1.1.7) variant was documented in the US on December 29, 2020 [22] , although the variant was likely introduced to the US prior to the detection. For the actual fitting of the model, we introduce 1, 000 cases of a variant strain (i.e., 2 (0) = 1, 000) on December 29, 2020 to estimate variant infections present prior to confirmation of its US detection. Furthermore, vaccination is introduced in the model (A-8) on January 22, 2021. Since our model does not explicitly consider the two-dose vaccine structure or account for any delay from inoculation to (partial) immunity, we introduce a delay from actual to modeled vaccination by constructing an inferred second dose time-series. That is, the available data only gives total doses delivered and does not disaggregate between first and second doses. Therefore, we assume that during the first 30 days of vaccine reporting, all administered doses are first Table 1 , illustrating the daily mortality generated from the model, in comparison to the observed daily mortality data for the US. For both panel, initial conditions used were (0) = 336 million, 1 (0) = 10 and all other state variables set at zero. The purple vertical line shows the time that 1,000 cases of the variant strain (strain 2) was introduced in the population (December 29, 2020), and assumed to be 1.56 times as infectious as strain 1 [25] . The brown vertical line shows the time at which vaccination was introduced into the model (January 22, 2021) using data from [10] .

doses. After 30 days, all first doses result in a second dose, and any remaining doses for that day are considered first doses. An inferred second dose time series can thus be constructed iteratively. Second doses were roughly 35% of the total doses administered each day following the first four weeks of vaccination [10] .

Furthermore, the following expression was used to determine the time-dependent vaccination rate, ( ):

This expression allows for the realistic possibility that individuals who have recovered from infection (and may have naturally-acquired immunity against future infection) will be vaccinated along with susceptible individuals. Thus, a proportion of the vaccination doses administered each day are delivered to wholly-susceptible individuals who have not developed prior COVID-19 immunity. The result of the data fitting of the model (A-8), using the fixed parameters in Table A .9, is depicted in Figure 2 . Figure 2 ). The contact rate for the wild-type strain, 1 , varies throughout the course of the pandemic due to implementation and compliance of non-pharmaceutical interventions (NPIs) and vaccination against it. Table 1 shows the values of the contact rate for the wild-type strain (strain 1), 1 , for each period of the pandemic. The first wave of the COVID-19 pandemic (January 22, 2020 to April 1, 2020) was the period with the highest contact rate, and occurred prior to most wide-scale NPI implementation (particularly the CDC's recommendation of face mask usage). The lockdown period (April 2, 2020 to June 15, 2020) was the period when numerous US states implemented stay-at-home policies, and was shown to have a lower contact rate. The second wave, which occurred during the early to mid summer of 2020 (June 16, 2020 to July 20, 2020) and was period when many states stay-at-home orders expired and businesses started returning to operation. This resulted in a higher contact rate, which in turn lead to more observed COVID-19 mortality. Many states began tightening restrictions on businesses again after seeing the rise in cases and deaths, which resulted in the post-second wave period (July 21, 2020 to September 18, 2020) having a decreased contact rate. As states lifted business restrictions and other COVID-19 related NPIs in mid-September, 2020 [78] , and many Americans traveled for the holidays, a large third wave of the pandemic was observed for September 19, 2020 to January 6, 2021, that had a similar contact rate as the second wave. As vaccination distribution commenced towards the final weeks of 2020, the contact rate once again decreased to a level similar to the post-second wave contact rate. We note that during the second and third waves, most US states still maintained a baseline level of NPI implementation and compliance (e.g., social distancing polices, face mask mandates). Thus, the contact rates for the second and third waves were not as high as the initial pandemic wave.

The model (A-8) has a unique disease-free equilibrium (DFE) given by

and all other components of 0 are zero. The asymptotic stability property of the disease-free equilibrium will be explored, for the special case of the model where ( ) is constant (i.e., we are considering the autonomous version of the model), using the next generation operator method [76, 26] . Using the notation in [76] , it follows that the non-negative matrix ( ) of new infection terms, and the matrix ( ) of new transition terms, are given, respectively, by:

where the matrices F 1 , F 2 , M 1 and M 2 are given in Appendix B, and 0 10×10 denotes the zero matrix of order 10. It follows that the vaccination reproduction number of the autonomous version of the model (A-8), denoted by  , is given by (where denotes the spectral radius of the next generation matrix FM −1 ):

where  [1] and  [2] (given in Appendix B) represent, respectively, the constituent reproduction numbers associated with the transmission of strains 1 and 2. The result below follows from Theorem 2 of [76] .

The threshold quantity,  , 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 implication of Theorem 3.1 is that a small influx if COVID-19 cases will not generate an outbreak in the community if the vaccination reproduction number ( ) is brought to (and maintained at) a value less than unity.

Herd immunity, which is a measure the fraction of susceptible individuals that need to be protected against infection in order to eliminate community transmission of an infectious disease, can be attained through two main ways: natural recovery from infection or from vaccination. However, the safest and fastest way to achieve herd immunity is through vaccination [7, 6] . Furthermore, at least 15% of the U.S. population currently cannot be vaccinated due to age restrictions [69] . Additionally, individuals who are pregnant, breastfeeding, have underlying medical conditions, as well as for other reasons, may be unable or unwilling to be vaccinated against COVID-19. Thus, it is critical to know what minimum proportion of the susceptible US population need to be vaccinated to achieve vaccine-derived herd immunity (so that the pandemic can be effectively controlled). In this section, we derive the expression for achieving vaccine-derived herd immunity against COVID-19 in the US based on widespread vaccination against the wild-type SARS-CoV-2 strain (i.e., strain 1).

In the absence of vaccination and other public health interventions (e.g., social-distancing, mask wearing, etc), the vaccination reproduction number ( ) reduces to the basic reproduction number (denoted by  0 ), given below:

where  1 and  1 are the constituent basic reproduction numbers for the wild-type (strain 1) and the variant (strain 2) SARS-CoV-2 strains, respectively, given by (where , with = 1, 2, ⋯ , 20 are given in Appendix B):

 1 = 1 [ 1 5 7 9 1 + 7 9 1 1 + 1 5 9 (1 − ) 1 .

For the case when all individuals in the population are vaccinated, the constituent reproduction numbers for the wildtype (strain 1) and variant (strain 2) SARS-CoV-2 strains become: The vaccination reproduction number for the (a) wild-type strain ( [1] , strain 1) with 1 = 0.24, (b) a variant strain that is 56% more transmissible than the wild-type strain ( [2] , strain 2), and (c) a variant strain that is 100% more infectious than the wild-type strain ( [2] , strain 2). Other parameter values used in this simulation are as given in Table A .9.

where expressions for ( = 1, 2, ⋯ , 8) are also given in Appendix B. Hence, for the case where all individuals in the population are vaccinated, the expression for the vaccination reproduction number of the model is given by:

It follows that the vaccination reproduction numbers of wild-type (strain 1,  [1] ) and variant (strain 2,  [2] ) can be re-written in terms of their basic reproduction numbers  1 and  2 , and reproduction numbers when the entire population is vaccinated, as below:

where = * * is the proportion of susceptible individuals in the community who have been fully vaccinated (at steady-state). The vaccine-derived herd immunity threshold for the wild-type and variant strains are computed by setting [1] = 1, [2] = 1 and solving for , respectively. We note that it is necessary to ensure, in this computation,

. Doing so gives:

It follows that the vaccine-derived herd immunity threshold for the COVID-19 pandemic in the US, in the presence of both the wild and variant SARS-CoV-2 strains, denoted by , is given by:

Using the fixed parameters in Table . 9 and a contact rate of 1 = 0.24 per day to account for transmission in the absence of nonpharmaceutical control measures (face mask usage, social distancing policies), it follows from Equation (7) that the vaccine-derived herd immunity threshold for the wild-type (strain 1) SARS-CoV-2 strain (in the absence of the variant strain) is = 1 = 0.59 if either the Pfizer or the Moderna vaccine (with efficacy = 0.94) is used. In other words, herd immunity against the wild-type strain can be achieved in the US, using either the Pfizer or Moderna vaccine, if at least 59% of the US population is fully-vaccinated with any of this vaccine. In this scenario, the herd immunity threshold rises to 64% if only the Johnson & Johnson vaccine (with = 0.67) is used. Thus, the wild-type SARS-CoV-2 strain can be effectively controlled in the US (in the absence of a co-circulating variant SARS-CoV-2 strain), using any of the three vaccines currently being administered in the US, if 59% (for Pfizer/Moderna vaccine) or 64% (for Johnson & Johnson vaccine) of the US populace is fully vaccinated.

For the case where the wild-type strain is co-circulating with the a dominant SARS-CoV-2 variant that is 56% more transmissible than the wild-type strain (such as the Alpha variant [25] ), our study shows that the vaccine-derived herd immunity threshold needed to effectively control both strains is 76% (i.e., = 0.76) using either the Pfizer or Moderna vaccine, or 93% when using only the Johnson & Johnson vaccine. On the other hand, if the wild-type strain is co-circulating with the Delta variant (which is 100% more transmissible than the wild-type strain), the herd immunity threshold needed to combat both strains increases to 82% if the Pfizer/Moderna vaccine is used. Figure 3 depicts contour plots of the vaccination reproduction number for the three scenarios with the wild-type strain only (Figure 3 (a) ) or in combination with a variant that is 56% (Figure 3 (b) ) or 100% ( (Figure 3 (c)) more transmissible than the wild-type strain, as a function of vaccine coverage and efficacy. This figure shows that the vaccination reproduction number decreases with increasing values of vaccine coverage and efficacy. This figure further shows that vaccine-derived herd immunity cannot be achieved using the Johnson & Johnson vaccine for the case where the wild-type strain is co-circulating with the Delta variant ( Figure 3 (c); since, for the 67% efficacy of this vaccine, the associated reproduction number cannot be brought to a value less than one, even with 100% vaccine coverage).

The model (A-8) contains numerous parameters, and uncertainties arise in the estimated values of some of the parameters (e.g., those related to vaccination and the dynamics of the variant strains) used in the numerical simulations. In this section, global sensitivity analysis is carried out, using Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) [55, 11, 56] , to determine the parameters that have the highest influence on the value of the chosen response function. Sensitivity analysis measures the extent that a response function changes with respect to variations in the input variables. For this study, the vaccination reproduction number for the variant strain ( [2] ) was used as the response function. One thousand bootstrap sample iterations were used for the LHS, and PRCC determined the correlation each individual parameter has on the response function ( [2] ). The baseline value of the contact rate for the variant strain, 2 , was chosen as 0.2495. The baseline values for the additional fixed parameters are as given in Table A .9. Each parameter was assumed to be uniformly distributed and range from 50% to 150% of the baseline value. The PRCC values generated for the parameters in the expression of the response function,  [2] , are tabulated in Table 2 . It follows from this table that our studies identifies the following four parameters that play a more dominant role in the size of the vaccination reproduction number for the variant strain ( [2] ): In other words, the parameters 2 , , , and have the greatest influence on the vaccination reproduction number of the variant strain ( [2] ). Hence, control measures against the variant strain (and, consequently, the disease in general) should be focused on these four parameters.

The model (A-8) is now simulated using the fixed parameters in Table A wild-type strain: Scenario 1 assumes the variant strain is the B.1.1.7 (Alpha) variant (which is 56% more transmissible than the wild-type strain) [25] ) and Scenario 2 assumes the variant strain is the Delta variant (which is 100% more transmissible than the wild-type strain) [19] . For mathematical tractability, it is assumed that recovered individuals do not acquire further SARS-CoV-2 infection. This assumption is consistent with other COVID-19 modeling studies [38, 40, 59, 58, 46, 29] . Furthermore, since less than 8% of fully-vaccinated individuals in the US have received the Johnson & Johnson vaccine [13] , simulations under Pfizer/Moderna vaccine administration can better estimate the future course of the pandemic.

Data from Korber et al. [50] suggests that the B.1.1.7. (Alpha) variant is 56% more transmissible than the wild-type strain in the US. Here, we simulate the dynamics of the two strains (wild-type strain and the Alpha variant) for various levels of the coverage level of fully-vaccinated individuals ( namely 50%, 66%, 75%) and cross-protective vaccination efficacy ( = 0.33, 0.5, 0.67, 0.75) offered by the Pfizer/Moderna vaccines. The cumulative COVID-19 mortality and cases, generated using the model for Scenario 1, are tabulated in Table 3 . This table shows that if 66% of the US population becomes fully-vaccinated with the Pfizer/Moderna vaccines, then under vaccination cross efficacy levels between 33% and 75% (i.e., ranges from 0.33 to 0.75), COVID-19 related deaths can be reduced between 27% and 44%, and cumulative cases can be reduced between 22% and 46%, when compared to the case where only 50% of the US population is fully-vaccinated. Furthermore, if as high as 75% of the US population becomes fully-vaccinated, then the reduction of COVID-19 related deaths and cases can increase to as high as 51% and 55%, respectively (when compared to the case when only 50% of the US population is fully-vaccinated). Table 4 tabulates the values and occurrence of the peak daily COVID-19 related deaths and new infections generated using the model. It can be seen from Table 4 that a future COVID-19 wave or surge can be averted if at least 66% of the US population is fully-vaccinated with either the Pfizer or Moderna vaccine, and the cross-protective efficacy the vaccine offers against the Alpha variant is at least 67%. Even if the cross-protective efficacy level is as low as 50%, a future pandemic wave can be avoided if at least 75% of the US populace is fully-vaccinated with either of the aforementioned vaccines.

Here, we simulated the case where the wild-type strain is co-circulating with the Delta variant (which is 100% more transmissible than the wild-type strain). For these simulations, the Pfizer and Moderna vaccines are used. Cumulative COVID-19 related deaths and cases for Scenario 2 are tabulated in Table 5 . This table shows that under 50% fullvaccination coverage and the same levels of cross-efficacy protection against the variant strain (i.e., ranging from 0.33 to 0.75), a 100% more infectious variant may result in 21% to 50% more COVID-19 related deaths and between 31% to 54% more cases, when compared to the 50% vaccination coverage level with the Alpha variant (i.e., a variant that is only 56% more infectious than the wild-type strain). Thus, the results in this table show a higher burden of the COVID-19 will be recorded in the US if the wild-type strain is co-circulating with the Delta strain, as against cocirculation with the Alpha variant (as expected). This table further shows that if the vaccination coverage is increased to 75%, then the COVID-19 related deaths and cases can be reduced by 24-56% and 20-58%, respectively, in comparison to only having 50% full-vaccination coverage, under this more infectious variant. Figure 5 shows the simulations of Scenario 2 when 75% of the US population is fully vaccinated with the Pfizer or Moderna vaccines. The values and occurrence of the peak daily COVID-19 related deaths and new infections are Table 6 . Table 6 shows that under this more infectious variant, a future pandemic wave can only be avoided if the US attains a full-vaccination coverage of 75% under the Pfizer and Moderna vaccines, and the cross-efficacy provided by the aforementioned vaccines against the Delta variant is at least 67%.

Since the emergence of COVID-19 in Wuhan City, China on December 31, 2019 [82, 86] , the pandemic has induced severe public health and socio-economic burden on an unprecedented scale. population has attained full-vaccination status as of August 3, 2021 [13] . However, the recent emergence of SARS-CoV-2 variants have posed major public health concerns in the US and globally. In particular, the B.1.1.7 (Alpha) and the B.1.617.2 (Delta) variants, that first emerged in the UK and India, respectively, have been circulating in many parts of the world, including the US [19] . Current data and studies show strong evidence that these variants are more easily transmissible, have higher rates of hospitalization and increased mortality rates, in comparison to the wild-type strain [45] . One of the greatest challenges associated with COVID-19 vaccination programs is the uncertainty of the level of cross-protection from currently-available COVID-19 vaccines offer against the variant strains. Numerous mathematical modeling studies, using varying modeling types and paradigms, have been conducted to provide insight into the transmission dynamics and control of COVID-19 in different populations. The modeling types used include statistical [48, 66, 64] , mechanistic/deterministic [58, 38, 40] , stochastic [42, 44, 51] , network [32, 68, 85] , and agent-based [23, 31] . The current study uses a deterministic system of 26 nonlinear differential equations to assess the population-level impact of the cross-protective efficacy of vaccination on SARS-CoV-2 variant strains. We used known parameter values for the wild-type strain and open-access COVID-19 death data to estimate the contact rate for infection transmission throughout the course of the pandemic in the United States. The disease-free equilibrium of the autonomous version (where the contact rate is assumed to be constant) of the resulting two-strain model (A-8) , which incorporates the vaccination of susceptible and naturally recovered individuals, is shown to be asymptoticallystable whenever the associated vaccination reproduction number of the model (denoted by  ) is less than unity. The threshold quantity  measures the average number of new SARS-CoV-2 cases generated by a typical infected individual introduced into a population that is partially-protected (via vaccination and natural recovery from infection). An explicit expression for the vaccine-derived herd immunity threshold (denoted by ) is derived. The herd immunity threshold is a measure of the minimum fraction of susceptible individuals in the population that need to be immunized Table 6 The and that the variant strain is 100% more infectious than the wild-type strain. A dashed line "-" indicates that the simulation did not produce a new peak of infections or deaths during the course of simulation. to ensure the protection of the individuals that cannot be vaccinated, and, consequently, leading to the suppression or halting of the COVID-19 pandemic.

Our study shows that if only the wild-type SARS-Cov-2 strain is circulating in the US, then at least 59% of the populace needs to be fully-vaccinated (with either the Pfizer or Moderna vaccine, each with estimated protective efficacy of at least 94%) to achieve the vaccine-derived herd immunity. In other words, if only the original wild-type SARS-CoV-2 strain is circulating in the US (which was the case prior to the emergence of the SARS-CoV-2 variants in the US around December 2020 [19] ), then at least 6 in every 10 people living in the US would need to be fully-vaccinated with the Pfizer/Moderna vaccines to achieve the vaccine-derived herd immunity. This herd value of vaccine-derived herd immunity, which is reasonably attainable in the US [13] , is consistent with the herd immunity threshold derived in previous studies, such as those in [40, 38] .

However, if the wild-type strain is co-circulating with the Alpha variant (which is 56% more transmissible than the wild-type strain), our study shows that the minimum threshold for achieving the vaccine-derived herd immunity (using either the Pfizer or Moderna vaccine) increases to 76%. The value increases to 82% if the wild-type is co-circulating with the Delta variant (which is 100% more transmissible than the wild-type strain). In other words, the emergence and circulation of the highly-transmissible Alpha and Delta variants contribute in raising the minimum value of the vaccinederived herd immunity threshold needed to effectively control, and mitigate the burden of, COVID-19 pandemic in the US. While achieving 60% vaccination coverage against the COVID-19 pandemic if only the wild-strain is circulating is certainly realistically-attainable in the US, achieving 76% or 82% coverage, if any of the two variants are circulating, is undoubtedly a tall order, for the numerous factors that contributed to the concerning level of vaccine hesitancy in the US [13] .

By implementing a global sensitivity analysis of the model, using the vaccination reproduction for the variant strain as the response function, this study identifies four parameters that have the highest influence on the SARS-CoV-2 variants (hence, of the COVID-19 pandemic) in the US. The identified parameters are the effective contact rate for the transmission of the SARS-CoV-2 variant strain contact rate ( 2 ), the cross-protective efficacy offered by the Pfizer/Moderna vaccine against the variant strain ( ), the coverage of fully-vaccinated individuals in the US ( ), and the relative infectiousness of fully-vaccinated infectious individuals, in comparison to the infectiousness of unvaccinated infectious individuals ( ).

We carried out extensive numerical simulations of the two-strain model to assess the population-level impact of various SARS-CoV-2 dynamic scenarios, based on varying the vaccination efficacy against the wild-type strain ( ), the vaccination cross-protective efficacy ( ) and the coverage of the fully-vaccinated population ( ) for the case where the wild-type strain is co-circulating with the Alpha or the Delta variant. These simulations show, as expected, that wide-scale vaccination is crucial to mitigating the burden and/or eliminating the COVID-19 pandemic in the US. For instance, for the Alpha variant, which was shown to be 56% more infectious than the wild-type strain [25] , we showed that increasing the coverage of fully-vaccinated individuals from the current 50% for the US to 75% can reduce anywhere from 33-51% of the cumulative mortality and 35-55% of cumulative cases (in comparison to when the vaccination coverage is 50%). For a SARS-CoV-2 variant that is more infectious than the Alpha variant, such as the Delta variant (which is 100% more transmissible than the wild-type strain), our simulations show that increasing the coverage of fully-vaccinated individuals from 50% to 75% could reduce anywhere from 24-56% of COVID-19 related deaths and 20-57% of cases, compared to the case with only 50% vaccination coverage.

Our study clearly shows that the vaccination coverage, cross-protective vaccination efficacy against variant strains, and the relative infectiousness of the variant strain, in comparison to the wild-type strain all greatly influence the dynamics of the COVID-19 pandemic in the US. Specifically, both the cross-protective vaccination efficacy against the variant and the vaccination coverage shift the peak of daily cases and deaths forward in time. This shifting increases with increasing values of the cross-protective efficacy and vaccination coverage. Our simulations show that, in general, the peak of daily deaths is expected to occur approximately 10 to 16 days after the peak of daily infections occurs.

One limitation of our study is the fact that our numerical simulations are based on the assumption that only the Pfizer or the Moderna vaccines are used in the population. In reality, three vaccines are being used in the US (the aforementioned two and the Johnson & Johnson vaccine), each with a presumably different level of cross-protective efficacy against the variant strains circulating in the population. Further, each of the three vaccines may differ in terms of their therapeutic properties (such as accelerating recovery, reducing risk of severe disease, hospitalization and deaths in breakthrough infections). Another limitation of our study is that the model is formulated to account for a single variant (i.e., we lumped all variants into a single dominant one). Although this assumption was necessarily made to enable mathematical tractability, extending the modeling framework to allow for the assessment of the dynamics of multiple SARS-CoV-2 variant strains (as is the case at the current time in the US, where both Alpha and Delta variants are co-circulating, but with the Delta clearly the prohibitive dominant one) would be quite useful. Since recent studies have shown that naturally-acquired immunity is less likely to produce an immune response to the existing variants than vaccination [62, 36] , another future extension of the modeling framework we developed in this study is to include the possibility that individuals who recovered from infection from the wild-type strain (i.e., those in the 1 class in our model) could acquire infection with the variant strain without losing their natural immunity against the wild-type strain.

The Alpha (B.1.1.7) variant was the most dominant SARS-CoV-2 variant circulating in the US throughout April and May 2021 [15, 20] , and was estimated to be 56% more transmissible than the wild-type strain [25] . The numerical simulations we conducted in this study showed that if a moderate level of vaccination coverage (≥ 66% of the US population) and moderately-high level of cross-protective vaccination efficacy is provided by the Pfizer and Moderna vaccines ( ≥ 0.67), then a surge of variant-induced COVID-19 infections and mortality will not occur (under the assumption that both the vaccine-derived and naturally-acquired immunity lasts at least two years). A recent study showed that cross-protective efficacy is as high as 93% against the Alpha variant when two doses of the Pfizer vaccine are administered, but drops to 33% when only one dose is received [37] . Thus, strictly following the two-dose regimen for the Pfizer and Moderna vaccines is essential to significantly reduce community transmission of COVID-19. The insights generated from this study align with the findings from other recent modeling studies for the dynamics of wild-type and variant strains of SARS-Cov-2 (also using SEIR-type modeling paradigm) [8, 34] .

Many US states relaxed control and mitigation measures against COVID-19, such as face mask mandates, once COVID-19 vaccination had become more widely available [27, 3] . Unfortunately, vaccination rates have since slowed down, since April 2021 [13] , and vaccinated individuals can still experience breakthrough cases that are infectious to others [14] . Since June 2021, the Delta variant has become the dominant SARS-Cov-2 variant in the US [20] , and is known to be more infectious than the Alpha variant [5] . For the Delta variant (which is 100% more transmissible than the wild-type strain), our simulations show that the COVID-19 pandemic can be effectively controlled in the US if 75% of the US population is fully-vaccinated (with either the Pfizer or Moderna vaccine) and the vaccine offers at least 67% cross-protective efficacy against the Delta variant. Thankfully, it was shown that the full recommended dosage of the Pfizer vaccine is 88% effective against the Delta variant [37] . However, if the vaccination coverage is lower, and a highly-transmissible variant (such as Delta) is co-circulating, the US still risks experiencing a severe COVID-19 wave in the future. As of August 2021, the US has experienced increased COVID-19 related cases and deaths [53] , with 97% and 99% of the COVID-19 related hospitalizations and deaths exclusively occurring in unvaccinated infected individuals, respectively [65] . Some US cities have since responded by resuming NPI (nonpharmaceutical intervention) mandates, such as mask wearing and social-distancing, in line with the new CDC recommendations [47] . Thus, it is critical to continue the practice of following NPIs, and encourage all eligible individuals to receive the full two-dose regimen of the Pfizer or Moderna vaccines. In summary, our study demonstrated that the prospect for the effective control and/or elimination of the SARS-CoV-2 variants (hence, the COVID-19 pandemic) in the US is promising provided the coverage of the fully-vaccinated individuals, and the level of cross-protection the vaccine offers against the SARS-CoV-2 variant strains, are high enough.

The authors declare no competing interests.

The two-strain, two-group vaccination model for the COVID-19 transmission dynamics and control in the US is given by the following deterministic system of nonlinear differential equations: 

The parameter Π is the recruitment rate (via birth or immigration) of individuals into the population. Individuals are vaccinated at a time-dependent rate ( ). The vaccine is assumed to produce protective efficacy 0 < < 1 against .

(A-10)

In (A-9) and (A-10), (with = 1, 2) is the effective contact rate for strain , ( ), ( ), and ( ) are, respectively, the modification parameters accounting for the assumed increase or reduction of contact rates of unvaccinated (vaccinated) pre-symptomatic, asymptomatic, and hospitalized individuals with strain , compared to symptomatic infectious individuals with strain . The modification parameter 0 ≤ ≤ 1 represents the assumed reduced infectiousness of individuals vaccinated against strain 1 (in comparison to unvaccinated infectious individuals with strain 1), while (0 ≤ ≤ 1) accounts for the assumed reduced infectiousness of strain 2 induced by the crossprotective of vaccination against strain 1.

Unvaccinated (vaccinated) individuals in the exposed classes (strain 1 or strain 2) progress to their respective strain's pre-symptomatic class at rate ( ) where = {1, 2} indicates the strain of infection. Non-vaccinated (vaccinated) individuals transition out of the pre-symptomatic classes at rate ( ). For strain 1, fraction 0 < ( ) < 1 of individuals develop symptoms while (1 − ) ((1 − )) individuals remain asymptomatic infectious. The fraction of unvaccinated (vaccinated) individuals infected with strain 2 who develop symptoms is 0 < ( ) < 1. Non-vaccinated (vaccinated) symptomatic infectious individuals recover at rate ( ), become hospitalized at rate ( ), and die from diseased-induced death at rate ( ). Non-vaccinated (vaccinated) asymptomatic infectious individuals recover at rate ( ). Non-vaccinated (vaccinated) Hospitalized individuals recover at rate ( ) and die from disease-induced death rate ( ). Non-vaccinated (vaccinated) individuals lose their naturally-acquired disease-induced immunity at rate ( ).

The state variables and parameters of the model are described in Tables A.7 and A.8, respectively. Vaccination efficacy against the wild-type strain (strain 1) 0 ≤ ≤ 1

Vaccination cross-protective efficacy against the variant strain (strain 2) ( ) Progression rate from unvaccinated (vaccinated) exposed (latently infected) individuals infected with strain to pre-symptomatic class 

where (noting that * = * + * , with * and * defined in Section 3): Using the approach in [39] , the vaccination reproduction numbers of the wild-type (strain 1) and variant (strain 2) SARS-CoV-2 strains are found from  respectively, where denotes the spectral radius. The expressions for  [1] and  [2] are given by: 

Executive Order GA 14

Executive Order GA 29

Executive Order GA 34

Effectiveness of the BNT162b2 Covid-19 Vaccine against the B.1.1.7 and B.1.351 Variants

Increased household transmission of COVID-19 cases associated with SARS-CoV-2 variant of concern B.1.617.2: a national case-control study

The concept of herd immunity and the design of community-based immunization programmes

Vaccination and herd immunity to infectious diseases

Could a new COVID-19 mutant strain undermine vaccination efforts? A mathematical modelling approach for estimating the spread of B.1.1.7 using Ontario, Canada as a case study

More Than 2.42 Billion Shots Given: Covid-19 Tracker

Bloomberg Covid-19 Vaccine Tracker Open Data (2021)

Sensitivity and Uncertainty Analysis of Complex Models of Disease Transmission: an HIV Model, as an Example

The Advisory Committee on Immunization Practices' Updated Interim Recommendation for Allocation of COVID-19 Vaccine -United States

COVID Data Tracker

COVID-19 vaccine breakthrough infections reported to CDC -United States

science/science-briefs/scientific-brief-emerging-variants.html?CDC_AA_refVal=https%3A%2F%2Fwww.cdc.gov% 2Fcoronavirus%2F2019-ncov%2Fmore%2Fscience-and-research%2Fscientific-brief-emerging-variants

How COVID-19 Spreads

Long-Term Effects

Pandemic Planning Scenarios

SARS-CoV-2 Variant Classifications and Definitions

Variant Proportions

COVID-19. Github repository

Gov. Polis and State Public Health Officials Announce First Case of COVID Variant of COVID-19 in Colorado

An agent-based model to evaluate the covid-19 transmission risks in facilities

BNT162b2 mRNA covid-19 vaccine in a nationwide mass vaccination setting

Estimated transmissible and severity of novel SARS-CoV-2 Variant of Concern 202012/01 in England

On the definition and the computation of the basic reproduction ratio 0 in models for infectious diseases in heterogenous populations

State of Arizona Executive Order 2021-06

State of Arizona Executive Order

To mask or not to mask: Modeling the potential for face mask use by the general public to curtail the COVID-19 pandemic

Theoretical Assessment of Public Health Impact of Imperfect Prophylactic HIV-1 Vaccines with Therapeutic Benefits

Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand. London: Imperial College COVID-19 Response Team

Using a real-world network to model localized covid-19 control strategies

LIST: Arizona cities with face mask requirements

Gatherings and Face Mask Order. Michigan Department of Health and Human Services

Comprehensive mapping of mutations in the SARS-CoV-2 receptor-binding domain that affect recognition by polyclonal human plasma antibodies

The Highly Contagious Delta Variant is On The Rise In The U.S. National Public Radio

A primer on using mathematics to understand COVID-19 dynamics: Modeling, analysis, and simulations

Existence of Multiple-Stable Equilibria for a Multi-Drug-Resistant Model of Mycobaceria

Towards achieving a vaccine-derived herd immunity threshold for COVID-19 in the U.S. Frontiers in Public Health

Effect of vaccination on household transmission of SARS-CoV-2 in England

Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts

The mathematics of infectious diseases

A stoachastic agent-based model of the sars-cov-2 epidemic in France

NERVTAG note on B.1.1.7 severity for SAGE. New and Emerging Respiratory Virus threads Advisory Group

Will an imperfect vaccine curtail the COVID-19 pandemic in the US?

These cities are reviving mask mandates with new CDC guidance, delta surge

Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator-days, and deaths by US state in the next 4 months

Florida Gov. Ron DeSantis issues statewide stay-at-home order

Tracking Changes in SARS-CoV-2 Spike: Evidence that D614G Increases Infectivity of the COVID-19 Virus

Early dynamics of transmission and control of covid-19: A mathematical modelling study. The Lancet Infectious Diseases

Initial report of decreased SARS-CoV-2 viral load after inoculation with the BNT162b2 vaccine

CDC says 7-day average of daily U.S. covid cases surpassed peak seen last summer

COVID-19 vaccines: Get the facts

A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code

Sensitivity and Uncertainty Analyses for a SARS Model with Time-Varying Inputs and Outputs

MDH lab testing confirms nation's first known COVID-19 case associated with Brazil P.1 variant

Mathematical assessment of the impact of non-pharmaceutical interventions on curtailing the 2019 novel Coronavirus

Could masks curtail the post-lockdown resurgence of covid-19 in the US?

COVID-19 MERS, & SARS. National Institute of Health

Arizona governor issues 'stay-at-home' order; will take effect close of business Tuesday

Vaccine-inuced immunity provides more robust heterotypic immunity than natural infection to emerging SARS-CoV-2 variants of concern

South Carolina Public Health Officials Detect Nation's First Known Cases of the COVID-19 Variant Originally Detected in South Africa

Understanding spatial heterogeneity of COVID-19 pandemic using shape analysis of growth rate curves

Experts call it a 'pandemic of the unvaccinated'. National Public Radio

Real-time monitoring the transmission potential of COVID-19 in Singapore

Prevention and Attenuation of covid-19 with the BNT162b2 and mRNA-1273 vaccines

A network-based explanation of why most covid-19 infection curves are linear

National Demographic Analysis Tables: 2020

Coronavirus (COVID-19) Update: FDA Issues Policies to Guide Medical Product Developers Addressing Virus Variants

United States Food and Drug Administration (2020)

FDA Briefing Document Pfizer-BioNTech COVID-19 Vaccine

FDA Issues Emergency Use Authorization for Third COVID-19 Vaccine

FDA Takes Key Action in Fight Against COVID-19 By Issuing Emergency Use Authorization for First COVID-19 Vaccine

United States Food and Drug Administration

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Transmission of SARS-CoV-2 Lineage B.1.1.7 in England: Insights from linking epidemiological and genetic data

Florida governor extends order suspending COVID-19 related enforcement fines

Executive Order No. 2020-42. State of Michigan Office of the Governor

Executive Order No. 2020-60. State of Michigan Office of the Governor

Executive Order No. 2020-110. State of Michigan Office of the Governor

WHO Diretor-General's opening remakrs at the mission briefing on COVID-19

mRNA-1273 vaccine induces neutralizing antibodies against spike mutants from global SARS-CoV-2 variants. biorxiv

Neutralization of SARS-CoV-2 spike 69/70 deletion, E484K, and N501Y variants by BNT162b2 vaccine-elicited sera

A data-driven network model for the emerging covid-19 epidemics in Wuhan

A novel coronavirus from patients with pneumonia in China

One of the authors (ABG) acknowledge the support, in part, of the Simons Foundation (Award #585022) and the National Science Foundation (grant number: DMS-2052363). The authors are grateful to the anonymous reviewers for their constructive comments.

☒ 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.☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:J o u r n a l P r e -p r o o f