Electronic Journal of Differential Equations, Vol. 2020 (2020), No. 80, pp. 1–24. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu OPTIMAL CONTROL APPLIED TO A VISCERAL LEISHMANIASIS MODEL BUDDHI PANTHA, FOLASHADE B. AGUSTO, IBRAHIM M. ELMOJTABA Communicated by Suzanne Lenhart Abstract. In this article, we developed a deterministic model for the trans- mission dynamics of visceral leishmaniasis in humans, canine reservoirs and sandflies, which is the only vector that transmits the disease parasite. The theoretical and epidemiological findings of this study indicates that the disease- free equilibrium of the model is locally and globally asymptotically stable when the associated reproduction number is less than unity. We perform sensitivity analysis on the model parameter to determine the parameter with the most impact on the reproduction number. Following the results obtained from the sensitivity analysis, we apply optimal control theory using three time depen- dent control variables representing personal protection, insecticide spraying and culling of infected canine reservoirs. Simulation results are presented for various outbreak scenarios which indicates that leishmaniasis can be eliminated from a region by the application of three time dependent controls represent- ing respectively, personal protection, insecticide spraying and culling infected canine reservoir. 1. Introduction Leishmaniasis is a vector borne parasitic disease caused by a protozoan parasite called Leishmania which is transmitted by a bite of an infected female phlebotomine sandflies [7, 46]. These sandflies species feed on animal blood which is required for egg development [45]. When a sandfly bites an infected human or animal, it picks up the parasites and can transmit them into a new uninfected host while feeding on (or biting) the host [13]. During the feeding, the infected sandfly produces saliva and parasitic proteins that interact with the skin of the vertebrate host, crucial for establishing Leishmania in the skin of a vertebrate [7, 46]. The species of the Leishmania parasite determines the symptoms of the disease such as cutaneous, mucocutaneous and visceral [7, 46]. Like many other disease, leishmanasis is more prevalent in the poorest communities in the developing world. The actual number of annual new cases of leishmaniasis is not known with certainty [7], however, the World Health Organization (WHO) estimate is about 350 million people at risk of contracting leishmanasis and occurrence of about 2 million new cases every year [46]. 2010 Mathematics Subject Classification. 35F21, 37N25. Key words and phrases. Visceral leishmanisis; PKDL; vaccination; canine reservoir; optimal control. c©2020 Texas State University. Submitted January 6, 2019. Published July 28, 2020. 1 2 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 There are over 30 sandflies species which cause human and animal leishmania- sis. The clinical symptoms and disease progression depend on Leishmania species, host immune system and factors such as environmental conditions, seasons, socio- economic status [7, 46]. Three major forms of human leishmaniasis infection are reported: cutaneous, mucocutaneous and visceral. The cutaneous leishmaniasis (CL), which causes skin sores and lesions, is most common form of leishmaniasis and is considered least fatal [7, 8, 46]. Most skin lesions due to CL infection heal without medical intervention, however, some infection spread to several parts of the body such as eyebrows, earlobes, limbs and throats. These CL infections may not heal by their own and may require medical attention [7, 46]. In the past 5 years, about 1 million CL cases were reported worldwide, with the majority oc- curring in Americas, the Mediterranean basin, the Middle East and Central Asia [7, 46]. Mucocutaneous leishmaniasis (MCL), causes partial or total destruction of mucous membranes of the nose, mouth and throat and usually occurs in patients infected with visceral leishmaniasis or post kala-azar dermal leishmaniasis (PKDL) [7, 46]. The patients co-infected with leishmaniasis and other disease such as HIV may also have MCL infection [7, 46]. MCL does not heals all spontaneously but it can be cured with timely treatment [7, 46]. Almost 90% of MCL cases occurs in the Plurinational State of Bolivia, Brazil and Peru [46]. Visceral leishmaniasis (VL), also known as kala-azar, is the most severe form of all leishmaniasis infections. It is fatal in almost all cases if left untreated [7, 13, 46]. Caused by either of L. donovani or L. infantum [39, 46] parasites, it can infect people of any ages but children under 10 years are more susceptible [7, 46]. The incubation period of visceral leishmanasis ranges from 10 days to several months and the onset of the symptoms is usually gradual [7, 46]. VL makes up an estimated 0.2 million to 0.4 million cases annually with over 0.2 million deaths. Over 90% new cases of VL occur in Bangladesh, Brazil, Ethiopia, India, South Sudan and Sudan [7, 46] annually. Most VL infections are asymptomatic but some victims eventually develop clinical visceral leishmaniasis. Following recovery from visceral leishmania, some individuals develop post-kala-azar dermal leishmaniasis (PKDL) which appears as macular, papular or nodular rash on the face, upper arms, trunks and other parts of the body [39, 46]. PKDL usually appears 6 months to 1 or more years after successful treatments of visceral leishmaniasis. In an endemic region, Leishmania species are sustained by the interactions be- tween sandflies species and the reservoir host species. The reservoirs can be zoonotic or anthroponotic. The zoonotic reservoir hosts may be wild or domestic animals, and the anthroponotic reservoir hosts are humans [12, 36, 46]. Some zoonotic reser- voirs may or may not show symptoms while others such as dogs may eventually die due to the disease. Such animals are perfect environment for Leishmania parasite spread since the parasites spread in their blood and on the dermis where they are bitten by sandflies [7, 46]. In an endemic area, an animal species is usually a prin- ciple reservoir host for a particular Leishmania species but other animal species in the same area may be an incidental reservoirs [7, 12, 36, 46]. In some VL endemic areas, seroprevalence test suggested that goats are reservoir host for VL [12, 36], however, many studies [12, 36, 46, 32] found that dogs are the principal reservoir host of Leishmania infantum, a major causative agent of visceral leishmaniasis, and more than 50% of all infected dogs are asymptotic carriers. EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 3 Because of technological advancement, many technologies for disease diagnosis, prevention and treatment for a wide range of diseases have been developed. How- ever, population growth, movement, environmental changes and urbanization have created favorable epidemiological conditions for disease transmission and persis- tence of some diseases such as leishmaniasis [12]. In fact, because of the rarity of functioning control programs, the new cases of leishmaniasis is in the increasing trend [46]. Since the disease transmission is maintained in a complex biological system involving human and animal reservoir, and vectors, the control measures should target interaction between these components [13, 39, 40, 41, 46]. Current prevention strategies include indoor residual spray and using door and window screens or insect repellent. These measures are only effective in reducing indoor interactions between human and sand flies but have no effect for the species of sand flies who feed on animals outdoors [46]. Also, it has been reported that the use of insect repellent can decrease sand flies and human interactions in a particular house but their biting rate increased in the nearby unprotected human and animals [7, 46]. In some cases, culling of seropositive reservoir such as dogs has been im- plemented to reduce this interaction. Some vaccine have recently been introduced but its efficacy in human is yet to be determined [12, 23, 24]. Canine vaccine has shown promising decrease in infection in some countries such as Brazil, France and Iran but no vaccines are available for humans [12, 46]. Numerous mathematical studies have been carried out on leishamaniasis. Chaves et al. [8] developed a parasite-reservoir-incidental host model consisting a system of ordinary differential equation (ODE) to study the threshold conditions for the persistence of the infection for American cutaneous leishmaniasis (ACL). Ribas et al. [34] developed an ODE model that includes constant control strategies (treat- ment, vaccination, culling and insecticide collar) applied to dogs. Their study found that the strategy of culling infected dogs is not the most efficient way to control zoonotic visceral leishmaniasis (ZVL) from both efficiency and ethical point of views. Other methods such as vector control and use of insecticide impregnated dog collars are more efficient in reducing ZVL in humans. Their study also noted that treating infected dogs is not effective to reducing the infections in human. Palatnik-de-Sousa et al. [32] developed an ODE model to study the efficacy of the culling of seropositive dogs. Their findings suggest that removing about 25% of the seropositive dogs from the endemic region significantly decreases the possibility of an epidemic outbreak in the region [32]. Some other studies such as Elmojtaba et al. [13], Stauch et al. [39, 40], Subramanian et al. [41] also explored the combina- tion of human drug treatment, vector and reservoir hosts control for the effectively curtailing the spread of VL in both human and reservoir populations. However, none of these studies incorporated the use of optimal control theory to investigate the impact of these intervention strategies and the cost of their application. In our model, we assume three different host populations: human host, sandflies (vector) and canine host; also, canine are assumed to be source of infection, i.e. sandflies may acquire the infection from dogs, as suggested by many studies; see for exam- ple [32]. Susceptible sandflies get infected after biting infected human or infected canine reservoir. Infected sandflies transmit the disease following their bites on susceptible human or canine hosts. We assume four compartments on human host populations, 2 compartments on vector population and 3 compartments on canine 4 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 reservoir populations. The detailed assumptions of our model are presented in Section 2. In this study, we aim to investigate the effect of optimal intervention strategies for visceral leishmanasis and the cost of applying these controls. In Section 2, we develop and analyze the dynamics of the mathematical model involving human- vector-reservoir populations and studied the stability of disease free equilibrium. In Section 3, we carry out the sensitivity analysis to identify parameters with the most impact in the disease transmission. We formulate the optimal control problem is Section 4, characterize the optimal controls, and present the numerical results in Section 5. Section 6 presents conclusions and discussions of our study. 2. Model formulation To formulate our model we consider three different populations: human host pop- ulation, NH(t), canine host population, NR(t), and vector population, NS(t). The human host population is divided into four sub-populations: susceptible individuals SH(t), visceral leishmaniasis infected individuals IH(t), those who develop PKDL after the treatment of visceral leishmaniasis, PH(t), and those who are recovered and have permanent immunity, RH(t). This implies that NH(t) = SH(t) + IH(t) + PH(t) +RH(t). Similarly, the canine host population is divided into three sub-populations: suscep- tible, SR(t), vaccinated, VR(t) and infected, IR(t), such that NR(t) = SR(t) + VR(t) + IR(t), and finally, the sandfly population have two sub-populations: susceptible, SS(t), and infected sandflies, IS(t), such that NS(t) = SS(t) + IS(t). It is assumed that susceptible individuals are recruited into the population at a constant rate ΛH , where recruitment is mainly by births. Also, assume that the susceptible animals acquire infection with leishmaniasis following contacts with infected sandflies at a per capita rate βHbS IS NH , where βH is the transmission prob- ability per bite per human (as the case for malaria, [28, 35]) and bS is the per capita biting rate of sandflies on humans. Infected humans die due to leishmaniasis at an average rate δH , or get treatment at an average rate γH , and a fraction α of them recover and acquire permanent immunity, and the complement fraction (1 − α) develop PKDL. Humans with PKDL get treated at an average rate τH , or recover naturally at an average rate σH , and acquire permanent immunity in both cases. There is a per capita natural mortality rate µH in all human sub-population. Susceptible sandflies are recruited at a constant rate ΛS , and acquire leishma- niasis infection following contacts with a leishmaniasis infected human or with a human having PKDL or leishmaniasis infected canine at an average rate equal to βSbS IH NH + βSbS PH NH + βSbSR IR NR , where bS is the per capita biting rate on human and bSR is the per capita biting rate on canine, and βS is the transmission proba- bility for sandfly infection after biting a human or a canine. Sandflies suffer natural mortality at a per capita rate µS regardless of their infection status. Susceptible canine reservoirs are recruited into the population at a constant rate ΛR, acquire infection with leishmaniasis following contacts with infected sandflies at a rate βRbSR IS NR where βR and bSR as described above; or vaccinated at an EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 5 average constant rate νR. Vaccinated canine reservoir acquire infection with leish- maniasis following contacts with infected sandflies at a rate βRbSR(1− ε) ISNR where ε is a modification parameter measures the efficacy of the vaccination, and it is assumed that the vaccine does not wane. Infected canine reservoirs die because of leishmaniasis at an average rate δR. Using the above description, the following mathematical model is proposed S′H = ΛH − βHbSISSH NH − µHSH I ′H = βHbSISSH NH − (γH + µH + δH)IH P ′H = (1− α)γHIH − (τH + σH + µH)PH R′H = αγHIH + (τH + σH)PH − µHRH S′S = ΛS − βSbS(IH + PH)SS NH − βSbSRIRSS NR − µSSS I ′S = βSbS(IH + PH)SS NH + βSbSRIRSS NR − µSIS S′R = ΛR − βRbSRISSR NR − (νR + µR)SR V ′R = νRSR − βRbSR(1− ε)ISVR NR − µRVR I ′R = βRbSRISSR NR + βRbSR(1− ε)ISVR NR − (µR + δR)IR (2.1) with initial conditions SH(0) = S0 H , IH(0) = I0 H , PH(0) = P 0 H , RH(0) = R0 H , SV (0) = S0 V , IV (0) = I0 V , SR(0) = S0 R, IR(0) = I0 V . and N ′H = ΛH − µHNH − δHIH , N ′S = ΛS − µSNS , N ′R = ΛR − µRNR − δRIR. Invariant region. All parameters of the model are assumed to be constant and nonnegative and all state variables are non-negative at time t = 0. Also, note that in the absence of disease induced death (i.e. δH = δR = 0), the total human population, NH → ΛH/µH as t → ∞. Similarly, NR → ΛR/µR and NS → ΛS/µS as as t→∞. This shows that the biologically-feasible region: Ω = { (SH , IH , PH , RH , SS , IS , SR, VR, IR) ∈ R9 + : SH , IH , PH , RH , SS , IS , SR, VR, IR ≥ 0, NH ≤ ΛH µH , NS ≤ ΛS µS , NR ≤ ΛR µr } is positively-invariant domain, and thus, the model is epidemiologically and math- ematically well posed, and it is sufficient to consider the dynamics of the flow generated by the system (2.1) in this positively-invariant domain Ω. The flow diagram of the system (2.1) is depicted in Figure 1 and the associated variables and parameters are described in Table 1. 6 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 IHSH RH PH βHbS IS NH αγ H (1 − α)γ H τHσH SS IS βSbS (IH+PH ) NH βSbSR IR NR SR IR VR νR βRbSR I S N R βRbSR(1− ε) IS NR ΛH µH µH δH µH µH ΛS µS µS ΛR µR µR µR δR Figure 1. Systematic flow diagram of model (2.1). EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 7 Table 1. Description of the variables and parameters for model (2.1) . Variable Description SH Population of susceptible humans IH Population of symptomatic humans PH Population of symptomatic humans RH Population of recovered humans SS Population of susceptible sandflies IS Population of infected sandflies SR Population of susceptible reservoir (canine) VR Population of vaccinated reservoir (canine) IR Population of infected reservoir (canine) Param. Description ΛH Recruitment rate of humans ΛS Recruitment rate of sandflies ΛR Recruitment rate of reservoir (canine) βH Transmission probability per contact for susceptible humans βS Transmission probability per contact for susceptible sandflies βR Transmission probability per contact for susceptible reservoir (canine) bS Sandflies biting rate in humans bSR Sandflies biting rate in reservoir (canine) µH Natural death rate of humans µS Natural death rate of sandflies µR Natural death rate of reservoir (canine) δH Disease induced death rate of humans δR Disease induced death rate of symptomatic reservoir (canine) γH Recovery rate of infected humans with VL dues to treatment α Fraction successfully treated for VL τH Recovery rate of humans with PKDL due to treatment σH Instantaneous recovery rate of humans with PKDL νR Vaccination rate in reservoir (canine) ε Vaccination efficacy in reservoir (canine) 2.1. Analysis of the Model. The disease-free equilibrium (DFE) of the system (2.1) is given by E0 = ( S∗H , I ∗ H , P ∗ H , R ∗ H , S ∗ S , I ∗ S , S ∗ R, V ∗ R, I ∗ R ) = (ΛH µH , 0, 0, 0, ΛS µS , 0, ΛR µR + νR , νRΛR µR(µR + νR) , 0 ) . To study the stability of the disease-free equilibrium, first we have to find the reproduction number which is defined as the number of secondary infections that occur when an infected individual is introduced into a completely susceptible pop- ulation [11, 16]. To calculate the effective reproduction number, we will use the next generation approach [11, 42]. First, we take the variables IH , PH , IS , IR as the infected compartments and then use the notation in [42], the Jacobian F and V matrices for new infectious terms and the remaining transfer terms, respectively, 8 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 are defined as: F =  0 0 βHbS 0 0 0 0 0 βSbSΛSµH µSΛH βSbSΛSµH µSΛH 0 βSbSRΛSµR µSΛR 0 0 βRbSR[µR+(1−ε)νR] µR+νR 0  V =  k1 0 0 0 −(1− α)γH k2 0 0 0 0 µS 0 0 0 0 k3  where k1 = γH+δH+µH , k2 = τH+σH+µH , k3 = µR+δR. Then the reproduction number is R0 = ρ(FV −1) = √ (RH +RR)RS where, ρ is the spectral radius and RH = b2SµHβH [(1− α)γH + k2)] k1k2ΛH , RR = b2SRβRµR[µR + (1− ε)νR] k3ΛR(µR + νR) , RS = βSΛS µ2 S . Furthermore, the quantityRH is the number of secondary infections in humans host by one infectious sandfly, RR is the number of secondary infections in the canine reservoir host by one introduced infectious sandfly, and lastly RS is the number of secondary infections in sandflies resulting from a newly introduced infectious human and canine reservoirs respectively. Using van den Driessche and Watmough [42, Theorem 2], the following result is established: Lemma 2.1. The disease-free equilibrium (E0) is locally asymptotically stable if R0 < 1, and unstable if R0 > 1. To study the global behavior of system (2.1), we use a theorem by Castillo- Chavez et al. [6], stated here for convenience. Theorem 2.2 ([6]). For the system dX dt = F (X,Z) dZ dt = G(X,Z) G(X, 0) = 0 , where the components of the column-vector X ∈ Rm denote the number of un- infected individuals and the components of vector Z ∈ Rn denote the number of infected individuals. U0 = (X∗, 0) denotes the disease-free equilibrium of this sys- tem. The fixed point U0 = (X∗, 0) is a globally asymptotically stable equilibrium for this system provided that R0 < 1 (locally asymptotically stable) and the following two conditions satisfied: (H1) For dX dt = F (X, 0), X∗ is globally asymptotically stable, EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 9 (H2) G(X,Z) = AZ − Ĝ(X,Z), Ĝ(X,Z) ≥ 0 for (X,Z) ∈ Ω, where A = DZG(X∗, 0) is an M-matrix (off-diagonal elements of A are non-negative) and Ω is the region where the model has biological meaning. We can rewrite our system (2.1) using above notation, where X = (SH , RH , SS , SR, VR), Z = (IH , PH , IS , IV ) (2.2) F (X, 0) = ( ΛH − µHSH 0 ΛS − µSSS ΛR − (νR + µR)SR νRSR − µRVR )T A =  −k1 0 βHbs 0 (1− α)γH −k2 0 0 βSbS SS NH βSbS SS NH −µS βSbSR SS NR 0 0 βSbSR SR+(1−ε)VR NR −k3  Ĝ(X,Z) =  βHbSIS ( 1− SH NH ) 0 βSbSSS ( 1− IH+PH NH ) + βSbSRSH ( 1− IR NR ) βRbSR ( 1− SR+(1−ε)VR NR )  it is clear that A is an M matrix and Ĝ(X,Z) > 0, and hence we have the following lemma. Lemma 2.3. The disease-free equilibrium is globally asymptotically stable when R0 < 1. 3. Sensitivity analysis The impact of leishmaniasis model (2.1) significant parameters are determined via sensitivity analysis [3, 29, 31] on the model outcome using Latin hyper-cubic sampling (LHS) technique and partial rank correlation coefficient (PRCC). LHS is a stratified sampling without replacement method which allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [3, 29, 30, 37]. PRCC measures the strength of the relationship between the model outcome and the parameters, stating the degree of the effect that each parameter has on the outcome [3, 29, 30, 37]. A total of 1,000 simulations of the model (2.1) per LHS run were carried out, using the ranges and baseline values tabulated in Table 2 (with the reproduction number, R0, as the response function). The values of the parameters are taken from published literature as mentioned in Table 2 and for each parameter value, an interval within 20% range of the parameter value is formed to test the sensitivity. Figure 2 depicts the PRCC values for each parameter of the models using the reproduction number, R0, as the response function. Parameters with the highest PRCC values have the largest impact onR0. Bars extending to the left (for negative PRCC values) or to the right (for positive PRCC values). The parameters γH , α, σH , νR, all have every low PRCC values with p-values > 0.01 (i.e., p-values 0.4390, 0.6311, 0.0215, 0.6777 respectively) indicating they have the least influence on the reproduction number, R0. Therefore, the key parameters influencing R0 are separated into those that decrease R0 when increased (those with significantly negative PRCC values) and those that causes R0 to increase when increased (those with significantly positive PRCC values). 10 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 H S R H S R b S b SR H S R H R H H H R Figure 2. PRCC values for the leishmaniasis model (2.1), using the reproduction number (R0) as the response function. Parameter values (baseline) and ranges used are as given in Table 2. The identification of these key parameters with significant impact on the re- production number R0 is vital to the formulation of effective control strategies necessary to combat the spread of the disease. The results from the sensitivity analysis therefore suggest that, to effectively curtail the spread of leishmaniasis transmission in a region, the control strategy to be implemented should decrease the transmission probability per contact for susceptible canine population (reduce βR), this should be followed by decreasing the transmission probability for sus- ceptible sandflies (reduce βS), the recruitment rate of sandflies (reduce ΛS), the sandflies biting rate on humans (reduce bS), the sandflies biting rate on the canine population (reduce bSR), and the transmission probability for susceptible humans (reduce βH). The result from the sensitivity analysis further suggests control strategy that increases the natural death rate of sand flies (increase µS), the vaccination efficacy in canine reservoir (increase ε), and the disease-induced death rate of infected canine population (reduce δR) will help to reduce the reproduction number and thereby reduce the spread of leishmaniasis transmission among humans and the canine population. Thus, in summary, the sensitivity analysis of the leishmaniasis transmission model (2.1) shows that the significant parameters are the natural death rate of sandflies (µS), sandflies biting rate on the canine population (bSR), transmission probability for susceptible sandflies (βS), recruitment rate of sandflies (ΛS), sand- flies biting rate in humans (bS), transmission probability for susceptible canine population (βR), transmission probability for susceptible humans (βH), vaccina- tion efficacy in canine population (ε), recruitment rate of canine population (ΛR), and the disease-induced death rate of infected canine population (δR). EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 11 Table 2. Parameters values of the leishmaniasis model (2.1). Param. Values Range References ΛH 4.16× 10−5 (1− 0.2)4.16× 10−5 - (1 + 0.2)4.16× 10−5 [44] ΛS 0.155 0.124 - 0.186 [20] ΛR 0.0027 0.00216 - 0.00324 assumed βH 0.5 0.4 - 0.6 [40] βS 0.22 0.176 - 0.264 [13] βR 0.0714 0.05712 - 0.08568 [13] bS 0.2856 0.22848 - 0.34272 [13] bSR 0.56 0.448 - 0.672 assumed µH 1.64× 10−5 1.312× 10−5 - 2.952× 10−5 [44] µS 0.056 0.0448 - 0.0672 [20] µR 2.11× 10−4 1.688× 10−4 - 2.532× 10−4 assumed δH 0.0014 0.00112 - 0.00168 [38] δR 0.0014 0.00112 - 0.00168 [2] γH 0.5 0.4 - 0.6 assumed α 0.40 0.32 - 0.48 [15] τH 0.033 0.02643 - 0.0396 [13] σH 0.00556 0.00448 - 0.006672 [13] νR 0.5 0.4 - 0.6 assumed ε 0.8 0.64 - 0.96 [12] 4. Optimal control problem Following the results obtained from the sensitivity analysis, we introduce three time dependent controls into the system (2.1). Let u1(t) be the time dependent rate of use of personal protection such as insect repellant, door and window screens against sand flies in the bid to reduce human contact with the flies. Let u2(t) be the time dependent indoor insecticide spraying rate (insecticide such as pyrethroids are effective against sandflies both indoors and outdoors [1]). Lastly, let u3(t) be the time dependent culling rate of the infected canine population [36]. Then the above system of ODEs (2.1) becomes S′H = ΛH − βHbS [1− u1(t)]ISSH NH − µHSH I ′H = βHbS [1− u1(t)]ISSH NH − (γH + δH + µH)IH P ′H = (1− α)γHIH − (τH + σH + µH)PH R′H = αγHIH + (τH + σH)PH − µHRH S′S = ΛS − βSbS [1− u1(t)](IH + PH)SS NH − βSbSRIRSS NR − [µS + u2(t)]SS I ′S = βSbS [1− u1(t)](IH + PH)SS NH + βSbSRIRSS NR − [µS + u2(t)]IS S′R = ΛR − βRbSRISSR NR − [µR + νR]SR 12 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 V ′R = νRSR − βRbSR(1− ε)ISVR NR − µRVR I ′R = βRbSRISSR NR + βRbSR(1− ε)ISVR NR − [µR + δR + u3(t)]IR. (4.1) We want to find the controls that minimizes the total infected humans (with VL and PKDL), infected canine population and the cost of controls. In other words, we want to find the optimal values of (u∗1, u ∗ 2, u ∗ 3) that minimizes the objective functional J(u1, u2, u3) where J(u1, u2, u3) = ∫ T 0 [A1IH +A2PH +A3IR +B1u2SS +B2u2IS +B3u3IR + C1u 2 1 + C2u 2 2 + C3u 2 3]dt. (4.2) subject to the differential equations (4.1), where T is the final time. This per- formance specification involves the total infected humans (with VL and PKDL), infected reservoir (canine), along with the cost of applying the controls u1(t), u2(t) and u3(t)). In this paper, a quadratic objective functional is implemented for mea- suring the control cost, such a cost has been frequently used [17, 18, 21, 22, 25, 47]. The positive coefficients, Ai, Bi, Ci; i = 1, · · · , 3 are balancing weight parameters. The controls, u1(t), u2(t) and u3(t), are bounded, Lebesgue integrable functions [18, 48]. And we seek to find optimal controls, u∗1, u ∗ 2 and u∗3, such that J(u∗1, u ∗ 2, u ∗ 3) = min (u1,u2,u3)∈U {J(u1, u2, u3)} (4.3) where the admissible set is U = {(u1, u2, u3) ∈ (L∞(0, T ))3 : 0 ≤ ui ≤Mi;Mi ∈ R+, i = 1, 2, 3}. The following theorem proves the existence of the solution of the system (4.1) as well as the non negativity and boundedness of the state variables. Theorem 4.1. Given controls u = (u1, u2, u3) ∈ U , there exist non-negative bounded solutions (SH , IH , PH , RH , SS , IS , SR, IR) to the state system (4.1) in the finite interval [0, T ] with given initial conditions. The existence and uniqueness of solutions for the state system (4.1) with a given control pair can be proven by using a result from Lukes [26]. The structure of system (4.1) gives the non-negativity and uniform boundedness of the state solutions. The next theorem proves the existence of the optimal controls. Theorem 4.2. There exists an optimal control tuple u∗ = (u∗1, u ∗ 2, u ∗ 3) ∈ U with cor- responding states (S∗H , I ∗ H , P ∗ H , R ∗ H , S ∗ S , I ∗ S , S ∗ R, V ∗ R, I ∗ R) that minimizes the objective functional J(u1, u2, u3). Proof. Since the controls and the state variables are uniformly bounded and non- negative on the finite interval [0, T ], there exists a minimizing sequence (un1 , u n 2 , u n 3 ) such that lim n→∞ J(un1 , u n 2 , u n 3 ) = inf (u1,u2,u3)∈U J(u1, u2, u3). Let us denote (SnH , I n H , P n H , R n H , S n S , I n S , S n R, V n R , I n R) = (SH , IH , PH , RH , SS , IS , SR, VR, IR)(uni , u n 2 , u n 3 ). EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 13 Since all of the state variables are bounded (Theorem 4.1), then their first deriva- tives are also bounded. Also, the control functions are assumed to be bounded. This implies that all state variables are Lipschitz continuous with the same Lipschitz con- stant. Thus the sequence (SnH , I n H , P n H , R n H , S n S , I n S , S n R, V n R , I n R) is uniformly equicon- tinuous in [0, T ]. Then by the Arzela-Ascoli Theorem [26], the state sequence has a subsequence that converges uniformly to (S∗H , I ∗ H , P ∗ H , R ∗ H , S ∗ S , I ∗ S , S ∗ R, V ∗ R, I ∗ R) in [0, T ]. The control sequence (un1 , u n 2 , u n 3 ) has a subsequence that converges weakly in L2(0, T ). Let (u∗1, u ∗ 2, u ∗ 3) ∈ U be such that uni ⇀ u∗i weakly in L2(0, T ) for i = 1, 2, 3. This result together with the uniform convergence of the state system implies the convergence of the terms like B1u n 2S n S . Using lower semi-continuity of norms in weak L2, we obtain ‖u∗i ‖2L2 ≤ lim inf n→∞ ‖uni ‖2L2 for i = 1, 2, 3. Hence, J(u∗1, u ∗ 2, u ∗ 3) ≤ lim inf n→∞ ∫ T 0 [ A1I n H +A2P n H +A3I n R +B1u n 2S n S +B2u n 2 I n S +B3u n 3 I n R + C1(un1 )2 + C2(un2 )2 + C3(un3 )2 ] dt = lim inf n→∞ J(un1 , u n 2 , u n 3 ) Thus, there exists a vector ~u = (u∗1, u ∗ 2, u ∗ 3) of controls that minimizes the objective functional J(u1, u2, u3). � Next, we characterize the optimal control triple (u1(t), u2(t), u3(t)). Pontryagin’s Maximum Principle [33] introduces adjoint functions that allow the state system (4.1) to be attached to the objective functional (4.2). Characterization of optimal controls. The necessary conditions that an op- timal control must satisfy come from the Pontryagin’s Maximum Principle [33]. This principle converts (4.1) and (4.2) into a problem of minimizing pointwise a Hamiltonian H, with respect to the controls u1, u2, u3 ∈ U . First we formulate the Hamiltonian [25] from the cost functional (4.2) and the governing dynamics (4.1) to obtain the optimality conditions. H = A1IH +A2PH +A3IR +B1u2SS +B2u2IS +B3u3IR + C1u 2 1 + C2u 2 2 + C3u 2 3 + λSH [ ΛH − βHbS [1− u1(t)]ISSH NH − µHSH ] + λIH [βHbS [1− u1(t)]ISSH NH − (γH + δH + µH)IH ] + λPH [ (1− α)γHIH − (τH + σH + µH)PH ] + λRH [ αγHIH + (τH + σH)PH − µHRH ] + λSS [ ΛS − βSbS [1− u1(t)](IH + PH)SS NH − βSbSRIRSS NR − [µS + u2(t)]SS ] + λIS [βSbS [1− u1(t)](IH + PH)SS NH + βSbSRIRSS NR − [µS + u2(t)]IS ] + λSR [ ΛR − βRbSRISSR NR − [µR + νR]SR ] 14 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 + λVR [ νRSR − βRbSR(1− ε)ISVR NR − µRVR ] + λIR [βRbSRISSR NR + βRbSR(1− ε)ISVR NR − [µR + δR + u3(t)]IR ] . Using the system of adjoint variables, λ′SH = − ∂H ∂SH = (λSH − λIH )βHbS [1− u1(t)]IS (IH + PH +RH) N2 H + µHλSH + (λIS − λSS )βSbS [1− u1(t)]SS (IH + PH) N2 H , (4.4) λ′IH = − ∂H ∂IH = −A1 + (λIH − λSH ) βHbS [1− u1(t)]ISSH N2 H − λPH (1− α)γH − λRH αγH + λIH (γH + δH + µH) + (λSS − λIS )βSbS [1− u1(t)]SS [SH +RH ] N2 H , λ′PH = − ∂H ∂PH = −A2 + (λIH − λSH ) βHbS [1− u1(t)]ISSH N2 H + λPH (τH + σH + µH) − λRH (τH + σH) + (λSS − λIS )βSbS [1− u1(t)]SS [SH +RH ] N2 H , λ′RH = − ∂H ∂RH = (λIH − λSH )βHbS [1− u1(t)]IS SH N2 H + λRH µH + (λIS − λSS )βHbS [1− u1(t)]SS (IH + PH) N2 H , λ′SS = − ∂H ∂SS = (λSS − λIS ) [ βSbS [1− u1(t)] (IH + PH) NH + βSbRIR NR ] −B1u2(t) + λSS [µS + u2(t)], λ′IS = − ∂H ∂IS = −B2u2(t) + (λSH − λIH ) βHbS [1− u1(t)]SH NH + λIS [µS + u2(t)] + (λSR − λIR) βRbSRSR NR + (λVR − λIR) βRbSR(1− ε)VR NR , EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 15 λ′SR = − ∂H ∂SR = (λIS − λSS ) βSbSRIRSS N2 R + (λSR − λIR) βRbSRIS(IR + VR) N2 R + λIR(µR + νR) − λVR νR + (λIR − λVR ) βRbSR(1− ε)ISVR N2 R , λ′VR = − ∂H ∂VR = (λIS − λSS ) βSbSRIRSS N2 R + (λIR − λSR ) βRbSRISSR N2 R + (λVR − λIR)βRbSR(1− ε)IS (SR + IR) N2 R + λVR µR, λ′IR = − ∂H ∂IR = −A3 −B3u3(t) + (λSS − λIS ) βSbSRSS(SR + VR) N2 R + (λIR − λSR ) βRbSRISSR N2 R + (λIR − λVR ) βRbSR(1− ε)ISVR N2 R + λIR(µR + δR + u3), with the transversality condition λi(tf ) = 0 for i = SH , IH , PH , RH , SS , IS , SR, VR, IR. (4.5) The optimality conditions are given as ∂H ∂uj = 0, j = 1, 2, 3. in the interior of the control set U . Thus, we have 2C1u1 + λ∗SH βHbSI ∗ SS ∗ H N∗H − λ∗IH βHbSI ∗ SS ∗ H N∗H + λ∗SS βSbS(I∗H + P ∗H)S∗S N∗H − λ∗IS βsbS(I∗H + P ∗H)S∗S N∗H = 0, B1S ∗ S +B2I ∗ S + 2C2u2 − λ∗SS S∗S − λ∗ISI ∗ S = 0, B3I ∗ R + 2C3u3 − λ∗IRI ∗ R = 0. Hence, the control characterization is u∗1 = min [ 1,max ( 1 2C1N∗H [βHbSI ∗ SS ∗ H(λ∗IH − λ ∗ SH ) + βSbSS ∗ S(I∗H + P ∗H)(λ∗IS − λ ∗ SS )], 0 )] , u∗2 = min [ 1,max ( 1 2C2 [−B1S ∗ S −B2I ∗ S + λ∗SS SS + λ∗ISI ∗ S ], 0 )] , 16 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 u∗3 = min [ 1,max ( 1 2C3 [−B3I ∗ R + λ∗IRI ∗ R], 0 )] . Setting bounds on the control. We set the bounds on the time dependent control, u1(t), for the personal protection rate between 0 and 1 (i.e, 0 ≤ u1(t) ≤ 1). The value of u1(t) = 1 means highly effective control with no contact between human and sand flies while the value of u1(t) close to 0 means high contact between human and sand-flies. For the the control u2(t), which is the rate of insecticide spray, the lower bound is assumed to be 0 which means no spraying. For the upper bound, we assume that in an ideal condition, at most 75% of whiteflies can be killed within 10 days. With this assumption, the upper bound of the spray rate u2(t) is 0.13 i.e, 0 ≤ u2(t) ≤ 0.13. Since the insecticide spray kills the sandflies that come in contact with the spray, we think killing coverage of 75% of the sandflies in 10 days is reasonable. Finally, for the infected dog culling rate u3(t), the lower bound for u3(t) is zero for no culling. For the upper bound, we assume that in an ideal condition, at most 50% of the infected dogs can be killed within a week. With this assumption, the upper bound of u3(t) is 0.1 i.e, 0 ≤ u3(t) ≤ 0.1. Next, we discuss the numerical solutions of the optimality system, the corre- sponding optimal control and the interpretations from various cases. 5. Numerical results The following algorithm was used to compute the optimal controls and state values using a Runge Kutta method of the fourth order in the interval [0,180] days i.e.(6 months). First, an initial estimate for the controls are made. Then the state variables are solved forward in time using the dynamics (4.1). The results obtained for the state variables are substituted into the adjoint equations (4.4). These adjoint equations with given final conditions (4.5) are then solved backwards in time, employing the backward fourth order Runge Kutta method. Both the state and adjoint values are then used to update the controls, and the process is repeated until the current state, adjoint, and controls values converge sufficiently [25]. To illustrate the optimal control strategies, the following initial conditions were used: SH(0) = 1000, IH(0) = 25, PH(0) = 25, RH(0) = 14, SS(0) = 500, IS(0) = 152, SR(0) = 224, VR(0) = 125, IR(0) = 114. Furthermore, the values A1 = A2 = A3 = 1.00, B1 = B2 = B3 = 0.01, C1 = C2 = C3 = 0.01 were chosen as the baseline weight parameters. It should be noted that the weights in the simulations here are only of theoretical sense to illustrate the control strategies proposed in this paper. Using parameter values in Table 2, the reproduction number is obtained as R0 = 2.33 > 1, thus indicating that the disease is endemic in the population. The results of the optimal control simulations of the leishmaniasis model (4.1) are depicted in Figures 3, where the total number of infectious humans, sandflies and canine reservoir in the absence of controls are denoted by blue curves and the corresponding numbers with optimal controls are denoted by pink curve. From the figure 3(a), we observed that the total number of infected humans (IH and PH) are reduced considerably with the applications of the optimal controls compared to those in the absence of controls. The total number of infected canine reservoirs and sandflies are shown in Figures 3(b) and 3(c). There are more infected in the absence of controls compared to the total number of infected with the applications of the optimal controls. The corresponding controls (u1, u2, u3) are depicted in Figures 4(a)-4(c). The time dependent control u1(t) is observed to be at its upper EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 17 Figure 3. Simulation of the leishmaniasis model (4.1) as a function of time without control (blue curves) and with optimal control (pink curves) for: (a) Total number of infected humans; (b) Total number of infected reserviors; (c) Total number of infected sandflies bound for about 2 months and then can be reduced to get the optimal results. Similarly, the control u2(t) is observed to be at its upper bound for about 3 and half months and is tapered down before gradually reducing and then turns off at 124 days. Finally, the control u3(t) is at the upper bound for about four and half months and can be reduced gradually. The total number of infectious humans (IH + PH) are at highest level (about 188) in the absence of time dependent control and decrease gradually but in the presence of control, the number of infectious humans start to decrease from the beginning and becomes less than one at the end of 95 days. The total number of infectious sandflies peaked at 9th day with value 181 and decreases with no control but the value decreases from the beginning in the presence of control and becomes zero in about a month. The number of infectious canine reservoir peaks at 48th day with value 142 and slightly decreases thereafter if no controls are applied. At the end of 6 months with no controls, the number of infected dogs is about 120. This implies that a huge number in the reservoir is present at the end of six months if no controls measures are applied. With optimal control the number of infectious reservoir decreases from the beginning and becomes zero on 50th day. The value of objective functional J without optimal control is 3.5128× 105 and with optimal control J is about 70% less which is 1.0701× 105. 18 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 Figure 4. Optimal controls of the leishmaniasis model (4.1) for: (a) Personal protection control; (b) Sandflies insecticide control; (c) Reser- voir culling control. Control within varying values of parameters ΛR, ΛS, bSR and τH . To carry out the sensitivity analysis test in Section 3, we set the values of the parameters to be within ±20% range of their baseline values and found the parameters τH , ΛR, ΛS and bSR among others to have a high impact on the reproduction number R0. However, the time dependent controls u1(t), u2(t) and u3(t) incorporated into system (4.1) did not directly act on these parameters. Hence, we vary these four parameters using a range of ± 50% from their baseline values and determine their impact on the system as well as the control variables. We have used the ± 50% range for emphasis purpose, the result is expected to be similar in the ± 20% range. All other parameters in the model 4.1 and the balancing parameters in the objective functional are taken to be in the baseline values from Table 2. The results of the optimal control simulations of the leishmaniasis model (4.1) in the uncertainty interval for these parameters ΛR,ΛS and bSR are depicted in Figures 5(a)- refs1(c) where the total number of infectious humans, sandflies and canine reservoir in the absence of controls are denoted by blue curves and the corresponding numbers with optimal controls are denoted by pink curve. These parameters are related to the canine reservoirs and sandflies. From Figure 5, we observed that the upper bounds of these three parameters lead to slightly bigger number of infectious humans, canine reservoirs and sandflies. However, these numbers are not very different and follows the same patterns as with the baseline values of these parameters. Also, as observed in 6(a)-6(c), the control profiles within the uncertainty interval considered are almost the same. Hence, regardless of the parameter values chosen within the uncertainty interval of EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 19 Figure 5. Simulation of the Leishmaniasis model (4.1) as a function of time without control and with optimal control for: (a) Total number of infected humans; (b) Total number of infected reserviors; (c) Total number of infected sandflies. The parameters ΛR,ΛS and bSR are set within ±50% interval of their baseline values given in Table 2 and all other parameters are taken in their baseline values. In all plots, the solid curves represent the plot for baseline values of all parameters, the dotted curves represent the plot for lower bound of the three parameters and the dashed curves represent the plot for upper bound of these three parameters. The blue curves are without controls and the pink curves are with the optimal controls. these parameters (ΛR,ΛS and bSR), the control profile remained the same although the actions of the controls are higher for the upper bound and lower for the lower bound. Note, that we have used very low weights, the profile is the same even for higher weights and tighter intervals as the ±20% uncertainty interval used for the sensitivity analysis. Next, we investigate the optimal control of system (4.1) within the uncertainty interval of parameter τH ; this parameter is the human recovery rate and we use the ±50% range from the baseline value. The results are depicted in Figures 7, where the blue curves represent total number of infectious humans, sandflies and canine reservoir in the absence of controls and the corresponding numbers with optimal controls are denoted by pink curves. We observed in Figure 7(a) that more infected humans at the lower bound and fewer at the upper bound of τH . Similarly for the optimal controls u1(t), u2(t) and u3(t) in Figures 7(b)-7(d). More control efforts are required when humans recover slowly, particularly in the use of insecticide spray rate (u1(t)) and dog culling rate (u3(t)). In summary, numerical simulations of the leishmaniasis control model (4.1) show that leishmaniasis can be reduced in the community by the application of the time dependent control triple (u1(t), u2(t) and u3(t)) which represent respectively the of personal protection(u1(t)) such as insect repellant, door and window screens against sand flies, insecticide spraying (u2(t)) and culling of infected canine reservoir (u3(t)). 20 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 Figure 6. Optimal controls of the leishmaniasis model (4.1) in the ±50% uncertainty interval of parameters ΛR,ΛS and bSR for: (a) Per- sonal protection control; (b) Sandflies insecticide control; (c) Reservoir culling control. The parameters ΛR,ΛS and bSR are set within ±50% interval of their baseline values given in Table 2 and all other parameters are taken in their baseline values. In all plots, the solid curves represent the plot for baseline values of all parameters, the dotted curves represent the plot for lower bound of the three parameters and the dashed curves represent the plot for upper bound of these three parameters. 6. Conclusions and discussion In this article, a new deterministic model is designed and used to study the transmission dynamics of leishmaniasis model. The model includes the transmission dynamics in humans and canine reservoirs. The study shows that the disease-free equilibrium of the model is locally- and globally-asymptotically stable whenever the associated reproduction number (R0, an epidemiological threshold quantity that measures the spreading capacity of the disease), is less than unity and unstable otherwise. This study identifies (via sensitivity analysis) the significant parameters using as model outcome the reproduction number. The parameters with the largest im- pact are the natural death rate of sandflies (µS), sandflies biting rate on the canine reservoir (bSR), transmission probability per contact for susceptible sandflies (βS), recruitment rate of sandflies (ΛS), sandflies biting rate in humans (bS), transmission probability per contact for susceptible canine reservoir (βR), transmission probabil- ity per contact for susceptible humans (βH), vaccination efficacy in canine reservoir EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 21 Figure 7. Optimal controls of the leishmaniasis model (4.1) in the ±50% uncertainty interval of τH for: (a) Total number of infected hu- mans; (b) Personal protection control; (c) Sandflies insecticide control; (d) Reservoir culling control. All of the other parameters are in their baseline values given in Table 2. In all plots, the solid curves represent the plot for baseline values of all parameters, the dotted curves represent the plot for lower bound of the three parameters and the dashed curves represent the plot for upper bound of these three parameters. (ε), recruitment rate of canine reservoir (ΛR), and the disease-induced death rate of infected reservoir (δR). The identification of these key parameters is vital to the formulation of effec- tive control strategies for combating the spread of the disease. Thus, following the results obtained from the sensitivity analysis, we introduce three time depen- dent controls into system (2.1). The time dependent controls represent the use of personal protection such as insect repellant, door and window screens against sand flies in the bid to reduce human contact with the flies, insecticide spraying with insecticide such as pyrethroids are effective against sandflies both indoors and outdoors and culling of infected canine reservoir. Results from the numerical simu- lations indicates that that leishmaniasis can be controlled in the community by the application of these time dependent controls. Hence, in this article, we formulated and analyzed a system of ordinary differen- tial equations for leishmaniasis in humans and canine reservoirs. Some of theoretical and epidemiological findings of this study are summarized below: (i) The leishmaniasis model (2.1) is locally-and globally-asymptotically stable (LAS) when R0 < 1 and unstable when R0 > 1; (ii) The sensitivity analysis of the model shows that the significant parameters are the natural death rate of sandflies (µS), sandflies biting rate on the canine reservoir (bSR), transmission probability per contact for susceptible sandflies (βS), recruitment rate of sandflies (ΛS), sandflies biting rate in humans (bS), transmission probability per contact for susceptible canine 22 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 reservoir (βR), transmission probability per contact for susceptible humans (βH), vaccination efficacy in canine reservoir (ε), recruitment rate of canine reservoir (ΛR), and the disease-induced death rate of infected reservoir (δR) (iii) Numerical simulations indicates that that leishmaniasis can be controlled in the community by the application of time dependent controls representing personal protection (u1(t)), insecticide spraying (u2(t)) and culling infected canine reservoir (u3(t)) respectively. (iv) The control profiles were approximately the same for all parameter values chosen within the uncertainty interval of these parameters (ΛR,ΛS and bSR). (v) In the uncertainty interval for parameter τH , more control efforts are re- quired when humans recover more slowly. Acknowledgments. This work was assisted through participation of the authors in an Investigative Workshop on Malaria-Leishmania Co-infection (May 26-28, 2015) at the National Institute for Mathematical and Biological Synthesis, an insti- tute sponsored by the National Science Foundation, the U.S. Department of Home- land Security, and the U.S. Department of Agriculture through NSF Award #EF- 0832858, with additional support from The University of Tennessee, Knoxville. The research of I.M. ELmojtaba is supported by Sultan Qaboos University grant number: IG/SCI/DOMS/16/16. References [1] Alexander, B.; Maroli, M.; Control of phlebotomine sandflies. Med Vet Entomol, 17 (2003)(1): 1–18. [2] Beneth, G.; Koutinas, A. F.; Laia, S.G.; Bourdeau, P.; Ferrer, L.; Canine leishmaniosis - new concepts and insights on an expanding zoonosis: part one Trends in Parasitology24 (2008) (7). [3] Blower, S. M.; Dowlatabadi, H.; Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int. Stat. Rev., 2 (1994), 229–243. [4] Bortoletto, D. V.; Utsunomiya, Y. T.; Perri, S. H. V.; Ferreira, F.; Nunes, C. M.; Age structure of owned dogs under compulsory culling in a visceral leishmaniasis endemic area. Cadernos de Sau de Publica, 32 (2016) (8). [5] Bradley, T.; King, R.; The dog economy is global but what is the world’s true canine capital? The Atlantic, http: //www.theatlantic.com/business/archive/2012/11/the-dog- economy-is-global-but-what-is-the-worlds-true-canine-capital/265155/ (2012). Accessed De- cember 30, 2015. [6] Castillo-Chavez, C.; Feng, Z.; Huang, W.; On the computation of R0 and its role on global sta- bility. In Mathematical approaches for emerging and re emerging infectious diseases. An In- troduction. Castillo-Chavez, C., Blower, S., van den Driessche, P., Kirschner, D. and Yakubu, A.-A. (Eds.) Springer-Verlag, New York (2002), 229-250. [7] Centers for Disease Control and Prevention (2013). Leishmaniasis. http://www.cdc.gov/parasites/leishmaniasis/ (2002). Accessed January 4, 2016. [8] Chaves, L. F.; Hernandez, M. J.; Mathematical modelling of American Cutaneous Leishma- niasis: incidental hosts and threshold conditions for infection persistence. ACTA TROPICA (Elsevier). 92 (2004), 255–252. [9] Costa, D. N. C. C.; Codeco, C. T.; Silva, M. A.; Werneck, G. L.; Culling Dogs in Scenarios of Imperfect Control: Realistic Impact on the Prevalence of Canine Visceral Leishmaniasis. PLOS: Neglected Tropical Disease . 7 (8) (2013), 1–8. [10] Deer, H. M.; Pesticide adsorption and half-life. AG/Pesticides, 15 (1999), 1. [11] Diekmann, O.; Heesterbeek, J. A. P.; Metz, J. A. J.; On the definition and computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol., 28 (1990), 365-382. EJDE-2020/80 A VISCERAL LEISHMANIASIS MODEL 23 [12] Diniz, S. A.; Silva, F. L.; Carvalho Neta, A. V.; Bueno, R., Guerra, R. M. S. N. C.; Abreu- Silva, A. L.; Santos, R. L.; Animal reservoirs for visceral leishmaniasis in densely populated urban areas. The Journal of Infection in Developing Countries. 2(1) (2008), 24–33. [13] Elmojtaba, I. M.; Mugisha, J. Y. T.; Hashim, M. H. A.; Mathematical analysis of the dynam- ics of visceral leishmaniasis in Sudan. Applied Mathematics and Computation 217(6) (2010), 2567-2578. [14] Ferrell, M. A.; Aagard, S. D.; Pesticide Adsorption and Half Life. University of Wyoming, Cooperative Extension Service, Department of Plant Sciences, College of Agriculture, 2003. [15] Gasim, S.; Elhasan, A. M.; Kharazmi, A.; Khalil, E. A. G.; Ismail, A.; Theander, T. G.; The development of post-kala-azar dermal leishmaniasis (PKDL) is associated with acquisi- tion of Leishmania reactivity by peripheral blood mononuclear cells (PBMC). Clinical and experimental Immunology, 119 (3) (2000), 523-529. [16] Hyman, J. M.; Li, J.; An intuitive formulation for the reproductive number for the spread of diseases in heterogeneous populations. Mathematical Biosciences, 167 (2000), 65-86. [17] Joshi, H. R.; Optimal Control of an HIV Immunology Model, Optim. Control Appl. Math, 23 (2002), 199–213. [18] Jung, E.; Lenhart, S.; Feng, Z.; Optimal Control of Treatments in a Two-Strain Tuberculosis Model, Discrete and Continuous Dynamical Systems-Series B. 2(4) (2002), 473–482. [19] Kamgang, J. C.; Sallet, G.; Computation of threshold conditions for epidemiological models and global stability of the disease-free equilibrium (DFE). Mathematical Biosciences, 213 (2008), 1-12. [20] Kasap, O. E.; Alten, B.; Comparative demography of the sand fly Phlebotomus papatasi (Diptera: Psychodidae) at constant temperatures. Journal of Vector Ecology, 31(2) (2006), 378-385. [21] Kern, D.; Lenhart, S.; Miller, R.; Yong, J.; Optimal Control Applied to Native-Invasive Population Dynamics. J Biol Dyn. 1 (4) (2007), 413–26. doi: 10.1080/17513750701605556 [22] Kirschner, D.; Lenhart, S.; Serbin, S.; Optimal Control of the Chemotherapy of HIV, J. Math. Biol. 35, (1997), 775–792. [23] Kumar, R.; Engwerda, C.; Vaccines to prevent leishmaniasis. Clinical & Translational Im- munology . 3(e13) (2014), 1–6. [24] Lemesre, J. L.; Holzmuller, P.; Gonçalves, R. B.; Bourdoiseau, G.; Hugnet, C.; Cavaleyra, M.; Papierok, G.; Long-lasting protection against canine visceral leishmaniasis using the LiESAp-MDP vaccine in endemic areas of France: double-blind randomised efficacy field trial. Vaccine. 25(21) (2007), 4223–4234. [25] Lenhart, S.; Workman, J. T.; Optimal Control Applied to Biological Models. Chapman and Hall/CRC. First Edition, 1998. [26] Lukes, D. L.; Differential Equations: Classical to Controlled. Mathematics in Science and Engineering. Academic Press, 1982. [27] Luo, Y.; Zhang, M.; Environmental modeling and exposure assessment of sediment-associated pyrethroids in an agricultural watershed. PloS one, 6(1) (2011), e15794. [28] Macdonald, G.; The Epidemiology and Control of Malaria, Oxford University Press, London, 1957. [29] Marino, S.; Hogue, I. B.; Ray, C. J.; Kirschner, D. E.; A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 254 (2008), 178–196. [30] Mckay, M. D.; Beckman, R. J.; Conover, W. J.; Comparison of 3 methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (1979), 239–245. [31] McLeod, R. G.; Brewster, J.F.; Gumel, A. B.; Slonowsky, D. A.; Sensitivity and uncertainty analyses for a sars model with time-varying inputs and outputs. Math. Biosci. Eng. 3(3) (2006), 527–544. [32] Palantik-de-sousa, C. B.; Batista-de-Melo, L. M.; Borja-Cabrera, G. P.; Plantik, M.; Lavor, C. C.; Improving methods for epidemiological control of canine visceral leishmaniasis based on a mathematical model. Impact on the incidence of the canine and human disease. Annals of Brazilian academy of Science. 76(3) (2004), 583–593. [33] Pontryagin, L. S.; Boltyanskii, V. G.; Gamkrelidze, R. V.; Mishchenko, E. F.; The mathe- matical theory of optimal processes. Wiley, New York, 1962. 24 B. PANTHA, F. B. AGUSTO, I. M. ELMOJTABA EJDE-2020/80 [34] Ribas, L. M; Shimozako, H. J.; Massad, E.; Estimating the Optimal Control of Zoonotic Visceral Leishmaniasis by the Use of a Mathematical Model. Hindawi: The Scientific World Journal. 2013 (2013), 1–6. [35] Ross, R. ; The Prevention of Malaria, John Murray, Oxford, London, 1911. [36] Singh, N.; Mishra, J.; Singh, R.; Singh S.; Animal Reservoirs of Visceral Leishmania in India. Journal of Parasitology, 99 (1) (2013), 64–67. [37] Sanchez, M.; Blower, S.; Uncertainty and sensitivity analysis of the basic reproduction rate. American Journal of Epidemiology, 145 (12) (1997), 1127–1137. [38] Standford University online resource on Visceral Leishmaniasis (2006). https://web.stanford.edu/class/humbio103/ParaSites2006/Leishmaniasis/visceral.htm . Ac- cessed on Sept. 19, 2015. [39] Stauch, A.; Duerr, H.-P.; Dujardin, J.-C.; Vanaerschot, M.; Sundar, S, et al.; Treat- ment of Visceral Leishmaniasis: Model-Based Analyses on the Spread of Antimony- Resistant L. donovani in Bihar, India. PLoS Negl Trop Dis 6 (12) (2012), e1973. doi:10.1371/journal.pntd.0001973 [40] Stauch, A.; Duerr, H. P.; Picado, A.; Ostyn, B.; Sundar, S.; Model-Based Investigations of Different Vector-Related Intervention Strategies to Eliminate Visceral Leishmaniasis on the Indian Subcontinent. PLoS Negl Trop Dis 8 (4) (2014), e2810. doi: 10.1371/jour- nal.pntd.0002810 [41] Subramanian, A.; Singh, V.; Sarkar, R. R.; Understanding Visceral Leishmaniasis Disease Transmission and its Control-A Study Based on Mathematical Modeling. Mathematics 3 (3) (2015), 913-944. doi:10.3390/math3030913 [42] van den Driessche, P.; Watmough, J.; Reproduction numbers and the sub-treshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180 (2002), 29-48. [43] Farrington, C. P.; Whitaker, H. J.; Estimation of effective reproduction numbers for infectious diseases using serological survey data. Biostatistics, 4 (4) (2003), 621-632. [44] World Bank; World Development Indicators: Population dynamics (2015) http://wdi.worldbank.org/table/2.1. Accessed January 20, 2016. [45] World Health Organization; Mosquitos and other biting Diptera. http://www.who.int/water sanitation health/resources/vector007to28.pdf (2015). Accessed January 20, 2016. [46] World Health Organization; Leishmaniasis (2015). http://www.who.int/mediacentre/factsheets/fs375/en/. Accessed January 4, 2016. [47] Yan, X.; Zou, Y.; Li, J.; Optimal quarantine and isolation strategies in epidemics control. World Journal of Modelling and Simulation. 3 3 (2007), 202–211. [48] Yan, X.; Zou, Y.; Optimal and sub-optimal quarantine and isolation control in SARS epi- demics. Mathematical and Computer Modelling. 47 (2008), 235–245. Buddhi Pantha Department of Science and Mathematics, Abraham Baldwin Agricultural College, Tifton, GA 31793, USA Email address: bpantha@abac.edu Folashade B. Agusto Department of Ecology and Evolutionary Biology, University of Kansas, Lawrence KS 66045, USA Email address: fbagusto@gmail.com Ibrahim M. Elmojtaba Department of Mathematics, Sultan Qaboos University, Al Khoudh 123 - Muscat, Oman Email address: elmojtaba@squ.edu.om 1. Introduction 2. Model formulation Invariant region 2.1. Analysis of the Model 3. Sensitivity analysis 4. Optimal control problem Characterization of optimal controls Setting bounds on the control 5. Numerical results Control within varying values of parameters R, S, bSR and H 6. Conclusions and discussion Acknowledgments References