A biological mathematical model of vector-host disease with saturated treatment function and optimal control strategies

: The aims of this paper to explore the dynamics of the vector-host disease with saturated treatment function. Initially, we formulate the model by considering three di ﬀ erent classes for human and two for the vector population. The use of the treatment function in the model and their brief analysis for the case of disease-free and endemic case are brieﬂy shown. We show that the basic reproduction number ( < or > ) than unity, the disease-free and endemic cases are stable locally and globally. Further, we apply the optimal control technique by choosing four control variables in order to maximize the population of susceptible and recovered human and to minimize the population of infected humans and vector. We discuss the results in details of the optimal controls model and show their existence. Furthermore, we solve the optimality system numerically in connection with the system of no control and the optimal control characterization together with adjoint system, and consider a set of di ﬀ erent controls to simulate the models. The considerable best possible strategy that can best minimize the infection in human infected individuals is the use of all controls simultaneously. Finally, we conclude that the work with e ﬀ ective control strategies.


Introduction
Mathematical modeling is used often for better understanding the infectious diseases dynamics. Mathematical models not only describe the mathematical process of the infectious diseases but also give useful information about the disease possible control and spread. There are a lot of infectious diseases in our world by providing many infected cases and death around the world. In which, vector-borne diseases are regarding a major threat to the human health causes many death and infection each year.
Vectors are the biological agents, which are the considered to be the main source of infection in human society. Dengue, malaria etc are the important vector-borne diseases that provide many infections and death cases to the human world. Vector-host disease mostly targeted the children especially in the developed countries. Some of the symptoms such as joint pains, headache, muscle, fever, and a skin rash similar to the measles. It is documented in infected cases few number of cases become the life-threatening dengue hemorrhagic fever. This results to bleeding, blood platelets with low level and with leakage of blood plasma, or the dengue shock syndrome in which a low blood pressure occurs. Approximately, one million deaths occurs per year and with over all 17% in all infectious diseases, so vector-borne diseases are considered to be responsible. Due dengue only, approximately 2.5 billion people in almost 100 countries of the world are currently at risk. Similarly, each year globally, 0.4 millions deaths are reported due to malaria, in all these cases most of the children are under the age of 5 years. Besides this, the infectious diseases such as leishmaniasis, chagas disease, and schistosomiasis provide millions of cases to the human population globally. Due to dengue, malaria, human African trypanosomiasis, yellow fever, schistosomiasis, onchocerciasis, and Japanese encephalitis contributed more than one billion cases and due to these a lot of deaths recorded/discovered globally [1].
Regarding the infectious diseases, the prevention is a useful protective tool to safe the society from infection whenever there is no vaccine or treatment. Dengue is a contagious disease caused by a virus, which is still epidemic in many regions such as tropical and sub-tropical areas of the world [2]. The disease is common in South Asia, Africa, USA and Western Pacific regions. There were 9 countries before 1970, which faced this problems but the increment were four times larger after 1995 [3]. The report of World Health Organization (WHO) suggested that dengue fever cases per year vary 50 to 100 million cases with approximately 10000 children death due to bleeding caused by dengue [4].
In order to understand the mechanism of vector-host diseases, the researchers developed numerous mathematical models in literature. For example, a mathematical model suggested by Ross [5] and then it is extended by the authors in [5] suggested the modeling and its analysis with optimal control analysis. The vector-borne disease transmission can be horizontally or vertically. A vector-borne disease model with time delay has been analyzed in [6]. The dynamics of vector-host model is studied in [7]. The analysis of vector-host disease with demographic structure is considered in [8]. The analysis of dengue dynamics with different mode of transmission is studied in [9]. The phenomenon of backward bifurcation analysis in dengue dynamics is considered in [10]. A vector-host disease with direct transmission is considered in [11]. Computer simulations and modeling formulation of dengue fever is analyzed in [12]. The dynamics of dengue infection in Pakistan with optimal control strategies has been proposed in [13]. Dengue dynamics with variable population is discussed in [14]. The authors considered in [15], the dynamics of malaria disease and presented the optimal control analysis with different control strategies. The dynamics of vector-host model with delay differential equation is studied in [16]. A dynamical model of vector-host disease with analysis of backward bifurcation and optimal control is considered in [17].
We aim here to formulate a mathematical model for vector-host dynamics through saturated treatment function. The use of the treatment function in mathematical models have been used by many authors, see [19,20,21,22]. The authors in [19] used the saturated treatment function in SIR model related to the network and presented the bifurcation analysis. The saturated treatment function for the age structured viral infection is analyzed in [20]. The rumors spread dynamics in social network is studied in [21]. The dynamics of SIR model with age dependent susceptibility with nonlinear incidence rate is investigated in [22]. We first develop the model and present mathematical results briefly. Then, we formulate a control problem and suggest a set of control combinations for possible control of infection. The rest of work is as follows: brief model formulation is given in Section 2. Model equilibria and its stability has been discussed in Section 3. Optimal control problem and its related results have been discussed in details in Section 4. The results are discussed briefly in Section 5 while the work is summarized in Section 6.

Model formulation
We present here briefly the dynamics of vector-host disease by denoting the total population of human by N h , subdividing further into three different classes, namely, the susceptible humans S h (t), infected humans I h (t) and the recovered humans R h (t) at any time t, thus N h (t) = S h (t) + I h (t) + R h (t). The population of susceptible human is increased by the recruitment of the individuals at a rate of Λ h . It is decreased by the effective contact with β 1 S h I v /(1 + α 1 I v ), where the disease contact rate between susceptible human and infected vector is represented by β 1 and α 1 is the saturation constant. It is further decreased by the natural death rate µ h . This rate of change can be represented through the following differential equation: The population of infected humans is generated by the effective contact rate β 1 S h I v /(1 + α 1 I v ) and decreased by the natural death rate µ h , the disease related death rate δ h and γuI h /(1 + buI h ). It can be seen when I or u is considered to be very small, then the treatment function converges to a near-zero value and whenever if the value of I is considered to be very large, then it approaches to a limit of finite value. The using of such type of function (treatment) will naturally reflect the epidemic system and thus, we consider it in our considered model. The term γ/b is defined to be the the maximal supply of medical resource per unit time while 1/(1 + buI h ) shows the reverse effect of infected people that are delayed for treatment and have an important effect on the disease spread, for details see [18]. We mathematically write the discussion in the form below: The individuals in the recovered class are generated by the treatment function γuI h /(1 + buI h ) while due to the natural death µ h it becomes decreasing. We mathematically obtain the following form: We denote the vector population by N v and distribute it into two subclasses, namely, S v susceptible vector and I v infected vector. Thus, we can write N v = S v + I v . The susceptible vector population is generated through the birth rate Λ v while decreased by the contact rate β 2 S v I h /(1+α 2 I h ) and the natural death rate µ v . This discussion leads to the differential equation given below: The infected vector population is generated through the contact rate β 2 S v I h /(1 + α 2 I h ) while decreased by the natural death rate µ v . This dynamics of the infected vector can be represented by the following differential equation: The equations (2.1-2.5) above can be written as a single system as follows: subject to the initial conditions (ICs) Let N h = S h + I h + R h , describes the dynamics of human population at time t and then it is given by According to the results that are given in Birkhoff and Rota [23], we have the following result: describes the total dynamics of vector at time t and then it is given by The exact solution of (2.10) is N v = Λ v µ v . The feasible region for the proposed model is is positively invariant.
Proof. To show the above result that is Ξ is positively invariant, we use standard comparison theorem As

Equilibria and local stability
In the present section, we examine the dynamics of the model (2.6) by the available fixed points. There exists two fixed points namely, the disease-free and the endemic equilibrium. We denote the disease-free equilibrium by P 0 and obtained the following: In order to find the stability analysis of the model (2.6), we need to calculate the basic reproduction number R 0 of the model (2.6) by considering the next generation method [24]. The desired matrices are computed as follows: The spectral radius R 0 = ρ(FV −1 ), that represents the basic reproduction number of (2.6), shown by .
Based on R 0 , the following are suggested: Theorem 3.1. If R 0 < 1, the disease-free equilibrium P 0 of the system (2.6) is locally asymptotically stable.

Endemic Equilibria
We obtain the endemic equilibria of the system (2.6) denoted by P * , and get, We have the following solution by using the above in first equation of the model (2.6): • If b or u are non-zero, then equation (3.1) becomes a quadratic equation with two roots for I h if C 1 < 0 and R 0 < 1. Also, if C 2 1 ≥ 4C 0 C 2 then there exists two positive roots, and namely the two • If b and u are both non-zero and R 0 < 1, then equation (3.1) has only one change of sign and so by the Descartes rule of sign it can be claimed that the system has a unique feasible equilibrium . Now, we have in the following the results for the local asymptotic stability of the model at P * 1 . Consider the theorem given below: Theorem 3.2. For R 0 > 1, then the vector-host system (2.6) at P * 1 is locally asymptotically stable. Proof. We obtain the Jacobian matrix below at P * 1 : where The coefficients involved in (3.4) are If the term (m 4 m 5 − m 2 m 3 ) > 0, then all coefficients k 1 , ..., k 4 become positive and then the Routh-Hurwitz condition k 1 > 0, k 3 > 0, k 4 > 0, and k 1 k 2 k 3 > k 2 3 + k 2 1 k 4 can be satisfied easily. The eigenvalues of the characteristics equation (3.3) then will have negative real parts if k i > 0 for i = 1, 2, 3, 4 > 0 and R 0 > 1 and the term under braces is positive. So, the result follows from Routh-Hurwitz criteria that the system (2.6) is locally asymptotically stable, if R 0 > 1 and the terms under braces is positive.

The existence of backward bifurcation
Here, we investigate the phenomenon of backward bifurcation for the system (2.6) using the center manifold theory described in [25]. Consider β 1 to be the bifurcation parameter and at R 0 = 1, we have Further, we make changes to the model variables by S h = y 1 , I h = y 2 , R h = y 3 , S v = y 4 , and I v = y 5 .
Using the vector notation y = (y 1 , y 2 , y 3 , y 4 , y 5 ) T , then, we write the model (2.6) in the form dy/dt = f , where f = ( f 1 , ..., f 5 ) is given by Evaluating the Jacobian matrix at P 0 with β 1 = β * 1 , we have It is obvious that a simple zero eigenvalues exists for the matrix J while the remaining have negative real part, so, it is possible now to apply the center manifold theory to the model (2.6). Next, we compute the left and right eigenvectors denoted by V = (v 1 , ..., v 5 ) and W = (w 1 , ..., w 5 ) and is given by Now, computing the values of a 1 and b 1 given by It can be seen that a 1 is negative while for the backward bifurcation a 1 and b 1 should be positive.

Global stability of DFE by Lyapunov method
We now consider the model (2.6) by obtaining its global stability at the disease-free and endemic case. First, defining the Lyapunov function for the model (2.6) at the disease-free case and present the result in the following theorem: The disease-free equilibrium of the model (2.6) is globally asymptotically stable, if R 0 < 1 and otherwise unstable.
Proof. In order to have the proof for the above result, we define the following Lyapunov function.
Now, by taking the time derivative of (3.6) and using the equations of the system (2.6), then we get and taking some arrangements of the terms, then we get Thus, by Principle [26], P 0 is globally asymptotically stable in Ξ.

Global stability of endemic equilibrium
We determine the global asymptotical stability of the model (2.6) by applying the geometric approach at P * 1 . To do this, we reduce the system (2.6) by using in the last equation of the model (2.6), and have the reduced system given by a new endemic equilibrium point P * 2 : subject to the non-negative initial conditions where g(x) : D → R n possesses a unique equilibrium x * and also a compact absorbing set exists for x * . Then, x * is globally asymptotically stable given that a function P(x) and a Lozinskii measure exists such that q = lim t−→∞ sup 1 t t 0 (H(x(s, x)))ds < 0 [27,28], where the symbols P , and H shall be defined in the result below.
Proof. The Jacobian matrix evaluated at P * 2 of the model (3.9) is given by Related to the matrix J, we define the following second additive compound matrix: where P f in the direction of vector field f shows the derivative of P. More precisely, we have: We obtain the matrix as follows: where Let the vector (û,v,ŵ) in R 3 and its norm . is defined as (û,v,ŵ) = max{|û|, |v|, |ŵ|}.
Let µH denote the Lozinski measure with the norm defined above. It follows from [27,28], we have |H 21 | and |H 12 | show the matrix norm related to the vector and µ, denote the Lozinski measure with respect to norm, then (3.11) Therefore, .

Now using system (3.9),
Then, we get Also, In f 2 above, we used the third equation of the system (3.9).
Then, we can get Then, This implies that q ≤ − µ 2 < 0. So, it follows from [27] that considered system is globally asymptotically stable.

Optimal control problem
This section investigates the application of the optimal control technique to the system (2.6) by modifying the birth rate of susceptible human and vector by the assumptions of the density effects Λ h → Λ h + cN h and Λ v → Λ v N v , while the constant c shows the density impact on the birth rate. Our main purpose is to formulate an optimal control problem and provide the best possible strategies of control for minimization of infection in human population. The use of optimal controls to the biological models with brief analysis are used by many researchers, see [29,30,31,32,33]. Here, in the optimal control system, we consider four controls, u i for i = 1, 2, 3, 4, which are defined as follows: u 1 is defined to be drugs or vaccine which can decrease the human and mosquitoes contacts such as insect repellents, the second control u 2 shows the level of larvicide and adulticide utilized in order to control mosquitoes breading places, the third control u 3 shows the minimization of human and mosquitoes contacts by the use of bed nets as a preventions and the control variable u 4 represents the control (through some specific prevention or treatment). The term (1 − u 1 ), is considered for the reduction of the force of infections in human population and b o is a positive rate constant. The factor (1 − u 2 ) is used for the reduction of the reproduction rate of mosquito population. The discussion above leads to the following control system: with the ICs (2.7) . In optimal control system (4.1), we considered the controls u(t) = (u 1 , u 2 , u 3 , u 4 ) ∈ U with a brief discussion. The control variables u(t) = (u 1 , u, u 2 , u 3 ) ∈ U are subjected to the state variables S h , E h , I h , S v and I v which are measured and bounded with We have the objective function for the vector-host control problem, given by, where U is defined in equation (4.2) and subjected to the system (4.1) with non-negative initial conditions. Consider the technique Pontryagin's Maximum Principle, to get the desired solution of the optimality system mathematically.

Existence of the optimality system
We use the results given in [34] for the control problem existence. The equations of the control (4.1) are bounded, which enable us to apply the result in [34] to our problem, if the following conditions are met: 1. O 1 : The state and control variables are nonempty.
To show these conditions (C 1 − C 4 ), we follow the results of [35] to find the result for the existence of (4.1). The controls and the state variables are clearly bounded, which confirms O 1 . The claim O 2 is confirmed because of bounded solution and convex. In order to fulfill C 3 , the model is bilinear in control variables. The last claim O 4 and their verification is given below, where D i , l 1 , l 2 > 0 and m > 1 for i = 1, ..., 6. Thus, we have 2) subject to the optimality (4.1), then there exists an optimal control u * = (u * i ) such that J(u * i ) = min U J(u i ) for i = 1, ..., 4. For the solution of an optimal control problem, the construction of the Lagrangian L and the Hamiltonian H is required, which are defined below: and X = (S h , I h , R h , S v , I v ), U = (u 1 , u 2 , u 3 , u 4 ) and λ = (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 ), to get: (4.6)

Optimal control problem solution
Using Pontryagin's Maximum Principle [36] for the solution of the optimality system as follows: suppose u * i for i = 1, ..., 4 denote the optimal solution of the optimality system (4.1), then the adjoint variables say, λ i for i = 1, ..., 5 exists which satisfy the conditions below, (4.7) Using these conditions to H, the following are obtained: For the controls u * i for i = 1, ..., 4 and S * h , I * h , R * h , S * v , I * v represent the solution of system of state, then there exists adjoint variables, say, λ i for i = 1, ..., 5, with transversality conditions Further, the control u * i for i = 1, ..., 4 are

10)
Proof. To obtain the results stated in above theorem, we solve the control system together with the Hamiltonian H (4.6) to have the results for the adjoint system (4.8) and the transversality conditions (4.9), with setting S h = S * h , I h = I * h , R h = R * h , S v = S * v , and I v = I * v and the derivative of H with respect to S h , I h , R h , S v , I v , we have (4.8). To get the equations of optimal control characterization in (4.10), we use ∂H ∂u i = 0, for i = 1, 2, 3, 4.

Discussion
This section describes the numerical results of the proposed model (2.6) and the optimal control problem (4.1), which are solved numerically. The optimal control solution is obtained through backward Runge-Kutta order four scheme. We denote the solution of the control system via dashed line and those without control by a bold line. The time unit considered in the numerical solution is per day. The numerical values for the parameters are presented in Table 1. The weight and balancing constants with their proposed values are D 1 = D 2 = 1000, D 3 = 10, D 4 = 0.005, D 5 = 0.03 and D 6 = 3. We choose different cases to investigate the optimal control solutions. We present the following cases: Case (i): In this case, we consider the control variable u 1 = 0 and make the rest of the controls u 2 = u 3 = u 4 0 and simulating the model. The resulting graphical results are depicted in Figure  2 with subfigures (a-f). In this case, the population of infected humans decreases and the recovered human increased. Also, the vector populations decrease sharply. This case effective for the infected population as it decreases very fast after day 14 and becomes steady.
Case (ii) In this case, we set u 2 = 0 and u 1 = u 3 = u 4 0. The resulting graphical results are presented through Figure 3 with subfigures (a-f). In this case, the population of susceptible human less increased compared to Case (i), but no decrease in the population of infected vector and susceptible vector. Although the population of the recovered and infected human the same as in Case case (ii). Thus, the strategy is not a good one.

Case (iii):
In this case, we set u 3 = 0 and u 1 = u 2 = u 4 = 0. The resulting graphical results are presented through Figure 4 with subfigures (a-f). In this case, the population of susceptible humans increases sharply compared to Case (i) and (ii). The population of infected humans, infected vector and susceptible vector is increasing more compared to previous strategies. Also, the population of recovered human increases.
Case (iv): In this case, we choose to set u 4 = 0 and u 1 = u 2 = u 3 0 and simulate the model and obtain the results graphically given in Figure 5 with subfigures (a-f). Comparing the control system with and without controls system, the population of susceptible individuals increases and decreases the population of infected but it can be seen that there is no increase in the population of recovered individuals. It can also be observed that this strategy minimizes the infection in the vector population.

Case (v):
In the above combinations of the controls and their simulations, which suggest the increase or decrease in different compartments of the humans and vector populations. In all these strategies from (i-iv) no one get the desired results for the humans and vectors population to be minimized as desired. So, we utilize all the controls active and simulate the model of control in connection with the model having no controls. We observe that this set of controls provide that the population of susceptible and recovered human increase sharply while the population of infected humans, susceptible vector and infected vector are decreasing better, see Figure 6. Comparing to the above cases, this strategy is comparatively better.

Conclusions
We presented a mathematical model for vector-host disease with saturated treatment function and presented its dynamical results with optimal controls analysis. Stability analysis of the model for the disease-free and endemic cases are obtained and discussed. The disease-free equilibrium is stable when the basic reproduction number R 0 < 1. When the basic reproduction number R 0 > 1, then the endemic equilibrium found to be stable both locally and globally. The optimal control problem with controls variables are formulated and the desired results are obtained and discussed briefly. The optimal control problem together with controls function, and with adjoint equations are simulated and the results of both the models with and without controls are showed. A set of different controls were used to obtain the graphical results and we found that the Case (v) is considered to be the best strategy to control the infection in humans. The use of saturated treatment function in the modeling of vector-host disease is a novel practice and could be useful for the mathematicians and scientists working on vector-host diseases research.