Sensitivity Analysis and Optimal Control of Anthroponotic Cutaneous Leishmania.

This paper is focused on the transmission dynamics and optimal control of Anthroponotic Cutaneous Leishmania. The threshold condition R0 for initial transmission of infection is obtained by next generation method. Biological sense of the threshold condition is investigated and discussed in detail. The sensitivity analysis of the reproduction number is presented and the most sensitive parameters are high lighted. On the basis of sensitivity analysis, some control strategies are introduced in the model. These strategies positively reduce the effect of the parameters with high sensitivity indices, on the initial transmission. Finally, an optimal control strategy is presented by taking into account the cost associated with control strategies. It is also shown that an optimal control exists for the proposed control problem. The goal of optimal control problem is to minimize, the cost associated with control strategies and the chances of infectious humans, exposed humans and vector population to become infected. Numerical simulations are carried out with the help of Runge-Kutta fourth order procedure.

Introduction

Anthroponotic cutaneous leishmaniasis (ACL) is an infectious disease caused by Leishmina Tropica. The parasite is transmitted by the sand fly Phlebotomus sergenti; about 2–3mm long, and is sandy-colored. These flies are blood-feeding and breed in forest areas, caves, or the burrows of small rodents. The sand fly latent period is assumed roughly to be 3–7 days [1]. Leishmaniasis is a group of infectious diseases. utaneous leishmaniasis(ACL), Muco-cutaneous leishmaniasis(MCL), Visceral leishmaniasis (VL) or kala-azar, and (PKDL) Post-kala-azar dermal leishmaniasis are four main clinical syndromes of Leishmaniasis [2, 3].

Cutaneous leishmaniasis (ACl, ZCL) is the most common form of leishmaniasis. The incubation period of Cutaneous leishmaniasis is two to eight weeks, however, in some cases longer periods have been reported. Leishmania major is the causative agent of human Cutaneous Leishmaniasis. Recovery from infection of Leishmania major, results in persistence of parasites at the primary infection site. Which causes long-lasting immunity to reinfection. Vaccination with killed parasites or recombinant proteins induces only short-term protection [4, 5]. Cross-immunity experiments between different species of leishmania is discussed in [6].

In Pakistan, ACL has the widest distribution, occurring in urban areas of Punjab (Multan), Balochistan (Quetta) and in the Northern Areas and Azad Kashmir. Cases of ACL are being increasingly reported from towns of Khyber Pakhtunkhwa. In 1997, a large outbreak of ACL was reported from Timargara refugee camp in Khyber Pakhtunkhw [7].

Mathematical model plays an important role in analyzing the transmission dynamics and control of infectious diseases. The formulation process of a mathematical model elucidates assumptions, parameters and variables. Chaves et al. [8] presented a mathematical model of American Cutaneous Leishmaniasis. They discussed the threshold conditions for infection persistence of the disease. In their recent work, Chaves et al. investigated the effect of deforestation on vector borne diseases [8, 9]. Das et al. studied the effect of delay on the dynamics of Cutaneous Leishmaniasis [10]. Calzada et al. investigated the compositional changes that occur in sand fly after the use of thermal fogging [11]. Chaves [12], Bacaer and Guernaoui [13] discussed the seasonal models of Cutaneous Leishmaniasis. Several strategies are adopted to control the transmission and eradication of different infectious diseases. However, it has been always a challenge to establish a balance between the cost and control of a strategy. Optimal control theory plays extra ordinary role in keeping the minimum cost for maximum efficiency of the control strategy. The optimal controls of several epidemic models are described in [14-17].

In this paper, we present a mathematical model for the transmission dynamics and optimal control of Anthroponotic cutaneous leishmaniasis. The threshold condition R0 for initial transmission of infection is obtained by using Next generation method. We investigated and discussed in detail different scenario of biological sense of the threshold condition. The sensitivity analysis of the reproduction number is presented and the most sensitive parameters are high lighted. On the basis of sensitivity analysis, some control strategies are introduced in the model. These strategies reduce the effect of the parameters with high sensitivity indices, on the initial transmission. The model will then be used to determine the cost-effective strategies for eradicating the disease transmission. In order to do this, an optimal control strategy is presented by taking into account the cost associated with control strategies. It is also shown that an optimal control exists for the proposed control problem. The goal of optimal control problem is to minimize, the cost associated with control strategies and the chances of infectious humans, exposed humans and vector population to become infected. Numerical simulations are carried out with the help of Runge-Kutta fourth order procedure. To address the effect of these parameters, we introduce some control variables in the model.

The paper is organized as follows. A mathematical model of the interaction between the human and vector is presented in Section 2. Section 3 represents mathematical and sensitivity analysis of reproduction number R0. Section 4 is concerned existence of control strategies. In section 5, the optimality system is numerically solved. The discussion and conclusion of the work is presented in section 6.

Model Formulation

In this section, we present a formulation of mathematical model which represents the Anthroponotic Cutaneous Leishmania infection in a community. The human population consists of four subclasses, S(t) represents the class of susceptible human, E(t) is the Cl-latent class of human, I(t) is the human class of infectious with Cl, R(t) is the human class recovered individuals from Cl at time t. Thus the total human population N(t) is N(t) = S(t) + E(t) + I(t) + R(t).

The vector population is divided into three sub-classes, S(t) represents the susceptible vector, E(t) is the exposed vector and I(t) infected vector. Total vector population N is N(t) = S(t) + E(t) + I(t). The flow chart shows the transmission dynamics of the each individual in the model given in Fig 1. Humans are recruited at a constant per capita rate Γ to the susceptible class S. The susceptible humans after being bitten by infected sand fly, are infected at the rate , where a is sand fly biting rate and b1 is the transmission probability of ACL to human from sand fly [18]. The exposed humans are recovered at the rate θ without becoming infected and the rest of the exposed humans after completing latency period, get infectious at the rate k1. Some of the infected humans are recovered naturally at the rate β. Sandflies after contact with infected humans are infected at the rate , where c1 is the transmission probability of ACL from human to sand fly. After completing incubation period, the exposed sand flies get infectious at the rate k2. μ and μ are natural mortality rates of human and sandflies, respectively. The dynamical system for human, reservoir and vector population is given by:

Fig 1

The graph represents the transmission dynamics of the each individual in the model.

Analysis of the Model

In this section, we study the stability analysis of the proposed model.

Invariant Region

We assume all the parameters to be nonnegative. The model is concerned with living population, therefore the state variables are assumed to be nonnegative at time t = 0. The dynamic of overall population is given by the following differential equations: The non-negative equilibrium of the above equations are Let us consider the biological feasible region Ψ is given by

By solving Eq (2), we obtain So Similarly we solve Eq (3) to obtain Hence Ψ is positively invariant domain, and the model is epidemiologically and mathematically well posed.

Reproduction Number

The number of secondary infection that occurs in completely susceptible population if we introduce an infectious individual to the population, is called reproduction number R0 [19]. This epidemiological measurement (number) is widely used to measure the transmission potential of the disease in susceptible population [20]. Reproduction number measures the initial disease transmission in population dynamics. The disease can invade the susceptible population and persist if R0 > 1. The disease dies out and does not spread if R0 < 1 [21]. To find the basic reproduction number, we use the Next generation method see for detail [22]. R0 = ρ(−FV−1), where ρ is spectral radius. In the absence of exposed human and vector, our model becomes and where

a1 = θ + k1 + μ, a2 = β + μ, a3 = k2 + μ.

The matrix say K = FV−1 obtained at disease-free equilibrium is referred as Next generation matrix for the given system. Each entry (m, n) of the Next generation matrix represents the expected number of secondary infections in compartment m caused by members of compartment n. The dominant nonnegative eigenvalue of the Next generation matrix is called Reproduction number, denoted by R0 [19]. The dominant eigen value of ρ(−FV−1) is: Which is the reproduction number

0.1 Biological Interpretation of Reproduction Number

In order to study the biological interpretation of R0 is given by Here a is the biting rate of sand fly, b1 is the transmission probability of Cl strain to human from sand fly and c1 is the transmission probability of Cl strain to sand fly from human.

If the human is susceptible and the sand fly is infected with Cl. The biting of human by sand fly, would result in the transmission of Cl stains to the human. The direction of transmission is denoted by the term ab1 in the R0. If the sand fly is not infected and the human is infected with Cl, then clearly ac1 is rightly indicating the secondary infections to sand fly from human. So R0 represents the transmission of Cl strains between human and sand fly. Thus R0 is biologically meaningful. Furthermore, we have emphasized ac1 and ab1, because we are interested in the direction of transmission. The rest of parameters used in R0 are concerned only with the magnitude and some others behavior in R0.

Sensitivity Analysis of R0

The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter. If the variable is a differentiable function of the parameter, the sensitivity index is then defined using partial derivatives [23].

Definition The normalized forward sensitivity index of a variable u that depends differentiability on a parameter p is defined as: To reduce the rate of disease transmission, we need to know the importance of different factors involved in its transmission. Initial disease transmission depends on reproduction number R0. So we investigate the sensitivity indices of the reproduction number R0, relative to the parameters involved. The sensitivity indices of the reproduction number R0 is given in Table 1. These indices allow us to measure the relative change in R0 with the change in a parameter. Using these indices, we find the parameters that highly effect R0, and need to be targeted by intervention strategies.

Table 1

The sensitivity indices of the reproduction number R0.

Parameter	Sensitivity index	Index at parameters value
a	+1	+1
β	-0.5ββ+μ_h	-0.4965
Γv	+0.5	+0.5
Γh	-0.5	-0.5
b1	0.5	+0.5
c1	0.5	+0.5
k1	0.5(Θ+μ_h)(k_1+Θ+μ_h)	+0.2345
k2	0.5μ_hk_2+μ_h	+0.00009
μv	-0.5(μ_v+2k_2)k_2+μ_v	-0.7571
μh	0.5(βk_1+βΘ-μh2)(β+μ_h)(k_1+Θ+μ_h)	+0.4958

Optimal Control Strategies

In this section, we use optimal control techniques to develop control strategies. In order to do this, we focus to control the transmission of infection like sand fly biting rate a has got highest sensitivity index 1. This means that decreases in biting rate by 10% would decrease R0 by 10%. The second highest index −0.7571 is that of sand fly’s mortality rate μ. That is increasing μ by 10% will decrease R0 by 7.5%. The parameters Γ, b1 and c1 have got sensitivity index of 0.5. By decreasing these parameters 10%, causes collective decrease of 15% in R0. Also Γ have got sensitivity index of -0.5. It is also to be noted that increase in human’s birth rate causes decrease in R0. The absolute sensitivity index of both treatment rate of human β and death rate of human μ is 0.49. So increasing the treatment rate of human by 10 percent will decrease R0 by 4.9%. Moreover, increase in treatment rate of human will cause decrease in human’s death rate, which will reduce R0.

Since the effect of all these parameters is coupled with three key parameters sand fly biting rate, treatment rate of infectious human and the mortality rate of sand fly. So instead of addressing all the parameters, we address the three key parameters; sand fly biting rate a, treatment rate of infectious human class and the mortality rate of sand flies which are main cause of transmission. Increase or decrease in these key parameters causes change in the rest of the parameters in the form of increase or decrease. For example decrease in sand fly biting rate a means decrease in the contact rate of human and sand fly. This creates difficulties for female sand fly to have human blood, which it needs for laying eggs. Consequently decrease occurs in sand fly birth rate Γ. While decrease in contact rate of sand fly and human, reduce the chances of sand fly to catch infection from human or to transmit infection to human. This will reduce the transmission probability b1 and c1 of ACL between humans and sand flies.

In order to develop optimal control strategies, we introduce the following three optimal control variables:

The control variable u1 represents the use of insecticide-treated bed nets and the use of sand fly’s repulsive lotions and electronic devices, to reduce sand fly biting rate.

The control variable u2 shows the use of effective medicines for the treatment of infectious humans.

The control variable u3 represents different measures like residual spraying of dwellings and animal shelters to kill sand flies at all stages.

By using these control variables our control problem becomes The goal of our optimal control strategies is to minimize the infectious and exposed human population, the vector population, sand fly biting rate and the cost of implementing the control by using possible minimal control variables u1(t), u2(t) and u3(t). In order to do this, we using the bounded Lebesgue measurable control to construct the objective functional is given by subject to the system (4).

In the objective functional, g1, g2 and g3 represent the weight constants of the exposed, infectious human and of vector population, respectively. d1, d2 and d3 are weight constants of human self protection, human treatment and vector control, respectively. The terms and described the cost of disease interventions. The cost associated with the first control strategy u1 comes from the cost of sand fly repellent lotions, electric mats and mosquito bed nets. The cost associated with the second control strategy u2(t), is the cost of expensive medication of human class. The cost associated with the third control strategy u3(t), could arise from applying different types chemical pesticides, to kill sand fly at any stage of its life. We have assumed the costs as proportional to the square of the corresponding control function. We aim to find control functions so that: subject to the system (4). The control set D is defined as:

Existence of the Control Problem

For existence of the control problem, we consider the control system (4) with initial conditions at t = 0. In case of bounded Lebesgue measurable controls and non-negative initial conditions, there exist non-negative bounded solution of the state system [14, 24]. For optimal solution of the system, first we find the Lagrangian and Hamiltonian. We define the Lagrangian of the control problem Eq (4) is given by where N(t) = S(t) + E(t) + I(t).

We define the Hamiltonian H, to find the minimal value of the Lagrangian, as follow:

Theorem There exists an optimal control such that J(u1, u2, u3) = min( J(u1, u2, u3) subject to the system (4) and initial conditions at t = 0.

Proof The control variables and the state variables are nonnegative. So we use the result in [14, 25], for the existence of an optimal control. The necessary convexity of the objective functional in u1, u2 and u3 are satisfied here. Also by definition, the set of the control variables (u1, u2, u3) ∈ D, is convex and closed. The compactness which we need for the existence of the optimal control is confirmed by boundedness of the optimal system. Also, the integrand in the functional Eq (5), is convex on the control set D. We can find a constant η > 1 and positive numbers ξ1 and ξ2 such that , as the state variables are bounded. Which completes the proof of existence of an optimal control.

For optimal solution, we apply Pontryagin’s Maximum Principle [26] as follow. If (x(t), u(t)) is an optimal solution of an optimal control problem, then there exists a non trivial vector function λ(t) = (λ1(t), λ2(t), …, λ(t)) satisfying the following:

We apply the necessary conditions to the Hamiltonian H in the system (7).

Theorem let be optimal state solutions with associated optimal control variables for the optimal control problem Eqs (4) and (5). Then there exist adjoint variables λ for i = 1, 2, …, 7 satisfying with transversality conditions (or boundary conditions) Furthermore, optimal controls and are given by

Proof: To find adjoint equations and the transversality conditions, we use Eq (7). Differentiating the Hamiltonian H with respect to each state variables, we obtain the system (8). Solving the equations , and on the interior of the control set and using the optimality conditions and the property of the control space D, we can derive Eqs (10)–(12). The optimal control and the state are found by solving the optimality system (5), the adjoint system , initial and boundary conditions (9), equations of the optimal control Eqs (10)–(12). Since the second derivatives of the Lagrangian with respect to u1, u2 and u3 are positive, so optimal problem is minimum at controls , , and . By substituting the values of , , and in the control system (7), we get the following system with

Numerical Simulations

In this section, we present numerical simulation. For numerical simulation, we use fourth order Runge-Kutta forward method. First, we use fourth order Runge-Kutta forward in time to solve the system (4), with estimated controls, over the simulated time. We then solve the system (8) with the help of backward method and the condition (9), utilizing the current iteration of the state equations. We apprise the controls u1, u2 and u3 through convex combination of both the previous iteration’s controls and the value of the characterization of Eqs (10)–(12). By keeping continue in this way until the values of unknowns at the consecutive iterations are too closed. For more detail see [27]. The value of weight constants used in the objective functional are g1 = 0.081, g2 = 0.01, g3 = 0.04, d1 = 0.08, d2 = 0.01 and d3 = 0.04. We used initial value of the state variables such as S(0) = 100, E(0) = 20, I(0) = 20, R(0) = 10, S(0) = 1000, E(0) = 20 and I(0) = 30.

Since the proposed model is composed of coupled system of nonlinear differential equations. Therefore applying control strategy to any targeted class say x of the model, will also effect the rest of the classes say x, in the model. In other words x is directly effected and x are indirectly effected by the applied strategy.

Discussion and Conclusion

First, we discuss results of numerical simulations. The control variable u1 reduces sand fly biting rate a, which causes decrease in the contact rate of humans and sand flies and consequently decrease in exposed and infectious classes of human as shown in Figs 2 to 8. The graph in Fig 2 represents the susceptible humans with and without control. The number of susceptible individuals approach to a small number due to optimal control. In Fig 3 The number of infectious individuals approach to a small number due to optimal control. The exposed and recovered humans with and without control are represented in Figs 4 and 5, respectively. The graph in Fig 6 shown the susceptible vectors with and without control. The number susceptible sand flies approach to a small number, due to optimal control. Figs 7 and 8 shown the graph the exposed and infected vectors with and without control, respectively. The number of exposed sand flies approach to a small number, due to optimal control. The control variable u1 indirectly effect the infectious class of human, where the control variable u2, directly effect the infectious class of humans. This is the reason that reduction in the infectious class is rapid as compared the reduction in exposed class of human. Fig 3 shows that for the initial period of few days, an increase occur in the infectious class of human. There are two reasons for this increase:

Fig 2

The graph represents the susceptible humans with and without control.

The number of susceptible individuals approach to a small number due to optimal control.

Fig 8

The graph shows infected vectors with and without control.

The number of infected sand flies approach to a small number, due to optimal control.

Fig 3

The graph represents the infectious humans with and without control.

The number of infectious individuals approach to a small number due to optimal control.

Fig 4

The graph shows the exposed humans with and without control.

Fig 5

The graph represents the recovered humans with and without control.

Fig 6

The graph represents the susceptible vectors with and without control.

The number susceptible sand flies approach to a small number, due to optimal control.

Fig 7

The graph represents the exposed vectors with and without control.

The number of exposed sand flies approach to a small number, due to optimal control.

Applying the control u1, minimize the number of new entries to the exposed human class. But those humans, whom are already exposed, get infectious after completing incubation period, which causes increase in the infectious class of human for initial few days.

The medication u2, takes some time to show its effect, because the infectious human does not recover immediate after the use of medicine.

Recovered human class depends on both infectious and exposed human classes. Therefore an increase is observed in the recovered human class, with the decrease in both infectious and exposed human classes, as shown Figs 4 and 5.

To control the spread of vector population we apply control variable u3. Figs 7 and 8 show that reduction in vector population is very rapidly. There are three reasons for this rapid decrease:

The control variable u1 minimize the biting rate of sand fly biting rate, which effects vector population in two different ways:

Female sand fly needs human blood before laying eggs. Therefore, reduction in sand fly biting rate causes decrease in the birth rate of sand fly and consequently a decrease in the number of susceptible vectors.

The control u1, reduces the contact rate of human and sand fly. Therefore the chances of susceptible vector to interact the infectious human decreases. Hence the rate of disease transmission from human to sand fly decreases, which causes reduction in the exposed and infectious classes of vectors.

The control u2 minimize the number of infectious humans, which reduces the chances of susceptible sand fly to catch infection from infectious human. This causes decrease in the number of exposed and infectious vectors represented in Fig 9.

Fig 9

The graph shows the control variables u1 and u2.

The control variable u3, directly effects the vector population. This control tries to kill sand flies in all compartments of the vector population.

Fig 9 shows that the efficiency of the control variable u1 is lowest for initial few days and remains maximum rest of the time. Similarly the efficiency of the control variable u2 is minimum for initial few days and remains maximum rest of the time. The reason is that both the strategies; medication of infectious human and the use of sand fly repulsive lotions, take time to show effect.

The efficiency of the control variable u3 is maximum all the times represented in Fig 10. Because the strategy is due to kill the sand flies with direct source, with immediate effect.

Fig 10

The graph shows the control variable u3.

We analyzed the mathematical model of Anthroponotic Cutaneous Leishmania using sensitivity indices of the reproduction number R0. With help of these sensitivity indices we found the relative importance of the role of different parameters, in the transmission of ACL. We addressed three key parameters; sand fly biting rate, mortality rate of sand flies and healing (recovery) rate of infectious humans. For this we introduced three control variables (strategies), in the optimal control problem, each for reducing sand fly biting rate, increasing mortality rate of sand fly and increasing recovery rate of humans. Our control strategies caused decrease in initial transmission rate R0. Since the control strategies are always effected by economics constraints, therefore, we have taken into account the constraints imposed by limited resources, in our objective functional. The control system was analyzed through Pontryagin’s Maximum Principle, to determine the necessary conditions for existence of optimal control problem. The results obtained from numerical simulations show that the control strategies are very effective, if implemented, on the same time, in the same area. And hence it is concluded that proposed optimal control works well in eradication and blocking further transmission of the Anthroponotic Cutaneous Leishmania.

We would like to extend this work by incorporating different biological behaviour in the future for instance:

The impact of temperature on the birth rate of sand flies and hence the transmission rate of Leishmaniasis, and to design a time and temperature dependent mathematical model. Since temperature of different regions are different, so there should be a particular mathematical model for a particular region.

The role of cross immunity in the control of Cutaneous Leishmaniasis.