Modeling the SARS-CoV-2 sublineages XBB and BQ.1 in Mexico, considering multiple vaccinations, booster dose, waning immunity and cross-immunity

1.   Introduction

Since the beginning of the pandemic of COVID-19, at least 771 million confirmed cases and 6.97 million deaths have been reported (November 2023) ^[1]. In particular, Mexico has reported 7.6 million cases and a death toll of 334,752 ^[2]. In 2020, Mexico entered a social distancing period called "Jornada de Sana Distancia" from March all the way through June 2020 as a non-pharmaceutical intervention that helped reduce the transmission of the virus ^[3]. Although social distancing requirements were relaxed by June 2020 through a system of traffic lights, the use of face masks in closed and open spaces was always "mandatory" ^[4]. Two years have passed and the use of face masks in Mexico is still being applied by the population. At the end of 2020, the "Comisión Federal para la protección contra riesgos sanitarios" (COFEPRIS) approved the use of a highly efficacious vaccine against COVID-19 (developed by Pfizer) ^[5], which was applied to our first line of defense, consisting of doctors and nurses. The strategy for vaccine application in Mexico was the same used by other countries, which focused on vaccinating the elderly to avoid more fatalities, as recommended by a series of mathematical models ^[6,7]. The strategy was divided into six stages: vaccines would be first applied to workers associated with healthcare, then teachers, followed by adults divided by age: older adults, 50–59 years old, 40–49 years old, and the rest of the population ^[8]. Notwithstanding the similarity of Mexico's vaccination plan against COVID-19, the difference with other countries is the heterogeneity of vaccine brands being applied. Mexico has signed and applied 7 types of vaccine brands with three types of technology for the making of the vaccine. We have the mRNA vaccines Pfizer and Moderna; the non-replicating viral vector vaccines CanSino, Sputnik V, AstraZeneca, and Johnson & Johnson; and the inactivated vaccine Sinovac ^[5].

There have been some mathematical models that evaluated vaccination programs ^[9,10,11,12]. In ^[9], Alvarez et al. used a simple system of ordinary differential equations (ODEs) to assess the importance of a high vaccination rate instead of a high coverage in reducing the number of infected individuals. In ^[10], Munguía-López et al. focused on an optimization strategy to guarantee a fair allocation of the distribution of the vaccines to all states in Mexico under different fairness schemes. Soria-Arguello et al. ^[11] analyzed how to distribute the vaccine based on several transportation modes to guarantee a correct distribution to all the states in Mexico. Saldaña et al. ^[12] developed a mathematical model where they included the mobility of a population, vaccine coverage, efficacy, and delivery time. The study found that by reducing the mobility between patches, the burden of infections will diminish the same as with a 30% immunization in the population. They also showed that by increasing vaccine coverage between 20% and 50%, the reduction of infection can be up to 30% and 50%, respectively. However, none of these works have focused on the heterogeneity of vaccines in Mexico based on their efficiency, which is one of our main concerns in this study.

One year after the start of the pandemic (2021), the dynamics of the spread of the virus changed because of the implementation of the vaccines and the emergence of new variants which has not stopped since 2021. Based on their varied levels of infectiousness, lethality, and response to the vaccine, the most threatening ones are identified as variants of concern (VOCs) ^[13]. The first VOC that appeared was the Alpha variant (B.1.1.7 under the Pango lineage) that gained prevalence in the first half of 2021 ^[14]; nevertheless, this variant failed to become the dominant variant in Mexico because it could not compete with the Mexican variant that circulated until May 2021 ^[15,16]. Instead, the Delta variant did not have the same behavior as Alpha, and it became the dominant variant in a short time until early December 2021 ^[17]. The Omicron variant was first detected in South Africa at the beginning of November 2021 ^[18,19]. This variant was the cause of the rise in new cases at the beginning of 2022 in the whole world. The Delta variant became dominant in most countries within two weeks since it was first detected due to its high transmissibility and higher immune escape from the protection provided by the vaccine or natural infection ^[20]. Since the variant emerged, around 5 sublineages have appeared, with BA.1 being responsible for the first wave of this variant. In the middle of 2022, the second wave of cases emerged with the sublineage BA.5 being responsible for this augmentation. This variant has mutations in the spike protein, developing a higher capacity of invading the host cell and evading the neutralizing function of antibodies ^[21,22]. Moreover, people who got infected in the first wave can still get re-infected with this strain ^[23]. In Mexico, the Omicron sublineage XBB.12.1 was dominant in the first two weeks of June 2022, and recent data shows that BQ.1 is the current strain that predominates the new infections ^[24]. At the end of 2022, two new sublineages: XBB and BQ.1 emerged and competed between them to see who will become the dominant variant.

Two new variables need to be incorporated to understand the dynamics of subsequent waves of Omicron: the waning immunity (either acquired through the vaccine or natural) and the existence of boosters. In Mexico, the first reinforcement or booster began to be administered at the beginning of December 2021 and the application was completed until February 2022 ^[25], and the second booster started in May 2022 ^[8]. It is important to notice that a second booster is the fourth dose for all those whose first series consists of two doses (Pfizer, Moderna, Sputnik V, AstraZeneca, and Sinovac) and the third dose for those whose first series consists of only one dose (Cansino and Johnson & Johnson). The purpose of this study is to implement a mathematical model to understand the impact of the boosters on the dynamics of new cases, keeping in mind the waning of said booster in the remainder of this year. We will still incorporate the widespread use of face masks in the model because in Mexico this NPI is still being used to curtail infections.

Recently, some mathematical models have been proposed to explain the dynamics of the Omicron wave using a deterministic system of nonlinear differential equations ^{[26,27,28,29]}. In ^[26], Benahmadi et al. used a modified susceptible-infected-recovered (SIR) model to include vaccination and fitted the parameters with the dynamics of Omicron in Morocco. They developed an optimized control analysis based on an imperfect vaccine to diminish the burden of new COVID-19 cases. In ^[27], Edholm et al. used a modified SIR model that included vaccination and breakthrough (individuals who get infected despite being vaccinated) exposed individuals in the province of Gauteng, South Africa. Their main result is that modifying the rates of the two different vaccines has the ability to affect the spread of the virus. Meanwhile, Gonzalez Parra and Arenas ^[28] developed an SIR model to assess the importance of vaccination and how the waning immunity and the efficacy of the vaccine are important parameters that affect the spread of Omicron in the United States. Finally, Safdar et al. ^[29] developed an extensive modified SIR model to acknowledge the presence of three types of vaccines applied in the US and how they can prevent more new COVID-19 infections. In this model, the authors also included the notion of a booster (third dose for mRNA vaccines and second dose for Johnson & Johnson) and studied how a high percentage of the booster is needed to prevent more infections at the beginning of 2023. In the present work, we developed a mathematical model of the same characteristics that incorporates the perks of the vaccination program in Mexico and the current knowledge of the changing epidemiology of COVID-19 ^[30], including parameters that describe the waning of vaccine-derived and natural immunity, as well as the use of face masks with their different efficiencies and mixed cross-immunity between Omicron sublineages. To our knowledge, there are no current mathematical models that describe the dynamics of the second Omicron wave in Mexico and provide short-term predictions of that wave.

2.   Materials and methods

In the following subsections, we describe the derivation of the mathematical formulas and parameter estimation for the model, including the parameters to explain the behavior of an imperfect vaccine. We also document the procedures for obtaining the control reproduction number, as well as the local sensitivity analysis.

2.1.   Derivation of the mathematical model of a partially vaccinated population of Mexico, considering the application of multiple vaccines and waning immunity for two sublineages of Omicron

Our model was developed by modifying a previously developed compartmental mathematical model that incorporated vaccination in Mexico ^[31] mixed with a model that incorporated multiple vaccines as an important variable ^[32]. Because we are modeling a vaccination campaign against COVID-19, we must incorporate concepts around an imperfect vaccine into our model, such as leakiness and the all-or-nothing protection, which are well explained in ^[33,34]. In view of the fact that the third Omicron wave in Mexico was first driven by the XBB and then by BQ.1, our mathematical model incorporates these two sublineages by adapting a previous compartmental model ^[35]. Finally, our model also incorporates the loss of immunity, either natural or acquired, and the existence of individuals with one or two boosters. The mathematical model contains two separate sets of equations for the Omicron variants XBB and BQ.1, which can affect the susceptible, vaccinated (regardless of the number of doses), and the partially protected susceptible individuals. Our compartmental models evaluate the dynamics of 16 populations at any given time $t$, which are denoted as $S(t)$, $F(t)$, etc. The flow between the compartments is represented in Figure 1.

Susceptible individuals will decrease when they interact with a symptomatic individual infected with BA.2 at a rate $\beta_1$ and with an asymptomatic individual at a rate $\beta_2$. Once exposed to the BA.2 sublineage, the infected individuals will become infectious after $1/w$ days. But only a proportion $p_1$ of these individuals will develop symptoms, whereas the rest of the exposed individuals will be asymptomatic. The symptomatic subpopulation can flow between two states: recovered or hospitalized. Symptomatic and asymptomatic individuals infected with BA.2 will recover after $1/\gamma_1$ days. However, in the case of symptomatic individuals, only a proportion $1-\varphi_1$ recovers from the disease; the rest of them become hospitalized after $1/\sigma$ days, and these individuals can recover at a rate $\omega_1(1-\chi_1)$ or die at a rate $\omega_1\chi_1$.

The second set of equations acknowledges the effects of individuals infected with BA.5, the flow between these compartments is the same as described for the BA.2. Partially protected susceptible individuals represent the compartment where immunity (either vaccine-acquired or naturally acquired) wanes and can be re-infected with the same or the other sublineage, it is depicted in Figure 1 as $S_p$. A full description of the model and the set of equations can be found in Appendix A1.

2.2.   Description of the parameters used as fixed values in the mathematical model of a partially vaccinated population of Mexico, considering the application of multiple vaccines and waning immunity for two sublineages of Omicron

To model the dynamics of the susceptible population, we use several parameters that need to be explained to better understand the dynamics of the spread of the two sublineages of Omicron. The first terms describe the dynamics of the population in question: $\Lambda$ and $\mu$ describe the natural recruitment and the natural death rate, respectively. The following parameters need to be included in our model because we are not considering a closed population: individuals migrate to the population ($\Lambda$), and they may die due to causes other than COVID-19 ($\mu$). Since we are considering two sublineages of Omicron, the susceptible population can be infected by both types of viruses; consequently, we have four transmission rates: $\beta_1$, $\beta_2$, $\beta_3$, and $\beta_4$. The first two parameters model the infection by sublineage one: $\beta_1$ represents when a susceptible individual encounters a symptomatic individual and $\beta_2$ represents when a susceptible individual interacts with an asymptomatic individual. Because asymptomatic people normally do not show symptoms, the probability of them infecting an individual is different than a symptomatic individual. Parameters $\beta_3$ and $\beta_4$ describe the same dynamics as the latter, the only difference is that the susceptible individuals can get infected by the sublineage two when they interact with a symptomatic and an asymptomatic individual, respectively. Because we are modeling an imperfect vaccine, as explained in the previous section, all individuals that get vaccinated at a rate $\rho$ may or may not be protected. We assume that only a proportion $\varepsilon_0\in[0, 1]$ of the people receiving the vaccine are protected. This is called the "all-or-nothing" protection and is defined by the expression

where $\varepsilon_{ai}$ and $w_i$ are the protection parameter and the relative vaccination rate, respectively, for each group of vaccines (high, medium, or low efficiency).

The next compartment consists of the fully vaccinated individuals, regardless of the manufacturer of the vaccine and the number of doses. Another concept related to vaccination is the vaccine effectiveness, which means that some individuals may get infected with either sublineage despite being vaccinated. This is what we will call "breakthrough infections". Because these individuals are vaccinated, they have developed a protection created by the vaccine, so if they encounter a symptomatic or asymptomatic individual infected with sublineage one, the probability of them getting infected is reduced by a rate $\mathcal{E}_1$. The same dynamic happens if they interact with an individual infected with sublineage two, the probability of infection is reduced by a rate $\mathcal{E}_2$. These parameters are defined as

where $\varepsilon_{1i}$ is the efficacy against sublineage one and $\varepsilon_{2i}$ is the efficacy against sublineage two, for vaccines in the $i$-th group ($i = 1, 2, 3$). Some fully vaccinated individuals will stay for a period in this compartment because of the protection provided by the vaccine; however, the acquired immunity will wane against time, and they transit to a partially protected susceptible compartment at a rate $\alpha_B$. In the fall of 2022, fully vaccinated individuals had the option of a booster dose, so this population may receive a booster at a rate $\rho_B$.

Booster individuals are the population that received an extra dose of vaccine to acquire more protection in the fall and winter of 2022; however, they can still get infected if they interact with infected individuals, regardless of their symptoms. Their probability of infection is reduced by a factor $1-\varepsilon_{B1}$ for sublineage one and $1-\varepsilon_{B2}$ for sublineage two. Eventually, the protection provided by the vaccine will wane over time, and these individuals will transit to a partially protected population at a rate $\alpha_B$.

Exposed individuals, regardless of the sublineage, are those that are exposed to the virus because they encountered an individual infected with either sublineage of Omicron. The flow of this compartment comes from the susceptible population, breakthrough infections from fully vaccinated and booster dose, and from partially immune susceptible individuals. Partially immune protected individuals can get infected from sublineage one of Omicron at a rate reduced by a factor $1-\psi_1$, and they can get infected by sublineage two at rate reduced by $1-\psi_2$. As the immunity wanes, people can still get infected by the virus; however, they already have a targeted immune response that will prevent hospitalization.

A proportion of infected individuals will develop symptoms at a rate $p_1w$ for sublineage one and $p_2w$ for sublineage two. Most infected individuals will recover at a rate $\gamma_1\left(1-\varphi_1\right)$; however, a portion $\varphi_1\in[0, 1]$ of these individuals will develop severe symptoms and can get hospitalized at a rate $\sigma\varphi_1$ for sublineage one. Meanwhile, for the sublineage two of Omicron, individuals will recover at a rate $\gamma_2\left(1-\varphi_2\right)$, and a portion $\varphi_2\in[0, 1]$ of them may develop symptoms and will need to be hospitalized at a rate $\sigma\varphi_2$.

In the hospitalized compartment for both sublineages of Omicron, individuals will eventually recover and transit to the recovery compartment for sublineage one at a rate $1-\chi_1$ multiplied by the average time of staying in hospital $\omega_1$. Unfortunately, a proportion $\chi_1\in[0, 1]$ of hospitalized individuals may pass away from the disease at a rate $\omega_1\chi_1$. Concurrently, hospitalized individuals from sublineage two will recover at a rate $1-\chi_2$, multiplied by the average time in the hospital $\omega_2$. In the same manner, a portion $\chi_2\in[0, 1]$ of hospitalized individuals from sublineage two may die from the disease at a rate $\omega_2\chi_2$.

For recovered individuals from a previous infection by either sublineage of Omicron, their natural immunity would wane off at rates $\eta_1$ and $\eta_2$ for sublineages one and two of Omicron, respectively. All baseline parameters are outlined with their respective value and units in Table 1. Meanwhile, parameters associated with infection ($\beta_1$, $\beta_2$, $\beta_3$, $\beta_4$), hospitalization ($\chi_1$ and $\chi_2$), recovery ($\varphi_1$ and $\varphi_2$), and vaccination rates ($\rho$ and $\rho_B$) are estimated based on data. The description of the optimization protocol is mentioned in the next subsection.

2.3.   Estimation of the parameters of the dynamics of the two most abundant Omicron sublineages in Mexico

In this section, we will estimate the parameter values for our model. First, we fit the model solutions to real data during the period from December 7, 2021 to April 20, 2022, which is when the first Omicron sublineage was dominant in Mexico. We assume that all reported infections during this period belong to this strain. Hence, the variables of our model corresponding to strain 2 will be set to zero. We set the fixed parameter values shown in Table 1. We estimate the values for $\beta_1$, $\beta_2$, $\varphi_1$, $\chi_1$, $\omega_1$, $\rho$, $\rho_B$, and the initial exposed population $E_1(0)$ by fitting our model to the daily COVID-19 infection, death and hospitalization data, as well as the cumulative numbers of vaccinations and booster doses. The optimization of parameters was performed in Matlab by minimizing the sum of squared errors (SSE), which was computed as follows. For a given vector of parameters ${\mathbf x}$, we obtain the numerical solutions of our model using the ode45 solver. Then, the sum of squared errors is given by

where $I_i^{mod}$, $A_i^{mod}$, $D_i^{mod}$, $H_i^{mod}$, $V_i^{mod}$, and $B_i^{mod}$ are the number of new symptomatic infections, new asymptomatic infections, new deaths, active hospitalizations, cumulative vaccinations, and booster doses obtained from the model for the $i$-th day, respectively, while $I_i^{rep}$, $A_i^{rep}$, $D_i^{rep}$, $H_i^{rep}$, $V_i^{rep}$, and $B_i^{rep}$ are the homologous quantities obtained from the reported data, respectively. Notice that the number of data points $n_1\ldots, n_6$ may be different for each variable. During this period, the number of boosters was reported only for the final date (April 20, 2022), so we used this single value to compute the SSE by setting $n_6 = 1$. The parameter vector that minimizes the SSE was found using three searches: a gradient-based method, a gradient-free algorithm, and again, a gradient-based method.

The data for COVID-19 cases and deaths were obtained from ^[36]. The data for hospitalized cases were taken from ^[37], and the data for the number of vaccinations and booster doses were taken from ^[38]. The initial values for all variables other than $E_1$ were chosen according to these data. The best fit values are shown in Table 2. Figure 2 shows a comparison between the solutions obtained with the model and the reported number of symptomatic/asymptomatic cases, hospitalizations, and number of COVID-19 deaths. Figure 3 shows a similar comparison for the vaccinated population, and Figure 4 depicts the model solutions for the other variables.

Next, we estimate the values of $\beta_3$, $\beta_4$, $\varphi_2$, $\chi_2$, $\omega_2$, and $\rho_B$, using the model simulations for both strains for two different periods (April 20 to July 15, 2022 and November 18 to December 31, 2022). The other parameter values are fixed as before. The optimization for the first period was carried out using the same method as before but considering the cumulative numbers of COVID-19 cases and deaths, as well as the daily hospitalizations and booster doses. The initial conditions were chosen using the values obtained in the previous simulations for April 20, 2022.

The best fit parameters for the two-strain model are shown in Table 3. Figure 5 depicts a comparison between the model solutions and the reported number of symptomatic/asymptomatic cases, hospitalizations, and number of COVID-19 deaths (for both strains). Figure 6 shows a comparison for the vaccinated population. Figure 7 depicts the model solutions for the susceptible, susceptible with variable immunity, and exposed populations. Figure 8 shows the simulated number of active infections and hospitalizations for each strain. These results were used to estimate the initial conditions for the optimization starting on November 18.

The comparison with reported data from November 18 to December 31, 2022 can be seen in Figure 9 (the number of vaccinations was not used in this case due to scarcity of reported data for this period), and simulations for both strains are depicted in Figure 10. The simulations show that the number of exposed, symptomatic, and asymptomatic individuals for strain 1 is almost zero for this period, while the number of hospitalized individuals is less than five for strain 1. This shows that strain 2 is dominating the pandemic dynamics from November 2022 onwards.

2.4.   Estimation of the control reproduction number

Our model always has a disease-free equilibrium, which is given by

where $S^0 > 0$, $F^0 \ge 0$, $B^0 \ge 0$, $S_p^0 \ge 0$ (the equalities hold in the case when there is no vaccination) and $E_1 = I_1 = A_1 = H_1 = E_2 = I_2 = A_2 = H_2 = 0$.

This equilibrium represents the scenario when the disease has become extinct, resulting in fixed fractions of the vaccinated, susceptible, and partially protected susceptible individuals. If the vaccination rate is positive, there will always remain a vaccinated portion and a partially protected portion of the population. Otherwise, if the vaccination program is discontinued, the entire population will become completely susceptible due to the waning rate of immunity.

We compute the control reproduction number $\mathcal{R}_c$ by applying the next-generation matrix method (see ^[39]) and using the equation $\mathcal{R}_c = \rho\left(FV^{-1}\right)$, where $F$ denotes the derivative of the matrix of new infections, $V$ is the transition matrix (flow between compartments), and $\rho$ is the spectral radius. The entire derivation of the reproduction number, as well as the theoretical analysis of the model can be found in Appendix A2.

2.5.   Local sensitivity analysis of the control reproduction number

The sensitivity analysis of the control reproduction number was carried out using the following definition, taken from ^[40]: the normalized forward sensitivity index of $\mathcal{R}_c$, which depends differentiably on a parameter $p$, is defined by

Using the estimated values for the parameters obtained in Section 2.3 as baseline values, we were interested in perturbing them to determine how their changes affect the reproduction number $\mathcal{R}_c$. We used the software Maple to compute the partial derivatives of the reproduction number with respect to each parameter and then evaluate the sensitivity index, following the formula above. A positive sign of this index correlates with an increase in $\mathcal{R}_c$, whereas a negative sign is associated with a decrease in $\mathcal{R}_c$.

3.   Results

3.1.   Modeling the SARS-CoV-2 Omicron sublineages variant spread in a multiple vaccinated population

In this section, we simulated our model to evaluate the behavior of the usage of face masks, vaccination and the importance of the booster and waning immunity in Mexico. We considered seven types of vaccines that are currently being applied in Mexico: Pfizer, Moderna, Sputnik V, AstraZeneca, Johnson & Johnson, Cansino, and Sinovac. To simplify our model, we grouped the seven types of vaccines into three groups based on their efficiency in preventing symptoms when a person gets infected with the virus. Pfizer, Moderna, and Sputnik V are grouped as high efficiency; AstraZeneca and Johnson & Johnson are considered medium efficiency, and finally, Cansino and Sinovac ^[41] consitute the low efficiency group. The simulations were carried out based on the estimated and fixed parameters presented in Tables 1–3.

Cases of infected individuals of the third wave of Omicron are starting to decrease, the peak was reached in the first week of January 2023 (Figure 11A). In the following months, cases of symptomatic individuals will start to descend, reaching a valley in the month of February for the XBB variant. In the last months of 2022 and early days of 2023, there will be a slight increase in cases, the peak may be reached in mid-March. Nevertheless, the peak will be 10 times higher than the peak reached between December and January 2023 with the XBB variant. The pattern of asymptomatic infections is the same as for symptomatic infections, the only difference is the number of cases. The peak was reached with less than 80,000 cases in the first week of January 2023, followed by a descent reaching a valley of zero active asymptomatic cases for the XBB variant, trailed to a peak in mid-December. The peak in the last month of the year will be 10 times lower compared with the peak that will be reached in March 2023 (Figure 11B). Notwithstanding an increase in cases, the number of hospitalizations will increase equally, but it will not be such a sudden increase that it exceeds the carrying capacity of medical supplies (Figure 11C). The peak was reached in the last week of December 2022 with a value of fewer than 2500 hospitalizations, the dynamics of descent will be the same as the infections and the peak of hospitalization may be until March of 2023 with roughly 250,000 hospitalizations, which will be more than 50% compared with the peak attained in the last week of the year (Figure 11C). Around 62% percent of the population of Mexico is fully vaccinated, and the trend of accomplishing the percentage to reach a higher population of vaccinated individuals is descending, which means that unvaccinated individuals will not be getting their dose of vaccine. Booster doses will continue to be applied in December with a descent in March 2023 (Figure 11D).

Now, we will review the behavior of the recovered, susceptible but partially protected individuals and the death toll by COVID-19. The recovery of individuals from XBB will continue to descend throughout the year, highlighting that it will never attain a value of zero-infected cases from this sublineage. In the meantime, the pattern of infections from BQ.1 will be decreasing for roughly two weeks and increasing in the following months, meaning that provided another variant or sublineage does not enter or emerge due to high transmission, BQ.1 will be circulating in the country going up and down in new infections (Figure 12A). In Figure 12B, we can acknowledge that the increase of new cases in the following months is caused by the dynamics of the susceptible but partially protected individuals, and that when this population decreases it will eventually be implicated in an increase of symptomatic or asymptomatic cases in the weeks following the decrease. When a period passes, the people who were not infected in the last wave of the year are the ones who will supply the next wave between the first months of 2023. The pattern of new infections will be dictated by this population, which derives from the fact that very few individuals will remain in the susceptible compartment without protection (Figure 12B). Because the cases of XBB are decreasing, the death toll from this sublineage will be zero; however, the pattern for BQ.1 is different (Figure 12C). By January 2023, there will be around 70,000 deaths because of their infection from BQ.1. If no other variant emerges in the remainder of the year, the death toll by BQ.1 sublineage may reach around 150,000 deaths in July 2023 (Figure 12C). This means that although the death rate is low, it will never be zero. Strategies must continue to be used to reduce transmission.

3.2.   Spread of the Omicron variant under heterogeneity of vaccination efficiency due to the application of multiple vaccines

Vaccines against COVID-19 have become the best strategy to decrease hospitalizations and deaths caused by this virus ^[42], and Mexico is not the exception. Less than 40% of the population in Mexico has not received a vaccine, which is still a significant amount, however, the population of susceptible but partially protected individuals is increasing because of the waning of natural or acquired immunity. The exposure rate of partially protected individuals is higher compared with the unvaccinated individuals (Figure S1), however, the probability to enter a hospital is greater in the unvaccinated population than in the partially protected population, regardless of the sublineage of Omicron (Figures 13 and S2). We separated the vaccinated population into two groups, those who received a booster and those who did not. The booster population has a higher exposure rate to the BQ.1 sublineage; however, their hospitalization rate is lower compared with the population that has not received any type of booster (Figure S3).

Finally, we wanted to evaluate how efficiently the three groups of vaccines that have been applied in Mexico have behaved in this second Omicron wave. The high-efficiency group had a lower exposed rate compared with the other efficiency groups, followed by the medium-efficiency group with a lower exposed rate compared with the low-efficiency group of vaccines applied in Mexico. This means that if a person is vaccinated with a low-efficiency vaccine, the probability of becoming exposed to the virus is higher than with other vaccines. Notwithstanding, this higher rate is not observed in the hospitalization rate: regardless of the vaccine, the proportion of hospitalization is relatively the same for the three groups of vaccines (Figures S2 and 13). For both sublineages, the hospitalization rate was higher for the susceptible non-protected individuals (Figures 13 and S2). Despite the heterogeneity of vaccines applied in Mexico, they have helped reduce the hospitalization rate, regardless of their vaccine efficiency. One flaw observed in the vaccines against COVID-19 is that they reduce the transmission between individuals, although not significantly.

3.3.   Assessing the impact of different combinations of COVID-19 vaccines

In this section, we will evaluate the combinations of vaccine efficiencies that can influence the infectivity of SARS-CoV-2 in Mexico. The vaccine with the highest number of doses applied in Mexico is AstraZeneca, followed by Pfizer ^[43]. Considering the three groups of vaccines applied in Mexico mixed with six combinations, we simulated the roll-out of new infections and hospitalization in this second Omicron wave. In Table 4, we list the combinations of the vaccines compared with the real strategy with their weighted average efficiency based on the administered doses. The first combination is the strategy applied by Mexico.

We evaluated the peak that occurred in July–August 2022 and the peak that may occur between November and December 2022. We evaluated whether the original vaccination strategy had been modified from the beginning based on the combinations proposed in Table 4. First, we evaluated the differences with the peak of the second wave in 2022. If Mexico had decided to vaccinate with only high efficacy, in general, there would be a decrease of 0.0259%, 0.03%, 0.0308%, and 0.0024% of symptomatic and asymptomatic cases, hospitalizations, and deaths, respectively (Figure S4). In contrast, using only a medium efficacy vaccine, there would be a slight increase of 0.0103% and 0.005% of symptomatic and asymptomatic cases, respectively. With a slight increase in hospitalizations, at 0.006%, deaths from COVID-19 would have remained the same. Vaccination with only low efficiency will produce an increase of roughly 0.5792% symptomatic and 0.0303% asymptomatic infections. In the case of hospitalizations, there will be an increase of 0.0352%, and deaths from COVID-19 will be associated with a 0.0024% increase for the second peak of the Omicron wave. For the other comparisons, there will be a slight increase in cases if they present no symptoms, as well as in hospitalizations (Figure S5). Only when we vaccinate more with a low-efficiency vaccine will there be an increase of 0.0012% in the death toll in mid-July. The peaks for each vaccine strategy are shown in Table 5.

Now, we will evaluate the behavior for the third Omicron wave in the last months of 2022. Overall, the behavior in terms of decay maintains the same pattern as the peak of mid-July. In relation to symptomatic cases, there will be a decay independently of the application of high, medium, or low-efficiency vaccines only and the decay will be 0.0119%, 0.0158%, and 0.0222% respectively. For hospitalizations, there will be a decay only when we apply high-efficiency vaccines, for medium and low efficiency there will be an increase of 0.004% and 0.0168%, respectively. Deaths from COVID-19 will remain the same using only a high-efficiency vaccine, while for medium and low efficiency there will be an increase of 0.0054%. Using all combinations is associated with an increase or decrease; however, the change of values is slight and we do not consider it significant. The peaks for each vaccine strategy are shown in Table 6.

3.4.   Evaluating the vaccine efficiency for the primary application and boosters in the rollout of infections, hospitalizations, and death from COVID-19

We evaluated the importance of different vaccine efficiencies for the primary set of doses and the booster, as well as how both sublineages compete based on the overall efficiency provided by the vaccine. For individuals who have not received a booster, there is a significant difference in the number of symptomatic individuals when we compare the low and high efficiency of the vaccines. If the vaccine has the highest efficiency, there will be a reduction of roughly 0.0013% of cases compared with the baseline efficiency. When we compare with low efficiency, there could be an increase of roughly 0.0016% in symptomatic infections (Figure 14A). The same pattern of increase and decrease occurs with the asymptomatic cases, hospitalizations, and death toll (Figure 14B–D). In terms of the competition between sublineages, vaccine efficiency helps in decelerating the increase of XBB infections, and this reduction was associated with an increase of BQ.1 infection. This means that XBB was the force for the third wave of Omicron in the last month of 2022 (Figure S6).

Booster doses behave in the same manner, that is, when the efficiency is high compared with the baseline, there will be a decrease in symptomatic cases (Figure S7A), asymptomatic cases (Figure S7B), hospitalizations, and deaths (Figure S7C, D). The opposite behavior occurs when we compare low efficiency with baseline efficiency. In terms of the competition between XBB and BQ.1, the same behavior occurs, XBB was the force to let BQ.1 thrive when the transmission was high. However, booster doses can reduce the peak that may occur in the first months of 2023 (Figure S8). This means that booster doses are more helpful in reducing the number of infected and hospitalized cases from the BQQ.1 Omicron variant.

3.5.   Assessing the importance of the use of a non-pharmaceutical strategy: face mask usage and efficiency are associated with shaping the behavior of the Omicron third wave

In this section, we evaluate the behavior of two important parameters: the proportion of the population that uses face masks ($C_m\in[0, 1]$) and the efficiency of the face masks ($\epsilon_m\in[0, 1]$) being used by said population. In Mexico, the use of face masks was not compulsory throughout the pandemic; however, most of the population used them and continue to wear them, specifically in closed areas. 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 S9, we evaluated three scenarios of the usage of cloth face masks based on the percentage of the population: from 25% to 75% of the population using face masks. If only 25% of the population uses a cloth face mask, the peak in March would have been roughly 466,130 symptomatic cases. In the meantime, if 50% of the population wears this type of face mask, there would have been a reduction of 18.74%, or roughly 87,340 fewer infected cases. Despite having a low efficiency in protection, the peak could have been reduced to only 298,570 active symptomatic cases if 75% of the population uses this face mask, a reduction of 35.95% of infected cases when compared with a 25% use of face masks. For the third wave in the last two months of the year, the difference in the percentage of use will dictate the time it takes to reach the peak. All three percentages will reach the same number of active symptomatic cases, and the difference is the length it takes to reach the peak. For asymptomatic individuals, if 75% of the population use face masks, the peak in March could have a 35.93% reduction and it will take roughly two months to achieve the peak in mid-April of 2023 (Figure S9B). Hospitalization is reduced in the same manner if only 50% of the population uses a face mask with low efficiency with a decrease of less than 14.16% admissions (Figure S9C). Deaths can be reduced by the beginning of 2023, with a minimum reduction of 15,401 deaths (28%) if 50% use this type face mask (Figure S9D).

If we reinforce our cloth face mask and add a tri-layer, which improves the efficiency to 50% of protection against the Omicron variant, the reduction is much more drastic. When 50% of the population uses it, the peak can be not only reduced but the time to reach the peak will be longer compared with a 25% use of this face mask. The symptomatic infections will be reduced with roughly 36.51% if only 50% of the population uses this face mask, compared with only 25% of usage. If 75% of the population uses the face mask, symptomatic active cases can be reduced by 363,575 infections, which is associated with a reduction of almost 90% compared with only 25% of usage (Figure S10A). The same pattern of reduction is observed in asymptomatic individuals (Figure S10B). Hospitalizations can be reduced by 99,495 ($-49.56\%$) or 180,743 ($-90.02\%$) if 50% or 75% of the population use this type of face mask, respectively, and the subsequent wave that may happen in the remaining of this year may be significantly reduced only if 50% or more people use this type of face mask (Figure S10C). The death toll is reduced if a minimum of 50% of the population uses this type of face mask (Figure S10D).

The most efficient face masks are the surgical and the N95, with 70% and 95% efficiency, respectively. However, the use of surgical face masks is slightly more accessible due to their cost compared with the N95. With only 25% of the population using this type of face mask, there would have been a significant amount of reduction than if no one uses it (Figure 15). By increasing to 50% the use of face masks, a considerable reduction of roughly 280,974 ($-89.99\%$) infections, which is a significant reduction of cases when the transmission is high (Figure 15A). For the upcoming wave of 2023, the peak with 50% of usage is lower than when 25% of the population uses a face mask, a reduction of 10,000 active symptomatic cases. If we compare the 75% use of this face mask, two important effects are seen: the peak is reduced and the time it takes to reach the peak is larger. The peak will be reduced by 331,722 ($-94.44\%$) active cases and it will be reached in roughly three months, enough time to reduce transmission and help gather medical supplies if needed. This pattern is observed in asymptomatic individuals with a reduction of the peak and 5 months to reach it (Figure 15B). For hospitalizations, the reduction is slightly more drastic. If 75% of the population uses a surgical mask, there will be 163,360 ($-95.94\%$) fewer hospitalized individuals in the peak for June 2023. In terms of the use of beds, this reduction is important to deal with other diseases that have been ignored because of this pandemic (Figure 15C). Deaths are more significantly reduced, with a decrease of roughly 30,477 deaths if 75% of the population uses a surgical mask and a decrease of 25,302 if 50% of the population uses it, compared with 25% of the population (Figure 15D).

Next, we will evaluate the use of the N95 face masks, but we consider that, due to their cost, a large part of the population does not use this type of face mask. However, the reduction of new infections is drastic. If 50% of the population of the population uses an N95 mask, the peak of mid-April could have been reduced to only 19,588 symptomatic cases, and if 75% of the population uses it, the peak would have been reduced by 274,470 infected symptomatic cases compared with 50% of usage (Figure S11A). With a 25% usage, the peak in the last month of the year will be only 280,460 cases; however, it would take only one month to reach this peak. Nevertheless, if 50% uses this type of face mask, the peak can be reduced by 150,000 infected cases, and it will take roughly two months to reach it, which is a considerable time to prepare for the subsequent wave in terms of gathering medical supplies. On the other hand, if 75% of the population uses this face mask, the subsequent waves can be avoided and new infections will approach a value close to zero, a value not seen since the pandemic began (Figure S11A). The same dynamics occur with asymptomatic cases: if more individuals use this type of face mask, the subsequent waves will smaller compared to the last one, and the time it takes to reach them will be longer (Figure S11B). Hospitalizations are reduced by 105,930 individuals if 50% use an N95, compared with a reduction of 107,430 hospitalized individuals if 75% of the population uses it (Figure S11C). Deaths are significantly reduced: if 75% of the population uses an N95, there will be no more deaths caused by COVID-19 since September, however, there is an important reduction if only 50% of the population uses this type of face mask, deriving in 17,748 fewer deaths compared with 25% of the population using this type of face mask. In conclusion, face masks are an important non-pharmaceutical strategy that still needs to be implemented to reduce transmission and the probability of the emergence of new more infectious variants.

3.6.   Varying the loss of acquired or natural immunity in a population that is continuously using face masks

In this section, we will evaluate how the loss of immunity once people get infected, regardless of vaccination status, will shape the behavior of the last wave in 2023. To better visualize the effect of this parameter, we divided the graphs into two periods, the wave in March (January–March 2023) and the wave that may happen in June (May–July 2023). If individuals recovering from COVID-19 developed lifelong protection by December 2022, the pandemic will not become extinct in that period: the number of infected symptomatic cases will be low, with only 296,110 cases (Figure S12, pink dotted line). However, as individuals can get re-infected from this virus, we varied the average immunity period between 50 days (red dotted line) and 150 days (black dotted line) with a baseline protection of 100 days (blue dotted line). The simulation was carried out assuming that a minimum of 50% of the population wears a face mask in the period evaluated, with an overall efficiency of 50% (improved cloth face mask). The dynamics of symptomatic and asymptomatic cases are highly influenced by the value of this parameter: if protection lasted only 50 days, this could result in a higher peak in February, with slightly more than 1,684,700 individuals with symptoms and 14,104,000 asymptomatic individuals (Figure S12A, B). The difference between the duration of 100 and 150 days is roughly 247,670 ($-22.97\%$) symptomatic and a 2,068,900 ($-2.29\%$) asymptomatic individuals (Figure S12A, B). The decay (at the end of March) followed by the peak in February has interesting dynamics for symptomatic, asymptomatic, and hospitalized individuals. A loss of protection in 50 days is results in a steady decay at the end of March 2023 with 468,870 symptomatic cases, 3,878,300 asymptomatic cases, and only 64,477 hospitalized individuals (Figure S12A–C). Instead, a loss of protection between 100 and 150 days will result in fewer than 279,410 active symptomatic cases at the end of March, 2,317,900 asymptomatic cases, and 47,165 hospitalized individuals (Figure S12A–C). Regardless of the duration of protection, all scenarios are associated with a stabilized behavior of infections at the end of March 2023. The death toll is affected by the re-infection parameter, where a 50-day period of protection yields a higher count with roughly 4458 more deaths compared with protection between 100 and 150 days (Figure S12D).

Now, we will assess the dynamics of cases between May and July 2023, the period where there may be a subsequent wave of Omicron. There is a difference of around two weeks in the time when the peak is reached when varying the loss of immunity: for 50 days, it will be reached in the second week of June; for 100 days, in the first two weeks of July 2023, and for 150 days, in the first two weeks of August 2023 (Figure 16). For a 50-day duration, symptomatic cases will peak at 985,000 and approximately 8,000,000 for asymptomatic cases (Figure 16A, B), and by the end of September 2023, there will be around 900,000 symptomatic and 5,000,000 asymptomatic cases (Figure 16A, B). In this graph, there is no decay regardless of the duration of the loss of immunity. For a 100-day duration, cases will peak at roughly 532,546 symptomatic and 5 million asymptomatic cases, followed by a slight descent to stabilize at approximately 500,000 symptomatic and 4.5 million asymptomatic active cases (Figure 16A, B). For a 150-day duration, cases will peak at 450,000 and 3.2 million for symptomatic and asymptomatic active cases, respectively. Once a peak is reached, there is a modest descent and a stabilization in 400,000 and 2.4 million symptomatic and asymptomatic cases, respectively (Figure 16A, B). There is a noticeable difference in the size of the subsequent wave of COVID-19 influenced by the waning of immunity. Despite a considerable increase of infected cases, hospitalized individuals will rise, but not by a considerable amount. For a 50-day duration in the peak, there will be roughly 620,000 individuals; for 100-days, 350,000 individuals; and for 150-days, 220,000 individuals (Figure 16C). This modest increase in hospitalizations is a product of a population with a level of protection from the virus that influences a re-infection, but the infection does not produce severe COVID-19 due to a combination of immunity acquired by the vaccine or by a prior infection. The death toll is affected as well, with a decay of around 100,000 deaths with 50 days of protection when we compared with 100-days, and a difference of 50,000 difference between the cases with 150 and 100 days (Figure 16D).

3.7.   Influence of the probability of re-infection prior to get infected with another Omicron sublineage on the pattern of new infections for the following months

Re-infection prior to get infected with another sublineage, or cross immunity, is defined as the immune response from one infection that will confer a protection from another infection. The following waves may be affected by the existence of cross immunity. The third wave of COVID-19 in Mexico was heavily influenced by the presence of XBB, which affected the peak of cases when BQ.1 was introduced in the population ^[24]. Since a prior infection with XBB could affect a sudden re-infection with BQ.1, an important portion of the population had protection and may be infected in this wave, but the story can be different in the subsequent wave. Between January and April 2023, the peak was influenced by the value of immunity prior to infection with XBB: there was a reduction of approximately 140,000 active symptomatic cases and 1,200,000 asymptomatic cases, because they had a prior infection with XBB (Figure S13A, B). Hospitalizations and death toll also are affected by this cross-immunity between sublineages with a reduction of 80,000 hospitalizations and 20,000 fewer deaths due to the influence of prior XBB infection (Figure S13C, D).

For the third Omicron wave in the last month of 2022, the size of the peak and the time it takes to reach that peak are affected by a prior infection of XBB and BQ.1. A lower protection is associated to a peak of around 18,000 new infections with symptoms and 200,000 asymptomatic cases (Figure 17A, B). With baseline protection, the peak is reduced by 100,000 infections, and it takes just about a one and a half months to reach that peak compared with a lower protection, when it only three weeks to reach it (Figure 17). For the case of asymptomatic infections, there may a be reduction of approximately 10,000 active cases that were infected in the prior wave (Figure 17B). Hospitalizations are affected as well: for lower protection, there may be roughly 20,000 new hospitalizations during the second wave of 2023; instead, baseline protection may be associated with a decrease of 2000 infections, reaching 18,000 new hospitalizations between May and July 2023 (Figure 17C). Finally, the death toll is influenced by the prior infection with another sublineage: for a low protection, there may be roughly 60,000 deaths due to COVID-19; alternatively, baseline protection can significantly decrease the death toll, with a reduction of 30,000 deaths (Figure 17D). The protection, or how well the new variant avoids elimination by the immune system, can influence the the dynamics of the following waves of COVID-19, unless another variant arises where a previous infection with another sublineage no longer results in some type of protection.

3.8.   Assessing the combined effect of face mask usage, face mask efficiency and efficiency of protection due to partial immunity in reducing the control reproduction number

In this section, we evaluate the importance of the use of face masks and the efficiency, combined with the efficiency of protection due to partial immunity, to reduce the control reproduction number. To take this into account, we rescale the transmission rates $\beta_1$, $\beta_2$, $\beta_3$, and $\beta_4$ by a factor $\left(1-C_m\epsilon_m\right)$, where $C_m$ denotes the proportion of the population that uses face masks and $\epsilon_m$ represents the efficiency of face masks.

Figure 18 depicts the values of the control reproduction number obtained using the baseline parameter values as before, for three different values of the protection due to partial immunity against sublineage BQ.1 ($\psi_2$). Since the reproduction number corresponding to sublineage XBB is lower than that of sublineage BQ.1, a variation in the value of the protection against the first sublineage ($\psi_1$) does not modify the value of the control reproduction number.

3.9.   Local sensitivity analysis of the control reproduction number

To evaluate the sensitivity of the control reproduction number as each of the model parameters is varied, we computed the elasticity indices of $\mathcal{R}_c$ with respect to each parameter. Figure 19 shows the corresponding values for the parameters that have an elasticity index greater than zero. The proportion of symptomatic cases for sublineage BQ.1 ($p_2$) and the transmission rate of the same sublineage by contact with symptomatic individuals ($\beta_3$) have the largest importance on increasing the reproduction number; in comparison, the transmission rate by contact with asymptomatic individuals ($\beta_4$) has a very small importance. On the other hand, the most important parameter associated with a decrease in $\mathcal{R}_c$ is the recovery rate of infection by sublineage BQ.1 ($\gamma_2$). The protection against this sublineage given by partial immunity (either natural or vaccine-acquired), denoted by $\psi_2$, has a smaller, although still significant importance on decreasing it, and it is followed by the face mask usage $C_m$ and face mask efficiency $\epsilon_m$. The next most important parameters are the hospitalization rate ($\sigma$) and the proportion of hospitalized cases for sublineage BQ.1 ($\varphi_2$). The rest of the parameters have only a very small influence on the value of the control reproduction number.

4.   Discussion

More than two years after the pandemic started, we are still using mathematical models to forecast the spread of SARS-CoV-2 and help inform public health decisions. In Mexico, the year 2022 was marked by the entry of Omicron into the population, mixed with the roll-out of boosters and the continuous use of face masks in open and closed places ^[44]. Despite the use of face masks in the population, the transmission was high when Omicron entered, which resulted in an exponential increase of newly infected individuals. As of January 2023, the population vaccinated against COVID-19 in Mexico was around 76% ^[45], meaning there is a sector of the population that was not vaccinated, with vaccine hesitancy being one of the most important reasons ^[46,47].

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 wave and subsequent waves, as well as the impact of vaccination campaigns in Mexico. Our compartmentalized model considered the spread in the unvaccinated and vaccinated subpopulations to better portray the ongoing COVID-19 pandemic. We also considered the particular vaccine campaigns in Mexico with seven different vaccines against COVID-19, as well as how this heterogeneity affected the global dynamics of the spread of COVID-19. Furthermore, we considered the transmission and immunity parameters associated with the XBB and BQ.1 variants. Our population was divided into four important subpopulations: unvaccinated (susceptible individuals), individuals vaccinated with only a two-dose scheme, individuals vaccinated with one or two boosters, and partially protected susceptible individuals (individuals with waning immunity, either acquired or natural). These subdivisions helped us understand which population will be the force for the subsequent Omicron waves if no other variant emerges and what public health decisions we need to face the following waves. Despite our model being constructed based on previous works and to our concern, there is no mathematical model applied in Mexico that introduces novelties such as vaccine heterogeneity, the inclusion of a hospitalized compartment, vaccine booster, the use of face masks and their related efficiencies, cross-immunity between Omicron sublineages, and modeling the waning of natural and acquired immunity as a two-stage process.

Face mask usage in the population and face mask efficiency are two important parameters with a large impact on the dynamics of the control reproduction number. We noted that N95 face masks are the most effective to reduce the spread and eliminate the pandemic, and this type of face mask is inexpensive based on the economy of this year ^[48]. The least effective face masks are cloth and improved cloth face masks, as expected. However, if 75% of the population used cloth face masks (30% efficiency), the new infections (symptomatic or not) can be reduced by 20%; while if 75% of the population uses an improved cloth face mask (50% efficiency), the transmission is reduced in 40% compared with no use of face masks. Surgical masks, which are a more suitable type of mask to purchase based on their cost, can reduce transmission by 64% compared with no use of face mask if 75% of the population wears them. If there are still people who do not want to wear face masks, with only 50% wearing surgical face masks, there is a 40% decrease in new cases, whether symptomatic or asymptomatic. Surgical face masks are the best option in terms of reduction in transmission and costs. Vaccines have reduced hospitalizations and therefore deaths ^[42], although they have not been able to reduce transmission between individuals. Our best strategy to reduce infections is the use of face masks ^[49,50]. We highly recommend that the use of face masks continues until the end of 2023 to reduce the spread, but this can change based on the dynamics of the spread in the following months, combined with the possible emergence of a new sublineage of Omicron or a variant.

The application of multiple vaccines in Mexico has significantly helped to reduce hospitalizations and deaths, even though they were developed to combat the original Wuhan strain. Even so, the cross-protection efficiency is waning over time for vaccine-acquired, natural, or hybrid immunity ^[51,52]. Our simulations showed an increase in new infections when we vary the immunity period between 50 days and 150 days, with the possibility of another Omicron BQ.1 wave between November and December 2022. The size of the following peak may be influenced by the loss of protection, either acquired or natural. If the loss of immunity were 50 days, there may be an increase of 50% in new infections. Furthermore, if we increase the loss of immunity to 150 days, the peak in November may be reduced by 10% compared with the baseline 100-day immunity period. The predicted peak will be 75% lower than the peak in July–August 2022. There is a possibility of another Omicron wave by May 2023, with a 40% reduction with respect to the wave in November–December 2022. The pattern of waning immunity will shape the behavior of following waves and how we will face the wave. Nevertheless, the application of a booster of bivalent vaccine that contains the ancestral Wuhan and the Omicron BA.1 in the fall may be a significant change in terms of loss of acquired immunity ^[53].

Cross-immunity is an important concept to be considered in a mathematical model that incorporates two sublineages of the same variant ^[54]. This means that if someone gets infected with a prior sublineage or variant, they may have protection because of the vaccine-acquired or natural immunity. Since XBB was the main driver for the second wave in Mexico, we incorporated how cross-immunity affected the BQ.1 wave in Mexico and how it may affect subsequent waves. We were able to observe that, as the cross-immunity varied, the size of the next wave changed. With greater protection, the peak decreases compared to the case with less protection, a reduction of almost 25%. In addition, if there was a higher protection, the time to reach the peak would increase. Less protection implies that it takes a month and a half to reach the peak, while with greater protection it takes two months and a week to reach the peak. We highly suggest monitoring cross-immunity protection on a long-term scale to visualize possible subsequent waves of Omicron sublineages or multiple rolling strains.

5.   Conclusions

In essence, our mathematical model provides a projection of the ongoing COVID-19 pandemic in Mexico, distinguishing new cases arising by the multiple vaccines applied in the vaccination campaign. We demonstrated that all vaccines applied, regardless of their vaccine efficiency against symptomatic COVID-19, have helped in reducing hospitalizations and deaths from COVID-19. We understand that this type of mathematical model is heavily influenced by social behavior, but it still can inform us to take important public health decisions. We recommend the continuous use of face masks in open and closed places with at least 50% efficiency protection against Omicron sublineage, at the minimum for the rest of the year. Finally, two important parameters will shape the behavior of the pandemic: waning immunity and cross-immunity between Omicron sublineages or with a possible emerging strain. We recommend monitoring these parameters and the way they affect the control reproduction number and herd immunity, as well as methods for controlling them using pharmaceutical or non-pharmaceutical strategies.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This article was supported in part by Mexican SNI under CVU 15284.

Data and materials availability

The original contributions presented in this study are included in the article–Supplementary material. All code used for analysis is available at https://github.com/UgoAvila/Mathematical-Model-of-COVID-19-in-Mexico-waning-and-cross-immunity.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A1.   Mathematical model with asymptomatic individuals, vaccination, waning immunity, and cross-immunity protection between sublineages

The model studied here is described by the following system of differential equations:

where

and the variables have the following meaning:

● $S$: unvaccinated susceptible individuals.

● $F$: individuals vaccinated with at least one dose but no booster doses.

● $B$: individuals vaccinated with at least one booster dose.

● $E_1$: individuals exposed to sublineage 1 Omicron.

● $I_1$: symptomatic individuals infected with sublineage 1 Omicron.

● $A_1$: asymptomatic individuals infected with sublineage 1 Omicron.

● $H_1$: hospitalized with sublineage 1 Omicron.

● $R_1$: recovered after infection with sublineage 1 Omicron.

● $S_p$: susceptible individuals with partial immunity.

● $E_2$: exposed individuals to sublineage 2 Omicron.

● $I_2$: symptomatic individuals infected with sublineage 2 Omicron.

● $A_2$: asymptomatic individuals infected with sublineage 2 Omicron.

● $H_2$: hospitalized with sublineage 2 Omicron.

● $R_2$: recovered after infection with sublineage 2 Omicron.

● $D_1$: cumulative number of deaths due to sublineage 1 Omicron.

● $D_2$: cumulative number of deaths due to sublineage 2 Omicron.

A2.   Mathematical analysis

A2.1.   Boundedness of the system

It is readily verified that every solution of system (A1.1) with non-negative initial conditions remains non-negative for all $t\ge0$. Let $X = S + F + B + E_1 + I_1 + A_1 + H_1 + R_1 + S_p + E_2 + I_2 + A_2 + H_2 + R_2$. Then,

which implies

Since $-(\omega_1\chi_1H_1+\omega_2\chi_2H_2)\leq 0$,

Using the integrating factor

we get

If we let $t \to \infty$, we obtain

hence, all the solutions of (A1.1) initiating in $\left\{ \mathbb{R}^{10}_+ \setminus \{0 \} \right\}$ are confined to the region

for any $\epsilon > 0$ and for $t \to \infty$.

A2.2.   Equilibria

We will proceed to find the equilibria of system (A1.1). First, let us notice that $D_1$ and $D_2$ are not in the other compartments, so we can omit them from the system. Now, we will equalize to zero $\frac{dS}{dt}, \frac{dF}{dt}, \frac{d B}{dt}, \frac{dE_1}{dt}, \frac{dI_1}{dt}, \frac{dA_1}{dt}, \frac{dH_1}{dt}, \frac{dR_1}{dt}, \frac{dS_p}{dt}, \frac{dE_2}{dt}, \frac{dI_2}{dt}, \frac{d A_2}{dt}, \frac{dH_2}{dt}, \frac{dR_2}{dt}$, obtaining

If we add all the equations, we obtain

and thus

A2.2.1.   The disease-free equilibrium

If $I_1 = I_2 = 0$, it follows from Eq (A2.5) that $p_1wE_1 = 0$, so $E_1 = 0$. In consequence, by (A2.6),

Then, $A_1 = 0$. It follows from (A2.7) that

Then, $H_1 = 0$. From (A2.8),

Thus, $R_1 = 0$. A similar argument, using $I_2 = 0$ with Eqs (A2.11)–(A2.14) tells us that $E_2 = A_2 = H_2 = R_2 = 0$.

In consequence, the system is reduced to

By (A2.16),

Substituting in (A2.17) and simplifying, we get

Using (A2.20) in (A2.18), we have

Last, substituting (A2.20) and (A2.21) in (A2.19) yields

Therefore, the disease-free equilibrium $E_0$ takes the form

where

A2.2.2.   Endemic equilibrium

From (A2.5), we have

Substituting (A2.22) in (A2.6), we get

From (A2.7),

Substituting (A2.22)–(A2.24) in (A2.8), we obtain

where

Analogously, we obtain

and

where

Substituting (A2.24) and (A2.26) in (A2.15), we obtain

where $D_0 = (\omega_1+\mu)(\omega_2+\mu)$, $D_2 = \omega_1\chi_1\varphi_1(\omega_2+\mu)$, and $D_3 = \omega_2\chi_2\varphi_2(\omega_1+\mu)$. Notice that, for the equilibrium to be positive, we require that $\mathcal{G}(I_1, I_2) > 0$, where $\mathcal{G}$ is defined by

From (A2.1),

Substituting (A2.23), (A2.25), and (A2.27), we obtain

where

From (A2.2),

Substituting (A2.23), (A2.25), (A2.27), and (A2.29) gives us

From (A2.3),

Substituting (A2.23), (A2.25), and (A2.27) gives us

From (A2.9), (A2.23), (A2.25), and (A2.27), we get

Using the expressions for $B$, $F$, $R_1$, and $R_2$, we can write the above equations in terms of $I_1$ and $I_2$. Finally, using (A2.4) and (A2.10), we can obtain a system of two quadratic equations in $I_1$ and $I_2$. However, since the resulting equations become too unwieldy, we will analyze them only in the case when the parameter $\rho$ is zero.

A2.2.3.   Endemic equilibrium without vaccination

In the case where there is no vaccination except for booster doses, that is, when $\rho = 0$, we obtain the following equations for the endemic equilibrium (the process is analogous to Section A2.2.2):

where

Furthermore, if $I_1, I_2 > 0$, we can see that the condition $\mathcal{G}(I_1, I_2) > 0$ becomes necessary and sufficient for all other variables to be nonnegative.

Substituting back in equations (A2.4) and (A2.10), we see that $I_1$ and $I_2$ satisfy

where

Therefore, we conclude the following result:

Theorem A2.1. When $\rho = 0$, system (A1.1) has an endemic equilibrium if and only if there exist $I_1, I_2 > 0$ satisfying (A2.33) such that

A2.3.   Basic reproduction number

The infected compartments of the population are $E_1, I_1, A_1, H_1, E_2, I_2, A_2$, and $H_2$. Let $m_1 = w+\mu$, $m_2 = \gamma_1(1-\varphi_1) + \sigma\varphi_1 + \mu$, $m_3 = \gamma_1+\mu$, $m_4 = \omega_1+\mu$, $m_5 = \gamma_2(1-\varphi_2) + \sigma\varphi_2 + \mu$, $m_6 = \gamma_2+\mu$, and $m_7 = \omega_2+\mu$. Using the notation in ^[39], it follows that the new infections and transition matrices are given by

Then,

where

Calculating the inverse of $V$, we get

Multiplying $F$ and $V^{-1}$, we get

which has the eigenvalues $\lambda_1 = \dfrac{\mathcal{M}_1wp_1m_3+\mathcal{M}_2w(1-p_1)m_2}{m_1 m_2 m_3} > 0$, $\lambda_2 = \dfrac{\mathcal{M}_3wp_2m_6 + \mathcal{M}_4(1-p_2)wm_5}{m_1 m_5 m_6} > 0$, and $\lambda_i = 0$, $i \in \{3, 4, 5, 6, 7, 8\}$.

Therefore, the basic reproduction number with vaccination, also known as the control reproduction number, is given by

where

In the particular case when $\rho = 0$, we can obtain the basic reproduction number with no vaccination, $\mathcal{R}_0$, which is given by the expression

Furthermore, by [39, Theorem 2], we can conclude the following result, which ensures the threshold dynamics of the disease-free equilibrium.

Theorem A2.2. The disease-free equilibrium $E_0$ of system (A1.1) is locally asymptotically stable if $\mathcal{R}_c < 1$ and unstable if $\mathcal{R}_c > 1$.