Toward Achieving a Vaccine-Derived Herd Immunity Threshold for COVID-19 in the U.S.

A novel coronavirus emerged in December of 2019 (COVID-19), causing a pandemic that inflicted unprecedented public health and economic burden in all nooks and corners of the world. Although the control of COVID-19 largely focused on the use of basic public health measures (primarily based on using non-pharmaceutical interventions, such as quarantine, isolation, social-distancing, face mask usage, and community lockdowns) initially, three safe and highly-effective vaccines (by AstraZeneca Inc., Moderna Inc., and Pfizer Inc.), were approved for use in humans in December 2020. We present a new mathematical model for assessing the population-level impact of these vaccines on curtailing the burden of COVID-19. The model stratifies the total population into two subgroups, based on whether or not they habitually wear face mask in public. The resulting multigroup model, which takes the form of a deterministic system of nonlinear differential equations, is fitted and parameterized using COVID-19 cumulative mortality data for the third wave of the COVID-19 pandemic in the United States. Conditions for the asymptotic stability of the associated disease-free equilibrium, as well as an expression for the vaccine-derived herd immunity threshold, are rigorously derived. Numerical simulations of the model show that the size of the initial proportion of individuals in the mask-wearing group, together with positive change in behavior from the non-mask wearing group (as well as those in the mask-wearing group, who do not abandon their mask-wearing habit) play a crucial role in effectively curtailing the COVID-19 pandemic in the United States. This study further shows that the prospect of achieving vaccine-derived herd immunity (required for COVID-19 elimination) in the U.S., using the Pfizer or Moderna vaccine, is quite promising. In particular, our study shows that herd immunity can be achieved in the U.S. if at least 60% of the population are fully vaccinated. Furthermore, the prospect of eliminating the pandemic in the U.S. in the year 2021 is significantly enhanced if the vaccination program is complemented with non-pharmaceutical interventions at moderate increased levels of compliance (in relation to their baseline compliance). The study further suggests that, while the waning of natural and vaccine-derived immunity against COVID-19 induces only a marginal increase in the burden and projected time-to-elimination of the pandemic, adding the impacts of therapeutic benefits of the vaccines into the model resulted in a dramatic reduction in the burden and time-to-elimination of the pandemic.

1. Introduction

The novel coronavirus (COVID-19) pandemic, which started as a pneumonia of an unknown etiology late in December 2019 in the city of Wuhan, became the most devastating public health challenge mankind has faced since the 1918/1919 influenza pandemic. The COVID-19 pandemic, which rapidly spread to essentially every nook and corner of the planet, inflicted devastating public health and economic challenges globally. As of January 24, 2021, the pandemic accounted for about 100 million confirmed cases and 2, 128, 721 cumulative mortality globally. Similarly, as of this date, the U.S., which reported its first COVID-19 case on January 20, 2020, recorded over 25, 123, 857 confirmed cumulative cases and 419, 204 deaths (1).

COVID-19, a member of the coronavirus family of RNA viruses, is primarily transmitted from human-to-human through inhalation of respiratory droplets from both symptomatic and asymptomatically-infectious humans (2) albeit there is limited evidence that COVID-19 can be transmitted via exhalation through normal breathing and aerosol (3). The incubation period of the disease is estimated to be between 2 and 14 days (with a mean of 5.1 days), and majority of individuals infected with the disease show mild or no clinical symptoms (4). The symptoms typically include coughing, fever and shortness of breadth (for mild cases) and pneumonia for severe cases (4). The people most at risk of dying from, or suffering severe illness from, COVID-19 are those with co-morbidities (such as individuals with diabetes, obesity, kidney disease, cardiovascular disease, chronic respiratory disease, etc.). Younger people, frontline healthcare workers and employees who maintain close contacts (within 6 feet) with customers and other co-workers (such as meat factory workers, retail store workers, etc.) are also at risk.

Prior to the approval of the three safe and effective vaccines (by AstraZeneca, Moderna, and Pfizer) for use in humans in December 2020 (5, 6), the control and mitigation efforts against COVID-19 have been focused on the use of non-pharmaceutical interventions (NPIs), such as quarantine, self-isolation, social (physical) distancing, the use of face masks in public, hand washing (with approved sanitizers), community lockdowns, testing, and contact tracing. Of these NPIs, the use of face masks in public was considered to be the main mechanism for effectively curtailing COVID-19 (4, 7–9). Furthermore, owing to its limited supply, the approved anti-COVID drug remdesivir is reserved for use to treat individuals in hospital who display severe symptoms of COVID-19. The U.S. started administering the Pfizer and Moderna vaccines by December 2020 (5, 6). Another vaccine by Janssen Biotech Inc., the vaccine division of Johnson & Johnson, received Emergency Use Authorization (EUA) in the U.S. in late February 2021 (10).

The Pfizer and Moderna vaccines, each offering a protective efficacy of about 95% (11–13), are genetic vaccines that trigger the immune system to recognize the coronavirus' spike protein and develop antibodies against it (11, 14). Each of these vaccines is administered in a two-dose structure (one dose to prime the immune system, and the second to boost it). For the Pfizer vaccine, the second dose is administered 19–42 days after the first dose, while that for the Moderna vaccine is administered 3–4 weeks after the first dose. Both vaccines need to be stored at appropriate refrigeration temperatures (15). The AstraZeneca vaccine, on the other hand, has estimated protective efficacy of 70% (11–13). It uses a replication-deficient chimpanzee viral vector that contains the genetic material of the SARS-CoV-2 virus spike protein (13). The AstraZeneca vaccine also requires two doses (1 month apart) to achieve immunity, and unlike the Pfizer and Moderna vaccines, does not have to be stored in super-cold temperatures (13). The Johnson & Johnson vaccine is an adenovirus single-dose vaccine. The efficacy of the Johnson & Johnson vaccine for preventing severe disease is 85% (16, 17). Like in the case of AstraZeneca vaccine, the Johnson & Johnson vaccine does not require extremely cold temperatures for storage.

An effective vaccine typically offers a range of protective and therapeutic benefits to the vaccinated individual, such as reducing the risk of acquiring infection and reducing the severity of disease, hospitalization, and mortality and accelerating recovery in breakthrough infections (for vaccines that offer strong therapeutic benefits) (18, 19). A vaccine that has protective and therapeutic efficacies, when introduced into a population during an epidemic, will play a major role in curtailing the epidemic, and its effective deployment would be dependent on the rollout strategy used. The goal of this study is to design a structured mathematical model that will allow for the realistic assessment of the population-level impact of vaccination programs based on using three of the aforementioned four COVID-19 vaccines (namely the AstraZeneca, Moderna, and Pfizer vaccines), with emphasis on determining the optimal coverage rate needed to achieve vaccine-derived herd immunity (which is required for eliminating the pandemic). A secondary objective of this study is to explore whether the prospect for eliminating the pandemic in the U.S. will be enhanced if the vaccination program is combined with NPIs, such as social-distancing at some level of compliance.

Numerous mathematical models, of various types, have been developed and used to provide insight into the transmission dynamics and control of COVID-19. The modeling types used include statistical (20), compartmental/deterministic [e.g., (4, 7–9, 21–23)], stochastic [e.g., (24, 25)], network [e.g., (26)], and agent-based [e.g., (27)]. A notable feature of the model to be developed in the current study is its multigroup nature. Specifically, the total population will be subdivided into two groups, namely those who habitually wear face mask in public and those who do not. Cumulative mortality data for COVID-19 pandemic in the U.S. will be used to parameterize the model. The expected outcome of the study is the determination of the minimum vaccine coverage level needed to effectively curtail (or eliminate) community transmission of COVID-19 in the U.S., and quantify the reduction in the required vaccine coverage if the vaccination program is supplemented with face masks usage (under various face masks efficacy and compliance parameter space). The rest of the paper is organized as follows. The novel multigroup model is formulated in section 2. The parameters of the model are also estimated, based on fitting the model with U.S. COVID-19 mortality data for the third wave of the pandemic. The model is rigorously analyzed, with respect to the asymptotic stability of the disease-free equilibrium of the model, in section 3. A condition for achieving community-wide vaccine-derived herd immunity is also derived. Numerical simulations of the model are reported in section 4. Discussions and concluding remarks are presented in section 5. It is worth stating that this study was carried out between December 2020 and January 2021, when the U.S. was experiencing the third wave of the pandemic (hence, what follows should be viewed in this context).

2. Formulation of Mathematical Model

In order to account for heterogeneity in face masks usage in the community, the total population of individuals in the community at time t, denoted by N(t), is split into the total sub-populations of individuals who do not habitually wear face mask in public (labeled “non-mask users”), denoted by N₁(t), and the total sub-populations of those who habitually wear face mask in public (labeled “mask users”), represented by N₂(t). That is, N(t) = N₁(t) + N₂(t). Furthermore, the sub-population N₁(t) is sub-divided into the mutually-exclusive compartments of unvaccinated susceptible [S_1u(t)], vaccinated susceptible [S_1v(t)], exposed [E₁(t)], pre-symptomatically-infectious [P₁(t)], symptomatically-infectious [I₁(t)], asymptomatically-infectious [A₁(t)], hospitalized [H₁(t)], and recovered [R₁(t)] individuals, so that

$N_1(t)=S_{1u}(t)+S_{1v}(t)+E_1(t)+P_1(t)+I_1(t)+A_1(t)+H_1(t)+R_1(t).$

Similarly, the total sub-population of the mask users, N₂(t), is stratified into the compartments for unvaccinated susceptible [S_2u(t)], vaccinated susceptible [S_2v(t)], exposed [E₂(t)], pre-symptomatically-infectious [P₂(t)], symptomatically-infectious [I₂(t)], asymptomatically-infectious [A₂(t)], hospitalized [H₂(t)], and recovered [R₂(t)] individuals. Hence,

$N_2(t)=S_{2u}(t)+S_{2v}(t)+E_2(t)+P_2(t)+I_2(t)+A_2(t)+H_2(t)+R_2(t).$

2.1. Infection Rates

In this section, the functional form of the infection rate (or effective contact rate) for a susceptible individual in group 1 or 2 will be derived. The model to be formulated has four infectious classes, namely the classes for pre-symptomatic (P_i), symptomatic (I_i), asymptomatic (A_i), and hospitalized (H_i) individuals (i = 1, 2). Hence, the rate at which an individual in group i acquires infection from an infectious individual in any of the four infectious classes is given by the average number of contacts per unit time (measured in days) for susceptible individuals (denoted by c_k; with k = {P_i, I_i, A_i, H_i} and i = 1, 2), times the sum (over all infectious compartments in group i) of the probability of transmission per contact with an infectious individual in group i (denoted by $\widehat{{\beta }_i}$) times the probability that a random infectious contact the susceptible individual makes is with an infectious individual in group i (denoted by ρ_k).

Let c_k be the average number of contacts an individual in epidemiological compartment k makes per unit time. It then follows that the probability that a random contact this individual makes is with someone else in epidemiological compartment k is given by the total number contacts made by everyone in that compartment, denoted by c_ki, divided by the total number of contacts for the entire population. That is, ${\rho }_k=\frac{c_{ki}}{c_{total}}$, where

$\begin{array}{l}{c}_{total}=c_{Siu}{S}_{1u}(t)+c_{Siv}{S}_{1v}(t)+c_{Ei}{E}_1(t)+c_{Pi}{P}_1(t)+c_{Ii}{I}_1(t)\\ +c_{Ai}{A}_1(t)+c_{Hi}{H}_1(t)+c_{Ri}{R}_1(t),i=1,2.\end{array}$

Based on the above definitions, it follows that the infection rate of a susceptible individual in group i, denoted by λ_i (i = 1, 2), is given by

$\begin{array}{l}{\lambda }_1=c_{S1u}\left[\frac{{\widehat{\beta }}_{P_1}{c}_{P1}{P}_1+{\widehat{\beta }}_{I_1}{c}_{I1}{I}_1+{\widehat{\beta }}_{A_1}{c}_{A1}{A}_1+{\widehat{\beta }}_{H_1}{c}_{H1}{H}_1}{c_{total}}\right]\\ +c_{S1v}(1-{\epsilon }_o)\left[\frac{{\widehat{\beta }}_{P_2}{c}_{P2}{P}_2+{\widehat{\beta }}_{I_2}{c}_{P2}{I}_2+{\widehat{\beta }}_{A_2}{c}_{A2}{A}_2+{\widehat{\beta }}_{H_2}{c}_{H2}{H}_2}{c_{total}}\right],& (2.1)\end{array}$

Similarly, the infection rate for a susceptible individual in group 2, denoted by λ₂, is given by:

$\begin{array}{l}{\lambda }_2=(1-{\epsilon }_i)c_{S2u}\left[\frac{{\widehat{\beta }}_{P_1}{c}_{P1}{P}_1+{\widehat{\beta }}_{I_1}{c}_{I1}{I}_1+{\widehat{\beta }}_{A_1}{c}_{A1}{A}_1+{\widehat{\beta }}_{H_1}{c}_{H1}{H}_1}{c_{total}}\right]\\ +c_{S2v}(1-{\epsilon }_i)(1-{\epsilon }_o)\\ \left[\frac{{\widehat{\beta }}_{P_2}{c}_{P2}{P}_2+{\widehat{\beta }}_{I_2}{c}_{P2}{I}_2+{\widehat{\beta }}_{A_2}{c}_{A2}{A}_2+{\widehat{\beta }}_{H_2}{c}_{H2}{H}_2}{c_{total}}\right].& (2.2)\end{array}$

In Equations (2.1) and (2.2), the parameters 0 < ε_o < 1 and 0 < ε_i < 1 represent the outward and inward protective efficacy, respectively, of face masks to prevent the transmission of infection to a susceptible individual (ε_o) as well as prevent the acquisition of infection (ε_i) from an infectious individual. For mathematical tractability (needed to reduce the number of parameters of the model to be developed), we assume that every member of the population has the same number of contacts. That is, we assume that c_S1u = c_S1v = ⋯ = c_R2 = k_c. Hence, c_total = k_cN(t). Let ${\beta }_k={\widehat{\beta }}_k{k}_c$. Using this definition of β_k and c_total = k_cN in Equations (2.1) and (2.2) gives, respectively,

$\begin{array}{l}{\lambda }_1=[\frac{{\beta }_{P_1}{P}_1+{\beta }_{I_1}{I}_1+{\beta }_{A_1}{A}_1+{\beta }_{H_1}{H}_1}{N}\\ +(1-{\epsilon }_o)\frac{{\beta }_{P_2}{P}_2+{\beta }_{I_2}{I}_2+{\beta }_{A_2}{A}_2+{\beta }_{H_2}{H}_2}{N}],& (2.3)\end{array}$

and,

$\begin{array}{l}{\lambda }_2=(1-{\epsilon }_i)[\frac{{\beta }_{P_1}{P}_1+{\beta }_{I_1}{I}_1+{\beta }_{A_1}{A}_1+{\beta }_{H_1}{H}_1}{N}\\ +(1-{\epsilon }_o)\frac{{\beta }_{P_2}{P}_2+{\beta }_{I_2}{I}_2+{\beta }_{A_2}{A}_2+{\beta }_{H_2}{H}_2}{N}].& (2.4)\end{array}$

2.2. Equations of Mathematical Model

Before giving the equations for the two-group vaccination model, it is important to recall that vaccination against COVID-19 in the U.S. is administered to individuals of a certain eligible age [e.g., as of April 2021, people 12 years of age and older for the Pfizer vaccine (28) and 18 years of age and older for the Moderna vaccine]. Consequently, in formulating a model that incorporates COVID-19 vaccines, it is important that demographic parameters (birth and natural death) are included to account for the new cohort of susceptible individuals that reach the minimum eligible age for receiving the vaccine. The equations for the rate of change of the sub-populations of non-mask users (i.e., individuals in group 1) is given by the following deterministic system of nonlinear differential equations (where a dot represents differentiation with respect to time t):

$\begin{array}{lll}{\dot{S}}_{1u}& =& \Pi +{\alpha }_{21}{S}_{2u}-{\lambda }_1{S}_{1u}-({\alpha }_{12}+{\xi }_v+\mu )S_{1u},\\ {\dot{S}}_{1v}& =& {\xi }_v{S}_{1u}+{\alpha }_{21}{S}_{2v}-(1-{\epsilon }_v){\lambda }_1{S}_{1v}-({\alpha }_{12}+\mu )S_{1v},\\ {\dot{E}}_1& =& {\lambda }_1{S}_{1u}+(1-{\epsilon }_v){\lambda }_1{S}_{1v}+{\alpha }_{21}{E}_2-({\alpha }_{12}+{\sigma }_1+\mu )E_1,\\ {\dot{P}}_1& =& {\sigma }_1{E}_1+{\alpha }_{21}{P}_2-({\alpha }_{12}+{\sigma }_P+\mu )P_1,\\ {\dot{I}}_1& =& r{\sigma }_P{P}_1+{\alpha }_{21}{I}_2-({\alpha }_{12}+{\varphi }_{1I}+{\gamma }_{1I}+\mu +{\delta }_{1I})I_1,\\ {\dot{A}}_1& =& (1-r){\sigma }_P{P}_1+{\alpha }_{21}{A}_2-({\alpha }_{12}+{\gamma }_{1A}+\mu )A_1,\\ {\dot{H}}_1& =& {\varphi }_{1I}{I}_1+{\alpha }_{21}{H}_2-({\alpha }_{12}+{\gamma }_{1H}+\mu +{\delta }_{1H})H_1,\\ {\dot{R}}_1& =& {\gamma }_{1I}{I}_1+{\gamma }_{1A}{A}_1+{\gamma }_{1H}{H}_1+{\alpha }_{21}{R}_2-({\alpha }_{12}+\mu )R_1.& (2.5)\end{array}$

where, λ₁ is as defined in (2.3).

In Equation (2.5), the parameter Π is the recruitment rate into the population (this parameter also captures the inflow of new susceptible individuals that have reached the minimum eligibility age for getting a vaccine). Furthermore, α₂₁ is the rate at which individuals in the habitual mask-wearing group 2 change their behavior and move to the non-masking group 1, and α₁₂ is the rate at which individuals in group 1 change their non-masking behavior and move to group 2. For mathematical tractability, we do not distinguish the change of behavior parameters (α₁₂ and α₂₁) for unvaccinated and vaccinated individuals, and we assume that all recruited individuals (at the rate Π) are initially in the non-masking group. The parameter ξ_v represents the per capita vaccination rate, and the vaccine is assumed to induce protective efficacy 0 < ε_v < 1 in all vaccinated individuals (i.e., the vaccine is imperfect). Natural deaths occurs in all epidemiological classes at a rate μ. Individuals in the E₁ class progress to the pre-symptomatic stage at a rate σ₁, and those in the pre-symptomatic class (P₁) transition out of this class at a rate σ_P [a proportion, q, of which become symptomatic, and move to the I class at a rate qσ_P, and the remaining proportion, 1 − q, move to the asymptomatically-infectious class at a rate (1 − q)σ_P]. Symptomatic infectious individuals are hospitalized at a rate ϕ_1I. They recover at a rate γ_1I and die due to the disease at a rate δ_1I. Hospitalized individuals die of the disease at the rate δ_1H.

Similarly, the equations for the rate of change of the sub-populations of mask users (i.e., individuals in group 2) is given by the following system of nonlinear differential equations:

$\begin{array}{lll}{\dot{S}}_{2u}& =& {\alpha }_{12}{S}_{1u}-{\lambda }_2{S}_{2u}-({\alpha }_{21}+{\xi }_v+\mu )S_{2u},\\ {\dot{S}}_{2v}& =& {\xi }_v{S}_{2u}+{\alpha }_{12}{S}_{1v}-(1-{\epsilon }_v){\lambda }_2{S}_{2v}-({\alpha }_{21}+\mu )S_{2v},\\ {\dot{E}}_2& =& {\lambda }_2{S}_{2u}+(1-{\epsilon }_v){\lambda }_2{S}_{2v}+{\alpha }_{12}{E}_1-({\alpha }_{21}+{\sigma }_2+\mu )E_2,\\ {\dot{P}}_2& =& {\sigma }_2{E}_2+{\alpha }_{12}{P}_1-({\alpha }_{21}+{\sigma }_P+\mu )P_2,\\ {\dot{I}}_2& =& q{\sigma }_P{P}_2+{\alpha }_{12}{I}_1-({\alpha }_{21}+{\varphi }_{2I}+{\gamma }_{2I}+\mu +{\delta }_{2I})I_2,\\ {\dot{A}}_2& =& (1-q){\sigma }_P{P}_2+{\alpha }_{12}{A}_1-({\alpha }_{21}+{\gamma }_{2A}+\mu )A_2,\\ {\dot{H}}_2& =& {\varphi }_{2I}{I}_2+{\alpha }_{12}{H}_1-({\alpha }_{21}+{\gamma }_{2H}+\mu +{\delta }_{2H})H_2,\\ {\dot{R}}_2& =& {\gamma }_{2I}{I}_2+{\gamma }_{2A}{A}_2+{\gamma }_{2H}{H}_2+{\alpha }_{12}{R}_1-({\alpha }_{21}+\mu )R_2,& (2.6)\end{array}$

with λ₂ defined in (2.4). Thus, Equations (2.5) and (2.6) represent the multi-group model for assessing the population-level impact of face masks usage and vaccination on the transmission dynamics and control of COVID-19 in a community. The flow diagram of the model {(2.5), (2.6)} is depicted in Figure 1 (the state variables and parameters of the model are described in Tables 1, 2, respectively).

FIGURE 1

Figure 1. Flow diagram of the model {(2.5), (2.6)}.

TABLE 1

Table 1. Description of the state variables of the model {(2.5), (2.6)}.

TABLE 2

Table 2. Description of the parameters of the model {(2.5), (2.6)}.

The multi-group model {(2.5), (2.6)} is an extension of the two-group mask-use model in (7) by, inter alia:

(i) allowing for back-and-forth transitions between the two groups (mask-users and non-mask-users), to account for human behavioral changes vis a vis decision to either be (or not to be) a face mask user in public;

(ii) incorporating an imperfect vaccine, which offers protective efficacy (0 < ε_v <1) against acquisition of COVID-19 infection;

(iii) allowing for disease transmission by pre-symptomatic and asymptomatically-infectious individuals.

2.3. Data Fitting and Parameter Estimation

In this section, cumulative COVID-19 mortality data for the U.S. (for the period October 12, 2020–January 20, 2021) will be used to fit the model (2.5)–(2.6) in the absence of vaccination. The fitting will allow us to estimate some of the key (unknown) parameters of the model. In particular, the parameters to be estimated from the data are the community transmission rate for individuals who do not wear face masks in public (β₁), the transmission rate for individuals who habitually wear face masks in public (β₂), the inward efficacy of masks in preventing disease acquisition by susceptible individuals who habitually wear face masks (ε_i), the outward efficacy of masks to prevent the spread of disease by infected individuals who habitually wear face masks (ε_o), the rate at which people who do not wear masks adopt a mask-wearing habit (α₁₂), the rate at which those who habitually wear face masks stop wearing masks in public (α₂₁), and the mortality rates of symptomatic infectious and hospitalized individuals (δ_i and δ_h, respectively). It should be mentioned that modification parameters η_P, η_I, η_A, , and η_H relating to disease transmission by pre-symptomatic infectious, symptomatic infectious, asymptomatic infectious, and hospitalized individuals, respectively, are introduced in the forces of infection λ₁ and λ₂, so that β_j = η_jβ_k (j ∈ {P_k, I_k, A_k, H_k}, k ∈ {1, 2}). The model was fitted using a standard nonlinear least squares approach, which involved using the inbuilt MATLAB minimization function “lsqcurvefit” to minimize the sum of the squared differences between each observed cumulative mortality data point and the corresponding cumulative mortality point obtained from the model (2.5)–(2.6) in the absence of vaccination (4, 29, 30). The choice of mortality over case data is motivated by the fact that mortality data for COVID-19 is more reliable than case data [see (8) for details]. The estimated values of the fitted parameters, together with their 95% confidence intervals, are tabulated in Table 3. The (fixed) values of the remaining parameters of the model are tabulated in Table 4. Figure 2 depicts the fitting of the model to the observed cumulative COVID-19 mortality data for the U.S. Furthermore, Figure 2 compares the simulations of the model using the fitted (estimated) and fixed parameters (given in Tables 3, 4) with the observed daily COVID-19 mortality for the United States.

TABLE 3

Table 3. Estimated (fitted) parameter values and their 95% confidence intervals for the model (2.5)–(2.6) in the absence of vaccination, using COVID-19 mortality data for the U.S. for the period from October 12, 2020 to January 20, 2021.

TABLE 4

Table 4. Baseline values of the fixed parameters of the model (2.5)–(2.6).

FIGURE 2

Figure 2. (A) Observed cumulative mortality (red dots), and the predicted cumulative mortality (blue curve) for the U.S. generated using the model (2.5)–(2.6) (in the absence of vaccination) for the period from October 12, 2020 to January 20, 2021. (B) Simulations of the model (2.5)–(2.6) (without vaccination) using the estimated (fitted) and fixed parameters tabulated in Tables 3, 4, respectively.

3. Mathematical Analysis

Since the model {(2.5), (2.6)} monitors the temporal dynamics of human populations, all of its state variables and parmeters are non-negative. Consider the following biologically-feasible region for the model:

$\begin{array}{l}\Omega =\{(S_{1u},S_{1v},S_{2u},S_{2v},E_1,E_2,P_1,P_2,I_1,I_2,A_1,A_2,H_1,H_2,R_1,R_2)\\ \in R_{+}^{16}:N(t)\le \frac{\Pi }{\mu }\}.& (3.1)\end{array}$

Theorem 3.1. The region Ω is positively-invariant with respect to the model {(2.5), (2.6)}.

Proof: Adding all the equations of the model {(2.5), (2.6)} gives

$\begin{array}{ll}\dot{N}=\Pi -\mu N-{\delta }_{1I}{I}_1-{\delta }_{1H}{H}_1-{\delta }_{2I}{I}_2-{\delta }_{2H}{H}_2.& (3.2)\end{array}$

Recall that all parameters of the model {(2.5), (2.6)} are non-negative. Thus, it follows, from (3.2), that

$\begin{array}{ll}\dot{N}\le \Pi -\mu N.& (3.3)\end{array}$

Hence, if $N>\frac{\Pi }{\mu }$, then Ṅ < 0. Furthermore, by applying a standard comparison theorem (37) on (3.3), we have:

$N(t)\le N(0)e^{-\mu t}+\frac{\Pi }{\mu }(1-e^{-\mu t}).$

In particular, $N(t)\le \frac{\Pi }{\mu }$ if $N(0)\le \frac{\Pi }{\mu }.$ If $N(0)>\frac{\Pi }{\mu }$ [i.e., N(0) is outside Ω], then $N(t)>\frac{\Pi }{\mu }$, for all t > 0 but with $\underset{t\to \infty }{\lim }N(t)=\frac{\Pi }{\mu }$ (and this type of solution trajectory strives to enter the region Ω). Thus, every solution of the model {(2.5), (2.6)} with initial conditions in Ω remains in Ω for all time t > 0. In other words, the region Ω is positively-invariant and attracts all initial solutions of the model {(2.5), (2.6)}. Hence, it is sufficient to consider the dynamics of the flow generated by {(2.5), (2.6)} in Ω (where the model is epidemiologically- and mathematically well-posed) (38). ⎕

3.1. Asymptotic Stability of Disease-Free Equilibrium

The model {(2.5), (2.6)} has a unique disease-free equilibrium (DFE), obtained by setting all the infected compartments of the model to zero, given by:

$\begin{array}{l}{\text{E}}_0:(S_{1u}^{*},S_{1v}^{*},S_{2u}^{*},S_{2v}^{*},E_1^{*},E_2^{*},P_1^{*},P_2^{*},I_1^{*},I_2^{*},A_1^{*},A_2^{*},H_1^{*},H_2^{*},R_1^{*},R_2^{*})\\ =(S_{1u}^{*},S_{1v}^{*},S_{2u}^{*},S_{2v}^{*},0,0,0,0,0,0,0,0,0,0,0,0),\end{array}$

where,

$\begin{array}{ll}{S}_{1u}^{*}& =\frac{\Pi ({\alpha }_{21}+{\xi }_v+\mu )}{({\xi }_v+\mu )({\xi }_v+{\alpha }_{12}+{\alpha }_{21}+\mu )},\\ S_{1v}^{*}& =\frac{\Pi ({\mu }^2{\xi }_v+2\mu \Pi {\alpha }_{21}{\xi }_v+{\alpha }_{21}^2{\xi }_v+\mu {\xi }_v^2+{\alpha }_{21}^2{\xi }_v^2)}{\mu (\mu +{\alpha }_{12}+{\alpha }_{21})(\mu +{\xi }_v)(\mu +{\xi }_v+{\alpha }_{12}+{\alpha }_{21})},\\ S_{2u}^{*}& =\frac{\Pi {\alpha }_{12}}{({\xi }_v+\mu )({\xi }_v+{\alpha }_{12}+{\alpha }_{21}+\mu )},\\ S_{2v}^{*}& =\frac{\Pi {\xi }_v{\alpha }_{12}(2\mu +{\alpha }_{12}+{\alpha }_{21}+{\xi }_v)}{\mu (\mu +{\alpha }_{12}+{\alpha }_{21})(\mu +{\xi }_v)(\mu +{\xi }_v+{\alpha }_{12}+{\alpha }_{21})}.\end{array}$

The local asymptotic stability property of the DFE (E₀) can be explored using the next generation operator method (39, 40). In particular, using the notation in (39), it follows that the associated non-negative matrix (F) of new infection terms, and the M-matrix (V), of the linear transition terms in the infected compartments, are given, respectively, by (where the entries f_i and g_i, i = 1, ⋯ , 8, of the non-negative matrix F, are given in Appendix I):

$\begin{array}{l}F=\left[\begin{array}{cccccccccc}0& f_1& f_2& f_3& f_4& 0& f_5& f_6& f_7& f_8\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& g_1& g_2& g_3& g_4& 0& g_5& g_6& g_7& g_8\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right],\end{array}$

and,

$\begin{array}{lll}V& =& \left[\begin{array}{cccccccccc}{K}_1& 0& 0& 0& 0& -{\alpha }_{21}& 0& 0& 0& 0\\ -{\sigma }_1& K_2& 0& 0& 0& 0& -{\alpha }_{21}& 0& 0& 0\\ 0& -r{\sigma }_p& K_3& 0& 0& 0& 0& -{\alpha }_{21}& 0& 0\\ 0& -(1-r){\sigma }_p& 0& K_4& 0& 0& 0& 0& -{\alpha }_{21}& 0\\ 0& 0& -{\varphi }_{1I}& 0& K_5& 0& 0& 0& 0& -{\alpha }_{21}\\ -{\alpha }_{12}& 0& 0& 0& 0& K_6& 0& 0& 0& 0\\ 0& -{\alpha }_{12}& 0& 0& 0& 0& K_7& 0& 0& 0\\ 0& 0& -{\alpha }_{12}& 0& 0& 0& -q{\sigma }_p& K_8& 0& 0\\ 0& 0& 0& -{\alpha }_{12}& 0& 0& -(1-q){\sigma }_p& 0& K_9& 0\\ 0& 0& 0& 0& -{\alpha }_{12}& 0& 0& -{\varphi }_{2I}& 0& K_{10}\end{array}\right],\end{array}$

where K₁ = α₁₂ + σ₁ + μ, K₂ = α₁₂ + σ_P + μ, K₃ = α₁₂ + ϕ_1I + γ_1I + μ + δ_1I, K₄ = α₁₂ + γ_1A + μ, K₅ = α₁₂ + γ_1H + μ + δ_1H, K₆ = α₂₁ + σ₂ + μ, K₇ = α₂₁ + σ_P + μ, K₈ = α₂₁ + ϕ_2I + γ_2I + μ + δ_2I, K₉ = α₂₁ + γ_2A + μ and K₁₀ = α₂₁ + γ_2H + μ + δ_2H. The theoretical analysis will be carried out for the special case of the model {(2.5), (2.6)} in the absence of the back-and-forth transitions between the no-mask and mask-user groups (i.e., the special case of the model with α₁₂ = α₂₁ = 0). This is needed for mathematical tractability. It follows that the control reproduction number of the model {(2.5), (2.6)} (with α₁₂ = α₂₁ = 0), denoted by ${{R}}_c$, is given by (where ρ is the spectral radius):

$\begin{array}{l}{ℛ}_c=\rho (F{V}^{-1})\\ =\frac{\begin{array}{l}{\sigma }_1[S_{1u}^{*}+(1-{\epsilon }_v)S_{1v}^{*}](({\overline{K}}_5[r{\overline{K}}_4{\beta }_{I_1}+(1-r){\overline{K}}_3{\beta }_{A_1}]+r{\overline{K}}_4{\varphi }_{1I}{\beta }_{H_1}){\sigma }_p\hfill \\ +{\overline{K}}_3{\overline{K}}_4{\overline{K}}_5{\beta }_{P_1})\hfill \end{array}}{N^{*}{\displaystyle \prod _{i=1}^5{\overline{K}}_i}},& (3.4)\end{array}$

where, ${\bar{K}}_1={\sigma }_1+\mu ,{\bar{K}}_2={\sigma }_P+\mu ,{\bar{K}}_3={\varphi }_{1I}+{\gamma }_{1I}+\mu +{\delta }_{1I},{\bar{K}}_4={\gamma }_{1A}+\mu ,{\bar{K}}_5={\gamma }_{1H}+\mu +{\delta }_{1H}$. The result below follows from Theorem 2 of (39).

Theorem 3.2. The DFE (E₀) of the model {(2.5), (2.6)}, with α₁₂ = α₂₁ = 0, is locally-asymptotically stable if ${{R}}_c<1$, and unstable if ${{R}}_c>1$.

The threshold quantity ${{R}}_c$ is the control reproduction number of the model {(2.5), (2.6)}. It measures the average number of new COVID-19 cases generated by a typical infectious individual introduced into a population where a certain fraction of the population is protected (via the use of interventions, such as face mask, social-distancing, and/or vaccination). The epidemiological implication of Theorem 3.2 is that a small influx of COVID-19 cases will not generate an outbreak in the community if the control reproduction number $({{R}}_c)$ is brought to, and maintained at a, value less than unity. In the absence of public health interventions (i.e., in the absence of vaccination, face mask usage and social-distancing), the control reproduction number (${{R}}_c$) reduces to the basic reproduction number (denoted by ${{R}}_0$), given by.

$\begin{array}{l}{ℛ}_0={ℛ}_c{|}_{{\epsilon }_0={\epsilon }_i={\epsilon }_v=S_{1v}^{*}=S_{2v}^{*}=0}\\ =\frac{{\sigma }_1(({\overline{K}}_5[r{\overline{K}}_4{\beta }_{I_1}+(1-r){\overline{K}}_3{\beta }_{A_1}]+r{\overline{K}}_4{\varphi }_{1I}{\beta }_{H_1}){\sigma }_p+{\overline{K}}_3{\overline{K}}_4{\overline{K}}_5{\beta }_{P_1})}{{\displaystyle \prod _{i=1}^5{\overline{K}}_i}}& (3.5)\end{array}$

3.2. Derivation of Vaccine-Induced Herd Immunity Threshold

Herd immunity is a measure of the minimum percentage of the number of individuals in a community that is susceptible to a disease that need to be protected (i.e., become immune) so that the disease can be eliminated from the population. There are two main ways to achieve herd immunity, namely through acquisition of natural immunity (following natural recovery from infection with the disease) or by vaccination. Vaccination is the safest and fastest way to achieve herd immunity (41, 42). For vaccine-preventable diseases, such as COVID-19, not every susceptible member of the community can be vaccinated, for numerous reasons (such as individuals with certain underlying medical conditions, infants, pregnant women, or those who opt out of being vaccinated for various reasons etc.) (9). So, the question, in the context of vaccine-preventable diseases, is what is the minimum proportion of susceptible individuals that we need to vaccinate in order to achieve herd immunity (so that those individuals that cannot be vaccinated will become protected owing to the community-wide herd-immunity). In this section, a condition for achieving vaccine-derived herd immunity in the U.S. will be derived.

Let $f_v=S_{1v}^{*}/N^{*}$, with N^* = Π/μ, be the proportion of susceptible individuals in Group 1 that have been vaccinated at the disease-free equilibrium (E₀). Using this definition in Equation (3.4) gives:

$\begin{array}{ll}{ℛ}_c=\frac{\begin{array}{l}{\sigma }_1(1-{\epsilon }_v{f}_v)(({\overline{K}}_5[r{\overline{K}}_4{\beta }_{I_1}+(1-r){\overline{K}}_3{\beta }_{A_1}]+r{\overline{K}}_4{\varphi }_{1I}{\beta }_{H_1}){\sigma }_p\hfill \\ +{\overline{K}}_3{\overline{K}}_4{\overline{K}}_5{\beta }_{P_1})\hfill \end{array}}{{\displaystyle \prod _{i=1}^5{\overline{K}}_i}}.& (3.6)\end{array}$

Setting ${{R}}_c$, in Equation (3.6), to unity and solving for f_v gives the herd immunity threshold (denoted by $f_v^c$) in terms of the basic reproduction number (9, 21):

$\begin{array}{ll}\begin{array}{c}{f}_v^c=\frac{1}{{\epsilon }_v}\left(1-\frac{1}{{}_0}\right)(\text{for }{R}_0>1).\end{array}& (3.7)\end{array}$

It follows from (3.6) and (3.7) that ${{R}}_c<(>)1$ if $f_v>(<)f_v^c$. Further, ${{R}}_c=1$ whenever $f_v=f_v^c$. This result is summarized below:

Theorem 3.3. Consider the special case of the model {(2.5), (2.6)} with α₁₂ = α₂₁ = 0. Vaccine-induced herd immunity can be achieved in the U.S., using an imperfect anti-COVID vaccine, if $f_v>f_v^c$ (i.e., if ${{R}}_c<1$). If $f_v<f_v^c$ (i.e., if ${ℛ}_c>1$, then the vaccination program will fail to eliminate the COVID-19 pandemic in the U.S.

The epidemiological implication of Theorem 3.3 is that the use of an imperfect anti-COVID vaccine can lead to the elimination of the COVID-19 pandemic in the U.S. if the sufficient number of individuals residing in the U.S. is vaccinated, such that $f_v>f_v^c$. The Vaccination program will fail to eliminate the pandemic if the vaccine coverage level is below the aforementioned herd immunity threshold (i.e., if $f_v<f_v^c$). Although vaccination, no matter the coverage level, is always useful (i.e., vaccination will always reduce the associated reproduction number, ${{R}}_c$, thereby reducing disease burden, even if the program is unable to bring the reproduction number to a value less than unity), elimination can only be achieved if the herd immunity threshold is reached (i.e., disease elimination is only feasible if the associated reproduction number of the model is reduced to, and maintained at, a value less than unity). The pandemic will persist in the U.S. if ${{R}}_c>1$.

Figure 3A depicts the cumulative mortality of COVID-19 in the U.S. for various steady-state vaccination coverage levels (f_v). This figure shows a decrease in cumulative mortality with increasing vaccination coverage. In particular, a marked decrease in cumulative mortality, in comparison to the baseline cumulative mortality (blue curve in Figure 3A), is recorded when herd immunity (i.e., when $f_v>f_v^c)$ is attained (green curve of Figure 3A). While a noticeable decrease in the cumulative mortality is also observed when the vaccine coverage equals the herd immunity threshold (gold curve of Figure 3A), the cumulative mortality dramatically increases (in comparison to the baseline, depicted by the blue curve of this figure) if the vaccine coverage is below the herd immunity threshold (magenta curve of Figure 3A).

FIGURE 3

Figure 3. Assessment of the effects of vaccine coverage (f_v) and efficacy (ε_v) on COVID-19 dynamics in the U.S. (A) Simulations of the model {(2.5), (2.6)}, with α₁₂ = α₂₁ = 0, showing the cumulative COVID-19 mortality in the U.S., as a function of time, for various values of vaccine coverage. Parameter values used are as given by the baseline values in Tables 3, 4, with α₁₂ = α₂₁ = 0 and various values of f_v. Magenta curve ($f_v=0.3021<0.5900=f_v^c$), blue curve (baseline parameter values, and baseline level of social-distancing compliance inherent in the cumulative mortality data, for the period October 12, 2020–January 20, 2021, used to fit the model), gold curve ($f_v=0.5900=f_v^c$) and green curve ($f_v=0.9216>0.5900=f_v^c$). The observed cumulative deaths data, fitted to the baseline scenario predicted by the model (blue curve), is shown in red dots. (B) Contour plot of the control reproduction number (${{R}}_c$) of the model {(2.5), (2.6)}, with α₁₂ = α₂₁ = 0, as a function of vaccine coverage (f_v) and vaccine efficacy (ε_v). Parameter values used are as given in Tables 3, 4, with α₁₂ = α₂₁ = 0.

The effect of vaccination coverage (f_v) and efficacy (ε_v) on the control reproduction number (${{R}}_c$) is assessed by generating a contour plot of ${{R}}_c$, as a function of f_v and ε_v. The results obtained (Figure 3B) show that the values of the control reproduction number for the U.S., during the simulation period (October 12, 2020–January 20, 2021), range from 0.4 to 2.2. Further, this figure shows that the control reproduction number decreases with increasing values of vaccination efficacy and coverage. For example, using the AstraZeneca vaccine (with efficacy ε_v = 0.7), about 80% of the U.S. population needs to be successfully vaccinated (with the two AstraZeneca doses) in order to bring the control reproduction number to a value less than unity. In other words, this figure shows that herd immunity can be achieved using the AstraZeneca vaccine in the U.S. if at least 80% of the populace received the two doses of the AstraZeneca vaccine. Using either the Pfizer or Moderna vaccine (each with efficacy of about 95%), on the other hand, the control reproduction number can be brought to a value less than unity (i.e., achieve herd immunity) if at least 60% of the U.S. populace received the two doses of either vaccine. Thus, this figure shows that the prospect of achieving vaccine-derived herd immunity using any of the three vaccines considered in this study (AstraZeneca, Pfizer, and Moderna) is promising if the coverage is moderately-high enough (with the prospect far more likely to be achieved using the Pfizer or Moderna vaccine, in comparison to using the AstraZeneca vaccine).

We also explored the potential impact of additional social-distancing on the minimum vaccination coverage needed to achieve herd immunity. It should, first of all, be stressed that, since our model was parameterized using the cumulative mortality data during the third wave of the pandemic in the U.S. (October 12, 2020–January 21, 2021), the effects of other non-pharmaceutical interventions, such as face masks usage and social-distancing, are already embedded into the results/data. In other words, the data (or the parametrization of our model) already includes some baseline level of these interventions. Specifically, we assume that the cumulative mortality data includes a baseline level of social-distancing compliance in the population (which is, clearly, quite high compared to what it was during the early stages of the pandemic in the U.S.) We now ask the question as to whether or not the minimum requirement for 80 and 60% coverage needed to achieve herd immunity, using the AstraZeneca or Pfizer/Moderna vaccine, respectively, can be reduced if the baseline social-distancing compliance is increased. In this study, we model social-distancing compliance by multiplying the effective contact rates (β₁ and β₂) with the factor 1−c_s, where 0 < c_s ≤ 1 is a measure of the additional social-distancing compliance (to the baseline social-distancing compliance achieved during the beginning of our simulation period; that is, by October 12, 2020).

We simulated the model {(2.5), (2.6)} using various values of c_s, and the results obtained are tabulated in Table 5. This table shows that if an additional 5% of the U.S. population observe social-distancing in public (in addition to the baseline social-distancing compliance achieved by October 12, 2020), the minimum vaccine coverages required to achieve herd immunity using the AstraZeneca and Pfizer/Moderna vaccines reduce, respectively, to 77 and 56.4%. Furthermore, if the increase in baseline social-distancing compliance is 10%, the minimum coverage needed to achieve herd immunity further reduce (but marginally) to 73 and 54%, respectively. However, when the increase in baseline social-distancing compliance is 30%, herd immunity can be achieved using the AstraZeneca vaccine by vaccinating only 53% of the U.S. population with this vaccine. For this scenario, only about 39% of the U.S. population needs to be vaccinated to achieve herd immunity if either the Pfizer or Moderna vaccine is used.

TABLE 5

Table 5. Vaccine-induced herd immunity threshold ($f_v^c$) for the U.S. for various levels of increases in baseline social-distancing compliance (c_s).

Thus, this study shows that the prospect of achieving vaccine-derived herd immunity in the U.S. using any of the three vaccines considered in this study is greatly enhanced if the vaccination program is complemented with an increased (and sustained) social-distancing strategy (from its baseline effectiveness and coverage level). In other words, if more people living in the U.S. will continue to observe social-distancing (e.g., additional 30% from the baseline social-distancing compliance), then COVID-19 elimination can be achieved if roughly only half the population is vaccinated using the AstraZeneca vaccine, or 2 in 5 vaccinated if either the Pfizer or Moderna vaccine is used instead. The U.S. is currently using the latter vaccines. Hence, with about 30% additional social-distancing compliance, we would only need to vaccinate about 2 in 5 residents of the U.S. to achieve vaccine-derived herd immunity (hence, eliminate the pandemic).

4. Numerical Simulations: Assessment of Control Strategies

The model {(2.5), (2.6)} will now be simulated to assess the population-level impact of the various intervention strategies described in this study. In particular, our objective is to assess the impact of social-distancing and face mask usage, implemented as sole interventions and in combination with any of the three vaccines (namely the AstraZeneca, Moderna, and Pfizer vaccines), on curtailing (or eliminating) the burden of the COVID-19 pandemic in the United States. Unless otherwise stated, the simulations will be carried out using the estimated (fitted) and fixed baseline values of the parameters of the model tabulated in Tables 3, 4. Furthermore, unless otherwise stated, the baseline initial size of the population of individuals who habitually wear face masks in public (assumed to be 30%), denoted by N₂(0), will be used in the simulations. The numerical simulation results for the baseline scenario (i.e., where baseline values of the parameters of the model, as well as the baseline initial size of the mask-wearing population, are used) will be illustrated in blue curves in the forthcoming figures. Furthermore, all numerical simulations will be carried out for the period starting from October 12, 2020 (which corresponds to the onset of the third wave of the pandemic in the United States).

4.1. Assessing the Impact of Initial Population of Face Mask Wearers

The model (2.5)–(2.6) is simulated to assess the community-wide impact of using face masks, as the sole intervention, in curtailing the spread of the pandemic in the United States. Specifically, we simulate the model using the baseline values of the parameters in Tables 3, 4 and various values of the initial size of the population of individuals who habitually wear face masks in public since the beginning of the pandemic in the United States [denoted by N₂(0)]. It should be noted that the parameters associated with other interventions (e.g., vaccination-related and social-distancing-related parameters) are kept at their baseline values given in Tables 3, 4. Although a sizable number of U.S. residents (notably individuals categorized in the first-tier priority group for receiving the COVID-19 vaccine, such as frontline healthcare workers, individuals at residential care facilities, the elderly etc.) had already been vaccinated using one of the two vaccines that received FDA Emergency Use Authorization in December 2020 [20.54 million vaccines doses of Pfizer and Moderna vaccines had already been administered in the U.S. as of January 23, 2021 (43)], these vaccines were not expected to be widely available to the general public until some time in March or April, 2021. Consequently, we set March 15, 2021 as our reference point for when we expect the vaccines to be widely available to the general public. Under this scenario (of vaccines expected to be widely available a few months after the initial starting point of our simulations, namely October 12, 2020), the objective of this set of simulations is to assess the impact of face masks usage, as a sole intervention, in controlling the spread of the pandemic in the U.S. before the two vaccines that received FDA EUA (Pfizer or Moderna) become widely available to the general U.S. public (to the extent that high vaccination coverage, such as vaccinating one million U.S. residents per day, can be realistically achieved). The new U.S. administration aims to vaccinate 100 million residents during its first 100 days.

The simulation results obtained, depicted in Figure 4, show (generally) that the early adoption of face masks control measures [as measured in terms of the initial proportion of the populace who choose to habitually wear face masks whenever they are out in the public, denoted by N₂(0)] play a vital role in curtailing the COVID-19 mortality in the U.S., particularly for the case when mask-wearers do not abandon their masks-wearing habit (i.e., α₂₁ = 0). For the case where the parameters associated with the back-and-forth transitions between the masking and non-masking sub-populations (i.e., α₁₂ and α₂₁) are maintained at their baseline values (given in Tables 3, 4), this figure shows that the size of the initial proportion of individuals who wear face masks has a significant impact on the cumulative COVID-19 mortality, as measured in relation to the cumulative mortality recorded when the initial proportion of mask wearers is at baseline level (blue curves in Figure 4). In particular, a 34% reduction in the cumulative mortality, in comparison to the cumulative mortality for the baseline scenario, will be recorded by March 15, 2021, if the initial proportion of mask-wearers is 40% (Figure 4A, magenta curve). Furthermore, the reduction in cumulative mortality by March 15, 2021 increases to 52% if the initial proportion of mask-wearers is 75% (Figure 4A, green curve). On the other hand, for the case when mask-wearers remain mask-wearers since the beginning of the simulation period (i.e., since October 12, 2020), so that α₂₁ = 0, while non-mask wearers (i.e., those in Group 1) can change their behavior and become mask-wearers (i.e., α₁₂≠0), our simulations show that the initial proportion of individuals who adopt masking only marginally affects the cumulative mortality (Figure 4B), in relation to the scenario in Figure 4A, where both α₁₂ and α₂₁ are nonzero). In particular, if 40% of the U.S. population adopted mask-wearing right from the aforementioned October 12, 2020, up to 37% of the baseline COVID-19 mortality can be averted (Figure 4B, magenta curve), in comparison to the baseline (Figure 4B, blue curve). Furthermore, the reduction in baseline cumulative mortality rises to 53% if three in every four Americans opted to wear face masks since the beginning of the simulation period (Figure 4B, green curve). This also represents a marginal increase in the cumulative deaths averted, in comparison to the scenario when α₁₂ ≠ 0 and α₂₁ ≠ 0 (Figure 4A, green curve).

FIGURE 4

Figure 4. Assessment of the impact of face mask usage, as a sole intervention, on COVID-19 pandemic in the U.S. Simulations of the model (2.5)–(2.6), showing cumulative mortality, as a function of time, for (A) face mask transition parameters (α₁₂ and α₂₁) maintained at their baseline values, (B) mask-wearers strictly adhere to wearing masks (α₂₁ = 0) and non-mask-wearers transit to mask wearing at their baseline rate (α₁₂ ≠ 0), and (C) non-mask wearers and mask-wearers do not change their behavior (i.e., α₁₂ = α₂₁ = 0). Mask use change is implemented in terms of changes in the initial size of the population of individuals who wear face masks (from the onset of simulations, on October 12, 2020). Blue curves (in each of the plots) represent the baseline scenarios where the initial size of the population of mask wearers is fixed at 30%, and the transition parameters, α₁₂ and α₂₁, are maintained at their baseline values. Parameter values used in the simulations are as given by the baseline values in Tables 3, 4, with different values of α₁₂ and α₂₁ (except for the blue curves, where α₁₂ and α₂₁ are fixed at their baseline values).

For the case when no back-and-forth transitions between the two (mask-wearing and non-mask-wearing) groups is allowed (i.e., when α₁₂ = α₂₁ = 0), our simulations show a far more dramatic effect of face mask usage on COVID-19 mortality (Figure 4C). For instance, this figure shows that higher cumulative mortality is recorded, in comparison to the baseline masks use scenario, when the initial size of the population of mask wearers is 40% (Figure 4C, magenta curve), in comparison to the blue curve of the same figure). Specifically, this represents a 55% increase, in comparison to the baseline cumulative mortality. This simulation result suggests that the 40% initial size of the populace wearing face masks, during the onset of the third wave of the pandemic in the U.S. (starting October 12, 2020), falls below the mask-use compliance threshold level needed to reduce the cumulative mortality during the third wave. On the other hand, if the initial size of the population of face masks wearers is increased to 50%, a decrease (and not an increase) in cumulative mortality is recorded, in comparison to the cumulative mortality for the baseline scenario (Figure 4C, gold curve, in comparison to the blue curve of the same figure). Further dramatic reduction (52%), in relation to the baseline scenario, will be achieved if the initial size of the mask-wearing population is increased to 75% (Figure 4C, green curve, in comparison to the blue curve of the same figure). Thus, these simulations show that, for the case when no change of mask-wearing behavior is allowed (i.e., everyone remains in their original group), there is a threshold value of the initial size of the population of mask wearers above (below) which the cumulative mortality is decreased (increased). Specifically, this simulation shows that (for this scenario with α₁₂ = α₂₁ = 0), at least half the population need to be wearing face masks right from the beginning of the epidemic to ensure greater reduction in cumulative mortality, in comparison to the baseline scenario (when the initial size of the mask-wearing sub-population is 30%).

In summary, comparing the same initial mask coverage (i.e., the same curve colors) in Figures 4A–C, it is clear that the scenario where individuals are allowed to change their behaviors from not wearing face masks to wearing face masks (i.e., α₁₂≠0), but masks wearers do not abandon masks wearing (i.e., α₂₁ = 0), depicted in Figure 4B, resulted in saving more lives (albeit only slightly), compared to the scenarios where no change of behavior is allowed for members of each group (Figure 4C) or members of both groups can change their behavior (Figure 4A). In other words, our study emphasizes the need for non-maskers to adopt a mask-wearing culture (i.e., α₁₂≠0) and for habitual mask-wearers not to abandon their mask-wearing habit (i.e., α₂₁ = 0).

4.2. Assessing the Impact of Additional Social-Distancing Compliance

In this section, we carry out numerical simulations to assess the potential impact of increases in the baseline social-distancing compliance (c_s) on the control of the pandemic. Specifically, the model {(2.5), (2.6)} will be simulated using the baseline parameter values tabulated in Tables 3, 4 with various values of c_s (corresponding to the various levels of the increase in baseline social-distancing compliance in the U.S., starting from October 12, 2020). It should be noted that, for these simulations, the baseline initial size of the masking population, N₂(0), is maintained. Furthermore, vaccine-related parameter values are maintained at their baseline levels in Tables 3, 4.

The simulation results obtained, depicted in Figure 5, show that, in the absence of additional increase in baseline social-distancing (i.e., c_s = 0, so that social-distancing compliance is maintained at the baseline level inherent in the cumulative mortality data by October 12, 2020), the U.S. would record about 500,000 cumulative deaths by March 15, 2021 (Figure 5A, blue curve). For this (baseline social-distancing) scenario, the U.S. would have recorded a peak daily mortality of about 3, 000 deaths on January 5, 2021 (Figure 5B, blue curve). The simulations in Figure 5 further show that the cumulative mortality (Figure 5A) and daily mortality (Figure 5B) decrease with increasing levels of the additional social-distancing compliance (c_s) in the population. For example, if the baseline social-distancing achieved during the onset of the third wave of the pandemic in the U.S. is further increased by only 5%, the simulation results show that up to a 19% of the cumulative mortality can be averted by March 15, 2021 (Figure 5A, magenta curve), in comparison to the baseline social-distancing scenario (Figure 5A, blue curve). Similarly, for this 5% increase in social-distancing (in relation to the baseline), up to 36% reduction in daily mortality can be achieved (Figure 5B, magenta curve), in comparison to the baseline scenario (Figure 5B, blue curve), and the pandemic would have peaked a week earlier (in late December 2020; the daily mortality at this peak would have been 1, 900), in comparison to the peak recorded in the baseline social-distancing scenario (Figure 5B, blue curve). More dramatic reduction in mortality will be recorded if the level of additional social-distancing compliance is further increased. For instance, if the baseline social-distancing compliance is increased by 10%, our simulations show that about 31% of the cumulative deaths recorded for the case with baseline social-distancing scenario (Figure 5A, blue curve) would have been averted (Figure 5A, gold curve). For this scenario, up to 59% of the daily deaths would have been prevented and the pandemic would have peaked in mid December 2020 (the daily mortality at this peak would have been 1, 229), as depicted in the gold curve of Figure 5B. Finally, if the baseline social-distancing compliance is increased by 30%, the pandemic would have failed to generate a major outbreak in the U.S. (Figure 5, green curves). In particular, the cumulative mortality for the U.S. by March 15, 2021 will be about 252, 400 (as against the nearly 400,000 fatalities that were recorded), as shown by the green curve of Figure 5A, in comparison to the blue curve of the same figure.

FIGURE 5

Figure 5. Assessment of the singular impact of increases in baseline social-distancing compliance on COVID-19 pandemic in the U.S. Simulations of the model (2.5)–(2.6) showing (A) cumulative mortality, as a function of time; (B) daily mortality, as a function of time, for various levels of increases in baseline social-distancing (SD) compliance (c_s) attained during the third wave of the pandemic in the United States. Parameter values used in the simulations are as given by the baseline values in Tables 3, 4, with β₁ and β₂ multiplied by (1 − c_s).

In summary, the results in Figure 5 show that COVID-19 could have been effectively suppressed in the U.S. if the baseline social-distancing compliance (recorded during the onset of the third wave of the pandemic in early October 2020) is increased by about 10–30%. These (recommended) increases in social-distancing compliance seem reasonably attainable. Hence, our study suggests that a moderate increase in the baseline social-distancing compliance will lead to the effective control of the COVID-19 pandemic in the U.S. This (increase in baseline social-distancing, as well as face masks usage) should be sustained until herd immunity is attained.

4.3. Assessment of Combined Impact of Vaccination and Social-Distancing

The model (2.5)–(2.6) will now be simulated to assess the community-wide impact of the combined vaccination and social-distancing strategy. Although the two vaccines received FDA's EUA by mid December 2020, we assume a hypothetical situation in which the vaccination started by mid October 2020 (the reason is to ensure consistency with the cumulative mortality data we used, which started from October 12, 2020 corresponding to the onset of the third wave of the pandemic in the United States). We consider the AstraZeneca vaccine (with estimated efficacy of 70%) and the two FDA-EUA vaccines (Moderna and Pfizer vaccines, each with estimated efficacy of about 95%). Simulations are carried out using the baseline parameter values in Tables 3, 4, with various values of the vaccination coverage parameter (ξ_v). For these simulations, parameters and initial conditions related to the other intervention (face mask usage) are maintained at their baseline values. Since the Moderna and Pfizer vaccines have essentially the same estimated efficacy (≈95%), we group them together in the numerical simulations for this section.

The simulation results obtained for the Moderna and Pfizer vaccine, depicted in Figures 6A–C, show that, in the absence of vaccination (and with social-distancing at baseline compliance level), approximately 511, 100 cumulative deaths will be recorded in the U.S. by April 10, 2021 (blue curves of Figures 6A–C). Furthermore, this figure shows a marked reduction in daily mortality with increasing vaccination coverage (ξ_v). This reduction further increases if vaccination is combined with social-distancing. For instance, with social-distancing compliance maintained at its baseline value on October 12, 2020 (i.e., c_s = 0), vaccinating at a rate of 0.00074 per day (which roughly translates to vaccinating 250, 000 people every day) resulted in a reduction of the projected cumulative mortality recorded by April 10, 2021 by 12%, in comparison to the case when no vaccination is used (magenta curve in Figure 6A, in comparison to the blue curve of the same figure). In fact, up to 31% of the projected cumulative mortality to be recorded by April 10, 2021 could be averted if, for this vaccination rate, the baseline social-distancing compliance is increased by 10% (i.e., c_s = 0.1; magenta curve in Figure 6C, in comparison to magenta curve in Figure 6A). If the vaccination rate is further increased to, for instance, ξ_v = 0.0015 per day (corresponding to vaccinating about 500, 000 people every day), while keeping social-distancing at its baseline compliance level (i.e., c_s = 0), our simulations show a reduction of 27% in the projected cumulative mortality by April 10, 2021, in comparison to the baseline social-distancing scenario (gold curve, Figure 6A, in comparison to the blue curve of the same figure). This reduction increases to 38% if the vaccination program is supplemented with social-distancing that increases the baseline compliance by 10% (gold curve, Figure 6C). If 1 million people are vaccinated per day (i.e., ξ_v = 0.003 per day), our simulations show that the use of the Moderna and Pfizer vaccines could lead to up to 36% reduction in the projected cumulative mortality by April 10, 2021 in the U.S. if the vaccination program is combined with a 10% increase in social-distancing compliance level (green curve of Figure 6C). Finally, compared to the Moderna and Pfizer vaccines, slightly lower reductions in the projected cumulative mortality are recorded when the AstraZeneca vaccine (with moderate to high vaccination coverage) is used (Figures 6D–F), particularly if combined with social-distancing. These results are summarized in Table 6.

FIGURE 6

Figure 6. Assessment of the combined impact of vaccination and social-distancing on cumulative mortality. Simulations of the model (2.5)–(2.6), depicting cumulative as a function of time, for the three vaccines considered in this study and various levels of increases in baseline social-distancing compliance starting from October 12, 2020 (c_s). (A–C) Pfizer or Moderna vaccine. (D–F) AstraZeneca vaccine. The vaccination rates ${\xi }_v=7.4\times 10^{-4}$, 1.5 ×10⁻³ per day, and 3.0 ×10⁻³ per day correspond, respectively, to vaccinating ~2.5 ×10⁵, 5.0 ×10⁵ and 1.0 ×10⁶ people per day. Other parameter values of the model used are as presented in Tables 3, 4.

TABLE 6

Table 6. Percentage reduction in projected cumulative COVID-19 mortality on April 10, 2021, in relation to the cumulative mortality in the absence of vaccination (511, 100 COVID-19 deaths on April 10, 2021), for the three vaccines considered in this study: AstraZeneca vaccine (efficacy ε_v = 0.7); Pfizer and/or Moderna vaccine (efficacy ε_v = 0.95), and various levels of increases in baseline social-distancing compliance attained on October 12, 2020 (c_s) and vaccination rate (ξ_v). SD, social-distancing compliance.

4.4. Combined Impact of Vaccination and Social-Distancing on Time-to-Elimination

The model (2.5)–(2.6) will now be simulated to assess the population-level impact of the combined vaccination and social-distancing interventions on the expected time the pandemic might be eliminated in the U.S. if the two strategies are implemented together. Mathematically, we define “elimination" to mean when the number of daily new cases is identically zero. As in section 4.3, we consider the three vaccines (AstraZeneca, Moderna, and the Pfizer vaccines), and assume that the vaccination program was started on October 12, 2020. The model is simulated to generate a time series of new daily COVID-19 cases in the U.S., for various vaccination rate (ξ_v) and levels of increases in baseline social-distancing compliance (c_s).

The results obtained, for each of the three vaccines, are depicted in Figure 7. This figure shows a marked decrease in disease burden (measured in terms of the number of new daily cases), with the possibility of elimination of the pandemic within 8–10 months from the commencement of the vaccination program. In particular, these simulations show that vaccinating 250, 000 people per day, with the Moderna or the Pfizer vaccine, will result in COVID-19 elimination in the U.S. by mid August of 2021, if the social-distancing compliance is kept at its current baseline compliance level (blue curve of Figure 7). For this scenario, the elimination will be reached in late August 2021 using the AstraZeneca vaccine. If the vaccination rate is further increased, such as to vaccinating 1 million people every day (and keeping social-distancing at its October 12, 2020 baseline), COVID-19 elimination is achieved much sooner in the United States. For instance, for this scenario (i.e., with ξ_v = 0.003 per day), the pandemic can be eliminated by late June of 2021 using the Moderna or the Pfizer vaccines (green curve of Figure 7A) and by mid July of 2021 using the AstraZeneca vaccine (blue curve of Figure 7D).

FIGURE 7

Figure 7. Effect of vaccination and social-distancing on time-to-elimination. Simulations of the model (2.5)–(2.6), depicting the impact of three vaccines against COVID-19 (the AstraZeneca vaccine, and the Pfizer or Moderna vaccine) and social-distancing, on time-to-elimination of the pandemic in the U.S. (A–C) Moderna or Pfizer vaccines. (D–F) AstraZeneca vaccine. The social-distancing compliance is baseline for (A,D), c_s = 0.05 for (B,E), and c_s = 0.10 for (E,F). The vaccination rates ${\xi }_v=7.4\times 10^{-4},1.5\times 10^{-3},3.0\times 10^{-3}$ per day correspond, respectively, to vaccinating approximately 2.5 ×10⁵, 5.0 ×10⁵, 1.0 ×10⁶ people per day. The values of the other parameters of the model used in the simulation are as given in Tables 3, 4.

Our simulations further show that if the vaccination program is combined with social-distancing that increases the baseline compliance by 10%, COVID-19 can be eliminated in the U.S. by as early as the end of May of 2021 using the Moderna or the Pfizer vaccine (green curve of Figure 7C), and by late June of 2021 using the AstraZeneca vaccine (green curve, Figure 7F). In conclusion, these simulations show that any of the three vaccines considered in this study will lead to the elimination of the pandemic in the U.S. if the vaccination rate is moderately-high enough. The time-to-elimination depends on the vaccination rate and the level of increases in the baseline social-distancing compliance attained by October 12, 2020. The pandemic can be eliminated as early as the end of May of 2021 if moderate to high vaccination rate (e.g., 1 million people are vaccinated per day) and social-distancing compliance (e.g., c_s = 0.1) is attained and maintained.

It is worth mentioning that two of the three vaccines that are currently being used in the U.S. were only approved by the FDA in December 2020 (the Pfizer vaccine was approved on December 11, 2020, while the Moderna vaccine was approved a week later), and their administration into the arms of Americans started late in December 2020. Therefore, as we noted earlier, the greater U.S. community might only be able to receive any of the vaccines by March or April 2021 (we chose March 15, 2021 as our reference point for simulation/comparative purposes). Thus, with a mass vaccination start date of mid March 2021 (i.e., if we can only achieve vaccinating 1 million or more people daily from mid March 2021), then COVID-19 elimination, assuming a 10% increase in baseline social-distancing compliance achieved on October 12, 2020, can be achieved by the end of October 2021 using the Moderna or the Pfizer vaccine (for the AstraZeneca vaccine, elimination will extend to November of 2021). It should be mentioned that the elimination can be achieved even earlier if large scale community vaccination in the U.S. is started earlier than our projected March 15, 2021, and particularly if this (early large scale vaccination before March 15, 2021) is also complemented with significant increase in baseline social-distancing compliance (such as increasing the baseline compliance by 10%).

In summary, our study clearly shows that the prospect of eliminating COVID-19 in the U.S. by the middle or early fall of 2021 is very much feasible if moderate level of coverage can be achieved using either of the two vaccines being used in the U.S., and if this vaccination coverage is complemented with a social-distancing strategy that increases the baseline compliance achieved by October 12, 2020 by a mere 10%. Our study certainly points to the fact that we will be seeing the back of this devastating Coronavirus beast, and socio-economic life may return to near normalcy, in 2021.

4.5. Assessing the Impacts of Waning Immunity, Mask Fatigue, and Relaxation of Mask Mandate for Fully-Vaccinated Individuals, and Therapeutic Benefits of Vaccines

In this section, the multi-group model (2.5)–(2.6) will be adapted and simulated to assess the population-level impact of three other factors that may significantly affect the effectiveness of the vaccination program against COVID-19, namely (a) waning natural and vaccine-derived immunity (44–46), (b) mask fatigue (and giving up masking) by fully-vaccinated individuals (47), and (c) therapeutic benefits of the vaccines (such as reducing development of severe disease, hospitalization and mortality in breakthrough infections, as well as in reducing transmissibility of infected vaccinated individuals) (18, 19, 48). Although the model (2.5)–(2.6) does not explicitly incorporate the aforementioned factors, it can readily be adapted to allow for their assessment. We describe below how the model can be adapted to achieve this objective, in addition to illustrating the effects of the factors via numerical simulations of the resulting adapted version of the model (2.5)–(2.6). For consistency, the simulations in this section will also be carried out from the beginning of the third pandemic wave in the U.S. (i.e., from October 12, 2020).

4.5.1. Waning Natural and Vaccine-Derived Immunity

Waning natural immunity can be incorporated in the model by allowing a transition from the compartment of recovered individuals (for each of the two groups) into the corresponding compartment for unvaccinated susceptible individuals (i.e., the immunity derived from natural recovery from COVID-19 infection ultimately wanes, and the recovered individuals subsequently become wholly-susceptible again). To adapt the model to account for this, we introduce a new parameter, ω_r, to represent the per capita rate at which recovered individuals revert to the corresponding unvaccinated susceptible compartment (i.e., the quantity ω_rR_i, with i = 1, 2) is subtracted from the equation for R_i and added to the corresponding equation for S_iu in the model (2.5)–(2.6).

Similarly, vaccine-derived waning immunity can be incorporated into the model (2.5)–(2.6) by allowing for transitions from the vaccinated susceptible compartments (S_iv; i = 1, 2) to the corresponding unvaccinated susceptible compartment (S_iu; i = 1, 2). We introduce a new parameter, ω_v, to represent the rate of waning of vaccine-derived immunity. To incorporate this into the model, the quantity ω_vS_iv (i = 1, 2) is subtracted from the equation for S_iv and added to the corresponding equation for S_iu in the model (2.5)–(2.6). For simulation purposes, we set ω_r and ω_v to be 1/270 per day and 1/180 per day (44, 46), respectively (corresponding to a 9 and 6 months duration for the waning of natural and vaccine-derived immunity, respectively).

The model (2.5)–(2.6) is now simulated, using the parameter values in Tables 3, 4, together with the above modifications (accounting for waning natural and vaccine-derived immunity, using the estimated values of ω_r and ω_v), to assess the potential impact of waning immunity on the COVID-19 dynamics in the United States. The results obtained, depicted in Figure 8A, show a slight increase in the peak number of new daily cases, in comparison to the results in Figure 7A, where the effect of waning immunity was not considered. In particular, if the vaccination rate is 250, 000 per day (i.e., if ξ_v is set at ${\xi }_v=7.4\times 10^{-4}$ per day), then the peak number of new cases increases by approximately 2% (in comparison to the case where no waning immunity is considered), and the time-to-elimination of the pandemic increases by about 13 days (compare blue curves in Figures 7A, 8A). If the daily vaccination rate is increased to one million per day (i.e., if ${\xi }_v=3.0\times 10^{-3}$ per day), then the peak new cases increases by up to 6% (in comparison to the case with no waning immunity) and the time-to-elimination increases by about a month (compare green curves in Figures 7A, 8A). The increases in burden and time-to-extinction in this case (with 1 million vaccinated daily, in comparison to the case with 250,000 people getting vaccinated daily) is due to the fact waning of both natural and vaccine-derived immunity causes a corresponding increase in the pool of susceptible individuals who can acquire infection (thereby increasing number of new cases and extending time-to-elimination). Thus, these simulations show that waning of natural and vaccine-derived immunity cause only a marginal increase in the burden and time-to-elimination of the pandemic.

FIGURE 8

Figure 8. Effect of (A) vaccine-induced and natural immunity waning, (B) unmasking by vaccinated individuals, and (c) therapeutic benefits of vaccines on the burden of the pandemic and the time-to-elimination. The vaccination rates ${\xi }_v=7.4\times 10^{-4},1.5\times 10^{-3},3.0\times 10^{-3}$, per day correspond, respectively, to vaccinating approximately 2.5 ×10⁵, 5.0 ×10⁵, 1.0 ×10⁶ people per day. The vaccine-induced and natural immunity waning rate parameters, ω_v and ω_r, are set to ω_v = 1/180 per day (44) and ω_r = 1/270 per day (46), respectively. The effect of therapeutic benefits of the vaccine depicted in (C) is incorporated by reducing the baseline values of the parameters (from Tables 3, 4) that are related to development of severe disease (r), hospitalization (ϕ_jI, j = 1, 2) and mortality (δ_jI, δ_jH, j = 1, 2) by 5%, in addition to increasing the baseline value of the parameters related to the recovery rate (γ_jI, γ_jA, γ_jH, j = 1, 2) by 5%. The values of the other parameters of the model used in the simulation are as given in Tables 3, 4.

4.5.2. Mask Fatigue and Relaxation of Mask Mandates for Fully-Vaccinated Individuals

To incorporate the effect of mask fatigue, or relaxation of mask mandates (47), in fully-vaccinated individuals into the model (2.5)–(2.6), we consider the worst-case scenario where all fully-vaccinated individuals opt to give up masking in public. To account for the worst case scenario of this (i.e., the case in which every fully-vaccinated individual abandons masking) in the model, we remove the state variable S_2v, for the vaccinated susceptible individuals in the mask-wearing group 2, from the model. Further, we re-direct all the new vaccinated individuals from group 2 into the vaccinated class of the non-masking group 1 (i.e., we add the term ξ_vS_2u from the equation for the rate of change of the S_2u population to that for the rate of change of the S_1v population, and the equation for S_2v is removed from the model) and also remove the term −α₁₂S_1v from the equation for the rate of change of the S_1v population (to ensure that vaccinated individuals in group 1 do not move to the mask-wearing group 2). Simulations of the model (2.5)–(2.6), under this setting (and using the parameter values in Tables 3, 4), depicted in Figure 8B, show a marginal change in the peak number of new cases and the time-to-elimination, in comparison to the case when fully-vaccinated individuals do not completely give up masking (i.e., compare Figure 8B with Figure 7A).

4.5.3. Therapeutic Benefits of COVID-19 Vaccines

Result from recent clinical trials have shown very promising therapeutic benefits for both the Pfizer and Moderna vaccines (18, 19). In this section, we seek to use the multi-group model (2.5)–(2.6) to assess the impact of such benefits on the dynamics of the disease in the United States. Since the model does not explicitly stratify the population of infected individuals according to whether they are vaccinated or not, a number of factors will come into play when estimating the overall impact of the therapeutic benefits, such as the high efficacy of the two vaccines (~95%, thereby significantly reducing the size of breakthrough infections), level of vaccine hesitancy in the community and the current daily infection rate in the community. Taking all these into account, we consider it plausible, as a first approximation, to estimate the overall therapeutic benefits in the U.S., at the beginning of the third wave characterized by low vaccination coverage (December 2020 until about February 2021), high disease burden (skyrocketing number of reported confirmed cases, hospitalizations and COVID-19 mortality), by a 5% reduction in severe or symptomatic illness, breakthrough transmission, hospitalization, and mortality, as well as a 5% increase in the rate of recovery from infection for vaccinated infected individuals. In other words, the effect of therapeutic benefits of the vaccine is incorporated into our model by reducing the baseline values of the parameters related to development of severe disease (r), hospitalization (ϕ_jI, j = 1, 2) and mortality δ_jI and δ_jH with j = 1, 2) by 5%, in addition to increasing the baseline value of the parameter related to the recovery rate (γ_jI, γ_jA and γ_jH, with j = 1, 2). The simulation results obtained, for this hypothetical scenario, show a marked reduction in disease burden and a decrease in time-to-elimination (Figure 8C), in comparison to the case where such therapeutic benefits are not accounted for (Figure 7A). In particular, if one million people are vaccinated daily (i.e., if the vaccination rate is set at ${\xi }_v=3.0\times 10^{-3}$ per day), up to 37% decrease in the peak number of new cases could be achieved. Further, the time-to-elimination decreases by 17 days (compare green curves in Figures 7A, 8A). Higher reductions in disease burden, and more accelerated time-to-elimination, will be achieved if higher percentages of therapeutic benefits are assumed. It should be mentioned that a more rigorous way to introduce the impact of therapeutic benefits into the multi-group model (2.5)–(2.6) will be to further restructure the infected compartments of the model in terms of whether they are vaccinated or not (doing so will result in a model with at least 28 nonlinear differential equations, which may be difficult to track mathematically and statistically).

In summary, it is shown in this section (based on the parameter values used in our simulations) that, while waning natural and vaccine-derived immunity generally induces a relatively small increase in the burden of the pandemic, together with a correspondingly marginal increase in the time-to-elimination (in comparison to the case when these effects are not incorporated into the model), the therapeutic benefits of the vaccines offer a dramatic impact on the trajectory of the disease (by significantly reducing both the burden and time-to-elimination of the pandemic, in comparison to the case when such benefits are not accounted for in the model). Finally, it is worth stating that, although the simulations carried out in section 4.5 are for the Pfizer and Moderna COVID-19 vaccines only (illustrated in Figure 7), similar simulations can also be carried out for the AstraZeneca and other vaccines with lower preventive effective efficacies. These simulations will, of course, show higher disease burden (owing to their reduced efficacy), in comparison to the case when Pfizer and Moderna vaccines are used.

5. Discussion and Conclusions

Since its emergence late in December of 2019, the novel Coronavirus pandemic has inflicted devastating public health and economic burden across the world. As of January 24, 2021, the pandemic accounted for over 100 million confirmed cases and 2.1 million fatalities globally (the United States accounted for 25, 123, 857 confirmed cases and 419, 204 deaths). Although control efforts against the pandemic have focused on the use of non-pharmaceutical interventions, such as social-distancing, face mask usage, quarantine, self-isolation, contact-tracing, community lockdowns, etc., a number of highly-efficacious and safe anti-COVID-19 vaccines have been developed and approved for use in humans. In particular, two of the three FDA-EUA vaccines (manufactured by Moderna Inc. and Pfizer Inc.) have estimated protective efficacy of about 95%. Furthermore the AstraZeneca vaccine, developed by the pharmaceutical giant, AstraZeneca and University of Oxford has protective efficacy of 70%. Mathematics (modeling, analysis, and data analytics) has historically been used to provide robust insight into the transmission dynamics and control of infectious diseases, dating back to the pioneering works of the likes of Daniel Bernoulli in the 1760s (on smallpox immunization), Sir Ronald Ross and George Macdonald between the 1920s and 1950s (on malaria modeling) and the compartmental modeling framework developed by Kermack and McKendrick in the 1920s (49–51). In this study, we used mathematical modeling approaches, coupled with rigorous analysis, to assess the potential population-level impact of the use of the three vaccines and data analytics in curtailing the burden of the COVID-19 pandemic in the U.S. We have also assessed the impact of other non-pharmaceutical interventions, such as face mask and social-distancing, implemented singly or in combination with any of the three vaccines, on the dynamics and control of the pandemic.

We developed a novel mathematical model, which stratifies the total population into two subgroups of individuals who habitually wear face masks in public and those who do not. The resulting two group COVID-19 vaccination model, which takes the form of a deterministic system of nonlinear ordinary differential equations, was initially fitted using observed cumulative COVID-induced mortality data for the U.S. Specifically, we fitted the model with the cumulative mortality data corresponding to the period when the U.S. was experiencing the third wave of the COVID-19 pandemic (estimated to have started around October 12, 2020). In addition to allowing for the assessment of the population-level of each of the three currently-available vaccines, the model also allows for the assessment of the initial size of the population of individuals who habitually wear face masks in public, as well as assessing the impact of increase in the baseline social-distancing compliance attained as of October 12, 2020. After the successful calibration of the model, we carried out rigorous asymptotic stability analysis to gain insight into the main qualitative features of the model. In particular, we showed that the disease-free equilibrium of the model is locally-asymptotically stable whenever a certain epidemiological threshold, known as the control reproduction number (denoted by ${{R}}_c$), is less than unity. The implication of this result is that (for the case when ${{R}}_c<1$), a small influx of COVID-infected individuals will not generate a alrge outbreak in the community.

The expression for the reproduction number (${{R}}_c$) was used to compute the nationwide vaccine-induced herd immunity threshold for a special case of the model where change of masking behavior is not allowed (i.e., where mask wearers remain wearers and non-wearers do not become wearers). The herd immunity threshold represents the minimum proportion of the susceptible U.S. population that needs to be vaccinated to ensure elimination of the pandemic. Simulations of our model showed, for the current baseline level of social-distancing in the U.S. (and baseline level of initial size of the population of face masks wearers), herd immunity can be achieved in the U.S. using the AstraZeneca vaccine if at least 80% of the susceptible population is fully vaccinated. The threshold herd immunity level needed when either the Pfizer or Moderna vaccine is used reduces to 59%. Our simulations further showed that the level of herd immunity needed to eliminate the pandemic decreases, for each of the three vaccines, with increasing levels of baseline social-distancing compliance. In particular, the baseline social-distancing achieved at the beginning of our simulation period (i.e., the level of social-distancing in the U.S. as of October 12, 2020) is increased by 10%, the herd immunity requirement for the AstraZeneca or Pfizer/Moderna vaccine reduced, respectively, to 73 and 54%. Furthermore, if the baseline social-distancing is increased by 30%, the herd immunity threshold needed to eliminate the pandemic using the AstraZeneca or Pfizer/Moderna vaccine reduced to a mere 53 and 39%, respectively. In other words, this study showed that the prospect of achieving vaccine-derived herd immunity, using any of the three vaccines considered in this study, is very promising, particularly if the vaccination program is complemented with increased levels of baseline social-distancing.

The multigroup nature of the model we developed in this study, where the total population is stratified into the two groups of those who habitually face mask in public and those who do not (with back-and-forth transitions between the two groups allowed), enabled us to assess the population-level impact of the initial sizes of the two groups in curtailing the spread of the pandemic in the United States. We assessed this by simulating the model during the beginning of the third wave of the pandemic in the U.S. (starting from October 12, 2020), and used the proportion of masks-wearers embedded in the cumulative mortality data we used to fit the model as the baseline. Our study emphasized the fact that early adoption of mask mandate plays a major role in effectively reducing the burden (as measured in terms of cumulative mortality) of the pandemic. This effect is particularly more pronounced when individuals in the face masks-wearing group do not change their behavior and transition to the non-mask wearing group (and non-mask wearers adopt a masks-wearing habit). Our study further showed that, for the case where the aforementioned back-and-forth transitions between the masks-wearing and the non-mask wearing groups are allowed, there is a threshold level of the initial size of the proportion of face masks-wearers above which the disease burden will be reduced, below which the disease burden actually increases. Our study estimated this threshold value of the initial size of the masks-wearing group to be about 50%. The epidemiological implication of this result is that the continued implementation of face masks use strategy (particularly at the high initial coverage level) will be highly beneficial in effectively curtailing the pandemic burden between now and the time when the two FDA-EUA vaccines become widely available to the general public in the U.S. (expected to be around mid March to mid April of 2021).

We further showed that the time-to-elimination of COVID-19 in the U.S., using a vaccine (and a non-pharmaceutical intervention), depended on the daily vaccination rate (i.e., number of people vaccinated per day) and the level of increase in baseline social-distancing compliance achieved at the onset of the third wave of the pandemic (October 12, 2020). Specifically, our study showed that the COVID-19 pandemic can be eliminated in the U.S. by early May of 2021 if we are able to achieve moderate level of daily vaccination rate (such as vaccinating 1 million people every day) and the baseline social-distancing compliance achieved on October 12, 2020 is increased by 10% (and sustained). It should, however, be mentioned that the time-to-elimination is sensitive to the level of community transmission of COVID-19 in the population (it is also sensitive to the effectiveness and coverage (compliance) levels of the other (non-pharmaceutical) interventions, particularly face mask usage and social-distancing compliance, implemented in the community). Specifically, our study was carried out between December 2020 and January 2021, when the United States was experiencing a devastating third wave of the COVID-19 pandemic (recording on average over 200, 000 confirmed cases per day, together with record numbers of hospitalizations and COVID-induced mortality). This explains the somewhat longer estimated time-to-elimination of the pandemic, using any of the three vaccines considered in this study, for the case where social-distancing compliance is kept at the baseline level. The estimate for the time-to-elimination (using any of the vaccines considered in this study) will be shorter if the community transmission is significantly reduced (as will be vividly evident from the reduced values of the transmission-related and mortality-related parameters of the re-calibrated version of our model).

It is worth emphasizing that at the time this study was carried out (between December 2020 and January 2021), it was unclear whether natural or vaccine-induced immunity to COVID-19 waned over time. It was also unclear whether the then new vaccines that received FDA's Emergency Use Authorization (Pfizer and Moderna) offer therapeutic benefits (such as reducing severe disease, hospitalization and deaths, in addition to accelerating recovery rate in vaccinated infected individuals). However, by the time we are reviewing the manuscript (June 2021), new data and studies have provided clarity on waning immunity to COVID-19 (45, 46) and on the therapeutic benefits of some of the COVID-19 vaccines (18, 19). Furthermore, the U.S. Centers for Disease Control and Prevention has modified its guidelines on masking, allowing fully-vaccinated individuals not to wear masks under certain circumstances (47). Consequently, we adapted the multi-group model we developed to allow for the assessment of the aforementioned new facts associated with the COVID-19 dynamics. Specifically, we adapted and simulated the model to assess the impact of waning immunity (both natural and vaccine-derived), mask fatigue and relaxation of mask mandates for fully-vaccinated individuals and the therapeutic benefits of the FDA-authorized vaccines on the disease burden (measured in terms of peak daily cases) and time-to-elimination of the pandemic in the United States. The simulations were carried out for the hypothetical scenario that the vaccination program was started at the beginning of the third wave of the pandemic in the U.S. (i.e., in October of 2020). Since the vaccines were not available until December 2020, and large scale vaccine rollout was only achieved some time in end of March 2021 or early April 2021, we adapted our conclusions appropriately to account for this time lag. Consequently, our simulation results, for these settings, show that, while waning natural and vaccine-derived immunity induces only a relatively marginal increase in both the burden and time-to-elimination of the pandemic, incorporating therapeutic benefits of the vaccine into the model causes a dramatic reduction in both the burden and time-to-elimination. If the impacts of therapeutic benefits are incorporated into the model from the very beginning of the third wave of the pandemic (October 2020), our simulations show that the pandemic could theoretically be eliminated in the U.S. by as early as late May, 2021 (note that, in the absence of such therapeutic benefits, our study estimated the time-to-elimination to be some time in October, 2021).

In summary, our study suggest that the prospects of COVID-19 elimination in the U.S. is very promising, using any of the three vaccines considered in this study. The elimination prospects are greatly enhanced if the therapeutic benefits of the FDA-EUA vaccines are incorporated into the multi-group model we developed and used in this study. While, for the baseline scenario, the AstraZeneca vaccine requires at least 80% of the U.S. population to be fully vaccinated to achieve herd immunity (needed for the elimination of the pandemic), such herd immunity can be achieved using any of the two FDA-EUA vaccines (Pfizer or Moderna) if only 59% of Americans are fully vaccinated. The prospects of eliminating COVID-19 using any of the three vaccines is greatly enhanced if the vaccination program is combined with a social-distancing strategy that increases the baseline compliance level of the social-distancing attained during the beginning of the third wave of the COVID-19 pandemic in the United States. In fact, our simulations strongly suggested that COVID-19 could have been eliminated in the U.S. as early as May 2021, depending on the level of increase in baseline social-distancing compliance. In other words, if we can continue to maintain social-distancing, while large scale vaccination is being implemented, COVID-19 can be history, and life can begin to return to normalcy or near-normalcy, in the summer or fall of 2021.

Some of the limitations of our model include not explicitly accounting for some important heterogeneities, such as age-structure, risk-structure, and vaccine dose structure, and the impacts of SARS-CoV-2 variants. Accounting for these may alter our results, especially during the early days of the vaccine administration (e.g., from December 2020 to April 2021) when the vaccine doses were generally in limited supply and needed to be prioritized to high-risk groups. Hence, our simulation results and conclusions should be interpreted with these limitations in mind. Further, while our multi-group model did not explicitly account for some other factors potentially relevant to COVID-19 dynamics, such as waning of natural and vaccine-derived immunity to COVID-19, mask fatigue and relaxation of mask mandates for fully-vaccinated individuals and the impacts of therapeutic benefits of the approved vaccines, our multi-group model was robust enough to allow for the assessment of these factors. We showed that, while incorporating waning of immunity and mask fatigue and relaxation of mask mandates in fully-vaccinated individuals in the model we developed only caused a marginal increase in disease burden and time-to-elimination, incorporating the impacts of the therapeutic benefits of the approved vaccines (even at a relatively low overall rate) resulted in a dramatic reduction in both the disease burden and time-to-elimination of the pandemic.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found at: https://coronavirus.jhu.edu/map.html.

Author Contributions

AG and CN conceived the study and developed the model. AG, EI, and CN carried out the analysis. All authors wrote and edited the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

One of the authors (AG) acknowledge the support, in part, of the Simons Foundation (Award #585022) and the National Science Foundation (Grant Number: DMS-2052363). CN acknowledges the support of the Simons Foundation (Award #627346). GN acknowledge the grants and support of the Cameroon Ministry of Higher Education, through the initiative for the modernization of research in Higher Education. All authors wish to express their deepest sympathy to the families of the victims of the SARS-CoV-2, and extend profound appreciation to the frontline healthcare workers for their heroic effort and sacrifices to save the lives of others. The authors are very grateful to the two reviewers for their very constructive comments, which have significantly enhanced the quality and clarity of the manuscript.

References

1. Center for Systems Science and Engineering at Johns Hopkins University. COVID-19; 2020. Available online at: https://github.com/CSSEGISandData/COVID-19 (accessed June 16, 2021).

Google Scholar

2. Centers for Disease Control and Prevention. Scientific Brief: SARS-CoV-2 Transmission. CDC Information. (2021) Available online at: https://www.cdc.gov/coronavirus/2019-ncov/science/science-briefs/sars-cov-2-transmission.html (accessed June 12, 2021).

Google Scholar

3. Centers for Disease Control and Prevention. Coronavirus Disease 2019 (COVID-19). National Center for Immunization and Respiratory Diseases (NCIRD), Division of Viral Diseases (2020). Available online at: https://www.cdc.gov/coronavirus/2019-ncov/index.html (accessed March 4, 2020).

Google Scholar

4. Ngonghala CN, Iboi E, Eikenberry S, Scotch M, MacIntyre CR, Bonds MH, et al. Mathematical assessment of the impact of non-pharmaceutical interventions on curtailing the 2019 novel Coronavirus. Math Biosci. (2020) 325:108364. doi: 10.1016/j.mbs.2020.108364

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Centers for Disease Control and Prevention. Different COVID-19 Vaccines. CDC Information. (2021) Available online at: https://www.cdc.gov/coronavirus/2019/ncov/vaccines/different/vaccines.html (accessed January 25, 2021).

Google Scholar

6. US Food and Drug Administration. FDA Takes Key Action in Fight Against COVID-19 By Issuing Emergency Use Authorization for First COVID-19 Vaccine. FDA Office of Media Affairs (2021). Available online at: https://www.fda.gov/news-events/press-announcements/fda-takes-key-action-fight-against-covid-19-issuing-emergency-use-authorization-first-covid-19 (accessed January 25, 2021).

Google Scholar

7. Eikenberry SE, Muncuso M, Iboi E, Phan T, Kostelich E, Kuang Y, et al. To mask or not to mask: modeling the potential for face mask use by the general public to curtail the COVID-19 pandemic. Infect Dis Model. (2020) 5:293–308. doi: 10.1016/j.idm.2020.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Ngonghala CN, Iboi E, Gumel AB. Could masks curtail the post-lockdown resurgence of COVID-19 in the US? Math Biosci. (2020) 329:108452. doi: 10.1016/j.mbs.2020.108452

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Iboi EA, Ngonghala CN, Gumel AB. Will an imperfect vaccine curtail the COVID-19 pandemic in the US? Infect Dis Model. (2020) 5:510–24. doi: 10.1016/j.idm.2020.07.006

CrossRef Full Text | Google Scholar

10. Food and Drug Administration. FDA Issues Emergency Use Authorization for Third COVID-19 Vaccine. FDA Office of Media Affairs (2021). Available online at: https://www.fda.gov/news-events/press-announcements/fda-issues-emergency-use-authorization-third-covid-19-vaccine (accessed June 25, 2021).

Google Scholar

11. Pfizer. fizer and BioNTech to Submit Emergency Use Authorization Request Today to the U.S. FDA for COVID-19 Vaccine. (2020). Available online at: https://www.pfizer.com/news/press-release/press-release-detail/pfizer-and-biontech-submit-emergency-use-authorization (accessed June 16, 2021).

Google Scholar

12. National Institute of Health. Promising Interim Results from Clinical Trial of NIH-Moderna COVID-19 Vaccine. (2020). Available online at: https://www.nih.gov/news-events/news-releases/promising-interim-results-clinical-trial-nih-moderna-covid-19-vaccine (accessed June 16, 2021).

Google Scholar

13. AstraZeneca. AZD1222 Vaccine Met Primary Efficacy Endpoint in Preventing COVID-19. (2020). Available online at: https://www.astrazeneca.com/media-centre/press-releases/2020/azd1222hlr.html (accessed June 16, 2021).

Google Scholar

14. Graham Lawton. Everything You Need to Know About the Pfizer/BioNTech COVID-19 Vaccine. (2020). Available online at: https://www.newscientist.com/article/2261805-everything-you-need-to-know-about-the-pfizer-biontech-covid-19-vaccine/ (accessed June 16, 2021).

Google Scholar

15. Moderna. Moderna Announces Longer Shelf Life for Its COVID-19 Vaccine Candidate at Refrigerated Temperatures. (2020). Available online at: https://investors.modernatx.com/news-releases/news-release-details/moderna-announces-longer-shelf-life-its-covid-19-vaccine (accessed June 16, 2021).

Google Scholar

16. Centers for Disease Control and Prevention (CDC). Different COVID-19 Vaccines. CDC Information. (2021). Available online at: https://www.cdc.gov/coronavirus/2019-ncov/vaccines/different-vaccines.html?s_cid=11304:johnson/and/johnson/covid/vaccine:sem.ga:p:RG:GM:gen:PTN:FY21 (accessed June 25, 2021).

Google Scholar

17. Pearson S. What is the Difference Between the Pfizer, Moderna, and Johnson & Johnson COVID-19 Vaccines? GoodRx (2021). Available online at: https://www.goodrx.com/blog/comparing-covid-19-vaccines/ (accessed June 25, 2021).

Google Scholar

18. Centers for Disease Control and Prevention. Benefits of Getting a COVID-19 Vaccine. CDC Information. (2021) Available online at: https://www.cdc.gov/coronavirus/2019-ncov/vaccines/vaccine-benefits.html (accessed June 11, 2021).

Google Scholar

19. Dagan N, Barda N, Kepten E, Miron O, Perchik S, Katz MA, et al. BNT162b2 mRNA COVID-19 vaccine in a nationwide mass vaccination setting. N Engl J Med. (2021) 384:1412–23. doi: 10.1056/NEJMoa2101765

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Srivastava A, Chowell G. Understanding spatial heterogeneity of COVID-19 pandemic using shape analysis of growth rate curves. medRxiv [Preprint]. (2020). doi: 10.1101/2020.05.25.20112433

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Gumel AB, Iboi EA, Ngonghala CN, Elbasha EH. A primer on using mathematics to understand COVID-19 dynamics: modeling, analysis and simulations. Infect Dis Model. (2021) 6:148–68. doi: 10.1016/j.idm.2020.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Schneider KA, Ngwa GA, Schwehm M, Eichner L, Eichner M. The COVID-19 pandemic preparedness simulation tool: CovidSIM. BMC Infect Dis. (2020) 20:1–11. doi: 10.2139/ssrn.3578789

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ngonghala CN, Goel P, Kutor D, Bhattacharyya S. Human choice to self-isolate in the face of the COVID-19 pandemic: a game dynamic modelling approach. J Theor Biol. (2021) 521:110692. doi: 10.1016/j.jtbi.2021.110692

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Hellewell J, Abbott S, Gimma A, Bosse NI, Jarvis CI, Russell TW, et al. Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts. Lancet Glob Health. (2020) 8:E488–96. doi: 10.1016/S2214-109X(20)30074-7

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kucharski AJ, Russell TW, Diamond C, Liu Y, Edmunds J, Funk S, et al. Early dynamics of transmission and control of COVID-19: a mathematical modelling study. Lancet Infect Dis. (2020) 20:553–8. doi: 10.1016/S1473-3099(20)30144-4

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Xue L, Jing S, Miller JC, Sun W, Li H, Estrada-Franco JG, et al. A data-driven network model for the emerging COVID-19 epidemics in Wuhan, Toronto and Italy. Math Biosci. (2020) 326:108391. doi: 10.1016/j.mbs.2020.108391

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Ferguson NM, Laydon D, Nedjati-Gilani G, Imai N, Ainslie K, Baguelin M, et al. Impact of Non-Pharmaceutical Interventions (NPIs) to Reduce COVID-19 Mortality and Healthcare Demand. London: Imperial College COVID-19 Response Team (2020).

PubMed Abstract | Google Scholar

28. Centers for Disease Control and Prevention (CDC). CDC Director's Statement on Pfizer's Use of COVID-19 Vaccine in Adolescents Age 12 and Older. CDC information (2021). Available online at: https://www.cdc.gov/media/releases/2021/s0512-advisory-committee-signing.html (accessed June 25, 2021).

Google Scholar

29. Banks HT, Davidian M, Samuels JR, Sutton KL. In: An Inverse Problem Statistical Methodology Summary. Dordrecht: Springer Netherlands (2009). p. 249–302. Available online at: https://doi.org/10.1007/978-90-481-2313-1_11

Google Scholar

30. Chowell G. Fitting dynamic models to epidemic outbreaks with quantified uncertainty: a primer for parameter uncertainty, identifiability, and forecasts. Infect Dis Model. (2017) 2:379–98. doi: 10.1016/j.idm.2017.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhou C. Evaluating new evidence in the early dynamics of the novel coronavirus COVID-19 outbreak in Wuhan, China with real time domestic traffic and potential asymptomatic transmissions. medRxiv [Preprint]. (2020). doi: 10.1101/2020.02.15.20023440

CrossRef Full Text | Google Scholar

32. Linton NM, Kobayashi T, Yang Y, Hayashi K, Akhmetzhanov AR, Jung Sm, et al. Incubation period and other epidemiological characteristics of 2019 novel coronavirus infections with right truncation: a statistical analysis of publicly available case data. J Clin Med. (2020) 9:538. doi: 10.3390/jcm9020538

PubMed Abstract | CrossRef Full Text | Google Scholar

33. World Health Organization. Coronavirus Disease 2019 (COVID-19): Situation Report. WHO (2020).

Google Scholar

34. Wu Z, McGoogan JM. Characteristics of and important lessons from the coronavirus disease 2019 (COVID-19) outbreak in China: summary of a report of 72 314 cases from the Chinese Center for Disease Control and Prevention. JAMA. (2020) 323:1239–42. doi: 10.1001/jama.2020.2648

CrossRef Full Text | Google Scholar

35. Kissler S, Tedijanto C, Goldstein E, Grad Y, Lipsitch M. Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period. Science. (2020) 368:860–8. doi: 10.1101/2020.03.04.20031112

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zou L, Ruan F, Huang M, Liang L, Huang H, Hong Z, et al. SARS-CoV-2 viral load in upper respiratory specimens of infected patients. N Engl J Med. (2020) 382:1177–9. doi: 10.1056/NEJMc2001737

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Lakshmikantham V, Vatsala A. Theory of differential and integral inequalities with initial time difference and applications. In: Rassias TM, Srivastava HM, editors. Analytic and Geometric Inequalities and Applications. Dordrecht: Springer (1999). p. 191–203. doi: 10.1007/978-94-011-4577-0_12

CrossRef Full Text | Google Scholar

38. Hethcote HW. The mathematics of infectious diseases. SIAM Rev. (2000) 42:599–653. doi: 10.1137/S0036144500371907

CrossRef Full Text | Google Scholar

39. van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Diekmann O, Heesterbeek JAP, Metz JA. On the definition and the computation of the basic reproduction ratio R₀ in models for infectious diseases in heterogeneous populations. J Math Biol. (1990) 28:365–82. doi: 10.1007/BF00178324

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Anderson RM, May RM. Vaccination and herd immunity to infectious diseases. Nature. (1985) 318:323–9. doi: 10.1038/318323a0

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Anderson RM. The concept of herd immunity and the design of community-based immunization programmes. Vaccine. (1992) 10:928–35. doi: 10.1016/0264-410X(92)90327-G

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ritchie H, Ortiz-Ospina E, Beltekian D, Mathieu E, Hasell J, Macdonald B, et al. Coronavirus (COVID-19) Vaccinations. Statistics and Research, Our World in Data. (2021) (accessed January 24, 2021).

Google Scholar

44. Curley B. How Long Does Immunity from COVID-19 Vaccination Last? Healthline. (2021) Available online at: https://www.healthline.com/health-news/how-long-does-immunity-from-covid-19-vaccination-last (accessed June 11, 2021).

Google Scholar

45. Seow J, Graham C, Merrick B, Acors S, Pickering S, Steel KJ, et al. Longitudinal observation and decline of neutralizing antibody responses in the three months following SARS-CoV-2 infection in humans. Nat Microbiol. (2020) 5:1598–607. doi: 10.1038/s41564-020-00813-8

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Dan JM, Mateus J, Kato Y, Hastie KM, Yu ED, Faliti CE, et al. Immunological memory to SARS-CoV-2 assessed for up to 8 months after infection. Science. (2021) 371:6529. doi: 10.1126/science.abf4063

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Centers for Disease Control and Prevention. Interim Public Health Recommendations for Fully Vaccinated People. CDC Information. (2021). Available online at: https://www.cdc.gov/coronavirus/2019-ncov/vaccines/fully-vaccinated-guidance.html (accessed June 11, 2021).

Google Scholar

48. Tenforde MW. Effectiveness of Pfizer-BioNTech and Moderna Vaccines AgainstCOVID-19 Among Hospitalized Adults Aged 65 Years-United States, January-March 2021. Morbidity and mortality weekly report. (2021) 70.

PubMed Abstract | Google Scholar

49. Bernoulli D. Essai d'une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l'inoculation pour la prévenir. Hist de l'Acad R Sci. (1760) 1:1–45.

Google Scholar

50. Ross R. The Prevention of Malaria. John Murray (1911).

Google Scholar

51. Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proc R Soc Lond Ser A. (1927) 115:700–21. doi: 10.1098/rspa.1927.0118

CrossRef Full Text | Google Scholar

Appendix I: Entries of the Non-Negative Matrix F

$\begin{array}{l}{f}_1={\beta }_{P_1}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],f_2={\beta }_{I_1}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],\end{array}$

$\begin{array}{l}{f}_3={\beta }_{A_1}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],f_4={\beta }_{H_1}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],\end{array}$

$\begin{array}{l}{f}_5=(1-{ϵ}_0){\beta }_{P_2}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],\\ f_6=(1-{ϵ}_0){\beta }_{I_2}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],\\ f_7=(1-{ϵ}_0){\beta }_{A_2}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],\\ f_8=(1-{ϵ}_0){\beta }_{H_2}\left[\frac{S_{1u}^{*}+(1-{ϵ}_v)S_{1v}^{*}}{N^{*}}\right],\end{array}$

$\begin{array}{l}{g}_1=(1-{ϵ}_i){\beta }_{P_1}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right],\\ g_2=(1-{ϵ}_i){\beta }_{I_1}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right],\\ g_3=(1-{ϵ}_i){\beta }_{A_1}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right],\\ g_4=(1-{ϵ}_i){\beta }_{H_1}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right],\\ g_5=(1-{ϵ}_i)(1-{ϵ}_0){\beta }_{P_2}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right],\\ g_6=(1-{ϵ}_i)(1-{ϵ}_0){\beta }_{I_2}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right],\\ g_7=(1-{ϵ}_i)(1-{ϵ}_0){\beta }_{A_2}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right],\\ g_8=(1-{ϵ}_i)(1-{ϵ}_0){\beta }_{H_2}\left[\frac{S_{2u}^{*}+(1-{ϵ}_v)S_{2v}^{*}}{N^{*}}\right].\end{array}$