Global Analysis and Optimal Control of a Periodic Visceral Leishmaniasis Model

In this paper, we propose and analyze a mathematical model for the dynamics of visceral leishmaniasis with seasonality. Our results show that the disease-free equilibrium is globally asymptotically stable under certain conditions when R0, the basic reproduction number, is less than unity. When R0 > 1 and under some conditions, then our system has a unique positive ω-periodic solution that is globally asymptotically stable. Applying two controls, vaccination and treatment, to our model forces the system to be non-periodic, and all fractions of infected populations settle on a very low level.


Introduction
Leishmaniasis is a vector borne, zoonotic disease transmitted to humans and animals through the bite of infected phlebotomine sandflies.The World Health Organization considers leishmaniasis as the sixth most important endemic disease in the world because it has the second-highest number of affected people after malaria and causes an important public health problem in several countries around the world [1].
Leishmaniaisis has many forms, the most important being cutaneous leishmaniasis and visceral leishmaniasis.The cutaneous form is rarely fatal, but it can lead to tissue destruction, scarring, and serious visible impairment.These can lead to social stigmas because of the deformation of the face with permanent scarring.The visceral form, which is almost 100% fatal if left untreated, requires a lengthy and very costly treatment that should be taken under medical supervision [2].In some cases, depending on the immunity, after the treatment of visceral leishmaniasis, another form of the disease may develop, the post-kala-azar dermal leishmaniasis (PKDL) form, which is not a fatal form; most cases will self-cure, but severe and chronic cases must be treated, and it plays an important role for the dynamics of the disease as it acts as a reservoir for the parasite [3].Naturally, a time delay could be inserted for PKDL development, as this would make the mathematics very rich and interesting.
Leishmaniasis is distributed in around 90 countries of the world, most of which are developing countries, and more than 90% of the cases are found in India, Nepal, Iran, Iraq and Sudan; there are approximately 350 million people at risk of leishmaniasis in these areas [4].Leishmaniasis is very hard to control, because it can be transmitted to many mammals, which play the role of a reservoir to the disease, and the transmission to humans can be restarted at any time, causing a high prevalence of the disease even after a long period of no disease or low prevalence of the disease.
Some studies have suggested that there is a correlation between seasonal variations and the abundance of sandflies, which leads to a link between seasonal variations and the prevalence of the disease [5][6][7][8][9][10].Tiwary et al. [9] found that sandflies are highly prevalent in the rainy season and that the lowest rate of infection with visceral leishmaniasis was observed in the rainy season.This might be due to the flushing effect of rainfall on immature, emerging sandflies.Turki et al. [10] found that the increase in temperature during summer can lead to an increase of cutaneous leishmaniasis cases.They also found an indication that rainfall has no impact or influence on the species of sandflies in the highlands; however, they found that in the lowlands, increased rainfall can lead to increased numbers of cases of cutaneous leishmaniasis, except during the winter season.There is not yet a vaccine for visceral leishmaniasis; however, this is biologically feasible, as many people before the age of vaccinology use a procedure called leishmanization to protect themselves from the infection.Leishmanization procedure is the process of introducing a live but weak form of the parasite to stimulate immune response in order protect from the real infection [11].
Some mathematical models have been developed to study the dynamics of visceral and cutaneous leishmanaisis.Dye [12] designed a mathematical model for the dynamics of zoonotic visceral leishmaniasis, and his model also focused on the dynamics of infection in dogs.His results show that insecticides are the most effective control method; the second-best strategy is to reduce susceptibility to leishmaniasis by vaccinating people or dogs, or by eliminating childhood malnutrition where it is common.Both killing vectors and reducing susceptibility (by whatever means) are more effective than killing dogs or treating them with drugs.Dye et al. [13] formulated an age-structured mathematical model to investigate the effects of visceral leishmaniasis on children, as well as the correlation between malnutrition and the susceptibility of leishmaniasis.His model consists of a system of partial differential equations for the human population, but it has no equation for the animal reservoir nor for the vector reservoir.Lysenko et al. [14] developed an epidemic model for the transmission of cutaneous leishmaniasis.However, the mechanisms that govern the dynamics of the transmission of leishmania parasites are not explicitly stated; therefore, no control strategies for the disease can be proposed.Recently, Rabinovich et al. [15] added some realism to this model by making the probability of developing cutaneous lesions a function of the number of infective bites by sandflies.Burattini et al. [16] formulated a mathematical model to describe the dynamics of the transmission of leishmaniasis, including populations of vector, human and animal (dogs) hosts.Assuming that all populations are both sources and sinks of leishmania parasites, they estimated R0: the infection persistence threshold condition from the sum of individual terms for transmission in humans and in dogs.However, their model did not include the role of PKDL humans in the disease transmission, who have a major role, as PKDL humans act as a source of infection for vectors [17].Chaves et al. [1] developed a mathematical model for the American cutaneous leishmaniasis.Their model describes the transmission of American cutaneous leishmaniasis in three different populations, a human host population, an animal reservoir population and a vector population, but they considered humans as sinks of the infection only; hence humans in their model are not a source of infection for the vectors.
Elmojtaba et al. [18] developed a mathematical model to understand the dynamics of visceral leishamniasis and determined the local and global stability of the disease-free equilibrium, as well as the local stability of the endemic equilibrium.Elmojtaba et al. [19] developed a model to study the role of cross-immunity between two different strains of leishmaniasis in the eradication of the disease; Elmojtaba et al. [20] also studied the dynamics of malaria and leishmaniasis co-infection to set up some control measures against the spread of these two diseases.Elmojtaba et al. [21] studied the role of infected immigrants and vaccination on the dynamics of visceral leishmaniasis under these settings.All these models used the same assumptions with constant coefficients; however, here in this paper, we consider a model with periodic coefficients.
This paper is organized as follows: the first section gives an introduction for the topic; in the second section, the model is built; the third section gives an intensive mathematical analysis of the model; the optimal control problem is proposed and analyzed in the fourth section; numerical simulations have been carried out and are presented in the fifth section; and a discussion is given in the last section.

Model Formulation
To formulate this model, we follow Elmojtaba et al. [18], and we consider the transmission of the disease between three different populations, a human host population, N H (t); a reservoir host population, N R (t); and a vector population, N V (t).We let the human host population be divided into four categories, susceptible individuals, S H (t); those who are infected with visceral leishmaniasis, 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 Similarly, we let the reservoir host population be divided into two categories, a susceptible reservoir S R (t), and an infected reservoir I R (t), such that We let the vector population have two categories, susceptiblesandflies S V (t), and infected sandflies It is assumed that susceptible individuals are recruited into the population at a constant rate Λ H and acquire infection with leishmaniasis following contact with infected sandflies at a per capita rate a(t)b I V N H , where a(t) is the per capita biting rate of sandflies on humans (or reservoirs) and b is the transmission probability per bite per human.The per capita biting rate of sandflies a is equal to the number of bites received per human from sandflies due to conservation of bites mechanism.Infected humans die as a result of leishmaniasis at an average rate δ or receive treatment at an average rate α 1 ; a fraction σ of those that are treated recover and acquire permanent immunity, and a fraction (1 − σ) develop PKDL.Humans with PKDL are treated at an average rate α 2 or recover naturally at an average rate β; they acquire permanent immunity in both cases.There is a per capita natural mortality rate µ h in all human subpopulations.
Susceptible reservoirs are recruited into the population at a constant rate Λ R and acquire infection with leishmaniasis following contact with infected sandflies at a rate a(t)b I V N H , where a(t) and b are as described above.We assume that the transmission probability per bite is the same for human and the reservoir because sandflies do not distinguish between humans and reservoirs.It is also assumed that reservoirs do not die as a result of the disease, but a per capita natural mortality rate µ r occurs in the reservoir population.
Susceptible sandflies are recruited at a constant rate Λ V and acquire leishmaniasis infection following contact with humans infected with leishmaniasis or having PKDL, or reservoirs infected with leishmaniasis, at an average rate equal to a(t , where a(t) is the per capita biting rate and c is the transmission probability for sandfly infection.Sandflies suffer natural mortality at a per capita rate µ v regardless of their infection status.
We assume that a(t) (the per capita biting rate of sandflies on humans) is a time-dependent ω-periodic function.With this assumption and the description of the terms, we obtain the following system of differential equations: All parameters of the model are assumed to be non-negative; furthermore, because Equation (1) monitors living populations, it is assumed that all the state variables are non-negative at time t = 0.It is noted that in the absence of the disease (δ = 0), the total human population size µ v } is a positively invariant domain, and thus the model is epidemiologically and mathematically well posed; it is sufficient to consider the dynamics of the flow generated by Equation (1) in this positively invariant domain Ω.

Analysis of the Model
In this section, we analyze Equation ( 1) to obtain the equilibrium points of the system and their stability.We consider the equations for the proportions by first scaling the subpopulations N H , N R and N V using the following set of new variables: and let m = N V N H be the female vector-human ratio defined as the number of female sandflies per human host [18].We note that the ratio m is taken as a constant because it is well known that a vector takes a fixed number of blood meals per unit time independent of the population density of the host.Similarly, we let n = N V N R be the female vector-reservoir ratio defined as the number of female sandflies per reservoir host.
Differentiating with respect to time t and simplifying gives the following reduced system of differential equations:

Mathematical Properties
We let C denote all continuous functions on the real line.Given f ∈ C, and if f is ω-periodic, then the average value of f on a time interval [0,ω] can be defined as The maximum and minimum values of f on a time interval [0,ω] are denoted as f M and f m , respectively, and are defined as Then the following result is valid: This proposition shows that the scaled system (Equation ( 2)) is epidemiologically and mathematically well-posed, and it is sufficient to consider the dynamics of the flow generated by Equation ( 2) in this positively invariant domain Ω.Here R 8 + denotes the non-negative cone of R 8 including its lower-dimensional faces.We denote the boundary and the interior of Ω by ∂Ω and • Ω, respectively.
Proof.The proof is straightforward.
We let (R n , R n + ) be an ordered n-dimensional Euclidean space with a norm ||.||, and we assume that int(R n + ) = ∅.We let x, y ∈ R n , and we say The following lemma will be useful for proving global stability of E 0 .
Lemma 1. (Lemma 2.1 in [22]).Let A(t) be a continuous, cooperative, irreducible, and ω-periodic n × n matrix function.Let Φ A(.) (t) be the fundamental matrix of the linear non-autonomous differential equation ẋ = A(t)x, where x is a n × 1 vector.Let µ = 1 ω ln r(Φ A(.) (ω)), where r(Φ A(.) (ω)) is the spectral radius of the monodromy matrix Φ A(.) (ω).Then there exists a positive ω-periodic function v(t) such that e µt v(t) is a solution of ẋ = A(t)x.Now let x = (i h , p h , i r , i v ) T , F (t, x) denotes the new infection rate and V (t, x) denotes the rate of transfer of individuals between infection compartment; therefore the vectors F (t, x) and V (t, x) given by Therefore, see [23], F(t) and V(t) given by and Now, F(t) is non-negative and −V(t) is co-operative (off-diagonal elements are non-negative).Thus, F(t) − V(t) is continuous and co-operative and irreducible.Therefore, Lemma 1 is applicable to the system Ẋ = (F(t) − V(t))X, where X = (i h , p h , i r , i v ) T .

Basic Reproduction Number
We let, for each s ∈ R, the 4 × 4 matrix Y(t, s), ∀t s satisfy the ω-periodic system: and Y(s, s) = I, where I is the 4 × 4 identity matrix.
We let C ω be the ordered Banach space of all ω-periodic functions from R to R 4 , which is equipped with maximum norm ||.|| and the positive cone C + ω = {φ ∈ C ω : φ(t) 0, ∀t ∈ R}.We define the linear operator L : C ω → C ω as follows: The definitions of the basic reproduction number (R 0 ) with seasonality in literature was given by [23].According to the definition, R 0 is defined as follows: In this paper, we follow the definition of [23] and analyze the model of Equation (2) in terms of the threshold parameter R 0 .
We use following lemma to infer stability of E 0 in terms of R 0 .
Proposition 3. If R 0 < 1 then the disease-free equilibrium E 0 of the system given by Equation ( 2) is globally asymptotically stable.

Existence and Permanence of the Endemic Periodic Solution
In this section, we study the existence and permanence of the endemic periodic state of the model of Equation (2).
Before giving the main result, we give following definitions from [26].
Definition 1.The system given by Equation ( 2) is uniformly persistent if ∃ an η > 0 (depending only on parameter values and not on initial conditions) such that for any initial value (s h (0), i h (0), p h (0), r h (0), s r (0), i r (0), ) such that every solution (s h (t), i h (t), p h (t), r h (t), s r (t), i r (t), s v (t), i − v(t)) of the system given by Equation (2) satisfies Definition 2. The system given by Equation ( 2) is said to be permanent if there exists a compact region Ω 0 ⊂ int(Ω) such that every solution of the system given by Equation ( 2) with initial condition (s h (0), i h (0), p h (0), r h (0), s r (0), i r (0), will eventually enter and remain in the region Ω 0 .
A natural concept of dissipation is to assume that there is a bounded set into which every orbit eventually enters and remains [27].Clearly, for a dissipative dynamical system, proving permanence is equivalent to proving uniform persistence.
We consider the following sets: We let f : X → X be a continuous map, and we define the following set: where f n is the nth iteration of the function.
The following lemma is used to show the uniform persistence of the system given by Equation (2).

Lemma 3. ([26]
) Let W s (X) be the stable set of X. Assume that 1. f (X 0 ) ⊆ X 0 and f has a global attractor A.
We claim the following result: Proposition 4. If R 0 > 1 then the solutions of Equation ( 2) are uniformly persistent, and the system admits at least one positive ω-periodic solution.

+
Let P 1 be the associated Poincaré map defined as P 1 := T(ω).We first show that P 1 is uniformly persistent with respect to (X 0 , ∂X 0 ).
It is clear that the sets X and X 0 are positively invariant for Equation (2).Now by Proposition 1, the bounded set Ω attracts every solution of Equation ( 2), and Ω is compact.Thus the Poincaré map P 1 is point-dissipative and compact on X.Therefore, it follows from Theorem 1.1.3in [26] that there is a global attractor A of P 1 that attracts each bounded set in X.

Global Stability of the Positive Periodic Solution
We claim the following result: Proposition 5.If R 0 > 1 and Λ H N H − 4δ > 0, then the system given by Equation ( 2) has a unique positive ω-periodic solution that is globally asymptotically stable.
Proof.If R 0 > 1, then we have by Proposition 4 that the system given by Equation ( 2) admits a positive periodic solution.
We construct the following Lyapunov function: We use the following formula: |x| = sgn(x)x to calculate the upper right-hand derivative (Dini's derivative) of L(t).

Constructing the Problem
We use two different interventions in the model of Equation ( 1) to reduce the visceral leishmaniasis, namely, vaccination and treatment using antibiotics.The vaccinated population is increased by a proportion u 1 (t)σ 1 of the susceptible individuals who are successfully vaccinated, where u 1 (t) is the per-week vaccination rate and σ 1 is the vaccine efficiency.The vaccinated population is decreased as a result of the waning of vaccine-based immunity (at a rate ) to become susceptible again and die (natural deaths) at a rate µ h .
The system of non-linear differential equations representing the effect of different interventions on our basic model from Equation ( 1) is given as follows: where the initial conditions S H (0), I H (0), P H (0), R H (0), V(0), S R (0), I R (0), S V (0), and I V (0) are given, and the definitions of the above model parameters are listed in Table 1.The control functions u 1 (t) and u 2 (t) are bounded and are Lebesgue integrable functions, where the control u 1 (t) measures the rate at which susceptible individuals are vaccinated in each time period, and the control u 2 (t) represents the effort of treatment of infected individuals.This performance specification involves the number of infection individuals as well as the cost of applying a vaccination control (u 1 (t)) and treatment control (u 2 (t)) in individuals with visceral leishmaniasis.The objective functional to be minimized is defined as where t f is the final time and A, C and D are positive constants to balance the cost factors.The costs can include funds needed for the cost of productive time loss for per premature death, the cost of vaccination per fully immunized person, the cost of vaccination per fully immunized person in high emergencies, the cost for medicine, and hospital admission per case.More often, the cost of implementing a control would be nonlinear.We thus seek to find an optimal control u * 1 (t) and u * 2 (t) such that where is the control set.

Analysis of the Optimal Control Problem
The necessary conditions that an optimal control must satisfy come from Pontryagin's Maximum Principle [35].This principle converts the system given by Equation ( 9) into the problem of minimizing pointwise a Hamiltonian, H, with respect to u 1 and u 2 : where g i is the right-hand side of the differential equation of the ith state variable and λ i are the adjoint variables by applying Pontryagin's Maximum Principle [35][36][37].
Proposition 6.Given an optimal control u * 1 and u * 2 and solutions S * H , I * H , P * H , R * H , V * H , S * R , I * R , S * V , and I * V of the corresponding state system given by Equation ( 9) that minimizes J(u 1 ,u 2 ) over U, then there exists adjoint variables λ i , 1 ≤ i ≤ 9, satisfying the following: with transversality conditions Furthermore, the control pair (u * 1 , u * 2 ) is given as Proof.The existence of an optimal control is guaranteed using the result by [38].Thus, the differential equations governing the adjoint variables are obtained by the differentiation of the Hamiltonian function evaluated at the optimal controls, giving the stated adjoint system from Equations ( 14) and (15).Furthermore, differentiating the Hamiltonian function with respect to the control variables in the interior of the control set U, where a i < u i < b i , i = 1, 2, we have Solving Equation ( 17) for controls (u * 1 , u * 2 ) and using the bounds on the controls, the characterization of Equation ( 16) can be derived.

Numerical Simulations
In this section, we present some numerical results obtained using data from [18], assuming that the force of infection is a periodic function with a period of 12 months, with one peak in each period.The parameter values are given in Table 1.To solve the optimal control problem, we used a forward-backward Runge-Kutta of order 4 (RK4) method.To solve the state system of ordinary differential equations, a simple forward RK4 method was applied, but to solve the adjoint ODE, the RK4 method had to be adapted to account for solving backwards in time.
In Figure 1, it is clear that without control, the fraction of infected humans decreases to low levels, then starts to increase and repeats this process periodically; however with optimal control, it decreases to a very low level and settles at that level.Figure 2 shows that the fraction of infected reservoirs without control increases at first, then decreases to a very low level, before increasing again; this process is repeated, and the fraction of infected reservoirs with control increases at first, then decreases again and settles on a medium level.This is because there is no control applied to the reservoir population.
From Figure 3, it is clear that the fraction of infected vectors without control increases at first, then decreases to a low level, before increasing again; this process is repeated, and the fraction of infected vectors with control increases slightly at first, then decreases again and settles on a very low level.The reason is that the control applied to the human population affects the vector population indirectly.
In Figure 4, we show that without control, the fraction of PKDL-infected humans decreases at first, then increases; this process repeats again until the end of the simulation time, while with control, this fraction decreases to a very low level.
The main conclusion drawn from the simulation is that when the optimal control is applied, the periodicity disappears from the model.

Discussion
Visceral leishmaniasis is associated with many different factors, such as complex vector dynamics, the role of an animal reservoir, and social and economical factors; one of the most important factors is the seasonal dynamics of the disease, which makes its control very difficult, as seasonality makes the recurrence of epidemics very common.
In this paper, we present and analyze a mathematical model for the dynamics of visceral leishmaniasis, assuming that the biting rate of the vector is a time-dependent function, to capture the seasonal behavior of the disease.Our results show that if R 0 is less than unity, then the disease-free equilibrium is globally asymptotically stable, and if R 0 > 1, then our system is uniformly persistent and admits at least one positive ω-periodic solution.Two time-dependent control variables, namely, the treatment of infected humans and the vaccination of susceptible individuals, were added to the model in order to investigate the best control strategy against the disease.The analysis and numerical simulation of the optimal control problem show that using an optimal combination of vaccination and treatment will make the system periodic-free by forcing it to reach the disease-free state.The optimal combination of vaccination and treatment means that an extensive use of treatment for all infected individuals and massive vaccination for susceptible humans should be applied.As shown previously, there is not yet a vaccine for visceral leishamaniasis; however, our results suggest that the research should continue to develop such as vaccine, as if it is applied perfectly with the treatment, it will lead to a disease-free state of the system, and more importantly, it will stop the recurrence of the disease.Our model assumes that vaccinated individuals could not be infected with visceral leishmaniasis and that the waning rate of the vaccine is very low; therefore the developed vaccine should be very effective to achieve these goals.Our results also suggest that the vaccine should be very cheap or otherwise freely distributed, as it should be applied widely, and almost all of the susceptible individuals should be vaccinated to eradicate the disease.

Figure 1 .
Figure 1.Simulation results showing the infected humans.

Figure 2 .
Figure 2. Simulation results showing the infected reservoir.

Figure 3 .
Figure 3. Simulation results showing the infected vector.

Table 1 .
Parameter values of the model.