Optimal control of dengue fever model with a logistically growing human population

1.   Introduction

Dengue is an infectious illness that is spread in almost all countries, including in Indonesia ^[1]. Severe dengue is a prominent leader of critical disease and death in some Asian and Latin American countries. Every year, dengue cases occur in East Java province, Indonesia. Based on the Ministry of Health of the Republic of Indonesia, dengue cases in Indonesia until July 2020 reached 71,633. East Java is one of the provinces in Indonesia with high endemic potential from year to year ^[2].

Dengue is a viral infection spread to humans through the bite of infected mosquitoes. Aedes aegypti mosquitoes are the major vectors that transmit the disease. The virus responsible for bringing dengue is called dengue virus (DENV). There are four DENV serotypes, and there is the potential to be infected four times. Currently, there is no particular treatment for dengue/severe dengue. Early detection of dengue disease progression and access to appropriate medical care reduce fatality rates of severe dengue to below 1%. Dengue prevention and regulation depends on effective vector control measures ^[1].

Mathematical modeling can play essential roles in helping understand the complex phenomenon of epidemic spread including the malaria ^[3,4], tuberculosis ^[5,6], measles ^[7], meningitis ^[8], monkeypox ^[9], zika ^[10], and COVID-19 ^[11,12]. With a mathematical model, the relationship between the disease spread and epidemiological parameters can be identified. Several mathematical models of dengue transmission with various types compartments have been studied by researchers ^[13,14]. In ^[15,16], the authors designed mathematical model to study the effect of temperature on the spread of dengue fever. The impact of reinfection with the same serotype on the dynamics of dengue transmission has been studied in ^[17]. Researchers at ^[18,19] developed a fractional order model to understand the dynamics of dengue fever transmission. The authors in ^[20] presented a dynamical model of dengue with hospitalization to analyze the infection in East Java province, Indonesia. Ndii ^[21] formulated a dengue model by incorporating seasonally varying mosquitoes, vaccination, and media. The authors of ^[22] discussed a dengue model that accommodates releasing Wolbachia into wild mosquitoes and vaccination. The authors in ^[23] presented the dengue model by considering the recruitment rate of the human population in terms of logistic growth. The proposed model is more realistic because the growth rate of the human population is not always constant.

This work aims to extend the previous study of ^[23] by formulating a dengue model that captures the control variables, such as fumigation and prevention, by applying the optimal control theory to the model. The model parameters will be estimated based on the dengue data in East Java Province, Indonesia. We employ efficiency and economic methods to examine the most cost-effective strategy among the implemented control scenarios. The remaining part of this paper proceeds as follows. Section 2 explains the model formulation. The analysis of the model is discussed in Section 3. The optimal control problem is studied in Section 4. Section 5 summarizes the main findings of this work.

2.   Model formulation

In this study, we consider the dengue transmission model with the assumption that the recruitment rate in the human population is logistic growth as proposed by ^[23]. The total human population $({N {}_{h} })$ is divided into three sub-populations: the susceptible $({S {}_{h} })$, the infectious $({I {}_{h} })$, and the recovered $({R {}_{h} })$. The total mosquitoes population $({N {}_{m} })$ is divided into two sub-populations: the susceptible mosquitoes $({S {}_{m} })$ and the infectious mosquitoes $({I {}_{m} })$. The flow diagram of the dengue model is displayed in Figure 1. As the diagram in Figure 1 shows, the dengue model is governed by an autonomous system of ODEs as follows:

with

and

The biologically feasible region of the model (2.1) is set out by

where

and

The description of the variables and parameters of the dengue model can be set out in Tables 1 and 2, respectively.

The change rate of the total human population can be represented by

The total of the mosquito population is constant. Hence, the system (2.1) can be reduced into the following system:

2.1.   Estimation of the parameters model

In this present section, we devote the parameters estimation of the model (2.3) using the fitting least-squares method ^[24,25]. The monthly cumulative data of dengue in East Java Province, Indonesia, from 2018 until 2020 ^[26] is utilized in this study. The parameters ${\mu }_h$, ${\mu }_m$, and $r$ are estimated taking the demographic information. The natural death rate of human $\left({\mu }_h\right)$ is observed from inverse of the average lifespan in East Java Province. The average of lifespan in East Java Province from 2018 until 2020 is 71.15 years. Therefore,

per year ^[27,28]. Next, the death rate of mosquito $\left({\mu }_m\right)$, is found from the inverse of the average lifespan of Aedes mosquitoes. Based on ^[29], the average of lifespan of Aedes mosquitoes is 25 days, that is

per day. Furthermore, for the growth rate of human $\left(r\right)$, is obtained as follows. The average population in East Java Province from 2018 to 2020 is $39,955,059$ ^[30]. Hence,

which is the total human population without the disease, so that $r =$ 561,560.91 per year. The other parameters, such as $B, \ \ {\beta }_{mh}, \ \ {\beta }_{hm}, \ \ {\eta }_h, \ \ K,$ and $N_m$, are estimated with the objective function that minimize

where $t_f$ is the end time of the cumulative dengue cases data $I^{data}_i$ and $I_i$, $\left(i = 0, \ 1, \ 2, \dots, \ t_f\right)$ are the cumulative numerical solutions of infected humans from system (2.3). Thus, the initial populations are

while the initial parameter values for estimation are

The lower bound parameter values are $\left({10}^{-3}, \ {10}^{-3}, \ {10}^{-3}, \ {10}^{-3}, \ 1.2{N_h}_0, \ 2N_{h_0}\right)$ and upper bound parameter values are $\left(1, \ 1, \ 1, \ 1, \ 2N_{h_0}, \ 3N_{h_0}\right)$. The result of the estimation and parameter value can be seen in Figure 2 and Table 3.

3.   Analysis model

In this section, we present the equilibria and the local stability of the DHF model.

3.1.   Equilibria of model

Model (2.3) has the dynamic characteristic of model (2.1), so to understand the dynamics behaviour of the model (2.1), we carry out the analysis of the model (2.3). Likewise, the biological domain in both models is the same. Based on ^[23], two equilibria are obtained, namely the disease-free equilibrium (DFE) and the endemic equilibrium. The DFE of the system (2.3) is provided by

The DFE $E_0$ is exist if $r > \mu _{h}$.

Meanwhile, the endemic equilibrium is

with

The endemic equilibrium exists if

or

Using the the next-generation matrix method ^[31], we yield the basic reproduction number as follows:

By using the formulation of $R_0$, the endemic equilibrium exists if

or $R_0 > 1$ and $r\mathrm{ > }\mu_h.$

3.2.   Local stability analysis

The local stability analysis of dengue model (2.3) can be expressed by the following theorem.

Theorem 3.1. The diseases-free equilibrium $E_{0}\left(N_{h}^{*}, S_{h}^{*}, 0, 0\right)$ of the model (2.3) is locally asymptotically stable if

Proof. The evaluation of Jacobian matrix of the system (2.3) at DFE $E_0$, we have

The two eigenvalues of Jacobian matrix (3.2) are $-\mu_h, -\left(r-\mu _{h} \right)$ that are obviously negative. The other two remaining eigenvalues can be found through the following characteristics equation:

for

and

Using Routh-Hurwitz criterion, the roots of the Eq (3.3) are negative if $a_{1}, a_{2} > 0\, .$ The coefficient $a_{1} > 0\,$, while $a_{2} > 0\,$ if

or

So, the theorem is proven. □

The endemic equilibrium E${}_{1}$ in model (2.3) will be asymptotically stable if

By Substituting the endemic equilibrium $E_1$ to Jacobian matrix of the system (2.3), we find that the eigenvalue is $-\left(r-\mu _{h} \right)$ and the others are the roots of the equation

Equation (3.4) can be expressed by

with

According to the Ruth-Hurwitz criteria, the roots of the characteristic equation will be negative or have a negative real part if

and

The coefficients b${}_{1}$–b${}_{3}$ containing $\tilde{N}_{h}, \tilde{S}_{h}, \tilde{I}_{h}, \tilde{I}_{m}$ so showing that the condition for the root must be positive will be difficult to do analytically. Therefore, the numerical simulation is implemented to analyze the stability of the endemic equilibrium. The parameter values used are given in Table 3.

The phase plane graph between the infected mosquito population and the infected human population is presented in Figure 3. Based on Figure 3, it can be seen that the three graphs converges to the endemic equilibrium E${}_{1}$ as time evolves when

3.3.   Sensitivity analysis

In this work, we discuss the sensitivity analysis to recognize the parameter that can influence the threshold $R_0$. To assign the sensitivity index, we refer to ^[32]. The calculation of the sensitivity index $R_0$ towards to some parameter, say $a$, is given by

Using the parameter values in Table 3, the sensitivity indexes of $R_0$ are summarized in Table 4.

The interpretation of the sensitivity index in Table 4 can be explained as follows. The positive sign reveals that when the parameter values are raised, the value of $R_0$ will also increase. Conversely, the negative sign reveals that when the parameter's value is raised, the value of $R_0$ will be decreased. For example, for

increasing the value of the biting rate of mosquito by $10\%$, causes $R_0$ to increase by $10\%$. Thus, for

increasing the carrying capacity by $10\%$, causes $R_0$ to decrease by $5\%$. Likewise, for the other indexes. The highly sensitive parameters should be considered carefully because small variations in these parameters will cause large quantitative changes. Table 4 shows that the parameter with the greatest influence on the threshold is parameter $B$. The implication is that an increase in the value of the biting rate of mosquitoes increases the spread of dengue disease in the population. Hence, we can anticipate it with preventive control by using a bed net and mosquito repellent lotion or fumigation control to reduce the mosquito population. Furthermore, the value of the biting rate of mosquitoes will decrease, which implies dengue transmission in the community also decreases.

Based on the contour plot in Figures 4 and 5, the basic reproduction number $\left(R_0\right)$ will increase in proportion to the results of the sensitivity analysis in Table 4. The parameters $B$ and ${\beta }_{mh}$ have a positive relation; however the parameters ${\mu }_m$ and ${\eta }_h$ have a negative relation.

4.   Optimal control

In this section, we apply the optimal control (OC) strategy to the dengue model. Based on the sensitivity analysis in the previous section, we can reduce the number of mosquito biting rates with fumigation and prevention efforts. Therefore, we define two control variables, namely $u_1$ as fumigation and $u_2$ as prevention. The dengue model with control variables can be written as follows

The control variables $u_1$ and $u_2$ are established on interval $\left[0, t_f\right]$, where

and $t_f$ represent the final time of the controls. We aim to minimize the number of infected human populations and mosquito populations and keep the cost of employing fumigation and prevention controls as low as possible. For this, we state the objective function

where $A_1$ and $A_2$ are weights of the objective function for $I_h$ and $I_m$, respectively, and $c_1$, and $c_2$ are weight parameters for fumigation and prevention respectively. We utilize the quadratic cost function for $J$ to depict the cost of control efforts. This quadratic function can explain a nonlinear cost increase associated with the performance of control attempts in the field ^[33,34].

Then, to solve this optimal control problem, we assumed that $N_h$ is constant. Using Pontryagin's maximum principle ^[35], the optimal solutions of $u_1$ and $u_2$ are provided by

The variables ${\lambda }_i$, $i = 1, 2, 3, 4, 5$, are represented as adjoint variables or co-state which satisfies the following equations

where the transversality conditions

In this optimal control problem, because $N_h$ is assumed to be constant,

as $N_h$ is the population in East Java Province in 2018. We assume that the values of the parameters are given by

and

and the other parameters are referred to Table 3. We solve the numerical optimal control simulation using the backward and forward iteration technique as stated in ^[36]. To determine which strategy or combination provides useful methods of controlling dengue, we addressed the following strategies for our simulation result, as enumerated below.

4.1.   Strategy 1. Single intervention: use of fumigation only

The result of the simulation of comparison $I_h$ and $I_m$ without and with control when fumigation usage only is given in Figure 6. We can observe that using this strategy, the infections of humans and mosquitoes tend to decrease compared to no control. The profile of the optimal control $u_1$ only is presented in Figure 7. It is apparent that fumigation control should be carried out intensively for a maximum of 50 months.

4.2.   Strategy 2. Single intervention: use of prevention only

The result of the simulation of comparison $I_h$ and $I_m$ without and with control when prevention usage only is given in Figure 8. From Figure 8, it can be seen that both infection of humans and mosquitoes significantly reduce using this strategy. The profile of the optimal control $u_2$ only is displayed in Figure 9. This figure shows that prevention control should be provided with full effort for almost 50 months.

4.3.   Strategy 3. Double intervention: combination of fumigation and prevention

The result of the simulation of comparison $I_h$ and $I_m$ without and with control when implementing both fumigation and prevention are given in Figure 10.

Using the last strategy, it is found that both infectious humans and mosquitoes diminish over to compare to without controls. The profile of the optimal control $u_1$ and $u_2$ simultaneously is demonstrated in Figure 11. As depicted in Figure 11, we found that fumigation and preventive control were both maintaining at maximum levels for nearly 50 months before decaying gradually to the lower bound.

Furthermore, to evaluate the most proper strategy, we assess and compare the merit and the costs related to the control measures by adopting the infection averted ratio (IAR), average cost-effectiveness ratio (ACER) and incremental cost-effectiveness ratio (ICER). These economic verify are mathematically established in ^[37,38] as follows:

The cumulative infected averted is presented by the sum of the difference between the total infectious individuals without and with control. More explicitly

where the symbol with subscript ${}^{*}$ is employed to represent the optimal solutions associated with the corresponding strategy. The most powerful strategy when taking IAR is the strategy with the highest ratio

It proposes a single intervention and evaluates it against its baseline desire. The aim of fumigation control $\left(u_1\right)$ is the reduction of mosquitoes, while prevention control $\left(u_2\right)$ is to protect susceptible humans by diminishing the contact rate. Therefore, we deal with

where $C_i$ represents the cost of the 2 possible interventions. The most cost-effective strategy when using ACER is the strategy with the smallest ratio.

Hence, ICER is utilized to compare two different strategies, namely $i$ and $j$. We calculate the cost of a strategy as undertaken for ACER. When comparing two or more competing intervention strategies incrementally, one intervention is compared with the next effective alternative in increasing order of total infected averted ^[37]. In computing ACER and ICER, we address the same cost for both interventions

The result of the calculation of the total infection and total cost is summarised in Table 5, while IAR and ACER is set out in Table 6.

Based on Table 5, Strategy 2 averts the greatest number of infections and yields the largest ratio

Hence, it is the most effective. This strategy is also the most cost-effective because it is associated with the smallest average cost-effectiveness ratio

Using the total infection and total cost in Table 5, the ICER indexes, as summarised in Table 7, are computed as follows:

Comparing Strategies 1 and 3, the utilization of Strategy 1 is cost-saving over Strategy 3. This means that Strategy 3 is less effective and more costly than the other strategy. Hence, Strategy 3 is eliminated. Furthermore, we recalculate the index of ICER as follows:

Comparing Strategies 1 and 2, the application of Strategy 2 is cost-saving over Strategy 1. This indicates that Strategy 1 is less effective and more costly than the other strategy. Hence, Strategy 1 is removed. Our results recommend that Strategy 2 is the most cost-effective intervention associated with the ICER.

5.   Conclusions

This study has investigated the dengue model by incorporating the logistic growth on the recruitment rate in the human population. The model analyses exhibited that the disease-free equilibrium is locally stable when the reproduction number is less than one, while the endemic equilibrium tends to be asymptotically stable when the threshold is greater than one. The parameters of the dengue model are estimated using data on cases of dengue in East Java Province, Indonesia, from 2018 to 2020. The findings indicate that the spread of dengue is persistent in the population.

Next, the sensitivity analysis shows that mosquito biting rate is the most influential parameter in the spread of dengue. Therefore, we carried out the control strategy in the form of fumigation and prevention. Pontryagin's maximum principle was adopted to accomplish the optimal control that minimizes the spread of dengue. The results of the optimal control problem show that the spread of dengue can be controlled by implementing sustainable control over a short period of time, which shows that the optimal control strategy is effective on humans and mosquitoes. Based on the cost-effectiveness evaluation, prevention is the most effective step to reduce dengue transmission using IAR, ACER, and ICER cost analysis.

Use of Generative-AI tools declaration

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

Acknowledgments

This research has been supported by the Ministry of Research, Technology and Higher Education, Republic of Indonesia, through the PDD 2020.

Conflict of interest

All authors declare that they have no conflicts of interest in this paper.