A deterministic transmission model for analytics-driven optimization of COVID-19 post-pandemic vaccination and quarantine strategies

1.   Introduction

The coronavirus disease of 2019 (COVID-19) is a disease caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus and was declared a pandemic by The World Health Organization (WHO) on March 11, 2020 ^[1]. Symptoms of COVID-19 include fever, cough, shortness of breath, fatigue, headache, nausea, vomiting, sore throat, and diarrhea. These symptoms typically appear about two to 14 days after the infection ^[2]. As of July 5, 2023, a total of 6,811,818 people in Indonesia have been reported infected with this virus, and 161,867 people have died from the infection ^[3]. Even though Indonesia has transitioned from the pandemic phase to the endemic phase of COVID-19 ^[4], it is important to note that the symptoms of COVID-19 can be exacerbated by the presence of comorbidities or underlying health conditions ^[5], and there are still reported cases of mortality due to COVID-19 ^[6]. Therefore, it remains crucial to continue practicing preventive measures against COVID-19.

Since the advent of COVID-19, which began affecting various parts of the world in late 2019, numerous scientific approaches have been implemented to comprehend the spread and control of the virus. These approaches encompass medical and public health ^[7,8] and social ^[9,10], mathematical ^{[11,12,13,14,15]}, and many other perspectives. In the field of mathematics, the approaches employed by the authors are quite diverse. These range from ordinary differential equations (ODEs) ^{[12,13,14,15,16]} to fractional differential equations ^[11,17,18], partial differential equations ^[19,20], stochastic models ^[21,22], and the utilization of big data analysis approaches ^[23,24].

Mathematical modeling is a crucial tool for examining and understanding the dynamics of COVID-19 transmission, particularly in the context of specific interventions. The most popular intervention in the early spread of COVID-19 is the implementation of social distancing. Some authors have introduced their mathematical model to analyze the impact of social distancing on controlling COVID-19. Aldila et al. ^[25] developed a model to assess the impact of social distancing policies in the city of Jakarta, Indonesia. They found that the intervention can reduce the spread of COVID-19, but the outbreak of COVID-19 will arise soon when social distancing is relaxed. Another model introduced by Gevertz et al. ^[26] in their paper considers nonsocially distant individuals. They found that implementing relaxation policies at regular intervals could effectively reduce the peak infectious load. However, it is crucial to note that this schedule is highly sensitive to both parameter values and the frequency of the schedule. Recently, Arik et al. ^[27] used a fractional-order differential equation to analyze the impact of the environment and social distancing on the spread of COVID-19.

Another intervention that many authors considered in their model is the vaccination strategy ^{[11,13,14,15]}. Saharan and Tee ^[28] used a simple susceptible-exposed-infectious-recovered (SEIR) model to analyze the effectiveness of a vaccine in controlling the spread of COVID-19 in Malaysia and Thailand. They used the Holt-Winters method to forecast their model dynamics. A fractional-order model was introduced by the authors in ^[29] to examine the impact of vaccines on the spread of COVID-19. A diffusion framework was used in their model to analyze the spatial spread of COVID-19. They found from their numerical experiments that it is imperative to implement strategic measures to ensure the timely and efficient distribution of vaccines. A mathematical model was introduced by the author in ^[30] to analyze the optimal allocation of COVID-19 vaccines in the Philippines. Using their model and data from the Philippines, they found that it is crucial to prioritize the primary vaccination series for the pediatric population (aged 5 to 11) within the first three months. Additionally, it is advisable to administer initial booster shots to individuals aged 12 and older during this timeframe. Please see ^{[12,31,32,33]} for more examples of the implementation of a mathematical model to assess the impact of vaccines on the spread of COVID-19.

Other models aim to demonstrate the significance of quarantine in suppressing the spread of COVID-19 ^{[34,35,36,37,38]}. These models emphasize the complementary role of quarantine, alongside vaccination, in effectively mitigating the transmission of the virus. Authors in ^[34,35] have highlighted that while quarantine alone can reduce the spread of COVID-19 in the long-term, it cannot completely eradicate the virus. Nevertheless, the quarantine program can be combined with other initiatives, such as vaccination ^[36,37] and interventions to enhance population awareness ^[38].

Based on this explanation, we can deduce that the long-term elimination of COVID-19 cannot be accomplished solely through vaccination or quarantine programs when executed in isolation from other factors. An essential element that can enhance the effectiveness of both vaccination and quarantine efforts is heightened public awareness. Consequently, certain models attempted to underscore the role of public awareness in the spread of COVID-19 ^{[38,39,40,41]}. Collectively, these models highlight that an increase in public awareness can effectively control the spread of COVID-19.

Hospital care and capacity play a crucial role in controlling COVID-19. Adequate hospital resources, such as beds, ventilators, medical staff, and equipment, are essential for managing and treating severe cases of COVID-19. Hospitals provide care for individuals with more severe symptoms, and having sufficient capacity is crucial to ensuring that everyone who needs medical attention can receive it promptly. Adequate and timely medical treatment provided in hospitals plays a crucial role in reducing the impact and transmission of diseases, as evidenced in countries such as the USA, Brazil, India ^[42] and Indonesia ^[43]. The availability and proper utilization of resources, such as masks and ventilators, also significantly contribute to the recovery of COVID-19 patients ^[44]. Moreover, several models highlight that hospitalization, when combined with other interventions like vaccination, can positively contribute to COVID-19 recovery ^[12]. Some models also indicate limited resources in treatment using a saturated function, leading to the conclusion that treatment resources are an important factor in COVID-19 recovery ^{[45,46,47,48,49]}.

From the above explanation, it can be broadly seen that maximizing the intensity of vaccine administration, implementing quarantine measures, enforcing social restrictions, and improving the quality of hospitals will increase the possibility of eliminating COVID-19 from the community. However, this comes with the consequence of the high cost associated with the required interventions. Therefore, an optimal intervention is needed to minimize the spread of COVID-19 in the population. One approach that can be taken is to tailor COVID-19 interventions based on field needs—in this case, depending on time. One mathematical approach that can be used is to construct the problem as an optimal control problem. The basic idea of this approach is to treat the intervention $u$ as a time-dependent variable $u(t)$. The aim of this implementation is to minimize the number of infected individuals at the lowest possible cost for interventions. Many authors have introduced mathematical models with optimal control for COVID-19 interventions. The authors in ^[50] developed an optimal control model to analyze the impact of virus protection on the dynamics of COVID-19. The optimal control model for how COVID-19 spread, considering lockdown measures and asymptomatic cases, is discussed in ^[51]. They found that implementing a lockdown is essential in the early stages of an outbreak. On the other hand, combining vaccination and treatment is effective in controlling and preventing the spread of COVID-19. A more complex model, introduced by the authors in ^[52], accommodates the low-high-risk population and three types of interventions, namely isolation, detection, and treatment. Their parameter values were estimated using infection incidence data from Shanghai (March–April 2022). They found that when the budget is limited, implementing isolation and detection is more effective. On the other hand, if the budget is sufficient, implementing treatment alone is more effective compared to other possible scenarios. For more examples of optimal strategies for COVID-19 control, refer to ^[53,54].

From the above explanation, we understand that although many references discuss the implementation of vaccination and quarantine as strategies to control COVID-19, there is still a lack of papers that consider the combined impact of population awareness on the effectiveness of these interventions. Therefore, the objective of this paper is to develop a deterministic model that incorporates vaccination and quarantine compartments while also considering the effects of saturation functions that represent awareness and resource limitations. Given that individuals who have contracted COVID-19 develop immunity to the disease for at least six months ^[55], this paper considers individuals with prior COVID-19 infections to be immune to the disease. An optimal control approach is also considered in our model to assess the optimal scenario for suppressing the spread of COVID-19 with a low cost of intervention through a combination of vaccination and quarantine.

This paper is organized as follows: The first section discusses the state of the art of the research. In Section 2, we encompass the construction of the model, where we delve into the assumptions employed, the differential equation system, and the transmission diagram. We engage in a comprehensive discussion on model analysis, which includes positivity analysis and the parameterization of model parameters in Section 3. In Section 4, we construct the nondimensionalized version of our model to examine the COVID-19-free equilibrium points, the endemic equilibrium, local stability analysis, and bifurcation analysis. Section 5 discusses a series of numerical experiments incorporating elasticity analysis, sensitivity analysis, and autonomous simulation. We define our optimal control model in Section 6. Several numerical experiments for various scenarios are conducted in this section. Lastly, some discussion and conclusions are provided in Section 7.

2.   Model construction

Before constructing the model, it is necessary to determine the assumptions or limitations of the model. In this paper, the following assumptions are made: We assume that the total population remains constant, there are no deaths caused by COVID-19, there is no waning immunity, and vaccination does not provide perfect protection against infection for vaccinated individuals. We assume that a vaccinated individual has better protection against contact with other individuals. Furthermore, we include medical interventions to accelerate the recovery rate of infected individuals. The human population is divided into five compartments based on their health status, namely susceptible $(S)$, vaccinated $(V)$, infected $(I)$, quarantined $(Q)$, and recovered $(R)$.

The model construction process is based on the transmission diagram shown in Figure 1. The modeling process is as follows: In our model, we incorporate the intervention of vaccination with a specific rate denoted as $u_1$. Under this intervention, vaccinated individuals are segregated into a separate class, denoted as $V$. It is important to note that, based on many reports ^[12], COVID-19 vaccines do not offer complete protection against infection. Therefore, vaccinated individuals can still potentially contract COVID-19, albeit at a reduced rate, represented by the parameter $\beta_2$. This rate is smaller than the infection rate for susceptible individuals in the $S$ class $(\beta_1)$. In several reports ^[56,57,58], it has been found that individuals who have received vaccination tend to exhibit a heightened awareness toward the spread of COVID-19. This is attributed to their enhanced knowledge of the dangers of COVID-19, whether acquired through media, social interactions, or other means. An increase in the number of COVID-19 infections tends to heighten public awareness of the disease spread. As the number of infected individuals rises, there is often greater media coverage, discussions, and public discourse surrounding the disease. This increased attention can contribute to a heightened awareness among the general population about the importance of preventive measures, the severity of the situation, and the potential risks associated with the disease ^[59]. Hence, we consider that the probability of infection should be a monotonically decreasing function of $I$. Hence, the infection term within the vaccinated class is modeled by a saturated function expressed as $\frac{{{\beta _2}}}{V}I1 + {b_1}I$.

Additionally, we assume that individuals who become infected can be detected at a rate denoted as $u_2$. Once detected, individuals receive medical treatment either in a hospital setting or are monitored from home, which we categorize as a quarantine intervention. As a result, their recovery rate increases. Drawing from various reports ^[60,61], it has been observed that the escalating number of COVID-19 infection cases during the pandemic may disrupt the quality of hospital services. This phenomenon is attributed, among other factors, such as hospital staff burnout and the presence of resource constraints, as well as challenges in adaptation and mitigation. Public health authorities and healthcare providers frequently collaborate to deploy strategic interventions aimed at effectively managing case surges and mitigating the adverse effects on the standard of healthcare. Based on these facts, we assumed that the recovery rate of quarantined individuals could become saturated due to the high number of individuals in quarantine in the hospital. Therefore, we model the recovery rate of quarantined individuals as a a saturated function $\gamma_0 + \dfrac{\gamma_1}{1 + b_2 Q}$. With this definition of the recovery function, we can see that the recovery rate of quarantined individuals reaches its maximum when $Q$ is relatively small ($Q \approx 0$), i.e., $\gamma_0+\gamma_1$. On the other hand, when the number of infected individuals is approximately large, the recovery rate for quarantined individuals reaches its minimum level, i.e., $\gamma_0$.

Based on the assumptions outlined above, the mathematical model for the spread of COVID-19 is formulated as follows:

with the parameter description given in Table 1. Since the total population $N$ is $S+V+I+Q+R$, then we have that the dynamic of total population is given by:

If the total of population is approximately constant, then we can approximate the value of $\Lambda_h$ by assuming $\frac{dN}{dt} = 0,$ which lead us to $\Lambda_h = \mu N$.

3.   Preliminary analysis

3.1.   Positive invariant region

In this section, a nonnegativity analysis of the solution obtained from Model (1) will be performed.

Theorem 1. Region $\Omega$ which is defined by:

is positively invariant with respect to Model (2.1).

Proof. Model (2.1) can be written as follows:

where:

while $d_1 = -u_1 - \mu - \beta_1 I$, $d_2 = -\dfrac{\beta_1}{1 + b_1 I} I - \mu$, $d_3 = \dfrac{\beta_2}{1 + b_1 I} I$, $d_4 = -u_2 - \mu -\gamma_0$, $d_5 = -\big(\gamma_0 + \dfrac{\gamma_1}{1 + b_2 Q}\big) - \mu$, and $d_6 = \gamma_0 + \dfrac{\gamma_1}{1 + b_2 Q}$. Due to the nonnegative entries outside the main diagonal, matrix $C$ is classified as a Metzler matrix. Since $D \geq 0$, then Model (2.1) possesses a positive invariant in $\mathbb{R}_+^{5}$ ^[67]. Hence, all trajectories of solution of Model (2.1), which originate from an initial state in $\mathbb{R}_+^5$, confine forever.

Summing up all equations in Model (2.1), we have:

Then the solution is $N(t) = \frac{\Lambda_h}{\mu} + \left(N(0) - \frac{\Lambda_h}{\mu}\right)e^{\mu t}.$ If $N(0) = \frac{\Lambda_h}{\mu}$, then we have $N(t) = \frac{\Lambda_h}{\mu}$ for all time $t$. If $N(0) < \frac{\Lambda_h}{\mu},$ then $N(t)$ will be monotonically increasing and tends to $\frac{\Lambda_h}{\mu}$ for $t\to \infty$. On the other hand, if $N(0) > \frac{\Lambda_h}{\mu},$ then $N(t)$ will be monotonically decreasing and tends to $\frac{\Lambda_h}{\mu}$ for $t\to \infty$. Therefore, all feasible solutions of the Model (2.1) enter the region $\Omega$, implying that the region is an attracting set.

3.2.   Parameter estimation

To estimate the parameters for Model (2.1), a data fitting process was conducted for the number of COVID-19 infections in West Java, Indonesia, from July 1, 2022, to October 1, 2022. In the context of our model, this data represents the dynamics of $I(t)$. Our aim is to fit the dynamic output of $I(t)$ from our model with the available data. Mathematically, this task can be expressed as minimizing the Euclidean distance between COVID-19 incidence data (denoted by $I^{\text{data}}$) and the solution of $I(t)$ from Model (2.1), using the best-fit parameters $\beta_1, \beta_2, u_1, u_2, b_1, b_2, \gamma_0, \gamma_1$, and the best initial conditions $S(0), V(0), I(0), Q(0), R(0).$ Hence, we define the following objective function:

where $T$ is the maximum number of existing data points. Our aim is to find the optimal parameters:

such that:

where $\Theta$ is the admissible value for the estimated parameters (fourth column of Table 1). The following constraints need to be mentioned:

● The infection rate of vaccinated individuals $(\beta_2)$ is smaller than that of susceptible individuals $(\beta_1)$. Hence, we have $\beta_2 < \beta_1$.

● Vaccination and quarantine is a non-stop intervention by the government. Hence, $u_1 > 0, u_2 > 0$.

● For biological interpretation purposes, we have $\gamma_0, \gamma_1, b_1, b_2,$ and all initial conditions are positive.

We approach our parameter estimation challenge as the task of minimizing a constrained, nonlinear, multivariable function. Specifically, our objective is to minimize $\mathcal{C}$ subject to the constraints imposed by the ODE model (2.1), while ensuring that the Model parameters and initial conditions adhere to biologically reasonable values. We use the fmincon toolbox in MATLAB to solve our problem. The result is shown in Figure 2. The best-fit parameter is provided in Table 1 (sixth column), while the best fit initial condition is as follows:

For more parameter estimation methods, readers may read ^{[68,69,70,71]} for further references.

From the results of parameter estimation in Table 1, we can see that the control reproduction number is still greater than one, which indicates that the number of infected individuals $(I(t))$ will continue to exist as $t \rightarrow \infty$ (see Section 4 for further analytical results on the COVID-19 endemic equilibrium point). We found that $u_1 = 0.0344$ and $u_2 = 0.1540$, which resulted in the findings shown in Figure 2. If the control interventions were stopped $(u_1 = u_2 = 0)$, then we would obtain $\mathcal{R}_0 = 9.316$. Any relaxation of the control variables will increase $\mathcal{R}_c$. On the other hand, if we increase the intensity of vaccination and quarantine, then $\mathcal{R}_c$ can be reduced. Figure 3 shows the impact of relaxation/improvement of vaccination and quarantine interventions on the dynamics of infected individuals $I$.

4.   Model analysis

4.1.   Nondimensionalization of the model

Let us assume the total population is constant. Hence we have:

which represents the proportion of each class to the total human population. Assuming $\tau = t \gamma_0$, $m = \dfrac{\mu}{\gamma_0}, a_1 = \dfrac{\beta_1 N}{\gamma_0}, a_2 = \dfrac{\beta_2}{\gamma_0}, g_1 = \dfrac{\gamma_1}{\gamma_0}, c_1 = \dfrac{b_1}{N}, c_2 = \dfrac{b_2}{N}, e_1 = \dfrac{u_1}{\gamma_0}, e_2 = \dfrac{u_2}{\gamma_0}$, and reducing the susceptible compartments, a nondimensionalized version of the Model (2.1) is obtained as follows:

where $x_5 = 1-x_1 - x_2 - x_3 - x_4.$

4.2.   COVID-19-free equilibrium and the control reproduction number

In this section, we initiate our analysis by determining the form of the COVID-19-free equilibrium, followed by the calculation of the control reproduction number, denoted as $\mathcal{R}_c$. Reproduction numbers have been utilized by various authors in epidemiological models ^[11,49] to assess the dynamic behavior of the proposed model, including the existence and stability of equilibrium points.

Theorem 2. The COVID-19-free equilibrium of system (4.1) is given by:

which always exists unconditionally.

Proof. The equilibrium points associated with Model (1) are obtained when the population changes in each compartment remain constant, namely, when $\dfrac{dx_1}{d\tau} = 0, \dfrac{dx_2}{d\tau} = 0, \dfrac{dx_3}{d\tau} = 0,$ and $\dfrac{dx_4}{d\tau} = 0.$ As the COVID-19-free equilibrium point represents a condition without COVID-19 in the population, we have $x_3^\dagger = 0$. By direct calculation, substituting $x_3^\dagger = 0$ into $\dfrac{dx_1}{dt}, \dfrac{dx_2}{dt},$ and $\dfrac{dx_4}{dt}$, we obtain $x_1^\dagger = \dfrac{m}{m + e_1}$, $x_2^\dagger = \dfrac{e_1}{m + e_1}$, and $x_4^\dagger = 0$. From the expressions of $x_1^\dagger, x_2^\dagger, x_3^\dagger,$ and $x_4^\dagger$, it is clear that they always exist in $\mathbb{R}_{\geq 0}^4$ without any existence criteria. Hence, the proof is complete.

From the direct consequences of the expression of $\mathcal{E}_0$, we have the following corollary.

Corollary 1. Without the intervention of vaccination, the COVID-19-free equilibrium of system (4.1) is given by:

Corollary 1. Without the intervention of vaccination, the COVID-19-free equilibrium of system (4.1) is given by:

Theorem 4.2. The control reproduction number from system (4.1) is $\mathcal{R}_c = \dfrac{a_{1} m +a_{2} e_{1}}{\left(e_{2}+1+m \right) \left(m +e_{1}\right) }$.

Proof. We use a next-generation matrix approach ^[72] to determine the control reproduction number of the system (4.1). At first, we determine the Jacobian matrix of the subsystem of (4.1) that only involves the infected compartments $x_3$ and $x_4$. The Jacobian matrix evaluated at $\mathcal{E}_0$ is given by:

The matrix $J$ can be written as the sum of the transmission matrix $T$ and the transition matrix $\Sigma$, where:

Since the second row of matrix T consists entirely of zeros, the next-generation matrix can be derived as:

where ${\bf E} = \left[\begin{array}{c} 1 \\ 0 \end{array}\right]$ and $\textbf{E}'$ is the transpose of $\textbf{E}$. Since the spectral radius of the matrix above is given by $\dfrac{a_{1} m +a_{2} e_{1}}{\left(e_{2}+1+m \right) \left(m +e_{1}\right)}$, we can conclude that $\mathcal{R}_c = \dfrac{a_{1} m +a_{2} e_{1}}{\left(e_{2}+1+m \right) \left(m +e_{1}\right)}$. Hence, the proof is complete.

4.3.   The existence of COVID-19 endemic equilibrium point

Theorem 4. Let $\mathcal{R}_c^*$ be $\mathcal{R}_c$ that satisfies $k_1^2-4k_0k_2 = 0$ and $\mathcal{K} = \left(e_{2}+1+m \right)^{2} \left(m +e_{1}\right)-a_{2} e_{1} \left(e_{2}+1+m \right) +a_{2} e_{1} \left(m c_{1}+a_{2}\right)$, then, system (4.1) has:

1) One endemic equilibrium if $\mathcal{R}_c > 1$,

2) No endemic equilibrium for $\mathcal{R}_c < 1$ if $\mathcal{K}\geq 0$,

3) No endemic equilibrium for $\mathcal{R}_c < \mathcal{R}_c^* < 1$ if $\mathcal{K} < 0$,

4) Two endemic equilibrium for $\mathcal{R}_c^* \leq \mathcal{R}_c < 1$ if $\mathcal{K} < 0$.

Proof. The endemic equilibrium represents a condition in which the disease persists within the system ^[73], specifically at the point $(x_1, x_2, x_3, x_4)$ with $x_3 \neq 0$. As there is no change in the system when the endemic equilibrium point is reached, we set $\dfrac{dx_1}{d\tau} = 0, \dfrac{dx_2}{d\tau} = 0, \dfrac{dx_3}{d\tau} = 0,$ and $\dfrac{dx_4}{d\tau} = 0$ with $x_3 \neq 0$. By direct calculation, we can find:

where $x_1^* = \dfrac{m}{a_{1} x_{3}^*+m +e_{1}}, x_2^* = \dfrac{e_{1} m \left(c_{1} x_{3}^*+1 \right)}{\left(a_{1} x_{3}^*+m +e_{1}\right) \left(m c_{1} x_{3}^*+m +a_{2} x_{3}^*\right)}$, $x_4^*$ is the positive root of $e_{2} x_{3}^*-\left(1+\dfrac{g_{1}}{c_{1} x_{4}^*+1}\right) x_{4}^*-m x_{4}^* = 0$, and $x_3$* is taken from the positive roots of the following polynomial:

where:

To analyze the existence of the endemic equilibrium, it is necessary to guarantee the existence of a positive root in the polynomial (4.4). If the roots are positive, then $x_1^*, x_2^*$, and $x_4^*$ will also be positive. From Descartes' rule of signs, it can be concluded that when $\mathcal{R}_c > 1$, the system will have a unique endemic equilibrium point. Furthermore, for $k_0 < 0$, it is possible for the system to have zero or two positive roots. Hence, further analysis is needed. Gradient analysis is conducted at $\mathcal{R}_c = 1, x_3 = 0.$ Using implicit derivative on polynomial (4.4), we have:

Hence, if $\mathcal{K} = \left(e_{2}+1+m \right)^{2} \left(m +e_{1}\right) - \left(e_{2}+1+m \right)a_{2} e_{1} +[a_{2} e_{1} \left(m c_{1}+a_{2}\right)] < 0$, then we will have a negative gradient of $x_3$ at $\mathscr{R}_c = 1, x_3 = 0$. This indicates the existence of positive roots $x_3$ for some intervals when $\mathcal{R}_c < 1$. Since $f(x_3^*)$ is a second-degree polynomial, then $f(x_3^*)$ will have a turning point at $\mathcal{R}_c = \mathcal{R}_c^*$, where $\mathcal{R}_c^*$ is the value of $\mathcal{R}_c$ that satisfies the discriminant of $f(x_3^*)$ equal to zero, which is $k_1^2-4k_0 k_2 = 0$. Hence, we have proved case 3 of the theorem. If $\mathcal{R}_c < \mathcal{R}_c^*$, then we will have no positive root (case 4 of the theorem proved). Lastly, if $\mathcal{K}\geq 0$, then we will have a positive gradient, which confirms case 2 in the proven theorem. Furthermore, we use Descartes' rule of signs to determine the possible number of positive roots of the polynomial $f(x_3)$. The results are shown in Table 2. It can be observed that $f(x_3)$ may have 0, 1, or 2 positive roots. Additionally, it is possible to have an endemic equilibrium even when $\mathcal{R}_0 < 1$.

4.4.   Local stability analysis

In the previous section, we have shown that the COVID-19-free equilibrium always exists, while the endemic equilibrium exists depending on the conditions of $\mathcal{R}_c$ and $\mathcal{K}$. In this section, we show how the stability of the COVID-19-free equilibrium can only be achieved if $\mathcal{R}_c < 1$.

Theorem 5. For $e_1 \neq 0$ and $e_2 \neq 0$, the point $\mathcal{E}_0$ in the system (1) is locally asymptotically stable if $\mathcal{R}_c < 1$ and unstable if $\mathcal{R}_c > 1$.

Proof. To demonstrate the stability of the COVID-19-free equilibrium point, we linearize system (4.1) around the COVID-19-free equilibrium point ^[73]. The Jacobian matrix is formulated as follows:

The eigenvalues of the matrix above can be obtained as follows:

$\mathcal{E}_0$ is locally asymptotically stable if all eigenvalues of matrix (4.6) have negative values. Since $\lambda_1$ to $\lambda_3$ are already negative and $\lambda_4 < 0$ for $\mathcal{R}_c < 1$, $\mathcal{E}_0$ is locally asymptotically stable for $\mathcal{R}_c < 1$. Moreover, $\mathcal{E}_0$ is unstable if any eigenvalue of matrix (4.6) is positive. Therefore, $\mathcal{E}_0$ is unstable for $\mathcal{R}_c > 1$.

4.5.   Bifurcation analysis

Here, we analyze the local stability of the COVID-19 endemic equilibrium point near the condition $\mathcal{R}_c = 1$ using the Castillo-Song bifurcation theorem ^[74].

Theorem 6. The COVID-19 model in system (4.1) undergoes backward bifurcation at $\mathcal{R}_c = 1$ if $\mathcal{K} < 0$, and forward bifurcation when $\mathcal{K} > 0$.

Proof. Using the Castillo-Song theorem ^[74], we examine bifurcation from system (2.1), which has already been nondimensionalized to become system (4.1). System (4.1) can be represented as follows:

For $\mathcal{R}_c = 1$, we take the bifurcation parameter $a_1 = a_1^{*}$ where $a_1$ represents the contact rate $\beta_1$. We have:

Furthermore, the linearization matrix around the COVID-19-free equilibrium point for $a_1 = a_1^{*}$ is obtained as follows:

Subsequently, we computed the right and left eigenvectors of the zero eigenvalues of $\mathcal{A}$, denoted as $\textbf{w} = (w_1, w_2, w_3, w_4)^T$ and $\textbf{v} = (v_1, v_2, v_3, v_4)$, respectively, yielding the following expression:

Using the theorem ^[74] and computing the nonzero partial derivatives of $f(x)$ at the disease-free equilibrium, the associated bifurcation coefficients are defined as follows:

Since $\mathcal{B} > 0$, according to ^[74], backward bifurcation occurs when $\mathcal{A} > 0 \iff \mathcal{K} < 0$, and forward bifurcation occurs when $\mathcal{A} < 0 \iff \mathcal{K} > 0$. Hence, the theorem is proved.

5.   Numerical experiments

In this section, we present the results of numerical experiments conducted on model (2.1) regarding the elasticity analysis, bifurcation diagram, and autonomous simulations. The control reproduction number of the original Model in (2.1) is given by:

Hence, without any intervention, the basic reproduction number of the original model in (2.1) is given by:

All model parameters that have been used in this section are given in Table 1, which shows that $\mathcal{R}_c = 1.827 > 1$. Without any improvement in the control variables, COVID-19 will continue to exist in West Java.

5.1.   Elasticity analysis

At this stage, an analysis will be carried out to examine the sensitivity of the control reproduction number ($\mathcal{R}_c$) to various parameters that affect its value, as well as the relationship between these parameters. By considering different values for these parameters, we can observe how the control reproduction number of the system (2.1) varies while keeping the values of the other parameters constant. To evaluate the sensitivity, we employ the forward sensitivity index of $\mathscr{R}_c$, which is denoted by:

where $\sigma$ is a parameter in the system (2.1).

We calculate the elasticity indices of $\mathscr{R}_c$ with respect to all parameter values to observe the effects of changes in parameters on $\mathscr{R}_c$. Using the formula and parameter values in Table 1, we obtained the following results:

To analyze the impact of changes in vaccination intensity and quarantine measures, we recalculate the elasticity index by adjusting the values of $u_1$ and $u_2$. The parameters used are the same as those shown in Figure 3. It is evident that reducing the intensity of vaccination and quarantine increases the elasticity of vaccine intervention while decreasing the elasticity of quarantine. Conversely, when we increase the intensity of vaccination and quarantine, the effectiveness of vaccination decreases while the elasticity of quarantine increases. Additionally, we observe that the quarantine intervention is more elastic in affecting $\mathcal{R}_0$ compared to vaccination. Another noteworthy point from Table 3 is that the infection rate of vaccinated individuals is much more sensitive compared to that of non-vaccinated individuals. Hence, it is essential to develop a high-quality vaccine that can provide maximum protection against the COVID-19 infection.

5.2.   Bifurcation diagram and autonomous simulation

In this section, we draw the bifurcation diagram using the best-fit parameter values obtained from the previous section. The result is presented in Figure 4, where the forward bifurcation occurs at $\mathscr{R}_c = 1$. The blue and red curves represent the equilibrium states of COVID-19-free and endemic conditions, respectively. The solid and dotted curves represent the stable and unstable equilibrium points. $BP$ denotes the branching point at $\mathscr{R}_c = 1$. It is clear to see that when $\mathscr{R}_c < 1$, the COVID-19-free equilibrium is stable, but it becomes unstable when it surpasses the branching point $BP$. At the same time, at the branching point $BP$, the endemic equilibrium starts to emerge as $\mathscr{R}_c$ increases.

To illustrate the dynamics of the Model (2.1), we use four sample points in Figure 4, denoted as $P_1, P_2, P_3$ and $P_4$. We observe that $P_1$ lies in the region where $\mathscr{R}_c < 1$, and the autonomous simulation is shown in Figure 5, panel (a), for the original model in the Model (2.1). It is evident that the system tends toward a stable COVID-19-free equilibrium. On the other hand, for sample points $P_2, P_3$ and $P_4$, when $\mathscr{R}_c > 1$, the system tends toward an endemic equilibrium. The numerical simulation results are given in panels (b)–(d), respectively.

6.   Optimal control model

6.1.   Optimal control characterization

Model (2.1) considers two control variables: $u_1$, representing vaccination intervention, and $u_2$, representing quarantine intervention. Trivially, from our sensitivity analysis in the previous section, we understand that larger values of $u_1$ and $u_2$ will yield better results for eradicating COVID-19. However, it will come at a high cost. Hence, our main objective is to minimize the number of infected individuals, $I$ and $Q$, with the lowest possible cost of $u_1$ and $u_2$. To address this problem, we constructed a cost function as follows:

which represents the total cost due to the high number of infected individuals and the level of control intensity. Note that $\omega_1$ and $\omega_2$ represent the weight parameters for $I$ and $Q$, respectively. On the other hand, $\varphi_1$ and $\varphi_2$ represent the weight costs for $u_1$ and $u_2$. Next, by applying Pontryagin's maximum principle ^[75], we obtain the Hamiltonian function as follows:

where $\lambda_1$ for $i = 1, 2, 3, 4, 5$ represents the adjoint variables for $S, V, I, Q,$ and $R$, respectively. Thus, taking the partial derivatives of $H$ with respect to each of the state variables yields the adjoint system, as shown below:

completed with the transversality condition $\lambda_i(t = T) = 0$ for $i = 1, 2, 3, 4, 5$. Furthermore, the optimal controls $(u_1^*, u_2^*)$ are given by:

where $u_i^{\min}$ and $u_i^{\max}$, for $i = 1, 2$, represent the lower and upper bounds of $u_1$ and $u_2$, respectively.

6.2.   Optimal control simulations

The solution of our optimal control problem in Eq (6.4) through analytical methods is unattainable due to the high dimensionality and nonlinearity of the ordinary differential equations system for the state variables in Model (2.1) and the adjoint variable in (6.3). Consequently, we employ a forward-backward sweep method to solve our problem numerically ^[73]. Additional examples demonstrating the implementation of this method can be found in ^{[76,77,78,79]}. The initial step involves providing an initial guess for the control variables $u_1$ and $u_2$ for all $t \in [0, T]$. By combining this guess with an initial condition for the state variables, we solve Model (2.1) in the time interval $t\in [0, T]$ forward in time. Subsequently, utilizing the obtained results, we solve the adjoint system in (6.3) backward in time, subject to the transversality condition $\lambda_i(T) = 0$. With the solutions for the state and adjoint variables in hand, we update the control variables using the formula in (6.4). This iterative process is repeated until the convergence criteria is met, with convergence defined as $|J^{\text{iteration-k+1}} - J^{\text{iteration-k}}| \leq 10^{-5}.$

We have divided the numerical experiments in this section into three different categories for three different purposes. The first simulation is conducted to determine the most cost-effective scenario based on the combination of controls. The second simulation aims to assess the impact of the initial population condition on the dynamics of control. The last simulation is conducted to evaluate the influence of the initial basic reproduction number on the effectiveness of the COVID-19 control strategy.

6.2.1.   Simulation for different strategies

In this section, we conducted a numerical experiment to address our optimal control problem, exploring various combinations of control interventions. The first scenario involves the simultaneous implementation of vaccination and quarantine measures to contain the spread of COVID-19. The outcomes are presented in Figure 6. It is evident that when both of these interventions are employed concurrently, there is a notable reduction in the number of infected individuals. The dynamics of vaccination control necessitate an initial high rate of application during the early stages of the simulation, followed by a gradual decrease over time. Conversely, the quarantine intervention exhibits a decreasing trend in response to a heightened number of infected individuals being placed under quarantine. The total cost of intervention for this scenario is $3.51 \times 10^6$.

The subsequent simulation aims to elucidate the population's response when vaccination is implemented as the sole intervention. The results are presented in Figure 7. It is evident that without quarantine measures, vaccination intervention necessitates a more intensive application, as indicated in panel (f). Consequently, the vaccinated population approaches nearly the entire total population, while the number of infected individuals, although exhibiting an increase, remains significantly lower in comparison to scenarios without any control measures. The overall cost of intervention in this scenario is only marginally higher than that of the first scenario, where vaccination and quarantine are implemented concurrently.

The final simulation was conducted to assess the impact of quarantine intervention when employed as the sole measure to mitigate COVID-19. The results are depicted in Figure 8. It is evident that without vaccination intervention, quarantine measures do not have a significant effect on curtailing the number of infected individuals, unlike the two preceding scenarios. To counterbalance the high number of infections, intensive vaccination interventions must be administered in the initial stages (refer to panel (f)). However, since our model assumes that the recovery of quarantined individuals follows a saturation function, it becomes increasingly challenging to increase the recovery rate as more individuals are placed under quarantine. This is because there are a high number of individuals who require treatment. Therefore, this intervention must be scaled back over time to minimize intervention costs. Consequently, quarantine interventions do not yield significant results in reducing infection numbers, as observed in the two previous scenarios.

6.2.2.   Simulation for different initial conditions

In contrast to the simulations in the previous subsection, which focused on intervention combinations, this subsection focuses on scenarios driven by distinct initial conditions. These simulations are conducted under the assumption of very low initial infection counts, and this scenario can be termed an endemic prevention scenario, where intervention is implemented before rapid COVID-19 transmission occurs. The initial condition is given by:

The simulation outcomes are illustrated in Figure 9. It is evident that when the number of infected individuals at the initial time is extremely low, the intervention administered is also relatively modest (refer to panel (f) in the blue curve). Vaccination intervention begins to escalate as the number of infected individuals increases (see panels (c)–(e) highlighted in red). On the other hand, the intensity of quarantine intervention diminishes over time, reflecting the reduced effectiveness of treatment intervention for individuals who are quarantined. This is due to the increasing number of individuals who are subject to quarantine. The total cost incurred for this scenario is $5.33 \times 10^6$, which is greater than when the initial condition involved a higher number of infected individuals. However, the number of new infections successfully averted is substantially larger, amounting to $6.51 \times 10^7$.

6.2.3.   Simulation for different initial basic reproduction numbers

The final simulation in this chapter was conducted to assess the influence of the endemic level, as indicated by the magnitude of $\mathscr{R}_c$ in the population, on the intensity of control measures required to manage COVID-19. To achieve this, the parameters used were consistent with those in Table 1, except for $\beta_1$ and $\beta_2$, which were reduced by 50%. As a result, we obtained $\mathcal{R}_c = 4.658$, which is 50% smaller than the estimate from the previous chapter. The results are displayed in Figure 10. Due to the lower value of $\mathscr{R}_c$, the intensity of vaccination control shown in Figure 10 (panel (f)) does not need to be as high as indicated in Figure 6 (panel (f)). Nonetheless, interventions of this nature have proven effective in significantly reducing the number of infected individuals (see panel (e)). The total cost incurred in this scenario is $9.43\times 10^6$, which is substantially smaller than the scenario in Figure 6. The number of infections successfully averted amounts to $3.4 \times 10^6$.

7.   Discussion and conclusions

In this paper, we delved into an analysis of the dynamic characteristics of the COVID-19 transmission model, specifically the susceptible-vaccinated-infected-quarantine-recovered $(SVIQR)$ model. Our investigation considers critical factors such as human awareness and resource constraints for treating quarantined individuals. To estimate model parameters, we utilized incidence data sourced from West Java, Indonesia. Our exploration encompassed mathematical analyses of equilibrium points, their local stability, and the control reproduction number ($\mathscr{R}_c$), all of which are presented analytically.

Our findings revealed a significant insight: The model may display forward or backward bifurcations at $\mathscr{R}_c = 1$. This suggests that achieving a reproductive number below one alone is not always sufficient to ensure the complete eradication of COVID-19 from a population. This result is consistent with another COVID-19 model introduced in ^[33,76]. In practice, backward bifurcation commonly occurs due to several possible factors, such as reinfection ^[80], relapse, imperfect vaccination ^[81,82], or limited treatment capacity ^[83]. In the context of our study on COVID-19, the existence of the backward bifurcation suggests the potential for decision-making challenges, as the expectation that COVID-19 will always disappear when $\mathscr{R}_c < 1$ is not consistently observed in real-world scenarios. The implications of this backward bifurcation highlight the possibility of a significant shift from a stable COVID-19-free state to a large endemic state, emphasizing the complex dynamics of disease control in public health. Another important result of our research is the sensitivity analysis. Sensitivity analysis underscores the potential effectiveness of minimizing contact between infected individuals who are not in quarantine and vaccinated individuals as a means of controlling the disease's spread.

To address budget limitations, we reconstructed our model as an optimal control model, introducing vaccination and quarantine interventions as time-dependent parameters. By employing Pontriagin's maximum principle, we rigorously formulated our optimal control problem and utilized the forward-backward sweep method for numerical solutions. Our simulations, conducted across various scenarios that may manifest in real-world situations, yielded noteworthy insights.

To begin, our numerical experiments underscore the paramount importance of vaccination as the primary intervention in curtailing infections. Vaccination has been proven effective in reducing the likelihood of vaccinated individuals contracting the virus, thereby minimizing the number of infected individuals. However, combining vaccination with quarantine measures still yields more promising results than relying solely on vaccination. A second critical finding emphasizes the necessity of implementing interventions at the onset of an outbreak rather than waiting for it to escalate. This early intervention strategy leads to a significantly higher number of preventable infections. {The results of our optimal control simulation study underscore the importance of considering vaccination and quarantine strategies in controlling the spread of COVID-19. Both interventions play distinct roles: Vaccination is highly effective in preventing the initial spread of COVID-19 during the early stages of the pandemic, while quarantine measures are best focused as the pandemic begins to escalate. The combination of both strategies is essential for effectively mitigating the outbreak of COVID-19. Implementing these strategies based on needs will yield optimal results, resulting in a reduced budget for intervention and preventing a larger number of COVID-19 infections.

While our study sheds considerable light on the significance of the COVID-19 vaccination, there are several possible future directions in COVID-19 research from this article and other references from other disciplines. In many instances, there exists a gap between genetics, environmental science, and epidemiology in the context of COVID-19. Potential novel insights lie in understanding the determinants of COVID-19 susceptibility and informing targeted interventions and public health strategies to mitigate the impact of the pandemic. Readers can explore a mathematical approach for this scenario in ^[84]. Another promising avenue in COVID-19 modeling involves bridging rumor spreading models with COVID-19 dynamics, leveraging insights from the Maki-Thompson model ^[85]. Readers can adapt the Maki-Thompson model to incorporate elements relevant to COVID-19 transmission dynamics, such as contact rates, transmission probabilities, and population demographics. Furthermore, readers can delve into research directions by exploring the concept of strategic decision-making and resource allocation in the context of pandemic response, utilizing a stochastic control approach ^[86]. Through this approach, we can gain a better understanding of the effectiveness of strategies for managing public health crises and ultimately contribute to improved preparedness and resilience in the face of future pandemics. In conclusion, by embracing interdisciplinary approaches and leveraging mathematical modeling techniques, we can pave the way for innovative solutions to combat COVID-19 and enhance our preparedness for future pandemics.

Use of AI tools declaration

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

Acknowledgments

This research is funded by Universitas Indonesia through the PUTI Research Grant Scheme, 2023 (ID Number: NKB-469/UN2.RST/HKP.05.00/2023).

Conflict of interest

The authors declare there is no conflict of interest.