OPTIMAL CONTROL OF VECTOR-BORNE DISEASES: TREATMENT AND PREVENTION

. In this paper we study the dynamics of a vector-transmitted disease using two deterministic models. First, we look at time dependent pre- vention and treatment eﬀorts, where optimal control theory is applied. Using analytical and numerical techniques, it is shown that there are cost eﬀective control eﬀorts for treatment of hosts and prevention of host-vector contacts. Then, we considered the autonomous counter part of the mode and we estab- lished global stability results based on the reproductive number. The model is applied to study the eﬀects of prevention and treatment controls on a malaria disease while keeping the implementation cost at a minimum. Numerical results indicate the eﬀects of the two controls (prevention and treatment) in lowering exposed and infected members of each of the populations. The study also highlights the eﬀects of some model parameters on the results.


1.
Introduction. Vector-borne diseases have been and are among the leading causes of death that remain challenges for many countries in the world, claiming lives of millions of people every year. The health as well as the socioeconomic impacts of emerging and re-emerging vector-borne diseases are significant. Thus, they are among the major concerns for several health organizations including, the World Health Organization, Center for Disease Control and National Health Institute [30,31,32]. Vectors are found in areas ranging from tropical to temperate zones and at different landscapes, thus, the geographic distribution of vector-borne diseases is vast and diverse.
Dengue fever, which is transmitted by Aedes aegypti mosquito, is common in tropical climates including the southern parts of the United States in the Texas and Florida regions [27]. Malaria is carried by female Anopheles mosquito, and it is prevalent mainly in Africa and some parts of Asia [11,12,19,25]. On the other The latent group in both populations is infected, but it does not transmit the disease. Members of the latent groups are at incubation stage, whereas the infected group is infectious in each case. We denote the total population size of the host by N = Σ 4 i=1 x i and the vector by P = Σ 3 i=1 y i . The transmission dynamics of the vector-borne disease is given by with initial conditions x i (0) ≥ 0, i = 1, · · · , 4 and y i (0) ≥ 0, i = 1, 2, 3. Note that 1 − u 1 (t) describes the failure rate of prevention efforts. The per capita recovery rate is r 0 u 2 (t), where 0 ≤ r 0 ≤ 1 is the proportion of effective treatment.
The rest of the functions and parameters in the model are defined as follows. δ 1 , δ 2 : per capita recruitment rate of the host and the vector populations, respectively; density dependent factors for the recruitment rates are g(N ) (for the host) and h(P ) (for the vector) populations; β : the probability that susceptible hosts become infected after contact with an infected vector; φ : the average number of contact made to a host by a single vector per unit time. Thus, φP N is the contact rate of vectors per host in a unit time. As y3 P is the proportion of infectious vectors, x 1 φP N y3 P describes the total number of contacts between infectious vectors and hosts per unit time. Finally, βx 1 φP N y3 P is the infection rate of hosts. However, the prevention effort u 1 (t) reduces the infection rate with a failure probability of (1 − u 1 (t)), which means if prevention controls are practiced, then the infection rate is βx 1 The immunity of treated hosts is lost at a rate of ψ per unit time, making treated individuals susceptible to the disease. Our model assumes that infectious individuals at acute stage could die at a rate of α. The natural recovery rate r is assumed to be larger than the disease induced death rate, d. Note that when there is no treatment, 1 r is the mean length that a host remains infectious before it naturally recovers to join the recovered group, x 4 .
The probability that susceptible vectors get infection after contacting infectious hosts is given by θ. Thus, since x3 N is the proportion of infectious hosts, assuming that vectors get infected by contacting only infectious hosts, the term θφx3y1 N (1 − u 1 (t)) describes the infection rate of vectors through contact with infectious hosts.
Since this model assumes that no genetic mutation or drug resistance of the disease takes place, the disease is caused by a single strain. Moreover, the community is assumed to be closed when it comes to immigration, but the general exit terms; µ in the host and γ in the vector, could include emigration rates as well as death rates. Therefore, 1 µ and 1 γ are respectively, the duration of time a host and a vector stay in the community until it dies or emigrates elsewhere. Further explanations of model parameters along with their estimates based on a malaria disease are given in Section 4.
It is biologically meaningful to assume that the recruitment functions g(N ) and h(P ) satisfy conditions g(0) = 0, h(0) = 0, g(N ) ≥ 0, and h(P ) ≥ 0. We further assume that g and h remain bounded on their domains with upper bounds k 0 and v 0 , respectively. Density dependent factors such as competition for resources could contribute to bounded recruitment functions. The following discussion could further imply that under these circumstances, the host and vector populations remain bounded. From (1) we have that Thus, for initial value P (0) ≤ δ2v0 γ , we have P ≤ δ2v0 γ . Let X = (x 1 , x 2 , x 3 , x 4 ) and Y = (y 1 , y 2 , y 3 ). Based on the above discussion, we define a set Ω as In the remaining part of this paper we restrict our state variables to this set. Proof. Clearly for H 1 = −(βφ + µ), H 2 = −(d + µ), H 3 = −(r + r 0 + α + µ) and Thus, solutions with initial value in Ω remain nonnegative for all t ≥ 0. Moreover, as (2) and (3) that N ≤ δ1k0 µ and P ≤ δ2v0 γ . Therefore, Ω is positively invariant under system (1). The analysis of the optimal control is carried out for time dependent N and P. However, we approximate the reproduction terms δ 1 g(N ) and δ 2 h(P ) by their average values throughout this paper.
Together with the mathematical model described by equation (1), we consider an optimal control problem with the objective (cost) functional given by A 1 and A 2 are weight constants of the latent and the infected group, respectively, whereas, B 1 and B 2 are weight constants for prevention and treatment efforts, respectively which regularize the optimal control. For technical purposes, it is assumed that the cost of prevention and the side effects of treatment are given in quadratic form in the cost function (5), that is, B 1 u 2 1 is the cost of prevention, and B 2 u 2 2 is the cost of the treatment effort. The cost of treatment could come from cost of drug, any cost associated with treating patients with other health complications, surveillance and follow up of drug management and fighting emergence of drug-resistant strains. Similarly, the cost of prevention is associated with costs of vaccination, pesticide sprays, educating the public about personal protection and supply of basic needs (treated bed screen, for example). Furthermore, l(x 1 (T ), x 4 (T )) is the fitness of the susceptible and the treated group at the end of the process as a result of the prevention and treatment efforts implemented.
We seek an optimal control pair (u * 1 , u * 2 ) such that J(u * 1 , u * 2 ) = min{J(u 1 , u 2 ) | (u 1 , u 2 ) ∈ Γ} subject to the system given by (1) and where the control set is Here a i and b i , i = 1, 2 are constants in [0,1]. The basic framework of an optimal control problem is to prove the existence of the optimal control and then characterize the optimal control through the optimality system.

2.1.
Existence of an optimal control. We note that the existence of an optimal control pair can be proved by using results from Fleming and Rishel, [14]. It is clear that the system of equations given by (1) is bounded from above by a linear system. The boundedness of solutions of system (1) for a finite time interval is used to prove the existence of an optimal control. We use the results on existence, [14,Theorem 4.1, in the special case of an optimal control for a free terminal point problem, where the initial time and state along with the final time are fixed, and there are no conditions on the final state. To use the theorem in [14], we first check the following properties.
1. The set of controls and corresponding state variables is non-empty. 2. The control set Γ is convex and closed. 3. The right hand side of the state system is bounded by a linear function in the state and control. 4. The integrand of the objective functional is convex on Γ and is bounded below by c 1 (|u 1 | 2 + |u 2 | 2 ) β 2 − c 2 , where c 1 , c 2 > 0 and β > 1. 5. The function l is continuous. An existence result in Lukes ([21], Theorem 9.2.1) for the state system (1) with bounded coefficients is used to give condition 1. The control set Γ is convex and closed by definition. The right hand side of the state system (1) satisfies condition 3 as the state solutions are a priori bounded (see Theorem 2.1). The integrand in the objective functional, , is clearly convex on Γ. Moreover, there are c 1 , c 2 > 0 and β > 1 satisfying , because, the state variables are bounded. Finally under assumption 5, there exists an optimal control pair (u * 1 , u * 2 ) that minimizes the objective functional J(u 1 , u 2 ). Note that we expect the fitness function l(x 1 (T ), x 4 (T )) to be maximized while the cost function J is minimum, for example, if we define the function l as follows; Then the function l is clearly continuous.
2.2. The optimality system. With the existence of the optimal control pair established, we now present the optimality system using a result from Lewis and Syrmos [20]. We discuss the theorem that relates to the characterization of the optimal control. The optimality system can be used to compute candidates for optimal control pair. To do this, we begin by defining a Lagrangian (which is the Hamiltonian augmented with penalty terms for the control constraints). Let where Ω is defined by (4) and Γ is defined by (6).
The Lagrangian for our problem consists of the integrand of the objective functional and the inner product of the right hand sides of the state equations and the adjoint variables Π = (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , λ 7 ). Penalty multipliers, w ij (t) ≥ 0, are attached to the control constraints. We define the Lagrangian as follows.
where w ij (t) ≥ 0 are the penalty multipliers satisfying Theorem 2.2. Given an optimal control pair, (u * 1 , u * 2 ), and solutions, x 1 , x 2 , x 3 , x 4 , y 1 , y 2 , and y 3 of the corresponding state system (1), there exist adjoint variables Π satisfying with the terminal conditions, Furthermore, u * 1 and u * 2 are represented by Proof. First of all, we differentiate the Lagrangian, L, with respect to states and then the adjoint system can be written aṡ The terminal conditions (9) of adjoint equations can is given by [20, pp.134] To obtain the optimality conditions (10), we also differentiate the Lagrangian, L, with respect to U = (u 1 , u 2 ) and set it equal to zero.
Solving for the optimal control, we obtain To determine an explicit expression for the optimal control without w 11 and w 12 , a standard optimality technique is utilized. We consider the following three cases.
This implies that (iii) On the set {t| u * 1 (t) = a 1 }, we have w 12 (t) = 0. Hence This implies that Combining these three cases, the optimal control u * 1 is characterized as Using similar arguments, we also obtain the second optimal control function We point out that the optimality system consists of the state system (1) with the initial conditions, Z(0), the adjoint (or costate) system (8) with the terminal conditions (9), and the optimality condition (10). Any optimal control pair must satisfy this optimality system. The graphical representations of the optimal control functions u 1 and u 2 are given in Figures 1. Parameter values are estimated based on a malaria disease as given in Table 1. 3. The autonomous model. In this section we analyze the stability of the diseasefree equilibrium by first calculating a key threshold quantity vital in epidemiology, and then study the existence and stability of fixed points results.
Without considering the control functions in (1) (u 1 = 0, u 2 = 0), the autonomous version of the model is given by If we approximate the reproduction function δ 1 g(N ) by a constant (possibly average) value Λ, then in the event when the community is disease free (x i = 0, i = 1, y i = 0, i = 1), we have from the first four equations in (11) Thus, for a low level of disease induced death rate (α ≈ 0), the population could eventually assume a steady-state value. Motivated by this, we consider a host population which assumes a steady state value Λ µ . Similarly, from the last three equations of (11), we see that the vector population equilibrates to ∆ γ as well, where δ 2 h(P ) is replaced by it average value ∆. We consider this value for the size of the vector population.
For technical purposes, in this section, we assume that both populations are at steady state. Clearly, the above system (11) has a unique disease-free equilibrium given by Using the next generation approach [7,10], we estimate a threshold quantity which is critical to the asymptotic stability of the disease-free equilibrium point. To this end, we consider the equations for the infectious group in the host and the vector populations, that is, the third and the seventh equations in (11), and let

KBENESH BLAYNEH, YANZHAO CAO AND HEE-DAE KWON
With X = (x 1 , x 4 , y 1 ), Y = (x 2 , y 2 ) and Z = (x 3 , y 3 ), the model (11) can be split into three equations as where f 3 = (h 1 , h 2 ) T (T is transpose) is given by (14). Let U 0 = (X * , 0, 0) ∈ R 7 denote the disease-free equilibrium of the system given by (11), where The constants used in the foregoing functions are: W = dβφ d+µ , E = ψr µ+ψ , F = r + α + µ, G = µN, H = kφθ∆ γ+k , J = γN. Then, the maximum eigenvalue, also known as the spectral radius of M D −1 is where N = Λ µ and P = ∆ γ . This threshold quantity is the basic reproduction number as defined in Diekmann and Heesterbreek [10]. To see this, note that the term kφβ γ(γ+k) describes the number of hosts that one vector infects (through contact) during the lifetime it survives as infectious, when all hosts are susceptible. On the other hand, the term dφθ (d+µ)(r+α+µ) describes the number of vectors that are infected through contact with one infectious host, while the host survives as infectious, assuming no infection among vectors.
Proof. The Jacobean of system (11) The third and the forth columns have diagonal entries. Therefore, the diagonal entries −(ψ +µ) and −γ are two of the eigenvalues of the Jacobian. Thus, excluding these columns and the corresponding rows, we calculate the remaining eigenvalues.
These eigenvalues are the solutions of the characteristic equation of the reduced matrix of dimension four which is given by To simplify the notation, we let where The Routh-Hurwitz conditions [22], which usually have different forms are the sufficient and necessary conditions on the coefficients of the polynomial (18). These conditions ensure that all roots of the polynomial given by (18) have negative real parts. For this polynomial, the Routh-Hurwitz conditions are: Clearly H 4 = B 0 H 3 . Since A i > 0, i = 1, 2, 3 we have B i > 0, i = 1, 2, 3. Moreover, if R 0 < 1, it follows that B 0 > 0. Thus, it is enough to prove that H 2 > 0 and Using Maple, it is easy to see that which is positive. Again, using Maple we can also see that which is clearly a positive quantity. Therefore, all of the eigenvalues of the Jacobian matrix have negative real parts when R 0 < 1. However, R 0 > 1 implies that B 0 < 0, and since all of coefficients (B 1 , B 2 andB 3 of the polynomial (18) are positive, not all roots of this polynomial can have negative real parts. This means, when R 0 > 1, the disease-free equilibrium point is unstable.
Note that the result in Lemma 3.1 is local, that is, we could only conclude that solutions with fairly small initial size in the invariant set Ω are attracted to the disease-free equilibrium point. It is possible to further reduce the dimension of the Jacobian in the proof of Lemma 3.1 by using y 1 = P − (y 2 + y 3 ) and x 1 = N − (x 2 + x 3 + x 4 ) without any technical difficulty.
Proof. We begin the proof by defining new variables and breaking the system given by (11) into subsystems. With Z = (x 2 , x 3 , y 2 , y 3 ) and X = (x 1 , x 4 , y 1 ), this system can be written as The two vector-valued functions G(X, Z) (dimension 4) and F (X, Z) (dimension 3) are given by where T is transpose. Now consider the reduced system: dX dt = F (X, 0) is a GAS equilibrium point for the reduced system dX dt = F (X, 0). To see this, solve the second equation in (24) to obtain x 4 (t) = x 4 (0)e −(µ+ψ)t ; it approaches zero as t → ∞. Similarly from the last equation we get y 1 (t) = ∆ γ + (y 1 (0) − ∆ γ )e −γt which approaches Λ γ as t → ∞. Finally, the first equation yields This asymptotic dynamics is independent of initial conditions in Ω. Hence, the convergence of solutions of (24) is global in Ω. Clearly G(X, Z) satisfies the following two conditions given as H2 in [7, pp. 246] namely, (i). G(X, 0) = 0 and (ii). G(X, Z) = AZ − G(X, Z), G(X, Z) ≥ 0 on Ω, where Note that since the host and the vector populations assume a steady-state value (N = Λ µ and P = ∆ γ ), the term φθx 3 ( ∆µ Λγ − y1 N ) in G(X, Z) is nonegative. Moreover, by Lemma 3.1 the disease-free equilibrium point is LAS for R 0 < 1. The global stability of the disease-free equilibrium point follows from the theorem in [7, pp. 246].
Theorem 3.3. The autonomous system given by (11) has a unique endemic equilibrium in Ω if R 0 > 1.
Proof. It is easy to see that the steady-state equations for x 1 , x 2 , x 3 and x 4 along with Also, from the equilibrium equations for y 1 , y 2 and y 3 , with y 1 = P − (y 2 + y 3 ) we get To simplify notations, let B = φθ N (r+α+µ) and C = γ+k γ . Then, the first equation in (26) yields Replacing y 3 by its new value (see the second equation in (26)) transforms equation (25) to It is easy to see that if P − Cy 2 = 0, then y 2 + y 3 = P which further leads to y 1 = 0 because, y 1 + y 2 + y 3 = P. As a result of this, dy1 dt = 0 implies that h(P ) = 0, which means that P = 0. But this is possible only if y 2 = 0 = y 3 , which is not of interest to us (there are vectors in the environment). Thus, y 1 > 0 and P > y 2 + y 3 = Cy 2 at steady-state. Now as P − Cy 2 > 0, we have x 2 > 0 as long as y 2 > 0, which we verify in the following. From the simplifying notations given above and (27), equation (25), which is also (28) gets its final form A unique nonzero solution of this equation satisfies Note that y 2 > 0 iff Γ = dφ 2 βθkP γN (r+α+µ)(d+µ)(γ+k) > 1, but R 2 0 = Γ. Thus, a unique endemic equilibrium point is possible only if R 0 > 1.

4.
Numerical results and discussion. The model considered in this paper is tested with data taken from a malaria disease of a single strain. This disease spreads via a cross-infection between vertebrate hosts including human and invertebrate vector which is the female anopheles mosquito. In the following two sub sections we discuss simulation results for the optimal control and the autonomous models for different parameters. The parameters used in solving the optimality system are estimates of data from malaria [28]-[31] as summarized in Table 1 which is presented at the end of this section.

4.1.
Simulation of the optimal control model. The optimality system is a twopoint boundary problem, because of the initial condition Z(0) of the state system (1), and the terminal condition Π(T ), (9) of the adjoint (or costate) system (8). We implemented a gradient method to solve the optimality system numerically. After making an initial guess for the control functions, we first solved the initial valued state system forward in time. Then, using the same initial guess for the control functions, we solve the adjoint system with terminal conditions backward in time. The controls are updated in each iteration using the optimality conditions (10). The iterations continue until convergence is achieved.
We balance the host populations and control functions by choosing weight constant values A 1 = A 2 = 1, B 1 = B 2 = 50, and Q 1 = Q 2 = 0.1 in the objective functional (5), because, the magnitudes of the host populations and the control functions are on different scales.
The numerical results are obtained for different values of φ, the average vectorhost contact rate, while keeping the remaining parameters unchanged in each simulation. At first we search for two optimal control functions u 1 and u 2 for prevention and treatment, respectively. These optimal control functions are designed in such a way that they minimize the cost function given by (5). Among outcomes of this, we find minimal latent and infectious host population levels. An increase in vector-host contact rate (φ) is an indication that prevention effort has lower efficacy. Therefore, in this case, to optimize the number of infectious and latent hosts, one has to increase the treatment rate. This can be noticed from comparison of φ = 3 and φ = 10 in Figure 1. After setting u 2 = 0 (no treatment), we search for optimal control function for prevention, u 1 . We repeat this process to search for optimal treatment control function, u 2 , setting u 1 = 0 (see Figure 2, φ = 3 and Figure 3, φ = 10). The cost of treatment is higher for larger contact rates.
As it can be seen in Figure 4, if we have to use only one control, then, prevention is more effective than treatment to diminish the size of latent and infectious hosts. This highlights the effectiveness of preventive measures in controlling vector-borne diseases. In perspective, one could conclude from the controls in Figures 1, 2 and 3 that we should give full prevention effort in the beginning of emergence of the disease while giving full treatment effort in the middle of the time interval when control efforts are practices. This means that prevention is more important in the beginning of disease outbreak. On the other hand, treatment is more important while the disease prevails.
As it is depicted in Figures 4, 5, 6 and 7, the latent and infectious population levels obtained using both controls (prevention and treatment) are lower than their counter parts which result from practicing only one control (prevention or treatment). Our conclusion from this is: intervention practices that involve both prevention and treatment controls yield a relatively better result.
Another notable feature in Figures 6 and 7 is that when the average host-vector contact rate (φ) is increased (preventive control has poor efficacy), then a combined control effort yields a result that is not remarkably different from those obtained using only one of the controls.

4.2.
Simulation of the autonomous model. We performed simulations for the autonomous model (11) for two cases; R 0 < 1 and R 0 > 1. R 0 is the basic reproduction number given by Equation (16). These simulations are performed for different values of φ and initial conditions. Some of the model parameters are considered based on a malaria disease (see Table 1). The simulation results are given in Figures 8,9 and 10. It is clear from Lemma 3.1, that the disease is endemic for R 0 > 1. Adding to this, the numerical simulation shows that the equilibrium levels of infectious and latent individuals in both populations increase when the contact rate, φ increases (See Figure 9).
The disease-free equilibrium point is globally stable for R 0 < 1 as it is given by Theorem 3.2. However, as we can see in Figure 8, although the size of all infected groups in each population dies out, its magnitude increases for a while with increased contact rates before it goes to zero. Besides the increased vector-host contact rate φ, increase in vector survival rate γ, is also a factor for the increased number of infectious hosts (see Figure 11). A vector-borne disease can be eliminated through a combination of effective vector control and personal protection from vectors. It is also important to note that if the contact rate is kept low, then an increase in infected (latent and infectious) vectors does not necessarily result in the spread of the disease among hosts. This is shown by the phase plane portrait for two different contact rates φ = 3 and φ = 10 in Figure 12.
Parameters for Malaria: There are different types of malaria causing viruses; what we consider in the simulations of the two models are basically the average values which encompass many of the features of the malaria disease including rate of infection, incubation period, length of infectious period in the vector and the host populations. Although estimates of some parameters are given in Table (1), here we give additional explanation and description of some of the parameters as they are given from reports by the Center for Disease Control and prevention and World Health Organization .
(1). The disease induced death rate α is allowed to vary. Most malaria parasites are not highly fatal, with the exception of chloroquine-resistant P. falciparum, thus, we keep α small. (2). The average life expectancy of adult mosquito, 1 γ , is estimated based on the range 15 to 20 days. For humans, the natural death rate is estimated based on average life of 60 years in developing countries. (3). It takes humans 10-30 days to develop malaria after a bite by infectious mosquito occurs (the average period is 20 days) [28]. Therefore, a reasonable estimate of the incubation rate d is 0.05 ≤ d ≤ 0.059. (4). The recovery rate could vary based on the gametocyte, and it may take an average of 30 days for the parasite to be cleared from bloodstream after treatment. We estimate the maximum recovery rate 1 r+r0 to be about 1/30. (5). It takes a female anopheles mosquito 10 to 26 days to develop the malaria parasite in its salivary gland after blood meal from infectious host [29]. Thus, a reasonable value for k could ranges from 1/26 to 1/10. Table 1. Estimated parametric values in the model for malaria case. The rates are given per day (6). After treatment, a person could take 1 to 20 years to be susceptible to the disease, therefore, we consider ψ to assume a value from 1 356 to 1 (20)(356) per day. In the simulation, we consider initial conditions x 1 (0) = 100, x 2 (0) = 20, x 3 (0) = 20, x 4 (0) = 10, y 1 (0) = 1000, y 2 (0) = 20, and y 3 (0) = 30. We bound prevention and treatment efficacy by a 1 = a 2 = 0 and b 1 = b 2 = 1.

5.
Conclusion. We wrap up our study by making some concluding remarks about our findings. As mentioned in Section 2, time dependent intervention strategies have been implemented to curtail a vector-borne disease on a finite time interval. This model allows intervention activities to be intensified in appropriate time intervals relevant to disease outbreak, which in most cases is in phase with climatic variations. Since the control functions are allowed to be piecewise continuous, the model incorporates possible scenarios including termination of activities in some parts of the year. The model also predicts that the time dependent efforts practiced to contain the disease outbreak could be cost effective: the disease control could be performed with minimal cost of prevention and treatment. Note that the same model could be used to predict minimal environmental side effects, such as emergence of drug resistant vectors, that result from preventive measures.
Prevention methods, which range from monitoring environmental changes to reduction of direct contact between vectors and hosts are among the best bets. While addressing the seasonality of the vector-borne disease, the results of numerical simulations of our model support the hypothesis that among intervention techniques, preventive practices are very effective in reducing the incidence of infectious hosts and vectors.
Although control theory results in optimal strategies to curtail vector-borne diseases, it is not capable of predicting the long term dynamics of the disease as a whole. The long term prediction of the disease dynamics is carried out, as in most epidemic studies, using the basic reproductive number of the autonomous counter part of our model. This is addressed in Section 3 assuming that the total population level of the vector and the host is constant. When disease-induced death rate is ignored in the host population, it is possible for the total population to eventually equilibrate.
In any case, assuming that the disease-induced death rate is at a minimum level, we considered the total populations of the host and the vector to be at steady state. This is used to simplify the technical details while analyzing the dynamics of the autonomous version, but no such assumption is made in the control model. Exploring the impacts of the size of A 1 , A 2 , B 1 , B 2 , Q 1 , and Q 2 , along with the upper and lower bounds of the control functions, that is, a i and b i , i = 1, 2 on the optimal solution is part of our future research.   Figure 2. The optimal value of the cost functional is calculated using only one control function. The first graph is an optimal control function for prevention with no treatment. The second graph is an optimal control function for treatment with no prevention. In both case, the average number of host-vector contacts per unit time is φ = 3.  Figure 3. The optimal value of the cost functional is calculated using only one control function. The first graph is an optimal control function for prevention with no treatment. The second graph is an optimal control function for treatment with no prevention. In both cases, the average number of host-vector contact per unit time is φ = 10.