Transmission Dynamics and Optimal Control of Malaria in Kenya

This paper proposes and analyses a mathematical model for the transmission dynamics of malaria with four-time dependent control measures in Kenya: insecticide treated bed nets (ITNs), treatment, indoor residual spray (IRS), and intermittent preventive treatment of malaria in pregnancy (IPTp). We first considered constant control parameters and calculate the basic reproduction number and investigate existence and stability of equilibria as well as stability analysis. We proved that if R 0 ≤ 1, the disease-free equilibrium is globally asymptotically stable in D. If R 0 > 1, the unique endemic equilibrium exists and is globally asymptotically stable. The model also exhibits backward bifurcation at R 0 = 1. If R 0 > 1, the model admits a unique endemic equilibrium which is globally asymptotically stable in the interior of feasible region D. The sensitivity results showed that the most sensitive parameters are mosquito death rate and mosquito biting rates. We then consider the time-dependent control case and use Pontryagin’s Maximum Principle to derive the necessary conditions for the optimal control of the disease using the proposed model. The existence of optimal control problem is proved. Numerical simulations of the optimal control problem using a set of reasonable parameter values suggest that the optimal control strategy for malaria control in endemic areas is the combined use of treatment and IRS; for epidemic prone areas is the use of treatment and IRS; for seasonal areas is the use of treatment; and for low risk areas is the use of ITNs and treatment. Control programs that follow these strategies can effectively reduce the spread of malaria disease in different malaria transmission settings in Kenya.


Introduction
Malaria is a leading cause of mortality and morbidity among the under-five group and the pregnant women in Sub-Saharan Africa [1].These groups are at high risk due to weakened and immature immunity, respectively.With the recent conversion of the Millennium Development Goals (MDGs) to Sustainable Development Goals (SDGs) as part of Global Malaria Action Plan for a malaria-free world by 2030, reducing malaria is critical to post-2015 malaria strategies and for achieving the SDGs such as ensuring healthy lives and promoting well-being for all at all ages.Most Kenyans are vulnerable to malaria because of poverty, inadequate health care infrastructures, and low income of the country.Prompt access to effective treatment for malaria is unacceptably low in Kenya due to the socioeconomic barriers to accessing health care.These challenges call for urgent development of effective and optimal strategies for preventing and controlling spread of malaria.Malaria transmission is highly variable across Kenya because of the different transmission intensities driven by climate and temperature.Kenya has four malaria epidemiological zones: the endemic areas, the seasonal malaria transmission, the malaria epidemic prone areas, and the low risk malaria areas [2].Malaria is caused by Plasmodium parasites and it is transmitted from one individual to another by the bite of infected female anopheline mosquitoes [1].Malaria is an entirely preventable and treatable disease, provided the currently recommended interventions are properly implemented.Controlling malaria transmission involves interrupting the malaria transmission for specific transmission settings and for the most at-risk groups of malaria.World Health organization (WHO) recommended malaria intervention strategies include the use of long-lasting insecticidetreated bed nets (LLINs), indoor residual spray (IRS), chemoprevention for the most vulnerable such as intermittent preventive treatment for pregnant women (IPTp), confirmation of malaria diagnostics through rapid diagnostics tests (RDTs) and microscopy for every suspected case, and timely treatment with artemisinin-based combination therapies (ACTs) [1,3].
Mathematical modelling has become an important tool in understanding the dynamics of disease transmission and in decision making processes regarding intervention programs for disease control.Mathematical models provide a framework for understanding the transmission dynamics for malaria and can be used for the optimal allocation of different interventions against malaria [4,5].Optimal control is a branch of mathematics developed to find optimal ways of controlling a dynamic system [6].Application of optimal control theory can be an important tool for estimating the efficacy of various policies and control measures and the cost of implementing them.Optimal control approaches have been successfully previously used in decision making for the infectious diseases such as tuberculosis [7,8], HIV [9], and influenza [10].Optimal control theory has also been applied in malaria control to assess the impact of antimalarial control measures by formulating the model as an optimal control problem.Most malaria models for analyzing effect of interventions in optimal control used the standard Susceptible-Exposed-Infectious-Recovered (SEIR) model for humans and Susceptible-Exposed-Infectious (SEI) model for mosquitoes.Okosun et al. [11] considered three control variables in assessing the optimal control and cost effectiveness of the interventions but not for different settings.Mwamtobe et al. [12] used three control variables and for only one region in Malawi.Kim et al. [13] used two control efforts in the optimal control model for malaria transmission in South Korea.Agusto et al. [14] used three system control variables.Silva and Torres [15] used one control variable.
IPTp is one of the WHO recommended prevention therapies for the pregnant women.IPTp has been shown to be effective in reducing maternal and infant mortality that are related to malaria for the most at-risk group for malaria [16][17][18][19].Very few studies have been carried out in applying optimal control theory to malaria transmission models for different transmission settings.The combined effect of ITNs, IRS, and natural death on reducing the mosquito population has not been demonstrated in optimal control theory in malaria control.The effect of IPTp which is WHO recommended preventive therapy for the most at-risk group for malaria (pregnant women) has not been studied in optimal control theory.No optimal control model for malaria interventions for different transmission settings exits for Kenya.No optimal control model for four control variables incorporating the IPTp malaria intervention studies exits for Kenya.No optimal control model has been stratified by the age group (under five) and specific categories (pregnant women).
In this paper, a model for malaria transmission dynamics with four control strategies is formulated and analyzed.We then formulate an optimal control problem and derive expressions for the optimal control for the malaria transmission model with four control variables and then use the optimal control theory to study the effectiveness of all possible combinations of four malaria preventive measures among the pregnant women and children under five years of age.

Model Formulation
The model is formulated by considering the human and mosquito subgroups.The considered model consists of population of susceptible  ℎ , Exposed humans  ℎ , infected humans  ℎ , recovered humans  ℎ , susceptible mosquitoes   , exposed mosquitoes,   and infected mosquitoes   .The total population sizes at time  for humans and mosquitoes are denoted by  ℎ () and   (), respectively.We employ the SEIRS type model for humans to describe a disease with temporary immunity on recovery from infection.Mosquitoes are assumed not to recover from the parasites so the mosquito population can be described by the SEI model.In the model we incorporate four time-dependent control measures simultaneously: (i) the use of treated bed nets  1 (), (ii) treatment of infective humans  2 (), (iii) spray of insecticides  3 (), and (iv) treatment to protect pregnant women and their newborn children: intermittent preventive treatment for pregnant women (IPTp)  4 (). ℎ () represents the number of individuals not yet infected with the malaria parasite at time ,  ℎ () represents individuals who are infected but not yet infectious,  ℎ () is the class representing those who are infected with malaria and are capable of transmitting the disease to susceptible mosquitoes, and  ℎ () represents the class of individuals who have temporarily recovered from the disease.
Figure 1 describes the dynamics of malaria in human and mosquito populations together with interventions.
The susceptible humans (pregnant women and children under the age of five) ( ℎ ) are recruited at the rate Λ ℎ .They either die from natural causes (at a rate  ℎ ) or move to the exposed class ( ℎ ) by acquiring malaria through contact with infectious mosquitoes at a rate (1 −  1 )(  / ℎ ) ℎ or (1 −  4 )(  / ℎ ) ℎ , where  is the transmission probability per bite,  is the per capita biting rate of mosquitoes,  is the contact rate of vector per human per unit time,  1 () ∈ [0, 1] is the preventive measure using ITNs,  4 () ∈ [0, 1] is the preventive measure using IPTp,   () is the infectious mosquitoes at time ,  ℎ () is the total number of individuals (pregnant women and children under the age of five), and  ℎ () is the total number of pregnant women.Susceptible class  ℎ is divided into whole population (children under the age of five and pregnant women) being exposed and the population for the pregnant women being exposed.Exposed individuals move to the infectious class after the development of clinical symptoms at the rate  ℎ .Infectious individuals are assumed to recover at a rate  +  2 , where  is the rate of spontaneous recovery,  2 () ∈ [0, 1] is the control on treatment of infected individuals, and  ∈ [0, 1] is the efficacy of treatment.Infectious individuals who do not recover die at a rate  ℎ +  ℎ .Individuals infected with malaria suffer a disease induced death at rate of  ℎ and natural death  ℎ .Infected individuals then progress to partially immune group where upon recovery the partially immune individual loses immunity at the rate  and becomes susceptible again.
Susceptible mosquitoes (  ) are recruited at the rate Λ  and acquire malaria infection (following contact with humans infected with malaria) at the rate   .They either die from natural causes (at a rate   ) or move to the exposed class by acquiring malaria through contacts with infected humans at a rate (1 −   1. Table 2 describes and shows parameters of the model.Table 3 describes and represents malaria prevention and control strategies practiced in Kenya. The following assumptions have been used in the formulation of the model: (i) Population for human and mosquito being constant (no immigrants).
(ii) No recovery for infected mosquitoes.
(iii) Mosquitoes not dying due to disease infection.
(iv) All parameters in the model being nonnegative.
Putting the above formulations and assumptions together gives the following system of nonlinear differential equations describing the dynamics of malaria in human and mosquito populations together with interventions:  ( With initial conditions =  ℎ / ℎ is the per capita incidence rate among mosquitoes (force of infection for susceptible mosquitoes),  ℎ =   / ℎ is the force of infection for susceptible humans (pregnant and under 5), and  ℎ =   / ℎ is the force of infection for susceptible pregnant humans.The total population size for the human is  ℎ =  ℎ +  ℎ +  ℎ +  ℎ and for mosquito is   =   +   +   and their differential equations are given by  ℎ / = Λ ℎ − ℎ  ℎ − ℎ  ℎ and   / = Λ  − (  +  1 +  3 )  , respectively.

Mathematical Analysis of the Malaria Model with
Intervention Strategies.We will assume that the control parameters are constant so as to determine the basic reproduction number, steady states, and their stability as well as the bifurcation analysis.
We describe the basic properties and analysis of the formulated malaria model with control strategies through mathematical analysis of the formulated model.

Basic Properties of the Model: Positivity and Invariant
Regions.All the state variables and parameters for model (1) are nonnegative for all  ≥ 0.
The feasible solutions set for model (1) given by is positively invariant and hence model ( 1) is biologically, epidemiologically meaningful and mathematically well posed in the domain .System (1) has always a disease-free equilibrium given by , 0, 0) . (4)

Basic Reproduction Number.
The matrices  and  for the new infection terms and the remaining transfer terms at disease-free equilibrium [26], respectively, are given by It follows that the basic reproduction of model (1), denoted by  0 , is given by

Stability Analysis of Disease-Free Equilibrium Point (1) Local Stability of Disease-Free Equilibrium Point
Theorem 1.The disease-free equilibrium point for system (1) is locally asymptotically stable if  0 < 1.
Proof.The Jacobian matrix () of the malaria model (1) at the disease-free equilibrium point is given by The eigenvalues of the Jacobian matrix are the solutions of the characteristic equation Expanding the determinant into a characteristic equation we have Hence we have where Thus, applying the Routh-Hurwitz criteria [27] to polynomial (10), we have that all determinants of the Hurwitz matrices are positive.Hence all the eigenvalues of the Jacobian have negative real part, implying that the DFE point is (at least) locally asymptotically stable ( 0 < 1).
Next, we study the global behavior of the disease-free equilibrium for system (1).
Next, we investigate the endemic equilibrium and its stability of system (1).

Stability Analysis of Endemic Equilibrium
Point.First we determine the existence of the endemic equilibrium points.
The endemic equilibrium ( 1 ) of the model is given by To derive  1 , we solve model ( 1) by equating it to zero.Substituting and solving for  * * ℎ (as an expression of parameters only) through some algebraic manipulation give or where We use the quadratic formula to find the roots of (17); that is, Hence From the quadratic equation ( 17) we analyze the possibility of multiple endemic equilibria.It is important to note that the coefficient  is always positive with  and  having different signs.All the negative terms in  are in the form of (1 − ) where  is the control that is bounded by 1.It follows that (i) there is a unique endemic equilibrium if  < 0 and  = 0 or  2 − 4 = 0, (ii) there is a unique endemic equilibrium if  < 0 (i.e., if  0 > 1), (iii) there are two endemic equilibria if  > 0,  < 0, and  2 − 4 > 0, (iv) there are no endemic equilibria otherwise.
Note that the hypothesis  > 0 is equivalent to  0 < 1.
Thus the results of this section can be summarized in the following theorem.Theorem 3. If  0 < 1,  0 is an equilibrium of system (1) and it is locally asymptotically stable.Furthermore, there exists an endemic equilibrium if conditions in (i) are satisfied or two endemic equilibria if conditions in (iii) are satisfied.If  0 > 1, then  0 is unstable and there exists a unique endemic equilibrium.
Item (iii) indicates the possibility of backward bifurcation in model (1) when  0 < 1.In the next section we will prove the occurrence of multiple equilibria for  0 < 1 comes from the backward bifurcation and this will give information on the local stability of the endemic equilibrium.
(1) Local Stability Analysis of Endemic Equilibrium Point.We use centre manifold theory [29] to investigate the stability of the endemic equilibria for model (1) where we carry out bifurcation analysis of system (1) at  0 = 1.We consider a transmission rate  as bifurcation parameter Ψ so that  0 = 1.
To apply the theory, we introduce dimensionless state variables into system (1).
The system of (1) can be written as We let and  7 =   .
That is, The Jacobian matrix of (1) calculated at Ψ * is given by ) .
Centre manifold approach [29] is then applied.
A right eigenvector associated with the eigenvalue zero is  = ( 1 ,  2 , . . .,  7 ).Solving the system we have the following right eigenvector: The left eigenvectors satisfying Solving the system, the left eigenvector will be given by , Computing for the sign of  and  as indicated in the theorem gives so that  is always positive.Therefore the following result is established.
Finally, we will investigate the global stability of the endemic equilibrium in the feasible region.
(2) Global Stability Analysis of Endemic Equilibrium Point.Global stability results for the endemic equilibrium can be obtained when it is unique and whenever it exists.We have established in Theorem 4 that if  0 > 1 this implies the existence and uniqueness of the endemic equilibrium.
The global behavior of the endemic equilibrium of model (1) when it exists is explored by proving that such an equilibrium is globally asymptotic stable in the interior of the feasible region .We will use the geometric approach to global stability as described by Li and Muldowney [30].The following conditions are required for the global stability of the endemic equilibrium,  1 : (i) the uniqueness of  1 in the interior of  (condition  1 ); (ii) the existence of an absorbing compact set in the interior of  (condition  2 ); and (iii) the fulfillment of a Bendixson criterion (i.e., inequality (A.6)).
Theorem 5.If  < 0 and  > 0 and  0 > 1, then the endemic equilibrium of the malaria model ( 1) is globally asymptotically stable in the interior of .
Proof.Following Li and Muldowney [30], for system (1), under the assumption of  0 > 1, satisfies conditions ( 1 )-( 2 ), the existence of the endemic equilibrium has also been shown, and the instability of DFE implies the uniform persistence [31]; that is, there exists a constant  > 0 such that any solutions ( ℎ (),  ℎ (),   ()) with ( ℎ (0),  ℎ (0),   (0)) in the interior of  satisfy min { lim The uniform persistence together with boundedness of  is equivalent to the existence of a compact set in the interior of  which is absorbing for (4) as described by Hutson and Schmitt [32].Thus, ( 1 ) is verified.Moreover,  1 is the only equilibrium in the interior of , so that ( 2 ) is also verified.
What remains is to find conditions for which the Bendixson criterion given by (A.6) is verified.To this aim, let us begin by observing that, from the Jacobian matrix associated with a general solution ( ℎ ,  ℎ ,   ) of reduced system (1), we get the second additive compound matrix  |2| : where where and the matrix  =    −1 +  |2|  −1 can be written in block form as where The vector norm | ⋅ | in R Therefore Therefore We rewrite the last two equations of system (1) for İℎ and İ as Substituting ( 39) into ( 37) and ( 40) into (38) we have For the uniform persistence constant  > 0, there exists a time  0 > 0 independent of (0) ∈ , the compact absorbing set, such that for  >  0 we have Relations (41) imply where Along each solution ( ℎ (),  ℎ (),   ()) to ( 1) with ( ℎ (0),  ℎ (0),   (0)) ∈  where  is the compact absorbing set, we have, for  >  0 , Which implies  2 < /2 < 0. This proves that the unique endemic equilibrium is globally asymptotically stable whenever it exist, thus completing the proof.

Sensitivity Analysis.
Sensitivity analysis is to assess the relative impact of each of the parameters of the basic reproductive number.The normalized forward sensitivity index of the reproduction number with respect to these parameters given in Table 4 is computed.The index measures the relative change in a variable with respect to relative changes in parameters.Definition 6.Following Chitnis et al. [33], the normalized forward sensitivity index of a variable, ℎ, that depends on a parameter, , is defined as   = (/ℎ) × (ℎ/).
The sensitivity index of the model parameters is given by Sensitivity indices for the control parameters are given by Through sensitivity analysis, it is observed that the most sensitive parameters to  0 across all the settings were the mosquito's natural death rate,   , and mosquito biting rate, ; this was followed by the by the mosquito contact rate with humans, , probability of mosquito getting infected, , the probability of humans getting infected, , and the recruitment rate of mosquitoes and humans (see Table 4).
Sensitivity Indices of  0 .Sensitivity indices are calculated using parameters in Table 5.

Analysis of Optimal Control of the Malaria Model with Intervention
Strategies.We consider the case of timedependent control variables.The malaria dynamics model is extended and an optimal control problem is formulated.We formulate an optimal control model for malaria disease in order to determine optimal malaria control strategies (ITNs, IRS, IPTp, and treatment) with minimal implementation cost.
For the optimal control problem of the given system, we consider the control variables () = ( 1 ,  2 ,  3 ,  4 ) ∈  relative to the state variables  ℎ ,  ℎ ,  ℎ ,  ℎ ,   ,   ,   where control variables are bounded and measured with We define the objective function as subject to Our aim with the given objective function is to minimize total number of mosquito population   , the number of exposed humans  ℎ (), and infected humans  ℎ () while minimizing the cost of control  1 (),  2 (),  3 (), and  4 ().We select to model the control efforts via a linear combination of quadratic terms and the constants which represents a measure of the relative cost of the interventions over [0, 1].
We seek an optimal control  *

Existence of Optimal Control
Problem.The existence of an optimal control can be proved by using the theorem given in Fleming and Rishel [34].It can be clearly seen that the system of equations given by ( 1) is bounded from above by a linear system.According to the result in Fleming and Rishel [34], the solution exists if the following hypotheses are met: (H In order to confirm the above hypotheses, a result given by Lukes [35] is used to establish the existence of solutions of state system (1).Since the coefficients are bounded, ( This concludes existence of an optimal control since the state variables are bounded.Hence we have the following theorem.(1) with initial conditions, then there exists an optimal control

Theorem 7. Given the objective functional
Lagrangian for a problem discusses how the techniques come and Hamiltonian helps in solving for the adjoint variable.In order to find an optimal solution, first we find the Lagrangian and Hamiltonian for the optimal control problem (51).The Lagrangian of the optimal problem is given by We seek to find the minimal value of the Lagrangian.To do this, we define the Hamiltonian  for the control problem which consists of the integrand of the objective functional (Lagrangian, ) and the inner product of the right hand sides of the state equations and the costate variables or adjoint variables ( 1 ,  2 ,  3 ,  4 ,  5 ,  6 ,  7 ) as Taking  = ( ℎ ,  ℎ ,  ℎ ,  ℎ ,   ,   ,   ),  = ( 1 ,  2 ,  3 ,  4 ), and  = ( 1 ,  2 ,  3 ,  4 ,  5 ,  6 ,  7 ) we obtain the Hamiltonian given by 3.3.The Optimality System.In order to find the necessary conditions for this optimal control, we apply Pontryagin's Maximum Principle [6] as described by Lenhart and Workman [36] as follows.
The state equation is The optimality condition is and the adjoint equation is Now, we apply the necessary conditions to the Hamiltonian.
with transversality conditions )} , Proof.To determine the adjoint equations and the transversality conditions we use the Hamiltonian .The Hamiltonian function, , is differentiated with respect to  ℎ ,  ℎ ,  ℎ ,  ℎ ,   ,   , and   .The adjoint/costate equation is given by with transversality conditions In order to minimize Hamiltonian, , with respect to the controls at the optimal controls,  is differentiated with respect to  1 ,  2 ,  3 , and  4 on the set , and the solution for the optimal control point is obtained after equating to zero.This is the optimality condition.
Hence the state and optimal control can be calculated using the optimality system.Hence using the fact that the second derivatives of the Lagrangian with respect to  1 ,  2 ,  3 , and  4 , respectively, are positive indicates that the optimal problem is a minimum at controls  * 1 ,  * 2 ,  * 3 , and  * 4 .The optimality system is solved using the forwardbackward fourth-order Runge-Kutta scheme in  statistical computing platform [37].The optimal strategy is obtained by solving the state and adjoint systems and the transversality conditions.First we start to solve state (1) using the Runge-Kutta fourth-order forward in time with a guess for the controls  1 ,  2 ,  3 , and  4 over the simulated time.Then, using the current iteration of the state equations with the initial guess for the controls, the adjoint equations in system (57) are solved by a backward method with the transversality conditions (62).Then the controls are updated by using a convex combination of the previous controls and the value from characterizations (63).This process is repeated and iterations stopped if the values of the unknowns at the previous iterations are very close to the ones at the present iterations [36].

Dynamics of Human Population with Intervention Strategies.
We simulated malaria model with intervention strategies to find the dynamics of human and mosquito populations.It is observed that the control strategy leads to decrease in the number of infected human ( ℎ ) as shown in Figure 2.
The figure shows steady decrease in susceptible human population at the initial period as the exposure of humans to disease increases.Thereafter graph of susceptible humans increases as the exposed and infected human population decreases due to positive effect of the intervention strategies being implemented.

Numerical Results on Optimal Control Analysis.
In this section, we investigate numerically the effect of the several optimal control strategies on the spread of malaria.We compare the numerical results from the simulation using one control and various combinations of two, three, and four control strategies.This was done by comparing when there were not any intervention strategies and when there were intervention strategies.There are 15 different control strategies for each of the four different epidemiological zones in Kenya that are explored.We use the case of endemic zone with the case of one control variable, two control variables, three control variables, and all the four control variables being in use for the illustration purpose.
The results in Figures 3-6 show a significant difference in  ℎ and   with the control strategy compared to  ℎ and   without the control strategy.It is observed that the control strategy leads to decrease in the number of infected humans ( ℎ ).The uncontrollable case leads to a decrease in the number of infected mosquitoes (  ), while the control strategy leads to decrease in the infected number.The control profiles show the upper bound time for each strategy for each setting before dropping to the lower bound.
There are several combinations for the different settings as described below.
Results of only one intervention strategy for the endemic epidemiological zones is as follows.
(a) Optimal Protection Using ITN.Only the control ( 1 ) on ITNs is used to optimize the objective function , while the control on treatment ( 2 ), the control on IRS ( 3 ), and control on IPTp ( 4 ) are set to zero.
(b) Optimal Treatment.Only the control ( 2 ) on treatment is used to optimize the objective function , while the control on ITNs ( 1 ), the control on IRS ( 3 ), and control on IPTp ( 4 ) are set to zero.
(c) Optimal IRS.Only the control ( 3 ) on IRS is used to optimize the objective function , while the control on treatment ( 2 ), the control on ITNs ( 1 ), and control on IPTp ( 4 ) are set to zero.
(d) Optimal IPTp.Only the control ( 4 ) on IPTp is used to optimize the objective function , while the control on treatment ( 2 ), the control on IRS ( 3 ), and control on ITNs ( 1 ) are set to zero.
Results of combining 2 intervention strategies for the endemic epidemiological zones are as follows.Results of combining the four intervention strategies for the endemic epidemiological zones are as follows.
Optimal ITN, Treatment, IRS, and IPTp.In this case all the control function ITNs control ( 1 ), treatment control ( 2 ), IRS control ( 3 ), and IPTp control ( 4 ) are used to optimize the objective function .
The same procedure is repeated for other combination for strategies and for other epidemiological zones (epidemic, seasonal, and low).Based on the simulation findings for the highest number of infections being inverted at a lower cost, the combined use of treatment and IRS reduces the infected human and mosquito population faster at a lower cost for the endemic settings (105 infections at $368.258).For the epidemic prone settings the use of treatment and IRS (111.03 infections at $388.6051) has more impact in reducing the infected human and mosquito population.For seasonal areas much impact will be felt when treatment is used (115.6983infections at $231.3967).For the low risk areas, just the use of ITNs and treatment (119.0659infections at 595.32) will be sufficient to reduce infected human and mosquito population.This is deduced from the intervention which takes shorter time to start reducing the number of infected mosquitoes and humans.

Discussion
In this paper, we formulated a mathematical model for the transmission dynamics of malaria with four time-dependent control measures in Kenya.We first consider control parameters to be constant and perform mathematical analysis of the model.The analysis showed that there exists a domain where the model is epidemiologically and mathematically well posed.Stability analysis of the model showed that the disease-free equilibrium is globally asymptotically stable if  0 ≤ 1 in .If  0 > 1, the unique endemic equilibrium exists and is globally asymptotically stable in .The model exhibited backward bifurcation backward bifurcation at  0 = 1 implying the existence of multiple endemic equilibria for  0 < 1.The existence of multiple endemic equilibria emphasizes the fact that  0 < 1 is not sufficient to eradicate disease from the population and the need to lower  0 much below one to make the disease-free equilibrium stable globally.This behavior has important public health implications because it means that bringing  0 below 1 is not enough to eradicate malaria.
We then consider the case of time-dependent control variable from where we formulated an optimal control problem and derived expressions for the optimal control for the malaria model with four control variables with an aim of minimizing the number of malaria infections in humans (derive optimal prevention and treatment strategies) while keeping the cost low.In the optimal control problem, use of one control at a time and the different combination of interventions can be explored to investigate and compare the effects of the control strategies on malaria eradication for different transmission settings.The analysis of the model showed that the state and optimal control can be calculated using the optimality system.The optimality system is comprised of the state system, the adjoint/costate system, initial conditions at  = 0, boundary conditions (transversality), and the characterization of the optimal control [11,12,14].
The results of the optimal control problem will be able to show which intervention or combination of the different intervention strategies has the highest impact on the control of the disease especially for different transmission settings.Our results shows that an effective IRS use and treatment will be beneficial to the community for the control of malaria disease (infected human and mosquito population) faster at a lower cost for the endemic settings.This is slightly different from the findings of Agusto et al. [14] who found that the combination of the personal protection, treatment, and insecticides spray had the highest impact on the control of the disease.This could be in endemic settings where both preventive and treatment measures work better which implies that the effect of protection using IRS is better.Griffin et al. [40] found that use of treatment, long-lasting insecticides treated bed nets (LLITNs), and IRS with high levels coverage would result in reducing malaria transmission for high settings though the study did not consider the cost aspect.
The findings show that, for the epidemic prone areas, the optimal control strategy for reducing the infected human and mosquito population was the use of treatment and IRS.This is slightly different from Agusto et al. [14] findings on resource limited settings in which the study recommended the use of personal protection and insecticides.This was further different from the findings of Mwamtobe et al. [12] who noted that the prevention strategies (use of ITNs and IRS) lead to the reduction of both the mosquito population and infected human individuals.This is because in epidemic areas emphasis is usually more placed on preventive strategies.
The results show that for seasonal areas much impact will be felt when treatment is used which is different from Mwamtobe et al. [12] who recommended IRS and ITNs.This is also comparable to Kim et al. [13] findings that mosquitoreduction strategies are more effective than personal protection.This is because in seasonal areas malaria transmission is usually not so high and hence if the mosquito-reduction strategies can be implemented then malaria transmission can be reduced.Griffin et al. [40] found that for the high seasonal transmission settings the use of LLITNs, IRS, and treatment would help reduce the transmission of malaria.
The results show that, for the low risk areas, just the use of ITNs and treatment will be sufficient to reduce infected human and mosquito population.This is comparable to Silva and Torres [15] who found the optimal use of ITNS would prevent malaria transmission the same as Kim et al. [13].The findings are comparable to those by Griffin et al. [40].In low transmission areas prevention strategies seem to be better because the population is not infected.
These findings support the WHO concerns on the capability of only one intervention strategy in reducing malaria transmission.The findings are however applicable to the designing of intervention strategies for malaria especially when costs aspects are of concern.This modelling approach also addresses effectiveness of the recommended intervention for at-risk group of malaria (pregnant women) by the WHO.The modelling approach has also been implemented in the  statistical computing platform which is free statistical software.
To the best of our knowledge, this is the first ever optimal control modelling and simulation of malaria intervention strategies in free  statistical computing platform; future testing and refinement of the model together with simulation with data from other representative settings should be done to improve the results and the model.These findings were based on the use of secondary data; a more designed study may be needed to ascertain the findings of these studies.Regardless, it does not invalidate the findings.

Conclusion
In this paper, we derived and analyzed a deterministic model for the transmission of malaria disease which incorporated the use of insecticide-treated bed nets (ITNs), treatment, indoor residual spray (IRS), and intermittent preventive treatment for pregnant women (IPTp) and performed optimal control analysis of the model.We first consider constant control parameters from where we investigate existence and stability of equilibria as well as stability analysis.We proved that if  0 ≤ 1, the disease-free equilibrium is globally asymptotically stable in .If  0 > 1, the unique endemic equilibrium exists and is globally asymptotically stable.The model also exhibits backward bifurcation at  0 = 1.
We then consider the time-dependent control case from where we derived and analyzed the necessary conditions for the optimal control of the disease.There were 15 different control strategies for each of the four different epidemiological zones in Kenya that were explored using control plots (numerical simulations) which compared when there were not any intervention strategies and when there were intervention strategies.Using the optimal control approach we found that the combined use of treatment and IRS would reduce the highest number of infected human and mosquito population faster at a lower cost for the endemic settings.For the epidemic prone settings the use of treatment and IRS has more impact in reducing the infected human and mosquito population.For seasonal areas much impact will be felt when treatment is used.For the low risk areas, just the use of ITNs and treatment will be sufficient to reduce infected human and mosquito population.This was deduced from the intervention which takes shorter time to start reducing the number of infected mosquitoes and humans.
We conclude that, according to our model, the optimal control strategy for malaria control in endemic areas was the combined use of treatment and IRS; in epidemic prone areas it was the use of treatment and IRS; in seasonal areas it was the use of treatment; and in low risk areas was the use of ITNs and treatment.Control programs that follow these strategies can effectively reduce the spread of malaria disease in different malaria transmission settings in Kenya.

Figure 2 :
Figure 2: Simulations showing the dynamics of human population with intervention strategies for the endemic setting.

Figure 3 :
Figure 3: Simulations of the model showing the effect of treatment on the spread of malaria for the endemic setting.

Figure 4 :
Figure 4: Simulations of the model showing the effect of ITNs and treatment on the spread of malaria for the endemic setting.

Figure 5 :
Figure 5: Simulations of the model showing the effect of ITNs, treatment, and IRS on the spread of malaria for the endemic setting.

Figure 6 :
Figure 6: Simulations of the model showing the effect of ITNs, treatment, IRS, and IPTp on the spread of malaria for the endemic setting.

Table 1 :
State variables of the malaria model.

Table 2 :
Description of parameter variables of the malaria model.

Table 4 :
Sensitivity indices (SI) of  0 to parameters for the malaria model.

Table 5 :
Parameter values for the full malaria model.

)
[11,12]objective function  is the final time and quantities  1 ,  2 , and  3 are weights constants of the total mosquito population, infected individuals, and exposed individuals, respectively, while  1 ,  2 ,  3 , and  4 are weight constants for use with ITNs, treatment effort, IRS, and IPTp efforts, respectively.The total mosquito population (  =   +   +   ) is part of the objective function because it is affected by the use of IRS and ITNs.In addition,  ℎ and  ℎ are included in the objective function because individuals in these classes are affected by the use of ITNs, IPTps, and treatment, respectively.A quadratic cost on the controls was chosen in line with what is known in the literature on epidemic optimal controls for malaria[11,12].The cost of implementing personal protection using ITNs is  1  2 1 , treatment of infected individuals is  2  2 2 , spraying of houses with IRS is  3  2 3 , and preventive method of IPTp is  4  2 4 .A linear function has been chosen for the cost incurred by exposed individuals  3  ℎ , infected individuals,  2  ℎ , and the mosquito population,  1   .A quadratic form is used for the cost on the controls  1  2 1 ,  2  2 2 ,  3  2 3 , and  4  2 4 , such that the terms (1/2) 1  2 1 , (1/2) 2  2 2 , (1/2) 3  2 3 , and (1/2) 4  2 4 describe the cost associated with the ITNs, treatment, mosquito control (IRS), and chemoprevention (IPTp), respectively.
H 3 ): right hand side of each equation of the state system in (1) is continuous, bounded by a linear functional in the state and control, and can be written as a linear function of  with coefficients depending on time and state.
1 ): the set of controls and corresponding state variables is nonempty.(H2 ): the control set  is convex and closed.( 1) is satisfied.The solutions are bounded and hence the control set satisfies ( 2 ) as well.The system is bilinear in  1 ,  2 ,  3 ,  4 and hence the right hand side of state system (1) satisfies condition ( 3 ) (since the solutions are bounded).Note that the integrand of the objective functional is convex.The last condition is also satisfied.The state and the control variables of system (1) are nonnegative values and nonempty.The control set  is closed and convex.The integrand of the objective cost function  expressed by (51) is a convex function of ( 1 ,  2 ,  3 ,  4 ) on the control set .The Lipschitz property of the state system with respect to the state variables is satisfied since the state solutions are bounded.It can easily be shown that there exist positive numbers  1 ,  2 and a constant  > 1 such that  1   +  2  ℎ +  3  ℎ where  1 ,  2 > 0,  1 ,  2 ,  1 ,  2 ,  3 ,  4 > 0, and  > 1.