Modeling the SARS-CoV-2 Omicron variant dynamics in the United States with booster dose vaccination and waning immunity

1.   Introduction

Throughout the entire COVID-19 pandemic, there has been an opportunity for the emergence of several SARS-CoV-2 variants due to the constant high transmission. These variants alter the transmission features of the virus, the vaccine effectiveness, and the severity of the infection. Multiple variants have emerged since December 2020, but not all have affected the dynamics of the pandemic. Said variants are catalogued by the World Health Organization (WHO) as variants of concern.

The B.1.1.529 variant, also known as Omicron, was first reported in South Africa on November 24, 2021 ^[1], but it was circulating in the United Kingdom on November 27, 2021 ^[2], which means that the variant was present in the population long before it was sequenced. This variant has over 50 mutations, and only 30 of them are in the Spike protein, which developed a higher capacity of transmission and affinity for the ACE2 binding receptor protein ^[3]. In addition, it has the capacity to evade the protection of the immune system. Neutralization by vaccination antibodies decreased significantly in the vaccinated community, which may increase breakthrough infections ^[4]. Moreover, the Omicron variant can evade neutralization antibodies of individuals that recovered from the infection and were not previously vaccinated ^[5]. These individuals can be re-infected with the Omicron variant, regardless of the previous variant infection.

The Omicron variant has three distinct features with respect to the former SARS-CoV-2 dominant variant (Delta): we have already discussed two of them, which are a higher transmission and evasion of the immune system. The other feature is severity: the Omicron variant develops a less severe infection than the Delta variant. This behavior can be explained by the fact that a booster shot was already applied prior to the Omicron wave in some countries, and prior infections of the unvaccinated in the Delta wave may explain the reduced severity ^[6].

Since the COVID-19 pandemic started, several methods have been applied to model the dynamics of this disease. Although many works have used systems of ordinary differential equations to this end, some authors have employed different methods, such as stochastic epidemic models ^[7,8], including some with delay ^[9]. Partial differential equation models were used in ^[10,11,12] to study the geographical spread of COVID-19. Furthermore, fractional differential equations have been employed in recent years to model several biological and epidemiological problems ^[13,14]. Some fractional models for COVID-19 have been studied in ^[15,16] using the Caputo–Fabrizio derivative and in ^[17] using the Atangana–Baleanu derivative. Gatto et al. ^[18] estimated the parameters of a spatially explicit SEIR model using mobility and epidemiological data from Italy. Bertozzi et al. ^[19] used three different parsimonious models: an exponential growth model, a self-exciting branching process and a compartmental SIR model to forecast the course of the pandemic in several countries. Other works, such as ^[20], have aimed to model the dynamics of SARS-CoV-2 inside the human body in the presence of humoral immunity and antiviral treatment.

Recently, some mathematical models have been proposed to explain the dynamics of the Omicron wave. Yu et al. ^[21] used an SEIHRD model with 11 equations to estimate the transmissibility of the Beta, Delta and Omicron variants using data from multiple COVID-19 waves in South Africa. Other studies, which included estimations for immune escape, were performed by Yang and Shaman ^[22] using a model-inference system, while Gozzi et al. ^[23] used a stochastic multi-strain epidemic model. Gowrisankar et al. ^[24] used the fractal interpolation function to predict the seven-day moving average of daily positive Omicron cases in six different countries, while Grabowski et al. ^[25] used a two-strain model, which allowed them to infer that the growth advantage of Omicron with respect to Delta may be attributed to immune evasion. Özköse et al. ^[26] employed fractional order modeling to examine the spread of the Omicron variant and its relationship with heart attacks using data from the United Kingdom. However, to our knowledge, there are no current mathematical models that describe the dynamics of the Omicron wave in the United States and provide short-term predictions of that wave.

In this article, we apply mathematical modeling to simulate the potential impact of the Omicron variant on the US population. We include parameters representing the use of face masks and their efficiency, the dynamics of the unvaccinated population and the population vaccinated with booster doses from three different manufacturers, as well as the waning of natural and acquired immunity.

The rest of this paper is structured as follows. We formulate the mathematical model in Section 2. In Section 3, we compute the basic reproduction number and develop a stability analysis for the disease-free equilibrium. In Section 4, we calibrate our mathematical model using daily cumulative data of infected, hospitalized and dead individuals, and the vaccination rates per vaccine manufacturer. In Section 5, we explore the simulations of the cases in the US. In Section 6, we perform a local and global sensitivity analysis for our model. Lastly, we provide a discussion and concluding remarks in Section 7.

2.   Formulation of the model

In this work, we will develop and apply a compartmental differential equation model to understand the dynamics of the spread of the Omicron variant in vaccinated and unvaccinated populations. We consider the importance of waning immunity for natural infection and acquired immunity. The mathematical model contains two separate sets of equations for the virus strain that affects susceptible and vaccinated subpopulations. Our $SV_1V_2V_3(EIAHW)^2RD$ model evaluates the dynamics of fifteen subpopulations at any given time $t$, which are denoted as $S(t)$, $E_1(t)$, etc. The biological interpretation of the model equations is described below.

Susceptible population $S(t)$: We consider natural recruitment at a rate $\Lambda$ and a natural death rate $d$. The susceptible population decreases at a rate $\beta_1$ by contact with people with symptomatic Omicron variant infection, and a rate $\beta_2$ by contact with people with asymptomatic Omicron variant infection. This subpopulation may increase at a rate $2K_v$ due to the waning immunity from the protection provided by the vaccine and a rate $2K_n$ due to the waning immunity from natural infection. Finally, this subpopulation decreases because they can be vaccinated by three types of vaccines at a rate $\rho_i\ge0$ ($i = 1$ denotes vaccination with Moderna, $i = 2$ represents Pfizer, and $i = 3$ is the Janssen (Johnson & Johnson) vaccine). The vaccination rates are multiplied by the all-or-nothing protection for the specific vaccine ^[27], which is provided at a rate $1-\varepsilon_{a, i}$ ($i = 1, 2, 3$). Hence,

where $\lambda = \beta_1\left(I_1 + \kappa I_B\right) + \beta_2\left(A_1 + \kappa A_B\right)$. The parameter $\kappa\in[0, 1]$ quantifies the reduction in infectivity for vaccinated individuals.

Unvaccinated population exposed to the Omicron variant $E_1(t)$: Since we are modelling the infection caused by the virus, we need to include the latent stage, i.e., the exposed population. It increases at a rate $\lambda$ (defined as above) due to symptomatic and asymptomatic infections from the Omicron variant. The rate of increment may be slowed down by the efficiency of face masks $\epsilon_m$ and usage of face masks by the population $C_m$. The exposed subpopulation represents newly infected individuals that have not received any type of vaccine. This subpopulation decreases at a rate $w$ as the exposed individuals become infectious, where $1/w$ denotes the average length of the latent period. Hence,

Symptomatic, unvaccinated population infected by the Omicron variant $I_1(t)$: Infected individuals will generate symptoms for the Omicron variant at a rate $p_1\in[0, 1]$ from the exposed unvaccinated subpopulation. They recover at a rate $\gamma_1$ and are hospitalized at a rate $\alpha_1$. Consequently,

Asymptomatic, unvaccinated population infected by the Omicron variant $A_1(t)$: This subpopulation is infected by the virus, but they do not develop any type of symptoms of COVID-19. It is important to include them in the model, because they are undetected and spreading the virus without knowing that they are infected. They are generated from the exposed compartment at a proportion $1-p_1$ and recover at a rate $\gamma_1$. Therefore,

Hospitalized, unvaccinated population infected by the Omicron variant $H_1(t)$: This population increases at a rate $\alpha_1$ from the symptomatic infections. The individuals in this subpopulation might die due to COVID-19 at a rate $\omega\chi$. The surviving individuals recover at a rate $\omega\left(1-\chi\right)$. The parameter $1/\omega$ represents the average hospitalization period, while $\chi\in[0, 1]$ denotes the probability of death by COVID-19. Hence,

Subpopulation with partial immunity (one dose) $V_1(t)$: This subpopulation increases at a rate $\rho_i\ge0$, multiplied by the all-or-nothing protection for the specific vaccine provided at a rate $1-\varepsilon_{a, i}$. It decreases at a rate $\left(1-\varepsilon_{L, i}\right)\lambda\left(1-C_m\epsilon_m\right)$ due to the imperfect vaccine effectiveness provided by one dose of each of the vaccines applied in the US, since people with partial immunity may still get infected by contact with symptomatic or asymptomatic individuals. This compartment includes also people with one dose of the J & J vaccine, which are considered fully vaccinated; it decreases at a rate $2K_{v1}$ due to the waning immunity of the vaccine. Finally, this population also decreases when individuals receive the second dose of the Moderna or Pfizer vaccines, which occurs at a rate $\theta$. Consequently,

Subpopulation with full immunity (two doses) $V_2(t)$: This population increases at a rate $\theta$ due to the application of second doses to the individuals that are part of the partial immunity subpopulation. Due to the imperfect vaccine, this population decays at a rate $\left(1-\varepsilon_{L2, i}\right)\lambda\left(1-C_m\epsilon_m\right)$ due to infections caused by the interaction with symptomatic and asymptomatic individuals infected with the Omicron variant. Finally, the population decays as well at a rate $2K_{v2}$ due to the loss of protection or waning immunity provided by the protection of the Pfizer and Moderna vaccines.

Vaccinated population exposed to the Omicron variant $E_B(t)$: This subpopulation is analogous to $E_1(t)$ but for vaccinated population, i.e., breakthrough infections. Vaccinated individuals who interact with symptomatic or asymptomatic infected individuals are incorporated at rates $1-\varepsilon_{L, i}$ ($i = 1, 2, 3$) and $1-\varepsilon_{L2, i}$ ($i = 1, 2$), considering the vaccine effectiveness of the vaccines to the Omicron variant for one and two doses, respectively. It decreases at a rate $w$ as exposed individuals become infectious.

Symptomatic, vaccinated population infected by the Omicron variant $I_B(t)$: For all individuals that are infected while being vaccinated, regardless of the number of doses received, we assume they develop symptoms at a proportion $p_2\in[0, 1]$ from the exposed vaccinated subpopulation. They recover at a rate $\gamma_2$ and are hospitalized at a rate $\alpha_B$. Consequently,

Asymptomatic, vaccinated population infected by the Omicron variant $A_B(t)$: We assume that individuals that are infected while being vaccinated have a probability $1-p_2$ of not developing symptoms, regardless of the number of doses received. They also recover at a rate $\gamma_2$. Therefore,

Hospitalized, vaccinated population infected by the Omicron variant $H_B(t)$: This population increases at a rate $\alpha_B$ from the vaccinated symptomatic infections. This subpopulation includes death from COVID-19 at a rate $\omega_B\chi_B$. The individuals that do not die recover at a rate $\omega_B\left(1-\chi_B\right)$, where $\chi_B\in[0, 1]$ denotes the probability of death due to COVID-19 for the vaccinated population. Hence,

Recovered population $R(t)$: All non-hospitalized individuals (with symptoms or not) will recover at a rate $\gamma_1$ for the unvaccinated population and $\gamma_2$ for the vaccinated population. Hospitalized unvaccinated individuals enter this subpopulation at a rate $\omega\left(1-\chi\right)$, and hospitalized vaccinated individuals recover at a rate $\omega_B\left(1-\chi_B\right)$. This population decreases at a rate $2K_n$ due to waning of natural immunity after infection with COVID-19. Thus,

Subpopulation with waning natural immunity $W(t)$: This subpopulation increases as the immunity provided by the natural infection of the recovered population wanes at a rate $2K_n$. This increment is associated with a decrement of the same rate, since these individuals become susceptible of getting infected again when they lose their protection completely. Parameters are multiplied by two, because we are considering the waning process in two stages ^[28], where the mean duration of natural immunity is $1/K_n$. Therefore,

Waning immunity of the vaccinated subpopulation $W_V(t)$: This population increases at two rates: the loss of protection from the J & J vaccine at a rate $2K_{v1}$ and the loss of protection provided by the two-dose vaccines at a rate $2K_{v2}$. Here, parameters are multiplied by two, since the waning of immunity due to vaccination also occurs in two stages. Vaccinated individuals can receive a booster dose (third dose for Moderna and Pfizer, second dose for J & J) to increase their protection at a rate $\rho_4$, which reduces the population with waning vaccine immunity. Vaccinated individuals who do not receive a booster shot lose their protection and may become susceptible to infection at a rate $2K_v$. Consequently,

Subpopulation with full immunity with a booster dose $V_3(t)$: This population increases at a rate $\rho_4$ when eligible individuals receive a third dose (for Moderna and Pfizer) or second dose (for J & J) to enhance their protection against COVID-19. Some individuals with booster dose can get infected at a rate $\lambda\left(1-C_m\epsilon_m\right)\mathcal{E}_3$ because of the vaccine leakiness. Individuals who receive their dose and do not get infected because of the imperfection of vaccine stay at this compartment and are fully protected from the Omicron variant. Hence,

Therefore, the model that will be studied here is defined by the system

where

Furthermore, we define an additional variable $D(t)$ that denotes the number of people deceased due to COVID-19, which is described by the equation

We assume that $\rho_1, \rho_2, \rho_3, \rho_4, K_{v1}, K_{v2}, K_n\ge0$ and $\kappa, p_1, p_2, \varepsilon_{a, i}, \varepsilon_{L2, i}, \mathcal{E}_3\in[0, 1]$ for $i = 1, 2, 3$, while all other parameters are strictly positive.

It is readily verified that every solution of (2.1) with nonnegative initial conditions remains nonnegative for $t\ge0$. If we define $N = S_1+E_1+I_1+A_1+H_1+V_1+V_2+E_B+I_B+A_B+H_B+R+W+W_V+V_3$, we obtain from (2.1) that

This implies that the set

is a positively invariant and attractive set for system (2.1).

3.   Theoretical analysis

In this section, we will prove some mathematical results for system (2.1). We will also perform a stability analysis of the diasease-free equilibrium in the case when the vaccination rates are zero and the booster dose is positive.

3.1.   The disease-free equilibrium

We will proceed to find the disease-free equilibrium (DFE) of (2.1). Now, we set the right-hand side of the system equal to zero and assume $I_1 = I_B = 0$.

If $I_1 = 0$, it follows from the third equation of the system that $p_1wE_1 = 0$, so $E_1 = 0$. In consequence, by the fourth equation, we have $-(\gamma_1+d)A_1 = 0$, so $A_1 = 0$.

A similar argument, using $I_B = 0$ with the ninth and tenth equations tells us that $E_B = A_B = 0$. Furthermore, from the fifth and eleventh equations, we obtain $H_1 = 0$ and $H_B = 0$. Substituting in the twelfth and thirteenth equations, we have $-2K_nR - dR = 0$ and $2K_nR - 2K_nW - dW = 0$. From this, it follows that $R = 0$ and $W = 0$. In consequence, the system is reduced to

Solving the above system, we obtain that the DFE is given by

with

The DFE of system (2.1) is non-negative if and only if

In the particular case when $\rho_1 = \rho_2 = \rho_3 = 0$, the DFE becomes

3.2.   Basic reproduction number

We will use the next-generation matrix method to determine the basic reproduction number of system (2.1), considering that the infected compartments of the population are $E_1$, $I_1$, $A_1$, $H_1$, $E_B$, $I_B$, $A_B$ and $H_B$.

Let $m_1 = w+d$, $m_2 = \alpha_1+\gamma_1+d$, $m_3 = \gamma_1+d$, $m_4 = \omega+d$, $m_5 = w+d$, $m_6 = \alpha_B+\gamma_2+d$, $m_7 = \gamma_2+d$ and $m_8 = \omega_B+d$.

Using the notation in ^[29], we obtain

Then, we have

and

where

Calculating the inverse of $V$, and then multiplying $F$ and $V^{-1}$, we get

The basic reproduction number with vaccination, also known as the control reproduction number, is given by $\mathcal{R}_c = \varrho\left(FV^{-1}\right)$, where $\varrho$ denotes the spectral radius. Hence,

The basic reproduction number with no vaccination, denoted by $\mathcal{R}_0$, is obtained by making $\rho_1 = \rho_2 = \rho_3 = 0$ in (3.1), which yields

3.3.   Analysis when the vaccination rates are zero

In this subsection, we will establish some local and global stability results for the DFE of model (2.1) in the case when the vaccination rates $\rho_1$, $\rho_2$ and $\rho_3$ are zero. However, we still allow the booster dose vaccination rate $\rho_4$ to be positive.

Theorem 1. The disease-free equilibrium $E_0$ of system (2.1) with $\rho_1 = \rho_2 = \rho_3 = 0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$, and it is unstable if $\mathcal{R}_0 > 1$.

Proof. This follows directly from the calculation of the basic reproduction number and an application of [29,Theorem 2].

Theorem 1 establishes a sufficient condition for the local stability of $E_0$. This implies that, if $\mathcal{R}_0 < 1$, the epidemic will become extinct provided that the initial conditions are sufficiently close to the disease-free equilibrium. However, we would like to obtain a condition that guarantees the eradication of the disease independently of the initial state. For this, we need to prove the following theorem, which gives a sufficient condition for the global stability of the DFE.

Theorem 2. Suppose that

Then, the disease-free equilibrium $E_0$ of system (2.1) with $\rho_1 = \rho_2 = \rho_3 = 0$ is globally asymptotically stable in $\Omega$.

Proof. Consider the following Lyapunov function:

where

Then, the time derivative of $\mathcal{L}$ evaluated at the solutions of system (2.1) is given by

After simplification, we obtain

Since the set $\Omega$ is positively invariant, we have $S_1 + V_3 \le \frac{\Lambda}{d}$ for all $t\ge0$. Then, the hypothesis (3.3) ensures that $\mathcal{L}'\le0$. Moreover, $\mathcal{L}' = 0$ if and only if

Substituting (3.4) in system (2.1) with $\rho_1 = \rho_2 = \rho_3 = 0$ shows that the solutions tend to the DFE as $t\to\infty$. Hence, the largest positively invariant set where $\mathcal{L}' = 0$ is $\{E_0\}$. Therefore, by LaSalle's invariance principle, we conclude that the DFE is globally asymptotically stable.

Theorem 3. System (2.1) with $\rho_1 = \rho_2 = \rho_3 = 0$ has a unique endemic equilibrium of the form

if and only if

where

Proof. To find the endemic equilibria of the model, we make the right-hand side of (2.1) equal to zero, which yields

From (3.11) and the positivity of parameters, it follows that $V_1 = 0$. Substituting successively in (3.12), (3.19), (3.20) and (3.13), we obtain $V_2 = 0$, $W_V = 0$, $V_3 = 0$ and $E_B = 0$. Similarly, Eqs (3.14)–(3.16) yield $I_B = 0$, $A_B = 0$ and $H_B = 0$.

From (3.8), we have

Substituting in (3.9), we get

From (3.10),

Substituting in (3.17), we obtain

From (3.18),

From (3.6),

Lastly, from (3.7), we get

Assuming $I_1 \ne 0$ and $2K_nG_1 - G_3 \ne 0$, it follows that

Therefore, we can see that condition (3.5) is necessary and sufficient for all variables to be nonnegative and $I_1$ to be positive. Hence, the theorem follows.

4.   Estimation of parameters

This section is devoted to the estimation of baseline values for the model parameters. The values obtained here will be used in Section 5 to make predictions for the dynamics of the Omicron wave in the United States, and in Section 6 to perform the sensitivity analysis of the model.

To estimate the parameter values related to the unvaccinated population, we will first fit the model solutions to real data during the period before vaccination started in the US. Then, we will obtain the parameter values related to the vaccinated population, using the data from December 13, 2020 to January 4, 2022, which includes the Alpha and Delta waves, and the onset of the Omicron wave.

4.1.   Model without vaccination

We used cumulative data for COVID-19 cases and deaths in the United States provided by the Johns Hopkins University repository ^[30], as well as daily hospitalization data obtained from ^[31], to fit the parameters of the model without vaccination. We considered the period from October 18, 2020 to January 14, 2021, which corresponds to the rising phase of the third wave in the US.

In the absence of vaccination, the model can be described by the following system:

If we define $\lambda_{10} = \beta_1\left(1-C_m\epsilon_m\right)$ and $\lambda_{20} = \beta_2\left(1-C_m\epsilon_m\right)$, then we can write $\lambda\left(1-C_m\epsilon_m\right) = \left(\beta_1I_1 + \beta_2A_1\right)\left(1-C_m\epsilon_m\right) = \lambda_{10}I_1 + \lambda_{20}A_1$.

We take as fixed the parameter values shown in Table 1, while the values for $\lambda_{10}$, $\lambda_{20}$, $\alpha_1$, $\chi$ and $\omega$ will be estimated by minimizing the sum of squared errors (SSE).

For model (4.1), we chose as initial conditions

The values for $H_1(0)$ and $D(0)$ were chosen as the reported number of hospitalized people and accumulated deaths due to COVID-19, respectively, for October 18, 2020. The number of active symptomatic cases $I_1(0)$ was taken as the number of infected cases reported in the previous 14 days (COVID-19 cases are usually not considered active if more than two weeks have passed since infection) minus the number of hospitalized cases on that date. Based on previous works with similar models ^[32], we assume that the number of asymptomatic cases $A_1(0)$ is eight times $I_1(0)$. The exposed population $E_1(0)$ was assumed equal to $I_1(0) + A_1(0)$. Since the starting date of simulations is about 9 months later than the date when COVID cases where first reported in the US, we assume that recovered population still has immunity against re-infection (see ^[33]); hence, we set $W(0)$ as zero. Finally, we computed the size of susceptible population $S_1(0)$ by subtracting the other initial values from the total population of the US (332.9 million people).

With the above assumptions, we obtained the best fit values for $\lambda_{10}$, $\lambda_{20}$, $\alpha_1$, $\chi$ and $\omega$ as shown in Table 2. Figure 1 shows a comparison between the solutions obtained with the model and reported cumulative number of symptomatic/asymptomatic infected cases, hospitalized people, and cumulative number of COVID-19 deaths.

4.2.   Estimation during vaccination period

In this section, we estimate the parameters of the full model with vaccination during the period from December 13, 2020 to January 4, 2022.

To model the vaccination rates, we will assume that the number of people per day who receive their first dose is increasing as the government accelerates the deployment of vaccines. This rate will reach a maximum and then start to decrease. As time tends to infinity, the vaccination rate for first doses will tend to zero since the vaccination program will focus instead on booster doses.

Hence, the time-dependent vaccination rates $\rho_i$ ($i = 1, 2, 3$) can be described by the following functions:

Here, $t$ represents the time since the date when vaccination began for Moderna and Pfizer vaccines (December 14, 2020). Vaccination with the Johnson & Johnson vaccine started 81 days later (on March 5, 2021), thus, the vaccination rate $\rho_3$ for Johnson & Johnson will be offset by a parameter $t_J = 81$ days.

We estimated the values of the parameters $\pi_i$, $\sigma_i$ ($i = 1, 2, 3$) using the reported number of vaccinations ^[37] for the period June 1, 2021 to January 13, 2022. For this, we denote the total number of people who have received at least one dose of Moderna, Pfizer and J & J vaccines by $T_1(t)$, $T_2(t)$ and $T_3(t)$, respectively, such that

The best fit values for each parameter, obtained by minimizing the SSE, are given in Table 3. The simulations for $T_1(t)$, $T_2(t)$ and $T_3(t)$ using these parameter values can be seen in Figure 2, and the graphs of the functions $\rho_1(t)$, $\rho_2(t)$ and $\rho_3(t)$ are shown in Figure 3.

To estimate the transmission rates $\beta_1$ and $\beta_2$, we use the values for $\lambda_{10}$ and $\lambda_{20}$ obtained in Section 4.1. This yields an estimate of $\beta_1 = 4.1979\times10^{-10}$ and $\beta_2 = 4.5929\times10^{-10}$ for the original variant of SARS-CoV-2. To account for the increased transmissibility of the variants that appeared subsequently in the US, we assume that the values for $\beta_1$ and $\beta_2$ become 2.8 times larger from April 18, 2021 to December 12, 2021 (the approximate time when the Alpha and Delta variants became dominant in the US), and five times larger from December 13, 2021 onwards (when the Omicron variant became dominant).

The vaccination rate for booster doses $\rho_4$ is set as

where $t_B = 242$ days corresponds to the time elapsed between the application of the first vaccines and the date when booster doses began to be applied (August 13, 2021).

The rest of the parameter values for system (2.1) were chosen as in Table 4. We consider the initial condition

where $H_1(0)$ and $D(0)$ are taken from the reported data for December 13, 2020, and the other values are taken from the simulations in Section 4.1.

Figure 4 depicts a comparison between the model solutions and the reported data from December 13, 2020 to January 4, 2022. The rest of the simulations for this period is shown in Figure 5.

5.   Results

5.1.   Simulation of the model without modifications

In this section, we simulated our model to evaluate the behavior of the usage of face masks, vaccination and the importance of the booster dose in the community of the US. We considered three types of vaccine that are currently being applied in the US: Pfizer, Moderna and Johnson & Johnson, mixed with the usage of improved cloth face masks. The simulations were carried out based on the estimated and fixed baseline parameters presented in Tables 1–4.

Cases of infected individuals from COVID-19 are starting to descend in mid-January, but the descent is slow in comparison with the increase of cases that started in December 2021. Between February and March, cases will descend until reaching a valley in the first week of April (Figure 6A, B). Between May and June, cases will increase, but the peak will be 2.88 times smaller than the peak from January 2022 (Figure 6A, B). Unvaccinated individuals are the carriers of this wave of Omicron cases: there will be approximately one breakthrough case (an infected vaccinated individual) for every 15 cases of unvaccinated individuals. The trend is maintained regardless of the presence or absence of symptoms (Figure 6A, B). In terms of hospitalizations, the peak was in the first week of February, reaching a peak of 140,000 hospitalized individuals. The unvaccinated population is the one that injects more hospitalized individuals compared to the vaccinated population (Figure 6C). Furthermore, the unvaccinated community has a higher hospitalization rate once they are exposed to the virus in comparison with the vaccinated community (Figure S1). In the following months, the number of people who are still protected by vaccination will descend, due to the waning immunity. The peak of booster doses was reached in the last two weeks of December 2021. The slight increase in breakthrough cases can be related to not getting a booster dose. The trend maintains for individuals getting their first dose or their second dose (Figure 6D). By simulating the behavior of vaccination based on the number of unvaccinated individuals, as well as people vaccinated with two doses and with a booster dose, we can see how the behavior of COVID-19 will be in the following months (Figure 6D).

The recovery of individuals from COVID-19 will increase in the following months regardless of the vaccination status (Figure S2A). Waning natural immunity will increase in the following months, and it will stabilize between April and May 2022 (red line in Figure S2B), the waning of acquired immunity (vaccinated individuals) is stabilized in the following months (blue dashed in line Figure S2B). The death toll of COVID-19 will increase in the following months, reaching 900,000 in early February 2022 (Figure S2C), and it will stabilize with less than one million deaths between March and April 2022. This stabilization is related to the decrease in new COVID-19 infections and hospitalizations (Figure 6A–C). By May 2022, deaths can increase based on another wave of COVID-19 that will start to increase between April and May 2022 (Figure 6A, B). Overall, 90% of the deaths of COVID-19 are in unvaccinated individuals, and this trend will continue in the following months (Figure S2D). Death of vaccinated individuals will continue to be low due to the protection provided by the vaccine (Figure S2D, blue dashed line).

5.2.   Assessing the importance of the usage of face masks

In this section, we evaluate the behavior of two important parameters: the proportion of the population that uses a face mask and the type of face mask being used by said population. We estimated the implementation of the use of face masks in the first week of January in the US. We evaluated four types of face masks: cloth masks, improved cloth masks, surgical face masks and N95 face masks. Cloth face masks have a 30% of protection against getting infected with the Omicron variant. In Figure S3, we evaluated five scenarios of the usage of cloth face masks based on the percentage of the population: from 0% to 100% of the population using face masks. For asymptomatic and symptomatic individuals, the usage of face masks for the reduction in new cases is important when cases are still increasing after reaching the peak. If 75% of individuals used a cloth mask the peak would have been slightly higher than one million compared with almost 1.25 million cases when individuals do not wear face masks, a reduction of 250,000 cases of COVID-19 (Figure S3A, B). The same trend is observed in hospitalized individuals (Figure S3C). 75% of the population using cloth face masks will be associated with roughly 20,000 fewer deaths compared with no usage of face masks (Figure S3D).

An improved cloth face mask gives you a 50% chance of not getting infected by COVID-19. As we increase the efficiency of face masks, the difference between in the number of new infections is clearly defined. The difference between 75% of the population using face masks and only 25% is roughly 400,000 symptomatic infections for the unvaccinated population (Figure 7A) and approximately 50,000 in the population of vaccinated individuals (Figure S4A) when the infections are increasing. Once the cases are decreasing, there is no difference between the usage of face masks, and the cases will be stabilized between May and June 2022. The same pattern of difference between the use of face masks occurs in asymptomatic individuals (Figures 7B and S4B). If only 25% of the population had used an improved cloth face mask, this behavior would have doubled hospital admissions compared to 75% of the population in the community of vaccinated individuals (Figure S4C). For the unvaccinated community where 75% used improved cloth face masks, hospitalization admission would have decreased in 40,000 individuals compared with only 25% using a face mask (Figure 7C). In the unvaccinated community, the use of face masks could have prevented around 50,000 deaths (Figure 7D). For the vaccinated population, roughly 10 deaths could have been prevented with respect to the case when only 25% of the population wore a face mask (Figure S4D).

For the usage of a surgical face mask and the N95 face mask, the pattern is completely different compared with the cloth and the improved cloth face mask, but it is very similar between them (Figures S5 and S7). For the use of a surgical face mask by 75% of the population, the cases in the community of unvaccinated individuals would have started to descend by January compared with the case with 25% usage, where the peak would have been reached at the end of January (Figure S5A). In the period of the valley of cases, the use of face masks only makes a difference in the minimum of infected cases and the height of the peak it can reach when cases increase between June and July 2022 (Figure S5A). The same pattern is observed in asymptomatic individuals (Figure S5B). For the vaccinated population, the same pattern is observed compared with the unvaccinated, the only difference is the number of cases is reduced due to the additional protection provided by the vaccine in symptomatic and asymptomatic individuals (Figure S6A, B). Hospital admission would have been highly reduced comparing the usage of face masks to a reduction of almost 4000 admissions between 75% and 25% of the use of a surgical face mask in the vaccinated individuals (Figure S5C), and a reduction of 1000 hospital admissions in the vaccinated community (Figure S6C). The overall death toll is drastically reduced in the vaccinated population by approximately 150,000 deaths compared with 75% usage and 25% usage of a surgical face mask (Figure S5D) and a reduction of 20 deaths in the vaccinated population (Figure S6D).

If only 25% of the population of unvaccinated individuals would have used an N95 face mask, 50,000 new infected cases could have been prevented. With 75% of the population using an N95, the infections would have decreased drastically compared with 25% of the population regardless of the symptomology (Figure S7A, B). The same pattern of prevention is seen in new infected breakthrough cases for the population of vaccinated individuals, but the difference between this population and the unvaccinated is the proportion of new cases (Figure S7A, B). Hospitalization as well would have decreased in 80,000 admissions if 75% of the unvaccinated population would have worn N95 face masks (Figure S7C), the same behavior at a different scale is seen for the vaccinated population (Figure S8C). Overall, 170,000 deaths could have been prevented in January if 75% of the unvaccinated community wore N95 face masks (Figure S7D), while 25 deaths could have been prevented in the vaccinated population (Figure S8D).

5.3.   Evaluating the impact of vaccination rate in a population wearing face masks

Our model incorporates the vaccination rate (VR) as a variable that changes with respect to time, thus we can modify this rate by increasing or decreasing it with respect to the baseline vaccination rate. Unfortunately, the VR for getting the first dose is decreasing, so the application of the second dose is decreasing as well (Figure S9A, B; blue dotted line). The application of the booster dose was increasing in the month of December 2021, but by January 2022 it is decreasing (Figure S9C; blue dotted line). In this section, we increased the VR to 200% (black dotted line) and decreased it to 50% (red dotted line) with respect to the baseline vaccination to evaluate if vaccination by itself is important in controlling the behavior of the Omicron variant (Figure S9A–C).

By increasing the VR, the unvaccinated compartment had a lower amount of new infected individuals regardless of the symptomology (Figure S10A, B). Hospitalizations can decrease by 1200 with respect to the baseline vaccination rate and by 1500 with respect to 50% VR (Figure S10C). The death toll is significantly lower for an increased vaccination rate (Figure S10D). However, the population of unvaccinated individuals is not significantly decreased, despite increasing the VR, because it is a large community.

The vaccinated population is substantially affected by increasing or decreasing the VR (Figure 8). Baseline VR will result in 320,000 new breakthrough cases, while a 200% VR will decrease the number of new infected symptomatic infections by 20,000. Decreasing the VR will develop an increase of roughly 45,000 new infected cases (Figure 8A). The asymptomatic infections behave in the same manner (Figure 8B). Hospitalizations can be reduced by 125 individuals if we increase the VR with respect to baseline VR. A decrease in the VR will be associated with 100 new hospitalizations with respect to the baseline VR (Figure 8C). Deceased individuals are affected as well by accelerating the VR by decreasing the overall deaths (Figure 8D).

On the other hand, varying the overall vaccine effectiveness of the three vaccines being applied in the US makes a difference only in the peak of symptomatic cases, asymptomatic cases and hospitalizations; when cases are decreasing there is no difference between vaccine effectiveness (Figure S11).

5.4.   Varying the vaccine-acquired and natural immunity waning rates in a population wearing face masks

It is known that the natural immunity derived from COVID-19 recoveries and the immunity derived by vaccines wane over time. In this section, we assess the natural and vaccine-acquired immunity. If individuals recovering from COVID-19 develop lifelong protection by March 2022, the pandemic would cease to exist if no other variant emerges in that period (Figure 9; magenta dotted line). We vary the average natural immunity period between seven months (red dotted line) and eleven months (black dotted line) with a baseline protection of nine months (blue dotted line). The simulations were carried out assuming that 50% of the population uses a 50% efficient face mask (improved cloth face mask). By March 2022, cases will decline until reaching a valley, the depth of the valley depends on the natural immunity. If the protection of natural immunity remains for seven months, it will trigger a subsequent wave of roughly less than 900,000 infected individuals between May and June 2022, reaching a peak in July 2022 and decreasing to a plateau of 500,000 infected individuals (Figure 9A). Asymptomatic infections behave in the same manner (Figure 9B). Hospitalizations will peak in early July 2022 with less than 90,000 individuals (Figure 9C) and deaths can reach one million by the end of March 2022 (Figure 9D). For waning immunity of nine months, the peak would occur in mid-July 2022 with approximately 700,000 symptomatic infections and 1,500,000 asymptomatic infections (Figure 9A, B; blue dotted line). Hospitalizations will peak until July 2022 with only 60,000 individuals, plateauing in 50,000 admissions by November 2022. For the scenario where immunity lasts up to nine months, the peak will be reached by mid-August 2022 with only 600,000 infected symptomatic and 9,000,000 asymptomatic individuals (Figure 9A, B, black dotted line). Hospitalizations will peak until August 2022 with 50,000 new admissions and plateauing to a value of 40,000 admissions (Figure 9C). Deaths can reach up to 950,000 by March 2022 (Figure 9D). The length and height of the subsequent wave, when natural immunity wanes, depends on whether another variant emerges, a scenario that has happened during this pandemic. Booster doses provide enough protection to slowly decrease breakthrough cases by mid-April 2022. This depends on the possible emergence of another variant and how it affects vaccine effectiveness (Figure S12).

5.5.   Evaluating the combined effect of face mask usage, vaccination and other non-pharmaceutical strategies (NPIs)

We incorporated an additional simulation to our model to assess whether there is any difference between using or not different NPIs in addition to the use of face masks. To evaluate this behavior, we multiplied the effective transmission rates $\beta_1$ and $\beta_2$ by the factor $1-C_m$, where the parameter $C_m$ ranges from 0 to 1 and represents an additional reduction in the transmission rate ^[44]. For this simulation, we added the variation of face mask efficiency, and the usage of said masks was fixed to 50% of the population.

The dynamics of using cloth face masks mixed with NPIs are depicted in Figure S13. For the unvaccinated community, the implementation of NPIs (black dotted line) yields a reduction of 100,000 symptomatic infections during the peak of the Omicron wave with respect to not implementing NPIs (red dotted line, Figure S13A). For the asymptomatic infections, the reduction is much clearer: the infections are curtailed by roughly 3,000,000 cases when compared with the absence of NPIs (Figure S13B). Hospital admissions are reduced by nearly 10,000 between mid-January 2022 (Figure S13C), and deaths can be reduced by 500 when NPIs are applied (Figure S13D). The vaccinated population is affected by the implementation of NPIs: symptomatic infections can be curtailed by 40,000 (Figure S14A), and asymptomatic infections are reduced by 200,000 in mid-January (Figure S14B). Hospitalizations are reduced by 400 compared with using NPIs (Figure S14C), while the death toll is reduced by 50 in mid-January and February. Due to the implementation of NPIs, the waning of natural infection is heavily influenced: by May 2021, there will be roughly 112 million individuals compared with 120 millions, and there will be 8,000,000 individuals whose protection has waned and are susceptible to infections by May 2021 (Figure S15A). Meanwhile, the number of vaccine-acquired infections is not influenced and decreases in the following months (Figure S15B).

For an improved face mask, which has an effectiveness of 50%, the peak comes forward for the unvaccinated community with NPIs, and it is reduced to nearly 18,000 symptomatic infections (Figure 10A). For asymptomatic infections, the peak occurs early with a reduction of roughly 2,000,000 infections by early February 2022 (Figure 10B). Hospitalizations are curtailed as well with a reduction of 2,000 admissions by the end of January 2022 (Figure 10C). The death toll is reduced by 450 deaths compared with not using NPIs (Figure 10D).

For the vaccinated population, the peak comes forward when we consider NPIs, and it has a reduction of nearly 30,000 infected asymptomatic individuals and 160,000 asymptomatic infections (Figure S16A, B). For hospital admissions, there was a reduction of 100 compared with no NPIs, and a global reduction of 15 deaths (Figure S16C, D). Waning immunity of the natural infections is influenced by implementation of NPIs: by May 2022, there will be 6,000,000 more individuals whose immunity (Figure 11A) has waned and have become susceptible; for the acquired immunity, using NPIs only slightly accelerates the decay and fewer individuals have waning immunity (Figure 11B).

Regarding the use of surgical masks and N95, the implementation of NPIs only accelerates the decline of cases. For the unvaccinated population, symptomatic cases will descend until reaching a value with 30,000 fewer infected compared with not using NPIs by Abril 2022 (Figure S17A). The same pattern of deceleration occurs for the asymptomatic infections, reaching a valley in April 2022 with 2,000,000 fewer cases than in the case with no NPIs (Figure S17B). Hospitalizations follow the same pattern: there is a reduction in April 2022 with a decay of 6,000 admissions, and a reduction of 250 deaths if NPIs are implemented by early January 2022 (Figure S17C, D). The same pattern of decay is observed in the vaccinated community: NPIs will curtail symptomatic infections (Figure 18A), as well as asymptomatic infections (Figure 18B). Global hospitalization admissions and deaths are reduced when NPIs are mixed with the use of a surgical mask (Figure S18C, D).

Waning immunity of natural infections would stabilize by April 2022 with roughly 25,000,000 fewer individuals if NPIs were implemented at the beginning of this year; these individuals are susceptible to re-infection if transmission is still high in the US. By July 2022, waning immunity will increase again for individuals that were not infected in the wave between June and August 2022 (Figure S19A). For the vaccinated population, individuals with waning immunity will descend, meaning that they will not become susceptible to re-infection in the period of simulation (Figure S19B). For the N95 face masks, the same pattern of decay is observed as for surgical face masks: the only difference is the amount of reduction of cases for the unvaccinated (Figure S20) and the vaccinated population (Figure S21). In the case of natural immunity, the dynamics are different with respect to the other types of face masks. Between January and May of 2022, the pattern is the same, reaching 14,000,000 less susceptible individuals if NPIs are implemented (Figure S22A). For all types of face masks, the waning community only augmented in the following months, while for N95 face masks the waning population decreased by June 2022 until it increased by August 2022. This decay is achieved by the combination of NPIs and the usage of face masks with higher efficiency (Figure S22A). For acquired immunity, we observe an acceleration of decay when NPIs are implemented (Figure S22B).

5.6.   Assessing the combined effect of face mask usage, face mask efficiency and efficiency of vaccination in reducing the basic reproduction number

In this section, we evaluate the importance of the use of face masks and their efficiency mixed with the vaccination rate to reduce the basic reproduction number. Since there is a reduction in the population with only one dose (Figure 6D), we do not contemplate the importance of this dose in the reduction of $\mathcal{R}_0$. We can observe a maximum for the basic reproduction number of $\mathcal{R}_0 = 3.2$ during the period of simulations between January 4th and November 2022, and a baseline value of $\mathcal{R}_0 = 2$, which means that, despite a reduction of cases, the spread of the disease will maintain a continuing trajectory if the vaccinated population maintains only a two-dose scheme (Figure S23). In the scenario of a low vaccine efficiency, the only face mask type that will help to reduce the value of $\mathcal{R}_0$ is N95. We need a combination of 84% of the population using face masks, but only a face mask with high efficiency (N95). As we increase the percentage of face mask usage, the value of $\mathcal{R}_0$ can descend to values less than one, meaning that the Omicron wave can be effectively controlled (Figure S23A). For baseline vaccine efficiency, we have $\mathcal{R}_0$ = 2 when 80% of the population uses face masks with a global face mask efficiency of 80%, which can be achieved if half of the 80% wears a surgical face mask and the other half an N95 (82.5% face mask efficiency). If we increase to 90% the use of face masks, we obtain $\mathcal{R}_0$ = 1.2, which means that the pandemic is slightly controlled (Figure S23B). A higher vaccine efficiency means that we need a smaller proportion of individuals wearing face masks to reduce the value of $\mathcal{R}_0$. In this case, we need roughly 75% of the population mixed with a face mask efficiency of 78% (Figure S23C). By increasing the use of face masks, the value of $\mathcal{R}_0$ will descend to a value less than 1, where the Omicron wave is controlled. The application of the booster dose in the population will significantly reduce the value of $\mathcal{R}_0$; moreover, when combining this with the use of face masks, the value will descend much faster. In a scenario of a low booster efficiency, it is enough to have roughly 60% of the population using face masks to obtain a value of $\mathcal{R}_0$ = 2 in combination with a 65% of face mask efficiency. This means that surgical face masks are enough to obtain a value of $\mathcal{R}_0$ equal to 2.

If the American population chooses to use N95 face masks, a face mask usage of at least 60% will be needed to reduce the value of $\mathcal{R}_0$ to 1.8 (Figure 12A). For a baseline booster efficiency, a compliance of 57% of the population will be enough to reduce the value of $\mathcal{R}_0$ to 2 in combination with a 58% of face mask efficiency. This means that individuals can use three types of face masks: N95, surgical and improved cloth face mask (lower proportion of individuals using this type) (Figure 12B). With a compliance of 91% or higher, the value of $\mathcal{R}_0$ will be less than one, which implies that this wave has been controlled. Finally, the higher booster efficiency is associated with less than 40% compliance to reduce the value of $\mathcal{R}_0$ to 2. If 80% of the population decided to use face masks in combination with a face mask efficiency of 38%, the value of $\mathcal{R}_0$ will be less than one and the Omicron wave will be eliminated. With a higher booster dose efficiency, all face mask types can be used to decrease the value of $\mathcal{R}_0$, but if higher mask efficiency is used, the descent will be much more significant (Figure 12C).

5.7.   Evaluating when we can start to reduce the use of face masks to control the wave of the Omicron variant

In this section, we evaluate when we can stop using face masks. We started the simulation by assuming that 50% of the population wears face masks (all of the same type); then, by April 2022 only 25% of the population wears face masks; finally, by October 2022 no one is wearing any type of face mask. We evaluated the dynamics of new infections (symptomatic or not), hospitalizations and deaths from the unvaccinated and vaccinated populations. We vary the efficiency of face masks and compare it with the basic model where the use of face masks is continuous throughout the timespan of simulation (black dotted line).

If the unvaccinated population only wears cloth face masks, the peak of symptomatic infections will be in mid-January 2022, and cases will descend until reaching a valley in April 2022 with less than 400,000 infections. When only 25% of the population uses cloth face masks, cases will increase at a higher pace compared with 50%. The peak of the other wave in this year will be the same, reaching 600,000 cases, but the difference with 50% is that it takes much longer to reach this peak compared with 25% of the population. By October 2022, if no individual is using a face mask, the number of cases will increase by roughly 100 000 symptomatic infections (Figure S24A). The same dynamics are observed in asymptomatic infections (Figure S24B). If 50% of the population wears a cloth face mask, the peak of 60,000 will be reached at a slower pace until July 2022. In the case of 25%, the peak will be the same amount, but it will be reached sooner due to the difference in cases between April and May 2022 (Figure S24C). Deaths are not significantly different until May 2022 (Figure S25). There is a significant reduction of almost 10,000 cumulative deaths than can be prevented in the population of unvaccinated individuals (Figure S25, left panel). For the vaccinated population, the trend is the same, the only difference is between April and June 2022 there will be a slight increase in cases when only 25% of the population wears a face mask compared with 50% in symptomatic or asymptomatic infection (Figure S26A and B). The same pattern is observed in hospitalizations: there will be a slight increase of 75 more admissions in the population with only 25% using a face mask (Figure S26C). For the death toll in the vaccinated population, there will be no difference until April (Figure S26D); once April passes, there is only a slight reduction (about two deaths) if Americans continue to use face masks (Figure S25, right panel). When only 25% of the population uses face masks, the waning of natural immunity community increases by more than 1,400,000 individuals (Figure S27A), while the waning of acquired immunity is not affected (Figure S27B).

For improved face masks, the peak will be around the second week of January 2022, and cases in the unvaccinated population will descend to a valley for both symptomatic and asymptomatic infections (Figure 13A, B). The peak is reached at a higher pace when only 25% of the population uses an improved cloth face mask compared with 50% of the population, for both symptomatic and asymptomatic infections (Figure 13A, B). In October 2022, when no individuals use face masks, cases will increase significantly (Figure 13A, B). Hospitalization admissions have the same acceleration pace when only 25% of the population use a face mask (Figure 13C). The difference in the death toll can be observed until May 2022 (Figure 13D), but there is a difference of roughly 10,000 preventable deaths if 50% of the population still wears an improved cloth face mask (Figure S28, left panel). For the vaccinated community, there is a slight increment between May and June 2022 for asymptomatic and symptomatic infections and new hospital admissions (Figure S29A–C). The death toll is reduced by only two deaths if 50% of the community of the vaccinated uses a face mask (Figure S28, right panel). Waning natural immunity is affected by the reduction in face mask usage, generating a higher number of individuals susceptible to re-infection in the following months (Figure S30A). The waning of acquired immunity is not affected by the use of face masks (Figure S30B).

When we compare with the use of surgical face masks, the dynamics keep the same pace, but the peak is slightly higher with respect to the case with 50% face mask usage, reaching a peak of 620,000 for symptomatic and 14,000,000 for asymptomatic infections (Figure S31A, B). Hospitalization increases at a higher pace when only 25% of the population uses a face mask, compared with 50% of the population (Figure S31C). The death toll is significantly different with 20,000 preventable deaths during May and August 2020 (Figure S32, left panel). For the vaccinated population, the dynamics are affected in an analogous manner: an increase of cases is observed in the second wave of this year between May and June 2022 when 25% of the population uses a surgical face mask (Figure S33A–C). Roughly five deaths will be prevented if the use of face masks is maintained at 50% (Figure S32, right panel).

The same dynamics are observed for the use of N95 face masks. The pattern for the unvaccinated population can be seen in Figure S34, and the increment of deaths when only 25% of the population use a face mask compared with 50% is shown in Figure S35. The dynamics of the vaccinated population for this case are observed in Figure S36.

For the surgical and N95 face masks, cases will start to descend between January and March 2022 (Figure 13A). We can observe that the descent of cases will take a period of three months until reaching a valley in April 2022, and there is a reduction of 30,000 infected cases compared with no variation in the valley. If only 25% of the community uses face masks by April 2022, there will be an increase of almost 32,000 cases between June and August 2022 compared with no variation (Figure 13A). For asymptomatic infections, the pattern is the same, the use of face masks is still going to be important despite the transmission being lower (Figure 13B). Hospitalizations are affected as well: by April 2022 the transmission will be low, and hospitalizations will be at their lower count. If only 25% of the population uses face masks by the beginning of April 2022, there will be an increment of roughly 3200 admissions (Figure 13C). It is important to notice the speed of the increment when only 25% of the population uses face masks regardless of the face mask efficiency compared with 50% of the population using face masks (Figure 13A–C; black dotted line). When there is no change in the use of face masks, the increment of cases will be high, but to a lesser extent compared with the case with variation, moreover, the increment will be at a slower pace (Figure 13A–C). The death toll is heavily affected by using face masks. By April 2022, the use of N95 will be associated with a significant decrease of 20,000 deaths with respect to the use of improved or not improved cloth face masks (Figure 13D). Wearing face masks when cases are increasing will still be important to reduce the increment of cases. For the vaccinated population, the same behavior is observed (Figure S24). The waning of natural immunity is heavily affected by the use of face masks, where the use of the N95 will have a lower number of individuals susceptible to re-infection compared with the other types of face masks (Figure S25A). The acquired waning immunity is not affected by the variation of the use in face masks, this population will decrease due to the protection provided by the booster dose (Figure S25B).

5.8.   Assessing the importance of using face masks at the end of 2022

In Section 5.2, we evaluated the importance of face masks at the beginning of the year and how they can reduce the spread of the virus. However, we just evaluated different efficiencies based on the face masks, but we still need to simulate which type of face masks is better and how the proportion of individuals using them affects the dynamics of the outbreak. The pandemic has not only had effects on the social behavior of humans but also on the economy ^[45,46]. Taking this into account, we compared the baseline use (absence of face masks) with an increase of 10% and 20% in three types of face masks: cloth, surgical, and N95.

As expected, the face masks with the largest reduction in cases in all simulations are N95, surgical, and lastly, the cloth face mask (see Figure 14). If 10% of the population used the cloth face mask (red dashed line), the cases would decrease by 5% around mid-December, when the peak occurs. When we compare an increase of 20% in cloth face masks compared with the baseline use, we obtain a decrease of 10%, as well as a slight reduction in cases at the beginning of 2023. However, we want to focus on the behavior of the moderate and the high-efficiency face masks: surgical and N95, respectively. If we increase the high-efficiency face mask by 10% (turquoise dashed line), the peak of mid-December will be reduced roughly by 11%, and the peak will move slightly with respect to the baseline (black line). If we increase the moderate-efficiency face masks, the peak will move around a week compared with the baseline usage of face masks, but the cases will drop by 10% (magenta dashed line) compared with the baseline line. This means that both types of face masks have almost the same behavior in reducing the number of cases in the last months of 2022. If 20% of the population wears moderate-efficiency face masks, the number of infected individuals will be reduced by almost 14% (magenta solid line). Simultaneously, if we increase in 20% the use of N95 face masks, the number of infected individuals will be reduced by 18% with respect to the baseline. However, taking into account how the economy has been affected globally by this pandemic, medium-efficiency masks are relatively cheaper than high-efficiency masks. Therefore, we can see in the figure that if more than 20% of the population uses a surgical mask, it has a greater decrease compared to the case when 10% uses a more effective mask. The dynamics of reduction due to face mask usage keep the same pace when we evaluate the dynamics in hospitalized people and the death toll of individuals by COVID-19. In summary, we highly recommend the use of face masks in the last months of 2022 and the first months of 2023, because they are still the best option to reduce the spread of the virus without affecting the social and economical behavior of humans. Further, we see that it is better to have more individuals using moderate-efficiency (surgical) face masks than fewer individuals using high-efficiency face masks. These results demonstrate the wrong assumption that surgical masks were not as efficient in reducing the spread of the virus. The more people use them, the more the dispersion is reduced.

5.9.   Determining the effect of booster vaccination rate, waning rate, and the use of face masks on the control reproduction number

In this subsection, we perform some simulations for the control reproduction number ($\mathcal{R}_c$) using the expression obtained in Subsection 3.2, considering that there is a certain population susceptible to the virus, vaccinated individuals with one, two, or three doses, and individuals whose vaccine protection has waned. We focus on two comparisons: booster vaccination rate vs use of face masks, and booster vaccination rate vs waning rate.

The first comparison means that the reproduction number is considered as a function of the booster vaccination rate $\rho_4$ and the use of face masks ($C_m$). We simulated the reproduction number based on three types of face masks: cloth, surgical, and N95 face masks, to optimize which percentage and booster rate will result in a decrease of the control reproduction to be less than unity. Before describing the results, it is important to mention that the application of the fall booster in the US is low when roughly only 1% of the global population has received this booster ^[47]. This being said, we need to depend on the use of face masks to slow down the spread of the virus at the end of 2022. In the left contour plot in Figure 15, we can visualize that the value of $\mathcal{R}_c$ can be less than one if a relatively high booster rate is delivered daily, and the use of cloth face masks must be high. We need to apply 481,255 daily doses of the booster dose to achieve a value less than one. If we apply 746,775 doses combined with the use cloth face masks, the value of $\mathcal{R}_c$ will be 0.5, which means that the spread of the virus will be very low, eventually meaning that the pandemic will be over. However, if we used a surgical face mask that is associated with moderate protection, the dynamics will be different. We need the same vaccination rate but with 48% of the population wearing face masks compared with 100% using cloth face masks (middle contour plot). If we enhance the use of face masks mixed with a booster vaccination rate of 481,255 daily doses, the value of $\mathcal{R}_c$ will be less than one, associated with a low spread of the virus in the community. Finally, when we compare the use of the N95 face masks (high protection), with the same vaccination rate but only 35% of the population using face masks, $\mathcal{R}_c$ will be less than one, implying the end of the pandemic. However, we highly recommend the use of surgical face masks because they are affordable and provide good protection to avoid getting infected and infecting others. Nevertheless, we need to continue with the application of the fall booster because individuals are still dying of COVID-19 ^[48], despite the fact that there is wide access to the vaccine. Vaccines are effective if the individual is up to date with the vaccination scheme: most of the associated deaths are individuals that only have one booster ^[48].

The waning rate of the protection provided by the vaccine is an important variable that will eventually shape the dynamics of the spread of the virus. To evaluate this variable, we calculated $\mathcal{R}_c$ based on the booster rate and the waning rate (Figure S37). The waning rate will shape the dynamics of the spread of the virus when we compare the use of three types of face masks in the contour plot in Figure S37 (left: cloth; middle: surgical and right: N95). The smaller the waning rate (bottom part of the contour plot), the more days the protection acquired by the vaccine lasts, whereas the larger the waning rate (upper part), the fewer days the protection acquired by the vaccine lasts. In the case of protection by cloth masks, if the waning rate is less than 0.002 (this means that the protection acquired by the vaccine lasts roughly 500 days), the value of $\mathcal{R}_c$ will be less than 1.5, regardless of the vaccination rate. If the waning rate is further reduced, the value of $\mathcal{R}_c$ will decrease to less than unity. However, the pandemic has demonstrated that the protection from the vaccine is lower (this may be because of new variants emerging ^[49]), and the only way to reduce $\mathcal{R}_c$ is by keeping a high booster vaccination rate. For the case of surgical masks, if the protection of the vaccine is roughly 333 days (which may be accurate) and is combined with a low booster rate, $\mathcal{R}_c$ may be reduced to a value less than 1.5. A higher booster rate and a low waning rate have the possibility of reducing its value to less than one. For the use of N95 face masks, if the waning rate is less than 200 days, the value of $\mathcal{R}_c$ will be reduced to 1.5, which is sufficiently low. Also, a higher booster rate is associated with a value of $\mathcal{R}_c$ less than unity. In conclusion, the waning rate of acquired protection is an important parameter that we need to watch closely in the following months of the COVID-19 pandemic, since this parameter is decreasing as variants emerge. The only way to slow down the emergence of new variants without lockdown is if a relatively high proportion of individuals use face masks, particularly in closed areas where a healthy distance is not possible.

6.   Sensitivity analysis

In this section, we focus on the sensitivity analysis for our model. The local sensitivity analysis helps us to evaluate the impact on model output when the parameters are perturbed but remain close to their baseline values. On the other hand, global sensitivity analysis allows us to investigate how the model solutions change when the parameters vary over a large range of values. In this work, we will apply the local analysis for the basic reproduction number, and a global technique to analyze the variation in deaths, hospitalizations and symptomatic infected cases.

6.1.   Local sensitivity of the basic reproduction number

We will now perform a local sensitivity analysis to determine the impact of each of the model parameters on the basic reproduction number $\mathcal{R}_0$.

If $\mathcal{R}_0$ is differentiable with respect to a given parameter $\vartheta$, then the normalized forward sensitivity index of $\mathcal{R}_0$ with respect to $\vartheta$ is defined as

Using the baseline parameter values given in Tables 1–4, we can compute the values of the sensitivity indices with respect to $\Lambda$, $d$, $C_m$, $\epsilon_m$, $w$, $\beta_1$, $\beta_2$, $p_1$, $\alpha_1$ and $\gamma_1$. Thus, we obtain

The parameters with the largest influence on the basic reproduction number are those that have the largest sensitivity indices in absolute value, that is, $d$ and $\Lambda$. However, these two are parameters related to the demographics of the model, so they are unlikely to change. Hence, the most important parameters to control the value of $\mathcal{R}_0$ are $\gamma_1$ (the recovery rate of unvaccinated population) and $\beta_2$ (the disease transmission rate of asymptomatic infectious individuals), and the least important is $w$, the rate of transfer from the latent to the intectious compartment (see Figure 16).

6.2.   Global sensitivity

Finally, we will use the extended Fourier amplitude sensitivity test (eFAST) to perform a global sensitivity analysis. For a parameter $i$, this method computes the first-order sensitivity index $S_i$ by varying only $i$ and leaving all other parameters constant. The index $S_i$ represents the fraction of model output variance explained by the input variation of the parameter $i$ ^[50].

We used the Matlab code available in ^[51] to perform eFAST for model (2.1) using 65 values for each parameter per search curve and repeated this process for 10 different search curves. We computed the solutions of the system using the variable values for January 4, 2022 as initial conditions. We assumed that the recruitment rate $\Lambda$, natural death rate $d$ and vaccination rates $\rho_1$, $\rho_2$ and $\rho_3$ have fixed values as in Section 4, while all other parameters are allowed to vary. A "dummy parameter'' was included in the eFAST so we can compute the p-values for the model parameters by comparing their sensitivity indices with that of the dummy variable.

Sensitivity indices were computed for the following time points: 20, 60,100 and 300 days after the initial date. Figure 17 shows the values of the first-order sensitivity indices ($S_i$) for the death toll $D(t)$. Only the most significant parameters (${\text{p-value}} < 0.05$) are shown. We can see that the recovery rate of unvaccinated infected population $\gamma_1$ is the most significant parameter at 60,100 and 300 days. However, at 20 days, the parameters $\chi_B$, $\alpha_B$, $p_1$ and $p_2$ have a larger sensitivity index.

Figure 18 depicts the values of $S_i$ for the total number of hospitalized people ($H_1(t) + H_B(t)$). In this case, $\omega$ (the inverse of the average hospitalization time) is the most significant parameter. The parameters $\gamma_1$, $\omega_B$ and $p_1$ were also significant at all time points, although with a smaller sensitivity index.

Lastly, the values of $S_i$ for the cumulative number of symptomatic infected cases can be seen in Figure 19. Here, the proportion of symptomatic cases for unvaccinated population $p_1$ has the largest significance. The transmission rate for asymptomatic people $\beta_2$ has the second largest significance at 20, 60 and 100 days. The immunity waning rate $K_n$ has little significance in the short term (20 and 60 days), but it becomes more significant at 100 and 300 days.

7.   Discussion and conclusions

The year 2021 was marked by the emergence of new variants of COVID-19, mixed with vaccination roll-out and reduction of non-pharmaceutical strategies, therefore increasing transmission and enhancing the appearance of new variants, as it occurred at the end of November 2021 ^[52]. Vaccination against COVID-19 has significantly reduced the spread of the virus, but it has not been sufficient to eliminate this pandemic, one of the main reasons being vaccine hesitancy in the population ^{[53,54,55,56]}. Booster doses of COVID-19 were applied to all individuals of age 18 and above in September 2021, but not all fully vaccinated individuals got their booster ^[57]. Since the first case of Omicron in the US, the population was divided in terms of vaccination into three major groups: unvaccinated, vaccinated with booster and vaccinated without booster, mixing this population with the non-usage of face masks creates insufficient protection against the Omicron wave ^[58].

There are a great number of mathematical models that explain the dynamics of Omicron in selected countries. For example, Oh et al. ^[59] developed a complete set of ODEs that includes three vaccination levels (one dose, full dose, and booster) and a compartment for a non-pharmaceutical intervention (NPI): quarantine. They evaluated multiple types of common SEIR models and obtained that the full model performed better in describing the population with respect to others. They also found that despite the vaccination and NPI strategy, COVID-19 cases will increase in South Korea at the beginning of 2023. Furthermore, the study in ^[60], which did not include vaccination nor NPIs, found that the force of infection, the probability of infectiousness with the Omicron variant and the period of incubation are the parameters that affect the reproduction number by increasing its values. A paper that includes vaccination and NPIs applied in the US found that NPIs are essential to slow down the spread of the Omicron variant, despite having a population with a suitable number of vaccinated individuals ^[61]. Recently, some researchers have studied the Omicron wave in the US ^[62,63,64]. Two of these models focused on the competition dynamics between two strains of the virus, where one of them is Omicron. In ^[62], the authors developed a mathematical model that recognized vaccination, waning immunity of vaccinated individuals and waning of hybrid immunity. The main result is that Omicron will have the ability to increase the death toll despite having a low fatality rate. Moreover, future new strains having a low fatality rate may behave in the same manner as Omicron. Meanwhile, Ngonghala et al. ^[63] developed an ODE model to evaluate the competition between Delta and Omicron in early 2022, and how increasing vaccination and use of face masks will be essential to slow down the spread, even with the capacity to disappear the pandemic. Treatment against COVID-19 has the same dynamics as vaccination and has a significant decay against hospitalization than daily cases. A recent Omicron model ^[64] stated that boosting vaccination will help to reduce new Omicron cases, and that waning immunity is important for future waves. In comparison with the previous works, our model includes a variable that most models do not include: the use of face masks as a more suitable approach for NPIs based on what we have learned throughout this pandemic. Further, we include the three types of vaccines applied in the US, with their respective immunity and how it decays. We acknowledge the importance of booster vaccination and the use of face masks to decrease the reproduction number and try to control the spread of this virus and lower the possibility of a new variant emerging due to high transmission.

In this work, we carried out a theoretical and numerical analysis for an epidemic model that can be used to analyze the dynamics of the Omicron variant and the impact of vaccination campaigns in the United States. Our model is related to some previously studied COVID-19 models, such as those in ^[65,66,67]. However, the present model introduces several novelties, such as the inclusion of a hospitalized compartment and vaccine booster doses, and modelling the waning of natural and acquired immunity as a two-stage process. We also introduced a recruitment rate and a natural death rate, which make the model more suitable to study the long-term dynamics of the epidemic. Furthermore, this is the first model incorporating all of these features that has been applied to predict the development of the Omicron wave and subsequent COVID-19 waves using parameters fitted to current US data.

We investigated the influence of face mask usage in the population and face mask efficiency, which are two parameters with a large impact on the basic reproduction number. We demonstrated that N95 face masks are the most effective to reduce the spread and eliminate the pandemic, but surgical masks are important as well. The least effective types of face masks are cloth and improved cloth face masks, as expected. The pandemic has affected the economy, and not everyone has access to buy an N95 due to its cost, so we simulated a case where most individuals use a surgical mask and only a few use the N95, which is effective to reduce the spread significantly compared with only a fraction of individuals using N95 face masks. We assessed a very important question: when can we stop wearing face masks? We highly recommend that the use of face masks continues until August 2022 to reduce the spread, but this can change based on the dynamics of the spread in the following months, mixed with a possible emergence of a new variant.

The application of the mRNA vaccines and the adenovirus vaccine in the US has helped significantly in reducing hospitalizations and deaths, even though they were developed to fight the original strain of COVID-19. However, the cross-protection efficiency is waning over time for vaccine-acquired and natural immunity ^[68,69]. Our simulations showed an increase in new infections when we vary the immunity period between 7 and 11 months, with the possibility of another Omicron wave between June and July 2022. The predicted peak will be 53% lower than the peak in January 2022. There is also a possibility of another Omicron wave by the end of November 2022, with a 25% reduction with respect to the wave in June-July 2022. The possible upcoming waves and the emergence of another SARS-CoV-2 variant may be decisive in the application of a fourth dose mixed with the use of face masks to end the pandemic, as Israel is already applying ^[70]. We must focus our attention on modelling the dynamics of the subvariant BA.2, which apparently has a higher transmission rate and immune evasion from the protection provided by the vaccine ^[71]. In addition, we must pay close attention to the sudden increase of COVID-19 cases in China that is happening in the last week of March 2022. This increase is usually associated with a new variant: it has happened in Great Britain, India and South Africa with Alpha, Delta and Omicron, respectively.

Acknowledgments

Ugo Avila Ponce de León is a doctoral student from Programa de Doctorado en Ciencias Biológicas of the Universidad Nacional Autónoma de México (UNAM). This paper was developed in the period of his Ph.D. studies. Ugo Avila Ponce de León also received a fellowship (CVU: 774988) from Consejo Nacional de Ciencia y Tecnología (CONACYT). This article was supported in part by Mexican SNI under CVU 15284.

Conflict of interest

The authors declare there is no conflict of interest.

Data and materials availability

The original contributions presented in the study are included in the Supplementary Material. All code used for analysis is available at https://github.com/UgoAvila/Mathematical-Model-SARS-CoV-2-booter-dose-waning-immunity-in-the-United-States.