Optimal control of visceral, cutaneous and post kala-azar leishmaniasis

This article focuses on the eradication of different strains of leishmaniasis with the help of almost nonpharmaceutical interventions (NPIs). A comprehensive mathematical model of the disease is formulated incorporating three types of populations: sandflies, humans and dogs (reservoirs), and 3-types of strains: Cl\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$C_{l}$\end{document}, cutaneous leishmaniasis; Vl\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$V_{l}$\end{document}, kala-azar; and PKDL, post kala-azar. We find R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0}$\end{document}, the basic reproduction number of the infection. On the basis of sensitivity test of R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0}$\end{document}, the most active/sensitive parameters are investigated. These active parameters are controlled with the help of control variables. In some cases different parameters depend on the same single parameter, like ovigenesis and biting rate, both of which are linked to the blood source. Therefore we introduce three nonpharmaceutical control variables in the proposed model to control the biting rate of sandflies, density of seropositive dogs, and density of vector population. Nonpharmaceutical interventions include bed nets, eradiation of infectious dogs, and residual sprays, and thus extend the proposed model to an optimal control model. Using Lagrangian and Hamiltonian, we minimize the densities infected classes in human, sandfly and vector populations. Adopting optimality approach, we check the existence of the optimal control for the system. Using Matlab, we produce numerical simulations for the validation of results of control variables.


Introduction
Leishmaniasis is a group of vector-borne diseases. Cutaneous and visceral strains are two main types of the disease. The disease is caused by protozoa which belong to Leishmania genus. The protozoa is transmitted/carried from one individual to another by a vector/sandfly. The three causative pathogens of visceral leishmaniasis are L. infantum, L. donovani, and L. chagasi. The incubation/latency period varies from three to six months with a minimum of 10 days, however, a longer period of one year has also been reported [1][2][3][4]. Each year 500,000 new cases of V l strain are reported worldwide [5]. The causative agents of cutaneous strain are L. major and L. tropica, however, in some cases L. infantum may also cause C l strain [6].
In very rare cases, very few individuals who seem to be recovered from visceral strain develop complications of the disease. The period in which complications may appear is from

Model formulation
In what follows, N h is the density of the human class and is divided into the following further subclasses: • S h , the susceptible human class; • E 1 and E 2 , the human classes exposed/latent for C l and V l infection, respectively; • E 12 , the human class recovered from C l and exposed to the infection of V l ; • E 21 , the humans recovered from V l and exposed to the infection of C l ; • E 23 , the dormant/exposed class of PKDL, the class which otherwise seems recovered from V l ; • E 123 , the exposed class of PKDL after being recovered from C l ; • I 1 and I 2 , the infectious classes of C l and V l ; • I 12 and I 21 , the infectious classes of V l and C l , developing from E 12 and E 21 ; • P 2 , the infectious class of PKDL; • P 12 , the infectious class of PKDL after being recovered from C l ; • R 1 , R 2 , and M, the classes recovered from C l strain, V l strain, and both strains, respectively. Next N r is the total reservoir population, divided into the following subclasses: • S r , the susceptible reservoir class; • I r , the infectious class of reservoir; • Z r , the reservoir class being recovered from infection. Then N v is the total vector population, divided in the following classes: • S v , the susceptible-vector class; • E v , the exposed-vector class, and • I v , the infectious-vector class. A susceptible human after catching infection follows Path-1 or Path-2, according to as the infection causing fly or the source fly was infected with C l or V l , respectively: • Path-1, In case of complication of V l strain, the transmission follows Path-3 or Path-4: • Path-3, I 2 − → P 2 − → R 2 ; • Path-4, I 12 − → E 123 − → P 12 . The infection follows the following path of transmission in dogs (reservoir): • S r − →I r − → Z r . The infection follows the following path of transmission in vectors: The overall path of transmission of leishmaniasis in all the three populations is shown in the following Fig. 1. The mathematical model of the 3-strains of leishmaniasis is given by the following set of coupled differential equations: E 123 = γ 5 I 12 -(k 6 + μ h )E 123 , P 12 = k 6 E 123 -(β 2 + δ 3 + μ h )P 12 , (1) The terms of interaction 1 , 2 , r , and v between the three populations are defined as follows: ρ(c 1 I 21 + c 2 I 12 + c 2 P 12 ) + c 1 I 1 + c 2 (I 2 + P 2 ) + cI r . Table 1 contains the values of the different parameters used in the model (1).

Model analysis
Here, in Sect. 3, the properties of the model, namely the disease-free equilibrium, invariant region, and basic reproduction number, are discussed.

Proposition 1 The region defined by
is a positively invariant domain, also the model is epidemiologically and mathematically well-posed as all the trajectories are forward-bounded.

Sensitivity analysis of R 0
Different parameters used in the model influence the evolution of the disease differently. The role of parameter K in the phenomenon Z is called the sensitivity of Z with respect to K and is given by [18,31] The sensitivity indices of different parameters are shown in Table 2

Control strategies
To control the further transmission of the disease in the community, we use sensitivity test to check the absolute value of the index of a parameter. The greater the absolute value of the index of a parameter, the higher would be the role of this parameter. The contact rate of the fly has sensitivity index of 1. This index is the highest, hence a plays the most significant role in the disease transmission. The probabilities of disease transmission from flies to dogs and from dogs to flies, denoted by b and c, respectively, have the indices of 0.47 each, while the index of sandfly birth rate, v has the index of 0.5 as shown in Table 2. So these three parameters need intervention. Alternatively, we address the biting rate, a, by introducing a control variable q 1 . This intervention consequently helps in controlling b, c, and v . Control variable q 1 is associated with insecticide-treated bed nets, sandfly repulsive lotion, and electric devices, etc. The sensitivity index of "elimination of seropositive/infected dogs" is -0.38. That is, an increase of 10% in the culling rate of infectious dogs will decrease R 0 by 3.8%. We introduce control variable q 2 for the control of seropositive dogs.
The mortality rate of flies has the sensitivity index of -1. This means that an increase in mortality of flies will cause a decrease in R 0 . To increase the death rate of flies, we introduce a control variable q 3 . As such, q 3 is associated with residual DDT sprays in animal and human shelters.
Based on the above control variables, we propose the following optimal control model: To reduce the density of sandflies, infected humans, and infected reservoir, we define an objective function (OF) as: where Q = g 1 I 1 + g 2 I 2 + g 3 P 2 + g 4 I 12 + g 5 P 12 + g 6 I 21 + g 7 I r .
Here N v denotes the total vector population and g i for i = 1, 2, 3, . . . , 8 denote the concerned weight constants; 1 ) denotes the cost associated with interventions.
Next we need to find the control function J 1 , subject to system (5), such that where D denotes the set of control variables and is defined as

Solution of the proposed optimal control (OC) model
Here we investigate existence of OC for the proposed model (5) at t = 0. Using [32], we claim that the solution of the state system (5) is bounded. This is because the control variables are bounded and Lebesgue measurable, and the initial conditions are nonnegative. Next we prove that the system has an optimal solution. For this, we introduce a Lagrangian and Hamiltonian as follows: where and Q = g 1 I 1 + g 2 I 2 + g 3 P 2 + g 4 I 12 + g 5 P 12 + g 6 I 21 + g 7 I r ; Theorem 1 The set (q * 1 , q * 2 , q * 3 ) of optimal controls minimizes J 1 over D, subject to the initial conditions specified at t = 0. Proof 1 Since both the state and control variables are nonnegative, using [33], the objective function is convex in the control variables q i , i = 1, 2, 3.
Clearly, from the definition, D is closed and convex. Further, since the optimal system is bounded, the set of optimal controls is compact.
And the control variables q 1 , q 2 , and q 3 enter as quadratic terms. Therefore the integrand in (6) is convex on the control set D.
The proof of existence of optimal control is completed.
Next we characterize the control variables of the proposed model. We use Pontryagin's maximum principle [34] stated here for convenience. Let be the states associated with control variables (q * 1 , q * 2 , q * 3 ). Let (y, q) be the solution of the control model/problem, for y ∈ and q = (q * 1 , q * 2 , q * 3 ), then there exists a vector function λ = (λ 1 , λ 2 , . . . , λ n ) satisfying the following: To establish the necessary condition on H, we prove the following result.

Numerical simulations and discussion
We solve the optimality system (14) numerically using RK4 method. The method solves the state system (5), forward in time, and the adjoint system (9), backward in the time, and system (10), the controls are updated continuously using systems (11) to (13). The process is continued until the results at the consecutive iterations (CI) are close. For details, see [35]. We assign the following values to weight constants (WC): g 1 = 0.5, g 2 = 0.5, g 3 = 0.5, g 4 = 0.5, g 5 = 0.5, g 6 = 0.5, g 7 = 0.5, g 8  Note Stage 1 infection means that a susceptible human at the first time is attacked by some strain. Stage 2 infection means that the individual has recovered form some strain of leishmaniasis and after recovery is attacked by another stain of leishmaniasis.
In the following figures we have presented the behavior of different subclasses of the total population involved in the model. Figure 2 represents two subclasses: S h , the susceptible human class, and E 1 , the C l -exposed human class. The figure shows that without intervention or control variables (presented by red lines), the density of the susceptible human class decreases because the susceptible humans are infected and they move from the susceptible class to the infected class E 1 . At the same time the density of the infected class increases. However, in case of proper intervention, the density of the exposed class decreases, and there is no notable decrease in the density of the susceptible human population. Figure 3 presents the behavior of E 2 , the V l -exposed human class, and I 1 , the infectious human class of C l strain. The gap between the red and blue lines represents the effectiveness of the control variables. Figure 4 presents the behavior of I 2 , the V l -infectious human class, and E 23 , the human class in the dormant period of developing PKDL. The figures show that the interventions have a notable effect on the control of I 2 and E 23 . Figure 5 represents the behaviors of P 2 , the human class infected with PKDL, and R 1 , the human class recovered from the cutaneous strain. Here too the effect of interventions is notable. The effect of interventions on R 2 is not so high. The reason might be the prolonged dormancy period of PKDL. All the above mentioned cases are concerned with Stage 1 infection and the effect of interventions is high.
From Fig. 6 to Fig. 9, these behaviors are associated with Stage 2 infection. Like in Fig. 6, E 12 represents the human class recovered from cutaneous leishmaniasis and infected/attacked by visceral leishmaniasis after the said recovery. In all Stage 2 infection cases, the effect of the interventions is not high. Figure 10 presents the behavior of S r , the susceptible reservoir, and I r , the infectious reservoir. The intervention involving "culling of seropositive dogs" directly addresses the infectious class of the reservoir population and hence the density of I r reduces rapidly. Similarly, the density of all the subclasses in the vector population reduces very rapidly because the use of residual sprays directly affects all the subclasses. Figure 11 represents the behavoir of susceptible sanflies and recovered resorvoirs. it takes 5-10 days to elimanate susceptible vectors using control variable. Figure 12 represents the behavoirs of exposed and infectious vectors with and without control. The density of exposed vectors reduces to zero in five days while that of infectious vectors reduces to zero in ten days. Figure 13 suggests that the interventions should be continue for a long time.

Conclusion
A model of three populations, namely sandflies, dogs, and humans, was considered for the epidemic of leishmaniasis. The main aim of the study was the minimization of the objective function, that is, the reduction of total sandflies population and the minimization of infected classes in the human and dog populations. Three control variables, q 1 , q 2 , and q 3 , were introduced to control the biting rate of sandflies and the densities of the infectious class of dog and vector populations. For this, first of all we investigated the sensitivity of the initial rate of transmission, R 0 , for different parameters involved in the transmission. Different parameters used in R 0 have different sensitivity indices. We needed to address the parameters with high sensitivity. However, some parameters are beyond the human control, like the birth and death rates of human and reservoir populations. The sensitivity index of a, the biting rate of sandfly, is 1. That is, a decrease of 20% in a would cause a decrease of 20% in R 0 , the initial rate of disease transmission. Similarly, μ v , v , μ r , and r have high sensitivity indices and hence a high impact on the disease transmission. The most sensitive parameters, namely a, the biting rate of sandfly, μ v , the mortality rate of sandflies, and μ r , the death rate of dogs, are directly influenced by using control variables q 1 , q 2 , and q 3 . Aso v , the birth rate of sandflies, is addressed indirectly by imposing restrictions on the collection of blood meal, and hence reducing ovigenesis. Numerical simulations show that as a result of these interventions cutaneous leishmaniasis can be eliminated in the period of 1200 days. However, I 2 , P 2 , and I 12 reduce to zero in 750, 1750, and 1150 days, respectively. Moreover, P 12 , I 21 , I r , and I v reduce to zero in 1700, 950, 7, and 7 days, respectively.