Since the emergence of coronavirus disease 2019 (COVID-19) in Wuhan, its causative agent, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has been rapidly evolving into new variants [1]. This emergence is expected as, like any virus, SARS-CoV-2 continues to mutate from time to time. So far, the World Health Organization (WHO) has identified five variants of concern, namely, Alpha, Beta, Gamma, Delta, and Omicron. These are the variants considered the most transmissible and dominant that are circulating the world [2]. While vaccine has been available, the campaign has suffered a series of setbacks and there are many issues leading to vaccine hesitancy [3, 4]. In the Philippines, the slow vaccine rollout, limited testing capacity, weak genomics surveillance, a fragile healthcare system, and a large informal economy contributed more to the many issues related to COVID-19. The latter is a social cost as few people can afford not to work but are forced to stay home in overcrowded housing. These factors exacerbated the situation as they gave the virus a perfect environment to mutate, feeding through a continuous supply of susceptibles.

Hence, a question of much importance is how to reduce COVID-19 transmission. In doing so, an approach of balancing the efforts on the vaccination campaign (booster shots amidst reinfection because of mutating variants and waning immunity), and preparedness of the healthcare system coupled with public health and social measures is a must. While the government is trying to strike a balance, the race between vaccination and fast mutating variants remains the biggest challenge.

Many computational modeling studies have investigated COVID-19 dynamics in the Philippines. For instance, in Arcede et al. [5], an SEIR-type model is constructed whose infected can either be symptomatic or not. The result shows that treating symptomatic alone does not reduce the spread. However, managing the number of susceptible does, which containment and vaccines have a significant impact role to play. Later, the same model was used to investigate the implemented non-pharmaceutical interventions (NPI) in the country [6]. Here, NPIs include lockdown, social distancing, mass testing, and strengthening the healthcare system. The study provided a choice for the government to implement the control by indicating economic cost (low, high) given no vaccine availability. In Bock et al. [7], an agent-based SIR model was used to investigate the prevalence of COVID-19 in two neighboring cities in Northern Mindanao, particularly in Iligan and Cagayan de Oro. The result shows that social distancing and age-specific quarantine can effectively slow down contagion. Furthermore, social distancing combined with an effective testing strategy can keep the epidemic at bay and prevent it from becoming a critical epidemic. In Arcede et al. [8], a regional COVID-19 model has been constructed in the cities mentioned earlier and in the Northern Mindanao region as a whole. The model is tailored to fit early transmission; hence appropriate models are suggested when laboratory-based disease reports are available. In [9], Mammeri et al. extended their SEIR-type model to account the spatial movement of individuals. Given five main islands and five main airports in the Philippines as nodes with index case assumed to start from Manila airport on day 1, their simulation show remarkably close similarity to what happened in the Philippines during its first 140 days. Studies mentioned do not deal with vaccination control. However, recent articles deal with vaccination strategies in the Philippine context. For instance, in [10], optimal control was used to investigate existing policy interventions, including vaccination rollouts, community quarantines, and simulated virus outbreaks. They found that early and effective implementation of precautionary measures such as community quarantines are crucial for containing outbreaks. They also found that even if vaccinations do not suffice, expanding the vaccine supply reduces the need for more resource-intensive interventions. Moreover, in Caga-anan et al. [11], a model with a delay on the vaccination compartment was constructed to study the impact of vaccination efforts on disease progression and herd immunity. The result shows that timely vaccination is preferred to maximize impact. They also assessed the performance of different vaccine brands in the model, showing Pfizer-BioNTech with the best results. Finally, some models were proposed for allocating resources. For instance, optimizing the location of vaccination sites implemented in San Juan Philippines [12] and distribution of COVID-19 testing kits in DOH-accredited testing centers in the country [13]. All studies mentioned earlier do not account explicitly for COVID-19 variants and randomness.

In this study, we evaluate the impact of vaccination and its waning induced immunity, as well as the uncertainties related to future mitigation policies and the further evolution of the virus, given the three variants that are circulating dominantly in the country.

Since we have a close system and N₀ = S + V + I₁ + I₂ + I₃ + C + R + D, we may just consider the system (1)–(7), i.e., without D. To find the disease-free equilibrium (DFE) (a steady-state solution of an epidemic model with all infected variables equals to zero), we equate Equations (1)–(7) to zero and the infective compartments I₁, I₂, I₃, and C equal to 0. We obtain a solution, which is the DFE given by (S^*, V^*, 0, 0, 0, 0, 0), where S*=λ1V*ν.

Let (Ω,F,{Ft}t≥0,ℙ) be a complete probability space with filtration {Ft}t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F0 contains all ℙ-null sets). Let ℝ+8={xi>0:i=1,2,...,8}. Let B(t) be a Brownian motion defined on the complete probability space Ω. Then, we have the following theorem showing that the stochastic system (9) has a unique non-negative global solution.

Theorem 4.1. For any given initial value x0∈ℝ+8, there is a unique solution x(t) of system (9) on t ≥ 0, and the solution will remain in ℝ+8 with probability 1, namely, x(t)∈ℝ+8 for all t ≥ 0 almost surely.

Proof. Since the coefficients of system (9) satisfy the local Lipschitz condition, it implies that for any given initial value x0∈ℝ+8, there is a unique local solution x(t) for every t ∈ [0, τ_e), where τ_e is the explosion time. To prove that the solution is global, we need to show that τ_e = ∞. To do so, we let s₀ ≥ 1 be sufficiently large so that all components of x₀ are contained in the interval [1/s₀, s₀]. For each integer s ≥ s₀, we define the stopping time by

Clearly, τ_s is increasing as s → ∞. Let τ∞=lims→∞τs, then τ_e ≥ τ_∞ almost surely. If we can show that τ_∞ = ∞ a.s., then τ_e = ∞ and x(t)∈ℝ+8 a.s. for all t ≥ 0.

Suppose τ_∞ < ∞, then there exists T > 0 such that ℙ{τ_∞ ≤ T} > ϵ for all ϵ ∈ (0, 1). Hence, there is an integer s₁ ≥ s₀ such that

Let a C2- function U:ℝ+8→ℝ+8 be defined by

Using Itô formula on Equation (15), we have

where

Note that K is a positive constant independent of the variables S, V, I₁,I₂, I₃, C, R, and D, and time t.

Thus,

Therefore, if t₁ ≤ T,

where τ_s ∧ t₁ = min{τ_s, t₁}. Taking expectations to both sides of (17), we obtain

By properties of Itô integral, we have

where U₀ = U(x₀). By Gronwall's inequality,

Let Ω_s = {τ_s ≤ T} for any s ≥ s₁. Then, by (14), we have ℙ(Ω_s) ≥ ϵ. Note that for every ω ∈ Ω_s, there is at least one of S(τ_s, ω), V(τ_s, ω), I₁(τ_s, ω), I₂(τ_s, ω), I₃(τ_s, ω), C(τ_s, ω), R(τ_s, ω), and D(τ_s, ω) that is equal to either s or 1/s. Consequently, U(x(τ_s, ω)) is no less than either

Thus,

It then follows from Equation (14) and Equation (21) that

where 1_Ωs is an indicator function of Ω_s. Letting s → ∞, then we have

which yields the contradiction. Therefore, we must have τ_∞ = ∞, almost surely.□

Data on confirmed cases used in parameter calibration are from the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU) [16]. The data are publicly available and so ethical approval is not required. The data are from 1 January 2021 (Day 0) to 28 June 2022 (Day 543). Values of some parameters and initial conditions are taken or estimated from sources as indicated in Tables 2, 3. We note that the parameter α varies over time as controls implemented by the government also vary. Hence, we considered α as a piecewise constant function, as shown in Table 2, where the dates correspond to the noticeable changes in the control measures implemented by the government.

To obtain the fitted values for the parameters β₁ and α, and the initial conditions I₁(0), I₂(0), and I₃(0), we minimized a non-linear least square function given by the sum of the square of the difference of the data and the model output. The optimization problem is solved using the approximate Bayesian computation approach combined with the Levenberg–Marquardt algorithm [23–25]. Visualization of the optimization result is given in Figure 2. In Figure 3, we plotted the evolution of the three variants based on the optimization result.

With our deterministic model, we have obtained parameters fitting observed data on confirmed cases in the Philippines. Our model has the added value that it could estimate the progression of the three main variants circulating the Philippines. Our simulations showed that the Omicron variant would be the dominant variant as we advance and the other variants die out (Figure 4). Moreover, it showed the significantly faster transmission of the Omicron variant compared with the other two, as reported [29] and validated by some studies [30, 31]. Our population is a closed model, and Theorem 3.1 showed that the Omicron variant would eventually die out. However, the time of its realization depends on the vaccine-induced immunity duration. With a duration of approximately 6 months, the simulation shows a noticeable presence of the Omicron variant by 31 December 2023. However, this presence is reduced significantly with vaccine-induced immunity durations of 12–18 months. We note that these vaccine-induced immunity durations can be achieved with booster shots in the population. Looking at the effective reproduction number given in (13), by Day 543 and vaccine-induced immunity duration of 6 months, we have R1e=0.1731, R2e=0.3463, and R3e=1.0390. By Day 5000, we only have R3e=0.9998. However, with a vaccine-induced immunity duration of 12 months and at Day 543, we already have R3e=0.5124.

We acknowledge that a lot of changes can happen in the future. For instance, the ongoing evolution of the virus and the ever-changing control measures being implemented by the government. Hence, incorporating these uncertainties in our simulations can prove beneficial. Our stochastic simulations (Figure 5) reveal that uncertainties, represented as noise in our model, can significantly affect future outcomes when vaccine-induced immunity duration is only 6 months, as seen in the spread of the stochastic model output. An increase in vaccine-induced immunity duration means a decrease in the variability of our model output and hence better projection. The government can then use better projection to design more robust and refined intervention strategies to control the virus effectively.

Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/CSSEGISandData/COVID-19.

RC-a and JA conceptualized the paper and drafted the manuscript. HC, MR, and RC-a did the mathematical analysis. JM did the biological part. RC-a performed the simulations. All authors contributed to manuscript revision, read, and approved the submitted version.