A mathematical model for policy of vaccinating recovered people in controlling the spread of COVID-19 outbreak

1.   Introduction

On March 9th, 2020, the World Health Organization (WHO) declared COVID-19 a global pandemic. The virus can be transmitted from person to person through droplets that are expelled when an infected individual coughs or sneezes. These droplets can contaminate objects or surfaces and can cause infection if an uninfected individual touches their eyes, nose, or mouth after touching the contaminated object or surface. The first case of COVID-19 was reported in Wuhan, China in December 2019 and had spread to 230 countries by November 2020 ^[1]. In Indonesia, the first case of COVID-19 was reported on March 2nd, 2020. As of November 15th, 2022, Indonesia has recorded 6,565,912 confirmed cases of COVID-19, with 6,356,794 people recovered and 159,158 patients who have died ^[2].

The Indonesian government has made various efforts to reduce the number of positive COVID-19 cases, including requiring the use of masks, implementing large-scale social restriction policies, and controlling transportation. Additionally, the government has administered various types of vaccines, including Sinovac, Astra Zeneca, Pfizer, Moderna, Janssen, and Sinopharm, to the community. The government also continues to boost Testing, Tracing, and Treatment (3T) and encourages the use of the PeduliLindungi application for digital tracing. Digital tracing is the government's effort to identify and detect individuals through tracking location data and information digitally. This application is connected and integrated with existing systems and databases at the Ministry of Health, and by using a QR code, it can track data on individuals who have been vaccinated and their test results. Tracing results with the PeduliLindungi application will make it easier for people to receive treatment or necessary handling if needed.

Several researchers in mathematics have modeled the spread of COVID-19 in Indonesia. Nuraini et al. used Richard's curve and the least squares method to represent the dynamics of COVID-19 patients in Indonesia ^[3]. Suwardi utilized the SEIR (susceptible, exposed, infectious, and recovery) model to describe the dynamics of COVID-19 spread in Indonesia ^[4]. Susanto used the SI (susceptible and infectious) model to depict the spread of COVID-19, assuming that there is a lockdown but citizens perform their activities as usual ^[5]. Additionally, Susanto et al. developed their research using the SIQRD (susceptible, infected, quarantine, recovery, and death) model to describe the administration of vaccines to a class of vulnerable populations ^[6]. Mukandavire, Nyabadza, and Malunguza, used the SEIR model to quantify early COVID-19 outbreak transmission in South Africa and explore vaccine efficacy scenarios ^[7]. Olivares and Staffetti used the SEI$_s$IAQR (susceptible, exposed, infected symptomatic, infected asymptomatic, quarantine, and recovery) model for uncertainty quantification of a mathematical model of COVID-19 transmission dynamics with a mass vaccination strategy. Model of spread and control of CCOVID-19 in China is studied in ^[8]. For a more general epidemic model, ^[9] provide a stochastic bifurcation analysis and stochastic delayed optimal control. Finally, Diagne used the SVEAIHR (susceptible, vaccinated, exposed, infected asymptomatic, infected symptomatic, hospitalized, and recovery) model to describe the administration of vaccines to a class of vulnerable populations ^[10]. Thus, researchers have developed various mathematical models to study the spread of COVID-19 with vaccines.

This paper is based on a previous study by Diagne ^[10] on constructing a model for the spread of COVID-19. The researcher added a compartment for recovered individuals who have been vaccinated to the model. In the previous research, infected individuals with no symptoms, infected individuals with symptoms, and quarantined individuals were assumed to be the ones who spread the virus. In this study, the model was modified to assume that the individuals who spread the virus are those without symptoms. Additionally, vulnerable individuals, recovered individuals, and vaccinated recovered individuals can be infected again if they come into close contact with infected individuals without symptoms. The study determined the basic reproduction number, which is a threshold parameter that indicates whether a disease is spreading or not. The next generation matrix method ^[11,12] was used to compute the basic reproduction number. The paper then conducted a sensitivity analysis on the basic reproduction number to identify the parameters that have the most significant impact on the spread of COVID-19. The results of this analysis can be used to inform policies regarding the control of the COVID-19 outbreak. The structure of this paper includes sections on the mathematical model and formulation, equilibrium point, basic reproduction number, numerical simulation, and sensitivity analysis of the basic reproduction number.

2.   Model formulation

The COVID-19-vaccine mathematical model presented in this study includes eight population classes. The vulnerable population class ($S$) represents individuals who are susceptible to COVID-19 infection. The population class ($V_s$) represents individuals who have been vaccinated from the vulnerable population. The population class ($E$) represents individuals who have been exposed to COVID-19. The asymptomatic infected population class ($A$) represents individuals who are infected with COVID-19 but do not exhibit symptoms. The infected population class with symptoms ($I$) represents individuals who are infected with COVID-19 and exhibit symptoms. The hospitalized population class ($H$) represents individuals who are infected with COVID-19 and have been hospitalized or are in quarantine. The recovered population class ($R$) represents individuals who have recovered from COVID-19. The vaccinated and recovered population class ($V_R$) represents individuals who have been vaccinated and have recovered from COVID-19. The total population is represented by $N$, where

In this model, the birth rate and death rate are different, which causes the population to be dynamic. Vulnerable individuals are assumed to be infected through contact with undetected infected individuals.

The assumptions in constructing the model of vaccination policy in controlling the spread of COVID-19 are as follow:

$(1)$ Infectious disease through a direct contact with infected individuals who show no symptoms. Infected individuals who show symptoms do not have a direct contact with vulnerable individuals and will perform self-quarantine or be hospitalized, so the individuals with symptoms will not spread the virus.

$(2)$ Individuals in the population are assumed to be not infected by COVID-19, therefore the increase of individuals in vulnerable population are from migration and birth rate.

$(3)$ The illness may cause death. The decrease caused by death from the illness happen in infected class $(I)$ and hospitalized $(H)$.

$(4)$ Population in hospitalized $(H)$ class do not interact with vulnerable and recovered population, therefore the spread of COVID-19 does not happen in hospitalized $(H)$ population.

$(5)$ The population who are exposed to COVID-19 will be the asymptomatic and symptomatic population.

$(6)$ The reduction in individuals in the population is caused by natural death, with a constant natural death rate and the same value for each class of the population.

$(7)$ The population in individuals who have been vaccinated and in the individual population who have recovered can be re-infected if interacting or making contact with individuals in the infected population is not detected.

$(8)$ Infected individuals with symptoms requiring more intensive care are admitted to the hospitalized$(H)$ class

The parameters used in the vaccination policy in controlling the spread of COVID-19 are shown in Table 1:

The transfer diagram model represents the spread of the Corona virus. The differential equation for the spread of the Corona virus model from the diagram is:

The solution to Eq (2.1) is in the domain:

Theorem 1. The solution to the system of Eq (2.1) is within the domain $\mathcal{D}$ in $\mathbb{R}^8$.

Proof. Consider the total population

We can derive the following equation by taking the derivative of $N$ with respect to time:

Integrating Eq (2.3), we obtain

By taking the limit as $t$ approaches infinity, we have

This implies that there is a maximum limit of $\dfrac{\Lambda}{\mu}$ for the population. Therefore, the solution to Eq (2.1) lies within the domain $\mathcal{D}$. □

3.   Model analysis

This part will discuss about equilibrium point, basic reproduction number, and stability analysis.

3.1.   Equilibrium point

The disease-free equilibrium point is a point where there is no disease in a population. In order to find the disease-free equilibrium point, we evaluate the system of Eq (2.1) when

Theorem 2. COVID-19 outbreak model (2.1) has only one equilibrium point

which represents the disease-free equilibrium point. At this point, there is no disease in the population.

Proof. There is only one equilibrium point for the COVID-19 outbreak model (2.1), which is the disease-free equilibrium point $x_{dfe}$. This point corresponds to the situation where there is no disease in the population. By substituting the values from $x_{dfe}$ to Eq (2.1), all derivatives become zero, which confirms that $x_{dfe}$ is indeed an equilibrium point for the model. Additionally, it can be concluded that the equilibrium point on $\mathcal{D}$ fulfills $E = A = I = H = 0$ condition, which is the condition for a disease-free equilibrium point. □

3.2.   Basic reproduction number

To understand the disease or infection spread level, the parameter to see how big is the potential of a disease spreading in a population is needed. The parameter is called basic reproduction number, whic is notated by $R_0$. The basic reproduction number $(R_0)$ is defines as the average sum of the secondary case caused by an infected individual during their infection span in the vulnerable population. The method used in obtaining the basic reproduction number is next generation matrix. Suppose the matrix $\mathcal{F}$ represents the occurence of new infection rate, matrix $\mathcal{V^{-}}$ represents the mobility of individuals who migrate from the first class to another class and matrix $\mathcal{V^{+}}$ represents the mobility of individuals who get into one class from another class. Equation (2.1) may be written into

with

defines matrix $\mathcal{F}$ and $\mathcal{V}$ with

and

The basic reproduction number is a spectral radius from next generation matrix $\mathcal{FV}^{-1}$. According to Eq (2.1), we get

and

therefore the basic reproduction number is

The basic reproduction number $(R_0)$ is defined as the number of new infections produced by an infected individual in a fully vulnerable population. The disease-free equilibrium point is asymptotically stable if $R_0 < 1$. When $R_0 < 1$, the number of infected individuals will decrease and, on average, an infected individual will infect less than one individual. This means that the infection cannot develop into a plague and the COVID-19 outbreak will eventually run out of the population. However, if $R_0 > 1$, the system will be unstable, leading to an increasing number of infected individuals. This enables the COVID-19 outbreak to develop and potentially cause significant harm to the population.

4.   Numerical simulation

This chapter will discuss the numerical simulation of vaccination distribution related to COVID-19 outbreak. The parameter values used in the simulation are shown in Table 2 below:

The initial condition is given to find the disease-free equilibrium point and the endemic equilibrium point. The condition is $S(0) = 100.000, \ V_s(0) = 50.000, \ E(0) = 150, \ A(0) = 100,$ $I(0) = 150, \ H(0) = 45, \ R(0) = 259$ and $V_R(0) = 35.$ After that, we get the value of $\beta_{1} = 0, 000815, \ \beta_2 = 0, 00000059, \ \beta_3 = 0.00000045, \ \beta_4 = 0.000000011.$

Figure 1 shows the disease-free equilibrium point. Class S will decrease until it converges to the disease-free equilibrium point and each class in $E$, $A$, $I$, $H$, and $R$ classes will initially increase, then decreases until it converges to zero point. Based on that condition, the solution for each class with a given initial value and parameter value will reach its equilibrium point. Furthermore, the values for endemic equilibrium point are $\beta_{1} = 0, 00815, \ \beta_2 = 0, 0000059, \ \beta_3 = 0.0000045, \ \beta_4 = 0.000000011$ and we get:

An endemic equilibrium point is an equilibrium point with the existence of a disease. The endemic equilibrium point has a basic reproduction number value of $R_0 = 3.782 > 1$, which means the outbreak will always happen. Figure 2 shows that many individuals in class $S$ undergo a decrease as time goes by and will converge to the equilibrium point. In addition, many individuals in classes $E$, $A$, $I$, $H$, $R$ undergo an initial increase, followed by a decrease until they converge into the equilibrium point.

The initial simulation of vaccine distribution is conducted using the given initial conditions: $S(0) = 100,000, \ E(0) = 150, \ A(0) = 50, \ I(0) = 100, \ H(0) = 100, \ R(0) = 249, \ V_s(0) = 0$, and $V_{R}(0) = 0$. Based on the parameter values in Table 2, the simulation is carried out.

The blue color in Figures 3–5 represents the scenario where individuals are not vaccinated, the yellow color represents the vaccination of the vulnerable population in class $S$, and the purple color represents the vaccination of the vulnerable populations in classes $S$ and $R$. The simulation results show that providing vaccines to both classes $S$ and $R$ can lead to the elimination of the disease, while vaccinating only class $S$ is not sufficient to prevent an outbreak.

Figure 6 shows the effect of vaccine rate $v$ and vaccine effectiveness $\theta$ on the populations of classes $A, I,$ and $H$. From the figure, it can be concluded that the vaccine rate $v$ and vaccine effectiveness $\theta$ significantly affect the populations of classes $A, I,$ and $H$. When both vaccine rate $v$ and vaccine effectiveness $\theta$ are equal to zero, the population of individuals in class $A$ is 175, the population in class $I$ is 650, and the population in class $H$ is 620. As the vaccine rate $v$ and vaccine effectiveness $\theta$ increase, the populations of classes $A, I,$ and $H$ decrease. The results show that the vaccination of the vulnerable group is more effective in suppressing the increase in the populations of classes $A, I,$ and $H$ compared to vaccination of class $R$.

Sensitivity analysis is conducted to understand the most influential parameter of the COVID-19 outbreak. The sensitivity to $p$ is defined by the equation

using Eq (4.1) based on Table 2, the sensitivity value of a parameter is shown in Table 3:

According to Table 3, the influential parameters of the COVID-19 outbreak are $\beta_{1}$, $\beta_{3}$, and $\gamma$. If we increase the value of these parameters, the $R_{0}$ value will also increase. However, if we decrease the parameter values of $v$, $\mu$, and $\alpha_1$, the value of $R_{0}$ will decrease. When $R_{0} > 1$, the parameters that positively affect the COVID-19 outbreak are $\beta_{1}$ and $\beta_{2}$. If the values of these parameters are decreased, it will reduce the value of $R_{0}$. Furthermore, a simulation of the connection between $\beta_{1}$ and $\beta_{2}$ to $R_0$ is conducted by modifying the parameter $v$ as follows:

Figure 7 shows the relationship between $\beta_{1}$ and $v$ parameters in the spread of the COVID-19 outbreak, where $R_0 = 1$. Eq (3.1) is used to form the function of $\beta_{1}$ in terms of $v$. Based on the parameters in Table 2, we obtain a linear function of $\beta_{1}$ in terms of $v$, which is $\beta_{1} = 0.000815-0.00000013v$. The blue line represents the linear function of $\beta_{1}$ in terms of $v$, and the red line represents the value of $R_0 = 1$. The intersection point between the blue and red line is the critical point, which is $v\approx0.023$. If the value of $v$ is greater than $0.023$, then the value of $R_0$ will be less than $1$, and the outbreak will not happen. However, if the value of $v$ is less than $0.023$, then the outbreak will occur. This simulation shows that the parameter $v$ plays a significant role in controlling the COVID-19 outbreak.

From Figure 8, we get that the area above the graph is an area with a basic reproduction number value of more than one, in other words, the area is an area of disease existence, while the area below the graph is an area with a basic reproduction number value of less than one, in other words, the area is a disease-free area. From the figure above it is also obtained that even though the vaccination rate is large, when the contact rate with an infected individual is large, the individual will be infected, for example from the figure the vaccine rate is $0.9$ but if the contact rate is $3.5 \times 10^{-3}$ then the spread of the disease in a population will still exist.

5.   Conclusions

In this paper, a mathematical model $(SV_s EIAHR V_R)$ for the spread of COVID-19 is constructed. The results show that vaccine distribution can suppress the spread of COVID-19. If the vaccine is distributed only once to the $S$ class, COVID-19 outbreak may still occur in the population. However, if the vaccine is distributed more than once, or distributed to the $R$ class, COVID-19 outbreak can be suppressed and the population can become disease-free. According to the sensitivity analysis, one of the ways to suppress COVID-19 is by increasing the $v$ parameter, which means increasing vaccine distribution, and the $\alpha_1$ parameter, which means accelerating the detection of undetected infected individuals who have recovered. In addition, decreasing the contact between vulnerable individuals and infected individuals will decrease the $\beta_{1}$ parameter, causing $R_0 < 1$. This means that an infected individual has the potential to infect less than one individual, leading to a disease-free population. Based on the sensitivity analysis, for a total susceptible population of 100,000 people, with 150 exposed, 50 undetected infected, 100 detected infected, 100 quarantined, and 259 recovered individuals, and based on the parameter values in Table 2, if the value of $v$ is set to 0.65 and $\beta_{1} < 3.8 \times 10^{-3}$ and $\beta_{2} < 0.25 \times 10^{-6}$, then the population will be disease-free.

Acknowledgments

This research was supported by the project Rekognisi Tugas Akhir Universitas Gadjah Mada (RTA-UGM) 1525/UN1/DITLIT/Dit-Lit/PT.01.05/2022.

Conflict of interest

The authors declare that they have no conflicts of interest.