Stability analysis of leishmania epidemic model with harmonic mean type incidence rate.

We discussed anthroponotic cutaneous leishmania transmission in this article, due to its large effect on the community in the recent years. The mathematical model is developed for anthroponotic cutaneous leishmania transmission, and its qualitative behavior is taken under consideration. The threshold number R 0 of the model is derived using the next-generation method. In the disease-free case, local and global stability is carried out with the condition that R 0 will be less than one. The global stability at the disease-free equilibrium point has been derived by utilizing the Castillo-Chavez method. On the other hand, at the endemic equilibrium point, the local and global stability holds with some conditions, and R 0 is greater than unity. The global stability at the endemic equilibrium point is established with the help of a geometrical approach which is the generalization of Lyapunov theory, by using the third additive compound matrix. The sensitivity analysis of the threshold number with other parameters is also taken into account. Several graphs of important parameters are discussed in the last section. © Società Italiana di Fisica and Springer-Verlag GmbH Germany, part of Springer Nature 2020.

Introduction

The main cause of leishmania type parasites is leishmaniasis. Four types of leishmaniasis have been studied in the literature [1, 2]. Among them is cutaneous leishmaniasis (CL) which is considered to be the most widespread form of leishmaniasis caused by L.tropica and L.major. The macrophagic cells of the host are attacked which results in lesion skin on uncovered parts of the body such as the legs, arm, and face. The parts of the world which are endemic to this disease are South America, Europe, and Asia. The disease is declared as one of the most dangerous diseases in the world by the WHO as the spreading rate of new cases per year is very fast (more than a million). The ten countries which together account for 70–75% of the global estimated CL incidence are Peru, CostaRica, North Sudan, Ethiopia, Syria, Iran, Brazil, Colombia, Algeria, and Afghanistan [3].

Bloodsucking and flies are the main transmitters of this disease. Throughout the world, 700 species have been spotted out, in which 37 species are recognized in Pakistan. Recently, in Pakistan cutaneous leishmaniasis (CL) spreads in the North-West part of the country, which in results killed a lot of people. Khyber Pakhtunkhwa (KPK) one of the provinces of Pakistan has been widely affected by CL, particularly the tribal areas. The same province was previously hit by CL in 1997 where an afghan refugee camp was situated [4]. At that time, Kabul (Afghanistan) was badly affected by the CL and due to cross-border movement the infected migrant carriers were the main cause of the epidemic in the refugee camp.

The transmission of parasites is carried out by the species of genera Phlebotomus type. The sand flies are the main transmitter of these parasites. The habitat of these flies enjoys a wide range i.e from desert to tropical rain forest. Not only this but also have several hosts in which dogs, chickens, humans, mammals, livestock, and vertebrates are considered to be the main hosts [5]. The color of these sand flies is sandy and normally they are 2–3 mm long. The latent period of these sand flies is considered to be between three and seven days [6]. From the ground level, the sand flies attain 2.51 m (8.3 ft) as it maximum height [7]. The incidental transference risk of blood-borne pathogens of humans is reduced by the use of animal blood and is considered to be cheaper than the conservation of animals and their preparation for feeding sand fly [8]. The fecundity of species is badly affected by the extreme temperatures [9]. Some recent studies on Cutaneous leishmaniasis epidemic models are described in [4, 5, 10–14].

In order to approach the infected population to extinct rapidly, we introduce the harmonic mean type of incidence rate between susceptible humans and infected vectors, and also between susceptible vectors and infected humans. Indeed mean of two values is the measure of centrality of a set of data. Geometric mean is primarily used to average ratios or rates of change of data. As for as harmonic mean is concerned, it is less sensitive for a few large values as compare to arithmetic or geometric means. It is sometimes used for highly skewed variables. The harmonic incidence rate shows the prospect of having the population approaches to extinction in a finite time but more rapidly as compared to other incidence rates. As compared to others, the harmonic mean is a better incidence rate when the number of individuals is defined in relation to the same unit. Moreover, in many situations, the harmonic incidence rate provides a true average of the rates and ratios. It is particularly sensitive to a single lower-than-average value.

The summary of our article is, to elaborate the CL epidemic model in Sect. 1. In Sect. 2, the CL epidemic model is formulated and the threshold value of the model is achieved by using next-generation method. If the human is susceptible and the sand fly is infected with CL, then the biting of humans by the sand flies would result in the transmission of CL stains to the humans. The direction of transmission is denoted by the term in the basic reproduction number. If the sand fly is not infected and the human is infected with CL, then clearly in the basic reproduction number is rightly indicating the secondary infections to sand fly from humans. So basic reproduction number represents the transmission of CL strains between humans and sand flies. This shows that the obtained basic reproduction number for our model is biologically meaningful. Some conditions are imposed on the threshold value to show the local and global stabilities of our proposed model in Sects. 3 and 4. The sensitivity analysis of the reproduction number is presented in Sect. 5, and the most sensitive parameters are highlighted. On the basis of sensitivity analysis, control strategies can be introduced in the model. These strategies will 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. Numerical simulations are carried out with the help of the Runge–Kutta fourth-order procedure in the last section.

Model formulation

Mathematical modeling of epidemic diseases have been widely studied by researchers [15-18]. In this section, CL epidemic model is presented consisting of seven classes. These classes are classified into four human population subclasses i.e , , , and , denoting the susceptible, exposed, infected, and recovered people. The three vector population subclasses i.e , , and represent susceptible, exposed, and infected vectors. represents the total human population i.e .

is the new recruitment of humans to susceptible class . The infected vector i.e sand fly bites the susceptible humans, and they are infected at the rate of , where transmission probability of CL to human from the sand fly is denoted by , and the sand fly biting rate is represented by a. shows the rate due to which the uninfected exposed humans are recovered and is the rate due to which the exposed humans get infectious. is the rate that shows the natural recovery of humans from infected class. The infected humans after being bitten by the sand flies infect the sand flies at the rate . is the natural death rate in humans, while shows the natural death rate in sand flies. The CL probability of transmission from humans to sand flies is represented by . is the rate which shows transmission of sand flies exposed class to infected class. The model is given below:withNext, for the system (1), we establish some basic results.

Basic properties of the model

Total population dynamics is represented by:The feasible region (biological) isFrom Eqs. (2) and (3), we obtainWhich shows that the model is well posed and is positively invariant domain.

Lemma 1

The orthant is invariant positively for the system described by (1).

Proof

Let and assumeSystem (43) is expressed in the following form:whereHere, we see that L is the Metzler matrix as it has nonnegative entries on its off diagonal and . Hence, it is concluded that the system (1) is positively invariant in .

Lemma 2

Solutions (if exist) of system (1) are positive under the initial conditions (2), for all .

Let us assume that the solutions exists in I, for all . Consider the second equation of (1) and the solution of it has the following formWe also take the third equation and solution of it has the following formClearly, it can be seen from the above solutions that these are strictly positive. In the same fashion, it can be shown that , and possess nonnegative solutions.

Basic reproductive number

The disease-free equilibrium point of system (1) is denoted by , i.eThe threshold value commonly known as basic reproduction number is very important and informative regarding the spread of infectious disease. First of all, we need to calculate this threshold value, for the sake of this let is our infected compartment, then it follows from system (1) that:Deploying the method of next-generation matrix, we have the following;whereThe inverse of V isThus,The dominant eigen value gives us , i.e

Local stability

We establish the local stability of the system (1) in this section at disease-free point as well as at endemic equilibrium point .

At disease-free equilibrium point

Theorem 3.1

The disease-free equilibrium point of the system (1) is stable locally asymptotically if .

The Jacobian matrix of the system (1) at is given byClearly, three eigenvalues of the Jacobian matrix of model (1) around the disease-free equilibrium are negative, i.e., , and . For the remaining, we taking the following reduced matrix:Using the elementary row operation, takes the following form:where , , and .

Clearly, the three eigenvalues of are negative, the last eigenvalue also has negative sign if implies that . The system (1) is locally stable at disease-free equilibrium point, if .

At endemic equilibrium point,

system (1) is rearranged to get and in terms of . Thus,

Theorem 3.2

If , then the endemic equilibrium of the system (1) is locally asymptotically stable.

The Jacobian matrix of the system (1) at is given byWhere .

Clearly, one eigenvalue of the Jacobian matrix of model (1) around the disease present equilibrium point is negative, i.e., . For the remaining six eigenvalues, we take the following reduced matrix:Using the elementary row operation, takes the following form:whereThe eigenvalues of takes the following form:The last eigenvalue also has negative real part ifthis implies that , finally we get . As we have shown that the real part of all eigenvalues are negative, hence endemic equilibrium point is locally asymptotically stable under the condition that .

Global asymptotic stability

For the model (1), the global stability at disease-free point is achieved by taking into account the Castillo–Chavez approach [19]. The method is summarized as, the proposed model (1) is reduced into the following two subsystems given byIn the system (26), and represent the number of uninfected and infected individuals, respectively, that is, and . The disease-free equilibrium is denoted by and define as Thus, the existence of the global stability at disease-free equilibrium point depends on the following two conditions At second condition, is an M-matrix that is the off diagonal entries are positive and is the feasible region. Then, the following statement holds.

If is globally asymptotically stable.

where for

Lemma 3

For , then the equilibrium point of the system (1) is said to be globally asymptotically stable, if the above conditions are satisfied [19].

Theorem 4.1

If , then the proposed model (1) is globally asymptotically stable at disease-free equilibrium and unstable otherwise.

Let and and define , whereBy using model system (1), we haveFor , and , we getThus, from equation (29) as , . So is globally asymptotically stable. Now,As and , hence . Clearly, B is M-matrix and hence both the conditions are proved, so by Lemma 1, the disease-free equilibrium point is stable globally asymptotically.

At endemic equilibrium point

For the global stability of (1) at endemic equilibrium , we use the geometrical approach [20]. The method is summarized as To investigate the sufficient condition through which the is globally asymptotically stable, consider the differential equationwhere the open set is simply connected and is a function such that . Assuming that is the solution of Eq. (32) and for , the following are true. The solution is said to be globally asymptotically stable in U,  if it is locally asymptotically stable and all trajectories in U converges to the equilibrium For a condition satisfied for f,  which precludes the existence of non-constant periodic solution of equation (32) known as Bendixson criteria. The classical Bendixson criteria for is robust under . Furthermore, a point is wandering for Eq. (32), if there exist a neighborhood N of and , such that is empty for all . Thus, the following global stability principle is established for autonomous system in any finite dimension.

There exist a compact absorbing set

System (32) has a unique equilibrium.

Lemma 4

If the conditions (a)–(b) and Bendixson criterion are satisfied for Eq. (32) [i.e., robust under local perturbation of f at all non-equilibrium, non-wandering point for Eq. (32)], then is globally asymptotically stable in U provided it is stable.

Define a matrix valued function P on U byEquation (33) is a matrix valued function on U. Further assume that exist and is continuous for Now, define a quantity define, such thatand be the third additive compound matrix of J, i.e., and . Let be the Lozinskii measure of the matrix B with respect to the norm in defined byHence, if , which shows that the presence of any orbit that give rise to a simple closed rectifiable curve, such as periodic orbits and heterocyclic cycles.

Lemma 5

Let U is simply connected and the condition (a)–(b) are satisfied, then the unique equilibrium of equation (32) is globally asymptotically stable in U, if [20].

Now, we apply the above techniques to prove the global stability of model (1) at endemic equilibrium. Thus, we have the following stability

Theorem 4.2

If , and , then the model (2) is globally asymptotically stable at endemic equilibrium and unstable otherwise.

Consider the subsystem of (1),For a matrix J,the third additive compound matrix is given byLet J be the Jacobian matrix of the system (36) given byThe third additive compound matrix of J is:whereConsider , such that and time derivative is .

Therefore,andSo that ,where;Consequently,and if , thenSimilarly,and if , thenNow,andThus, combining the above four inequalities, we get the following inequality:where is the Lozinskii measure. The subsystem of model (1) containing four nonlinear differential equations is globally asymptotically stable. Solving the remaining three linear differential equations results in , and as . Hence, is globally asymptotically stable.

Sensitivity analysis

Determining the parameters which 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 model. Sensitivity analysis of has received much attention from the ecologist and epidemiologist.

Definition 1

The normalized forward sensitivity index of the 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 (iii) by linearizing system (1), and then solving the obtain set of linear algebraic equations. We will apply the direct differentiation method as it gives analytical expressions for the indices. The indices not only shows us the influence of various aspects associated with the spreading of infectious disease but also gives us important information regarding the comparative change between and different parameter. Consequently, it helps in developing control strategies. Table 1 shows that the parameters a, , , , and have a positive influence on the reproduction number , which describe that the growth or decay of these parameters say by 10% will increase or decrease the reproduction number by 10%, 4.9%, 1.0%, 5.0%, 0.17%, respectively. But on the other hand, the index for parameters , , and illustrates that increasing their values by 10% will decrease the values of reproduction number by 0.04%, 3.6%, 5.1%, and 2.2%, respectively. and have no impact on the reproduction number.

Table 1

Comparison of sensitivity indices of the reproduction number against mentioned parameters

Parameter	S.Index	Value	[21]
a	S_a	1.000000000	1.0
θ	S_θ	-0.004792332271	-0.021
c_1	S_c_1	0.5000000001	0.5
μ_h	S_μ_h	-0.3685739184	0.4958
β	S_β	-0.2272727272	-0.4965
Γ_v	S_Γ_v	0.0	0.5
b_1	S_b_1	0.4999999999	5.0
k_1	S_k_1	0.1006389778	0.2345
k_2	S_k_2	0.01754385961	0.00009
μ_v	S_μ_v	-0.5175438600	-0.7571
Γ_h	S_Γ_h	0.0	-0.5

Sand fly biting rate a has got the highest sensitivity index, i.e., 1, while has the second-highest sensitivity index, i.e., . The third highest sensitivity index is of which is 0.5. Since the effect of all parameters is coupled with these three key parameters i.e sand fly biting rate, transmission of CL from humans to sand flies, and the mortality rate of sand flies. So instead of addressing all the parameters, one can address only the three key parameters, which are the main cause of transmission. An increase or decrease in these key parameters causes a change in the rest of the parameters in the form of an increase or decrease. For example, decrease in the sand fly biting rate a means decrease in the contact rate of humans and sand flies. This creates difficulties for female sand flies to have human blood, which they needs for laying eggs. Consequently, a decrease occurs in the sand fly birthrate . While a decrease in contact rate of sand fly and human reduces the chances of sand fly to catch infection from human or to transmit the infection to human. This will reduce the transmission probability and of CL between humans and sand flies (Fig. 1).

Fig. 1

Sensitivity analysis of different parameters

Numerical simulations and discussion

The RK4 method is a fourth-order method, meaning that the local truncation error is on the order of , while the total accumulated error is on the order of , where h is step-size. Estimating the error has little or negligible computational cost compared to a step with the higher-order method. In most situations of interest, a fourth-order Runge–Kutta integration method represents an appropriate compromise between the competing requirements of a low truncation error per step and a low computational cost per step.for usingWe used initial condition of the state variables as and The Runge–Kutta method is popular because of its simplicity and efficiency. It is one of the most powerful predictor-correctors methods, following the form of a single predictor step and one or more corrector steps. The fourth-order Runge–Kutta approximation is given by

In this paper, we established global dynamics and sensitivity analysis of the anthroponotic cutaneous leishmania epidemic model. The sharp threshold parameter i.e basic reproduction number totally establishes the global stability of the proposed model in Theorems 4.1 and 4.2. Theorem 4.1 guarantees the global stability in the disease-free equilibrium case with the condition that the threshold parameter will be less than or equal to one. On the other hand, Theorem 4.2 proves the global stability at the endemic equilibrium point with the condition that the threshold parameter will greater than one. The global dynamics have been studied by utilizing the Castillo-Chavez method and geometrical approach. Now, we provide some simulations of anthroponotic cutaneous leishmania epidemic model to attest our findings. Choosing parameters as mentioned in Table 2, the basic reproduction number , hence by theorem 2.1 the disease-free equilibrium is globally asymptotically stable, and , , , , , , and tend to their disease-free equilibrium point (see Figs. 2, 3).

Table 2

Parametric values of model (1) used for simulation

Parameter	Value	Parameter	Value
a	0.012	b_1	0.071
θ	0.32	k_1	0.23
c_1	0.41	k_2	0.19
μ_h	0.21	μ_v	0.75
β	0.49	Γ_h	0.4
Γ_v	0.6

Fig. 2

, , , , , , and tend to their disease-free equilibrium point when

Fig. 3

, , , , , , and are unstable when

Conclusions

Recently, Covid-19 which was initially an epidemic but quickly become pandemic broke out in China. Covid-19 is locally as well as globally unstable till March, 2020. Similarly, from time to time cutaneous leishmaniasis (CL) also adopts the shape of the epidemic, particularly in Pakistan. That is why it is very necessary to model cutaneous leishmaniasis and discuss its local as well as global stability. Important parameters are highlighted which is sensitive to a threshold value commonly known as basic reproductive number. For endemic stability analysis, we consider the generalization of the Lyapunov method called a geometrical approach in which the third additive compound matrix is taken into account. The feasibility of our result is verified by numerical simulations.

Extending our work, one can use harmonic mean type incidence rate to reformulate the visceral leishmaniasis epidemic model. One can check its stability analysis, the sensitivity of parameters, bifurcation analysis, and optimal control.