ABSTRACT
Multiple effective vaccines are currently being deployed to combat the COVID-19 pandemic (caused by SARS-COV-2), and are viewed as the major factor in marked reductions in disease burden in regions around the world with moderate to high coverage of these vaccines. 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 and potentially more virulent than the wild (resident) strain, may at least partially evade existing vaccines. A new two-strain (one resident, the other wild) 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. Analysis of the model, which is fitted using COVID-19 mortality data for the US, shows that vaccine-induced herd immunity can be achieved if 61\% of the American population is fully vaccinated with the Pfizer or Moderna vaccines. Parameter sensitivity analysis suggests three main factors that significantly affect the COVID-19 burden in the US, namely (a) daily vaccination rate, (b) the level of cross-protection the vaccines offer against the variant, and (c) the relative infectiousness of the dominant variant relative to the wild strain. This study further suggests that a new variant can cause a significant disease surge in the US if (i) the vaccine coverage against the wild strain is low (roughly $<50\%$), (ii) the variant is much more transmissible (e.g., twice 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 70\%). A new variant will not cause such surge in the US if it is only moderately more transmissible (e.g., $1.56\%$ more transmissible) than the wild 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.
1. Introduction
Since its emergence late in December 2019, the novel coronavirus pandemic (COVID-19), 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 (50)), COVID-19 represents the greatest global public health challenge since the 1918-1919 influenza pandemic. As of June 18, 2021, there have been over 176 million confirmed COVID-19 cases and over 3.8 million deaths worldwide, with the US the most affected county (with over 33 million confirmed cases and now over 600,000 COVID-19 deaths), although cases and deaths may be dramatically under-counted in some communities or jurisdictions (17).
The most common symptoms of COVID-19 disease include dry cough, shortness of breath, fever, fatigue, muscle aches, and loss of taste or smell. Symptoms develop roughly 2-14 days after contacting the SARS-CoV-2 virus, likely via respiratory droplets or aerosols, and typically last for up to two weeks, although a minority experience “long COVID-19,” when symptoms endure for several weeks or months, and interfere with day-to-day life (13). The severity of symptoms varies widely from wholly asymptomatic disease to severe multi-organ system failure and death. Risk of a poor outcome increases with comorbidities, such as immunocompromised status, heart or lung disease, diabetes, etc., while the mortality rate increases exponentially with age. Transmission to susceptible individuals peaks early in the natural history of the disease, with infected individuals infectious several days before the onset of symptoms (pre-symptomatic stage); those who never develop symptoms may also still pass the disease as asymptomatic carriers (14).
Until late December 2020, control and mitigation measures against COVID-19 in the US were exclusively non-pharmaceutical interventions (NPIs), such as the quarantine of those suspected of being exposed to the disease, isolation of individuals with disease symptoms, community “lockdowns,” social (physical)-distancing, and the use of face masks (24; 49; 39; 26; 32; 33). Numerous states and jurisdictions within the US have implemented lockdown and mask mandates of varying strictness and duration to limit COVID-19 spread (23; 51; 28; 67; 68; 69; 29; 41; 1; 2).
The US Food and Drugs Administration (FDA) gave emergency use authorization (EUA) for the first vaccine against COVID-19 on December 11, 2020 (60; 62), developed by Pfizer-BioNtech and based on delivery of RNA (mRNA) encoding the SARS-CoV-2 spike protein (44). Two doses of the vaccine administered 21 days apart showed 95% efficacy against symptomatic COVID-19 in clinical trials including 44,000 participants (62; 60). A week later (December 18, 2020), the FDA authorized a vaccine developed by Moderna Inc. (59). This vaccine is also based on mRNA technology. Delivered in two doses four weeks apart, the vaccine showed efficacy of 94.5% and in initial clinical trials (59). The most recent vaccine to receive EUA (February 27, 2021) is the Janssen vaccine, developed by Johnson & Johnson and administered in a single dose (61). Unlike the mRNA-based Pfizer and Moderna vaccines, the Johnson & Johnson vaccine was developed an using adenovirus vector encoding the SARS-CoV-2 spike protein. Results from clinical trials have shown the Johnson & Johnson vaccine to be 67% effective in preventing moderate to severe COVID-19 14 days post-vaccination (61). Although the vaccines have generally only been available to those 16 and older (Pfizer) or 18 and older (Moderna, J&J), but the FDA recently extended EUA for the Pfizer vaccine to be administered to children of ages 12–15 (63).
Following EUA in December 2020, limited vaccine doses were made available according to prioritization schedules that varied somewhat by state, but generally prioritized the elderly, healthcare workers, and long-term care residents (10). Vaccinations rates generally increased nearly linearly since then for several months, but have since peaked and are on the decline as of June 2021. Even with all individuals aged 12 and above now eligible to receive the Pfizer vaccine in the US, increasing vaccine coverage is now widely viewed as a problem of limited demand rather than supply. While vaccines are highly effective against the original SARS-CoV-2 strain, sustainable control of the COVID-19 epidemic may hinge on the characteristics of emerging variant strains.
Viruses are subject to mutations when replicating, and mutations that increase the fitness of the virus tend to persist in the population. In fact, by the middle of March 2020, the dominant SARS-CoV-2 strain in the US has already mutated from the original strain that emerged in Wuhan, China three months earlier (42). Since then, several additional variants from the original SARS-CoV-2 strain have appeared around the globe (12) with many likely more transmissible than the original SARS-CoV-2 strain (20; 65), and possibly more deadly. For instance, the B.1.1.7 variant first appeared in the United Kingdom during late summer of 2020 (65), and may be associated with an increased risk of disease-induced death (12). The mutations of this variant were found on the spike protein of the SARS-CoV-2 virus (15). This variant was first detected in the US on December 29, 2020 in Colorado (18). The B.1.1.7. variant was shown to be 56% more transmissible than the wild-type (original) SARS-CoV-2 strain (20). Individuals with the B.1.1.7. variant were 1.4 times more likely to be hospitalized, and 1.58-1.71 times more likely to die from COVID-19 complications (38). Another variant, B.1.351, emerged independently from the B.1.1.7. variant, but shares some similar mutations (12). The B.1.351 variant was first identified in South Africa in October 2020, and was detected in the US at the end of January 2021 (53). Furthermore, the P.1. variant of the SARS-CoV-2 virus has 17 unique mutations, with origins tracing back to Brazil. This variant also appeared in the US at the end of January 2021 (47). Finally, the most recent VOC to have emerged in the US is the B.1.617.2 that first appeared in India.
Other SARS-CoV-2 variants have appeared in RNA sequencing, but the B.1.1.7, P.1, and B.1.351 variants represent the greatest proportion of variants currently circulating in the US and are currently classified as variants of concern (VOC) (16; 15). VOCs have shown evidence of (i) more disease rapid spread, (ii) more severe disease, (iii) evading current diagnostic tests, and/or (iv) resistance to naturally-acquired or vaccine-induced immunity when compared to the wild-type SARS-CoV-2 strain, which does not contain any major mutations (12; 15; 70). These concerns are compounded by the fact that, as more individuals are vaccinated, mask mandates and social-distancing measures are being relaxed, or fully-lifted altogether (3; 22; 66). 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 variant B.1.617.2 (16). It can be reasonably assumed that the reported cases of COVID-19 caused by variants underestimates the number of true cases. The B.1.617.2 variant, for example, is estimated to be responsible for 6% of US cases, and as much as 18% for some states (31). This is due to the limited genomic sequencing capabilities and strain surveillance, when compared to the magnitude of infection transmission. Alhough several government-operated, university affiliated, and independent laboratories throughout the US sequence several thousand strains of the SARS-CoV-2 virus each week, this still amounts to less than 1% of total infections reported (12).
Recent studies regarding the cross-protective immunity of existing vaccines on the SARS-CoV-2 VOCs have shown varying results (52; 70; 30; 71). 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 (52; 71). However, for those who acquired immunity through infection recovery (convalescent immunity) or only received a single dose of the vaccine, less neutralizing antibodies were produced, particularly for the B.1.351 variant (52; 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 (70), 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 (58; 12). However, it is of major public health interest to understand how new variants can, in principle, impact the community, and if the current vaccination program will be enough to eliminate the pandemic in the United States, under a plausible range of scenarios with respect to variant transmissibility and cross-protection.
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 (35; 37; 43; 54; 55; 27; 72; 56; 26; 19; 40; 24; 39; 49; 32; 33). 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 cross-protective immunity from vaccines protect the community from the emerging COVID-19 variants. The current study uses a novel mathematical model to assess the impact of vaccination and vaccine-induced cross-protection 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 was estimated from the fitting. Further, basic qualitative features of the model, including the derivation of vaccine-derived herd immunity threshold for the US and parameter sensitivity analysis, are discussed in Section 3.1. Numerical simulations are presented in Section 3.4.
2. Methods
2.1. Model Formulation
In this section, the new two-strain, two-group mechanistic model for the transmission dynamics of two SARS-CoV-2 strains (the wild-type and variant strains) in the presence of vaccination, is formulated. The total population of individuals in the community at time t, denoted by N(t), is stratified into two groups based on vaccination status, namely the total unvaccinated group, denoted by NU (t) and the total vaccinated group denoted by NV (t), so that N(t) = NU (t) + NV (t). The total unvaccinated population (NU (t)) is further stratified into the mutually-exclusive compartments of individuals that are unvaccinated susceptible (SU), latent or “exposed” (Ei), pre-symptomatically infectious (Pi), symptomatically infectious (Ii), asymptomatically infectious (Ai), hospitalized (Hi), and recovered (Ri), with i = 1, 2 (represents the SARS-CoV-2 strain), so that Similarly, the total vaccinated population (NV (t)) is sub-divided into the sub-populations of individuals that are vaccinated susceptible (SV (t)), exposed (EV i), pre-symptomatically infectious (PV i), symptomatically infectious (V i), asymptomatically infectious (AV i), hospitalized (HV i), and recovered (RVi), with i = 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 in 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: strain 1, the wild-type (original) strain that is assumed to be the predominant strain circulating the US throughout 2020 (42) and strain 2, a new dominant, variant strain assumed to be circulating in the US since the end of 2020 (18)). For mathematical tractability and convenience, we assume that susceptible individuals can be infected with either of the two strains, but not with both. Further, we have no 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 εv, and are assumed to induce cross-protective efficacy (εc) against the dominant variant (strain 2) (52; 71; 30; 70). 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 the following:
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 (36).
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 (60)), we include demographic (birth and death) parameters to account for the inflow of new susceptible individuals that are eligible for vaccination.
Vaccination is only offered to wholly-susceptible individuals or those who naturally recovered from COVID-19 infection (but not for currently-infected individuals).
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.
The vaccines are imperfect against infection (i.e, breakthrough infection can occur) (60; 59; 61). The vaccines offer strong therapeutic benefits, vis a vis reducing severe disease, hospitalization, and death (60; 59; 61).
Vaccine-derived and natural immunity may wane over time in individuals, in which they revert to the wholly-susceptible class.
2.2. Data Fitting and Parameter Estimation
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) (17) and daily vaccination data from Bloomberg (from December 21, 2020 to March 6, 2021) (8) for the US. 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 (17; 8). 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) (14).
The first case of the B.1.1.7 variant was documented in the US on December 29, 2020 (18), although the variant was likely introduced to the US prior to detection. For the actual fitting, we introduce 1,000 cases of a variant strain (i.e.,I2(0) = 1, 000) on December 29, 2020 to estimate variant infections present prior to confirmation of its US detection. Further, vaccination is introduced in the model (A-8) on January 21, 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 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 (8).
Furthermore, the following expression was used to determine the time-dependent vaccination rate, αv(t): 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. The left panel of this figure shows the model (A-8) fit using the observed cumulative US mortality data. The right panel shows the simulations of the model (A-8) with the fixed and fitted parameters, compared to the observed daily mortality data for the US. The purple vertical dotted line in each panel indicates the time that a variant strain was introduced into the US population (December 29, 2020). The contact rate for the variant strain, (β2) was assumed to be 1.56 times greater than the contact rate for the wild-type strain (β1) (20). The brown dotted lines indicate the time when vaccination was introduced into the population (January 21, 2021). The model (A-8) with the chosen fixed parameter values fits well with the observed cumulative and daily US COVID-19 death data.
The contact rate, β1, varies throughout the course of the pandemic due to implementation and compliance of non-pharmaceutical interventions (NPIs) and vaccination. 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 (66), 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.
3. Results
3.1. Mathematical Analysis: Computation of Reproduction Number
In this section, the reproduction number of the model (A-8) is computed for the special case where the vaccination ate, αv(t), is treated as constant. That is, We consider the autonomous version of the model with αv(t) = αv, a constant. The autonomous version of the model (A-8) has a unique disease-free equilibrium (DFE) given by where and all other components of 𝔼0 are zero. The next generation operator method (64; 21) can be used to analyze the local asymptotic stability of the DFE (𝔼0). Noting that , and using the notation in (64), it follows that the non-negative matrix (F) of new infection terms, and the matrix (M) of new transition terms, are, respectively, given by: where the matrices F1, F2, M1 and M2 are given in Appendix B, and 010×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 ℛvac, is given by (where ρ denotes the spectral radius of the next generation matrix FM−1) where and (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 (64).
The DFE (𝔼0) of the autonomous version of the model (A-8) is locally-asymptotically stable if ℛvac < 1, and unstable if ℛvac > 1. The threshold quantity, ℛvac, 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 (ℛvac) is brought to (and maintained at) a value less than unity.
3.2. Computation of the Herd Immunity Threshold
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 (6; 5). Furthermore, at least 15% of the U.S. population currently cannot be vaccinated due to age restrictions (57). 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 herd immunity (so that the pandemic can be effectively controlled). In this section, we derive the expression for achieving 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 (ℛvac) 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 ki, with i = 1, 2, ., 20 are given in Appendix B): For the case when all individuals in the population are vaccinated, the constituent reproduction numbers for the wild-type (strain 1) and variant (strain 2) SARS-CoV-2 strains become: where expressions for Qi, i = 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, ) and variant (strain 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 vaccinated (at steady-state). The herd immunity threshold is computed by setting and solving for fv (noting that it is necessary to ensure, in this computation, that ). Doing so gives: Using the fixed parameters in Table A.9 and a contact rate of β1 = 0.24 to account for transmission in the absence of nonpharmaceutical control measures (face mask usage, social distancing policies), the herd immunity threshold for the wild-type (strain 1) SARS-CoV-2 strain is 0.61 for the Pfizer or Moderna vaccines (εv = 0.94), and 0.86 using the Johnson & Johnson vaccine (εv = 0.67). Thus, 61% of community members must be vaccinated if all individuals use the Pfizer or Moderna vaccines, or 86% of community members if the Johnson & Johnson vaccine is used.
3.3. Parameter Sensitivity Analysis
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 dynamics of variants) used in the numerical simulations. In this section, sensitivity analysis is carried out, using Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) (45; 9; 46), 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 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 . The baseline value of the contact rate, β2 was chosen as 0.24, with additional baseline values for fixed parameters from Table A.9. The PRCC values generated for the parameters in the expression of the response function, , are tabulated in Table 2. The contact rate of the variant strain (β2), the fraction of non-vaccinated and vaccinated individuals who experience symptomatic infection from the variant strain (q, qv), and the modification parameter (ηA2) accounting for relative infectiousness of non-vaccinated asymptomatic individuals infected with the variant strain compared to non-vaccinated asymptomatic individuals infected with the wild-type strain have the greatest PRCC magnitude. That is, β2, q, qv, and ηA2 most affect the vaccination reproduction number of the variant strain , and, thus, the dynamics of the model (A-8).
3.4. Numerical Simulations: Assessment of Vaccination Program
The model (A-8) is now simulated using the fixed parameters in Table A.9 to assess the population-level impact of the emerging variant on the dynamics of COVID-19. Simulations for the US are carried out for the period March 6,2021 to March 6, 2023 for two levels of variant strain transmissibility in relation to the wild-type strain: Scenario 1 assumes the variant strain is 1.56 times transmissible (i.e., B.1.1.7 variant (42)), and Scenario 2 assumes the variant strain is twice as transmissible. The values of each model compartment on the last point of the fitting period on March 6, 2021 (see Section 3) were used as initial conditions to simulate projections of future COVID-19 dynamics in the U.S. The contact rate for the wild-type strain was chosen to be 0.0805 from March 6, 2021 to March 30, 2021 (i.e., the result of fitting the post-third wave period, see Table A.9) and increased to 0.161 on April 1, 2021 to account for the alleviation of face mask mandates and social distancing policies in various states (22; 3). Each Scenario is assessed for three levels of vaccination coverage: 50%, 66%, and 75% (i.e., fv = 0.5, 0.66, and 0.75). Simulations are carried out assuming vaccinated individuals receive both doses of the Pfizer/Moderna vaccine (εv = 0.94).
Due to inconclusive data surrounding the length of vaccine and naturally-acquired immunity, we assume for the period of simulation that previously infected individuals remain in their respective recovered classes. This assumption is consistent with other COVID-19 modeling studies (32; 33; 49; 48; 39; 24). Thus, results using administration of the Johnson & Johnson vaccine are not shown, since the assumption that individuals will not become re-infected during the simulation period results in an overall lower COVID-19 compared to administering the Pfizer/Moderna vaccines. In other words, under administration of the less efficacious Johnson & Johnson vaccine (εv = 0.67), more individuals will become infected with the less deadly wild-type strain rather than the more deadly variant strain. Furthermore, since less than 8% of fully-vaccinated individuals in the US have received the Johnson & Johnson vaccine (11), simulations under Pfizer/Moderna vaccine administration can better estimate the future course of the pandemic.
3.4.1. Scenario 1: Variant is 1.56 Times More Transmissible Than Wild-type Strain
Data from (42) suggests that the B.1.1.7. variant is 1.56 times more transmissible (i.e., 1.56 times more infectious) than the wild-type strain in the US. Here, we simulate the dynamics of the two strains under this assumption, for various levels of full-vaccination coverage (50%, 66%, 75%) and cross-protective vaccination efficacy (εc = 0, 0.25, 0.5, 0.7) offered by the Pfizer/Moderna vaccines. Cumulative COVID-19 related deaths and cases for Scenario 1 are tabulated in Table 3, and values and occurrence of the peak daily COVID-19 related deaths and new infections are included in Table 4.
The absolute worst-case-scenario for the population would be if the vaccines do not offer any cross-protection (i.e, if εc = 0). Fortunately, the three EUA vaccines currently used in the US do offer a significant level cross-protective efficacy against existing variants (58; 12). However, it is helpful to quantify how much morbidity and mortality burden from the variant strain could be averted through cross-protection. Table 3 shows that, under the worst-case scenario, nearly half of the US population would have acquired infection, causing over two million deaths. Roughly 73% of cases and deaths would be caused by the variant strain. The peak daily burden would be experienced during mid September 2021, with over 2 million new cases per day.
If 50% of the US population became fully-vaccinated, then as much as 192,000 to 918,000 lives can be saved under 25% to 70% cross-protective vaccination efficacy. If vaccination coverage is increased to 75%, then 3 to 14.8 million deaths can be averted. At a moderately-high level of cross-protective efficacy (εc = 0.7), then even vaccination coverage levels as low as 66% can avoid another pandemic wave. Figure 4 shows the simulations of Scenario 1 when 75% of the US population is fully vaccinated with the Pfizer/Moderna vaccines.
3.4.2. Scenario 2: Variant is Twice as Transmissible as Wild-type Strain
If we assume now that the variant is, instead, twice as transmissible as the wild-type strain in the US, the worst-case-scenario (εc = 0) would cause an additional 214,000 cumulative deaths and 17.6 million cases than the worst case from Scenario 1. The peak of daily cases and deaths would shift 9 weeks earlier to late July 2021 and exceed 3.76 million new cases a day. If 50% of the US population is fully vaccinated, 23.5 to 35.5 million more cases and 300,000 to 432,000 more deaths may occur when compared to the same vaccination coverage under Scenario 1. An additional 17.6 to 44.1 million cumulative cases and 214,000 to 537,000 cumulative deaths will occur under 75% vaccination coverage when compared to Scenario 1. With a more infectious variant strain, even high levels of vaccination coverage (75%) and cross-protective vaccination efficacy (εc = 0.7) will result in a future pandemic wave. Figure 5 shows the simulations of Scenario 2 when 75% of the US population is fully vaccinated with the Pfizer/Moderna vaccines. Cumulative COVID-19 related deaths and cases for Scenario 2 are tabulated in Table 5, and values and occurrence of the peak daily COVID-19 related deaths and new infections are included in Table 6.
4. Discussion & Conclusions
Since the emergence of COVID-19 in Wuhan, China late 2019, the COVID-19 pandemic has ravaged the world by infecting every country and altered life as we know it. Fortunately, scientists and researchers have swiftly developed several vaccines to combat the virus using mRNA (Pfizer and Moderna) and adenovirus (Johnson & Johnson) technology. Over 2.23 billion shots have been administered worldwide as of June 18, 2021 and over 43% of the US population is fully vaccinated (7; 11). However, recent variants of SARS-CoV-2 have emerged and become a public health concern, particularly: the B.1.1.7 variant first discovered in the UK, the B.1.351 variant initially found in South Africa, the P.1 variant that originated from Brazil, and the B.1.617.2 variant from India have infected numerous countries, including the US. Current data and studies show strong evidence that these variants are more easily transmissible and have higher rates of hospitalization and disease-induced death. The greatest concern of the variants is the inconclusive cross-protective efficacy from the current COVID-19 vaccines available.
Numerous studies have investigated the impact of COVID-19 on different populations through mathematical modeling techniques. These studies include statistical (40; 55; 54), mechanistic/deterministic (48; 25; 32; 33), stochastic (35; 37; 43), network (27; 56; 72), and agent based (19; 26) modeling types. This study presents a novel deterministic system of 26 nonlinear differential equations to assess 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 two-strain model (A-8) includes vaccination of susceptible and naturally recovered individuals, and its autonomous equivalent is shown to have an asymptotically stable disease free equilibrium (DFE) when the vaccination reproduction number (ℛvac) is less than unity. The vaccination reproduction number measures the total amount of new infections that result from a single infection in a partially protected population, and an explicit expression for ℛvac is derived, along with the expression for the herd immunity threshold . The herd immunity threshold measures the fraction of susceptible community members that must acquire immunity through vaccination, in the U.S. population to halt the COVID-19 pandemic. Additionally, the parameters most sensitive to the vaccination reproduction number of the variant strain were found to be the variant strain contact rate (β2), the proportion of individuals (non-vaccinated and vaccinated) who experience symptoms from infection with the variant strain (q, qv), and the modification parameter (ηA2) that accounts for relative infectiousness of non-vaccinated asymptomatic individuals infected with the variant strain compared to non-vaccinated asymptomatic individuals infected with the wild-type strain. Finally, the model (A-8) was numerically simulated for four different Scenarios with varying vaccine efficacy (εv) of the wild-type strain, cross-protective vaccination efficacy (εc) of the variant strain, successful vaccination coverage level (fv), and relative transmission rate of the variant strain compared to the wild-type strain.
The worst case occurs when there is no level of cross-protective vaccination efficacy against the variant (i.e., if εc = 0). The worst case predicts that 49.85% of the U.S. population will have acquired infection by the March 6, 2023, over 2 million individuals would have lost their lives from COVID-19, and 72.5% of cases and COVID-19 related deaths resulting from a variant strain that is 1.56 times as infectious as the wild-type strain. If the variant ends up being twice as infectious as the wild-type strain, then as much as 55% of the U.S. population will have been infected with COVID-19 by March 6, 2023, with now more than 2.2 million COVID-19 related deaths, and 75.4% of the pandemic burden being caused by the variant strain.
Wide-scale vaccination is essential to eliminating the pandemic and reducing the overall COVID-19 burden. The vaccination coverage level, cross-protective vaccination efficacy against variant strains, and the relative infectiousness of the variant strain when compared to the wild-type strain all greatly influence the dynamics of the model (A-8). Both cross-protective vaccination efficacy against the variant and vaccination coverage shift the peak of daily cases and deaths forward in time with increasing levels. The peak of daily deaths is expected to occur approximately 10 to 16 days after the peak of daily infections occurs.
No data is currently available about if, or how much infectiousness is reduced by vaccinated individuals. Therefore, this study assumes that vaccinated individuals, when infected, are equally likely to transmit infection as non-vaccinated individuals (θv = θc = 1). This study also assumes that exclusive deployment of the Pfizer or Moderna vaccines in the simulations. In reality, there’s a combination of individuals receiving the Johnson & Johnson versus the Pfizer or Moderna vaccines, with each possibly having a different level of cross-protective efficacy against variant strains and therapeutic effects (e.g., increased recovery rates, decreased disease-related death rates). Moreover, another limitation of the current framework of the model (A-8) is that it can only analyze one variant strain of SARS-CoV-2. Extension of this model could include multiple variant strains to determine the individual and collective impact on the COVID-19 pandemic. Since current studies show that naturally-acquired immunity is less likely to produce an immune response to the existing variants than vaccination (52; 30), another future extension of the model (A-8) could include the possibility that individuals who recovered from the wild-type strain (R1 class) could be directly infected by the variant strain without losing developed immunity from the wild-type strain. This would alter the vaccination reproduction number and basic reproduction number expressions of the model (A-8).
The B.1.1.7 is currently the most dominant variant circulating in the US population (12; 16), and estimated to be 56% more infectious than the wild-type strain (20). The numerical simulations showed that if a moderate level of vaccination coverage (≥ 66% of the US population) and moderately-high level of cross-protective vaccination efficacy (εc ≥ 0.7) is attained, then a surge of variant-induced COVID-19 infections and mortality will not occur (under the assumption that both vaccine and naturally-acquired immunity lasts at least two years). One study showed that cross-protective efficacy is as high as 93% against the B.1.1.7 variant when two doses of the Pfizer vaccine are administered, but drops to 33% when only one dose is received (31). If the variant is much more transmissible than the wild-type strain, then even at moderately high levels of cross-protective efficacy (εc = 0.7) and vaccine coverage (75% of the US population), the U.S. could still experience a severe wave of COVID-19 cases caused by variant strains that will be detrimental to the community. Numerical simulations show that if the US experiences another wave of infections, the peak daily COVID-19 related deaths and new cases are almost exclusively from the variant strain. This could result in another period of lockdown and exacerbate morbidity and mortality from the virus. Despite many states relieving COVID-19 related restrictions and face mask mandates now that vaccination has become more widely available (22; 3), it imperative that we continue to follow non-pharmaceutical intervention (NPI) protocol. Further, the risk of the B.1.617.2 variant that recently emerged in the US has had devastating consequences in South Asia, and we are wait to see how the transmission of this variant unfolds in the US. Current studies show that the full recommended dosage of the Pfizer vaccine is 88% against the B.1.617.2 variant (31). Thus, it is critical to encourage all eligible individuals to receive the full two dose regimen of the Pfizer or Moderna vaccines.
Data Availability
Data used for simulations are publicly available at Johns Hopkins University COVID-19 repository
Declaration of Competing Interests
The authors declare no competing interests.
Conflict of Interest
Declaration of interests
☒ 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:
Appendix A Equations of the Mathematical Model
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:
Description of Parameters of the Model
The parameter ∏ is the recruitment rate (via birth or immigration) of individuals into the population. Individuals are vaccinated at a time-dependent rate αv(t). The vaccine is assumed to produce protective efficacy 0 < εv < 1 against strain 1, potential cross-protective efficacy 0 ≤ εc < 1 against strain 2, and wane at rate ωv. Natural death occurs in all epidemiological classes at rate μ. Susceptible individuals acquire infection with strain 1 at a rate λ1 or with strain 2 at a rate, λ2, where
and, In (A-9) and (A-10), βi (with i = 1, 2) is the effective contact rate for strain i, ηPi (ηPvi), ηAi (ηAvi), and ηHi (t7Hvi) 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 i, compared to symptomatic infectious individuals with strain i. The modification parameter 0 ≤ θv ≤1 represents the assumed reduced infectiousness of individuals vaccinated against strain 1 (in comparison to unvaccinated infectious individuals with strain 1), while (0 ≤ θc ≤1) accounts for the assumed reduced infectiousness of strain 2 induced by the cross-protective 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 σi (σvi) where i = {1, 2} indicates the strain of infection. Non-vaccinated (vaccinated) individuals transition out of the pre-symptomatic classes at rate, ξi (,ξvi). For strain 1, fraction 0 < r (rv) < 1 of individuals develop symptoms while (1 − r) ((1 − rv)) individuals remain asymptomatic infectious. The fraction of un-vaccinated (vaccinated) individuals infected with strain 2 who develop symptoms is 0 < q (qv) < 1. Non-vaccinated (vaccinated) symptomatic infectious individuals recover at rate γIi (γIvi), become hospitalized at rate ϕi (ϕvi), and die from diseased-induced death at rate δIi (δIvi). Non-vaccinated (vaccinated) asymptomatic infectious individuals recover at rate γAi (γAvi). Non-vaccinated (vaccinated) Hospitalized individuals recover at rate γHi (γHvi) and die from disease-induced death rate δHi (δHvi). Non-vaccinated (vaccinated) individuals lose their naturally-acquired disease-induced immunity at rate Ψi (Ψvi).
Table of State Variables and Parameters of the Model
The state variables and parameters of the model are described in Tables A.7 and A.8, respectively.
Values of Fixed Parameters
Appendix B Computation of ℛVac
The associated next generation matrices of the model (A-8) are given by: with entries
and with entries with entries and, with entries Using the approach in (34), 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 and are given by: Where,
Footnotes
Suggested Reviewers: Enahoro Iboi, PhD, Assistant Professor, Spelman College, enahoroiboi{at}spelman.edu, Has extensively worked on infectious disease modeling for COVID-19, Calistus Ngonghala, PhD, Assistant Professor, University of Florida, calistusnn{at}ufl.edu, Has extensively worked on infectious disease modeling for COVID-19
References
- [1].↵
- [2].↵
- [3].↵
- [4].↵
- [5].↵
- [6].↵
- [7].↵
- [8].↵
- [9].↵
- [10].↵
- [11].↵
- [12].↵
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].↵
- [29].↵
- [30].↵
- [31].↵
- [32].↵
- [33].↵
- [34].↵
- [35].↵
- [36].↵
- [37].↵
- [38].↵
- [39].↵
- [40].↵
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].↵
- [47].↵
- [48].↵
- [49].↵
- [50].↵
- [51].↵
- [52].↵
- [53].↵
- [54].↵
- [55].↵
- [56].↵
- [57].↵
- [58].↵
- [59].↵
- [60].↵
- [61].↵
- [62].↵
- [63].↵
- [64].↵
- [65].↵
- [66].↵
- [67].↵
- [68].↵
- [69].↵
- [70].↵
- [71].↵
- [72].↵