Stability analysis of five-grade Leishmania epidemic model with harmonic mean-type incidence rate

In this paper, we discuss the Anthroponotic Cutaneous Leishmania transmission. In addition, we develop a mathematical model for the Anthroponotic Cutaneous Leishmania transmission and consider its qualitative behavior. We derive the threshold number R0 of the model using the next generation method. In the disease-free case, we carry out the local and global stability under the condition R0 < 1. Moreover, we derive the global stability at the disease-free equilibrium point by utilizing the Castillo-Chavez method. On the other hand, at the endemic equilibrium point, we show the local and global stability to be held under specific conditions and R0 > 1. We also establish the global stability at the endemic equilibrium point with the help of a geometrical approach, which is a generalization of Lyapunov theory, by using a second additive compound matrix. Finally, we take into account the sensitivity analysis of the threshold number with other parameters. We also discuss several graphs of important parameters.


Types of leishmaniasis
I. Mucocutaneous Leishmaniasis (MCL) is a less common form of Leishmaniasis, which is the partial or complete destruction of mucous membranes in the nose, throat, or mouth. It is either CL skin lesion to mucosa membrane or carried by direct bite of a sand fly on mucus membranes.
II. Visceral Leishmaniasis (VL) is the most lethal disease, which appears from days to year after a bite of an infected sand fly. This disease affects the internal organ like the liver, spleen, and bone marrow. Patients are symptomatically present with enlargement of the spleen and liver, weight loss, and episodes of fever. In Pakistan, both MCL and VL are rare.
III. Cutaneous Leishmaniasis (CL) is an ulcer appearance on the skin, which is usually developed in a couple of days after a bite of the sand fly. If the immunity of the patient is strong enough, then mostly these sores may heal spontaneously within a few days. According to the World Health Organization (WHO), 1-1.5 million new cases of CL are reported every year worldwide. Pakistan is one of the countries where the incident of CL is increasing yearly, and its main reason is poverty, negligible planning for vector control, inadequate medical facilities, armed conflicts, and mass migration of humans and cattle.

Treatment methods
Standard drugs: Meglumine antimoniate and sodium stibogluconate are used as a first-line conventional treatment worldwide to treat CL. These drugs are used oral and injectable. Besides their hard administration, it offers many side effects, and drugs are costly as well.
Conventional treatment: Some conventional domestic methods help in healing in CL wound: applying garlic water on the lesion, cleaning wound repeatedly with antiseptics, and placing a hot stone on wound have been reported by people.
Photodynamic therapy: Photodynamic therapy (PDT), the use of light of specific wavelength and energy, is a rapidly evolving therapeutic option for the treatment of CL. The FDA has approved indications of PDT for various skin diseases, actinic keratosis, basal cell carcinoma, superficial squamous cell carcinoma, psoriasis, acne vulgaris, and sarcoidosis, verruca Vulgaris, and condyloma acuminatum. Treatment protocols of PDT are easy, patients' comfort and lesion response to PDT make this the best option to treat CL.
Treatment protocols: Herein we present a number of cases of CL lesions at different sites of the body treated with PDT from Pakistan. These lesions had been diagnosed clinically as cutaneous Leishmaniasis (CL). The patients had received treatment (i.e., pentavalent antimoniate) with no satisfactory outcome. Clinical assessment showed that the patients had developed resistance to pentavalent antimoniate. The patients were planned for photodynamic therapy (PDT). Written and informed consent was obtained from the patient. To perform PDT, the lesions and the adjacent skin were cleaned and scrubbed to remove the necrotic layer and exudates dried up. The photosensitizer, methyl aminolevulinate (MAL), was applied locally on the lesions under an adhesive covering, followed by an incubation period of 3 hours for the absorption of the photosensitizer. After the incubation period, the cream was removed, and the lesion was washed with normal saline. The lesion was irradiated with red laser light (wavelength 635 nm); a light energy dose of 75 J/cm 2 was delivered as per the institutional protocol. The patient received three sessions of PDT; the first two sessions were repeated at two weeks intervals, whereas the third session was given at one-month intervals. The patient showed a complete response to PDT, being followed up regularly, and has remained disease-free for over two years.

Model formulation
Mathematical modeling of epidemic diseases has been widely studied by researchers [19][20][21]. The model is a modified version of Khan et al. [22], in which the total population N(t) is partitioned into the following subpopulations. These classes are classified into four human population subclasses S h (t), E h (t), I 1h (t), I 2h (t), and R h (t), denoting the susceptible, exposed, infected people from grades I-III, infected people in grades IV and V, and recovered people. The three vector population subclasses S v (t), E v (t), and I v (t) represent susceptible, exposed, and infected vectors; the total human population is rep- The model is given as follows: with

Basic properties of the model
The total-population dynamics is represented by The feasible region (biological) is From Eqs. (2) and (3) we obtain which shows that the model is well posed and is a positively invariant domain.
System (1) is expressed in the form where Here we see that L is the Metzler matrix as it has nonnegative entries on its of-diagonal and B ≥ 0. Hence we conclude that system (1) is positively invariant in R 8 + .
Lemma 2 If solutions of system (1) exist, then they are positive under the initial conditions (2) for all t > 0.
Proof Let us assume that the solutions exist in I for all t ∈ I ⊂ [0, ∞). Consider the second equation of (1). Its solution has the following form: We also take the third equation, and its solution has the following form: We can clearly see from the above solutions that these are strictly positive. In the same fashion, we can show that S h , R h , S v , E v , I 2h , and I v possess nonnegative solutions.

Basic reproductive number R 0
The disease-free equilibrium point of system (1) is Let (E h , I 1h , I 2h , E v , I v ) be our infected compartment. Then from system (1) we have: Using the next generation matrix approach [23,24], the Jacobian matrix J for this system at the disease-free equilibrium point is given by Now decomposing the matrix J in terms of F and V as J = F -V , we get The characteristics equation of [-FV -1 ] becomes The dominant eigenvalue gives us R 0 , that is, .

Local stability
In this section, we establish the local stability of system (1) in this section at the diseasefree point E 0 and endemic equilibrium point E * .

Geometric interpretation of stability
For the DFE, the local stability implies that if there is a small perturbation of the system (e.g., a small number of infected individuals is introduced into the population), then after some time, the system will return to the DFE. However, a big perturbation may be able to drive the system to a different behavior, for example, convergence to the EE. If the stability of the DFE is global, then no matter the size of the perturbation, the disease will not be able to persist in the population.
On the other hand, the local stability of an equilibrium point means that if you put the system somewhere nearby the point, then it will move itself to the equilibrium point in some time. The global stability means that the system will come to the equilibrium point from any possible starting point (i.e., there is no "nearby" condition). Even in more physical interpretation, it could be like this: If DFE or EE is locally stable, then all epidemiological situations not-so-much different from the given stable equilibrium will (with time) evolve to (or transform into) the equilibrium point. Also, this means that the equilibria are stable to small perturbations: if you push the situation a bit out from the equilibrium point, then the situation will return back on its own (from the physicist's point of view, this means that the equilibrium may be a stable situation in real life because the real world always is somewhat noisy). The global stability of an equilibrium point in this case maybe described as "the inevitable fate of the epidemic process regardless of its starting situation". However, a caveat should be put that this "inevitability" holds as long as the world strictly follows the underlying mathematical model of the epidemic process. For still even more lay (or medical) audience, all these things may be easily interpreted in terms of "adding or removing a small number of infectious people", but we believe it would make the central idea less clear.

At disease free equilibrium point
Proof For system (1), the Jacobian matrix at E 0 is given by We can clearly see that three eigenvalues are negative, namely, λ 1 = -μ h , λ 2 = -μ h , and λ 3 = -μ v . For the rest of eigenvalues, we consider the following reduced matrix: After simplification, we get: The eigenvalues of J |0| 2 take the following form: The last eigenvalue also has the negative real part since implies that Thus under the condition R 0 < 1, the real part of all eigenvalues are negative. Hence, for R 0 < 1, system (1) is locally asymptotically stable.

At endemic equilibrium point
, , then the endemic equilibrium E * of system (1) is locally asymptotically stable. (1), the Jacobian matrix at E * is given by

Proof. For system
where , a 66 = acI * 1h + adI * 2h + μ v , and a 76 = acI * 1h + adI * 2h . Clearly, one eigenvalue of the Jacobian matrix J | * | of model (1) around the disease present equilibrium point E * is negative, that is, λ 1 = -μ h < 0. For the remaining seven eigenvalues, we take the following reduced matrix: After simplification we get: where The eigenvalues of J | * | 2 take the following forms:

Global asymptotic stability
Here we discuss the global stability analysis of model (1) for both disease-free and endemic equilibria. We use the method of Castillo-Chavez et al. [25] to establish the global stability for disease-free equilibrium, whereas for the global stability of endemic equilibrium, we use the generalization of Lyapunov theory [26].

At disease-free equilibrium point
For model (1), the global stability at the disease-free point is achieved by taking into account the Castillo-Chavez approach [25]. The method is summarized by the reduction of the proposed model (1) to the following two subsystems: In system (20), χ 1 and χ 2 represent the numbers of uninfected and infected individuals, respectively, that is, The disease-free equilibrium is denoted by E 0 and defines as E 0 = (χ 0 1 , 0). Thus the existence of the global stability at the disease-free equilibrium point depends on the following two conditions: 1. If dχ 1 dt = G(χ 1 , 0), then χ 0 1 is globally asymptotically stable.
By using model system (1) we have Thus from equation (23) we have that χ 1 → χ 0 1 as t → ∞. So χ 1 = χ 0 1 is globally asymptotically stable. Now As the total population is bounded by , which implies thatH(χ 1 , χ 2 ) ≥ 0. Clearly, B is an M-matrix, and hence both the conditions are proved, so by Lemma 1 the disease-free equilibrium point E 0 is stable globally asymptotically.

Endemic equilibrium (global stability)
"For the global stability of (1) at the endemic equilibrium E * , we use the geometrical approach [26]. Thus we investigate the sufficient condition through which the E * is globally asymptotically stable. Therefore consider the differential equatioṅ where the open set U ⊂ R n is simply connected, and f : U → R n is a function such that f ∈ C 1 (U). Assuming that f (x * ) = 0 is the solution of equation (25), for x(t, x 0 ), the following are true: 3. There exist a compact absorbing set K ∈ U. 4. System (25) has a unique equilibrium. The solution x * is said to be globally asymptotically stable in U if it is locally asymptotically stable and all trajectories in U converge to the equilibrium x * . For n ≥ 2, a condition satisfied for f , which precludes the existence of nonconstant periodic solution of equation (25), is known as the Bendixson criterion. The classical Bendixson criterion div f (x) < 0 for n = 2 is robust under C 1 [26]. Furthermore, a point x 0 ∈ U is wandering for equation (25) if there exist a neighborhood N of x 0 and τ > 0 such that N ∩ x(t, N) is empty for all t > τ . Thus the following global stability principle established for autonomous system in any finite dimension. (25) (i.e., robust under C 1 local perturbation of f at all nonequilibrium nonwandering point for equation (25)), then x * is globally asymptotically stable in U, provided that it is stable.

Lemma 4 If the conditions (3)-(4) and Bendixson criterion are satisfied for equation
Define a matrix-valued function P on U by Further, assume that P -1 exists and is continuous for x ∈ K . Now define the quantitȳ where J [2] is the second additive compound matrix of J, that is, J(x) = Uf (x) and B = P f P -1 + PJ [2] P -1 . Let (B) be the Lozinski measure of the matrix B with respect to the norm · in R n [27] defined by Hence ifq < 0, then the presence of any orbit gives rise to a simple closed rectifiable curve such as a periodic orbit and heterocyclic cycle. (25) is globally asymptotically stable in U ifq < 0.

Lemma 5 Let U be simply connected, and let conditions (3)-(4) be satisfied. Then the unique equilibrium x * of equation
Now we apply the above techniques to prove the global stability of model (1) at the endemic equilibrium. We have the following stability. 2 , and R 0 > 1, then model (29) is globally asymptotically stable at the endemic equilibrium E * and unstable otherwise.

Theorem 4.2 If
Proof To prove the global asymptotic stability of the proposed model (1) at the endemic equilibrium E * , let us consider the subsystem of (1) Let us start with the Jacobian matrix of system (29) The third additive compound matrix is given by Let J be the Jacobian matrix of system (29) given by The third additive compound matrix of J is denoted by J |3| and defined as where Let us choose a function p( A direct computation shows that B = P f P -1 + PJ |3| P -1 , which becomes so that B = P f P -1 + PJ |3| P -1 , where Consequently, The integration of the Lozinski measure μ(B) and taking the limits as t → ∞ lead to the following equations: Thus combining these four inequalities, we get the following equation: The system containing the first four equations of model (1) is globally asymptotically stable around its interior equilibrium (S (34) Now rewrite the system in the form The integrating factors for the system are e t(β+k 2 +μ h ) , e t(μ h +θ 1 ) , e t(μ h ) , and e t(μ v ) .
Using the integrating factors, we solve the system. So for large time t, that is, t → ∞, I 1h → I * 1h , I 2h → I * 2h , R h → R * h , and I v → I * v , which means that the endemic equilibrium point E * is globally asymptotically stable.

Sensitivity analysis
Determining the parameters that are helpful in decreasing the spread of infectious disease is carried out by sensitivity analysis. Forward sensitivity analysis is considered a vital component of disease modeling although its computation becomes tedious for complex biological models. Sensitivity analysis of R 0 has received much attention from the ecologists and epidemiologists. The behavior of the R 0 with some suitable values of the parameters are shown in Fig. 1.

Definition 1
The normalized forward sensitivity index of R 0 that depends differentiably on a parameter ω is defined as Three methods are normally used to calculate the sensitivity indices: (i) by direct differentiation, (ii) by a Latin hypercube sampling method, and (iii) by linearizing system (1) and then solving the obtained set of linear algebraic equations. We will apply the direct differentiation method as it gives analytical expressions for the indices. The indices not only show us the influence of various aspects associated with the spreading of infectious disease but also gives us important information regarding the comparative change between R 0 and different parameters. Consequently, it helps in developing the control strategies.

Optimal control
By using these control variables our control problem becomes The goal of our optimal control strategies is minimizing the infectious and exposed human population, the vector population, sandfly biting rate, and the cost of implementing the control by using possible minimal control variables u 1 (t), u 2 (t), and u 3 (t). To do this, we use the bounded Lebesgue-measurable control to construct the objective functional subject to system (37). In the objective functional, g 1 , g 2 , g 3 , and g 4 represent the weight constants of the exposed, infectious human, and of vector population, respectively, d 1 , d 2 , and d 3 are weight constants of human self-protection, human treatment, and vector control, respectively. The terms (1/2)d 1 u 2 1 (t), (1/2)d 2 u 2 2 (t), and (1/2)d 3 u 2 3 (t) describe the costs of disease interventions. The cost associated with the first control strategy u 1 (t) comes from the cost of sandfly repellent lotions, electric mats, and mosquito bed nets. The cost associated with the second control strategy u 2 (t) is the cost of expensive medication of human class. The cost associated with the third control strategy u 3 (t) can arise from applying different types chemical pesticides to kill sand flies at any stage of their life. We have assumed the costs as proportional to the square of the corresponding control function. We aim to find control functions such that subject to system (4). The control set D is defined as is Lebesgue measurable on [0, 1], 0 ≤ u i (t) < 1, i = 1, 2, 3 .

Existence of the control problem
For the existence of the control problem, we consider the control system (37) with initial conditions at t = 0. In case of bounded Lebesgue-measurable controls and nonnegative initial conditions, there exists a nonnegative bounded solution of the state system [28,29]. For optimal solution of the system, we first find the Lagrangian and Hamiltonian. We define the Lagrangian of the control problem (37) as where We define the Hamiltonian H to find the minimal value of the Lagrangian as follows: Theorem 6.1 There exists an optimal control u * = (u * 1 , u * 2 , u * 3 ) ∈ D such that J(u 1 , u 2 , u 3 ) = min (u 1 ,u 2 ,u 3 )∈D J(u 1 , u 2 , u 3 ) subject to system (37) and initial conditions at t = 0.
Proof The control variables and the state variables are nonnegative. So we use the result in [28,29] for the existence of an optimal control. The necessary convexity of the objective functional in u 1 , u 2 , and u 3 is satisfied here. Also, by definition the set of the control variables (u 1 , u 2 , u 3 ) ∈ D is convex and closed. The compactness, which we need for the existence of an optimal control, is confirmed by the boundedness of the optimal system. Also, the integrand in the functional (38), g 1 E h (t) + g 2 I 1h (t) + g 3 I 2h (t) + g 4 N v (t) + 1 2 (d 1 u 2 1 (t) + d 2 u 2 2 (t) + d 3 u 2 3 (t)), is convex on the control set D. We can find a constant η > 1 and positive numbers ξ 1 and ξ 2 such that J(u 1 , u 2 , u 3 ) ≥ ξ 1 (|u 1 | 2 , |u 2 | 2 , |u 3 | 2 ) η 2ξ 2 , as the state variables are bounded. This completes the proof of the existence of an optimal control.
We apply the necessary conditions to the Hamiltonian H in system (40).
(41) Figure 2 Plots of different classes showing the difference between "with control" and "without control" recommended that proper record keeping and documentation systems for leishmaniasis be initiated by health authorities at the local, provincial, and national levels and be well maintained to identify leishmaniasis outbreaks so that control measures can be started well in time. Further, IDP camps must be monitored regularly to minimize the risk that nonendemic areas will be exposed to the disease by infected IDPs.