<figure><img src="image_1.png" /></figure>

HAL
open science

# Effect of Vaccination to COVID-19 Disease Progression and Herd Immunity

Randy Caga-Anan, Michelle Raza, Grace Shelda G. Labrador, Ephrime Metillo, Pierre Del Castillo, Youcef Mammeri

► To cite this version:

Randy Caga-Anan, Michelle Raza, Grace Shelda G. Labrador, Ephrime Metillo, Pierre Del Castillo, et al.. Effect of Vaccination to COVID-19 Disease Progression and Herd Immunity. Computational and Mathematical Biophysics, 2021, 9 (1), pp.262-272. 10.1515/cmb-2020-0127. hal-03636448

HAL Id: hal-03636448

https://hal.science/hal-03636448

Submitted on 28 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effect of Vaccination to COVID-19 Disease Progression and Herd Immunity

Randy Caga-Anan, Michelle Raza, Grace Shelda G. Labrador, Ephrime Metillo, Pierre Del Castillo, Youcef Mammeri

► To cite this version:

Randy Caga-Anan, Michelle Raza, Grace Shelda G. Labrador, Ephrime Metillo, Pierre Del Castillo, et al.. Effect of Vaccination to COVID-19 Disease Progression and Herd Immunity. Computational and Mathematical Biophysics, 2021, 9 (1), pp.262-272. 10.1515/cmb-2020-0127. hal-03636448

HAL Id: hal-03636448

https://hal.science/hal-03636448v1

Submitted on 28 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Research Article

Open Access

Randy L. Caga-anan*, Michelle N. Raza, Grace Shelda G. Labrador, Ephrime B. Metillo, Pierre del Castillo, and Youcef Mammeri

# Effect of Vaccination to COVID-19 Disease Progression and Herd Immunity

https://doi.org/10.1515/cmb-2020-0127

Received July 17, 2021; accepted December 3, 2021

**Abstract:** A mathematical model of COVID-19 with a delay-term for the vaccinated compartment is developed. It has parameters accounting for vaccine-induced immunity delay, vaccine effectiveness, vaccination rate, and vaccine-induced immunity duration. The model parameters before vaccination are calibrated with the Philippines' confirmed cases. Simulations show that vaccination has a significant effect in reducing future infections, with the vaccination rate being the dominant determining factor of the level of reduction. Moreover, depending on the vaccination rate and the vaccine-induced immunity duration, the system could reach a disease-free state but could not attain herd immunity. Simulations are also done to compare the effects of the various available vaccines. Results show that Pfizer-BioNTech has the most promising effect while Sinovac has the worst result relative to the others.

**Keywords:** COVID-19; vaccination; delay differential equation

**MSC:** 92D30; 37N25; 34D20

## 1 Introduction

The COVID-19 pandemic, which was reported to have originated from China [18], is a rapidly evolving public health problem that immobilized nearly the entire world, causing many social disruptions, and brought havoc on many nations' economies. Considering the health threat brought by the disease, affected countries implemented control measures and strategies to minimize and prevent the spread of the virus. Before the availability of vaccines, countries worldwide practiced and applied non-pharmaceutical interventions (NPIs), e.g., wearing of face masks and face shields, lockdowns, social distancing, contact tracings, quarantines, isolation, and other public health measures and strategies. Since the start of the pandemic, various dilemmas affected the progression estimation of the disease and designing effective interventions. Various mathematical models have been created to understand and estimate the effect of these conjectural factors, and to evaluate and design intervention strategies. For instance, at the early part of the pandemic, the role of the asymptomatic infected are not yet clearly studied and so mathematical models are used to examine its likely impact on disease progression and intervention. Among the very large literature, we can cite [1, 3, 8, 14]. Models were also used to find optimal policies to mitigate the disease, e.g. [4, 10] dealt with the Philippines.

<figure><img src="image_4.png" /></figure>

<footnote>*Corresponding Author: Randy L. Caga-anan: Mindanao State University-Iligan Institute of Technology, Iligan City, Philippines</footnote>

<footnote>Mathematical Biology Research Cluster, Complex Systems Group, PRISM, MSU-IIT, E-mail: randy.caga-anan@g.msuiit.edu.ph</footnote>

<footnote>Michelle N. Raza: Institute of Mathematics, University of the Philippines-Diliman, Quezon City, Philippines</footnote>

<footnote>Grace Shelda G. Labrador: Mindanao State University-Iligan Institute of Technology, Iligan City, Philippines</footnote>

<footnote>Mathematical Biology Research Cluster, Complex Systems Group, PRISM, MSU-IIT</footnote>

<footnote>Ephrime B. Metillo: Mindanao State University-Iligan Institute of Technology, Iligan City, Philippines</footnote>

<footnote>Pierre del Castillo, Youcef Mammeri: Laboratoire Amiénois de Mathématique Fondamentale et Appliquée, CNRS UMR 7352, Université de Picardie Jules Verne, 80069 Amiens, France</footnote>
<figure><img src="image_2.png" /></figure>

Research Article

Open Access

Randy L. Caga-anan*, Michelle N. Raza, Grace Shelda G. Labrador, Ephrime B. Metillo, Pierre del Castillo, and Youcef Mammeri

# Effect of Vaccination to COVID-19 Disease Progression and Herd Immunity

https://doi.org/10.1515/cmb-2020-0127

Received July 17, 2021; accepted December 3, 2021

**Abstract:** A mathematical model of COVID-19 with a delay-term for the vaccinated compartment is developed. It has parameters accounting for vaccine-induced immunity delay, vaccine effectiveness, vaccination rate, and vaccine-induced immunity duration. The model parameters before vaccination are calibrated with the Philippines' confirmed cases. Simulations show that vaccination has a significant effect in reducing future infections, with the vaccination rate being the dominant determining factor of the level of reduction. Moreover, depending on the vaccination rate and the vaccine-induced immunity duration, the system could reach a disease-free state but could not attain herd immunity. Simulations are also done to compare the effects of the various available vaccines. Results show that Pfizer-BioNTech has the most promising effect while Sinovac has the worst result relative to the others.

**Keywords:** COVID-19; vaccination; delay differential equation

**MSC:** 92D30; 37N25; 34D20

## 1 Introduction

The COVID-19 pandemic, which was reported to have originated from China [18], is a rapidly evolving public health problem that immobilized nearly the entire world, causing many social disruptions, and brought havoc on many nations' economies. Considering the health threat brought by the disease, affected countries implemented control measures and strategies to minimize and prevent the spread of the virus. Before the availability of vaccines, countries worldwide practiced and applied non-pharmaceutical interventions (NPIs), e.g., wearing of face masks and face shields, lockdowns, social distancing, contact tracings, quarantines, isolation, and other public health measures and strategies. Since the start of the pandemic, various dilemmas affected the progression estimation of the disease and designing effective interventions. Various mathematical models have been created to understand and estimate the effect of these conjectural factors, and to evaluate and design intervention strategies. For instance, at the early part of the pandemic, the role of the asymptomatic infected are not yet clearly studied and so mathematical models are used to examine its likely impact on disease progression and intervention. Among the very large literature, we can cite [2, 3, 15, 11, 17]. Models were also used to find optimal policies to mitigate the disease, e.g. [4, 10] dealt with the Philippines.

<figure><img src="image_3.png" /></figure>

<footnote>*Corresponding Author: Randy L. Caga-anan: Mindanao State University-Iligan Institute of Technology, Iligan City, Philippines</footnote>

<footnote>Mathematical Biology Research Cluster, Complex Systems Group, PRISM, MSU-IIT, E-mail: randy.caga-anan@g.msuiit.edu.ph</footnote>

<footnote>Michelle N. Raza: Institute of Mathematics, University of the Philippines-Diliman, Quezon City, Philippines</footnote>

<footnote>Grace Shelda G. Labrador: Mindanao State University-Iligan Institute of Technology, Iligan City, Philippines</footnote>

<footnote>Mathematical Biology Research Cluster, Complex Systems Group, PRISM, MSU-IIT</footnote>

<footnote>Ephrime B. Metillo: Mindanao State University-Iligan Institute of Technology, Iligan City, Philippines</footnote>

<footnote>Pierre del Castillo, Youcef Mammeri: Laboratoire Amiénois de Mathématique Fondamentale et Appliquée, CNRS UMR 7352, Université de Picardie Jules Verne, 80069 Amiens, France</footnote>
With the availability of vaccines, other uncertainties and questions emerged and various mathematical models have been developed and studied to answer some of them. In Jentsch et al. [7], a coupled social-epidemiological model of SARS-CoV-2 transmission to address the problem of vaccine prioritization was developed. The authors provided answers whether to prioritize vaccinating individuals who cause the most transmission or those who are at the highest risk of death. One of their main results is that if vaccines become available sufficiently late in the pandemic, it is relevant to prioritize vaccinating individuals who caused the most transmission than those at higher risk of death (like those aged 60 years and older). In Moore et al. [14], a mathematical model structured by age and the UK region was used to conduct an analysis of a vaccination programme together with the relaxation of NPIs. Their modelling study estimated that vaccination alone is insufficient to contain the outbreak. In [1], Acuña-Zegarra et al. formulated an optimal control problem with mixed constraints to study optimal vaccination policies. Their solution identifies vaccination policies that minimize the burden of COVID-19 quantified by the number of disability-adjusted years of life lost.

In this study, we develop a mathematical model of COVID-19 that is used to investigate the now yet unclear factors related to vaccination. These factors are vaccine-induced immunity delay, vaccination rate, vaccine effectiveness, and vaccine-induced immunity duration. We will investigate on the impact of these factors in mitigating the spread of COVID-19 and achieving the goal of herd immunity in the Philippines. We note that the model can estimate the number of undetected infections which may represent the naturally immune members of the population. There are mathematical models developed to study the effect of delays in vaccination for other infectious diseases, e.g. [16, 13]. However, none of these models have an explicit vaccinated compartment with a delay term.

We arranged the rest of the paper as follows: Section 2 describes the model formulation. The analysis of the model is in Section 3. In Section 4, we present the calibrated parameters and show the results of the simulations. And lastly, in Section 5 we have a brief discussion of the results and its implications.

## 2 Mathematical Model

We focus on five compartments dividing a population infected by COVID-19. The compartments are Susceptible (S), Vaccinated (V), Infected (I), Confirmed (C), and Removed (R). Removed here means either recovered or dead. Figure 1 shows the dynamics of the model. A susceptible individual (S) could become infected after contact with infected individuals (I). Infected individuals (I) is considered as a confirmed case (C) after being detected. In either case, the individual will end up in the removed compartment (R). We assume that infected individuals will have full and lasting immunity. But, we assume a waning vaccine-induced immunity. Moreover, we note that vaccinated individuals will not become immune instantaneously but it will involve a delay time $\tau$.

<figure><img src="image_5.png" /><figcaption>Figure 1: The flowchart of the model with a delay term to account for vaccine-induced immunity delay.</figcaption></figure>

Here, we adapted a closed population model, i.e. we do not consider natural birth and death rates. We denoted by $N_0$ the total number of population that is assumed constant. The model is written as frequency dependent, whereas it can be rewritten as density dependent provided that $\beta$ is replaced by $\beta/N_0$. We assume that vaccinated individuals will survive in the span of the needed delay time and become fully immune
With the availability of vaccines, other uncertainties and questions emerged and various mathematical models have been developed and studied to answer some of them. In Jentsch et al. [7], a coupled social-epidemiological model of SARS-CoV-2 transmission to address the problem of vaccine prioritization was developed. The authors provided answers whether to prioritize vaccinating individuals who cause the most transmission or those who are at the highest risk of death. One of their main results is that if vaccines become available sufficiently late in the pandemic, it is relevant to prioritize vaccinating individuals who caused the most transmission than those at higher risk of death (like those aged 60 years and older). In Moore et al. [14], a mathematical model structured by age and the UK region was used to conduct an analysis of a vaccination programme together with the relaxation of NPIs. Their modelling study estimated that vaccination alone is insufficient to contain the outbreak. In [1], Acuña-Zegarra et al. formulated an optimal control problem with mixed constraints to study optimal vaccination policies. Their solution identifies vaccination policies that minimize the burden of COVID-19 quantified by the number of disability-adjusted years of life lost.

In this study, we develop a mathematical model of COVID-19 that is used to investigate the now yet unclear factors related to vaccination. These factors are vaccine-induced immunity delay, vaccination rate, vaccine effectiveness, and vaccine-induced immunity duration. We will investigate on the impact of these factors in mitigating the spread of COVID-19 and achieving the goal of herd immunity in the Philippines. We note that the model can estimate the number of undetected infections which may represent the naturally immune members of the population. There are mathematical models developed to study the effect of delays in vaccination for other infectious diseases, e.g. [16, 13]. However, none of these models have an explicit vaccinated compartment with a delay term.

We arranged the rest of the paper as follows: Section 2 describes the model formulation. The analysis of the model is in Section 3. In Section 4, we present the calibrated parameters and show the results of the simulations. And lastly, in Section 5 we have a brief discussion of the results and its implications.

## 2 Mathematical Model

We focus on five compartments dividing a population infected by COVID-19. The compartments are Susceptible (S), Vaccinated (V), Infected (I), Confirmed (C), and Removed (R). Removed here means either recovered or dead. Figure 1 shows the dynamics of the model. A susceptible individual (S) could become infected after contact with infected individuals (I). Infected individuals (I) is considered as a confirmed case (C) after being detected. In either case, the individual will end up in the removed compartment (R). We assume that infected individuals will have full and lasting immunity. But, we assume a waning vaccine-induced immunity. Moreover, we note that vaccinated individuals will not become immune instantaneously but it will involve a delay time $\tau$.

<figure><img src="image_4.png" /><figcaption>Figure 1: The flowchart of the model with a delay term to account for vaccine-induced immunity delay.</figcaption></figure>

Here, we adapted a closed population model, i.e. we do not consider natural birth and death rates. We denoted by $N_0$ the total number of population that is assumed constant. The model is written as frequency dependent, whereas it can be rewritten as density dependent provided that $\beta$ is replaced by $\beta/N_0$. We assume that vaccinated individuals will survive in the span of the needed delay time and become fully immune
afterwards. The system of differential equations with the delay term is given by

$$ \frac{dS}{dt} = -\alpha\beta I \frac{S}{N_0} - vS(t-\tau) + \eta V \quad (1) $$

$$ \frac{dV}{dt} = vS(t-\tau) - \eta V \quad (2) $$

$$ \frac{dI}{dt} = \alpha\beta I \frac{S}{N_0} - (\gamma + \rho)I \quad (3) $$

$$ \frac{dC}{dt} = \rho I - \delta C \quad (4) $$

$$ \frac{dR}{dt} = \gamma I + \delta C. \quad (5) $$

The parameter $\beta$ denotes the transmission rate of the disease, while the parameter $\alpha$ represents the level of control measures to limit the transmission. The parameter $v$ denotes the vaccination rate, which is further refined in (12). The parameter $\eta$ denotes the vaccine-induced immunity waning rate, while the parameter $\rho$ denotes the proportion of infections detected and confirmed through testing. Lastly, the parameters $\delta$ and $\gamma$ denote the removal rates from $C$ to $R$ and $I$ to $R$, respectively.

# 3 Qualitative analysis

The initial value problem requires to know the past history. We assume that for $t \in [-\tau, 0]$

$$ S(t) = S_0(t), V(t) = V_0(t), I(t) = I_0(t), C(t) = C_0(t), R(t) = R_0(t), $$

where $S_0, V_0, I_0, C_0, R_0$ are continuous functions on $[-\tau, 0]$. It is now standard to check that the problem is well-posed, and the domain $\{(S, V, I, C, R) \in \mathbb{R}_+^5; 0 \le S + V + I + C + R \le N_0\}$ is positively invariant.

The only equilibrium is the Disease Free Equilibrium of the form $(S^*, V^* = \frac{\gamma}{\eta}S^*, I^* = 0, C^* = 0, R^*)$ with $S^* + V^* + R^* = N_0$.

## 3.1 Basic Reproduction Number

Since the infected individuals are in $I$ and $C$, the rate of new infections in each compartment ($\mathcal{F}$) and the rate of other transitions between compartments ($\mathcal{V}$) can be rewritten as

$$ \mathcal{F} = \begin{pmatrix} \alpha\beta I \frac{S}{N_0} \\ 0 \end{pmatrix}, \quad \mathcal{V} = \begin{pmatrix} (\gamma + \rho)I \\ \delta C - \rho I \end{pmatrix}. $$

Thus,

$$ J_{\mathcal{F}} = \begin{pmatrix} \alpha\beta \frac{S}{N_0} & 0 \\ 0 & 0 \end{pmatrix}, J_{\mathcal{V}} = \begin{pmatrix} \gamma + \rho & 0 \\ -\rho & \delta \end{pmatrix}, J_{\mathcal{V}}^{-1} = \begin{pmatrix} \frac{1}{\gamma+\rho} & 0 \\ \frac{\rho}{\delta(\rho+\gamma)} & \frac{1}{\delta} \end{pmatrix}. $$

Therefore, the next generation matrix is

$$ J_{\mathcal{F}} J_{\mathcal{V}}^{-1} = \begin{pmatrix} \frac{\alpha\beta}{\gamma+\rho} \frac{S}{N_0} & 0 \\ 0 & 0 \end{pmatrix}. $$

We deduce as in Diekman et. al [5] and Driessche and Watmough [6] that the basic reproduction number

$$ \mathcal{R}_0 := \frac{\alpha\beta}{\gamma + \rho} \frac{S^*}{N_0}. $$

The term $\frac{\alpha\beta}{\gamma+\rho}$ represents the contact rate with infected individuals during the infection period. Herd immunity is then obtained when $\mathcal{R}_0 < 1$, i.e. if the proportion of non-immune individuals $1 - \frac{S^*}{N_0}$ is smaller than $1 - \frac{\gamma+\rho}{\alpha\beta}$.
afterwards. The system of differential equations with the delay term is given by

$$ \frac{dS}{dt} = -\alpha\beta I \frac{S}{N_0} - \nu S(t - \tau) + \eta V \quad (1) $$

$$ \frac{dV}{dt} = \nu S(t - \tau) - \eta V \quad (2) $$

$$ \frac{dI}{dt} = \alpha\beta I \frac{S}{N_0} - (\gamma + \rho)I \quad (3) $$

$$ \frac{dC}{dt} = \rho I - \delta C \quad (4) $$

$$ \frac{dR}{dt} = \gamma I + \delta C. \quad (5) $$

The parameter $\beta$ denotes the transmission rate of the disease, while the parameter $\alpha$ represents the level of control measures to limit the transmission. The parameter $\nu$ denotes the vaccination rate, which is further refined in (12). The parameter $\eta$ denotes the vaccine-induced immunity waning rate, while the parameter $\rho$ denotes the proportion of infections detected and confirmed through testing. Lastly, the parameters $\delta$ and $\gamma$ denote the removal rates from C to R and I to R, respectively.

# 3 Qualitative analysis

The initial value problem requires to know the past history. We assume that for $t \in [-\tau, 0]$

$$ S(t) = S_0(t), V(t) = V_0(t), I(t) = I_0(t), C(t) = C_0(t), R(t) = R_0(t), $$

where $S_0, V_0, I_0, C_0, R_0$ are continuous functions on $[-\tau, 0]$. It is now standard to check that the problem is well-posed, and the domain $\{(S, V, I, C, R) \in \mathbb{R}_+^5; 0 \le S + V + I + C + R \le N_0\}$ is positively invariant.

The only equilibrium is the Disease Free Equilibrium of the form $(S^*, V^* = \frac{\nu}{\eta}S^*, I^* = 0, C^* = 0, R^*)$ with $S^* + V^* + R^* = N_0$.

## 3.1 Basic Reproduction Number

Since the infected individuals are in I and C, the rate of new infections in each compartment ($\mathcal{F}$) and the rate of other transitions between compartments ($\mathcal{V}$) can be rewritten as

$$ \mathcal{F} = \begin{pmatrix} \alpha\beta I \frac{S}{N_0} \\ 0 \end{pmatrix}, \quad \mathcal{V} = \begin{pmatrix} (\gamma + \rho)I \\ \delta C - \rho I \end{pmatrix}. $$

Thus,

$$ J_{\mathcal{F}} = \begin{pmatrix} \alpha\beta \frac{S}{N_0} & 0 \\ 0 & 0 \end{pmatrix}, J_{\mathcal{V}} = \begin{pmatrix} \gamma + \rho & 0 \\ -\rho & \delta \end{pmatrix}, J_{\mathcal{V}}^{-1} = \begin{pmatrix} \frac{1}{\gamma + \rho} & 0 \\ \frac{\rho}{\delta(\rho + \gamma)} & \frac{1}{\delta} \end{pmatrix}. $$

Therefore, the next generation matrix is

$$ J_{\mathcal{F}} J_{\mathcal{V}}^{-1} = \begin{pmatrix} \frac{\alpha\beta}{\gamma+\rho} \frac{S}{N_0} & 0 \\ 0 & 0 \end{pmatrix}. $$

We deduce as in Diekman et. al [5] and Driessche and Watmough [6] that the basic reproduction number

$$ \mathcal{R}_0 := \frac{\alpha\beta}{\gamma + \rho} \frac{S^*}{N_0}. $$

The term $\frac{\alpha\beta}{\gamma+\rho}$ represents the contact rate with infected individuals during the infection period. Herd immunity is then obtained when $\mathcal{R}_0 < 1$, i.e. if the proportion of non-immune individuals $1 - \frac{S^*}{N_0}$ is smaller than $1 - \frac{\gamma+\rho}{\alpha\beta}$.
## 3.2 Stability of the equilibrium points with delay

Since the population is constant, we restrict ourselves to the system $(S, I, C, R)$, since $V = N_0 - (S + I + C + R)$. Considering a delay $\tau > 0$, the characteristic equation then is given by

$$ \det(\lambda I_d + J_0 + e^{-\lambda \tau} J_1) = 0, $$

with

$$ J_0 = Jac \begin{pmatrix} -\alpha \beta I_{\frac{S}{N_0}} + \eta(N_0 - (S + I + C + R)) \\ \alpha \beta I_{\frac{S}{N_0}} - (\gamma + \rho)I \\ \rho I - \delta C \\ \gamma I + \delta C \end{pmatrix} (S^*, 0, 0, R^*) = \begin{pmatrix} -\eta & -\alpha \beta \frac{S^*}{N_0} - \eta & -\eta & -\eta \\ 0 & \alpha \beta \frac{S^*}{N_0} - (\gamma + \rho) & 0 & 0 \\ 0 & \rho & -\delta & 0 \\ 0 & \gamma & \delta & 0 \end{pmatrix}, $$

and

$$ J_1 = Jac \begin{pmatrix} -\nu S \\ 0 \\ 0 \\ 0 \end{pmatrix} = \begin{pmatrix} -\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}. $$

It comes

$$ \lambda(\lambda + \delta) \left( \lambda - \alpha \beta \frac{S^*}{N_0} + (\gamma + \rho) \right) \left( \lambda + \eta + \nu e^{-\lambda \tau} \right) = 0. $$

Real eigenvalues are $0, -\delta, (\gamma + \rho)(\mathbb{R}_0 - 1)$. One can note that there is no real solution of the equation $(\lambda + \eta + \nu e^{-\lambda \tau}) = 0$. Let $\lambda = a + ib$, we deduce

$$ e^{-a\tau} \cos(b\tau) = -\frac{a+\eta}{\nu} \qquad (6) $$

$$ e^{-a\tau} \sin(b\tau) = \frac{b}{\nu}. \qquad (7) $$

Studying the sign of $a$, we obtain the following result.

**Theorem 3.1.** 1. If $\mathbb{R}_0 < 1$ and $0 \le \tau < \frac{\pi}{2\nu}$, then the DFE is asymptotically stable.

2. If $\tau > \frac{\pi}{2\nu}$, then the DFE is unstable.

# 4 Simulations

## 4.1 Parameter values

From [8], the infection period has been estimated at 12 days and so we set $\gamma = 1/12$. Symptoms usually appear after 5 days, and we assume the detection is then performed. The parameter $\delta$ is set to 1/7.

Parameters $\alpha, \beta$ and $\rho$, are estimated by fitting the model with the cumulative confirmed cases of the Philippines from January 30, 2020 to February 28, 2021, the last day before the start of the vaccination campaign of the Philippines. The data set can be downloaded from https://data.gov.ph. These data are available to the public and so ethical approval is not required.

We acknowledge that the three parameters are varying through time due to the varying controls being implemented and of the evolution of the virus. Thus, our parameter estimation are piecewise functions as follows:

$$ \alpha = \begin{cases} \alpha_1, & \text{for } 0 \le t < 46 \ (\text{January 30, 2020} - \text{March 15, 2020}) \\ \alpha_2, & \text{for } 46 \le t < 123 \ (\text{March 16, 2020} - \text{May 31, 2020}) \\ \alpha_3, & \text{for } 123 \le t < 320 \ (\text{June 1, 2020} - \text{Dec. 14, 2020}) \\ \alpha_4, & \text{for } 320 \le t < 340 \ (\text{Dec. 15, 2020} - \text{Jan. 3, 2021}) \\ \alpha_5, & \text{for } 340 \le t < 396 \ (\text{Jan. 4, 2021} - \text{Feb. 28, 2021}) \\ \alpha_6, & \text{for } t \ge 396, \ (\text{March 1, 2021 onwards}) \end{cases} \qquad (8) $$
## 3.2 Stability of the equilibrium points with delay

Since the population is constant, we restrict ourselves to the system $(S, I, C, R)$, since $V = N_0 - (S + I + C + R)$. Considering a delay $\tau > 0$, the characteristic equation then is given by

$$ \det(\lambda I_d + J_0 + e^{-\lambda \tau} J_1) = 0, $$

with

$$ J_0 = Jac \begin{pmatrix} -\alpha \beta I_{\frac{S}{N_0}} + \eta(N_0 - (S + I + C + R)) \\ \alpha \beta I_{\frac{S}{N_0}} - (\gamma + \rho)I \\ \rho I - \delta C \\ \gamma I + \delta C \end{pmatrix} (S^*, 0, 0, R^*) = \begin{pmatrix} -\eta & -\alpha \beta \frac{S^*}{N_0} - \eta & -\eta & -\eta \\ 0 & \alpha \beta \frac{S^*}{N_0} - (\gamma + \rho) & 0 & 0 \\ 0 & \rho & -\delta & 0 \\ 0 & \gamma & \delta & 0 \end{pmatrix}, $$

and

$$ J_1 = Jac \begin{pmatrix} -\nu S \\ 0 \\ 0 \\ 0 \end{pmatrix} = \begin{pmatrix} -\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}. $$

It comes

$$ \lambda(\lambda + \delta) \left( \lambda - \alpha \beta \frac{S^*}{N_0} + (\gamma + \rho) \right) \left( \lambda + \eta + \nu e^{-\lambda \tau} \right) = 0. $$

Real eigenvalues are $0, -\delta, (\gamma + \rho)(\mathbb{R}_0 - 1)$. One can note that there is no real solution of the equation $(\lambda + \eta + \nu e^{-\lambda \tau}) = 0$. Let $\lambda = a + ib$, we deduce

$$ e^{-a\tau} \cos(b\tau) = -\frac{a+\eta}{\nu} \qquad (6) $$

$$ e^{-a\tau} \sin(b\tau) = \frac{b}{\nu}. \qquad (7) $$

Studying the sign of $a$, we obtain the following result.

**Theorem 3.1.** 1. If $\mathbb{R}_0 < 1$ and $0 \le \tau < \frac{\pi}{2\nu}$, then the DFE is asymptotically stable.

2. If $\tau > \frac{\pi}{2\nu}$, then the DFE is unstable.

# 4 Simulations

## 4.1 Parameter values

From [8], the infection period has been estimated at 12 days and so we set $\gamma = 1/12$. Symptoms usually appear after 5 days, and we assume the detection is then performed. The parameter $\delta$ is set to 1/7.

Parameters $\alpha, \beta$ and $\rho$, are estimated by fitting the model with the cumulative confirmed cases of the Philippines from January 30, 2020 to February 28, 2021, the last day before the start of the vaccination campaign of the Philippines. The data set can be downloaded from https://data.gov.ph. These data are available to the public and so ethical approval is not required.

We acknowledge that the three parameters are varying through time due to the varying controls being implemented and of the evolution of the virus. Thus, our parameter estimation are piecewise functions as follows:

$$ \alpha = \begin{cases} \alpha_1, & \text{for } 0 \le t < 46 \ (\text{January 30, 2020} - \text{March 15, 2020}) \\ \alpha_2, & \text{for } 46 \le t < 123 \ (\text{March 16, 2020} - \text{May 31, 2020}) \\ \alpha_3, & \text{for } 123 \le t < 320 \ (\text{June 1, 2020} - \text{Dec. 14, 2020}) \\ \alpha_4, & \text{for } 320 \le t < 340 \ (\text{Dec. 15, 2020} - \text{Jan. 3, 2021}) \\ \alpha_5, & \text{for } 340 \le t < 396 \ (\text{Jan. 4, 2021} - \text{Feb. 28, 2021}) \\ \alpha_6, & \text{for } t \ge 396, \ (\text{March 1, 2021 onwards}) \end{cases} \qquad (8) $$
where the dates correspond to the noticeable changes in the control measures of the country and we put in $\alpha_6$ our estimate for the future transmission reduction control,

$$ \beta = \begin{cases} \beta_1, & \text{for } 0 \le t < 343 \text{ (January 30, 2020 - Jan. 6, 2021)} \\ \beta_1(1 + \beta_2), & \text{for } t \ge 343, \text{ (Jan. 7, 2021 onwards)} \end{cases} \qquad (9) $$

where $\beta_2$ corresponds to the increase of the transmission rate due to the introduction of the more transmissible new variant of the virus in the population,

$$ \rho = \begin{cases} \rho_1, & \text{for } 0 \le t < 184 \text{ (January 30, 2020 - July 31, 2020)} \\ \rho_2, & \text{for } t \ge 184, \text{ (Aug. 1, 2020 onwards)} \end{cases} \qquad (10) $$

where $\rho_2$ corresponds to the start of the current testing rate of the country.

The cost function to be minimized for the parameter estimation consists of a nonlinear least square function given by

$$ F(\theta) = \sum_{i=1}^{N} (T_d(t_i) - T_s(t_i, \theta))^2, \qquad (11) $$

where $\theta = (\alpha, \beta, \rho)$ and constrained by $0 \le \theta \le 1$. $T_d(t_i)$ and $T_s(t_i, \theta)$ denote the cumulative confirmed cases from data and model at time $t_i$, respectively. The optimization problem is solved using the Levenberg-Marquardt algorithm [9, 12].

The obtained parameter values are given in Table 1 and the fit is shown in Figure 2. With the estimated parameters, the relative cost function value is $1.3058 \times 10^{-3}$.

<figure><img src="image_6.png" /><figcaption>**Figure 2:** The model's fit with data and the infection curve including the undetected ones. With $N_0 = 108, 116, 615, I_0 = 2, C_0 = 1, S_0 = N_0 - I_0 - C_0, R_0 = 0$, the model predicted a cumulative infection (including the undetected ones) of 2,387,393 against the confirmed cases of 576,352 by Feb. 28, 2021.</figcaption></figure>

The vaccination rate is given by

$$ v(t) = \begin{cases} v_e \left( \frac{v_p}{N_0} \right), & \text{for } t > t_v \\ 0, & \text{for } t < t_v, \end{cases} \qquad (12) $$

where $v_e$ is the parameter for vaccine effectiveness and $\frac{v_p}{N_0}$ is the proportion of vaccinated susceptibles per day, and $t_v$ is the start time of vaccination. The parameters are summarized in Table 1.

## 4.2 Effect of vaccination on disease progression

We study the possible impact of vaccination, with its related factors, in controlling COVID-19 in the country. Considering that some travel restrictions are being lifted and easing of community quarantines are already
where the dates correspond to the noticeable changes in the control measures of the country and we put in $\alpha_6$ our estimate for the future transmission reduction control,

$$ \beta = \begin{cases} \beta_1, & \text{for } 0 \le t < 343 \text{ (January 30, 2020 - Jan. 6, 2021)} \\ \beta_1(1 + \beta_2), & \text{for } t \ge 343, \text{ (Jan. 7, 2021 onwards)} \end{cases} \qquad (9) $$

where $\beta_2$ corresponds to the increase of the transmission rate due to the introduction of the more transmissible new variant of the virus in the population,

$$ \rho = \begin{cases} \rho_1, & \text{for } 0 \le t < 184 \text{ (January 30, 2020 - July 31, 2020)} \\ \rho_2, & \text{for } t \ge 184, \text{ (Aug. 1, 2020 onwards)} \end{cases} \qquad (10) $$

where $\rho_2$ corresponds to the start of the current testing rate of the country.

The cost function to be minimized for the parameter estimation consists of a nonlinear least square function given by

$$ F(\theta) = \sum_{i=1}^{N} (T_d(t_i) - T_s(t_i, \theta))^2, \qquad (11) $$

where $\theta = (\alpha, \beta, \rho)$ and constrained by $0 \le \theta \le 1$. $T_d(t_i)$ and $T_s(t_i, \theta)$ denote the cumulative confirmed cases from data and model at time $t_i$, respectively. The optimization problem is solved using the Levenberg-Marquardt algorithm [9, 12].

The obtained parameter values are given in Table 1 and the fit is shown in Figure 2. With the estimated parameters, the relative cost function value is $1.3058 \times 10^{-3}$.

<figure><img src="image_5.png" /><figcaption>**Figure 2:** The model's fit with data and the infection curve including the undetected ones. With $N_0 = 108, 116, 615, I_0 = 2, C_0 = 1, S_0 = N_0 - I_0 - C_0, R_0 = 0$, the model predicted a cumulative infection (including the undetected ones) of 2,387,393 against the confirmed cases of 576,352 by Feb. 28, 2021.</figcaption></figure>

The vaccination rate is given by

$$ v(t) = \begin{cases} v_e \left( \frac{v_p}{N_0} \right), & \text{for } t > t_v \\ 0, & \text{for } t < t_v, \end{cases} \qquad (12) $$

where $v_e$ is the parameter for vaccine effectiveness and $\frac{v_p}{N_0}$ is the proportion of vaccinated susceptibles per day, and $t_v$ is the start time of vaccination. The parameters are summarized in Table 1.

## 4.2 Effect of vaccination on disease progression

We study the possible impact of vaccination, with its related factors, in controlling COVID-19 in the country. Considering that some travel restrictions are being lifted and easing of community quarantines are already
<table><caption>**Table 1:** The parameters of the model.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>ν<sub>p</sub></td><td>number of vaccinated individuals</td><td>persons/day</td><td>0 – 500000</td></tr><tr><td>ν<sub>e</sub></td><td>vaccination effectiveness</td><td>dimensionless</td><td>0.5 – 1</td></tr><tr><td>η</td><td>vaccine-induced immunity waning rate</td><td>1/day</td><td><sup>1</sup>&frasl;<sub>600</sub> – <sup>1</sup>&frasl;<sub>90</sub></td></tr><tr><td>τ</td><td>vaccine-induced immunity delay</td><td>days</td><td>14 – 98</td></tr><tr><td rowspan="4">α</td><td rowspan="4">transmission reduction</td><td rowspan="4">dimensionless</td><td>α<sub>1</sub> | 1</td></tr><tr><td>α<sub>2</sub> | 0.48181</td></tr><tr><td>α<sub>3</sub> | 0.41553</td></tr><tr><td>α<sub>4</sub> | 0.60158</td></tr><tr><td>α<sub>5</sub></td><td></td><td></td><td>0.41712</td></tr><tr><td rowspan="2">β</td><td rowspan="2">transmission rate</td><td rowspan="2">1/day</td><td>β<sub>1</sub> | 0.26467</td></tr><tr><td>β<sub>2</sub> | 3.8981 × 10<sup>-7</sup></td></tr><tr><td>ρ</td><td>proportion of infections detected</td><td>1/day</td><td></td></tr><tr><td>ρ<sub>1</sub></td><td>0.01154</td><td></td><td></td></tr><tr><td>ρ<sub>2</sub></td><td>0.036766</td><td></td><td></td></tr><tr><td>δ</td><td>removal rate from C to R</td><td>1/day</td><td>1/7</td></tr><tr><td>γ</td><td>removal rate from I to R</td><td>1/day</td><td>1/12</td></tr></tbody></table>

being done in many parts of the country, we let $α_6 = α_5 \times 1.2$. With the vaccination campaign starting at day 396, Figure 3 shows the impact of vaccination considering vaccination rate ($ν_p$), vaccine effectiveness ($ν_e$), vaccine-induced immunity delay ($τ$), and vaccine-induced immunity duration ($1/η$). We measure this impact by the maximum value of the Confirmed compartment.

In the simulations in Figure 3, the ranges are [0, 500000], [0.5, 1], [14, 98] (days), and [90, 600] (days) for $ν_p$, $ν_e$, $τ$, and $1/η$, respectively, which correspond to the probable values of the parameters for the Philippines. These values are represented in percentile. When these parameters are needed to be fixed, we used the values $ν_p = 100000$, $ν_e = 0.9$, $τ = 56$, and $η = 0$. From the start of vaccination, we let the simulations run for an additional 671 days or up to Dec. 31, 2022.

<figure><img src="image_7.png" /><figcaption>**Figure 3:** The effect of each of the vaccine-related parameters $ν_p$, $ν_e$, $τ$, and $η$ to the maximum value of the Confirmed (C) compartment. The x-axis is in terms of the percentile of the ranges [0, 500000], [0.5, 1], [14, 98] (days), and [90, 600] (days), for $ν_p$, $ν_e$, $τ$, and $1/η$, respectively.</figcaption></figure>
<table><caption>**Table 1:** The parameters of the model.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>ν<sub>p</sub></td><td>number of vaccinated individuals</td><td>persons/day</td><td>0 – 500000</td></tr><tr><td>ν<sub>e</sub></td><td>vaccination effectiveness</td><td>dimensionless</td><td>0.5 – 1</td></tr><tr><td>η</td><td>vaccine-induced immunity waning rate</td><td>1/day</td><td><sup>1</sup>&frasl;<sub>600</sub> – <sup>1</sup>&frasl;<sub>90</sub></td></tr><tr><td>τ</td><td>vaccine-induced immunity delay</td><td>days</td><td>14 – 98</td></tr><tr><td rowspan="5">α</td><td rowspan="5">transmission reduction</td><td rowspan="5">dimensionless</td><td>α<sub>1</sub> | 1</td></tr><tr><td>α<sub>2</sub> | 0.48181</td></tr><tr><td>α<sub>3</sub> | 0.41553</td></tr><tr><td>α<sub>4</sub> | 0.60158</td></tr><tr><td>α<sub>5</sub> | 0.41712</td></tr><tr><td rowspan="2">β</td><td rowspan="2">transmission rate</td><td rowspan="2">1/day</td><td>β<sub>1</sub> | 0.26467</td></tr><tr><td>β<sub>2</sub> | 3.8981 × 10<sup>-7</sup></td></tr><tr><td>ρ</td><td>proportion of infections detected</td><td>1/day</td><td></td></tr><tr><td>ρ<sub>1</sub></td><td>0.01154</td><td></td><td></td></tr><tr><td>ρ<sub>2</sub></td><td>0.036766</td><td></td><td></td></tr><tr><td>δ</td><td>removal rate from C to R</td><td>1/day</td><td>1/7</td></tr><tr><td>γ</td><td>removal rate from I to R</td><td>1/day</td><td>1/12</td></tr></tbody></table>

being done in many parts of the country, we let $α_6 = α_5 \times 1.2$. With the vaccination campaign starting at day 396, Figure 3 shows the impact of vaccination considering vaccination rate ($ν_p$), vaccine effectiveness ($ν_e$), vaccine-induced immunity delay ($τ$), and vaccine-induced immunity duration ($1/η$). We measure this impact by the maximum value of the Confirmed compartment.

In the simulations in Figure 3, the ranges are [0, 500000], [0.5, 1], [14, 98] (days), and [90, 600] (days) for $ν_p$, $ν_e$, $τ$, and $1/η$, respectively, which correspond to the probable values of the parameters for the Philippines. These values are represented in percentile. When these parameters are needed to be fixed, we used the values $ν_p = 100000$, $ν_e = 0.9$, $τ = 56$, and $η = 0$. From the start of vaccination, we let the simulations run for an additional 671 days or up to Dec. 31, 2022.

<figure><img src="image_6.png" /><figcaption>**Figure 3:** The effect of each of the vaccine-related parameters $ν_p$, $ν_e$, $τ$, and $η$ to the maximum value of the Confirmed (C) compartment. The x-axis is in terms of the percentile of the ranges [0, 500000], [0.5, 1], [14, 98] (days), and [90, 600] (days), for $ν_p$, $ν_e$, $τ$, and $1/η$, respectively.</figcaption></figure>
## 4.3 Effect of vaccination on herd immunity

One of the goals of vaccination is for the country to reach herd immunity against the virus. We compute the day when the country will reach herd immunity considering the various factors related to vaccination. As shown in Section 3.1, the population has achieved herd immunity if more than $1 - \frac{\gamma+\rho}{a\beta}$ of its susceptible population is already immune to the virus. Note that we assume here that individuals who recovered from the infection will develop permanent immunity and that vaccination will be given only to the susceptible individuals.

Considering that $v_p$ is the dominant factor among the vaccine-related parameters, we pair it with the other vaccine-related parameters to produce the heatmaps in Figures 4, 5, and 6, showing the day when the country will have reached the herd immunity threshold.

<figure><img src="image_8.png" /><figcaption>**Figure 4:** A heatmap showing the days after Day 0 (Jan. 30, 2020) when the country reaches the herd immunity threshold, considering the combination of values of $v_p$ and $1/\eta$. The parameters $v_e$ and $\tau$ are fixed at 0.9 and 56 days, respectively. Black boxes mean herd immunity is not achieved at the end of simulation (Dec. 31, 2022).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5:** A heatmap showing the days after Day 0 (Jan. 30, 2020) when the country reaches the herd immunity threshold, considering the combination of values of $v_p$ and $v_e$. The parameters $\tau$ and $\eta$ are fixed at 56 days and 0, respectively. Black boxes mean herd immunity is not achieved at the end of simulation (Dec. 31, 2022).</figcaption></figure>
## 4.3 Effect of vaccination on herd immunity

One of the goals of vaccination is for the country to reach herd immunity against the virus. We compute the day when the country will reach herd immunity considering the various factors related to vaccination. As shown in Section 3.1, the population has achieved herd immunity if more than $1 - \frac{\gamma+\rho}{a\beta}$ of its susceptible population is already immune to the virus. Note that we assume here that individuals who recovered from the infection will develop permanent immunity and that vaccination will be given only to the susceptible individuals.

Considering that $v_p$ is the dominant factor among the vaccine-related parameters, we pair it with the other vaccine-related parameters to produce the heatmaps in Figures 4, 5, and 6, showing the day when the country will have reached the herd immunity threshold.

<figure><img src="image_7.png" /><figcaption>**Figure 4:** A heatmap showing the days after Day 0 (Jan. 30, 2020) when the country reaches the herd immunity threshold, considering the combination of values of $v_p$ and $1/\eta$. The parameters $v_e$ and $\tau$ are fixed at 0.9 and 56 days, respectively. Black boxes mean herd immunity is not achieved at the end of simulation (Dec. 31, 2022).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5:** A heatmap showing the days after Day 0 (Jan. 30, 2020) when the country reaches the herd immunity threshold, considering the combination of values of $v_p$ and $v_e$. The parameters $\tau$ and $\eta$ are fixed at 56 days and 0, respectively. Black boxes mean herd immunity is not achieved at the end of simulation (Dec. 31, 2022).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6:** A heatmap showing the days after Day 0 (Jan. 30, 2020) when the country reaches the herd immunity threshold, considering the combination of values of $v_p$ and $\tau$. The parameters $v_e$ and $\eta$ are fixed at 0.9 and 0, respectively. Black boxes mean herd immunity is not achieved at the end of simulation (Dec. 31, 2022).</figcaption></figure>

## 4.4 Comparing the different vaccines

We consider the different available vaccines distinguish by vaccine effectiveness and vaccine-induced immunity delay. We fixed $v_p = 100000$ and $\eta = 0$. The parameters $v_e$ and $\tau$ corresponding to the different vaccines are given in Table 2. We wanted to see the effect of the different vaccines to the number of susceptible, the fully immune individuals through vaccination, and the total number of infections. The results are shown in Figures 7, 8, and 9. The values at the end of simulation are tabulated in Table 2.

<table><caption>**Table 2:** Comparing the different vaccines. $S(t_f)$, $V(t_f)$, and $I_{cum}(t_f)$ are the number of susceptible, fully immune individuals through vaccination, and cumulative infections, respectively, at the final time of simulation, which is Day 1066 (Dec. 31, 2022).</caption><thead><tr><td>Vaccine</td><td>$v_e$</td><td>$\tau$</td><td>$S(t_f)$</td><td>$V(t_f)$</td><td>$I_{cum}(t_f)$</td></tr></thead><tbody><tr><td>Pfizer-BioNTech</td><td>0.95</td><td>42</td><td>58,589,579</td><td>45,498,723</td><td>4,030,627</td></tr><tr><td>Moderna</td><td>0.94</td><td>56</td><td>59,126,976</td><td>44,705,577</td><td>4,286,529</td></tr><tr><td>Johnson</td><td>0.66</td><td>28</td><td>69,501,815</td><td>34,244,861</td><td>4,372,459</td></tr><tr><td>AstraZeneca</td><td>0.7</td><td>56</td><td>68,238,781</td><td>35,081,867</td><td>4,798,741</td></tr><tr><td>Sputnik</td><td>0.91</td><td>42</td><td>60,036,593</td><td>43,992,393</td><td>4,089,979</td></tr><tr><td>Sinovac</td><td>0.51</td><td>28</td><td>75,746,331</td><td>27,434,886</td><td>4,938,258</td></tr><tr><td>Bharat</td><td>0.78</td><td>56</td><td>65,107,767</td><td>38,419,780</td><td>4,591,718</td></tr></tbody></table>

## 5 Discussion

Even though various non-pharmaceutical interventions have been implemented worldwide, it is widely believed that vaccination is the main tool to end the COVID-19 pandemic and bring back normalcy to the world. Every country's vaccination campaign aims at achieving herd immunity. This immunity could be achieved by being infected of the virus (natural immunity) or being vaccinated (vaccine-induced). Natural immunity thus has a positive effect on the outlook of herd immunity in a country. In Figure 2, the calibrated model estimated that the naturally-immune individuals in the Philippines could be around four times the recorded cases. A test before vaccination could reveal these individuals and the vaccines intended for them could instead be used to the completely susceptible individuals. Figure 3 shows us that the number of vaccines given per day
<figure><img src="image_9.png" /><figcaption>**Figure 6:** A heatmap showing the days after Day 0 (Jan. 30, 2020) when the country reaches the herd immunity threshold, considering the combination of values of $v_p$ and $\tau$. The parameters $v_e$ and $\eta$ are fixed at 0.9 and 0, respectively. Black boxes mean herd immunity is not achieved at the end of simulation (Dec. 31, 2022).</figcaption></figure>

## 4.4 Comparing the different vaccines

We consider the different available vaccines distinguish by vaccine effectiveness and vaccine-induced immunity delay. We fixed $v_p = 100000$ and $\eta = 0$. The parameters $v_e$ and $\tau$ corresponding to the different vaccines are given in Table 2. We wanted to see the effect of the different vaccines to the number of susceptible, the fully immune individuals through vaccination, and the total number of infections. The results are shown in Figures 7, 8, and 9. The values at the end of simulation are tabulated in Table 2.

<table><caption>**Table 2:** Comparing the different vaccines. $S(t_f)$, $V(t_f)$, and $I_{cum}(t_f)$ are the number of susceptible, fully immune individuals through vaccination, and cumulative infections, respectively, at the final time of simulation, which is Day 1066 (Dec. 31, 2022).</caption><thead><tr><th>Vaccine</th><th>v<sub>e</sub></th><th>τ</th><th>S(t<sub>f</sub>)</th><th>V(t<sub>f</sub>)</th><th>I<sub>cum</sub>(t<sub>f</sub>)</th></tr></thead><tbody><tr><td>Pfizer-BioNTech</td><td>0.95</td><td>42</td><td>58,589,579</td><td>45,498,723</td><td>4,030,627</td></tr><tr><td>Moderna</td><td>0.94</td><td>56</td><td>59,126,976</td><td>44,705,577</td><td>4,286,529</td></tr><tr><td>Johnson</td><td>0.66</td><td>28</td><td>69,501,815</td><td>34,244,861</td><td>4,372,459</td></tr><tr><td>AstraZeneca</td><td>0.7</td><td>56</td><td>68,238,781</td><td>35,081,867</td><td>4,798,741</td></tr><tr><td>Sputnik</td><td>0.91</td><td>42</td><td>60,036,593</td><td>43,992,393</td><td>4,089,979</td></tr><tr><td>Sinovac</td><td>0.51</td><td>28</td><td>75,746,331</td><td>27,434,886</td><td>4,938,258</td></tr><tr><td>Bharat</td><td>0.78</td><td>56</td><td>65,107,767</td><td>38,419,780</td><td>4,591,718</td></tr></tbody></table>

## 5 Discussion

Even though various non-pharmaceutical interventions have been implemented worldwide, it is widely believed that vaccination is the main tool to end the COVID-19 pandemic and bring back normalcy to the world. Every country's vaccination campaign aims at achieving herd immunity. This immunity could be achieved by being infected of the virus (natural immunity) or being vaccinated (vaccine-induced). Natural immunity thus has a positive effect on the outlook of herd immunity in a country. In Figure 2, the calibrated model estimated that the naturally-immune individuals in the Philippines could be around four times the recorded cases. A test before vaccination could reveal these individuals and the vaccines intended for them could instead be used to the completely susceptible individuals. Figure 3 shows us that the number of vaccines given per day
<figure><img src="image_11.png" /><figcaption>**Figure 7:** Comparing the different vaccines by the decrease in susceptible.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8:** Comparing the different vaccines by the number of fully immune individuals through vaccination.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9:** Comparing the different vaccines by the number of cumulative infections.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7:** Comparing the different vaccines by the decrease in susceptible.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8:** Comparing the different vaccines by the number of fully immune individuals through vaccination.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9:** Comparing the different vaccines by the number of cumulative infections.</figcaption></figure>
is the dominant determining factor in using vaccination to curb the disease. It also shows that the greater the delay in having vaccine-induced immunity the worse it will be with respect to the number of infections.

Assuming 90% vaccine effectiveness, unlimited vaccine-induced immunity duration, and 56 days vaccine-induced immunity delay, Figures 4, 5, and 6 show that per day vaccination of 150,000 individuals may not be enough for the country to achieve herd immunity by Dec. 31, 2022. This stresses again the need for the country to strengthen its vaccine procurement strategy.

A dark horse in this vaccination campaign is the vaccine-induced immunity duration. Figure 4 shows its impact to herd immunity. Natural immunity could cause herd immunity. However, due to the various non-pharmaceutical interventions to control infections, it is clear that vaccine-induced immunity will be the dominant factor in achieving herd immunity. Figure 10 shows a scenario where the country will have already eliminated the disease but still could not achieve herd immunity even with a sustained constant vaccination. This is due to the system going to a disease-free equilibrium where the number of susceptible is still above the herd immunity threshold. Thus, studies determining the immunity duration cause by the different available vaccines is very important so that countries could design effective vaccination strategies now and for the next years. For example, if achieving herd immunity seems impossible with the current rate, then the vaccination campaign should focus on protecting the vulnerable - minimizing deaths and hospitalization, instead of aiming for herd immunity.

Finally, from Figures 7, 8, and 9, we see that Sinovac has the worst effect among the vaccines. Due to its low vaccine effectiveness, in the simulations it only achieved a little more than half of what Pfizer-BioNTech has achieved with respect to the number of fully immune individuals through vaccination, as can be seen in Table 2. Pfizer-BioNTech is the best vaccine for decreasing the number of susceptible and infections, while increasing the fully immune individuals. It is followed closely by Moderna and Sputnik. The other vaccines are in between.

<figure><img src="image_14.png" /><figcaption>**Figure 10:** Simulation with Pfizer-BioNTech with 300,000 vaccinations per day, but with vaccine-induced immunity for only 9 months. At the end of simulation (Day 1797 or Dec. 31, 2024), $S = 61,286,308$, $I = 0$, $I_{cum} = 3,219,785$, $V = 43,615,996$, and $R = 3,218,590$.</figcaption></figure>

**Acknowledgement:** RLC, MNR, GSGL and EBM are funded by DOST's PCHRD and SEI, and MSU-IIT's OVCRE through PRISM, while PdC and YM are funded by the Agence National de la Recherche and Région Hauts-de-France, projet Space-covid ANR Résilience.

**Author's Statement:** Conflict of interest: Authors state no conflict of interest.
is the dominant determining factor in using vaccination to curb the disease. It also shows that the greater the delay in having vaccine-induced immunity the worse it will be with respect to the number of infections.

Assuming 90% vaccine effectiveness, unlimited vaccine-induced immunity duration, and 56 days vaccine-induced immunity delay, Figures 4, 5, and 6 show that per day vaccination of 150,000 individuals may not be enough for the country to achieve herd immunity by Dec. 31, 2022. This stresses again the need for the country to strengthen its vaccine procurement strategy.

A dark horse in this vaccination campaign is the vaccine-induced immunity duration. Figure 4 shows its impact to herd immunity. Natural immunity could cause herd immunity. However, due to the various non-pharmaceutical interventions to control infections, it is clear that vaccine-induced immunity will be the dominant factor in achieving herd immunity. Figure 10 shows a scenario where the country will have already eliminated the disease but still could not achieve herd immunity even with a sustained constant vaccination. This is due to the system going to a disease-free equilibrium where the number of susceptible is still above the herd immunity threshold. Thus, studies determining the immunity duration cause by the different available vaccines is very important so that countries could design effective vaccination strategies now and for the next years. For example, if achieving herd immunity seems impossible with the current rate, then the vaccination campaign should focus on protecting the vulnerable - minimizing deaths and hospitalization, instead of aiming for herd immunity.

Finally, from Figures 7, 8, and 9, we see that Sinovac has the worst effect among the vaccines. Due to its low vaccine effectiveness, in the simulations it only achieved a little more than half of what Pfizer-BioNTech has achieved with respect to the number of fully immune individuals through vaccination, as can be seen in Table 2. Pfizer-BioNTech is the best vaccine for decreasing the number of susceptible and infections, while increasing the fully immune individuals. It is followed closely by Moderna and Sputnik. The other vaccines are in between.

<figure><img src="image_13.png" /><figcaption>**Figure 10:** Simulation with Pfizer-BioNTech with 300,000 vaccinations per day, but with vaccine-induced immunity for only 9 months. At the end of simulation (Day 1797 or Dec. 31, 2024), $S = 61,286,308$, $I = 0$, $I_{cum} = 3,219,785$, $V = 43,615,996$, and $R = 3,218,590$.</figcaption></figure>

**Acknowledgement:** RLC, MNR, GSGL and EBM are funded by DOST's PCHRD and SEI, and MSU-IIT's OVCRE through PRISM, while PdC and YM are funded by the Agence National de la Recherche and Région Hauts-de-France, projet Space-covid ANR Résilience.

**Author's Statement:** Conflict of interest: Authors state no conflict of interest.
References

[1] Acuña-Zegarra, M. A., Díaz-Infante, S., Baca-Carrasco, D., & Liceaga, D. O. (2021). COVID-19 optimal vaccination policies: a modeling study on effectiveness, natural and vaccine-induced immunity responses. Mathematical Biosciences, 108614.

[2] Ali, M., Shah, S. T. H., Imran, M., & Khan, A. (2020). The role of asymptomatic class, quarantine and isolation in the transmission of COVID-19. Journal of biological dynamics, 14(1), 389-408.

[3] Arcede, J. P., Caga-anan, R. L., Mentuda, C. Q., & Mammeri, Y. (2020). Accounting for Symptomatic and Asymptomatic in a SEIR-type model of COVID-19. Mathematical Modelling of Natural Phenomena, 15, 34.

[4] Buhat, C. A. H., Duero, J. C. C., Felix, E. F. O., Rabajante, J. F., & Mamplata, J. B. (2021). Optimal allocation of COVID-19 test kits among accredited testing centers in the Philippines. Journal of healthcare informatics research, 5(1), 54-69.

[5] Diekmann, J. Heesterbeek, & J. Metz, On the definition and the computation of the basic reproductive ratio $r_0$ in models for infectious diseases in heterogeneous populations, *J Math Biol* **28** (1990), 365–383.

[6] V. Driessche, & J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, *Math Biosci*, **180** (2002), 29-48.

[7] Jentsch, P. C., Anand, M., & Bauch, C. T. (2021). Prioritising COVID-19 vaccination in changing social and epidemiological landscapes: a mathematical modelling study. The Lancet Infectious Diseases.

[8] S. A. Lauer, K. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R. Meredith, A. S. Azman, N. G. Reich, and J. Lessler. The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmed cases: Estimation and application. Annals of Internal Medicine, 03 2020. URL https://doi.org/10.7326/M20-0504.

[9] Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares. Quarterly of applied mathematics, 2(2), 164-168.

[10] Macalisang, J. M., Caay, M. L., Arcede, J. P., & Caga-anan, R. L. (2020). Optimal Control for a COVID-19 Model Accounting for Symptomatic and Asymptomatic. Computational and Mathematical Biophysics, 8(1), 168-179.

[11] Mammeri, Y. (2020). A reaction-diffusion system to better comprehend the unlockdown: Application of SEIR-type model with diffusion to the spatial spread of COVID-19 in France. Computational and Mathematical Biophysics, 8(1), 102-113.

[12] Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. Journal of the society for Industrial and Applied Mathematics, 11(2), 431-441.

[13] Medlock, J., & Meyers, L. A. (2009). Optimizing allocation for a delayed influenza vaccination campaign. PLoS currents, 1.

[14] Moore, S., Hill, E. M., Tildesley, M. J., Dyson, L., & Keeling, M. J. (2021). Vaccination and non-pharmaceutical interventions for COVID-19: a mathematical modelling study. The Lancet Infectious Diseases.

[15] Peirlinck, M., Linka, K., Costabal, F. S., Bhattacharya, J., Bendavid, E., Ioannidis, J. P., & Kuhl, E. (2020). Visualizing the invisible: The effect of asymptomatic transmission on the outbreak dynamics of COVID-19. Computer Methods in Applied Mechanics and Engineering, 372, 113410.

[16] Pesco, P., Bergero, P., Fabricius, G., & Hozbor, D. (2015). Mathematical modeling of delayed pertussis vaccination in infants. Vaccine, 33(41), 5475-5480.

[17] Viguerie, A., Veneziani, A., Lorenzo, G., Baroli, D., Aretz-Nellesen, N., Patton, A., ... & Auricchio, F. (2020). Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study. Computational Mechanics, 66(5), 1131-1152.

[18] Zhu N, Zhang D, Wang W et al. A novel coronavirus from patients with pneumonia in China, 2019. N Engl J Med 2020; 382:727-33.
References

[1] Acuña-Zegarra, M. A., Díaz-Infante, S., Baca-Carrasco, D., & Liceaga, D. O. (2021). COVID-19 optimal vaccination policies: a modeling study on effectiveness, natural and vaccine-induced immunity responses. Mathematical Biosciences, 108614.

[2] Ali, M., Shah, S. T. H., Imran, M., & Khan, A. (2020). The role of asymptomatic class, quarantine and isolation in the transmission of COVID-19. Journal of biological dynamics, 14(1), 389-408.

[3] Arcede, J. P., Caga-anan, R. L., Mentuda, C. Q., & Mammeri, Y. (2020). Accounting for Symptomatic and Asymptomatic in a SEIR-type model of COVID-19. Mathematical Modelling of Natural Phenomena, 15, 34.

[4] Buhat, C. A. H., Duero, J. C. C., Felix, E. F. O., Rabajante, J. F., & Mamplata, J. B. (2021). Optimal allocation of COVID-19 test kits among accredited testing centers in the Philippines. Journal of healthcare informatics research, 5(1), 54-69.

[5] Diekmann, J. Heesterbeek, & J. Metz, On the definition and the computation of the basic reproductive ratio $r_0$ in models for infectious diseases in heterogeneous populations, *J Math Biol* **28** (1990), 365-383.

[6] V. Driessche, & J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, *Math Biosci*, **180** (2002), 29-48.

[7] Jentsch, P. C., Anand, M., & Bauch, C. T. (2021). Prioritising COVID-19 vaccination in changing social and epidemiological landscapes: a mathematical modelling study. The Lancet Infectious Diseases.

[8] S. A. Lauer, K. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R. Meredith, A. S. Azman, N. G. Reich, and J. Lessler. The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmed cases: Estimation and application. Annals of Internal Medicine, 03 2020. URL https://doi.org/10.7326/M20-0504.

[9] Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares. Quarterly of applied mathematics, 2(2), 164-168.

[10] Macalisang, J. M., Caay, M. L., Arcede, J. P., & Caga-anan, R. L. (2020). Optimal Control for a COVID-19 Model Accounting for Symptomatic and Asymptomatic. Computational and Mathematical Biophysics, 8(1), 168-179.

[11] Mammeri, Y. (2020). A reaction-diffusion system to better comprehend the unlockdown: Application of SEIR-type model with diffusion to the spatial spread of COVID-19 in France. Computational and Mathematical Biophysics, 8(1), 102-113.

[12] Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. Journal of the society for Industrial and Applied Mathematics, 11(2), 431-441.

[13] Medlock, J., & Meyers, L. A. (2009). Optimizing allocation for a delayed influenza vaccination campaign. PLoS currents, 1.

[14] Moore, S., Hill, E. M., Tildesley, M. J., Dyson, L., & Keeling, M. J. (2021). Vaccination and non-pharmaceutical interventions for COVID-19: a mathematical modelling study. The Lancet Infectious Diseases.

[15] Peirlinck, M., Linka, K., Costabal, F. S., Bhattacharya, J., Bendavid, E., Ioannidis, J. P., & Kuhl, E. (2020). Visualizing the invisible: The effect of asymptomatic transmission on the outbreak dynamics of COVID-19. Computer Methods in Applied Mechanics and Engineering, 372, 113410.

[16] Pesco, P., Bergero, P., Fabricius, G., & Hozbor, D. (2015). Mathematical modeling of delayed pertussis vaccination in infants. Vaccine, 33(41), 5475-5480.

[17] Viguerie, A., Veneziani, A., Lorenzo, G., Baroli, D., Aretz-Nellesen, N., Patton, A., ... & Auricchio, F. (2020). Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study. Computational Mechanics, 66(5), 1131-1152.

[18] Zhu N, Zhang D, Wang W et al. A novel coronavirus from patients with pneumonia in China, 2019. N Engl J Med 2020; 382:727-33.