Optimal control on COVID-19 eradication program in Indonesia under the effect of community awareness

1.   Introduction

An outbreak of the novel coronavirus was first identified in Wuhan, China, and was quickly transmitted worldwide. As of September 7, 2020, approximately $27, 475, 333$ individuals have been infected ^[1]. In Indonesia, the first confirmed case was reported on March 02, 2020, and as of September 7, 2020, approximately $196, 989$ cases have been detected ^[2]. DKI Jakarta, West Java, and East Java are the three provinces with the highest number of COVID-19 cases. With no vaccine becoming available in the near future, it is predicted that the number of COVID-19 infections may still increase. Therefore, the implementation of current strategies, such as social and physical distancing, medical masks, and other nonpharmaceutical interventions, should be investigated further.

Among these strategies that have been implemented to minimize the risk of infection are social distancing and medical mask use. The implementation of these strategies can be more effective if there is a high level of individual awareness. If people become aware of the importance of such interventions, they would decide to implement these strategies. To deeply understand the transmission dynamics of the disease and intervention effectiveness, extensive research using mathematical models could be used as an alternative solution.

Mathematical models have a long history of application to help humans understand how the dynamics of a disease spread in a population, for example in dengue ^[3,4,5], malaria ^[6,7,8], tuberculosis ^[9,10], and many more. These models try to accommodate various essential factors in the spread of a disease, such as the presence of a disease vector, the phenomenon of relapse and reinfection, symptomatic and asymptomatic cases, analysis of the success of interventions with limited costs, and others. The existence of the influence of community awareness has also been discussed by several previous authors, including in ^{[10,11,12,13]}. This awareness effect has been applied to the model as a variable that has its dynamics, as well as a parameter that affects the behavior of other parameters such as the rate of infection, treatment, and others.

Since COVID-19 rapidly spread in many countries in the world, numerous mathematical models have been conducted ^{[14,15,16,17,18]} to give a better understanding on the transmission mechanism of COVID-19, and how to prevent it. Similar to the mentioned mathematical models in the previous paragraph, this mathematical model tries to divide the human population based on their health status, then a dynamic analysis is carried out in-depth. Kurchaski et al. ^[15] formulated a stochastic mathematical model to assess the variations in disease transmission and measure the probability of newly introduced cases triggering outbreaks in other regions based on the transmission variation in Wuhan, China. Prem et al. ^[16] used mathematical models to assess the effects of changes in the population interacting with each other during the outbreak. They found that restricting activities would help to minimize the number of infections. Furthermore, research showed that the use of masks could reduce the disease transmission dynamics ^[19]. Based on the model which considers the detected/undetected cases and the symptomatic/asymptomatic cases, Giordano et al. ^[20] found that combination of restrictive social-distancing, rapid testing and contact trace should be implemented partially to reduce the spread of COVID-19. A similar result also given by Aldila et al. ^[21].

Although many studies have been conducted to understand COVID-19 transmission dynamics and intervention, only a few studies on the COVID-19 disease transmission in Indonesia have been performed ^[17,18,21]. Using the SEIR model, Soewono ^[17] calculated the basic reproduction number of COVID-19 in DKI Jakarta, which was found to be 2.5148. In addition, Ndii et al. ^[18] formulated deterministic and stochastic mathematical models and determined the probability of disease extinction as well as calculating the basic reproduction number for Indonesia. They found that reducing the contact rate by approximately 70% can minimize the number of infected individuals. Aldila et al. ^[21] investigated the effects of rapid testing and social distancing in controlling the spread of COVID-19 in Jakarta–Indonesia and found that a massive rapid-test intervention should be implemented if strict social distancing is relaxed. However, none of these works analyzed the effects of the individual's awareness or included a cost-effectiveness analysis. According to ^[23,24], developing human awareness to COVID-19 could help to prevent the spread of COVID-19. This campaign help to disseminate and help to dispel misinformation on COVID-19, and in the same way to promote precautionary measures like washing hands, physical distancing, and many other prevention strategies. Hence, it is important to discuss how the awareness of the community effect the spread of COVID-19.

In the present study, we formulated a mathematical model with nonpharmaceutical and pharmaceutical interventions for COVID-19 control programs and performed a cost-effectiveness analysis. The interventions considered were the use of medical masks, rapid test interventions, and pharmaceutical interventions, which are all affected by community awareness. As the purpose of the study is to determine the best intervention strategies, the optimal control approach was implemented in the model. The effectiveness of all the intervention schemes was analyzed using the infection averted ratio (IAR) and average cost-effectiveness ratio (ACER) methods. To test our model, we estimated the values of the parameters using incidence data from Jakarta, West Java, and East Java in Indonesia.

Compared with other established COVID-19 transmission model ^{[14,15,16,17,18,19,20,21]}, the novelty of this paper lies in several important features. Compared to the model in ^[20], our model accommodates the latent stage before an individual is ready to spread COVID-19. Compared to ^[15,16,17], our model accommodates detected/undetected cases. The first novelty is that our model accommodates how human awareness of the detected cases may reduce the infection rate. Furthermore, our model also considers the use of medical masks massively, not only for infected individuals but also for susceptible and recovered individuals. Our parameters were estimated using incidence data from West Java and East Java, which has never been used previously, and also incidence data from DKI Jakarta.

The remainder of the paper is organized as follows. Section 2 presents the basic model and the model's analysis regarding the equilibrium, basic reproduction number, and bifurcation analysis. Section 3 extends the model from Section 2 to include the effects of community awareness on the infection rate. The characterization of the optimal control problem is also conducted in this section using Pontryagin's maximum principle. The parameter estimations are presented in Section 4, followed by numerical experiments for the autonomous system. The numerical experiments described in Section 5 consist of a sensitivity analysis, optimal control simulations, and a cost-effectiveness analysis. Finally, conclusions for our research are presented in the last section.

2.   The basic model

2.1.   Model formulation

We divide the human population into seven categories based on their health status: susceptible individuals $(S)$, exposed/latent individuals $(E)$, asymptomatic individuals $(A)$, symptomatic individuals $(I)$, quarantined individuals $(Q)$, hospitalized individuals $(H)$, and recovered individuals $(R)$. Hence, the total human population at time $t$ is given by

A flowchart of the model is given in Figure 1. The susceptible population is generated by recruitment through births at a constant rate of $\Lambda$. Due to the direct contact with an infected individual $A$ and $I$, susceptible individuals will become infected by COVID-19 and transferred to the exposed compartment. For the infection process, the following assumptions have been used. a) Exposed individuals cannot spread COVID-19 as they are still in the incubation period. b) Infected individuals who undergo quarantine or hospitalization cannot spread the disease. c) As asymptomatic individuals do not show any symptoms, which in this case means they do not spread pathogens from sneezing as often as symptomatic individuals. Li et al. ^[22] showed that the infection rate of asymptomatic cases is lower than that of symptomatic cases and hence we set that the infection rate of asymptomatic individuals is lower than that of symptomatic individuals, with the reduction parameter $\sigma \in (0, 1)$. Therefore, the infection term in our model is given by

We assume that an intervention to reduce the transmission rate has been implemented in the population that forces all individuals to use a medical mask. Assume that a $q$ proportion of $N$ use a medical mask, while the rest $1-q$ do not. Hence, the total mass contact between each $S$ and $I$ is given by

If we assume that an individual who uses a medical mask reduces the chance to spread/receive pathogens from other individuals with the reduction factor $\phi$, the total of newly infected individuals caused by contact between susceptible and symptomatic individuals is given by

Similarly, the total of newly infected individuals caused by contact between susceptible and asymptomatic individuals is given by

After the incubation period $\alpha^{-1}$, a portion of exposed individuals $p$ progress to be symptomatic infected individuals, while the remainder, $1-p$, are determined to be asymptomatic. Note that there is a possibility that exposed individuals could recover without progressing to be infected individuals due to their immune systems. We symbolize this possibility as $\gamma_1$. The number of infected individuals $I$ then decreases due to the recovery rate $\gamma_3$ or rapid assessment at a rate of $\eta_1$. Detected symptomatic individuals will be forced to be hospitalized. Hence, the $H$ compartment increases due to the transition from $I$ at a rate of $\eta_1$. Similarly, the number of asymptomatic individuals decreases due to a natural recovery rate of $\gamma_2$ and rapid assessment $\eta_2$, which transfers them into a $Q$ compartment. Asymptomatic individuals may progress into symptomatic individuals. Hence, we have a transition from $A$ to $I$ at a rate of $\xi_1$, and from $Q$ to $H$ with the rate of $\xi_2$. Finally, the number of infected individuals in the quarantine and hospitalized compartments will decrease due to the recovery rates of $\gamma_4$ and $\gamma_5$, respectively. According to ^[25], from more than 23 million cases in the world on August 25, 2020, the reported death cases has been only 816, 469 cases (Case fatality rate is around $\pm 3\%$). Hence, since our model's simulation only conducted for a short period, it is sufficient to assume that the death rate induced by COVID-19 to be neglected.

Hence, the full model equation in terms of the rate of change for each sub-population is given as follows.

where $\mu$ is the natural death rate. Notably, the model (2.3) is supplemented with non-negative initial conditions

Note that all parameters are positive, with their descriptions given in Table 1.

2.2.   Model analysis

2.2.1.   Basic properties

For the biological significance of model (2.3), it is necessary to ensure that all variables in our model are non-negative at time $t\geq0$ when the initial condition is also non-negative. Hence, we have the following lemma.

Lemma 1. For a non-negative initial condition for system (2.3) as follows

the solutions of $S(t), E(t), A(t), Q(t), I(t), H(t),$ and $R(t)$ of model (2.3) are non-negative for all time $t > 0$.

Proof. Please see Appendix A for the proof.

Lemma (1) guarantee each variable's positivity for all time $t > 0$, which is needed for epidemiological interpretation. In our cases, all our variables in system (2.3) represent the number of the human population. Furthermore, the following lemma completes the well-posed biological properties of our model, which guarantee that each variable in system (2.3) is bounded for $t \to \infty$.

Lemma 2. The feasible region $\Omega$ defined by

is positively invariant under system (2.3).

Proof. Please see B for the proof.

2.2.2.   Equilibrium and the basic reproduction number

With a straightforward calculation, system (2.3) has a unique disease-free equilibrium as follows.

Using the next-generation matrix ^[40], the basic reproduction number $(\mathscr{R}_0)$ of system (2.3) is given by

It is possible to rewrite the expression of $\mathscr{R}_0$ to account for the source of infection as follows

where

This expression shows infections resulting from two sources, namely, asymptomatic and symptomatic individuals. The first component $\left(\dfrac{\Lambda}{\mu}\right)$ in each term represents the total number of individuals who can be infected by infected individuals at an early stage of the infection being spread.

The second expression in each component is $\left(\dfrac{ (1-q+q\phi)^2 \beta}{\mu+\alpha+\gamma_1} \right)$, which represents the infection that is affected by the use of a medical mask during the lifetime period of an individual in the exposed compartment. Hence, increasing the value of the denominator will reduce this component. However, we cannot increase the value of $\alpha$ as it leads to the acceleration of the incubation period, and it is not possible to shorten the life expectancy of humans. Therefore, the only possible way is by increasing the recovery rate of $\gamma_1$. Before the virus becomes active in the human body, there is a short time period during which the immune system will fight against the virus. If the immune system successfully kills the virus, the human will recover. Therefore, increasing $\gamma_1$ is highly related to increasing the effectiveness of the immune system in the human body, such as increasing its endurance through additional supplements/vitamins, exercise, and the consumption of healthy food.

The difference between $\mathscr{R}_{\text{Asymptomatic}}, \mathscr{R}_{\text{Symptomatic-1}},$ and $\mathscr{R}_{\text{Symptomatic-2}}$ is presented in their third component. In $\mathscr{R}_{\text{Aymptomatic}}$, the infection term is multiplied by the ratio between the transition from being exposed to the infection period of the asymptomatic class. To become a symptomatic individual from the exposed compartment, two paths can be taken. The first path is $E \to I$ directly, which contributes to $\mathscr{R}_{\text{Symptomatic-1}}$, while the second path is $E \to A \to I$, which contributes to $\mathscr{R}_{\text{Symptomatic-2}}$. Notably, the third component for each $\mathscr{R}$ is multiplied by $\alpha$, which represents the incubation period of coronavirus. Hence, we understand that a shorter incubation period of COVID-19 will increase $\mathscr{R}_0$.

An important result is stated in the following theorem.

Theorem 1. The disease-free equilibrium $\mathcal{E}_0$ of system (2.3) is locally asymptotically stable if $\mathscr{R}_0 < 1$, and unstable if $\mathscr{R}_0 > 1$.

This theorem has been reviewed by the author in ^[41]. Hence, we do not show it in this article. The theorem implies that it is possible to eradicate COVID-19 if this threshold is less than unity. The basic reproduction number is defined as an expected number of secondary cases due to infection from one primary case during its infection period in a completely susceptible population ^[42]. This means that the basic reproduction number picturing the number of new COVID-19 which is produced by one infected individual when the initial condition of the population is at the COVID-19 free state. Hence, it is understandable that the number of an infected individual will increase and tends to the endemic state whenever the basic reproduction number is larger than unity. Many epidemiological models generate the same results (see ^{[43,44,45,46]} for some examples). However, not always $\mathscr{R}_0 < 1$ indicates the disease may not persist. When backward bifurcation appears, another stable equilibrium, which in this case is the endemic equilibrium, is locally stable. Please refer to ^{[21,47,48,49]} for examples. Hence, it is important to understand the bifurcation type of our proposed model in (2.3).

The next equilibrium is the endemic equilibrium point, $\mathcal{E}_1$, which is given by

where

while $I_1$ is taken from the positive roots of

with

As $c_1 < 0$ and $c_0 > 0 \iff \mathscr{R}_0 > 1$, we have the following theorem regarding the existence of the endemic equilibrium.

Theorem 2. There exists a unique endemic equilibrium $\mathcal{E}_1$ of system (2.3) whenever $\mathscr{R}_0 > 1$.

Based on Theorems 1 and 2, our model exhibits a change in stability and the existence of an equilibrium at $\mathscr{R}_0 = 1$. Hence, we investigate the stability of the endemic equilibrium using the well-known Castillo-Song theorem ^[50] at $\mathscr{R}_0 = 1$.

Theorem 3. System (2.3) undergoes a forward bifurcation at $\mathscr{R}_0 = 1$.

Proof. Please refer to Appendix C for the proof of this theorem.

Theorems (2) and (3) indicate that COVID-19 will exist if the threshold number $\mathscr{R}_0$ larger than unity. Hence, whenever the COVID-19 free equilibrium $\mathcal{E}_0$ exist and stable, then the COVID-19 equilibrium does not exist, and vice versa. Our analysis in this section suggests the importance of paying attention to the size of $\mathscr{R}_0$. From the expression of $\mathscr{R}_0$ in (2.6), it can be seen that $q$ and $\eta_1$ are inversely proportional to $\mathscr{R}_0$. Hence, increasing both of these parameters will reduce $\mathscr{R}_0$, as shown in Figure 2. It is easy to verify that the minimum proportion/rate of medical mask/rapid testing required to reach $\mathscr{R}_0 < 1$ will decrease whenever $\beta$ decreases. Hence, it is important to reduce the transmission rate simultaneously with other interventions. Naturally, the transmission rate will be decreased when the community becomes aware of the existence of the disease. The massive amount of information about COVID-19 provided through social media, TV, or other sources could lead to community awareness of the disease. Hence, we improve our model to consider the community awareness effect on the transmission rate in the following section.

3.   Improvement of the model for a COVID-19 control program under the awareness effect

3.1.   Model formulation

Here, we improve our proposed model in (2.3) by considering two factors: population awareness, which will decrease the infection rate, and the introduction of time-dependent control variables.

The first improvement is to include the population awareness of the danger of COVID-19. We assume that the more aware the community is, the more readily the infection rate will decrease. Let $m$ describe the level of awareness of the community. A larger $m^{-1}$ indicates a high level of community awareness, while a low $m^{-1}$ indicates a low awareness level. Next, we assume that the community awareness depends on the number of reported cases, which in our model is the individuals in $H$ and $Q$. Hence, instead of treating $\beta$ as a constant, it should be a function that depends on $m, H$, and $Q$, namely, $\bar{\beta}(m, H, Q)$. This function should be a monotonic decreasing function with respect to $H$ and $Q$, i.e., $\dfrac{\partial \bar{\beta}}{\partial H} < 0$, $\dfrac{\partial \bar{\beta}}{\partial Q} < 0$ and $\dfrac{\partial \bar{\beta}}{\partial m} < 0$. Furthermore, $\lim_{H+Q\to \infty } \bar{\beta} = \bar{\beta}_{\min}$ and $\lim_{m^{-1}\to \infty } \bar{\beta} = \bar{\beta}_{\min}$, which represents the transmission rate, tends to its minimum values whenever the number of reported cases tends to infinity or when the awareness level is high (large $m^{-1}$). Conversely, $\bar{\beta}$ should also fulfill $\lim_{H+Q \to 0} \bar{\beta} = \bar{\beta}^{\max}$ and $\lim_{m^{-1} \to 0} \bar{\beta} = \bar{\beta}^{\max}$ which represent the transmission rate tends to it maximum value whenever the number of reported cases and awareness are very low. In this article, we chose the following infection function:

where $\beta_0$ is the maximum of $\bar{\beta}$ if the level of awareness is low ($m^{-1}\to 0$) or $H+Q\to0$. With this chosen function, we have:

Next, we consider the intervention parameters to be time-dependent variables. First, we introduce the control variable as the proportion of individuals who use a medical mask, denoted by $u_1(t)$. The second and third controls are for rapid test interventions for asymptomatic and symptomatic individuals, respectively. Hence, we change $\eta_1$ and $\eta_2$ into $u_2(t)$ and $u_3(t)$, respectively. The last control variable is the pharmaceutical interventions dedicated to accelerating the recovery rate for hospitalized individuals, denoted by $u_4(t)$.

Hence, the model for COVID-19 transmission, considering the level of awareness and interventions (pharmaceutical and nonpharmaceutical), is given as follows.

supplemented with nonnegative initial conditions

3.2.   Model analysis

For the model analysis, let us assume that all control variables are constant parameters; hence, $u_i(t) = u_i,$ for $i = 1, 2, 3, 4$. Therefore, the COVID-19 model in (3.2) can now be written as follows.

Using a similar approach as with Lemma (1) and Lemma (2), we also have the following two properties for model in (3.3)

Lemma 3. For the non-negative initial conditions for system (3.3) as follows

the solutions of $S(t), E(t), A(t), Q(t), I(t), H(t),$ and $R(t)$ of model (3.3) are non-negative for all time $t > 0$.

Lemma 4. The feasible region $\Omega$ defined by

is positively invariant under system (3.3).

Similar to the previous model in (2.3), the awareness-based model in (3.3) has a COVID-19-free equilibrium given by:

Using the next-generation matrix ^[40], the basic reproduction number $(\mathscr{R}_0)$ of system (3.3) is given by

The local stability of $\mathscr{E}_0$ is also determined by $\mathscr{R}_0^*$ for the case of the awareness-based model in (3.3). This is stated in the following theorem.

Theorem 4. The disease-free equilibrium $\mathscr{E}_0$ of system (3.3) is locally asymptotically stable if $\mathscr{R}_0^* < 1$, and unstable if $\mathscr{R}_0^* > 1$.

From the form of $\mathscr{R}_0^*$, it can be seen that $m, \gamma_4, \gamma_5,$ and $u_4$ do not determine the type of stability of $\mathscr{E}_0$ as they do not appear in $\mathscr{R}_0^*$. However, these parameters play an important role in determining the size of the epidemic when the endemic equilibrium occurs and determine the speed to reach the COVID-19-free equilibrium when $\mathscr{R}_0 < 1$. We will provide the analysis using numerical simulations in Section 5 for these claims.

The endemic equilibrium of system (3.3) is given by

where

where $S_1^*$ is a function of other variables and has a considerably long expression, which could not be included in this paper, while $A_1^*$ and $I_1^*$ are taken from the intersection of the following polynomials:

Note that $k_i$ for $i = 1, 2, \dots 10$ has an extremely long expression, which could not be included in this paper. From the above expression of $\mathscr{E}_1$, which depends on the intersection between $P_1$ and $P_2$, it is difficult to determine the number of possible endemic equilibria. Thus, we plot the polynomials $P_1$ and $P_2$ in Figure 3 using the set of parameters given in Table (1).

We close our dynamical analysis in this section with the following theorem that states the bifurcation type of model in (3.3) at $\mathscr{R}_0^* = 1$.

Theorem 5. The awareness-based COVID-19 transmission model in (3.3) undergoes a forward bifurcation at $\mathscr{R}_0^* = 1$.

Proof. Please refer to Appendix A for the proof of this theorem.

Our result for the COVID-19 model which considers awareness of the population in system (3.2) shows a similar qualitative behavior with the standard model in system (2.3), where the related basic reproduction number becomes the unique threshold to guarantee the existence/persistence of COVID-19 from the population for each model. Hence, we can conclude that COVID-19 could be eradicated as long as we can reduce the basic reproduction number less than unity. Otherwise, COVID-19 will persist. Furthermore, it can be seen that $m$ does not appear in $\mathscr{R}_0^*$. This result indicates that when the community awareness only appears in the transmission term, then the community awareness does not change the condition such that $\mathscr{R}_0^* = 1$. However, it does change the size of the endemic equilibrium for each variable and the time to reach the outbreak of COVID-19. We discuss this result in more detail using numerical experiments in Section 4, Figure 5.

3.3.   Optimal control characterization

As we already stated in an earlier section of this manuscript, four different control variables will be implemented, namely, the use of medical masks to reduce the infection probability $u_1(t)$, hospitalization rate for symptomatic individuals $u_2(t)$, rapid testing to detect asymptomatic individuals and push them to conduct a self-quarantine $u_3(t)$, and increases in the medical treatment quality to accelerate recovery rate $u_4(t)$.

The optimal control problem seeks to minimize the number of people infected by COVID-19, while keeping the cost for control implementations as low as possible. To do this, let us consider the following objective function.

where $\omega_j$ for $j = 1, \dots 5$ and $\upsilon_i$ for $i = 1, \dots 4$ are positive constants that will balance the relative purposes for each term in the objective function $\mathcal{J}$, and $T_f$ is the final time of control implementation. Furthermore, let the set of admissible control sets $\mathcal{U}$ be given by:

The first component of $\mathcal{J}$, i.e., $\int_0^{T_f} \left[\omega_1 E(t) + \omega_2 A(t) + \omega_3 Q(t) + \omega_4 I(t) + \omega_5 H(t)\right] dt$, represents the cost related to the existing number of infected individuals in the field. This cost is not related to the control variables. For example, this term is related to the economic cost of the pandemic. Term $\int_0^{T_f} \left[\sum_{i = 1}^4 \upsilon_i u_i^2 \right] dt$ is related to the total cost of control implementations to achieve the eradication of COVID-19 from the community. We choose a quadratic cost function to model the cost for interventions as has already done by many authors, such as in ^{[44,45,52,53,54]}. Biologically, this quadratic function represents a condition in which a larger intervention that needs to be implemented will be more costly, which means it is not linear. For example, the increment for the implementation of medical masks from $u_1 = 0.1$ to $u_1 = 0.2$ is not difficult to implement as the number of medical masks is easy to find. In contrast, the increment from $u_1 = 0.8$ to $u_1 = 0.9$ is more difficult to implement because of the limitation of medical masks available in the field.

The sufficient condition to determine the optimal control $u_i^*$ for $i = 1, 2, 3, 4$ in $\mathcal{U}$ such that

with the constraints of our COVID-19 model in (3.3) can be obtained using Pontryagin's maximum principle ^[51]. The principle of this method is to transfer our system, which involved the state system (3.3), cost function (3.8), and minimization problem (3.10), into minimizing the Hamiltonian function $\mathcal{H}$ problem with respect to $u_i$ for $i = 1, 2, 3, 4$, that is

where $\lambda_i$ for $i = 1, 2, \dots, 7$ are the adjoint variables for the state system $S, E, A, I, Q, H, R$, respectively. These adjoint variables satisfy the following system of ODEs.

supplemented with the transversality condition $\lambda_i(T_f) = 0$ for $i = 1, 2, \dots 7$. In characterizing the optimal controls, the Hamiltonian function $\mathcal{H}$ is differentiated partially with respect to each control variable $u_i$ for $i = 1, \dots 4$, which gives us:

Hence, considering the upper and lower bound of each control parameter by $u_i^{\max}$ and $u_i^{\min}$, respectively, we can characterize the optimal controls as

for $i = 1, \dots, 4$, where $u_i^{\dagger}$ is the solution of $\frac{\partial \mathcal{H}}{\partial u_i} = 0$ for $i = 1, \dots, 4$.

4.   Parameter estimation from COVID-19 incidence in Indonesia and the autonomous simulation

In this section, the transmission rate $(\beta)$ is estimated against COVID-19 data from the DKI Jakarta, East Java, and West Java Provinces, and the other parameters are obtained from the literature and presented in Table 1. The total populations of DKI Jakarta, West Java, and East Java are 10, 374, 235, 49, 316, 712 and 39, 501, 000, respectively.

The sum of the squared error between the model and data is minimized, which is given by

where $H_t$ and $Q_t$ is the number of active cases of $H$ and $Q$ up to day $t$, respectively, while $f_t(x)$ and $g_t(x)$ is the number of active cases for $H$ and $Q$ up to day $t$ from the model's solution, respectively. The transmission rate, $\beta_0$ and $\beta_1$, is then estimated using the "lsqnonlin'' built-in function in MATLAB.

The aim of the estimation is to obtain a general insight regarding the transmission rates of the disease in these three provinces during the early incidence period. The incidence data are taken from ^[56] for Jakarta, ^[57] for East Java, and ^[58] for West Java. Each dataset was obtained during a one-month period from the beginning of the recorded incidents. The fitted values for the transmission rate of DKI Jakarta are $\beta_0 = 2.015 \times 10^{-7}$ and $\beta_1 = 0.94\times10^{-7}$. Hence, the lowest infection rate in Jakarta was $1.075 \times 10^{-7}$, which is 47% less than the initial infection rate of $\beta_0$. On the other hand, the transmission rates for East Java are $\beta_0 = 4.64 \times 10^{-8}$ and $\beta_1 = 1.856 \times 10^{-8}$. The maximum reduction in infection rate due to community awareness in East Java was 60%. Finally, the transmission rates for West Java are $\beta_0 = 4.9\times10^{-8}$ and $\beta_1 = 4.41135\times10^{-8}$. Therefore, the maximum reduction in the transmission rate was 90% in West Java. Using these parameter values, the basic reproduction number during the early spread of COVID-19 for Jakarta, East Java, and West Java are $4.18, \, 3.67,$ and $4.84$, respectively, which indicates that COVID-19 will persist in the community if no further intervention is implemented. The results of the comparison between the actual and simulated data are given in Figure 4. Although there is some systematic bias in the data and model simulations for West Java in the middle of the outbreak, the result is sufficient to meet the purpose of the estimation, which is to obtain a general insight into the transmission rate in the early period of the outbreak.

Various interventions have been conducted by policymakers in different countries, such as social distancing to reduce the contact rate, hospitalization to medicate infected individuals, contact tracing with rapid test intervention, and medical masks, the most popular intervention. In this section, we will determine how the control interventions affect the dynamics of the proposed model. The first simulation was conducted with various values for the awareness level of the community $(m)$, and other parameters were kept constant with $\beta_0$ and $\beta_1$ from the Jakarta data. The results are given in Figure 5. It can be seen that although the awareness level did not affect the size of $\mathscr{R}_0$, which in this case means it also did not determine the equilibrium stability type, it is clear that a high level of community awareness could reduce the level of the outbreak and delay the occurrence time. The second simulation was conducted to determine the effect of the medical mask intervention. To perform the simulation, we kept all parameters constant, while $u_1$ varied. The results are given in Figure 6. It can be seen that the effect of medical mask use was significant in reducing the outbreak and could delay the outbreak occurrence time. This result confirms the reason policymakers suggest the use of a medical mask, not only for the infected individual but also for the susceptible population. The application of medical mask use is based on the difficulty of finding infected people, especially those who do not show symptoms. Therefore, protecting all individuals using a medical mask is a reasonable option.

The next simulation was conducted to determine the effect of rapid testing to trace the existence of infected individuals in the community. To run this simulation, we used various values of $u_2$ and $u_3$ simultaneously, while the other parameters remained constant. From this simulation, it can be seen that rapid testing succeeded in reducing the total number of infected individuals, detected individuals, and undetected individuals. The implementation of rapid tests was successful in reducing the number of undetected cases, which in this case, transferred to the quarantined or hospitalized compartments. When these individuals were moved to these compartments, the recovery duration from COVID-19 was improved, which reduced the number of detected cases.

5.   Numerical experiment on a sensitivity analysis and optimal controls

This section presents the results of a sensitivity analysis and optimal control simulations.

5.1.   Sensitivity analysis

A global sensitivity analysis was performed using Latin hypercube sampling (LHS) in conjunction with the partial rank correlation coefficient (PRCC) multivariate analysis ^[59,60]. When the PRCC value of the parameter closes to positive or negative one, it indicates that the parameters are influential. The sign (positive or negative) indicates the relationship between the parameters and the output of interest. For example, when the sign is negative, it means that an increase in the parameter values results in a decrease in the output (in our case, it is the number of infections or the reproduction number). We measured the increasing number of infected individuals and the reproduction number (Eq (3.6)). The increasing number of infected individuals is given by the following equation

where $C_I$ is the cumulative number of infectious individuals. The aim is to investigate the influential parameters on the increasing number of infected individuals and the reproduction number. The PRCC values for each parameter measured against the increasing number of infected individuals are given in Figure 8.

For the sensitivity analysis, 2000 simulations were performed to assess the model's sensitivity to the parameters, and the results are given in Figure 8. This figure shows that the transmission rate $\beta_0$ and control rates, which are the use of masks $(u_1)$, are the most influential parameters and affect the increasing number of infectious individuals. Parameter $\beta_0$ has a positive relationship, and the control rate ($u_1$) has a negative relationship, to the model's outcome. This means that an increase in this control rate would reduce the number of infected individuals. Interestingly, the effects of the other controls are not as strong as this control. The results imply that to reduce the number of infected individuals, the use of masks should be strongly implemented. We also ran 2000 samples to determine the most influential parameters on the reproduction number and found the same results in which the transmission rate $(\beta_0)$ and control rate $(u_1)$ were the most influential parameters. The first rate has a positive relationship and the other rate has a negative relationship. The results are given in Figure 9. The control rate ($u_2$) also affects the reproduction number, although it is not as strong as the effect of medical mask use ($u_1$.)

5.2.   Optimal control simulation

The computations for the optimal control problem were performed numerically using the Runge–Kutta method of the fourth order with MATLAB. The algorithm is summarized as follows. First, an initial guess of the control variables is made and used to solve the state system (3.2) forward in time. The results for the state variables and initial guess of $u_i$ are then substituted into the adjoint system (3.12), which is solved backward in time with the transversality condition $\lambda(T_f) = 0$. Both the state and adjoint values are then used to update the control (3.14), and the process is repeated until the current state and adjoint and control values converge sufficiently ^[61]. Please see ^{[44,45,52,53,54,55]} for more examples of this method implemented in epidemiological models.

To illustrate the optimal control strategies, we used parameter values that were fitted to the COVID-19 incidence data in Jakarta and initial conditions as follows: $S(0) = 10\, 374\, 231, E(0) = 10, A(0) = 10, I(0) = 10, Q(0) = 0, H(0) = 4,$ and $R(0) = 0$. For the weight factors, we chose $\omega_i = 1$ for $i = 1, \dots, 5$ and $\upsilon_j = 10^6$ for $j = 1, \dots, 4$. It should be pointed out that these values are theoretical as they were chosen only to illustrate the control strategies proposed in this article.

In this section, we analyzed the effect of community awareness on the dynamics of the infected population and on how the control variables responded to each scenario. Based on our numerical simulations in Figure 5, smaller values of $m$ (high levels of community awareness) could reduce and delay the outbreak occurrence. To perform the simulation, we used three values of $m^{-1}$, i.e., $m^{-1} = 0.1$, which represents a high awareness level, and $m^{-1} = 0.01$ and $m^{-1} = 0.001$ to represent the medium and low awareness levels, respectively. The simulation results are shown in Figure 10. It can be seen that all simulations show a similar behavior of the control trajectories in which medical mask use provides a high intensity early in the simulation and tends to its lower bound after the outbreak has passed. As a response, the rapid testing $u_2$ and $u_3$ and hospitalization $u_4$ should start to increase to balance the decreasing of medical mask use in the community. The cost function for the cases of $m^{-1} = 0.1$, $m^{-1} = 0.01$, and $m^{-1} = 0.001$ are $2.49 \times 10^6$, $2.5\times10^6$, and $2.53 \times 10^6$. However, the numbers of infections that can be avoided in the case when $m^{-1} = 0.1$ provides the highest results, i.e., $34\, 980$, while those for $m^{-1} = 0.01$ and $m^{-1} = 0.001$ are $525$ and $74$, respectively.

To analyze the effectiveness of the optimal control simulations under the effect of community awareness, we used two cost-effective analysis techniques. The first technique is the infection averted ratio (IAR), with the formula given by

The simulation with the highest ratio is the most cost-effective. With this formula, we calculate the IAR for each scenario as follows:

Hence, as IAR $(m^{-1} = 0.1) >$ IAR $(m^{-1} = 0.01) >$ IAR $(m^{-1} = 0.001)$, we can conclude that the most cost-effective is the case in which $m^{-1} = 0.1$, which represents the highest level of community awareness.

The second analysis is the average-cost-effectiveness ratio (ACER) technique with the formula given as follows.

The lowest ACER ratio is the most effective strategy. Using the above formula, the ACER value for each scenario is given as follows.

Hence, as ACER $(m^{-1} = 0.1) <$ ACER $(m^{-1} = 0.01) <$ ACER $(m^{-1} = 0.001)$, we can conclude that the most effective strategy is when $m^{-1} = 0.1$, followed by $m^{-1} = 0.01$ and $m^{-1} = 0.001$. From both of these simulations, we can confirm the importance of community awareness to increasing the chance of a successful COVID-19 eradication program.

6.   Conclusion

Although specific medicine to cure infected individuals or a vaccine to protect the susceptible population from COVID-19 have not yet been found, various interventions have been implemented by the government in many countries, such as social/physical distancing, rapid testing, the use of medical masks, quarantines, and the improvement of hospitalization services. In this work, we presented two deterministic mathematical models in the form of systems of ordinary differential equations to describe the transmission dynamics and consider several interventions (medical masks, rapid testing, and improvement of medical treatment in hospitals), with and without community awareness. The first model was constructed by dividing the human population into susceptible, exposed, asymptomatic, symptomatic, quarantined, hospitalized, and recovered groups. The second model used a similar separation of the population but involved community awareness, which decreased the infection rate whenever the number of hospitalized and quarantined individuals increased.

Mathematical analyses showed that both models have a COVID-19-free equilibrium point, which is locally asymptotically stable if the basic reproduction number is less than unity, and unstable otherwise. The endemic equilibrium for the first model was shown analytically, and we found that it existed whenever the basic reproduction number was larger than unity. The endemic equilibrium for the second model (model with awareness), the COVID-19 endemic equilibrium, was shown numerically. The center-manifold theorem was applied to both models to analyze the bifurcation type at a basic reproduction number equal to unity. We found that both models undergo forward bifurcation when the basic reproduction number is equal to unity.

Furthermore, we used our model to fit the incidence data in the three provinces in Indonesia that have the highest recorded COVID-19 incidence, namely, Jakarta, East Java, and West Java. Our results suggest that West Java has the largest basic reproduction number, followed by Jakarta and East Java. However, we found that West Java has the highest reduction in the infection rate due to "community awareness" as the infection rate could reduce to 90% whenever the number of reported cases increased. In contrast, Jakarta has the lowest effect, with the reduction in infection rate at only 47%, compared with East Java, which has a 60% infection reduction rate.

The results of the global sensitivity analysis showed that the infection rate and control rate ($u_1$) are the most influential parameters on the increasing number of new infections. The first rate has a positive relationship, and the other rate has a negative relationship. This means that increasing the proportion of individuals who use masks results in a decrease in the number of COVID-19 infections. The same influential parameters were also found when we measured the basic reproduction number. Furthermore, the control rates, $u_2$ and $u_3$, which are rapid testing for asymptomatic and symptomatic individuals, affect the basic reproduction number and have a negative relationship. This means that increasing these control rates would reduce the basic reproduction number.

From the numerical simulations of the autonomous simulation, we found that increasing community awareness not only succeeded in suppressing the level of the COVID-19 outbreak but also delayed the occurrence time of the outbreak. Hence, we analyzed these results further using the optimal control approach. The optimal control problem was characterized using Pontryagin's maximum principle and solved numerically using the forward-backward sweep method with MATLAB. Our optimal control simulation suggests that time-dependent intervention is effective in reducing the spread of COVID-19. Furthermore, the implementation cost for the COVID-19 eradication program is more efficient when the community has a high level of awareness.

Acknowledgments

The $1^{st}$ author is financially supported by Universitas Indonesia with QQ Research Grant scheme, 2019 (Grant No. NKB-0268/UN2.R3.1/HKP.05.00/2019).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.

A.   Appendix: Proof of Lemma 1

Under the given initial conditions, from $\frac{dS}{dt}$ in system (2.3), we have

This can be re-written as

where $P = \int_0^t \left(\phi q + 1-q\right)^2 \beta (I(\tau)+\sigma A(\tau)) d \tau + \mu t.$ Therefore,

Hence,

Therefore,

In a similar way, it can be shown that $E(t)\geq 0, A(t) \geq0, Q(t)\geq0, I(t)\geq0, H(t)\geq0,$ and $R(t)\geq0$. Thus, the solution of $S(t), E(t), A(t), Q(t), I(t), H(t),$ and $R(t)$ of model (2.3) are non-negative for all time $t > 0$.

Proof of Lemma 2

From model (2.3)

Solving this, we have $N(t) = \frac{\Lambda}{\mu} + \left(N(0)-\frac{\Lambda}{\mu}\right) e^{-\mu t}$, where $N(0)$ represents the initial conditions of the total population. Thus, we have $N(t) = \frac{\Lambda}{\mu}$ as $t \to \infty$. Hence, all feasible solutions of system (2.3) enter the region

Therefore, it is a positively invariant set for system (2.3).

Proof of Theorem 3

For the system, it is assumed that

Therefore, it can be written as

Parameter $\beta$ as the bifurcation parameter is obtained by solving $\mathscr R_0 = 1$ respect to $\beta$. Next, $\mathcal{E}_0$ and the bifurcation parameters are substituted into the Jacobian matrix of system (A.1). Thus, the eigenvalue of the Jacobian matrix is obtained. Since one zero eigenvalues appear while other eigenvalues are negative, we can proceed to use the center-manifold theorem to analyze the bifurcation type of system (2.3). Furthermore, we attempt to determine the right and left eigenvectors. To find the right eigenvector, we use vector $\vec w = (w_1, w_2, w_3, w_4, w_5, w_6, w_7)$. The right eigenvector $\vec w$ is obtained as follows:

Then, we look for the left eigenvector using vector $\vec v = (v_1, v_2, v_3, v_4, v_5, v_6)$. The left eigenvector $\vec v$ is obtained as follows:

We find that the value of the eigenvector $v_1 = 0, v_4 = 0, v_6 = 0$, and $v_7 = 0$; thus, there is no need to determine partial derivatives of $g_1, g_4, g_6$, and $g_7$. Therefore, we will determine derivatives of $g_2$, $g_3$, and $g_5$ to obtain the values of $\mathcal A$ and $\mathcal B$. The non-zero $g_2$, $g_3$, and $g_5$ derivatives are as follows:

Therefore, $\mathcal A$ and $\mathcal B$ are obtained as follows:

Because $\mathcal A < 0$ and $\mathcal B > 0$, there is a forward bifurcation at $\mathscr{R}_0 = 1$ for model (2.3).

Proof of Theorem 5

For the system, it is assumed that

Therefore, it can be written as

Parameter $\beta_0$ as the bifurcation parameter is obtained by solving $\mathscr R_0^* = 1$ respect to $\beta_0$. Next, $\mathscr{E}_0$ and the bifurcation parameters are substituted into the Jacobian matrix of system (A.5). Thus, the eigenvalue of the Jacobian matrix is obtained. Since the zero eigenvalues appear while other eigenvalues are negative, we can proceed to analyze the bifurcation type of our model using the center manifold theorem. Furthermore, we determine the right and left eigenvectors. To find the right eigenvector, we use the vector $\vec w = (w_1, w_2, w_3, w_4, w_5, w_6, w_7)$. The right eigenvector $\vec w$ is obtained as follows:

Then, we determine the left eigenvector using vector $\vec v = (v_1, v_2, v_3, v_4, v_5, v_6)$. The left eigenvector $\vec v$ is obtained as follows:

We find the values of the eigenvector $v_1 = 0, v_4 = 0, v_6 = 0$, and $v_7 = 0$; thus, there is no need to determine a partial derivative of $g_1, g_4, g_6$, and $g_7$. Therefore, we will determine derivatives of $g_2$, $g_3$, and $g_5$ to get the values $\mathcal A$ and $\mathcal B$. The non-zero $g_2$, $g_3$, and $g_5$ derivatives are as follows:

Thus, $\mathcal A$ and $\mathcal B$ are obtained as follows:

Because $\mathcal A < 0$ and $\mathcal B > 0$, there is a forward bifurcation appear for model (3.3) at $\mathscr{R}_0^* = 1$.