Vaccination and variants: A COVID-19 multi-strain model evolution for the Philippines

Coronavirus disease 2019 (COVID-19) management and response is a challenging task due to the uncertainty and complexity of the nature surrounding the virus. In particular, the emergence of new variants and the polarizing response from the populace complicate government efforts to control the pandemic. In this study, we developed a compartmental model that includes (1) a vaccinated compartment, (2) reinfection after a particular time, and (3) COVID-19 variants dominant in the Philippines. Furthermore, we incorporated stochastic terms to capture uncertainty brought about by the further evolution of the new variants and changing control measures via parametric perturbation. Results show the importance of booster shots that increase the vaccine-induced immunity duration. Without booster shots, simulations showed that the dominant strain would still cause significant infection until 31 December 2023. Moreover, our stochastic model output showed significant variability in this case, implying greater uncertainty with future predictions. All these adverse effects, fortunately, can be effectively countered by increasing the vaccine-induced immunity duration that can be done through booster shots.


. Introduction
Since the emergence of coronavirus disease 2019  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 agentbased 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.
. Model formulation . . Deterministic model In this study, we divided the population into eight compartments, namely, susceptible (S), vaccinated (V), infected by the original strain (I 1 ), infected by the Delta variant (I 2 ), infected by the Omicron variant (I 3 ), confirmed (C), recovered (R), and dead (D). Figure 1 shows the dynamics of the model. The model considers the following facts. First, vaccines do not provide lasting immunity. Hence, vaccinated people will be back to being susceptible after some time. Second, reported confirmed cases do not show the full extent of the infection. Third, estimates of unreported cases are not being accounted for. Finally, recovered individuals (from natural infection) also do not have lasting immunity against the virus; hence, both vaccinated and recovered individuals may return to the susceptible population after some time. In the model, we assume density-dependent transmission rates.
An infected individual by any of the variants can infect people in compartment S. We assume that vaccinated people in V are immune to any of the variants. The parameter β i denotes the transmission rate of the disease caused by the ith strain. The parameter α represents the level of control measures implemented to limit the transmission. The parameter ν denotes the vaccination rate of susceptible people. The parameters λ 1 and λ 2 denote the vaccine immunity waning rate and natural immunity waning rate, respectively. The parameter ǫ denotes the proportion of infections detected and confirmed through testing also known as detection rate. We also have the removal rate from C given by δ. Finally, the recovery rate from unconfirmed infections caused by the ith strain is denoted as γ i and the probability to recover from the infection is denoted as ρ. The parameters are summarized in Table 1.
In this study, we adapted a closed population model, i.e., we did not consider natural birth and death rates. We denote by N 0 the total number of population at the start. The dynamics of our model is governed by the following system of ordinary differential equations (ODE): Frontiers in Applied Mathematics and Statistics frontiersin.org . /fams. . .

. Stochastic model
As will be seen in Section 5, the third strain will become dominant as time goes on. We acknowledge that we need to accommodate future uncertainties in our simulation. Hence, we modify our ODE model to add stochasticity. We apply parametric perturbation to the reduced transmission rate αβ 3 . The resulting system with stochastic differential equations (SDE) is as follows: where dB/dt is the white noise, i.e., the derivative of the standard Brownian motion B(t), and σ > 0 denotes the intensity of that noise.
. Qualitative analysis of the ODE model Since we have a close system and N 0 = S + V + I 1 + I 2 + I 3 + 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 1 , I 2 , I 3 , and C equal to 0. We obtain a solution, which is the DFE

. . Reproduction number
By definition, the basic reproduction number R 0 denotes the average number of individuals directly infected by a single infected individual over the duration of its infectious period in a population without any deliberate intervention to stop its spread. We will compute R 0 using the next generation matrix method defined by Diekman et al. [14] and van den Driessche and Watmough [15]. Let X be the vector of the infected classes and Y be the vector of the other classes. Let F (X, Y) be the vector of new infection rates (flows from Y to X) and let V(X, Y) be the vector of all other rates (not new infections). Then, we have Evaluating the derivatives of F and V at the DFE, we are led to the following matrices Hence, the next generation matrix is given by The eigenvalues of K are the following: , and η 4 = 0. The eigenvalue η 1 is associated with strain 1 and gives rise to the basic reproduction number . Similarly, eigenvalue η 2 associated with strain 2 corresponds to the basic , and eigenvalue η 3 associated with strain 3 corresponds to R 3 = λ 1 αβ 3 ν(γ 3 + ǫ) Finally, we take R 0 = max{R 1 , R 2 , R 3 }.
. . Stability analysis of the DFE Proof. Consider that and This is finite, since D is bounded. Hence, δ(1 − ρ) C(t) → 0 as t → +∞. Since δ(1 − ρ) > 0, we have C(t) → 0 as t → +∞. Using the same argument and noting that we also have I 1 (t), I 2 (t), I 3 (t) → 0 as t → +∞. Similar deduction can be used to show that R(t) → 0 as t → +∞, using Equation (7). 2 We may define the effective reproduction number associated with strain i by Compared to the basic reproduction number R 0 , the effective reproduction number R e i (t) denotes the average number of new infections associated with strain i, at time t, caused by a single infected individual, considering that in the population at this time, there are already some individuals who are no longer susceptible.

. Existence of solution for the SDE Model
Let ( , F , {F t } t≥0 , P) be a complete probability space with filtration {F t } t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F 0 contains all P-null sets). Let R 8 + = {x i > 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 x 0 ∈ R 8 + , there is a unique solution x(t) of system (9) on t ≥ 0, and the solution will remain in R 8 + with probability 1, namely, x(t) ∈ R 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 x 0 ∈ R 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 0 ≥ 1 be sufficiently large so that all components of x 0 are contained in the interval [1/s 0 , s 0 ]. For each integer s ≥ s 0 , we define the stopping time by τ s = inf{t ∈ [0, τ e ) : at least one of S, V, I 1 , I 2 , I 3 , C, R, or D ∈ (1/s, s)}.
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.2

. Numerical simulations
Data on confirmed cases used in parameter calibration are from the COVID-19 Data Repository by the Center for Systems Science Frontiers in Applied Mathematics and Statistics frontiersin.org . /fams. .

FIGURE
Output of the model with calibrated parameters compared with data. The gray curves represent the other outputs of the approximate Bayesian computation approach, while the red curve represents the best fit model output (i.e., using the parameters in Table ). The green dots are the data from the COVID-data repository of JHU CSSE. 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 β 1 and α, and the initial conditions I 1 (0), I 2 (0), and I 3 (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.

FIGURE
Projection on the e ects of di erent levels of vaccine-induced immunity duration. The green vertical line is at Day , the end of the data used. One may interpret the case when λ = / to be the case when individuals only opt to be fully vaccinated. The cases λ = / and λ = / may correspond to the cases when individuals are also getting first and second booster shots, respectively.

. . Deterministic simulations on waning vaccine-induced immunity
The Philippines started its vaccination campaign last 1 March 2021, with the two-dose Sinovac vaccine [26]. We estimated that it would need 4 weeks or 28 days to be fully protected from the vaccine [11], so we started the parameters ν and λ 1 by Day 87 (29 March 2021). It is estimated that the vaccine-induced immunity wanes after 6 months and that booster shots are recommended after that time interval [27]. Some people are unwilling to take booster shots due to vaccine hesitancy. Accounting for this social behavior in the model, we consider three different vaccine-induced immunity duration through the parameter λ 1 by setting it to 1/180 (6 months). This accounts for the case when the population only takes full vaccination but no booster shots. On the contrary, we set λ 1 = 1/365 (12 months) when the people finished taking the first booster shot, while λ 1 = 1/540 (18 months) when they completed the second booster shot. We simulate up to 31 December 2023 (Day 1094). The result is shown in Figure 4.

. . Stochastic simulations on the dominant transmission rate
As shown in Figure 4, our simulation suggested that the original and Delta variants are to die out even with the minimum vaccineinduced immunity duration of 6 months. But it is not the case with regard to the Omicron variant. However, we acknowledge that much uncertainty can affect the reduced transmission rate αβ 3 -for instance, the changing control measures of the government and the further evolution of the variant. Hence, we also want simulations incorporating noise on αβ 3 . In Figure 5, we show the result of our stochastic simulations. The curves are for the confirmed compartment. The figure can be viewed as a 3 × 2 matrix where the rows correspond to the cases corresponding to the different vaccine-induced immunity duration, given by the value of λ 1 . The columns show the results concerning different noise intensities, given by the value of σ . The stochastic simulations start after Day 543 and run up to 31 December 2023 (Day 1094). The numerical simulations are implemented using the Euler-Maruyama scheme [28].

. Discussion
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 vaccineinduced 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  we only have R e 3 = 0.9998. However, with a vaccine-induced immunity duration of 12 months and at Day 543, we already have R e 3 = 0.5124. We acknowledge that a lot of changes can happen in the future. For instance, the ongoing evolution of the virus and the everchanging 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.

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