Abstract

An existing model is extended to assess the impact of some antimalaria control measures, by re-formulating the model as an optimal control problem. This paper investigates the fundamental role of three type of controls, personal protection, treatment, and mosquito reduction strategies in controlling the malaria. We work in the nonlinear optimal control framework. The existence and the uniqueness results of the solution are discussed. A characterization of the optimal control via adjoint variables is established. The optimality system is solved numerically by a competitive Gauss-Seidel-like implicit difference method. Finally, numerical simulations of the optimal control problem, using a set of reasonable parameter values, are carried out to investigate the effectiveness of the proposed control measures.

1. Introduction

Malaria is a tropical infectious disease, transmitted from infected person to susceptible one through Anopheles female mosquito each time the infected mosquito takes a blood meal, and it is prevalent mainly in Africa and some parts of Asia [1]. Clinical symptoms such as fever, pain, chills, and sweats may develop a few days after an infected mosquito bite. Malaria is the cause of the death of over 1 million people each year and many more are infected with the disease. Scientists of several countries are trying to create an effective vaccine to prevent malaria, but it has been difficult. Yet it still poses a major problem throughout much of the world. One of the most important public health programs in many countries is the program to control or to eliminate malaria. This is because malaria is regarded as a very dangerous disease that may lead to death [2].

In the absence of an effective vaccine [3], current control programs for malaria have focused on personal protection and mosquito reduction strategies (such as larviciding, adulticiding and elimination of mosquito breeding sites). A way to prevent the malaria epidemic is to control the growth of mosquito population. Intensive use of insecticides has been one of the major efforts in many years. Most of the mosquitoes use favorable climatic conditions to flourish [4]. Thus, combating efforts of malaria are more effective and economical if it is in phase with climatic changes. Consequently, we consider an optimal control model with three time-dependent controls, prevention ??1, treatment ??2 and mosquito reduction ??3, respectively.

The use of mathematical modeling is increasing influence the theory and practice of disease transmission and control. The mathematical modeling can help in figuring out decisions that are of significant importance on the outcomes and provide complete examinations that enter into decisions in a way that human reasoning and debate cannot. Several health reports and studies in the literature address that malaria is increasing in severity, causing significant public health and socioeconomic burden [5, 6]. Malaria remains the world’s most prevalent vector-borne disease. Despite decades of global eradication and control efforts, the disease is reemerging in areas where control efforts were once effective and emerging in areas thought free of the disease. The global spread necessitates a concerted global effort to combat the spread of malaria. The present study illustrates the use of mathematical modeling and analysis to gain insight into the transmission dynamics of malaria in a population, with main objective on determining optimal control measures. In order to manage the disease, one needs to understand the dynamics of the spread of the disease. Some health scientists have tried to obtain some insight in the transmission and elimination of malaria using mathematical modeling [7].

Recently, some theoretical studies and mathematical models have used control theory [8]. The optimal control efforts are carried out to limit the spread of the disease, and in some cases, to prevent the emergence of drug resistance. The authors in [9] studied a model to assess the impact of some anti-WNV control measures, by reformulating the model as an optimal control problem with density-dependent demographic parameters. They have used two control functions, one for mosquito-reduction strategies and the other for personal (human) protection. In [10], time-dependent prevention and treatment efforts are investigated, where optimal control theory is applied. Using analytical and numerical techniques, they have shown that there are cost-effective control efforts for treatment of hosts and prevention of host-vector contacts. Three types of control functions, one for vector-reduction strategies and the other two for personal (human) protection and blood screening, respectively, are considered in [11]. They have investigated that there are cost-effective control efforts for prevention of direct and indirect transmission of vector-borne disease. In this paper, we use three control functions one for mosquito-reduction strategies and other two for prevention and treatment efforts, respectively. The goal of this paper is to use optimal control theory to evaluate the effectiveness of the control functions. We want to minimize the exposed, infected human and susceptible, exposed, infected mosquito populations with minimum implementation cost.

The paper is organized as follows. Section 2 describes a mathematical model of malaria with three control terms. The analysis of optimization problem is presented in Section 3. The numerical implementation and the strategy used to solve the problem is given in Section 4. Finally, the conclusions are summarized in Section 5.

2. Malaria Model with Controls

The model presented in this section will be a continuation of ideas from recent vector-host models [911]. We will begin with the presentation of the optimal control problem for the transmission dynamics of malaria in order to derive optimal prevention, treatment and mosquito reduction strategies with minimal implementation cost. Our aim is to show that it is possible to implement time-dependent anti-malaria control techniques while minimizing the cost of implementation of such measures. In this paper, we introduce three control functions, ??1, ??2 and ??3. The control functions ??1, ??2 represent time-dependent successful efforts of prevention and treatment, respectively. The well-known practices of prevention efforts include, surveillance, use of mosquito nets, treating mosquito-breading ground and reducing contacts between human and mosquitoes. On the other hand, treatment efforts are carried out by screening patients, administering drug intake. The control function ??3(??) represents the level of larvacide and adulticide used for mosquito control administered at mosquito breeding sites.

We consider a compartmental model that divide the human and mosquito populations into different classes. For humans, the four compartments represent the total population of humans at a given time ??. Let ??h denotes the number of members of a population susceptible to the disease, ??h is the number of infective members of a population, ??h is the number of exposed members of a population and ??h is the number of members who have been removed from the possibility of infection. For mosquitoes the three compartments represent susceptible mosquitoes ????, exposed mosquitoes ???? and infectious mosquitoes ???? populations at a given time ??. The immune class in the mosquito population does not exist, since the mosquito once infected never recover.

For the human population: ?h is the human input rate of new individuals entering the population (assumed susceptible). ??h and ??h are the natural and disease induced death rates respectively. ??h is the progression rate of the exposed class to infectious class. ?? is treatment rate, ??0 and ??0 are rate constants. ????h is the inoculation rate, where ??h is the probability that a bite by an infectious mosquito results in transmission of the disease to the susceptible human and ?? is the contact rate between the two. In the model, the term ????h??h???? denotes the rate at which the human hosts ??h get infected by infected mosquitoes ????.

For the mosquito population: ??? represent the per capita birth rate of mosquito and ???? is the natural death rate in mosquito and ???? is the progression rate of the exposed class in mosquito to infectious class. ?????? is the rate of transmission where ???? is the probability that a bite results in transmission of the parasite to a susceptible mosquito. The term ????????????h refers to the rate at which the susceptible mosquito ???? are infected by the infected human hosts ??h. Parameter definitions and assumptions lead to the following model which involves a system of coupled nonlinear differential equations and three controls.????h????=?h-????h??h(??)?????(??)1-??1?(??)-??h??h(??),????h????=????h??h(??)?????(??)1-??1?(??)-??h??h-??h??h,????h????=??h??h-???+??0??2???h???(??)-h+??h???h(??),????h=???????+??0??2???h(??)-??h??h(??),??????????=????????1-??1?-??????????(??)??h?(??)1-??1?(??)-????????(??)-??0??3(??)????(??),??????????=??????????(??)??h?1-??1?(??)-????????(??)-????????(??)-??0??3(??)????(??),??????????=????????(??)-????????(??)-??0??3(??)????(??),(2.1) with initial conditions given at ??=0. The associated force of infections is reduced by factor of (1-??1(??)), where ??1(??) measures the level of successful prevention (personal protection) efforts. The control variable ??1(??) represents the use of drugs or vaccine which are alternative preventive measures to minimize or eliminate mosquito human contacts (such as the use of insect repellents). The per capita recovery rate is proportional to ??2(??), where ??0>0 is a rate constant. Finally, we describe the role of the third control variable ??3(??). Most of the mosquito use favorable climatic conditions to flourish, particularly during hot and wet seasons. These problems are less pressing during cold seasons. Therefore, we can use a time-dependent mosquito control, preferably applied in seasons favorable for mosquito outbreak. The control variable ??3(??) represents the level of larvicide and adulticide used for mosquito control administered at mosquito breeding sites to eliminate specific breeding areas. It follows that the reproduction rate of the mosquito population is reduced by a factor of (1-??3(??)). It is assumed that under the successful control efforts the mortality rate of mosquito population increases at a rate proportional to ??3(??), where ??0>0 is a rate constant. The impact of the controls is explored via simulations.

3. Mathematical Analysis of the Model

3.1. Positivity and Boundedness of Solutions

Since the model (2.1) represent human and vector population, it is important to prove that all solutions with nonnegative initial data will remain nonnegative for all time.

Theorem 3.1. If ??h(0), ??h(0), ??h(0), ??h(0), ????(0), ????(0), ????(0) are nonnegative, then ??h(??), ??h(??), ??h(??), ??h(??), ????(??), ????(??), ????(??) remain nonnegative for all ??>0.

Proof. To see this, let ??*?=sup??>0:??h>0,??h=0,??h>0,??h=0,????>0,????=0,?????=0.(3.1) Thus, ??*>0. Then, from the first equation of the system (2.1) we have ????h????=?h-??h??h?????1-??1?-??h??h,=?h-???h?????1-??1?+??h???h.(3.2) Letting ??(??)=??h??h????(1-??1), the above equation can be written as ?????????h??exp??0??(??)????+??h????=?h??exp??0??(??)????+??h???,(3.3) integrating both sides from ??=0 to ??=??*, ??h???*???exp??*0??(??)????+??h??*?-??h?(0)=??*0?h??exp??0??(??)????+??h???????,(3.4) multiplying both sides by ?exp{-??*0??(??)????-??h??*}, ??h???*?=??h?-?(0)exp??*0??(??)????-??h??*??-?+exp??*0??(??)????-??h??*?×???*0?h??exp??0??(??)????+??h???????>0.(3.5)
Thus, ??h(??*) being the sum of positive terms is positive. By the same argument, it can be proved that the quantities ??h,??h,??h,????,????, and ???? are positive for all times ??>0. This completes the proof.

The objective functional F formulates the optimization problem of interest, namely, that of identifying the most effective strategies. The overall preselected objective involves the minimization of the number of exposed, infected human and minimizing the susceptible, exposed and infective mosquitoes at a minimal cost over a finite time interval [0,??].

Define the objective functional F,F???1,??2,??3?=???0???1??h+??2??h+??3????+12???1??21+??2??22+??3??23??????.(3.6) The goal is to minimize the cost functional (3.6). This functional includes the number of exposed, infectious and the total number of mosquito population, respectively, as well as the social costs related to the resources needed for, personal protection ??1??21, treatment ??2??22, and spraying of insecticides operations, ??3??23. In words, we are minimizing the number of exposed, infectious human and susceptible, exposed and infectious mosquito populations as well as the cost based on the implementation of the control functions. We choose to model the control efforts via a linear combination of quadratic terms, ??2??(??)(??=1,2,3). The constants ???? and ????(??=1,2,3) represent a measure of the relative cost of the interventions over [0,??]. The objective of the optimal control problem is to seek optimal control functions (??*1(??),??*2(??),??*3(??)) such that F???*1,??*2,??*3??F???=min1,??2,??3?,???1,??2,??3?????,(3.7)

where the control set is defined as ??????=??=1,??2,??3?|????(??)isLebesguemeasurable,0=????[]?(??)=1,???0,??for??=1,2,3(3.8) subject to the system (2.1) and appropriate initial conditions. Pontryagin’s Maximum Principle is used to solve this optimal control problem and the derivation of the necessary conditions. First we prove the existence of an optimal control for system (2.1) and then derive the optimality system.

3.2. Existence of an Optimal Control

Theorem 3.2. Given the objective functional F(??1,??2,??3)=???0(??1??h+??2??h+??3????+1/2(??1??21+??2??22??+????3??23))????, where the control set U given by (3.8) is measurable subject to system (2.1) with initial conditions given at t=0, then there exists an optimal control ??*=(??*1(??),??*2(??),??*3(??)) such that F(??*1,??*2,??*3)=min{F(??1,??2,??3),(??1,??2,??3)???}.

Proof. The integrand of the objective functional F given by (3.6) is a convex function of (??1,??2,??3) and the state system satisfies the Lipschitz property with respect to the state variables since state solutions are bounded. The existence of an optimal control follows [12].
In order to find an optimal solution, first we should find the Lagrangian and Hamiltonian for the optimal control problem (2.1)–(3.6). The Lagrangian of the control problem is given by??=??1??h+??2??h+??3?????+????+?????+12???1??21+??2??22+??2??23?.(3.9) We seek for the minimal value of the Lagrangian. To do this, we define the Hamiltonian function ?? for the system, where ????, ??=1,,7 are the adjoint variables: ??=??1??h+??2??h+??3????+12???1??21+??2??22+??2??23?+??1??h-????h??h(??)?????(??)1-??1?(??)-??h??h?(??)+??2?????h??h(??)?????(??)1-??1?(??)-??h??h-??h??h?+??3???h??h-???+??0??2???h???(??)-h+??h???h?(??)+??4????+??0??2???h(??)-??h??h?(??)+??5?????????1-??3?(??)-??????????(??)??h?(??)1-??1?(??)-????????(??)-??0??3(??)?????(??)+??6???????????(??)??h?1-??1?(??)-????????(??)-????????(??)-??0??3(??)?????(??)+??7?????????(??)-????????(??)-??0??3(??)?????.(??)(3.10)
In order to derive the necessary conditions, we use Pontryagin’s maximum principle [13] as follows.
If (??,??) is an optimal solution of an optimal control problem, then there exists a non trivial vector function ??=(??1,??2,,????) satisfying the following conditions.????=????????(??,??,??,??),????0=????(??,??,??,??),????????=????????(??,??,??,??).????(3.11)

We now derive the necessary conditions that optimal control functions and corresponding states must satisfy. In the following theorem, we present the adjoint system and control characterization.

Theorem 3.3. Given an optimal control ??*=(??*1,??*2,??3*) and a solution ??*=(??*h,??*h,??*h,??*h,??*??,??*??,??*??) of the corresponding state system (2.1), there exists adjoint variables ????, ??=1,,7 satisfying????1(??)????=????h???1-??2??1-??1?????+??h??1,????2(??)????=??h???2-??3?+??h??2-??1,????3(??)=???????+??0??2????3-??4?+???h+??h???3+?????????5-??6??1-??1?????-??2,????4(??)????=??h??4,????5(??)????=-?????5?1-??3?+?????????5-??6??1-??1???h+??????5+??0??5u3-??3,????6(??)????=-?????5?1-??3?+???????6-??7?+??????6+??0??6??3-??3,????7(??)????=-?????5?1-??3?+????h???1-??2??1-??1???h+??????7+??0??7??3-??3,(3.12) with transversality conditions ???????end?=0,,??=1,2,,7.(3.13) Furthermore, the control functions ??*1, ??*2, and ??*3 are given by ??*1????=maxmin1??,??,1,0*2????=maxmin2??,??,1,0*3????=maxmin3??,,1,0(3.14) where ??1=????h???2-??1???*h??*??+?????????6-??5???*????*h??1,??2=???3-??4???0??*h??2,??3=?????5??*??+??0???5??*??+??6??*??+??7??*?????3.(3.15)

Proof. To determine the adjoint equations and the transversality conditions we use the Hamiltonian (3.10). The adjoint system results from Pontryagin’s Maximum Principle [13]. ????1(??)????=-????????h,????2(??)????=-????????h,,????7(??)????=-??????????,(3.16) with ????(??)=0.
To get the characterization of the optimal control given by (3.14), solving the equations,????????1=0,????????2=0,????????3=0,(3.17) on the interior of the control set and using the property of the control space ??, we can derive the desired characterization (3.14).

4. Numerical Results and Discussion

The numerical algorithm presented below is a semi-implicit finite difference method.

We discretize the interval [??0,????] at the points ????=??0+????(??=0,1,,??), where ?? is the time step such that ????=????. Next, we define the state and adjoint variables ??h(??), ??h(??), ??h(??), ??h(??), ????(??), ????(??), ????(??), ??1(??), ??2(??), ??3(??), ??4(??), ??5(??), ??6(??), ??7(??) and the controls ??1(??), ??2(??), ??3(??) in terms of nodal points ????h, ????h, ????h, ????h, ??????, ??????, ??????, ????1, ????2, ????3, ????4, ????5, ????6, ????7, ????1, ????2 and ????3. Now a combination of forward and backward difference approximation is used as follows.

The Method, developed by [14] and presented in [15, 16], to adapt it to our case as following: ??h??+1-????h??=?h-????h??h??+1???????1-????1?-??h??h??+1,??h??+1-????h??=????h??h??+1???????1-????1?-???h+??h???h??+1,??h??+1-????h??=??h??h??+1-???+??0????2???h??+1-???h+??h???h??+1,??h??+1-????h??=???+??0????2???h??+1-??h??h??+1,??????+1-????????=??????????+1+??????+????????1-????3?-????????????+1??h??+1?1-????1?-?????+??0????3???????+1,??????+1-????????=????????????+1??h??+1?1-????1?-?????+????+??0????3???????+1,??????+1-????????=??????????+1-?????+??0????3???????+1.(4.1) By using a similar technique, we approximate the time derivative of the adjoint variables by their first-order backward-difference and we use the appropriated scheme as follows??1??-??-??1??-??-1??=????h???1??-??-1-??2??-????1-????1???????+1+??h??1??-??-1,??2??-??-??2??-??-1??=??h???2??-??-1-??3??-???-??1+??h??2??-??-1,??3??-??-??3??-??-1??=???+??0????2????3??-??-1-??4??-???+???h+??h???3??-??-1+?????????5??-??-??6??-????1-????1???????+1-??2,??4??-??-??4??-??-1??=??h??4??-??-1,??5??-??-??5??-??-1??=-?????5??-??-1?1-????3?+?????????5??-??-1-??6??-????1-????1???h??+1+?????+??0????3???5??-??-1-??3,??6??-??-??6??-??-1??=-?????5??-??-1?1-????3?+???????6??-??-1-??7??-???+?????+??0????3???6??-??-1-??3,??7??-??-??7??-??-1??=-?????5??-??-1?1-????3?+????h???1??-??-1-??2??-??-1??1-????1???h??+1+?????+??0????3???7??-??-1-??3.(4.2) The algorithm describing the approximation method for obtaining the optimal control is the following.

Algorithm 4.1. Step 1. Consider ??h(0)=??h0, ??h(0)=??h0, ??h(0)=??h0, ??h(0)=??h0, ????(0)=S??0, ????(0)=????0, ????(0)=????0?,? ????(????)=0 (??=1,…, 5), and ??1(0)=??2(0)=??3(0)=0.Step 2. For ??=1,,??-1, do??h??+1=????h+???h?1+??????h???????1-????1?+??h?,??h??+1=????h+??????h??h??+1???????1-????1????1+??h+??h?,??h??+1=????h+????h??h??+1?1+????+??0????2+??h+??h?,??h??+1=????h?+????+??0????2???h??+11+????h,??????+1=??????+???h???????+????????1-????3??1+??-????1-????3?+????????h??+1?1-????1?+????+??0????3?,??????+1=??????+??????????????+1??h??+1?1-????1????1+????+????+??0????3?,??????+1=??????+????????????+1???1+????+??0????3?,??1??-??-1=??1??-??+????2??-??????h?1-????1???????+1?1+??????h?1-????1???????+1+??h?,??2??-??-1=??2??-?????+??h??3??-??+??1????1+??h+??h?,??3??-??-1=??3??-??+??????+??0????2???4??-??-?????????5??-??-??6??-????1-????1???????+1+??2??1+????+??0????2+??h+??h?,??4??-??-1=??4??-??1+????h,??5??-??-1=??5??-???+??????????6??-???1-????1???h??+1+??3??1+??-????1-????3?+???????1-????1???h??+1+????+??0????3?,??6??-??-1=??6??-????+??????5??-??-1?1-????3?+??????7??-??+??3????1+????+????+??0????3?,??7??-??-1=??7??-????+??????5??-??-1?1-????3?-????h???1??-??-1-??2??-??-1??1-????1???h??+1+??3????1+????+??0????3?,??1??+1=???2??-??-1-??1??-??-1?????h??h??+1??????+1+???6??-??-1-??5??-??-1?????????????+1??h??+1??1,??2??+1=???3??-??-1-??4??-??-1???0??h??+1??2,??3??+1=?????5??-??-1???????+1+??????+1+??????+1?+??0???5??-??-1??????+1+??6??-??-1??????+1+??7??-??-1??????+1???3,??1??+1????=min1,max1??+1,??,0??2??+1????=min1,max2??+1,??,0??3??+1????=min1,max3??+1,,0??(4.3) end for. Step 3. For??=1,…, ??-1, write ??*h(????)=????h, ??*h(????)=????h, ??*h(????)=????h, ??*h(????)=????h, ??*??(????)=??????, ??*??(????)=??????, ??*??(????)=????h, ??*1(????)=????1, ??*2(????)=????2, ??*3(????)=????3.
end for.

We have plotted susceptible, infected, exposed and recovered human population with and without control by considering parameter values as given in Table 1, with initialization ??h(0)=100, ??h(0)=20, ??h(0)=20, ??h(0)=10. The control individuals are marked by solid blue line while the individual without control are marked by black line. Similarly, we have plotted susceptible, infected, exposed and total number of mosquito population with and without control by considering parameter values as given in Table 1, with initial conditions ????(0)=1000, ????(0)=20, ????(0)=30. The control individuals are marked by solid blue line while the individual without control are marked by black line. The weight constant values in the objective functional are ??1=1, ??2=1, ??3=1, ??1=150, ??2=50 and ??3=300. Figure 1, represents the population of susceptible, exposed, infected and recovered humans with and without control. The solid black line denotes the population of individuals in system (2.1) without control while the blue line denotes the population of exposed individuals in system (2.1) with control. We see that the population of the exposed, infected human with control is more sharply decreased after 12 or 13 days than the individuals without control. Figure 2 represents the population of susceptible, exposed, infected and the total number of mosquito population in the system (2.1) with and without control. The population of susceptible, exposed and infected mosquito with control is more sharply decreased than without control and becomes very small. In Figure 2, we see that if there are no control susceptible, exposed and infected mosquito without control constantly increases, but if there is control the population of susceptible, exposed and infected mosquito begins to decrease from the very beginning day of the control measure. Figures 3, 4, 5, 6 represent the optimal controls ??1, ??2 and ??3, the balancing constants ??1, ??2, ??3, ??1, ??2 and ??3 are given in each of the figure captions and parameters are given in Table 1. Parameter values used in the numerical simulations are estimated based on malaria disease as given in Table 1 and are taken from [10]. The constant ??0 is arbitrarily chosen with ??0=0.06. For illustration purposes, we consider the parameter values in Table 1 for numerical simulation.

5. Conclusion

A comprehensive, continuous model for the transmission dynamics of malaria has been presented. We sought to determine optimal control strategies that would minimize not only the exposed, infected human, and susceptible, exposed, infected mosquito but also the cost of implementation of the control as well. Our model incorporates three control measures. We analyzed the optimal control using the functional F in terms of quadratic forms. Minimizing the cost we obtained the optimal controls ??1, ??2 and ??3 where F was minimized. The model is analyzed for the existence of control. The proposed optimal control shows the result of optimally controlling the disease using three controls, personal protection, treatment and mosquito reduction strategies, respectively. We have developed the necessary conditions for the optimal control using the Pontryagin’s Maximum Principle. Using the state and adjoint system together with the characterization of the optimal control, we solved the problem numerically via an efficient numerical method based on optimal control with parameter values estimated based on malaria. A comparison between optimal control and without control dynamics is presented. It is easy to see that the optimal control has a very desirable effect upon the population for reducing the number of infected individuals. In order to illustrate the overall picture of the disease, the numbers of susceptible, exposed, infected, recovered human population and susceptible, exposed, infected mosquito population under the optimal control and without control are shown in figures.

Acknowledgments

This research was partially supported by the Deanship of Scientific Research, King Khalid University, Saudi Arabia (KKU-SCI-11-025) and the work of Prof. Il Hyo Jung was supported by the Cooperative Research Program for Agriculture Science and Technology Development (Project No. PJ007420) Rural Development Administration, Republic of Korea. A. A. Lashari acknowledges the support from the National University of Sciences and Technology Islamabad, Pakistan.