\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2012 (2012), No. 81, pp. 1--22.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2012 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2012/81\hfil Application of optimal control] {Application of optimal control to the epidemiology of malaria} \author[F. B. Agusto, N. Marcus, K. O. Okosun\hfil EJDE-2012/81\hfilneg] {Folashade B. Agusto, Nizar Marcus, Kazeem O. Okosun } % in alphabetical order \address{Folashade B. Agusto \newline Department of Mathematics and Statistics, Austin Peay State University, Clarksville, Tennessee 37044, USA} \email{fbagusto@gmail.com} \address{Nizar Marcus \newline Department of Mathematics and Applied Mathematics, University of the Western Cape, South Africa} \email{nmarcus@uwc.ac.za} \address{Kazeem O. Okosun \newline Department of Mathematics, Vaal University of Technology, Private Bag X021, Vanderbijlpark, 1900, South Africa} \email{kazeemoare@gmail.com} \thanks{Submitted July 2, 2010. Published May 22, 2012.} \subjclass[2000]{92B05, 93A30, 93C15} \keywords{Malaria; optimal control; insecticide treated bed nets; \hfill\break\indent mosquito insecticide} \begin{abstract} Malaria is a deadly disease transmitted to humans through the bite of infected female mosquitoes. In this paper a deterministic system of differential equations is presented and studied for the transmission of malaria. Then optimal control theory is applied to investigate optimal strategies for controlling the spread of malaria disease using treatment, insecticide treated bed nets and spray of mosquito insecticide as the system control variables. The possible impact of using combinations of the three controls either one at a time or two at a time on the spread of the disease is also examined. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction} Malaria is a common and serious disease. It is reported that the incidence of malaria in the world may be in the order of 300 million clinical cases each year. Malaria mortality is estimated at almost 2 million deaths worldwide per year. The vast numbers of malaria deaths occur among young children in Africa, especially in remote rural areas. In addition, an estimated over 2 billion people are at risk of infection, no vaccines are available for the disease \cite {MAL, WHO}. Malaria is transmitted to humans through the bite of an infected female Anopheles mosquito, following the successful sporozoite inoculation, \emph{plasmodium falciparum} is usually first detected 7-11 days. This is followed after few days of the bites, by clinical symptoms such as sweats, shills, pains, and fever. Mosquitoes on the other hand acquire infection from infected human after a blood meal. Although malaria is life-threatening it is still preventable and curable if the infected individual seek treatment early. Prevention is usually by the use of insecticide treated bed nets and spraying of insecticide but according to the World Health Organization position statement on insecticide treated mosquito nets \cite{WHO3}, the insecticide treated bed nets(ITNs), long-lasting insecticide nets (LLINs), indoor residual spraying (IRS), and the other main method of malaria vector control, may not be sufficiently effective alone to achieve and maintain interruption of transmission of malaria, particularly in holo-endemic areas of Africa. Many studies have been carried out to quantify the impact of malaria infection in humans \cite{BRW,Green,KWM,MAC,ROSS,DLS}. Many of these studies focuses only on the transmission of the disease in human and the vector populations but recently, Chiyaka et.al \cite{chi} formulated a deterministic system of differential equations with two latent periods in the non-constant host and vector populations in order to theoretically assess the potential impact of personal protection, treatment and possible vaccination strategies on the transmission dynamics of malaria. Blayneh et al \cite{KYH}, used a time dependent model to study the effects of prevention and treatment on malaria, similarly Okosun \cite{oko} used a time dependent model to study the impact of a possible vaccination with treatment strategies in controlling the spread of malaria in a model that includes treatment and vaccination with waning immunity. Thus, following the WHO position statement \cite{WHO3} it is instructive to carry out modeling studies to determine the impact of various combinations of control strategies on the transmission dynamics of malaria. In this paper, we use treatment of symptomatic individuals, personal protection and the straying of insecticide as control measures and then consider this time dependent control measures using optimal control theory. Time dependent control strategies have been applied for the studies of HIV/AIDS disease, Tuberculosis, Influenza and SARS \cite{ABM, Agu, CEC, KARRA,Ki, SUPS,van, XIE}. Optimal control theory has been applied to models with vector-borne diseases \cite{KYH,rafi,Tho, WICK}. Our goal is to develop mathematical model for human-vector interactions with control strategies, with the aim of investigating the role of personal protection, treatment and spraying of insecticides in malaria transmission, in line with concerns raised WHO \cite{WHO3}; in order to determine optimal control strategies with various combinations of the control measures for controlling the spread of malaria transmission. The paper is organized as follows: in Section \ref{form}, we give the description of the human-vector model, stating the assumptions and definitions of the various parameters of the model. The analysis of the equilibrium points are discussed in Sections \ref{dfe} and \ref{stab}. In Section \ref{cont}, we state the control problem as well as the objective functional to be minimized, we then apply the Pontryagin's Maximum Principle to find the necessary conditions for the optimal control. In Sections \ref{num}, we show the simulation results to illustrate the population dynamics with preventative measures and treatment as controls. \section{Model formulation} \label{form} The model sub-divides the total human population at time t, denoted by $N_h(t)$, into the following sub-populations of susceptible individuals $(S_h(t))$, those exposed to malaria parasite $(E_h(t))$, individuals with malaria symptoms $(I_h(t))$, partially immune human $(R_h(t))$. So that $$ N_h(t) = S_h(t) + E_h(t) + I_h(t)+ R_h(t). $$ The total vector (mosquito) population at time $t$, denoted by $N_v(t)$, is sub-divided into susceptible mosquitoes $(S_v(t))$, mosquitoes exposed to the malaria parasite $(E_v(t))$ and infectious mosquitoes$(I_v(t))$. Thus, $$ N_v(t) = S_v(t) + E_v(t) + I_v(t). $$ It is assumed that susceptible humans are recruited into the population at a constant rate $\Lambda_h$. Susceptible individuals acquire malaria infection following contact with infectious mosquitoes (at a rate $\beta\varepsilon_h\phi$), where $\beta$ is the transmission probability per bite and $\varepsilon_h$ is the biting rate of mosquitoes, $\phi$ is contact rate of vector per human per unit time. Susceptible individuals infected with malaria are moved to the exposed class $(E_h)$ at the rate $\beta\varepsilon_h\phi$ and then progress to the infectious class, following the development of clinical symptoms (at a rate $\alpha_h$). Individuals with malaria symptoms are effectively treated (at a rate $\tau$) where ($0\leq\tau\leq 1$). Human spontaneous recovery rate is given by $b$, where $0\leq b<\tau$. And individuals infected with malaria suffer a disease-induced death (at a rate $\psi$). Infected individual then progress to the partially immuned group. Upon recovery, the partially immuned individual losses immunity (at the rate $\kappa$) and becomes susceptible again. Susceptible mosquitoes ($S_v$) are generated at the rate $\Lambda_v$ and acquire malaria infection (following effective contacts with humans infected with malaria) at a rate $\lambda\phi\varepsilon_v (I_h + \eta R_h)$, where $\lambda$ is the probability of a vector getting infected through the infectious human and $\varepsilon_v$ is the biting rate of mosquitoes. We assume that humans in the $R_h(t)$ class can still transmit the disease, thus, the modification parameter $\eta \in [0,1)$ gives the reduced infectivity of the recovered individuals \cite{Cov,Ros}. Mosquitoes are assumed to suffer natural death at a rate $\mu_v$, regardless of their infection status. Newly-infected mosquitoes are moved into the exposed class ($E_v$ ), and progress to the class of symptomatic mosquitoes ($I_v$) following the development of symptoms (at a rate $\alpha_v$). Thus, putting the above formulations and assumptions together gives the following human-vector model, given by system of ordinary differential equations below as \begin{equation} \begin{gathered} \frac{dS_h}{dt} = \Lambda_h + \kappa R_h - \beta\varepsilon_h\phi I_vS_h - \mu_h S_h, \\ \frac{dE_h}{dt} = \beta\varepsilon_h\phi I_vS_h - (\alpha_h+\mu_h )E_h,\\ \frac{dI_h}{dt} = \alpha_hE_h -(b + \tau) I_h - ( \psi + \mu_h)I_h ,\\ \frac{dR_h}{dt} = (b + \tau) I_h - (\kappa+ \mu_h) R_h,\\ \frac{dS_v}{dt} = \Lambda_v - \lambda\phi\varepsilon_v(I_h + \eta R_h) S_v - \mu_v S_v,\\ \frac{dE_v}{dt} = \lambda\phi\varepsilon_v(I_h + \eta R_h)S_v - (\alpha_v +\mu_v) E_v,\\ \frac{dI_v}{dt} = \alpha_vE_v - \mu_v I_v, \end{gathered} \label{eM13} \end{equation} The associated model variables and parameters are described in Table \ref{p1}. \subsection{Basic properties of the malaria model} \subsubsection{Positivity and boundedness of solutions}\label{Pos} For the malaria transmission model \eqref{eM13} to be epidemiologically meaningful, it is important to prove that all its state variables are non-negative for all time. In other words, solutions of the model system \eqref{eM13} with non-negative initial data will remain non-negative for all time $t > 0$. \begin{theorem} \label{thm0} Let the initial data $S_h(0) \geq 0$, $E_h(0)\geq 0$, $I_h(0) \geq 0$, $R_h(0) \geq 0$, $S_v(0) \geq 0$, $E_v(0)\geq 0$, $I_v(0) \geq 0$. Then the solutions $(S_h, E_h, I_h, R_h$, $S_v, E_v, I_v)$ of the malaria model \eqref{eM13} are non-negative for all $t > 0$. Furthermore $$ \limsup_{t\to\infty} N_h(t)\leq\frac{\Lambda_h}{\mu_h},\quad \limsup_{t\to\infty} N_v(t)\leq\frac{\Lambda_v}{\mu_v}, $$ with $N_h = S_h+ E_h+ I_h+R_h$ and $N_v = S_v+E_v+I_v$. \end{theorem} \begin{proof} Let $t_1=\sup \{t>0:S_h(t)>0,E_h(t)>0,I_h(t)>0,R_h(t)>0, S_v(t)>0,I_v(t)>0,E_v(t)>0\}$. Since $S_h(0)>0,E_h(0)>0,I_h(0)>0,R_h(0)>0,S_v(0)>0,E_v(0)>0,I_v(0)>0 $, then, $t_1>0$. If $t_1<\infty $ , then $S_h$, $E_h$, $I_h$, $R_h$, $S_v$, $E_v$ or $I_v$ is equal to zero at $t_1$. It follows from the first equation of the system \eqref{eM13}, that $$ \frac{dS_h}{dt} = \Lambda_h - \beta\varepsilon_h\phi I_vS_h - \mu_h S_h + \kappa R_h $$ Thus, $$ \frac{d}{dt}\big\{S_h(t) \exp[(\beta\varepsilon_h\phi I_v + \mu_h )t] \big\} = (\Lambda_h+ \kappa R_h)\exp [ (\beta\varepsilon_h\phi I_v + \mu_h )t] $$ Hence, $$S_h(t_1)\exp[(\beta\varepsilon_h\phi I_v + \mu_h )t ] -S_h(0)= \int^{t_1}_0(\Lambda_h+ \kappa R_h) \exp [ (\beta\varepsilon_h\phi I_v + \mu_h )p]dp $$ so that \begin{align*} S_h(t_1) &= S_h(0)~ exp[-(\beta\varepsilon_h\phi I_v + \mu_h )t_1 ] +\exp[-(\beta\varepsilon_h\phi I_v + \mu_h )t_1]\\ &\quad\times \int^{t_1}_0(\Lambda_h+ \kappa R_h)\exp [ (\beta\varepsilon_h\phi I_v + \mu_h )p ]dp>0. \end{align*} and \begin{align*} R_h(t_1)&=R_h(0)\exp[-(\mu_h + \kappa)t_1] + \exp[(\mu_h + \kappa)t_1]\int^{t_1}_0 (b+\tau) I_h \exp[(\mu_h + \kappa)p ]dp\\ &>0. \end{align*} It can similarly be shown that $ E_h > 0$, $I_h > 0$, $S_v > 0$, $E_v > 0$ and $I_v > 0$ for all $t>0$. For the second part of the proof, it should be noted that $00$). The rate of change of the humans and mosquitoes populations is given in equation \eqref{eDW}, it follows that \begin{equation} \label{e3} \begin{gathered} \frac{dN_h(t)}{dt}\leq \Lambda_h -\mu_hN_h(t),\\ \frac{dN_v(t)}{dt}\leq \Lambda_v -\mu_vN_v(t). \end{gathered} \end{equation} A standard comparison theorem \cite{Lak} can then be used to show that $N_h(t)\leq N_h(0)e^{-\mu_ht}+\frac{\Lambda_h}{\mu_h}(1-e^{-\mu_ht})$ and $N_v(t)\leq N_v(0)e^{-\mu_vt}+\frac{\Lambda_v}{\mu_v}(1-e^{-\mu_vt})$. In particular, $N_h(t)\leq\frac{\Lambda_h}{\mu_h}$ and $N_v(t)\leq\frac{\Lambda_v}{\mu_v}$ if $N_h(0)\leq\frac{\Lambda_h}{\mu_h}$ and $N_v(0)\leq\frac{\Lambda_v}{\mu_v}$ respectively. Thus, the region $\mathcal{D}$ is positively-invariant. Hence, it is sufficient to consider the dynamics of the flow generated by \eqref{eM13} in $\mathcal{D}$. In this region, the model can be considered as been epidemiologically and mathematically well-posed \cite{Het}. Thus, every solution of the basic model \eqref{eM13} with initial conditions in $\mathcal{D}$ remains in $\mathcal{D}$ for all $t > 0$. Therefore, the $\omega$-limit sets of the system \eqref{eM13} are contained in $\mathcal{D}$. This result is summarized below. \begin{lemma} \label{lem1} The region $\mathcal{D}= \mathcal{D}_h\cup \mathcal{D}_v \subset \mathbb{R}^4_+\times \mathbb{R}^3_+$ is positively-invariant for the basic model \eqref{eM13} with non-negative initial conditions in $\mathbb{R}^7_+$ \end{lemma} \subsection{Stability of the disease-free equilibrium (DFE)} \label{dfe} The malaria model \eqref{eM13} has a DFE, obtained by setting the right-hand sides of the equations in the model to zero, given by $$ \mathcal{E}_0 = (S^*_h,E^*_h,I^*_h,R^*_h,S^*_v,E^*_v,I^*_v) = \Big(\frac{\Lambda_h}{ \mu_h},0,0,0,\frac{\Lambda_v}{\mu_v},0,0\Big). $$ The linear stability of $\mathcal{E}_0$ can be established using the next generation operator method \cite{van} on the system \eqref{eM13}, the matrices $F$ and $V,$ for the new infection terms and the remaining transfer terms, are, respectively, given by \begin{gather*} F=\begin{pmatrix} 0&0 &0 & 0 &\beta\varepsilon_h\phi S^*_h\\ 0 & 0 & 0 &0 & 0\\ 0 & 0 & 0 &0 & 0\\ 0&\lambda\varepsilon_v\phi S^*_v &\lambda\varepsilon_v\phi\eta S^*_v &0 &0 \\ 0 & 0 & 0 &0 & 0 \end{pmatrix}, \\ V=\begin{pmatrix} k_1& 0 & 0 & 0 &0\\ -\alpha_1 & k_2 & 0 &0 & 0\\ 0& -(b+\tau) & k_3 &0 & 0 \\ 0&0 & 0 & k_4 & 0 \\ 0 &0 & 0 & -\alpha_2 & \mu_v \end{pmatrix}, \end{gather*} where $k_1= \alpha_h+\mu_h$, $k_2=b+\tau+\psi+\mu_h$, $k_3= \kappa+\mu_h$, $k_4= \alpha_v+\mu_v$. It follows that the reproduction number of the malaria system \eqref{eM13}, denoted by $\mathcal{R}_0$, is \begin{equation} \mathcal{R}_0 = \sqrt{\frac{\alpha_1\alpha_2\lambda \beta[k_3 +\eta(b+\tau)]\phi^2\epsilon_h\varepsilon_vS^*_hS^*_v}{k_3 k_4 k_2 k_1\mu_v}}, \label{R1} \end{equation} Further, using \cite[Theorem 2]{van}, the following result is established. \begin{theorem} \label{thm1} The DFE of the model \eqref{eM13}, given by $\mathcal{R}_{0}$, is locally asymptotically stable (LAS) if $\mathcal{R}_0 < 1$, and unstable if $\mathcal{R}_0 > 1$. \end{theorem} \section{Existence of endemic equilibrium point (EEP)} \label{stab} Next conditions for the existence of endemic equilibria for the model \eqref{eM13} is explored. Let $$ \mathcal{E}_1 = \big(S^{**}_h,E^{**}_h,I^{**}_h,R^{**}_h,S^{**}_v,E^{**}_v, I^{**}_v\big), $$ be the arbitrary endemic equilibrium of model \eqref{eM13}, in which at least one of the infected components of the model is non-zero. Let \begin{gather} \lambda^{**}_h = \beta\phi\varepsilon_hI_v , \label{eLH} \\ \lambda^{**}_v = \lambda\phi\varepsilon_v (I_h + \eta R_h) \label{eLV} \end{gather} be the force of infection in humans and in the vector. Setting the right-hand sides of the equations in \eqref{eM13} to zero gives the following expressions (in terms of $\lambda^{**}_h$ and $\lambda^{**}_v$) \begin{equation} \label{eEEP} \begin{gathered} S^{**}_h = \frac{\Lambda^{**}_h k_1 k_2 k_3}{(\lambda_h +\mu_h ) k_1 k_2 k_3 -\kappa\lambda^{**}_h \alpha_h (b+\tau)},\\ E^{**}_h = \frac{k_2\lambda^{**}_h\Lambda_h k_3}{(\lambda_h +\mu_h ) k_1 k_2 k_3 -\kappa\lambda^{**}_h \alpha_h (b+\tau)}, \\ I^{**}_h = \frac{\lambda^{**}_h\Lambda_h k_3\alpha_1}{(\lambda_h +\mu_h ) k_1 k_2 k_3 -\kappa\lambda^{**}_h \alpha_h (b+\tau)},\\ R^{**}_h = \frac{(b+\tau)\lambda^{**}_h\Lambda_h\alpha_1}{(\lambda_h +\mu_h ) k_1 k_2 k_3 -\kappa\lambda^{**}_h \alpha_h (b+\tau)}, \\ S^{**}_v = \frac{\Lambda_v }{(\lambda^{**}_v +\mu_v )}, \quad E^{**}_v = \frac{\lambda^{**}_v\Lambda_v}{k_4(\lambda^{**}_v +\mu_v ) } ,\quad I^{**}_v = \frac{\alpha_v\lambda^{**}_v\Lambda_v}{k_4\mu_v(\lambda^{**}_v +\mu_v )} \end{gathered} \end{equation} Substituting \eqref{eEEP} and \eqref{eLV} into \eqref{eLH}, gives $a_0 \lambda^{**}_h + b_0=0$, where \begin{gather*} a_0 = k_4\mu_v \{\lambda\phi\varepsilon_v\Lambda_h\alpha_h [k_3 +\eta (b+\tau)]+\mu_v [k_3 k_1 k_2 -\kappa\alpha_h (b+\tau)]\} \\ b_0 = \mu_h\mu_v^2 k_4 k_3 k_1 k_2(1-\mathcal{R}^2_0). \end{gather*} The coefficient $a_0$ is always positive, the coefficient $b_0$ is positive (negative) if $\mathcal{R}_0$ is less than (greater than) unity. Furthermore, there is no positive endemic equilibrium if $b_0\geq 0$. If $b_0 < 0$, then there is a unique endemic equilibrium (given by $\lambda_h = b_0/a_0$). This result is summarized below. \begin{lemma} \label{lem2} The model \eqref{eM13} has a unique positive endemic equilibrium whenever $\mathcal{R}_0 > 1$, and no positive endemic equilibrium otherwise. \end{lemma} \subsection{Global stability of endemic equilibrium for a special case} In this section, we investigate the global stability of the endemic equilibrium of model \eqref{eM13}, for the special case when $\kappa=0$, that there is no lost of immunity. Using the approach in the proof of Lemma \ref{lem1}, it can be shown that the region $$ \tilde{\mathcal{D}}= \tilde{\mathcal{D}}_h \cup \tilde{\mathcal{D}}_v\subset \mathbb{R}^4_+\times \mathbb{R}^3_+, $$ where \begin{gather*} \tilde{\mathcal{D}}_h = \big\{ (S_h,E_h,I_h,R_h) \subset \mathcal{D}_h~\colon S_h \leq S^*_h\big\}, \\ \tilde{\mathcal{D}}_v = \big\{ (S_v,E_v,I_v) \subset \mathcal{D}_v : S_v\leq S^*_v\big\}. \end{gather*} is positively-invariant for the special case of the model \eqref{eM13} described above. It is convenient to define $$ \tilde{\mathcal{D}} =\big\{ (S_h,E_h,I_h,R_h,S_v,E_v,I_v) \in \mathcal{D}~\colon E_h=I_h=R_h=E_v=I_v=0\big\}. $$ \begin{theorem} \label{impeep} The unique endemic equilibrium, $\tilde{\mathcal{E}}_1$, of the model \eqref{eM13} is GAS in $\tilde{\mathcal{D}}\backslash \tilde{\mathcal{D}}_0$ whenever $\tilde{\mathcal{R}}_0{_{|_{\kappa=0}}} > 1$. \end{theorem} \begin{proof} Let $\tilde{\mathcal{R}}_0 > 1$, so that the unique endemic equilibrium ($\tilde{\mathcal{E}}_1$) exists. Consider the non-linear Lyapunov function \begin{align*} \mathcal{F} &= S_h^{**}\Big(\frac{S_h}{S_h^{**}} - \ln \frac{S_h}{S_h^{**}} \Big) + E_h^{**}\Big(\frac{E_h}{E_h^{**}} - \ln \frac{E_h}{E_h^{**}} \Big) + \frac{k_1}{\alpha_h}I_h^{**}\Big(\frac{I_h}{I_h^{**}} - \ln \frac{I_h}{I_h^{**}}\Big) \\ &\quad + \frac{k_2k_1}{\alpha_h\gamma}R_h^{**}\Big(\frac{R_h}{R_h^{**}} - \ln \frac{R_h}{R_h^{**}}\Big) + S_v^{**}\Big(\frac{S_v}{S_v^{**}} - \ln \frac{S_v}{S_v^{**}} \Big) + E_v^{**}\Big(\frac{E_v}{E_v^{**}} - \ln \frac{E_v}{E_v^{**}} \Big) \\ &\quad + \frac{k_4}{\alpha_v}I_v^{**}\Big(\frac{I_v}{I_v^{**}} - \ln \frac{I_v}{I_v^{**}}\Big) , \end{align*} where $\gamma = b+\tau$ and the Lyapunov derivative is \begin{align*} \dot{\mathcal{F}} &= \Big(1 - \frac{S_h^{**}}{S_h} \Big)\dot{S_h} + \Big(1 - \frac{E_h^{**}}{E_h} \Big)\dot{E}_h + \frac{k_1}{\alpha_h}\Big(1 - \frac{I_h^{**}}{I_h}\Big)\dot{I_h} + \frac{k_2k_1}{\alpha_h\gamma}\Big(1 - \frac{R_h^{**}}{R_h}\Big)\dot{R_h}\\ &\quad + \Big(1 - \frac{S_v^{**}}{S_v} \Big)\dot{S_v} + \Big(1 - \frac{E_v^{**}}{E_v} \Big)\dot{E}_v + \frac{k_4}{\alpha_v}\Big(1 - \frac{I_v^{**}}{I_v}\Big)\dot{I_v}. \end{align*} Substituting the expressions for the derivatives in $\dot{\mathcal{F}}$ (from \eqref{eM13} with $\kappa=0$) gives \begin{align*} \dot{\mathcal{F}} &= \Lambda_h - \lambda_hS_h - \mu_h S_h - \frac{S_h^{**}}{S_h} \Big(\Lambda_h - \lambda_hS_h - \mu_h S_h \Big)\\ &\quad + \lambda_hS_h - k_1E_h - \frac{E_h^{**}}{E_h} \Big( \lambda_hS_h - k_1E_h \Big)\\ &\quad + \frac{k_1}{\alpha_h} \Big(\alpha_hE_h -k_2I_h\Big) - \frac{k_1}{\alpha_h} \frac{I_h^{**}}{I_h}\Big(\alpha_hE_h -k_2I_h \Big)\\ &\quad + \frac{k_2k_1}{\alpha_h\gamma}\Big( \gamma I_h - k_3 R_h\Big) - \frac{k_2k_1}{\alpha_h\gamma}\frac{R_h^{**}}{R_h} \Big(\gamma I_h - k_3 R_h\Big)\\ &\quad + \Lambda_v - \lambda_v S_v - \mu_v S_v -\frac{S_v^{**}}{S_v} \Big( \Lambda_v - \lambda_v S_v - \mu_v S_v \Big)\\ &\quad + \lambda_v S_v -k_4 E_v + \frac{E_v^{**}}{E_v} \Big( \lambda_v S_v -k_4 E_v \Big) \\ &\quad + \frac{k_4}{\alpha_v}\Big( \alpha_vE_v - \mu_v I_v \Big) + \frac{k_4}{\alpha_v}\frac{I_v^{**}}{I_v}\Big(\alpha_vE_v - \mu_v I_v\Big). \end{align*} %\label{eFP} so that \begin{equation} \label{eFPF} \begin{split} \dot{\mathcal{F}} &= \lambda_h S_h^{**}\Big(1-\frac{S_h^{**}}{S_h}\Big) + \mu_hS^{**}_h\Big(2-\frac{S_h}{S^{**}_h}-\frac{S^{**}_h}{S_h}\Big) + \lambda_hS^{**}_h - \frac{E^{**}_h}{E_h}\lambda_h S_h \\ &\quad + k_1 E^{**}_h - k_1\frac{I^{**}_h}{I_h}E_h + \frac{k_2k_1}{\alpha_h}I^{**}_h -\frac{k_2k_1}{\alpha_h}\frac{R^{**}_h}{R_h}I_h + \frac{k_3k_2k_1}{\alpha_h\gamma}R_h -\frac{k_3k_2k_1}{\alpha_h\gamma}R_h \\ &\quad + \lambda_v S_v^{**}\Big(1-\frac{S_v^{**}}{S_v}\Big) + \mu_vS^{**}_v\Big(2-\frac{S_v}{S^{**}_v}-\frac{S^{**}_v}{S_v}\Big) + \lambda_vS^{**}_v\\ &\quad - \frac{E^{**}_v}{E_v}\lambda_v S_v + k_4 E^{**}_v - k_4\frac{I^{**}_v}{I_v}E_v + \frac{k_4\mu_v}{\alpha_v}I^{**}_v - \frac{k_4\mu_v}{\alpha_v}I_v \end{split} \end{equation} Finally, equation \eqref{eFPF} can be further simplified to give \begin{equation} \begin{split} {\dot{\mathcal{F}}} &= \mu_h S_h^{**} \Big(2- \frac{S_h^{**}}{S_h} -\frac{S_h}{S_h^{**}} \Big) +k_1E_h^{**} \Big( 5 - \frac{S_h^{**}}{S_h} - \frac{E_h^{**}}{E_h}- \frac{E_h}{E^{**}_h}\frac{I_h^{**}}{I_h}\\ &\quad - \frac{I_h}{I_h^{**}}\frac{R_h^{**}}{R_h} - \frac{R_h}{R_h^{**}} \Big) + \ \mu_v S_v^{**} \Big(2- \frac{S_v^{**}}{S_v} -\frac{S_v}{S_v^{**}} \Big) \\ &\quad +k_4E_v^{**} \Big( 4 - \frac{S_v^{**}}{S_v}- \frac{E_v^{**}}{E_v} - \frac{E_v}{E^{**}_v}\frac{I_v^{**}}{I_v}- \frac{I_v}{I_v^{**}} \Big). \end{split}\label{eF2} \end{equation} Since the arithmetic mean exceeds the geometric mean, it follows that \begin{gather*} 2 -\frac{S_h^{**}}{S_h} -\frac{S_h}{S_h^{**}} \leq 0,\quad 2 - \frac{S_v^{**}}{S_v} -\frac{S_v}{S_v^{**}}\leq 0, \\ 4 - \frac{S_v^{**}}{S_v}- \frac{E_v^{**}}{E_v} - \frac{E_v}{E^{**}_v}\frac{I_v^{**}}{I_v}- \frac{I_v}{I_v^{**}} \leq 0,\\ 5 - \frac{S_h^{**}}{S_h}- \frac{E_h^{**}}{E_h} - \frac{E_h}{E^{**}_h}\frac{I_h^{**}}{I_h} - \frac{I_h}{I_h^{**}}\frac{R_h^{**}}{R_h} - \frac{R_h}{R_h^{**}} \leq 0 \end{gather*} Since all the model parameters are non-negative, it follows that $\dot{\mathcal{F}}\leq 0 $ for $\tilde{\mathcal{R}}_{0}{_{|_{\kappa=0}}}>1$. Thus, it follows from the LaSalle's Invariance Principle, that every solution to the equations in the model \eqref{eM13} (with initial conditions in $\tilde{\mathcal{D}}\backslash \tilde{\mathcal{D}}_0$) approaches the EEP ($\tilde{\mathcal{E}}_1$) as $t\to \infty$ whenever $\tilde{\mathcal{R}}_{0}{_{|_{\kappa=0}}}>1$. \end{proof} \section{Analysis of optimal control} \label{cont} We introduce into the model \eqref{eM13}, time dependent preventive ($u_1,u_3$) and treatment ($u_2$) efforts as controls to curtail the spread of malaria. The malaria model \eqref{eM13} becomes \begin{gather} \frac{dS_h}{dt} = \Lambda_h + \kappa R_h - (1-u_1)\beta\varepsilon_h\phi I_vS_h - \mu_h S_h, \notag \\ \frac{dE_h}{dt} = (1-u_1)\beta\varepsilon_h\phi I_vS_h - (\alpha_h+\mu_h )E_h, \notag\\ \frac{dI_h}{dt} = \alpha_hE_h -(b + u_2) I_h - ( \psi + \mu_h)I_h , \notag\\ \frac{dR_h}{dt} = (b + u_2) I_h - (\kappa+ \mu_h) R_h, \notag\\ \frac{dS_v}{dt} = \Lambda_v - (1-u_1)\lambda\varepsilon_v\phi (I_h + \eta R_h) S_v - u_3(1-p)S_v - \mu_v S_v, \label{eM12}\\ \frac{dE_v}{dt} = (1-u_1)\lambda\varepsilon_v\phi (I_h + \eta R_h) S_v - u_3(1-p)E_v - (\alpha_v +\mu_v) E_v, \notag\\ \frac{dI_v}{dt} = \alpha_vE_v - u_3(1-p)I_v - \mu_v I_v. \notag \end{gather} The function $0\leq u_1\leq 1$ represent the control on the use of mosquitoes treated bed nets for personal protection, and $0 \leq u_2 \leq a_2$, the control on treatment, where $a_2$ is the drug efficacy use for treatment. The insecticides used for treating bed nets is lethal to the mosquitoes and other insects and also repels the mosquitoes, thus, reducing the number that attempt to feed on people in the sleeping areas with the nets \cite{CDC,WHO3}. However, the mosquitoes can still feed on humans outside this protective areas, and so we have included the spraying of insecticide. Thus, each mosquitoes group is reduced (at the rate $u_3$ $(1 - p)$), where $(1 - p)$ is the fraction of vector population reduced and $0 \leq u_3 \leq a_3$, is the control function representing spray of insecticide aimed at reducing the mosquitoes sub-populations and $a_3$ represent the insecticide efficacy at reducing the mosquitoes population. This is different from what was implemented in \cite{KYH}, where only two control measures of personal protection and treatment were used. With the given objective function \begin{equation} J(u_1, u_2, u_3) = \int^{t_f}_0[ mI_h + nu_1^2 + cu_2^2 + du_3^2 ]dt \label{eJ} \end{equation} where $t_f$ is the final time and the coefficients $m,n,c,d$ are positive weights to balance the factors. Our goal is to minimize the number of infected humans $I_h(t)$, while minimizing the cost of control $u_1(t),~u_2(t),~u_3(t)$. Thus, we seek an optimal control $u^*_1, u_2^*, u_3^*$ such that \begin{equation} J(u^*_1, u_2^*, u_3^*)=\min_{u_1,u_2,u_3}\{J(u_1,u_2,u_3)|u_1,u_2,u_3\in\mathcal{U}\} \end{equation} where the control set $$ \mathcal{U}=\{(u_1,u_2,u_3)\mid u_i:[0,t_f]\to [0,1], \text{ Lebesgue measurable } i=1,2,3\}. $$ The term $mI_h$ is the cost of infection while $nu^2_1$, $cu_2^2$ and $du^2_3$ are the costs of use of bed nets, treatment efforts and use of insecticides respectively. The necessary conditions that an optimal control must satisfy come from the Pontryagin's Maximum Principle \cite{Pon}. This principle converts \eqref{eM12}-\eqref{eJ} into a problem of minimizing pointwise a Hamiltonian $H$, with respect to $(u_1,u_2,u_3)$ \begin{equation} \begin{split} H&= m I_h+n u^2_1+cu_2^2 +du_3^2 + \lambda_{S_h}\{\Lambda_h + \kappa R_h - (1-u_1)\beta\varepsilon_h\phi I_vS_h - \mu_h S_h\}\\ &\quad + \lambda_{E_h}\{(1-u_1)\beta\varepsilon_h\phi I_vS_h - (\alpha_h+\mu_h )E_h \}\\ &\quad + \lambda_{I_h}\{\alpha_hE_h -(b + u_2) I_h - ( \psi + \mu_h)I_h \}\\ &\quad + \lambda_{R_h}\{(b + u_2) I_h - (\kappa+ \mu_h) R_h \}\\ &\quad + \lambda_{S_v}\{\Lambda_v - (1-u_1)\lambda\varepsilon_v\phi (I_h + \eta R_h) S_v - u_3(1-p)S_v - \mu_v S_v \}\\ &\quad + \lambda_{E_v}\{ (1-u_1)\lambda\varepsilon_v\phi (I_h + \eta R_h) S_v - u_3(1-p)E_v - (\alpha_v +\mu_v) E_v\}\\ &\quad + \lambda_{I_v}\{\alpha_v2E_v - u_3(1-p)I_v - \mu_v I_v \} \end{split}\label{eH} \end{equation} where the $\lambda_{S_h}$, $\lambda_{E_h}$, $\lambda_{I_h}$, $\lambda_{R_h}$, $\lambda_{S_v}$, $\lambda_{E_v}$, $\lambda_{I_v}$ are the adjoint variables or co-state variables. \cite[Corollary 4.1]{Fl} gives the existence of optimal control due to the convexity of the integrand of $J$ with respect to $u_1,~u_2$ and $u_{3}$, a \textit{priori} boundedness of the state solutions, and the \textit{Lipschitz} property of the state system with respect to the state variables. Applying Pontryagin's Maximum Principle \cite{Pon} and the existence result for the optimal control from \cite{Fl}, we obtain the following theorem. \begin{theorem} \label{thm4} Given an optimal control $u^*_1$, $u^*_2$, $u^*_3$ and solutions $S_h^*$, $E_h^*$, $I_h^*$, $R_h^*$, $S_v^*$, $E_v^*$, $I_v^*$ of the corresponding state system \eqref{eM12} that minimizes $J(u_1,u_2,u_3)$ over $\mathcal{U}$. Then there exists adjoint variables $\lambda_{S_h}$, $\lambda_{E_h}$, $\lambda_{I_h}$, $\lambda_{R_h}$, $\lambda_{S_v}$, $\lambda_{E_v}$, $\lambda_{I_v}$ satisfying \begin{equation} \begin{gathered} - \frac{d\lambda_{S_h}}{dt}= -[(1-u_1)\beta\varepsilon_h\phi I_v + \mu_h] \lambda_{S_h} + (1-u_1)\beta\varepsilon_h\phi I_v\lambda_{E_h} \\ - \frac{d\lambda_{E_h}}{dt}= -(\mu_h + \alpha_h)\lambda_{E_h} + \alpha_h\lambda_{I_h} \\ \begin{aligned} - \frac{d\lambda_{I_h}}{dt} &= m - [ (b+u_2)+ (\mu_h + \psi)]\lambda_{I_h} + (b+u_2)\lambda_{R_h} \\ &\quad + (1-u_1)\lambda\varepsilon_v\phi S_v(\lambda_{E_v} -\lambda_{S_v}) \end{aligned} \\ - \frac{d\lambda_{R_h}}{dt}= \kappa\lambda_{S_h} - (\mu_h + \kappa)\lambda_{R_h} + (1-u_1)\lambda\varepsilon_v\phi\eta S_v(\lambda_{S_v}-\lambda_{E_v}) \\ \begin{aligned} - \frac{d\lambda_{S_v}}{dt}&= - [(1-u_1)\lambda\varepsilon_v\phi(I_h+\eta R_h) + u_3(1-p) + \mu_v]\lambda_{S_v} \\ &\quad + (1-u_1)\lambda\varepsilon_v\phi(I_h+\eta R_h)\lambda_{E_v} \end{aligned}\\ - \frac{d\lambda_{E_v}}{dt}= -[ u_3(1-p)+ \alpha_v + \mu_v]\lambda_{E_v} + \alpha_v\lambda_{I_v} \\ - \frac{d\lambda_{I_v}}{dt}= -(1-u_1)\beta\varepsilon_h\phi S_h\lambda_{S_h} + (1-u_1)\beta\varepsilon_h\phi S_h\lambda_{E_h} - [u_3(1-p)+ \mu_v ]\lambda_{I_v} \end{gathered}\label{eA} \end{equation} and with transversality conditions \begin{equation} \lambda_{S_h}(t_f)=\lambda_{E_h}(t_f)=~\lambda_{I_h}(t_f) =\lambda_{R_h}(t_f)=\lambda_{S_v}(t_f)=\lambda_{E_v}(t_f)=\lambda_{I_v}(t_f)=0 \label{eT} \end{equation} and the controls $u_1^*, u_2^*$ and $u_3^*$ satisfy the optimality condition \begin{equation} \begin{gathered} u^*_1 = \max\Big\{0,\min\Big(1, \frac{ \beta\varepsilon_h\phi I^*_v(\lambda_{E_h} - \lambda_{S_h}) S^*_h+ \lambda\varepsilon_v\phi (I^*_h+\eta R^*_h)(\lambda_{E_v} - \lambda_{S_v}) S^*_v }{2n}\Big)\Big\}, \\ u^*_2 = \max\Big\{0,\min\Big(1, \frac{ (\lambda_{I_h} - \lambda_{R_h} )I_h^* }{2c}\Big)\Big\}\\ u^*_3 = \max\Big\{0,\min\Big(1, \frac{ (1-p)(S_v^*\lambda_{S_v} + E_v^*\lambda_{E_v } + I_v^*\lambda_{I_v} ) }{2d}\Big)\Big\} \end{gathered} \label{eC} \end{equation} \end{theorem} \begin{proof} The differential equations governing the adjoint variables are obtained by differentiation of the Hamiltonian function, evaluated at the optimal control. Then the adjoint system can be written as \begin{gather*} - \frac{d\lambda_{S_h}}{dt}= \frac{\partial H}{\partial S_h} = -[(1-u_1)\beta\varepsilon_h\phi I_v + \mu_h]\lambda_{S_h} + (1-u_1)\beta\varepsilon_h\phi I_v\lambda_{E_h} \\ - \frac{d\lambda_{E_h}}{dt}= \frac{\partial H}{\partial E_h} = -(\mu_h + \alpha_h)\lambda_{E_h} + \alpha_h\lambda_{I_h} \\ \begin{aligned} - \frac{d\lambda_{I_h}}{dt}= \frac{\partial H}{\partial I_h} &= m - [ (b+u_2)+ (\mu_h + \psi)]\lambda_{I_h} + (b+u_2)\lambda_{R_h} \\ &\quad + (1-u_1)\lambda\varepsilon_v\phi S_v(\lambda_{E_v}-\lambda_{S_v}) \end{aligned}\\ - \frac{d\lambda_{R_h}}{dt}= \frac{\partial H}{\partial R_h} = \kappa\lambda_{S_h} - (\mu_h + \kappa)\lambda_{R_h} + (1-u_1)\lambda\varepsilon_v\phi\eta S_v(\lambda_{S_v}-\lambda_{E_v}) \\ \begin{aligned} - \frac{d\lambda_{S_v}}{dt}= \frac{\partial H}{\partial S_v} & = - [(1-u_1)\lambda\varepsilon_v\phi(I_h+\eta R_h) + u_3(1-p) + \mu_v]\lambda_{S_v} \\ &\quad + (1-u_1)\lambda\varepsilon_v\phi(I_h+\eta R_h)\lambda_{E_v} \end{aligned}\\ - \frac{d\lambda_{E_v}}{dt}= \frac{\partial H}{\partial E_v} = -[ u_3(1-p)+ \alpha_v + \mu_v]\lambda_{E_v} + \alpha_v\lambda_{I_v} \\ \begin{aligned} - \frac{d\lambda_{I_v}}{dt}= \frac{\partial H}{\partial I_v} &= -(1-u_1)\beta\varepsilon_h\phi S_h\lambda_{S_h} + (1-u_1)\beta\varepsilon_h\phi S_h\lambda_{E_h} \\ &\quad - [u_3(1-p)+ \mu_v ]\lambda_{I_v} \end{aligned} \end{gather*} with transversality conditions \begin{equation} \lambda_{S_h}(t_f)=\lambda_{E_h}(t_f)=~\lambda_{I_h}(t_f) =\lambda_{R_h}(t_f)=\lambda_{S_v}(t_f)=\lambda_{E_v}(t_f)=\lambda_{I_v}(t_f)=0 \label{eT1} \end{equation} On the interior of the control set, where $0