\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{tikz}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2020 (2020), No. 80, pp. 1--24.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2020 Texas State University.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2020/80\hfil A visceral leishmaniasis model]
{Optimal control applied to a visceral leishmaniasis model}
\author[B. Pantha, F. B. Agusto, I. M. Elmojtaba \hfil EJDE-2020/80\hfilneg]
{Buddhi Pantha, Folashade B. Agusto, Ibrahim M. Elmojtaba}
\address{Buddhi Pantha \newline
Department of Science and Mathematics,
Abraham Baldwin Agricultural College, \newline
Tifton, GA 31793, USA}
\email{bpantha@abac.edu}
\address{Folashade B. Agusto \newline
Department of Ecology and Evolutionary Biology,
University of Kansas,
Lawrence KS 66045, USA}
\email{fbagusto@gmail.com}
\address{Ibrahim M. Elmojtaba \newline
Department of Mathematics,
Sultan Qaboos University,
Al Khoudh 123 - Muscat, Oman}
\email{elmojtaba@squ.edu.om}
\dedicatory{Communicated by Suzanne Lenhart}
\thanks{Submitted January 6, 2019. Published July 28, 2020.}
\subjclass[2010]{35F21, 37N25}
\keywords{Visceral leishmanisis; PKDL; vaccination; canine reservoir;
\hfill\break\indent optimal control}
\begin{abstract}
In this article, we developed a deterministic model for the transmission 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 dependent 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 representing respectively, personal protection,
insecticide spraying and culling infected canine reservoir.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks
\section{Introduction}
Leishmaniasis is a vector borne parasitic disease caused by a protozoan parasite
called \emph{Leishmania} which is transmitted by a bite of an infected female
phlebotomine sandflies \cite{CDC2016, WHO2016}.
These sandflies species feed on animal blood which is required for egg development \cite{WHO}.
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
\cite{Elmojtaba2010}. 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 \cite{CDC2016, WHO2016}.
The species of the Leishmania parasite determines the symptoms of the disease such
as cutaneous, mucocutaneous and visceral \cite{CDC2016, WHO2016}.
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 \cite{CDC2016}, 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 \cite{WHO2016}.
There are over 30 sandflies species which cause human and animal leishmaniasis.
The clinical symptoms and disease progression depend on Leishmania species,
host immune system and factors such as environmental conditions, seasons,
socio-economic status \cite{CDC2016,WHO2016}. 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 \cite{CDC2016, Chaves2004,WHO2016}.
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 \cite{CDC2016,WHO2016}. In the past 5 years, about 1 million CL
cases were reported worldwide, with the majority occurring in Americas,
the Mediterranean basin, the Middle East and Central Asia \cite{CDC2016,WHO2016}.
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) \cite{CDC2016,WHO2016}.
The patients co-infected with leishmaniasis and other disease such as HIV may also have MCL
infection \cite{CDC2016,WHO2016}. MCL does not heals all spontaneously but it can be
cured with timely treatment \cite{CDC2016,WHO2016}. Almost 90\% of MCL cases occurs in
the Plurinational State of Bolivia, Brazil and Peru \cite{WHO2016}.
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
\cite{CDC2016,Elmojtaba2010, WHO2016}. Caused by either of \emph{L. donovani}
or \emph{L. infantum} \cite{Stauch2012,WHO2016} parasites, it can infect people of
any ages but children under 10 years are more susceptible \cite{CDC2016,WHO2016}.
The incubation period of visceral leishmanasis ranges from 10 days to several months
and the onset of the symptoms is usually gradual \cite{CDC2016,WHO2016}.
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 \cite{CDC2016, WHO2016} 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
\cite{Stauch2012, WHO2016}. 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 between
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 \cite{Diniz2008,Singh2013,WHO2016}.
Some zoonotic reservoirs 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 \cite{CDC2016, WHO2016}. In an endemic area, an animal species
is usually a principle reservoir host for a particular Leishmania species but other animal
species in the same area may be an incidental reservoirs
\cite{ CDC2016,Diniz2008,Singh2013, WHO2016}. In some VL endemic areas, seroprevalence
test suggested that goats are reservoir host for VL \cite{Diniz2008, Singh2013}, however,
many studies \cite{Diniz2008,Singh2013,WHO2016,Palatnik-de-Sousa2004} found that dogs
are the principal reservoir host of \emph{Leishmania infantum}, a major causative agent
of visceral leishmaniasis, and more than 50\% of all infected dogs are asymptotic carriers.
Because of technological advancement, many technologies for disease diagnosis, prevention and
treatment for a wide range of diseases have been developed. However, population growth,
movement, environmental changes and urbanization have created favorable epidemiological
conditions for disease transmission and persistence of some diseases such as
leishmaniasis \cite{Diniz2008}. In fact, because of the rarity of functioning control programs,
the new cases of leishmaniasis is in the increasing trend \cite{WHO2016}.
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 \cite{Elmojtaba2010,Stauch2012,Stauch2014,Subramanian2015,WHO2016}.
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 \cite{WHO2016}. 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
\cite{CDC2016, WHO2016}. In some cases, culling of seropositive reservoir such as dogs
has been implemented to reduce this interaction. Some vaccine have recently been introduced
but its efficacy in human is yet to be determined \cite{Diniz2008, Kumar2014, Lemesre2007}.
Canine vaccine has shown promising decrease in infection in some countries such as Brazil,
France and Iran but no vaccines are available for humans \cite{Diniz2008,WHO2016}.
Numerous mathematical studies have been carried out on leishamaniasis. Chaves et al.\
\cite{Chaves2004} 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.\
\cite{Ribas2013} developed an ODE model that includes constant control strategies
(treatment, 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.\ \cite{Palatnik-de-Sousa2004}
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 \cite{Palatnik-de-Sousa2004}. Some other studies such as Elmojtaba et al.\
\cite{Elmojtaba2010}, Stauch et al.\ \cite{Stauch2012, Stauch2014}, Subramanian et al.\
\cite{Subramanian2015} also explored the combination 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 example \cite{Palatnik-de-Sousa2004}. 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 reservoir populations. The detailed assumptions of our model
are presented in Section \ref{Formulation}.
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 \ref{Formulation}, 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 \ref{Sensitivity}, 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 \ref{Control}, characterize the
optimal controls, and present the numerical results in Section \ref{Numerical}.
Section \ref{Conclusion} presents conclusions and discussions of our study.
\section{Model formulation} \label{Formulation}
To formulate our model we consider three different populations: human host population,
$N_H(t)$, canine host population, $N_R (t)$, and vector population, $N_S (t)$.
The human host population is divided into four sub-populations: susceptible individuals
$S_H (t)$, visceral leishmaniasis infected individuals $I_H (t)$, those who develop PKDL
after the treatment of visceral leishmaniasis, $P_H (t)$, and those who are recovered
and have permanent immunity, $R_H (t)$. This implies that
$$
N_H (t) = S_H (t) + I_H (t) + P_H (t) + R_H (t).
$$
Similarly, the canine host population is divided into three sub-populations:
susceptible, $S_R (t)$, vaccinated, $V_R (t)$ and infected, $I_R (t)$, such that
$$
N_R (t) = S_R (t) + V_R (t) + I_R (t),
$$
and finally, the sandfly population have two sub-populations: susceptible, $S_S (t)$,
and infected sandflies, $I_S (t)$, such that
$$
N_S (t) = S_S (t) + I_S (t).
$$
It is assumed that susceptible individuals are recruited into the population at a constant
rate $\Lambda_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 $\beta_H b_S \frac{I_S}{N_H}$, where $\beta_H$ is the transmission
probability per bite per human (as the case for malaria, \cite{Macdonald1957,Ross1911})
and $b_S$ is the per capita biting rate of sandflies on humans. Infected humans die due
to leishmaniasis at an average rate $\delta_H$, or get treatment at an average rate
$\gamma_H$, and a fraction $\alpha$ of them recover and acquire permanent immunity,
and the complement fraction $(1 - \alpha)$ develop PKDL. Humans with PKDL get treated
at an average rate $\tau_H$, or recover naturally at an average rate $\sigma_H$, and
acquire permanent immunity in both cases. There is a per capita natural mortality rate
$\mu_H$ in all human sub-population.
Susceptible sandflies are recruited at a constant rate $\Lambda_S$, and acquire
leishmaniasis infection following contacts with a leishmaniasis infected human or with
a human having PKDL or leishmaniasis infected canine at an average rate equal to
$\beta_S b_S\frac{I_H}{N_H} + \beta_S b_S \frac{P_H}{N_H} + \beta_S b_{SR} \frac{I_R}{N_R}$,
where $b_S$ is the per capita biting rate on human and $b_{SR}$ is the per capita biting
rate on canine, and $\beta_{S}$ is the transmission probability for sandfly infection
after biting a human or a canine. Sandflies suffer natural mortality at a per capita
rate $\mu_S$ regardless of their infection status.
Susceptible canine reservoirs are recruited into the population at a constant rate
$\Lambda_R$, acquire infection with leishmaniasis following contacts with infected
sandflies at a rate $\beta_R b_{SR} \frac{I_S}{N_R}$ where $\beta_R$ and $b_{SR}$
as described above; or vaccinated at an average constant rate $\nu_R$.
Vaccinated canine reservoir acquire infection with leishmaniasis following contacts
with infected sandflies at a rate $\beta_R b_{SR} (1 - \epsilon) \frac{I_S}{N_R}$
where $\epsilon$ 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 $\delta_R$.
Using the above description, the following mathematical model is proposed
\begin{equation} \label{eL}
\begin{aligned}
S_H' &= \Lambda _H - \frac{\beta_Hb_S I_SS_H}{N_H} - \mu_H S_H \\
I_H' &= \frac{\beta_Hb_S I_SS_H}{N_H}- (\gamma_H + \mu_H + \delta_H) I_H \\
P_H' &= (1-\alpha) \gamma_H I_H - (\tau_H + \sigma_H + \mu_H) P_H \\
R_H' &= \alpha\gamma_H I_H + (\tau_H + \sigma_H) P_H - \mu_H R_H \\
S_S' &= \Lambda_S - \frac{\beta_Sb_S(I_H+P_H) S_S}{N_H} - \frac{\beta_Sb_{SR} I_R S_S}{N_R}
- \mu_S S_S \\
I_S' &= \frac{\beta_Sb_S(I_H+P_H) S_S}{N_H} + \frac{\beta_Sb_{SR}I_R S_S}{N_R} - \mu_S I_S\\
S_R' &= \Lambda_R - \frac{\beta_Rb_{SR}I_SS_R}{N_R} - (\nu_R +\mu_R) S_R \\
V_R' &= \nu_R S_R - \frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R} - \mu_R V_R \\
I_R' &= \frac{\beta_Rb_{SR}I_SS_R}{N_R} + \frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R} - (\mu_R + \delta_R) I_R
\end{aligned}
\end{equation}
with initial conditions
\begin{gather*}
S_H(0)=S_H^0, \quad I _H(0)=I_H^0, \quad P_H(0)=P^0_H,\quad R_H(0)=R^0_H, \\
S_V(0)=S_V^0,\quad I_V(0)=I_V^0,\quad S_R(0)=S^0_R,\quad I_R(0)=I^0_V.
\end{gather*}
and
\[
N_H' = \Lambda_H - \mu_H N_H - \delta_H I_H ,\quad
N_S' = \Lambda_S - \mu_S N_S,\quad
N_R' = \Lambda_R - \mu_R N_R - \delta_R I_R.
\]
\subsection*{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.\ $\delta_H = \delta_R = 0)$, the total human population,
$N_H\to \Lambda_H/ \mu_H$ as $t\to \infty$. Similarly, $N_R\to\Lambda_R/ \mu_R$ and
$N_S \to \Lambda_S/ \mu_S$ as as $t\to \infty$. This shows that the biologically-feasible
region:
\begin{align*}
\Omega = \big\{&(S_H,I_H,P_H,R_H,S_S,I_S,S_R,V_R,I_R) \in \mathbb R^{9}_{+} :\\
&S_H,I_H,P_H,R_H,S_S,I_S,S_R,V_R,I_R \geq 0, \;
N_H \leq \frac{\Lambda_H}{\mu_H}, N_S \leq \frac{\Lambda_S}{\mu_S},
N_R \leq \frac{\Lambda_R}{\mu_r} \big\}
\end{align*}
is positively-invariant domain, and thus, the model is epidemiologically and
mathematically well posed, and it is sufficient to consider the dynamics of the
flow generated by the system \eqref{eL} in this positively-invariant domain $\Omega$.
The flow diagram of the system \eqref{eL} is depicted in Figure \ref{p1} and the
associated variables and parameters are described in Table \ref{t1}.
\begin{figure}[ht]
\begin{center}
\tikzstyle{block}=[rectangle,draw,fill=orange!40,text width=2em, text centered, rounded corners, minimum height=3em, minimum width=3em]
\tikzstyle{nblock}=[fill={white},text width={}, rounded corners]
\begin{tikzpicture}[node distance=2.0cm, auto, >=stealth]
% Humans
\coordinate(IH) at (0,0);
\node[block,fill = brown!50] (IH) {\Large$I_H$};
\node[block,fill = yellow!50](SH)[left of = IH, node distance=4cm] {\Large$S_H$};
\node[nblock](S)[right of = IH, node distance=4cm] {};
\node[block,fill = blue!20](RH)[below of=S, node distance=3cm] {\Large$R_H$};
\node[block](PH)[above of =S, node distance=3cm] {\Large$P_H$};
\draw[->,line width=0.37mm](SH)--node[pos=.5,sloped,above]{{\scriptsize$\beta_Hb_S \frac{I_S}{N_H}$}}(IH);
\draw[->,line width=0.37mm](IH)--node[pos=.5,sloped,above]{{\scriptsize$\alpha\gamma_H$}}(RH);
\draw[->,line width=0.37mm](IH)--node[pos=.5,sloped,below]{{\scriptsize$(1-\alpha)\gamma_H$}}(PH);
\draw[->,line width=0.37mm](PH)--node[pos=.5,right]{{\scriptsize$\tau_H$}}(RH);
\draw[->,line width=0.37mm](PH)--node[pos=.5,left]{{\scriptsize$\sigma_H$}}(RH);
% sandflies
\node[block,fill = yellow!50](SS) [below of=SH, node distance=5.5cm] {\Large$S_S$};
\node[block,fill = brown!50](IS)[below of=RH, node distance=2.5cm] {\Large$I_S$};
\draw[->,line width=0.37mm](SS)--node[pos=.5,sloped,above]{{\scriptsize$\beta_Sb_S \frac{(I_H+P_H)}{N_H}$}}(IS);
\draw[->,line width=0.37mm](SS)--node[pos=.5,sloped,below]{{\scriptsize$\beta_Sb_{SR} \frac{I_R}{N_R}$}}(IS);
% Reservoirs
\node[block,fill = yellow!50](SR)[below of=SS, node distance=4.5cm] {\Large$S_R$};
\node[block,fill = brown!50](IR)[below of=S, node distance=12.5cm] {\Large$I_R$};
\node[block,fill = red!10](VR)[below of=IS, node distance=2.5cm] {\Large$V_R$};
\draw[->,line width=0.37mm](SR)--node[pos=.5,sloped,above]{{\scriptsize$\nu_R$}}(VR);
\draw[->,line width=0.37mm](SR)--node[pos=.5,sloped,above]{{\scriptsize$\beta_Rb_{SR} \frac{I_S}{N_R}$}}(IR);
\draw[->,line width=0.37mm](VR)--node[pos=.5, right]{{\scriptsize$\beta_Rb_{SR}(1-\varepsilon)\frac{IS}{N_R}$}}(IR);
% Humans
\node [] (birth) [left of=SH, node distance=1.5cm]{};
\draw[->] (birth) edge node[pos=.5,above]{{\scriptsize$\Lambda_H$}}(SH) ;
%
\node [] (death) [below of=SH, node distance=1.5cm]{};
\draw[<-] (death) edge node[pos=.5,left]{{\scriptsize$\mu_H$}}(SH) ;
%
\node [] (death) [below of=IH, node distance=1.5cm]{};
\draw[<-] (death) edge node[pos=.5,left]{{\scriptsize$\mu_H$}}(IH) ;
%
\node [] (death) [below of=IH, node distance=1.5cm]{};
\draw[<-] (death) edge node[pos=.5,right]{{\scriptsize$\delta_H$}}(IH) ;
%
\node [] (death) [right of=RH, node distance=1.5cm]{};
\draw[<-] (death) edge node[pos=.5,above]{{\scriptsize$\mu_H$}}(RH) ;
%
\node [] (death) [right of=PH, node distance=1.5cm]{};
\draw[<-] (death) edge node[pos=.5,above]{{\scriptsize$\mu_H$}}(PH) ;
%
%
% sandflies
\node [] (birthS) [left of=SS, node distance=1.5cm]{};
\draw[->] (birthS) edge node[pos=.5,above]{{\scriptsize$\Lambda_S$}}(SS) ;
%
\node [] (deathS) [below of=SS, node distance=1.5cm]{};
\draw[<-] (deathS) edge node[pos=.5,left]{{\scriptsize$\mu_S$}}(SS) ;
%
\node [] (deathS) [right of=IS, node distance=1.5cm]{};
\draw[<-] (deathS) edge node[pos=.5,above]{{\scriptsize$\mu_S$}}(IS) ;
%
%
% Reservoirs
\node [] (birthM) [left of=SR, node distance=1.5cm]{};
\draw[->] (birthM) edge node[pos=.5,above]{{\scriptsize$\Lambda_R$}}(SR) ; %
%
\node [] (deathM) [below of=SR, node distance=1.5cm]{};
\draw[<-] (deathM) edge node[pos=.5,left]{{\scriptsize$\mu_R$}}(SR) ;
%
\node [] (deathM) [right of=VR, node distance=1.5cm]{};
\draw[<-] (deathM) edge node[pos=.5,above]{{\scriptsize$\mu_R$}}(VR) ;
%
\node [] (death) [right of=IR, node distance=1.5cm]{};
\draw[<-] (death) edge node[pos=.5,above]{{\scriptsize$\mu_R$}}(IR) ;
%
\node [] (death) [right of=IR, node distance=1.5cm]{};
\draw[<-] (death) edge node[pos=.5,below]{{\scriptsize$\delta_R$}}(IR) ;
\end{tikzpicture}
\end{center}
\caption{\small Systematic flow diagram of model \eqref{eL}.}
\label{p1} %\label{flowchart}
\end{figure}
\begin{table}[htb]
\caption{Description of the variables and parameters for model \eqref{eL}}.
\label{t1}
\centering
\begin{tabular}{l l} \hline
Variable & Description \\ \hline
$S_H$ & Population of susceptible humans\\
$I_{H}$& Population of symptomatic humans \\
$P_{H}$ & Population of symptomatic humans \\
$R_H$ & Population of recovered humans \\
$S_S$ & Population of susceptible sandflies\\
$I_{S}$ & Population of infected sandflies \\
$S_R$ & Population of susceptible reservoir (canine)\\
$V_R$ & Population of vaccinated reservoir (canine)\\
$I_{R}$ & Population of infected reservoir (canine) \\[3pt]
\hline
Param. & Description \\ \hline
$\Lambda_H$ & Recruitment rate of humans\\
$\Lambda_S$ & Recruitment rate of sandflies\\
$\Lambda_R$ & Recruitment rate of reservoir (canine)\\
$\beta_H$ & Transmission probability per contact for susceptible humans\\
$\beta_S$ & Transmission probability per contact for susceptible sandflies\\
$\beta_R$ & Transmission probability per contact for susceptible reservoir (canine)\\
$b_S $ & Sandflies biting rate in humans\\
$b_{SR} $ & Sandflies biting rate in reservoir (canine) \\
$\mu_H$ & Natural death rate of humans\\
$\mu_S$ & Natural death rate of sandflies\\
$\mu_R$ & Natural death rate of reservoir (canine)\\
$\delta_{H}$ & Disease induced death rate of humans \\
$\delta_{R}$ & Disease induced death rate of symptomatic reservoir (canine) \\
$\gamma_H$ & Recovery rate of infected humans with VL dues to treatment\\
$\alpha$ & Fraction successfully treated for VL\\
$\tau_H$ & Recovery rate of humans with PKDL due to treatment\\
$\sigma_H$ & Instantaneous recovery rate of humans with PKDL\\
$\nu_R$ & Vaccination rate in reservoir (canine) \\
$\varepsilon$ & Vaccination efficacy in reservoir (canine) \\
\hline
\end{tabular}
\end{table}
\subsection{Analysis of the Model}
The disease-free equilibrium (DFE) of the system \eqref{eL} is given by
\[
\begin{aligned}
\mathcal{E}_0
&= \Bigl(S_H^{*},I_H^{*},P_H^{*},R_H^{*},S_S^{*},I_S^{*},S_R^{*},V_R^{*},I_R^{*}
\Bigr) \\
&= \Bigl(\frac{\Lambda_H}{\mu_H},0,0,0,\frac{\Lambda_S}{\mu_S},0,
\frac{\Lambda_R}{\mu_R + \nu_R},\frac{\nu_R\Lambda_R}{\mu_R(\mu_R + \nu_R)},0\Bigr).
\end{aligned}
\]
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 population
\cite{Diekmann1990,Hyman2000}. To calculate the effective reproduction number,
we will use the next generation approach \cite{Diekmann1990,VanDen2002}.
First, we take the variables $I_H,P_H,I_S,I_R$ as the infected compartments and then
use the notation in \cite{VanDen2002}, the Jacobian $F$ and $V$ matrices for new
infectious terms and the remaining transfer terms, respectively, are defined as:
\begin{gather*}
\mathbf{F} =
\begin{pmatrix}
0 & 0 & \beta_H b_S & 0\\
0 & 0 & 0 & 0\\
\frac{\beta_S b_S \Lambda_S \mu_H}{\mu_S\Lambda_H} & \frac{\beta_S b_S \Lambda_S \mu_H}{\mu_S\Lambda_H} & 0 & \frac{\beta_S b_{SR} \Lambda_S \mu_R}{\mu_S \Lambda_R}\\
0 & 0 & \frac{\beta_R b_{SR}[\mu_R + (1 - \varepsilon)\nu_R]}{\mu_R + \nu_R} & 0\\
\end{pmatrix}
\\
\mathbf{V} =
\begin{pmatrix}
k_1 & 0 & 0 & 0\\
-(1 - \alpha)\gamma_H & k_2 & 0 & 0\\
0 & 0 & \mu_S & 0\\
0 & 0 & 0 & k_3 \\
\end{pmatrix}
\end{gather*}
where
$k_1 = \gamma_H + \delta_H + \mu_H$, $k_2 = \tau_H + \sigma_H + \mu_H$,
$k_3 = \mu_R + \delta_R$.
Then the reproduction number is
\[
\mathcal{R}_0 = \rho(FV^{-1}) = \sqrt{(\mathcal{R}_H +\mathcal{R}_R)\mathcal{R}_S }
\]
where, $\rho$ is the spectral radius and
\begin{gather*}
\mathcal{R}_H = \frac{b_S^2\mu_H\beta_H[(1-\alpha)\gamma_H+k_2)]}{k_1k_2\Lambda_H},\\
\mathcal{R}_R = \frac{b_{SR}^2\beta_R\mu_R[\mu_R+(1-\varepsilon)\nu_R]}{k_3\Lambda_R
(\mu_R+\nu_R)},\\
\mathcal{R}_S = \frac{\beta_S\Lambda_S}{\mu_S^2}.
\end{gather*}
Furthermore, the quantity $\mathcal{R}_H $ is the number of secondary infections in humans
host by one infectious sandfly, $\mathcal{R}_R$ is the number of secondary infections
in the canine reservoir host by one introduced infectious sandfly, and lastly $\mathcal{R}_S$
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 \cite[Theorem 2]{VanDen2002}, the following result is established:
\begin{lemma} \label{lem1}
The disease-free equilibrium ($\mathcal{E}_0$) is locally asymptotically stable if
$\mathcal{R}_0 < 1$, and unstable if $\mathcal{R}_0 >1$.
\end{lemma}
To study the global behavior of system \eqref{eL}, we use
a theorem by Castillo-Chavez et al.\ \cite{Carlos1}, stated here
for convenience.
\begin{theorem}[\cite{Carlos1}] \label{thm1}
For the system
\begin{gather*}
\frac{dX}{dt} = F(X,Z) \\
\frac{dZ}{dt} = G(X,Z)\\
G(X,0) = 0\,,
\end{gather*}
where the components of the column-vector $X \in R^{m}$ denote the number of uninfected
individuals and the components of vector $Z \in R^{n}$ denote the number of infected
individuals. $U_0 = (X^{\ast},0)$ denotes the disease-free equilibrium of this system.
The fixed point $U_0 = (X^{\ast},0)$ is a globally asymptotically
stable equilibrium for this system provided that $\mathcal{R}_0 < 1$
(locally asymptotically stable) and the following two conditions satisfied:
\begin{itemize}
\item[(H1)] For $\frac{dX}{dt} = F(X,0)$, $X^{\ast}$ is globally asymptotically stable,
\item[(H2)] $G(X,Z) = AZ - \hat{G}(X,Z),\;\hat{G}(X,Z)\geq 0$ for
$(X,Z)\in\Omega$, where $A = D_{Z}G(X^{\ast},0)$ is an M-matrix
(off-diagonal elements of $A$ are non-negative) and
$\Omega$ is the region where the model has biological meaning.
\end{itemize}
\end{theorem}
We can rewrite our system \eqref{eL} using above notation, where
\begin{equation} \label{GS2}
X = (S_H,R_H,S_S,S_R,V_R),\quad
Z = (I_H,P_H,I_S,I_V)
\end{equation}
\begin{gather*}
F(X,0) = \begin{pmatrix}\Lambda_H - \mu_H S_H & 0 &
\Lambda_S - \mu_S S_S & \Lambda_R -(\nu_R + \mu_R) S_R & \nu_R S_R - \mu_R V_R
\end{pmatrix}^T \\
\mathbf{A}= \begin{pmatrix}
-k_1 & 0 & \beta_H b_s & 0\\
(1-\alpha)\gamma_H & -k_2 & 0 & 0\\
\beta_S b_S \frac{S_S}{N_H} & \beta_S b_S \frac{S_S}{N_H}
& -\mu_S & \beta_S b_{SR} \frac{S_S}{N_R}\\
0 & 0 & \beta_S b_{SR} \frac{S_R + (1-\epsilon)V_R}{N_R} & -k_3\\
\end{pmatrix}
\\
%\label{Ghat}
\hat{G}(X,Z) =\begin{pmatrix}
\beta_H b_S I_S \bigl(1 - \frac{S_H}{N_H}\bigr)\\
0 \\
\beta_S b_S S_S \bigl(1 - \frac{I_H + P_H}{N_H}\bigr)
+ \beta_S b_{SR} S_H \bigl(1 - \frac{I_R}{N_R}\bigr)\\
\beta_R b_{SR} \bigl(1 - \frac{S_R + (1-\epsilon)V_R}{N_R}\bigr)
\end{pmatrix}
\end{gather*}
it is clear that A is an M matrix and $\hat{G}(X,Z) > 0$, and hence we have the
following lemma.
\begin{lemma} \label{lem2}
The disease-free equilibrium is globally asymptotically
stable when $\mathcal{R}_0 < 1$.
\end{lemma}
\section{Sensitivity analysis} \label{Sensitivity}
The impact of leishmaniasis model \eqref{eL} significant parameters are determined
via sensitivity analysis \cite{Blower1994, Marino2008, McLeod2006}
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
\cite{Blower1994, Marino2008, Mckay1979,Sanchez1997}.
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 \cite{Blower1994, Marino2008, Mckay1979,Sanchez1997}.
A total of 1,000 simulations of the model \eqref{eL} \emph{per} LHS run were carried out,
using the ranges and baseline values tabulated in Table \ref{t2}
(with the reproduction number, $\mathcal{R}_0$, as the response function).
The values of the parameters are taken from published literature as mentioned in
Table \ref{t2} and for each parameter value, an interval within 20\% range of the
parameter value is formed to test the sensitivity.
Figure \ref{PRCC1} depicts the PRCC values for each parameter of the models using
the reproduction number, $\mathcal{R}_0$, as the response function. Parameters with
the highest PRCC values have the largest impact on $\mathcal{R}_0$.
Bars extending to the left (for negative PRCC values) or to the right
(for positive PRCC values). The parameters $\gamma_H$, $\alpha$, $\sigma_H$, $\nu_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, $\mathcal{R}_0$. Therefore, the key parameters
influencing $\mathcal{R}_0$ are separated into those that decrease $\mathcal{R}_0$
when increased (those with significantly negative PRCC values) and those that causes
$\mathcal{R}_0$ to increase when increased (those with significantly positive PRCC values).
\begin{figure}[htb]
\begin{center}
\includegraphics[width= 0.9\textwidth]{fig2} % PRCCR0.eps
\end{center}
\caption{\small PRCC values for the leishmaniasis model \eqref{eL}, using the reproduction
number ($\mathcal{R}_0$) as the response function. Parameter values (baseline) and
ranges used are as given in Table \ref{t2}.}
\label{PRCC1}
\end{figure}
The identification of these key parameters with significant impact on the reproduction
number $\mathcal{R}_0$ 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 \emph{per} contact for susceptible canine population
(reduce $\beta_R$), this should be followed by decreasing the transmission probability
for susceptible sandflies (reduce $\beta_S$), the recruitment rate of sandflies
(reduce $\Lambda_S$), the sandflies biting rate on humans (reduce $b_S$), the sandflies
biting rate on the canine population (reduce $b_{SR})$, and the transmission probability
for susceptible humans (reduce $\beta_H$).
The result from the sensitivity analysis further suggests control strategy that increases
the natural death rate of sand flies (increase $\mu_S$), the vaccination efficacy in canine
reservoir (increase $\varepsilon$), and the disease-induced death rate of infected canine
population (reduce $\delta_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 \eqref{eL}
shows that the significant parameters are the natural death rate of sandflies ($\mu_S$),
sandflies biting rate on the canine population ($b_{SR})$, transmission probability for
susceptible sandflies $(\beta_S)$, recruitment rate of sandflies $(\Lambda_S)$,
sandflies biting rate in humans $(b_S)$, transmission probability for susceptible
canine population $(\beta_R)$, transmission probability for susceptible humans $(\beta_H)$,
vaccination efficacy in canine population $(\varepsilon)$, recruitment rate of canine
population $(\Lambda_R)$, and the disease-induced death rate of infected canine population
$(\delta_R$).
\begin{table}[htb]
\caption{Parameters values of the leishmaniasis model \eqref{eL}.}
\label{t2}
\centering
\begin{tabular}{llll} \hline
Param. & Values & Range & References \\\hline
$\Lambda_H$ & $4.16\times10^{-5}$ & $(1-0.2)4.16\times 10^{-5}$ - $(1+0.2)4.16\times 10^{-5}$ & \cite{WB2015} \\
$\Lambda_S$ & $0.155$ & $0.124$ - $0.186$ & \cite{Kasap2006} \\
$\Lambda_R$ &$0.0027$ & $0.00216$ - $0.00324$ & assumed \\
$\beta_H$ & $0.5$ & $0.4$ - $0.6$ & \cite{Stauch2014} \\
$\beta_S$ & $0.22$ & $0.176$ - $0.264$ & \cite{Elmojtaba2010} \\
$\beta_R$ & $0.0714$ & $0.05712$ - $0.08568$ & \cite{Elmojtaba2010} \\
$b_S$ & $0.2856$ & $0.22848$ - $0.34272$ & \cite{Elmojtaba2010} \\
$b_{SR}$ & $0.56$ & $0.448$ - $0.672$ & assumed \\
$\mu_H$ & $1.64\times10^{-5}$ & $1.312\times10^{-5}$ - $2.952\times 10^{-5}$ & \cite{WB2015} \\
$\mu_S$ & $0.056$ & $0.0448$ - $0.0672$ & \cite{Kasap2006} \\
$\mu_R$ & $2.11\times10^{-4}$ & $1.688\times10^{-4}$ - $2.532\times10^{-4}$ & assumed \\
$\delta_H$ & $0.0014$ & $0.00112$ - $0.00168$ & \cite{standford2006} \\
$\delta_R$ & $0.0014$ & $0.00112$ - $0.00168$ & \cite{{Beneth2008}} \\
$\gamma_H$ & $0.5$ & $0.4$ - $0.6$ &assumed \\
$\alpha$ & $0.40$ & $0.32$ - $0.48$ & \cite{Gasim2000} \\
$\tau_H$ & $0.033$ & $0.02643$ - $0.0396$ & \cite{Elmojtaba2010} \\
$\sigma_H$ & $0.00556$ & $0.00448$ - $0.006672$ & \cite{Elmojtaba2010} \\
$\nu_R$ & $0.5$ & $0.4$ - $0.6$ & assumed \\
$\varepsilon$ & $0.8$ & $0.64$ - $0.96$ & \cite{Diniz2008} \\
\hline
\end{tabular}
\end{table}
\section{Optimal control problem} \label{Control}
Following the results obtained from the sensitivity analysis, we introduce three
time dependent controls into the system \eqref{eL}. Let $u_1(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 $u_2(t)$ be the time dependent indoor insecticide spraying rate (insecticide such
as pyrethroids are effective against sandflies both indoors and outdoors
\cite{Alexander2003}). Lastly, let $u_3(t)$ be the time dependent culling rate of the
infected canine population \cite{Singh2013}. Then the above system of ODEs \eqref{eL}
becomes
\begin{align}
S_H' &= \Lambda _H - \frac{\beta_Hb_S[1-u_1(t)]I_SS_H}{N_H} - \mu_H S_H \nonumber \\
I_H' &= \frac{\beta_Hb_S [1-u_1(t)]I_SS_H}{N_H}- (\gamma_H + \delta_H + \mu_H) I_H \nonumber \\
P_H' &= (1-\alpha) \gamma_H I_H - (\tau_H + \sigma_H + \mu_H) P_H \nonumber \\
R_H' &= \alpha\gamma_H I_H + (\tau_H + \sigma_H) P_H - \mu_H R_H \nonumber \\
S_S' &= \Lambda_S- \frac{\beta_Sb_S[1-u_1(t)](I_H+P_H) S_S}{N_H}
- \frac{\beta_Sb_{SR} I_R S_S}{N_R} - [\mu_S+u_2(t)]S_S \nonumber \\
I_S' &= \frac{\beta_Sb_S[1-u_1(t)](I_H+P_H) S_S}{N_H} + \frac{\beta_Sb_{SR}I_R S_S}{N_R}
- [\mu_S + u_2(t)]I_S \nonumber \\
S_R' &= \Lambda_R - \frac{\beta_Rb_{SR}I_SS_R}{N_R} - [\mu_R + \nu_R]S_R \nonumber \\
V_R' &= \nu_R S_R - \frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R} - \mu_R V_R \nonumber \\
I_R' &= \frac{\beta_Rb_{SR}I_SS_R}{N_R} + \frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R}
- [\mu_R + \delta_R + u_3(t)]I_R. \label{eCL}
\end{align}
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(u_1,u_2,u_3)$ where
\begin{equation} \label{eJ}
\begin{aligned}
J(u_1,u_2,u_3)&=\int_0^{T}[A_1I_H+A_2P_H+A_3I_R+B_1u_2S_S+B_2u_2I_S+B_3u_3I_R \\
&\quad + C_1u_1^2+C_2u_2^2+C_3u_3^2]dt.
\end{aligned}
\end{equation}
subject to the differential equations \eqref{eCL}, where $T$ is the final time.
This performance specification involves the total infected humans (with VL and PKDL),
infected reservoir (canine), along with the cost of applying the controls $u_1(t),u_2(t)$
and $u_3(t)$). In this paper, a quadratic objective functional is implemented for
measuring the control cost, such a cost has been frequently used
\cite{Joshi2002, Jung2002, Kern2007, Kirschner1997, Lenhart2007, Xiefei2007}.
The positive coefficients, $A_i,B_i,C_i; i=1,\cdots,3$ are balancing weight parameters.
The controls, $u_1(t),u_2(t)$ and $u_3(t)$, are bounded, Lebesgue integrable functions
\cite{Jung2002,Xiefei2008}. And we seek to find optimal controls, $u_1^*,u_2^*$ and $u^*_3$,
such that
\begin{equation}
J(u_1^*,u_2^*,u_3^*)=\min_{(u_1,u_2,u_3)\in\mathcal U}\{J(u_1,u_2,u_3)\}
\end{equation}
where the admissible set is
\[
\mathcal{U}=\{(u_1,u_2,u_3)\in (L^\infty(0,T))^3:
0\le u_i\le M_i; M_i\in \mathbb R^+,i=1,2,3 \}.
\]
The following theorem proves the existence of the solution of the system
\eqref{eCL} as well as the non negativity and boundedness of the state variables.
\begin{theorem} \label{thm2}
Given controls $u=(u_1,u_2,u_3)\in \mathcal{U}$, there exist non-negative bounded solutions
$(S_H,I_H,P_H,$ $R_H,S_S,I_S,S_R,I_R)$ to the state system \eqref{eCL} in the finite
interval $[0,T]$ with given initial conditions.
\end{theorem}
The existence and uniqueness of solutions for the state system \eqref{eCL} with a given
control pair can be proven by using a result from Lukes \cite{Lukes1982}.
The structure of system \eqref{eCL} gives the non-negativity and uniform boundedness of
the state solutions. The next theorem proves the existence of the optimal controls.
\begin{theorem} \label{thm3}
There exists an optimal control tuple $u^*=(u_1^*,u_2^*,u_3^*)\in \mathcal{U}$ with
corresponding 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(u_1,u_2,u_3)$.
\end{theorem}
\begin{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 $(u_1^n,u_2^n, u_3^n)$ such that
$$
\lim_{n\to \infty}J(u_1^n,u_2^n,u_3^n)=\inf_{(u_1,u_2,u_3)\in\mathcal U}J(u_1,u_2,u_3).
$$
Let us denote
\begin{align*}
&(S_{H}^n,I_{H}^n,P_{H}^n,R_{H}^n,S_{S}^n,I_{S}^n,S_{R}^n,V_{R}^n,I_{R}^n)\\
&=(S_H,I_H,P_H,R_H,S_S,I_S,S_R,V_R,I_R)(u_i^n,u_2^n,u_3^n).
\end{align*}
Since all of the state variables are bounded (Theorem \ref{thm2}), then their first derivatives
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
constant. Thus the sequence
$(S_{H}^n,I_{H}^n,P_{H}^n,R_{H}^n,S_{S}^n,I_{S}^n,S_{R}^n,V_{R}^n,I_{R}^n)$
is uniformly equicontinuous in $[0,T]$. Then by the Arzela-Ascoli Theorem \cite{Lukes1982},
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 $(u_1^n,u_2^n,u_3^n)$ has a subsequence that converges weakly in
$L^2(0,T)$. Let $(u_1^*,u_2^*,u_3^*)\in U$ be such that $u_i^n\rightharpoonup u_i^*$
weakly in $L^2(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 $B_1u_2^nS_S^n$.
Using lower semi-continuity of norms in weak $L^2$, we obtain
$$
\| u_i^*\|^2_{L^2}\leq \liminf_{n\to \infty} \|u_i^n\|^2_{L^2}
\quad \text{for } i=1,2,3.
$$
Hence,
\begin{align*}
J(u_1^*,u_2^*,u_3^*)
&\leq \liminf_{n\to\infty} \int_0^T\big[A_1I_H^n+A_2P_H^n+A_3I_R^n+B_1u_2^nS_S^n
+B_2u_2^nI_S^n \\
&\quad +B_3u_3^nI_R^n + C_1(u_1^n)^2+C_2(u_2^n)^2+C_3(u_3^n)^2\big]dt\\
&=\liminf_{n\to\infty} J(u_1^n,u_2^n,u_3^n)
\end{align*}
Thus, there exists a vector $\vec u=(u_1^*,u_2^*,u_3^*)$ of controls that minimizes
the objective functional $J(u_1,u_2,u_3)$.
\end{proof}
Next, we characterize the optimal control triple ($u_1(t),u_2(t),u_3(t)$).
Pontryagin's Maximum Principle \cite{Pontryagin1962} introduces adjoint functions
that allow the state system \eqref{eCL} to be attached to the objective functional \eqref{eJ}.
\subsection*{Characterization of optimal controls}\label{Con}
The necessary conditions that an optimal control must satisfy come from the Pontryagin's
Maximum Principle \cite{Pontryagin1962}. This principle converts \eqref{eCL}
and \eqref{eJ} into a problem of minimizing pointwise a Hamiltonian $H$, with respect
to the controls $u_1,u_2,u_3\in \mathcal U$. First we formulate the Hamiltonian
\cite{Lenhart2007} from the cost functional \eqref{eJ} and the governing dynamics
\eqref{eCL} to obtain the optimality conditions.
\begin{align*}
H&= A_1I_H+A_2P_H+A_3I_R+B_1u_2S_S+B_2u_2I_S+B_3u_3I_R + C_1u_1^2+C_2u_2^2+C_3u_3^2\\
&\quad + \lambda_{S_H}\Big[\Lambda _H - \frac{\beta_Hb_S[1-u_1(t)]I_SS_H}{N_H}
- \mu_H S_H \Big]\\
&\quad + \lambda_{I_H}\Big[\frac{\beta_Hb_S [1-u_1(t)]I_SS_H}{N_H}- (\gamma_H
+ \delta_H + \mu_H) I_H\Big]\\
&\quad + \lambda_{P_H}\Big[(1-\alpha) \gamma_H I_H - (\tau_H + \sigma_H + \mu_H) P_H\Big]\\
&\quad + \lambda_{R_H}\Big[\alpha\gamma_H I_H + (\tau_H + \sigma_H) P_H - \mu_H R_H\Big]\\
&\quad + \lambda_{S_S}\Big[\Lambda _S - \frac{\beta_Sb_S[1-u_1(t)](I_H+P_H) S_S}{N_H}
- \frac{\beta_Sb_{SR} I_R S_S}{N_R} - [\mu_S+u_2(t)]S_S\Big]\\
&\quad + \lambda_{I_S}\Big[ \frac{\beta_Sb_S[1-u_1(t)](I_H+P_H) S_S}{N_H}
+ \frac{\beta_Sb_{SR}I_R S_S}{N_R} - [\mu_S + u_2(t)]I_S\Big]\\
&\quad + \lambda_{S_R}\Big[\Lambda_R-\frac{\beta_Rb_{SR}I_SS_R}{N_R} - [\mu_R + \nu_R]S_R\Big]\\
&\quad + \lambda_{V_R}\Big[\nu_R S_R - \frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R}
- \mu_R V_R\Big]\\
&\quad + \lambda_{I_R}\Big[\frac{\beta_Rb_{SR}I_SS_R}{N_R}
+ \frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R}- [\mu_R + \delta_R + u_3(t)]I_R\Big].
\end{align*}
Using the system of adjoint variables,
\begin{equation} \label{eA}
\begin{aligned}
\lambda'_{S_H}&=-\frac{\partial H}{\partial S_H} \\
&=(\lambda_{S_H}-\lambda_{I_H})\beta_H b_S[1-u_1(t)]I_S\frac{(I_H+P_H+R_H)}{N_H^2}
+\mu_H\lambda_{S_H}\\
&\quad+(\lambda_{I_S}-\lambda_{S_S})\beta_S b_S[1-u_1(t)]S_S\frac{(I_H+P_H)}{N_H^2},
\end{aligned}
\end{equation}
\begin{align*}
\lambda'_{I_H}&=-\frac{\partial H}{\partial I_H}\\
&=-A_1+ (\lambda_{I_H}-\lambda_{S_H})\frac{\beta_Hb_S[1-u_1(t)]I_SS_H}{N_H^2}-\lambda_{P_H}
(1-\alpha)\gamma_H-\lambda_{R_H}\alpha\gamma_H \\
&\quad +\lambda_{I_H}(\gamma_H+\delta_H+\mu_H)+(\lambda_{S_S}
-\lambda_{I_S})\beta_S b_S[1-u_1(t)]S_S\frac{[S_H+R_H]}{N_H^2},
\end{align*}
\begin{align*}
\lambda'_{P_H}
&=-\frac{\partial H}{\partial P_H}\\
&=-A_2+(\lambda_{I_H}-\lambda_{S_H})\frac{\beta_H b_S[1-u_1(t)]I_SS_H}{N_H^2}
+ \lambda_{P_H}(\tau_H+\sigma_H+\mu_H) \\
&\quad -\lambda_{R_H}(\tau_H+\sigma_H)+(\lambda_{S_S}
-\lambda_{I_S})\beta_S b_S[1-u_1(t)]S_S\frac{[S_H+R_H]}{N_H^2},
\end{align*}
\begin{align*}
\lambda'_{R_H}
&=-\frac{\partial H}{\partial R_H}\\
&= (\lambda_{I_H}-\lambda_{S_H})\beta_H b_S[1-u_1(t)]I_S\frac{S_H}{N_H^2}+\lambda_{R_H}\mu_H \\
&\quad +(\lambda_{I_S}-\lambda_{S_S})\beta_H b_S[1-u_1(t)]S_S\frac{(I_H+P_H)}{N_H^2},
\end{align*}
\begin{align*}
\lambda'_{S_S}
&=-\frac{\partial H}{\partial S_S}\\
&=(\lambda_{S_S}-\lambda_{I_S})\Big[\beta_S b_S[1-u_1(t)]\frac{(I_H+P_H)}{N_H}
+\frac{\beta_Sb_RI_R}{N_R}\Big] \\
&\quad -B_1u_2(t)+\lambda_{S_S}[\mu_S+u_2(t)],
\end{align*}
\begin{align*}
\lambda'_{I_S}
&=-\frac{\partial H}{\partial I_S} \\
&=-B_2u_2(t)+(\lambda_{S_H}-\lambda_{I_H})\frac{\beta_Hb_S[1-u_1(t)]S_H}{N_H}
+\lambda_{I_S}[\mu_S+u_2(t)] \\
&\quad +(\lambda_{S_R}-\lambda_{I_R})\frac{\beta_Rb_{SR}S_R}{N_R}
+(\lambda_{V_R}-\lambda_{I_R})\frac{\beta_Rb_{SR}(1-\varepsilon)V_R}{N_R},
\end{align*}
\[
\begin{aligned}
\lambda'_{S_R}
& =-\frac{\partial H}{\partial S_R} \\
&=(\lambda_{I_S}-\lambda_{S_S})\frac{\beta_Sb_{SR}I_RS_S}{N_R^2} \\
&\quad + ( \lambda_{S_R}-\lambda_{I_R})\frac{\beta_Rb_{SR}I_S(I_R+V_R)}{N_R^2}
+\lambda_{I_R}(\mu_R+\nu_R) \\
&\quad -\lambda_{V_R}\nu_R+(\lambda_{I_R}
-\lambda_{V_R})\frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R^2},
\end{aligned}
\]
\begin{align*}
\lambda'_{V_R}
&=-\frac{\partial H}{\partial V_R} \\
&= (\lambda_{I_S}-\lambda_{S_S})\frac{\beta_Sb_{SR}I_RS_S}{N_R^2}
+(\lambda_{I_R}-\lambda_{S_R})\frac{\beta_Rb_{SR}I_SS_R}{N_R^2} \\
&\quad +(\lambda_{V_R}-\lambda_{I_R})\beta_Rb_{SR}(1-\varepsilon)I_S\frac{(S_R+I_R)}{N_R^2}
+\lambda_{V_R}\mu_R,
\end{align*}
\begin{align*}
\lambda'_{I_R}
&=-\frac{\partial H}{\partial I_R} \\
&= -A_3-B_3u_3(t)+(\lambda_{S_S}-\lambda_{I_S})\frac{\beta_Sb_{SR}S_S(S_R+V_R)}{N_R^2} \\
&\quad +(\lambda_{I_R}-\lambda_{S_R})\frac{\beta_Rb_{SR}I_SS_R}{N_R^2}
+(\lambda_{I_R}-\lambda_{V_R})\frac{\beta_Rb_{SR}(1-\varepsilon)I_SV_R}{N_R^2} \\
&\quad +\lambda_{I_R}(\mu_R+\delta_R+u_3),
\end{align*}
with the transversality condition
\begin{equation}
\lambda_i(t_f)=0\quad \text{for }i=S_H,I_H,P_H,R_H,S_S,I_S,S_R,V_R,I_R. \label{eC}
\end{equation}
The optimality conditions are given as
$$
\frac{\partial H}{\partial u_j}=0,\quad j=1,2, 3.
$$
in the interior of the control set $ \mathcal U$.
Thus, we have
\begin{gather*}
\begin{aligned}
&2C_1u_1 + \lambda^*_{S_H}\frac{\beta_Hb_SI^*_SS^*_H}{N^*_H}-\lambda^*_{I_H}
\frac{\beta_Hb_SI^*_SS^*_H}{N^*_H} \\
&+\lambda^*_{S_S}\frac{\beta_Sb_S(I^*_H+P^*_H)S^*_S}{N^*_H}
-\lambda^*_{I_S}\frac{\beta_sb_S(I^*_H+P^*_H)S^*_S}{N^*_H}=0,
\end{aligned} \\
B_1S^*_S + B_2I^*_S+2C_2u_2-\lambda^*_{S_S}S^*_S-\lambda^*_{I_S}I^*_S=0,\\
B_3I^*_R + 2C_3u_3-\lambda^*_{I_R}I^*_R=0.
\end{gather*}
Hence, the control characterization is
\begin{gather*}
\begin{aligned}
u_1^*= \min\Big[&1,\max\Big(\frac{1}{2C_1N^*_H}[\beta_Hb_SI^*_SS^*_H(\lambda^*_{I_H}
-\lambda^*_{S_H}) \\
& +\beta_Sb_SS^*_S(I^*_H+P^*_H)(\lambda^*_{I_S}
-\lambda^*_{S_S})],0\Big)\Big],
\end{aligned}\\
u_2^*= \min\Big[1,\max\Big(\frac{1}{2C_2}[-B_1S^*_S-B_2I^*_S
+\lambda^*_{S_S}S_S+\lambda^*_{I_S}I^*_S], 0\Big)\Big],\\
u_3^*= \min\Big[1,\max\Big(\frac{1}{2C_3}[-B_3I^*_R+\lambda^*_{I_R}I^*_R], 0\Big)\Big].
\end{gather*}
\subsection*{Setting bounds on the control}
We set the bounds on the time dependent control, $u_1(t)$, for the personal protection
rate between 0 and 1 (i.e, $0\le u_1(t)\le 1$). The value of $u_1(t)=1$ means highly
effective control with no contact between human and sand flies while the value of $u_1(t)$
close to 0 means high contact between human and sand-flies. For the the control $u_2(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 $u_2(t)$ is 0.13 i.e, $0\le u_2(t)\le 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 $u_3(t)$,
the lower bound for $u_3(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 $u_3(t)$ is 0.1 i.e, $0\le u_3(t)\le 0.1$.
Next, we discuss the numerical solutions of the optimality system, the corresponding optimal
control and the interpretations from various cases.
\section{Numerical results}\label{Numerical}
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 \eqref{eCL}. The results obtained for the state
variables are substituted into the adjoint equations \eqref{eA}. These adjoint equations
with given final conditions \eqref{eC} 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 \cite{Lenhart2007}.
To illustrate the optimal control strategies, the following initial conditions were used:
$S_H(0) = 1000$, $I_H(0) = 25$, $P_H(0) = 25$, $R_H(0) = 14$, $S_S(0) = 500$,
$I_S(0) = 152$, $S_R(0) = 224$, $V_R(0) = 125$, $I_R(0) = 114$.
Furthermore, the values $A_1=A_2=A_3=1.00$, $B_1=B_2=B_3=0.01$,
$C_1=C_2=C_3=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 \ref{t2}, the reproduction number is obtained as
$\mathcal{R}_0=2.33>1$, thus indicating that the disease is endemic in the population.
\begin{figure}[htb]
\begin{center}
\includegraphics[width= 0.96\textwidth]{fig3} % InfectedHuman.png
\end{center}
\caption{\small Simulation of the leishmaniasis model \eqref{eCL} 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}
\label{f1}
\end{figure}
The results of the optimal control simulations of the leishmaniasis model \eqref{eCL}
are depicted in Figures \ref{f1}, 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 \ref{f1}(a), we observed that the total number of infected humans
($I_H$ and $P_H$) 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 \ref{f1}(b) and \ref{f1}(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
($u_{1}, u_{2}$, $u_{3}$) are depicted in Figures \ref{f2}(a)-\ref{f2}(c). The time
dependent control $u_1(t)$ is observed to be at its upper bound for about 2 months and
then can be reduced to get the optimal results. Similarly, the control $u_2(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 $u_3(t)$ is at
the upper bound for about four and half months and can be reduced gradually.
The total number of infectious humans ($I_H+P_H$) 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\times10^5$ and
with optimal control $J$ is about 70\% less which is $1.0701\times 10^5$.
\begin{figure}[htb]
\begin{center}
\includegraphics[width= 0.96\textwidth]{fig4} % u110.png
\end{center}
\caption{\small Optimal controls of the leishmaniasis model \eqref{eCL} for:
(a) Personal protection control; (b) Sandflies insecticide control;
(c) Reservoir culling control.}
\label{f2}
\end{figure}
\subsection*{Control within varying values of parameters
$\Lambda_R$, $\Lambda_S$, $b_{SR}$ and $\tau_H$}
To carry out the sensitivity analysis test in Section 3, we set the values of the
parameters to be within $\pm 20\%$ range of their baseline values and found the parameters
$\tau_H$, $\Lambda_R$, $\Lambda_S$ and $b_{SR}$ among others to have a high impact on
the reproduction number $\mathcal{R}_0$. However, the time dependent controls
$u_1(t)$, $u_2(t)$ and $u_3(t)$ incorporated into system \eqref{eCL} did not directly
act on these parameters. Hence, we vary these four parameters using a range of $\pm~50\%$
from their baseline values and determine their impact on the system as well as the control
variables. We have used the $\pm~50\%$ range for emphasis purpose, the result is expected
to be similar in the $\pm~20\%$ range. All other parameters in the model \ref{eCL} and the
balancing parameters in the objective functional are taken to be in the baseline values
from Table \ref{t2}.
The results of the optimal control simulations of the leishmaniasis model \eqref{eCL}
in the uncertainty interval for these parameters $\Lambda_R,\Lambda_S$ and $b_{SR}$
are depicted in Figures \ref{s1}(a)- ref{s1}(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.
\begin{figure}[htb]
\begin{center}
\includegraphics[width= 0.96\textwidth]{fig5} % InfectedHumanSensitivity.png
\end{center}
\caption{\small Simulation of the Leishmaniasis model \eqref{eCL} 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 $\Lambda_R,\Lambda_S$ and $b_{SR}$ are set within $\pm 50\%$
interval of their baseline values given in Table \ref{t2} 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.}
\label{s1}
\end{figure}
From Figure \ref{s1},
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 \ref{s4}(a)-\ref{s4}(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
these parameters ($\Lambda_R,\Lambda_S$ and $b_{SR}$),
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 $\pm20\%$
uncertainty interval used for the sensitivity analysis.
\begin{figure}[htb]
\begin{center}
\includegraphics[width= 0.96\textwidth]{fig6} % u110Sensitivity.png
\end{center}
\caption{\small Optimal controls of the leishmaniasis model \eqref{eCL} in the $\pm 50\%$
uncertainty interval of parameters $\Lambda_R,\Lambda_S$ and $b_{SR}$ for:
(a) Personal protection control; (b) Sandflies insecticide control;
(c) Reservoir culling control. The parameters $\Lambda_R,\Lambda_S$ and $b_{SR}$ are
set within $\pm 50\%$ interval of their baseline values given in Table \ref{t2} 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.}
\label{s4}
\end{figure}
Next, we investigate the optimal control of system \eqref{eCL} within the uncertainty
interval of parameter $\tau_H$; this parameter is the human recovery rate and we use the
$\pm 50\%$ range from the baseline value. The results are depicted in Figures \ref{s7},
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 \ref{s7}(a) that more infected humans
at the lower bound and fewer at the upper bound of $\tau_H$.
Similarly for the optimal controls $u_1(t),~u_2(t)$ and $u_3(t)$ in Figures
\ref{s7}(b)-\ref{s7}(d). More control efforts are required when humans recover
slowly, particularly in the use of insecticide spray rate ($u_1(t)$) and dog culling
rate ($u_3(t)$).
\begin{figure}[htb]
\begin{center}
\includegraphics[width= 0.96\textwidth]{fig7} % InfectedHumanSensitivity_tauH.png
\end{center}
\caption{\small Optimal controls of the leishmaniasis model \eqref{eCL} in the $\pm 50\%$
uncertainty interval of $\tau_{H}$ for:
(a) Total number of infected humans;
(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 \ref{t2}. 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.}
\label{s7}
\end{figure}
In summary, numerical simulations of the leishmaniasis control model \eqref{eCL}
show that leishmaniasis can be reduced in the community by the application of
the time dependent control triple ($u_{1}(t), u_{2}(t)$ and $u_{3}(t)$) which represent
respectively the of personal protection($u_1(t)$) such as insect repellant, door
and window screens against sand flies, insecticide spraying ($u_2(t)$) and culling of
infected canine reservoir ($u_3(t)$).
\section{Conclusions and discussion} \label{Conclusion}
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
(${\mathcal R}_0$, 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 impact are the
natural death rate of sandflies ($\mu_S$), sandflies biting rate on the canine reservoir
($b_{SR})$, transmission probability \emph{per} contact for susceptible sandflies
$(\beta_S)$, recruitment rate of sandflies $(\Lambda_S)$, sandflies biting rate in humans
$(b_S)$, transmission probability \emph{per} contact for susceptible canine reservoir
$(\beta_R)$, transmission probability \emph{per} contact for susceptible humans $(\beta_H)$,
vaccination efficacy in canine reservoir $(\varepsilon)$, recruitment rate of canine reservoir
$(\Lambda_R)$, and the disease-induced death rate of infected reservoir $(\delta_R$).
The identification of these key parameters is vital to the formulation of effective control
strategies for combating the spread of the disease. Thus, following the results obtained
from the sensitivity analysis, we introduce three time dependent controls into system \eqref{eL}. 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 simulations 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 differential
equations for leishmaniasis in humans and canine reservoirs. Some of theoretical
and epidemiological findings of this study are summarized below:
\begin{itemize}
\item[(i)] The leishmaniasis model \eqref{eL} is locally-and globally-asymptotically
stable (LAS) when $\mathcal{R}_0 < 1$ and unstable when $\mathcal{R}_0 > 1$;
\item[(ii)] The sensitivity analysis of the model shows that the significant parameters
are the natural death rate of sandflies ($\mu_S$), sandflies biting rate on the canine
reservoir ($b_{SR})$, transmission probability \emph{per} contact for susceptible
sandflies $(\beta_S)$, recruitment rate of sandflies $(\Lambda_S)$, sandflies biting rate
in humans $(b_S)$, transmission probability \emph{per} contact for susceptible canine
reservoir $(\beta_R)$, transmission probability \emph{per} contact for susceptible humans
$(\beta_H)$, vaccination efficacy in canine reservoir $(\varepsilon)$, recruitment rate of
canine reservoir $(\Lambda_R)$, and the disease-induced death rate of infected reservoir
$(\delta_R$)
\item[(iii)] Numerical simulations indicates that that leishmaniasis can be controlled in
the community by the application of time dependent controls representing personal protection
($u_1(t)$), insecticide spraying ($u_2(t)$) and culling infected canine reservoir ($u_3(t)$)
respectively.
\item[(iv)] The control profiles were approximately the same for all parameter values chosen
within the uncertainty interval of these parameters ($\Lambda_R,\Lambda_S$ and $b_{SR}$).
\item[(v)] In the uncertainty interval for parameter $\tau_H$, more control efforts are
required when humans recover more slowly.
\end{itemize}
\subsection*{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 institute sponsored by the National
Science Foundation, the U.S. Department of Homeland 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.
\begin{thebibliography}{00}
\bibitem{Alexander2003} Alexander, B.; Maroli, M.;
Control of phlebotomine sandflies. \emph{Med Vet Entomol}, \textbf{17} (2003)(1): 1--18.
\bibitem{Beneth2008} Beneth, G.; Koutinas, A. F.; Laia, S.G.; Bourdeau, P.; Ferrer, L.;
Canine leishmaniosis - new concepts and insights on an expanding zoonosis:
part one \emph{Trends in Parasitology}\textbf{24} (2008) (7).
\bibitem{Blower1994} Blower, S. M.; Dowlatabadi, H.;
Sensitivity and uncertainty analysis of complex models of disease transmission:
an HIV model, as an example. \emph{Int. Stat. Rev.}, \textbf{2} (1994), 229--243.
\bibitem{Bortoletto2016} 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).
\bibitem{Bradley2012} Bradley, T.; King, R.;
The dog economy is global but what is the world's true canine capital?
\emph{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 December 30, 2015.
\bibitem{Carlos1} Castillo-Chavez, C.; Feng, Z.; Huang, W.;
On the computation of $R_0$ and its role on global stability. In
\emph{Mathematical approaches for emerging and re emerging infectious diseases.
An Introduction}. Castillo-Chavez, C., Blower, S., van den Driessche, P., Kirschner,
D. and Yakubu, A.-A. (Eds.) Springer-Verlag, New York (2002), 229-250.
\bibitem{CDC2016} Centers for Disease Control and Prevention (2013).
Leishmaniasis.\\
http://www.cdc.gov/parasites/leishmaniasis/ (2002). Accessed January 4, 2016.
\bibitem{Chaves2004} Chaves, L. F.; Hernandez, M. J.;
Mathematical modelling of American Cutaneous Leishmaniasis:
incidental hosts and threshold conditions for infection persistence.
\emph{ACTA TROPICA (Elsevier)}. \textbf{92} (2004), 255--252.
\bibitem{Costa2013} 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. \emph{PLOS: Neglected Tropical Disease }. \textbf{7} (8) (2013),
1--8.
\bibitem{Deer1999} Deer, H. M.;
Pesticide adsorption and half-life. AG/Pesticides, 15 (1999), 1.
\bibitem{Diekmann1990} Diekmann, O.; Heesterbeek, J. A. P.; Metz, J. A. J.;
On the definition and computation of the basic reproduction ratio $R_0$ in
models for infectious diseases in heterogeneous populations. \emph{J.
Math. Biol.}, \textbf{28} (1990), 365-382.
\bibitem{Diniz2008} 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.
\emph{The Journal of Infection in Developing Countries}. \textbf{2}(1) (2008), 24--33.
\bibitem{Elmojtaba2010} Elmojtaba, I. M.; Mugisha, J. Y. T.; Hashim, M. H. A.;
Mathematical analysis of the dynamics of visceral leishmaniasis in Sudan.
\emph{Applied Mathematics and Computation} \textbf{217}(6) (2010), 2567-2578.
\bibitem{Ferrell2003} 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.
\bibitem{Gasim2000} 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
acquisition of Leishmania reactivity by peripheral blood mononuclear cells (PBMC).
\emph{Clinical and experimental Immunology}, \textbf{119} (3) (2000), 523-529.
\bibitem{Hyman2000} Hyman, J. M.; Li, J.;
An intuitive formulation for the reproductive number for the spread of diseases in
heterogeneous populations. \emph{Mathematical Biosciences}, \textbf{167} (2000), 65-86.
\bibitem{Joshi2002} Joshi, H. R.;
Optimal Control of an HIV Immunology Model, \emph{Optim. Control Appl. Math},
\textbf{23} (2002), 199--213.
\bibitem{Jung2002} Jung, E.; Lenhart, S.; Feng, Z.;
Optimal Control of Treatments in a Two-Strain Tuberculosis Model,
\emph{Discrete and Continuous Dynamical Systems-Series B}. \textbf{2}(4) (2002), 473--482.
\bibitem{Kamgang2008} Kamgang, J. C.; Sallet, G.;
Computation of threshold conditions for epidemiological models and global stability of
the disease-free equilibrium (DFE). \emph{Mathematical Biosciences}, \textbf{213}
(2008), 1-12.
\bibitem{Kasap2006} Kasap, O. E.; Alten, B.;
Comparative demography of the sand fly Phlebotomus papatasi (Diptera: Psychodidae)
at constant temperatures. \emph{Journal of Vector Ecology}, \textbf{31}(2) (2006), 378-385.
\bibitem{Kern2007} Kern, D.; Lenhart, S.; Miller, R.; Yong, J.;
Optimal Control Applied to Native-Invasive Population Dynamics.
\emph{J Biol Dyn.} {\bf 1} (4) (2007), 413--26. doi: 10.1080/17513750701605556
\bibitem{Kirschner1997} Kirschner, D.; Lenhart, S.; Serbin, S.;
Optimal Control of the Chemotherapy of HIV, \emph{J. Math. Biol.} 35, (1997), 775--792.
\bibitem{Kumar2014} Kumar, R.; Engwerda, C.;
Vaccines to prevent leishmaniasis. \emph{Clinical \& Translational Immunology }.
\textbf{3}(e13) (2014), 1--6.
\bibitem{Lemesre2007} 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.
\emph{Vaccine}. \textbf{25}(21) (2007), 4223--4234.
\bibitem{Lenhart2007} Lenhart, S.; Workman, J. T.;
Optimal Control Applied to Biological Models.
\emph{Chapman and Hall/CRC}. First Edition, 1998.
\bibitem{Lukes1982} Lukes, D. L.;
Differential Equations: Classical to Controlled.
\emph{Mathematics in Science and Engineering}. Academic Press, 1982.
\bibitem{Luo2011} Luo, Y.; Zhang, M.;
Environmental modeling and exposure assessment of sediment-associated pyrethroids
in an agricultural watershed. PloS one, 6(1) (2011), e15794.
\bibitem{Macdonald1957} Macdonald, G.;
The Epidemiology and Control of Malaria, Oxford University Press, London, 1957.
\bibitem{Marino2008} Marino, S.; Hogue, I. B.; Ray, C. J.; Kirschner, D. E.;
A methodology for performing global uncertainty and sensitivity analysis in systems biology.
\emph{J. Theor. Biol.} {\bf 254} (2008), 178--196.
\bibitem{Mckay1979} 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. \emph{Technometrics} \textbf{21} (1979), 239--245.
\bibitem{McLeod2006} 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.
\emph{Math. Biosci. Eng}. {\bf3}(3) (2006), 527--544.
\bibitem{Palatnik-de-Sousa2004} 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.
\emph{Annals of Brazilian academy of Science}. \textbf{76}(3) (2004), 583--593.
\bibitem{Pontryagin1962} Pontryagin, L. S.; Boltyanskii, V. G.; Gamkrelidze, R. V.;
Mishchenko, E. F.;
The mathematical theory of optimal processes. Wiley, New York, 1962.
\bibitem{Ribas2013} Ribas, L. M; Shimozako, H. J.; Massad, E.;
Estimating the Optimal Control of Zoonotic Visceral Leishmaniasis by the Use of a
Mathematical Model. \emph{Hindawi: The Scientific World Journal}. \textbf{2013} (2013), 1--6.
\bibitem{Ross1911} Ross, R. ;
\emph{The Prevention of Malaria}, John Murray, Oxford, London, 1911.
\bibitem{Singh2013} Singh, N.; Mishra, J.; Singh, R.; Singh S.;
Animal Reservoirs of Visceral Leishmania in India. \emph{Journal of Parasitology},
{\bf 99} (1) (2013), 64--67.
\bibitem{Sanchez1997} Sanchez, M.; Blower, S.;
Uncertainty and sensitivity analysis of the basic reproduction rate.
\emph{American Journal of Epidemiology}, {\bf 145} (12) (1997), 1127--1137.
\bibitem{standford2006} Standford University online resource on Visceral Leishmaniasis (2006).\\
https://web.stanford.edu/class/humbio103/ParaSites2006/Leishmaniasis/visceral.htm .
Accessed on Sept. 19, 2015.
\bibitem{Stauch2012} Stauch, A.; Duerr, H.-P.; Dujardin, J.-C.; Vanaerschot, M.;
Sundar, S, et al.;
Treatment of Visceral Leishmaniasis: Model-Based Analyses on the Spread of
Antimony-Resistant L. donovani in Bihar, India.
\emph{PLoS Negl Trop Dis} \textbf{6} (12) (2012), e1973.
doi:10.1371/journal.pntd.0001973
\bibitem{Stauch2014} 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.
\emph{PLoS Negl Trop Dis} \textbf{8} (4) (2014), e2810. doi: 10.1371/journal.pntd.0002810
\bibitem{Subramanian2015} Subramanian, A.; Singh, V.; Sarkar, R. R.;
Understanding Visceral Leishmaniasis Disease Transmission and its Control-A Study
Based on Mathematical Modeling. \emph{Mathematics} \textbf{3} (3) (2015), 913-944.
doi:10.3390/math3030913
\bibitem{VanDen2002} van den Driessche, P.; Watmough, J.;
Reproduction numbers and the sub-treshold endemic equilibria for compartmental models
of disease transmission. \emph{Mathematical Biosciences}, \textbf{180} (2002), 29-48.
\bibitem{Farrington2003} Farrington, C. P.; Whitaker, H. J.;
Estimation of effective reproduction numbers for
infectious diseases using serological survey data. \emph{Biostatistics}, \textbf{4} (4)
(2003), 621-632.
\bibitem{WB2015} World Bank; World Development Indicators: Population dynamics (2015) \\
http://wdi.worldbank.org/table/2.1. Accessed January 20, 2016.
\bibitem{WHO} World Health Organization;
Mosquitos and other biting Diptera.\\
http://www.who.int/water\_sanitation\_health/resources/vector007to28.pdf (2015).
Accessed January 20, 2016.
\bibitem{WHO2016} World Health Organization;
Leishmaniasis (2015). \\
http://www.who.int/mediacentre/factsheets/fs375/en/. Accessed January 4, 2016.
\bibitem{Xiefei2007} Yan, X.; Zou, Y.; Li, J.;
Optimal quarantine and isolation strategies in epidemics control.
\emph{World Journal of Modelling and Simulation}. \textbf{3} 3 (2007), 202--211.
\bibitem{Xiefei2008} Yan, X.; Zou, Y.; Optimal and sub-optimal quarantine and
isolation control in SARS epidemics. \emph{Mathematical and Computer Modelling}.
\textbf{47} (2008), 235--245.
\end{thebibliography}
\end{document}