Control strategies and sensitivity analysis of anthroponotic visceral leishmaniasis model

ABSTRACT This study proposes a mathematical model of Anthroponotic visceral leishmaniasis epidemic with saturated infection rate and recommends different control strategies to manage the spread of this disease in the community. To do this, first, a model formulation is presented to support these strategies, with quantifications of transmission and intervention parameters. To understand the nature of the initial transmission of the disease, the reproduction number is obtained by using the next-generation method. On the basis of sensitivity analysis of the reproduction number , four different control strategies are proposed for managing disease transmission. For quantification of the prevalence period of the disease, a numerical simulation for each strategy is performed and a detailed summary is presented. Disease-free state is obtained with the help of control strategies. The threshold condition for globally asymptotic stability of the disease-free state is found, and it is ascertained that the state is globally stable. On the basis of sensitivity analysis of the reproduction number, it is shown that the disease can be eradicated by using the proposed strategies.


Introduction
Leishmaniasis is a family of infectious diseases. This disease is transmitted by different species of phlebotomize sandflies. Anthroponoticvisceral leishmaniasis (AVL) is caused by Leishmania donovani and is transmitted by sandfly. Human is the main reservoir of the virus. The clinical symptoms of AVL include fatigue, prolonged fever, losing weight, bleeding tendency and enlargement of both spleen and liver. The average incubation period of visceral leishmaniasis (VL) is 2-6 months; however, longer and shorter periods (from 10 days to 1 year) [1,26]. The latency period of sandfly is assumed roughly to be 3-7 days [21,22]. Post-kala azar dermal leishmaniasis (PKDL) is the complication of VL in a patient who has recovered from VL. The interval between PKDL and VL is observed from 0 to 6 months in Sudan and from 2 to 3 years in India [27]. Somalia, Yemen, Saudi Arabia, Ethiopia, Kenya and Uganda are the countries highly suffering with AVL [14,19]. For more details and related literature of Leishmaniasis, the readers can refer to [2,3,8,16,17,25,26].
In this paper, we consider the work of Stauch et al. [23] by incorporating some new biological features related to develop a mathematical model. To do this, we introduce the system of nonlinear differential equations to represent the dynamics of the disease. We show that the model is epidemiologically and mathematically well posed. To understand the nature of the initial transmission of the disease, the reproduction number R 0 is obtained by using the next-generation method. On the basis of sensitivity analysis of R 0 , we propose control strategies for disease transmission. To understand the dynamical behaviour and stability analysis, we use the Routh-Hurwitz criteria and stability analysis theory of nonlinear systems of differential equations. The threshold condition for globally asymptotic stability of the disease-free state is presented and shown that this state is globally stable. Numerical simulations are carried out to justify the effect of control strategies on the prevalence period of AVL.
The paper is organized as follows. Some basic ideas of the related problem are presented in Section 2. A mathematical model of the interaction between the human and the vector is presented in Section 3. Section 4 presents stability and sensitivity analysis of the proposed model. Section 5 is concerned with numerical simulations of the control strategies, with quantifications of transmission and intervention parameters. Finally, we present our conclusion.

Preliminaries
In this section, we present some definitions related to the present work [5,12,23].
where (H 1 ) : The system is defined on a positively invariant set X of the nonnegative orthant. The system is dissipative on X. That X is the region, where the model has biological sense, well posed and all its trajectories are forward bounded.

Mathematical formulation
Stauch et al. [23] studied the biological behaviour of the disease, while we present a mathematical formulation of the model. Our newly developed model represents the dynamics of visceral strains of Leishmania in two different populations, human population N h and sandfly population N v . The total population is distributed in 14 compartments:           S v : The susceptible sandfly population.
E v : The exposed sandfly population.
The infectious sandfly population.

The total human population
The susceptible humans become latent at the rate λ h , after contact with infected sand fly. After completing sojourn time, the patient develop second latent stage E 2 , at the rate ξ 1 . The fraction γ 1 ξ 2 of the patients at E 2 develops symptomatic Kala Azar, the fraction γ 2 ξ 2 enter the dormant period of PKDL and the fraction γ 3 ξ 2 recover and get PCR − . The fraction α 3 of these recovered individuals become (P − , D − , L + ). Some of these recovered individuals lose their cellular immunity at the rate α 1 , getting (P − , D − , L − ) and become susceptible again.
After diagnosis, the symptomatic infectious individuals at I 1 are put on first-line treatment at the rate ξ 3 . After the first-line treatment, these infectious individuals can be divided into six subgroups: • Fraction μ k of these infectious human dies due to toxicity of medicines.
• Fraction μ h dies due to natural death. • Fraction p 1 τ 1 does not show +ve response to first-line treatment and is put on secondline treatment. • Fraction p 2 τ 1 enters the dormant period before the development of PKDL. • Fraction p 3 τ 1 is recovered.
After completing a sojourn time at R 3 , the person further improves and become (P − , D − , L + ), and enters the class R 2 . After second-line treatment, infectious individuals can be divided into three subgroups: • A group passes away due to, toxicity of second-line treatment, μ T2 , natural mortality μ h and disease death μ k . • A group p 5 τ 2 is recovered and enters the recovered class R 3 .
• A group p 4 τ 2 enters the dormancy period of PKDL.
After completing the sojourn period (dormancy period) α 2 at E 3 , the victimized individuals develop PKDL, subject to their survival. Some of these individuals are recovered due to treatment (180 days = τ 3 ) and some recover naturally at the rates β 2 and β 1 , getting (P − , D − , L + ). After completing the sojourn period α 4 , these recovered individuals enter the recovered class R 2 . A susceptible vector, after contact with a person in latent or infectious states, gets infected and enters the exposed class E v at the rate λ v . After completing the incubation period σ v at E v , the vector becomes infectious and enters the class I v . Figure 1 represents the flow of the disease in susceptible population.
The dynamical system for human and vector population is given bẏ Here, The description of the parameters is given in the following table:

Replace model analysis by model analysis
In this section, we discuss invariant region, the disease-free equilibrium point and reproduction number R 0 of the system (1).

Invariant region
The proposed model is concerned with living population; therefore the state variables are nonnegative. The dynamic of overall population is obtained by adding all the classes concerned with humans (Ṅ h ) and adding all the classes concerned with vector (Ṅ f ) and is given by the following differential equations: where δ 1 = μ k , δ 2 = μ k + μ T1 and δ 3 = μ k + μ T2 . If the human population is diseasefree, then This shows that the biological feasible region is given by which is a positively invariant domain. The model is epidemiologically and mathematically well posed [6] and all the trajectories are forward bounded.

Disease-free equilibrium and reproduction number
For the disease-free equilibrium, we equate the right-hand sides of all the equations in Equation (1) to zero; also, we assume that initially, there is no infection. Then, the diseasefree equilibrium of the model (1) is The number of secondary infections caused in completely susceptible population by introducing an infectious individual to the population is called reproduction number R 0 [4].
In order to find the basic reproduction number, we use the next-generation method [5].

Biological interpretation and sensitivity of R
where a is the sandfly biting rate, b is the transmission probability of VL infection to human from sandfly and c 1 is the transmission probability of VL infection to sandfly from human in state E 1 . If human is susceptible and the sandfly is infected with VL, then the term ab confirms the transmission of VL infection from sandfly to human. If human is in the latency period, stage E 1 and the sandfly is susceptible, then ac 1 confirms the transmission of VL infection from human to sandfly. So, the term R a denotes VL transmission, between human and sand fly.
where c 2 is the transmission probability of VL infection to sandfly from human in state E 2 . If human is susceptible and the sandfly is infected with VL, then the term ab confirms the transmission of VL infection from sandfly to human. If human is in the latency period, stage E 2 and the sandfly is susceptible, then ac 2 confirms the transmission of Vl infection from human to sandfly. So, the term R b denotes Vl transmission, between human and sandfly. Similarly, the term R c denotes Vl transmission, with involvement of state I 4 . Thus, R 0 is biologically meaningful.

Sensitivity analysis of R 0 Definition 4.1 ([18]):
The normalized forward sensitivity index of a variable, x, that depends differentiably on a parameter y is defined as To reduce the rate of disease transmission, it is important to know the role of different parameters involved in its transmission. Since initial disease transmission depends on basic reproduction number R 0 . Therefore, we find the sensitivity indices of the parameters involved in reproduction number R 0 . These indices allow us to measure the relative change in R 0 with the change in a parameter. With the help of these indices, we find the parameters that are highly effective in disease transmission, and need to be targeted by intervention strategies. Table 1 shows the sensitivity indices of the parameters involved in the initial disease transmission. Sandfly biting rate a and duration of feeding cycle β have got highest sensitivity indices 1. This means that the decrease in biting rate by 10% would decrease R 0 by 10%. The second highest index −1 is that of sandfly 's mortality rate μ v . That is, increasing μ v by 10% will decrease R 0 by 10%. The fraction of exposed human in the class E 2 which develop symptomatic kala azar is denoted by ω 1 ·ω 1 has got a sensitivity index of −0.65. The transmission probabilities of infection between human and sandfly have got a sensitivity index of 0.5. The parameter v ; birth rate of sandflies, have got a sensitivity index of 0.5. Decrease of 10% in the birth rate of sandflies will decrease R 0 by 5%. h have got a sensitivity index of −0.5. Increase in human's birth rate causes a decrease in R 0 . The sensitivity index of death rate of human μ h is 0.49. Increase in treatment rate of human will cause a decrease in human's death rate, which will reduce R 0 . We develop four different control strategies which touch, directly or indirectly, the parameters effecting the initial transmission rate R 0 of ACL.

Stability of disease-free state
For global stability of the disease-free equilibrium, we proceed as follows: Let

Proof:
The above sub-system in accordance with Theorem 2.1 is equivalent tȯ The system can be written as follows, on the domain G = {X ∈ ; X I = 0, X s = 0} : The system (4) is a linear system. This system is globally asymptotically stable at the equilibrium ( h /μ h , 0, 0, 0, v /μ v ), corresponding to the disease-free equilibrium where the hypotheses H 1 and H 2 are satisfied.

Theorem 4.2:
For the sub-system (5), A 2 is Metzler and irreducible ∀X ∈ , and there exists a matrixĀ 2 such that where α is the stability modulus ofĀ 2 .

Proof:
We can write the sub-system (5) aṡ Here, A 2 (X) is given by the following matrix: All the off-diagonal entries of the matrix A 2 (X) are nonnegative, on the domain G. Hence, A 2 (X) is Metzler and irreducible ∀X ∈ G.
The upper bond matrix of the matrix A 2 (X) is denoted byĀ 2 and is given bȳ This upper bound is not attained in , and particularly not realized for the Jacobian at the disease-free equilibrium. Thus, we can have sufficient condition. So H 4 of Theorem 2.1 equivalent to Equations (6) and (7) holds.
Finally, we show that H 5 or Equation (8) holds; we state the following theorem.

Theorem 4.3: The Metzler matrix satisfy the axiom 'α(Ā
where is the additional threshold number given by Proof: We decompose the matrixĀ 2 in four blocks, such that where B,C,D and E are 7 × 7, 7 × 2, 2 × 7 and 2 × 2 sub-matrices, respectively. We know that the matrixĀ 2 is stable if and only if B and E − DB −1 C are Metzler stable [12]. Here, B is Metzler stable because all its off-diagonal entries are nonnegative, and all the eigenvalues are negative. To show that E − DB −1 C is stable, we proceed as follows: Let From [15], α(Ā 2 ) ≤ 0 only if Here n 1 = ab, n 2 = ac 1 , n 3 = ac 2 , n 4 = ac 3 , n 5 = ac 3 , n 6 = ac 3 , n 7 = ac 4 . By putting these values in Equation (9), we have a 2 bc 1 σ v a 1 a 13 μ v + a 2 bc 2 ξ 1 σ v a 1 a 2 a 13 μ v + a 2 bc 4 a 8 k 2 σ v ξ 1 a 1 a 2 a 11 a 12 a 13 μ v < 1.
Thus, we have shown that axioms H 1 . . . H 5 of Theorem 2.1 do hold. Now, we are in a position to claim the following result.

Theorem 4.4:
If the parameters used in the model satisfy the condition α(Ā 2 ) ≤ 0, then the disease-free equilibrium of the system (1) is globally asymptotically stable.

Simulation results of the model
We use four different control strategies and generate numerical simulation for each, using Matlab software. In these strategies, we focus on treatment of human class and control of sandfly class. τ 1 , τ 2 , τ 3 , ξ 3 are interventions by treatment of human class, where ν 1 , ν 2 are vector-related interventions ( Table 2). The ratio of human and sandfly is taken as 100:527.
In the following numerical simulations, E 1 denotes early latent period of VL in human population, E 2 is the second stage of VL-latency in human population, E 3 is the dormant period before development of PKDL, I 1 is the early symptomatic and diagnosis stage of VL in human, I 2 is the second stage of symptomatic kala azar, I 3 is the third stage of symptomatic kala azar and I 4 denotes the human class with PKDL. E f is the class of exposed       sandflies and I f denotes the infectious class of sandflies. We apply control strategies to these infectious classes to eliminate the disease. Figures 2-5 present the time spent in eradication of each infectious class.
In Table 3, we compile the summary of the results of the control strategies, obtained from numerical simulations. The table presents the time spent (TS) in elimination of each infectious class. (E 1 , TS) means the time spent (unit days) in elimination of exposed class E 1 .

Discussion and conclusion
In this work, a mathematical model of Anthroponotic leishmania transmission was developed. On the basis of sensitivity analysis of the reproduction number, we presented four control strategies. For quantification of prevalence period of the disease, we performed numerical simulations. The results shown that the disease can be eradicated by using the proposed strategies. Control strategy 2 and strategy 3 take comparatively small time in the elimination of the disease. Since the prolonged prevalence of PKDL may cause the new outbreaks of VL, we recommend control strategy 2, where the prevalence period of PKDL is comparatively low. The reproduction number of the model is most sensitive to a, sandfly biting rate; β, the feeding period of sandfly and μ v , the mortality rate of sandfly. So along with treatment of human's infectious classes, we need to focus control variables ν 1 and ν 2 , using different measures to control phlebotomize sandflies and its biting rate. This includes residual spraying of dwellings, insecticide-treated nets and application of repellents/insecticides to skin. Sandfly is susceptible to all the major insecticidal groups. The main intervention for vector control is reduction in sandfly life expectancy, which lowers the size of vector population and hence reducing the biting pressure of vector on humans. By reducing sandfly life expectancy, the vector is less likely to survive long enough to bite twice -once to acquire infection and again to infect a host, and the vector spends small time as infected. The extinction of PKDL states, I 4 and E 3 , takes comparatively a long time. Although, both, theorems and numerical simulations agree with the global stability of disease-free equilibrium yet the long prevalence of state I 4 in the population works as a reservoir and hence cannot be neglected.