Understanding Visceral Leishmaniasis Disease Transmission and its Control — A Study Based on Mathematical Modeling

Understanding the transmission and control of visceral leishmaniasis, a neglected tropical disease that manifests in human and animals, still remains a challenging problem globally. To study the nature of disease spread, we have developed a compartment-based mathematical model of zoonotic visceral leishmaniasis transmission among three different populations—human, animal and sandfly; dividing the human class into asymptomatic, symptomatic, post-kala-azar dermal leishmaniasis and transiently infected. We analyzed this large model for positivity, boundedness and stability around steady states in different diseased and disease-free scenarios and derived the analytical expression for basic reproduction number (R0). Sensitive parameters for each infected population were identified and varied to observe their effects on the steady state. Epidemic threshold R0 was calculated for every parameter variation. Animal population was identified to play a protective role in absorbing infection, thereby controlling the disease spread in human. To test the predictive ability of the model, seasonal fluctuation was incorporated in the birth rate of the sandflies to compare the model predictions with real data. Control scenarios on this real population data were created to predict the degree of control that can be exerted on the sensitive parameters so as to effectively reduce the infected populations. OPEN ACCESS Mathematics 2015, 3 914


Introduction
Leishmaniasis is an endemic infectious disease that has a wide potential to infect large human populations on a global scale.Visceral leishmaniasis (VL), sometimes usually referred to as "kala-azar" is the most severe form of leishmaniasis, a systemic disease that is fatal, if left untreated.Though important, the World Health Organization (WHO) classifies leishmaniasis as a neglected tropical disease [1] posing a serious but important concern for researchers and epidemiologists to explore the disease, and provide practically feasible solutions to limit its spread through populations.The overall prevalence of leishmaniasis worldwide is between 10 and 12 million; out of which there is an estimated incidence of 500,000 VL cases alone per year [1][2][3].WHO reports suggest that around 90% of visceral leishmaniasis cases are highly concentrated mainly, in six countries: Bangladesh, Brazil, Ethiopia, India, South Sudan and Sudan [1].Reports also suggest that visceral leishmaniasis alone claims around 20,000 human lives per year from India, pointing out the severity of disease and the serious threat that it poses to several developing countries [4,5].
Epidemiologically, visceral leishmaniasis can be transmitted through either the anthroponotic medium, where human is the sole source of infection and disease transmission is between human to human, or the zoonotic medium, where an animal is involved in the transmission cycle and the transmission is between animal and human [6] and is usually transmitted indirectly between hosts by sandflies of the genera Phlebotomus and Lutzomyia [7].Further, the disease is also seasonal, with the infected cases fluctuating in accordance with varying vector populations [8,9].On the other hand, clinically, the VL infected human cases can be categorized as asymptomatic, symptomatic and post-kala-azar dermal leishmaniasis (PKDL)-infected representing different stages of disease diagnosis.Asymptomatic infection represents that proportion of individuals which do not display any detectable symptoms of VL infection and clinically constitute the most crucial class of individuals responsible for spreading infection as they remain untreated until appearance of symptoms.Symptomatic infected individuals display symptoms within a defined period of parasite incubation within them.After continuous treatment and subsequent recovery from a VL infection, patients start developing a less severe complication of visceral leishmaniasis, termed as post-kala-azar-dermal leishmaniasis.It is largely characterized by a macular, maculopapular, and nodular rash in a patient, who has recovered from VL and who is otherwise well [6].The PKDL-infected cases act as a constant source of VL infection spread among human populations.Both the impact of the mentioned diagnostic complications along with the diversity in mechanisms of disease spread among host populations impose the need of a deliberate effort to combat visceral leishmaniasis.Several conventional control strategies have been developed for control of leishmaniasis spread.They are broadly classified into-(a) vector-based control methods: Methods like insecticide and residual sprays, sandfly repellents, barrier for sandfly bites through appropriate netting and clothing provide control through either elimination of vector populations or barrier for sandfly bites [10]; (b) parasite-based control methods: Treatment of leishmaniasis patients with drugs, like sodium stibogluconate, amphotericin B, miltefosine, paromomycin, etc. as a primary line of treatment [11].Even though the above control strategies have been largely used separately to curb the spread of the disease, optimizing and combining different control strategies can further help in large-scale elimination of the disease [12].
A number of experimental and theoretical approaches are available to detect the spread of VL infection and provide solutions to control disease spread in the human population [11].From the 19th century onwards, mathematical and statistical models of infection transmission have provided a greater understanding of infection persistence in a population and have played an imperative role in providing solutions to curb the spread of infection [13].A number of mathematical and statistical modeling studies have been published [5,[14][15][16][17][18][19][20][21][22][23][24][25][26] to explore the dynamics of leishmaniasis disease transmission and its incidence in populations.Most of these models have focused on transmission of cutaneous leishmaniasis through different populations [14,[21][22][23][24][25].Only a few studies that explicitly model visceral leishmaniasis transmission exist [15,16,19,20,27,28].These models have attempted to study zoonotic VL disease transmission among populations.A general model for zoonotic transmission of visceral leishmaniasis through populations was developed and analyzed by Burattini et al. [17], where the sandfly, animal and human populations were considered as a SEIR type model.A region-specific SIR type model for modeling visceral leishmaniasis disease spread in Sudan was developed by ELmojtaba et al. [15], where they included one more class for PKDL individuals in addition to the infected human class considered by Burratini et al. [17].Even though these models were formulated as standard SIR/SEIR type models, they did not consider the contribution of differential clinical manifestations of the disease in the human population to the infection spread.To address this issue, two models have been proposed to understand the VL disease dynamics in the Indian subcontinent via an anthroponotic medium [19,20].Both these studies have explicitly divided the human infection class as asymptomatic, symptomatic and PKDL-infected to detect the spread of infections in human populations.However, none of the available studies have considered both the role of clinically distinct infected human classes and the animal population in VL disease dynamics.
With this background, we propose a detailed compartment-based mathematical model that studies the VL transmission process in three distinct populations-the human and the animal as hosts, and sandfly as the vector, further considering the detailed division of the human infected population into its clinically distinct classes-symptomatic, asymptomatic and PKDL infected.With respect to these population groups, the model comprises a total of 14 compartments (modelled as variables) and the rate of disease transmission as a flow from one compartment to another.The model was analyzed both analytically and numerically for this large system to demonstrate the effect of disease transmission between the heterogeneous populations and to estimate the possible control mechanism.
As a first step, positivity and boundedness analysis of the model was performed.The VL infection dynamics was investigated in the presence or absence of the constituent compartments to understand their individual contribution in disease spread.For this, the disease-free equilibrium points for different scenarios under the presence or absence of each population in the infection was calculated analytically and conditions of their stability were elucidated.The expression for key epidemiological parameter determining the persistence of infection in the population, that is the basic reproduction number (R0) [29], was derived for this large system.
The model was also verified numerically in these different scenarios to compare it with the analytical results.The model dynamics not only agrees with previously published models [17], but also provides a detailed explanation of disease transmission through the host populations.We also performed local parameter sensitivity analysis to identify the most sensitive parameters for epidemiologically important variables in the model.The sensitive parameters were then varied to understand the effect of these parameters in the steady state of each variable.The epidemic thresholds were identified from the model by calculating R0 of the system for every parameter change.With respect to these results, we discuss important control measures that could be taken in order to the curb the disease spread.Most importantly, from our analysis, we suggest that combinatorial control of the identified sensitive parameters in specific ranges is the most optimal and effective way to eliminate VL disease spread.In addition, from our model predictions, we could hypothesize that animal populations recovering from infection can be effective absorbers of infection thereby, reducing infection spread through human populations.Finally, as an assertive point, we consider the birth rate of sandfly (density) as a binormal distribution function of rainfall and attempted to fit the model outcomes to the number of cases occurred in Muzaffarpur district, Bihar, India for the year 2005 [26].The model predictions were comparable with the behavior of real population data, thereby substantiating the importance of our detailed model in capturing the disease dynamics.We further extended our parameter control analysis to real population data where we observe the effect of variation in the identified sensitive parameters on real population data substantiating their role in control of infection.

Model Construction
Visceral leishmaniasis demonstrates a complex disease dynamics by infecting heterogeneous host populations [6].Obtaining a realistic perspective of VL infection dynamics among heterogeneous populations is a challenging task.In order to model VL infection dynamics, we considered three different population groups possibly required for the infection to persist-the human, the animal and the sandfly.The model equations were derived with respect to few assumptions.The total population of each group was assumed to be constant (each group normalized by its total population) [17].Similarly, the total human infected class consisted of symptomatic infected, asymptomatic infected and the PKDL-infected individuals considering each as a separate sub-compartment [19,20].The asymptomatic infected human population represents the set of individuals that do not present any clinical/physiological symptoms of leishmaniasis [19].With respect to the immunogenic potential of individual, the asymptomatic individuals may show some symptoms adding to the symptomatic infected population or they might get recovered gradually.Accordingly, a transient stage between the asymptomatic infected and the recovered human classes, was introduced as a separate compartment to capture the delay in conversion of asymptomatic infected into the recovered human population separately.The transmission of exposed class to infected class is absolutely assumed to be dependent on the number of infectious bites provided by the sandfly to the human or animal and ability to induce an infection in human or animal.In addition, the transition of the population from one compartment to another depends on the ability of the human, animal and sandfly population to fight back the infection.
Considering the aforementioned assumptions, we constructed a compartment based model integrating the three heterogeneous population groups.Let NH (t) represent the total human population, NF (t) represent the total sandfly (vector) population and NA(t) the total animal population.To observe the changes in the relative proportion of individuals in each stage/compartment for each population group during transmission of disease, we considered NH (t) = NF (t) = NA(t) = 1.In addition, the model parameter values have been scaled accordingly to capture changes in relative population fractions.To represent a biologically realistic complex scenario among populations, we have further divided the total human population into seven compartments namely, susceptible SH(t), exposed EH(t), asymptomatic infected IHa(t), symptomatic infected IHs(t), PKDL infected IHp(t), transient stage population TH(t) and recovered populations, those who are permanently immune to the infection RH(t).Similarly, the vector population is divided in to three compartments, susceptible sand-flies SF(t), latent sand-flies EF(t) and infected sand-flies IF(t), and the animal or reservoir populations into four compartments, susceptible animal SA(t), latent animal EA(t), infected animal IA(t) and recovered animal RA(t).The TH or the transient class represents a fraction of asymptomatic individuals that display symptoms after a delayed period of time and are subsequently treated to become recovered [19].The individuals belonging to this class demonstrate recovery solely due to the primary and secondary lines of treatment that attempt to eliminate the parasite population within each host.The flow of infection between these compartments in the human and animal populations have been modeled using the standard SEIR model whereas the flow of infection through the vector has been modeled using the SEI carrier type model [30].Figure 1 represents the model structure encompassing various components considered in the model, each arrow between the compartments denoting the flow of infection from one compartment to another.The arrows providing to each susceptible population denote the birth rate of individuals αH, αA, αF and arrows leaving from each compartments denote the mortality rate µH, µA, µF of each individual compartment considered in the model.pH and pA denote the probabilities that any sandfly will bite a susceptible human or animal, with respect to the availability of the corresponding host blood meal.Further, the interactions of susceptible animals and humans with the fly populations have been modeled using standard incidence function, considering the probability that an infected sandfly would transmit infection to susceptible human or animal host (pHi, pAi).A latent period of infection (incubation period) is considered between the exposed and infected individuals of each population (σHs, σHa, σHp, σA, and σF) to differentiate between the infectious (exposed) and the infected.We also assume that a fraction of exposed individuals become either symptomatic (bHs) or asymptomatic (1-bHs).As mentioned before, the asymptomatic individuals become symptomatic and the symptomatic individuals after treatment might become PKDL-infected.To include this, we assume that a fraction of asymptomatic individuals (bHt) remain dormant representing the transient population and the remaining become symptomatic (1-bHt).Similarly, a fraction of symptomatic individuals after treatment become PKDL-infected (bHp) and the remaining recover from the infection (1-bHp).The corresponding rates of transition between the asymptomatic to symptomatic and symptomatic to PKDL-infected are given by ρHs and ρHp respectively.Further, each human and animal infected population is assumed to recover from the infection with a specific recovery rate (γHt, γHs, γHp, γA).The transition rate of the infection from the asymptomatic infected to the transient stage (γHs) was also considered separately.As we consider the division of the infected human class into its 3 constituent sub-classes (asymptomatic, symptomatic and PKDL), we also assume that the weight of infection gained by the vector from each of these three infected subclasses (fHa, fHs, fHp) is unequal, assuming a higher weight of infection from the symptomatic human populations.With respect to the assumptions reflected, by formulating a set of coupled ordinary differential equations (Equations ( 1)-( 14)) that considers these compartments as variables, we attempted to understand the transmission dynamics of VL infection through these heterogeneous populations.

Results and Discussion
The set of ordinary differential equations (ODEs) were analyzed both analytically and numerically to interpret the biological significance of the VL infection.Primarily, positivity and boundedness conditions of the system of equations (Equations ( 1)-( 14)) were proven.Equilibrium points and epidemic threshold were calculated for the set of formulated ODEs.With respect to these calculated values the stability of the system was analyzed.
We observe that the right hand side of Equations (1) ( 14) are smooth functions of the variable , , , , these quantities are non-negati es, hence local existence and uniqueness properties hold in R .
It is also true for the other variables , , Therefore, the positivity condition holds.For boundedness of the solutions we propose the following theorem: Theorem 1.All the solutions of Equations ( 1)-( 14) which initiate in 14 0

R  are uniformly bounded if
; .
Proof: We define a function,

Deriving the Equilibrium Points from the Model
The model equations (Equations ( 1)-( 14)) were analyzed to calculate the equilibrium points for the disease free scenario and the endemic scenario.Each scenario was incorporated in the model as distinct mathematical conditions and the model was thoroughly analyzed.The equilibrium points in each situation were calculated using substitution method [13].

Disease Free Equilibrium
Disease-free equilibria were calculated for three distinct conditions-first, when there was no infection that means there were no infected sandflies; second, when animal populations were not infected and other two populations (human and sandfly) were infected and third, when animal and sandfly populations are infected and human populations were not infected.These conditions were created to substantiate the usability of the model in predicting disease-free equilibria for any single host and vector pair and eventually to compare the behavior in different conditions of specific host availability.
(a) Absence of infected sandflies: In order to create the situation during which there were no infected sandflies, for calculations of equilibrium point (SH      (1 ) *** , ( 1){ } )

A AA and A A A A
) -, where , and , are defined in the Equations ( 32)- (37).

Basic Reproduction Number
The basic reproduction number (R0) suggests the nature of the disease spread through a population [29].The basic reproduction number (R0) is an important parameter that gives the average number of susceptible individuals each infected individual would infect in a population over a period of infection.In theoretical epidemiology, R0 is calculated as a dominant eigenvalue r(K) of next-generation operator K [29].For our model, R0 can also be calculated as R0 = r (FV −1 ), where the matrices, for the new infection terms, F, and, of the transition terms, V, are given, respectively, by

Ht Ha H Ht
Hence for our system, the basic reproduction number (R0) can be calculated by substituting the expressions of the equilibrium points in the matrix F and subsequently finding the dominant eigenvalue of FV −1 through symbolic evaluation of the matrix inverse, multiplication and eigenvalue determination using MATLAB Mathworks R2012a, and performing algebraic manipulations to obtain a simplified expression as,

F H Hp Hp Ha H Hi Hp Hs
From our R0 expression it is clear that the explicit relationship between the parameters and the epidemiological threshold, basic reproduction number, is difficult to obtain.Since, R0 represents the average number of secondary cases generated when one infected case is introduced into a completely susceptible population, hence from the expression it can be observed that the number of sandfly infections generated by infectious human and animal (near the equilibria) is given by the sum of the product of the infection rate of infected humans and animals ( ) ,

Hi Ai H H A H s H p A A H H i A A i H H i
calculated at the equilibrium value [37], where the susceptible sand flies can acquire infection through blood meal from an infectious human (different groups) and animals.However, due to the complex expression of R0, numerical estimation of R0 and its variations due to different parametric conditions, which eventually can alter the disease endemicity, is presented subsequently.It is to be noted here that while studying the stability of the system around each of the equilibrium points (different disease-free and endemic scenarios) as well as variations of R0 through changes of parameters to observe endemicity, the respective equilibrium values are substituted in the R0 expression (or in the F matrix) and the results were obtained numerically.

Stability Analysis
The local asymptotic stability of the system around disease free and endemic equilibrium points were analyzed.Stability analysis was performed by finding the eigenvalues of the Jacobian matrix of this large system.Since from our model, the Jacobian matrix is of size 14 × 14, we have used the following definition [38,39] and Theorems 2 and 3 [40,41] related to the stability of matrices of higher orders.

 
T ; there exits a positive diagonal matrix such that J positive-definite (positive) stable matrices, Theorem 2. Let J be a 14 × 4 matrix of this model then J will be reducible to the form where blocks J1 and J3 are square matrices.Then J  ϑ, if J1  ϑ; J3  ϑ and J0 is a null matrix.
Theorem 3. Let J be a square real matrix.Then all the eigenvalues of a matrix J have negative (positive) real part if only if there exists a symmetric positive-definite matrix H such that, HJ + J T H is negative (positive) definite.Such a matrix J is said to be negative (positive) stable.

Stability Analysis around Disease-Free Equilibrium
The model stability was analyzed for disease free situations of each population group.Under disease free scenario, the possible situations for the model are: first when sandfly were not infected, secondly when animal populations were non-infected and last when human populations were non-infected.Therefore, using the matrices J1 and J3 found from Theorem 2 and Theorem 3 on ODEs (Equations ( 1)-( 14)), we analyzed the stability of the system (Equations ( 1)-( 14)).
(a) Absence of infected sandflies: The behavior of system could be understood by finding eigenvalues of J1 and J2 matrix.The eigenvalues were found after substituting the equilibrium points which were calculated for the sandfly uninfected case.As there was no infection in sandfly population, implies that there would be no infection through other population group.We obtained eigenvalues for this condition from the model as -µH, -µA, -µF, 0,0,0,0,0,0,0,0,0,0,0.When there is no disease transmission throughout the population, the system would be stable if Theorem 4. The system (Equations ( 1)-( 14)) is locally, asymptotically stable around the disease free equilibrium points (SH * , EH * , IHa * , IHs * , IHp * , TH * , RH * , SA * , EA * , IA * , RA * , SF * , EF * , IF * ) if < 1, otherwise unstable.
Proof: From the above eigenvalue analysis using Theorems 2 and 3, the proof of Theorem 4 is obvious if the inequality, < 1 holds true.

Stability Analysis around Endemic Equilibrium
To observe the behavior of the system when all three group populations were infected, the system was analyzed on the basis of eigenvalues of matrix J1 and J3.The eigenvalues were determined after substituting the critical points which were calculated for the endemic scenario.The eigenvalues calculated were, ( --), -, ( --), ( --), Ht .

Hi Ai F H Ha Ha Hp Hp Hs Hs
Hence we propose the following theorem.

Proof:
The above theorem can be proved easily by obtaining the eigenvalues of the jacobian matrix at the endemic equilibrium point by using the definition [37,38] and Theorems 2 and 3 [39,40] related to the stability of matrices of higher orders.The analytical expressions of some eigenvalues are not explicit as they contain the steady state values of the endemic equilibrium point, but if we use the positivity conditions for the endemic equilibrium point to exist (Lemma 1), then we can easily show that the inequality stated in Theorem 6 will hold true to obtain real and negative eigenvalues, which ensures the stability of the endemic equilibrium point.
The conditions obtained from local stability analysis and the system around each of the equilibrium points provides a strong biological basis to interpret the system behavior under different parametric conditions.The dependency of the biologically important parameters on other system parameters for obtaining the stable situation under different disease scenario can also be estimated from the above analysis.Moreover, the estimation of R0 and its variations due to different parametric conditions, which eventually can alter the disease endemicity, can also be tested.This further helps to get a coherent estimation of the control of disease conditions.All the above analytical results were further verified by numerical simulations.The existence of different equilibria (disease free or endemic) and the conditions stated in Theorems 4, 5 and 6 for the stability of the system around different equilibria are verified numerically using the biologically feasible parameter values in Table 1.The inequalities in Theorems 4, 5 and 6 are not satisfied and the existence conditions are violated if parameter values are not biologically feasible.

Numerical Analysis of the Model
The model was also analyzed numerically to observe the behavior of the infection transmission in the considered heterogeneous populations.Parameter values are given in Table 1.The numerical simulations were performed using the 4th order Runge-Kutta method implemented as a C program.The simulated results were plotted to visualize the behavior of model using MATLAB R2012a, Mathworks.In addition, the most sensitive parameters for infection to occur were identified from the model by performing local parameter sensitivity.The equations were incorporated as reaction rate laws in COPASI [42] for performing parameter sensitivity analysis.These sensitive parameters were varied in appropriate and wide ranges to capture the changes in the steady state behavior of some important model variables.

Dynamics of Visceral Leishmaniasis Transmission:
To compare the numerical simulations with the analytical solutions of the model, the previously mentioned disease free and disease-persistent scenarios were re-created and analyzed numerically.The model was simulated for 1000 days and dynamics of disease transmission was monitored.Simulations were performed using parameters, mostly obtained from literature and few are assumed to be biologically feasible, and are provided in Table 1.
The following situations in the model represent the different disease-persistent and disease-free scenarios (a) Vector populations remain uninfected: In order to create this situation, IF was kept to be zero.Hence, it is evident that absence of infected sandflies in the population would lead to no disease transmission.This could be captured from the model simulations as well (vector uninfected, Figure 2).The susceptible population stabilized around 1 and other populations stabilized around 0.
As no infection persists in the population, the equilibrium prevalence of all population fractions remains same as the initial value justifying our analytical observation.(b) No infections through the animal population: When there were no infected animals, the disease transmission was restricted to humans and sandflies alone (animal uninfected, Figure 2).In order to create this situations, the preference of sandfly blood meal to animal (pA) and human population (pH) was changed from their default values to 0 and 1.The infected fly population shows a transient increase within a few days as disease load is totally dependent on a single host population that is the human population, avoiding delay in infection spread as opposed to distribution within two hosts where the infection would be lost by clearance of infection from each host.Hence the corresponding infected human populations also show a transient increase.
As the susceptible human population (SH) becomes continually infected, SH shows a rapid decrease until it is totally transported.In addition, with respect to immunity of individual the recovered class shows a similar increase.(c) No infections through the human population: When there were no infected humans, the disease transmission was restricted to animals and sandflies alone (human uninfected, Figure 2).Similar to the non-infected animal scenario, the preference of sandfly blood meal to animal (pA) and human population (pH) was changed from their default values to 1 and 0 respectively to maintain infection only through animals.The infected fly population shows an increase similar to the animal uninfected case for the initial few days after which a rapid decline is observed.However, comparing both the situations, it can be observed that steady state of IF is around 100 fold higher in the situation when only animals are present as compared to situation when they are absent.This suggests that the animal population act as the major contributor for continual infection.(d) All populations were infected: Under the influence of both the host populations (human and animal), the dynamics of infected sandfly population demonstrated a combinatorial effect where the initial population rise resembled the infection under presence of only the human infected population and the latter persistent infected population resembled the infection under presence of only the animal infected population (all infected, Figure 2).Another observation was that the steady state of recovered animal population decreased further in presence of both animals and human populations as compared to situation where only animals were present suggesting an increased flow of infection from animals to humans via the increase in proportion of infected sandflies.A similar behavior of each population considered in our model was also observed in a previous study [17].

Figure 2.
Numerical simulations of the model were performed in the vector uninfected, animal uninfected, human uninfected and all infected scenarios for every variable (in the human, animal and sandfly compartments) considered in the model and the temporal behavior of the system was monitored for each population group.
An interesting observation in this context (i.e., for endemic steady state) is that for a wide range of parameter variations (10 5 fold increase and decrease from the basal values mentioned in Table 1, for all the parameters involved in the inequality stated in Theorem 6) the condition of Theorem 6 holds true, i.e., the system is stable around the endemic steady state.Due to the largeness of the system, we have not been able to prove explicitly the global stability of the system.However, from our numerical simulation and biologically feasible parameter variations within a wide range of values we can expect globally stability around the endemic equilibrium point.

Animal Populations Act as a Reservoir of VL Infection
The equilibrium prevalence of infections in each population under the above created scenarios was compared (Figure 3).It can be observed that the equilibrium value of infected sandfly population (IF) was high under the human non-infected scenario (Figure 3A) as compared to the animal non-infected scenario (Figure 3B).Similarly, while considering spread of infection in all the three populations (Figure 3C), the equilibrium value of animal infected population (IA) was high as compared to every infected human compartment as well as the total infection in human.Both these observations suggest that the animal population contributes more to the VL infection than the human, clearly indicating that sandfly prefers animal over human.Further, R0 values were calculated in each of the aforementioned scenarios (Figure 3D).R0 for the non-infected animal scenario was higher than non-infected human scenario confirming that infection spread is highly modulated by the animal population.It has been previously reported that when a reservoir host is removed from the system, R0 increases [43].This confirms that indeed the animal population acts as a reservoir for VL infection.

Designing Disease Control Strategies from Model Predictions
The most important aim of our study is to model the complex disease dynamics of visceral leishmaniasis and to understand important parameters that can control disease transmission among different populations.By identifying the most important parameters that govern VL infection, it can be possible to develop optimal strategies that can lead to the control of these parameters and hence, cause reduction in spread of the infection through populations.For identification of the most important parameters, a two-way analysis was performed where the most important parameters that are sensitive to infected population groups in our model was identified, followed by variation of these sensitive parameters in feasible ranges to understand the dependency of equilibrium infected populations on these parameters.

Identification of Important Parameters that Control VL Infection
Local sensitivity analysis is a mathematical procedure that identifies the most sensitive parameters in a system of ODEs for a given output [44].Local sensitivity analysis was performed in our model using COPASI 4.11 (Build 65) [42] to identify the most sensitive parameters that influence VL disease dynamics.The importance of local sensitivity analysis is that it suggests the importance of parameter that is affecting the output of the model when perturbed in a very small local range.According to the standard procedure, every parameter was increased or decreased up to 0.1% of its original value (keeping the other parameters constant) and then its sensitivity in this local range to steady state of the output as checked.The sensitivity of variables is measured by the values of the sensitivity coefficients.The scaled sensitivities of each of the model parameters to steady state of infected human, animal and fly populations (IHa, IHs, IHp, IA, and IF) are plotted in Figure 4.The negative values of sensitivity coefficients indicate negative effect of that particular parameter on steady state and positive values of sensitivity coefficients indicate that the particular parameter has a positive effect on the steady state of all infected population groups considered in the model.As some parameters are common for different events occurring in the model, the variable affected during the event has been indicated in brackets along with the parameter.For example, μF, death rate of a sandfly is a parameter that is common for susceptible, exposed and infected flies.However, µF (SF) indicates the death rate of susceptible flies.
There are few parameters αH, αF, αA, μH, μA and μF, which are common for different populations in the model, the population affected in the reaction has been indicated in brackets along with the parameter.From Figure 4, it can be observed that the biting rate β as expected has a positive effect on the infected individuals of each population.Similarly, the birth rates of each population (αH, αF, and αA) have a positive effect on the corresponding variables.In addition, the mortality rates (µH, µA) of each population except IF have very small negative effect on their corresponding variable.In Figure 4A, we could observe that ρHs has negative effect on IHa suggesting that increase in duration between asymptomatic infection and symptoms lead to reduction in asymptomatic infection.Another trivial result was that the increase in recovery rate of infection (γHs, γHp, γA) reduces rate of infections (Figure 4A-D), which depends on immunogenic potential of the corresponding population to the parasite infection.Other trivial results include the highly negative effect of µF (IF) and on IF (Figure 4E).In addition, the death rate of susceptible fly µF (SF) has negative effect on IF indicating cessation in transmission of the infection from susceptible to infected under increased death rate of susceptible flies.
A unique result was that ρHs (recovery rate of symptomatic individuals) had a positive effect on IHp probably suggesting the flow of infection through an alternative route, i.e., IHp (Figure 4C).Conversely, γHs demonstrates a negative effect on IHp (Figure 4C) signifying the reduction in the rate of conversion of IHs into IHp.Further, it was further observed that the recovery rate of animal from infection (γA) had a high negative effect on IHp and IF (Figure 4E) further suggesting the role of the animal population in absorbing the infection towards itself.

Variation of Parameters to Control VL Infection Spread
The parameters identified through sensitivity analysis were varied in appropriate ranges to study equilibrium dependence of the infected populations on each identified parameter.By performing extensive parameter variations, we also compare the choice of parameter values and identify ranges where disease transmission would be considerably reduced.Further, in each case, change in R0 was also monitored with corresponding change in the parameter value to identify epidemic thresholds for each parameter that can lead to termination of the disease.
(a) Infected sandfly biting rate vs. fly mortality rate: As indicated by parameter sensitivity (Figure 4), the infected sandfly biting rate (β) has a positive effect whereas the fly death rate has a negative effect on equilibrium infected populations.Varying β between 0-1/person/day, it can be observed that all the equilibrium infected populations demonstrate a sigmoidal increase with variation in biting rate (Figure 5A).When β < 0.2/person/day, there are no or low proportion of infected individuals in each population.But for a β > 0.2/person/day, the proportion of infected individuals switches from low numbers to considerably high numbers, representing the spontaneous spread of infection through all the populations.Change in infection spread can also be monitored by R0.Hence, varying the infected sandfly biting rate, it can be observed that R0 < 1 when, β ≤ 0.168 suggesting no infections in the population.Similarly, varying the fly death rate (μF) between 0-0.5/cycle, a sigmoidal decrease in proportion of all infected humans can be observed with increase in fly death rate (IHa, IHs, IHp in Figure 5B).At the same time, an exponential decrease can also be observed in the infected animal and fly populations (IA, IF in Figure 5B).In addition, R0 exponentially decreases with increase in μF.An R0 < 1 can be observed when μF > 0.1.Thus, maintaining the infected fly biting rate β ≤ 0.168/person/day or mortality rate of flies μF > 0.1/cycle, can lead to cessation of the disease.
Apart from varying either of the parameters, we also perform variations in the above two parameters simultaneously to identify combinatorial ranges of parameters where the infection dies out.It can be observed that for mortality rate of flies (μF) > 0.23, R0 is always less than 1 (Figure 5C).This suggests that μF is more influential as compared to β in reducing the infection.However, when μF < 0.12, β contributes more in spreading the infection among different populations.This can be observed by the high R0 values (orange colored regions, Figure 5C) under the presence of increasing β.This indicates that combinatorial strategies that can control both the infected sandfly biting rate (β) and mortality rate of flies (μF) in the predicted ranges can lead to better reduction of disease spread.(b) Recovery rate of animal vs. symptomatic humans: As indicated by parameter sensitivity (Figure 4), both the recovery rate of animal and recovery rate of symptomatic humans negatively affect the infected populations.Varying recovery rate of animal (γA) between 0-0.06/day, it can be observed that the fraction of individuals in each infected class of humans demonstrate a delayed sigmoidal decrease (IHa, IHs, IHp in Figure 6A).Whereas, infected animals and flies demonstrate a rapid exponential decay with increase in γA (IA, IF in Figure 6A).Varying γA further leads to an exponential decrease in R0.However, in the given parameter ranges of γA, the infection did not die out i.e., R0 was never less than 1.This suggested that the animal recovery rate does reduce the occurrence of VL disease among populations but does not lead to cessation of disease.Further, the infection that persists in the population is because of the infected human populations acting as a sole source of infection.Similarly, varying the recovery rate of symptomatic individuals (γHs), it can be observed that there is an exponential decrease in relative proportions of all the infected populations considered in the model (Figure 6B).In addition, it can be observed that the rate of decrease in the steady state asymptomatic and PKDL-infected populations is slower than the rate of decrease in symptomatic populations.A similar behavior can be observed for the transient class (TH) also.This suggests the role of IHa, IHp and TH in acting as important sources for VL infection for a certain period of time, even after increasing the recovery rate of symptomatic individuals.Further, R0 also demonstrates an exponential decrease with increase in γHs.However, in the given parameter ranges of γHs, the infection did not die out i.e., R0 was never less than 1.As increasing the recovery rate of animals and humans individually did not lead to cessation of disease spread, it is necessary to vary both the parameters and observe its effect on persistence of infection.
Apart from varying either of the parameters, we also perform variations in the above two parameters simultaneously to identify combinatorial ranges of parameters where the infection dies out (Figure 6C).We monitor R0 while varying the recovery rates γA and γHS.It can be observed that for γA > 0.05 and γHS > 0.04, R0 is less than 1.Thus, for optimal disease control, it is required to develop strategies that can maintain high rates of recovery within both animal and host.This also indicates that healthy and rapidly recovering animals contribute towards reduction of VL infection by absorbing the infection towards itself and thereby reducing the load towards the recovering human population.

Testing the Model with Real Data
The monsoon season with moderate temperature and relative humidity, is the most conducive environment for the breeding of sandflies.A significant positive correlation between sandfly density and rainfall is known for the sandfly species Phlebotomus papatasi [45].Further, sandfly abundance is known to the highest in monsoon and post-monsoon seasons [46].With respect to this information, it is palpable that seasonal fluctuations like rainfall affect the existence of effective sandfly population.To reinforce our model predictions with real data, we intended to include this seasonal effect on sandfly density.
Thus, in order to test the predictive capability of the model with real population data, we replaced the birth rate of sandfly as a binormal distribution function F(t) in the model, which was fitted for the bimodal nature of the rainfall distribution in the Muzaffarpur district, Bihar, India for the year 2005 [24].
Here, t denotes the time of the year (in days); a1, a2 gives the peak value of mosquito density in that year; c1 and c2 gives the width of the two normal distributions; and b1, b2 denotes the position of the peaks.As a case study, the values of a1, a2, b1, b2, c1, and c2 were adjusted so as to appropriately predict sandfly density from reported rainfall data [26].This function was then introduced as sandfly birth-rate in the model.To perform this simulation, the default parameter values were retained.For the purpose of validation, the behavior of this modified model was scaled and fitted with the number of VL infected symptomatic human cases reported for Muzzafarpur district, Bihar, India in the year 2005 [26].While performing simulations, the initial values of the human compartments were varied in proportion so as to relatively match the predicted symptomatic cases with the monthly average of January 2005.Using these initial values for the human compartments, the model was simulated for the whole year.The average number of predicted symptomatic VL cases for each month were then calculated from the model simulations, scaled with respect to total number of observed VL infected cases in that year and compared with the average number of known cases in each month (for the year 2005).A recent study has reported the evidence of zoonotic transmission of visceral leishmaniasis in the state of Bihar, India [47].Hence, to create the actual scenario, we retained the effect of animal populations in the transmission of the disease while performing the model simulations.Data in Figure 7B was the average number of symptomatically infected humans recorded in that region for each month in the year 2005.The goodness of fit of model prediction with real observed data was reasonable (R 2 = 0.69) substantiating the strength of model to predict the dynamics of VL disease transmission.Since most of the environmental, epidemiological and entomological parameters for this specific region were not available, hence we restricted ourselves to predict for only one year using the global parameter values for our general model.It is worthy to note that even if the model is general in nature it is still able to capture the trend of the disease for a specific region and the predictions are also reasonable.This confers the model with wide applicability to different epidemiological situations and ability to predict disease incidence if all the realistic parameter values are available along with the information of rainfall for a specific region.
Further, to exhibit the strength of the model (in which sandfly density was fitted for rainfall distribution) to predict control scenarios, the sensitive parameter values (β, μF, γA and γHs) identified in the previous analysis were fixed to regimes where R0 was predicted to be less than 1 in the previous analysis, and the simulations were performed for 12 months (Figure 7C).The population proportions were then scaled according to the VL infected cases in the year 2005 to understand the control effects on real population.It can be observed that the given combination of β and μF parameter values imposed the greatest degree of control in VL infected cases (two-fold reduction in the peak number of cases as compared to the default situation) as compared to their singular controls.Similarly, the given combination of γA and γHs parameters reduced the VL infected cases to relatively low values as compared to their singular controls (around 1.7-fold reduction in the peak number of cases as compared to the default situation).This extended analysis further demonstrates the usage of the model in predicting the control measures for yearly VL infection.

Conclusions
In this study, we propose a compartment-based mathematical model that explains the visceral leishmaniasis transmission process in three distinct populations-the human, sandfly and the animal considering the detailed division of the human infected population into the symptomatic, asymptomatic and PKDL-infected classes.Our model is unique as compared to other models of visceral leishmaniasis as it integrates the detailed classification of the human infected population along with other heterogeneous populations (vector and animal) to study their combined effect on disease transmission.Further, the proposed model is generic for VL disease transmission and can be applied to any region/conditions by introduction of region-specific parameter values and conditions.However, the model of this scale posed a challenging task for analytical study and estimation of different epidemic thresholds, including basic reproduction number.We could analyze the model both analytically and numerically in different diseased and disease-free scenarios, to understand the equilibrium prevalence of infection in different populations of the model under various situations.In addition, performing parameter sensitivity, biting rate of the sandfly, birth rate of sandflies, recovery rates of symptomatic humans and animals were identified to be the most sensitive parameters.One of the interesting outcomes of this analysis was to establish the fact that rapidly recovering animal populations may act as absorbers of VL infection.In addition, through our model parameter variation analyses, we demonstrate the potency of the model in enumerating equilibrium thresholds for different model parameters and hence, suggest solutions for curbing disease spread among populations.The predictions from the model suggests that measures should be taken to either control infected sandfly populations by eradicating their populace through external perturbations, like insecticide sprays or by avoiding contact of human or animal populations with the infected sandflies.Further, the drain of the infected human or animal populations into the recovered state can be accelerated by discovery of vaccines or drugs that can increase the rate of their transition from infected to recovered class.Through our parameter variation analysis, the role of asymptomatic, PKDL and transiently infected classes as constant sources for infection has been established and further suggest the requirement of a mass-identification programme to identify the total number of individuals belonging to these classes and thereby design control strategies to target each of these individual classes.Through our analysis, we prominently stress on the control of a combination of parameters (in ranges where R0 is less than 1), as the most optimal method to control VL disease spread.As an important measure to effectively reduce the VL disease spread, a combinatorial control strategy that can reduce biting rate of sandflies and increase mortality rate of sandflies needs to be implemented.In contrast to zoonotic disease control for visceral leishmaniasis, like culling of animals [48,49], we observed that VL spread among human population decreases in the presence of animal population.Thus, a combinatorial strategy to target both animal and human recovery rates can highly contribute to the disease elimination programme.Model predictions were also comparable to real infected population data demonstrating the applicability of our model to predict the region-specific nature of the VL infection.This also emphasizes the importance of our model in predicting the expected number of asymptomatic, PKDL, transient and recovery classes of human hosts while comparing the predicted number of asymptomatic cases with actual reported VL symptomatically infected cases.Through our generic model analysis and predictions, we conjecture that our generic model can be used to understand the disease dynamics and spread as well as may be useful to provide practical solutions to effectively control visceral leishmaniasis in the human population.

Figure 1 .
Figure 1.Schematic representing the flow of infection between different population groups considered in the model.Each directed arrow represents the flow of infection.The corresponding parameters are given above each arrow.Arrows leaving a compartment but not connected to any other compartment represent the mortality (death) of individuals in that compartment.The arrows directed towards the susceptible compartments and not connected with any other compartment represent the birth rates of that population.

Figure 3 .
Figure 3.Comparison of equilibrium and R0 values in different scenarios created from the model.(A) Comparison of equilibria of infected populations in the non-infected animal scenario; (B) Comparison of equilibria in the non-infected human scenario; (C) Comparison of equilibria in the default, all infected scenario; (D) R0 values in the above three scenarios.

Figure 4 .
Figure 4. Local parameter sensitivity analysis for different model parameters and its effects on infected population groups.(A) Local parameter sensitivity for the variable IHa; (B) Local parameter sensitivity for the variable IHs; (C) Local parameter sensitivity for the variable IHp; (D) Local parameter sensitivity for the variable IA; (E) Local parameter sensitivity for the variable IF.As some parameters are common for different events occurring in the model, the variable affected during the event has been indicated in brackets along with the parameter.For example, μF, death rate of a sandfly is a parameter that is common for susceptible, exposed and infected flies.However, µF (SF) indicates the death rate of susceptible flies.

Figure 5 .
Figure 5. Parameter variations in β and μF-(A) Effect of β on equilibrium infection dynamics and R0; (B) Effect of μF on equilibrium infection dynamics and R0; (C) Effect of simultaneous variations in β and μF on R0.

Figure 6 .
Figure 6.Parameter variations in γA and γHs-(A) Effect of γA on equilibrium infection dynamics and R0; (B) Effect of γH on equilibrium infection dynamics and R0; (C) Effect of simultaneous variations γA and γHs on R0.

Figure 7 .
Figure 7. Testing model with real data-(A) Monthly average sandfly density versus monthly average rainfall in 2005; (B) Observed (actual) VL-infected cases in the Muzzafarpur district, Bihar, India (2005) vs. predicted VL infected cases from the model.R 2 indicates the goodness of fit of the predicted cases to the actual population data; (C) Control situations for previously considered parameters in the model (modified to accommodate the sandfly density function) and their effects on VL-infected populations.

Table 1 .
As mentioned before, the parameter values have been scaled to per individual in a population so as to compute dynamic changes in relative fraction of individuals in each compartment.For simplicity of the model, we have also assumed that αH = μH, αF = μF and αA = μA.

Table 1 .
Model parameters and their numeric values used for further simulations.
Therefore total animal population would be equal to susceptible animal population.The equilibriums of all other variables were found by substituting value of IF**.Hence, Absence of infected animal population: In this situation, only disease transmission between human and sandfly populations was considered but animal populations remained non-infected.To calculate equilibrium point (SH**, EH**, IHa**, IHs**, IHp**, TH**, RH**, SA**, EA**, IA**, RA**, SF**, EF**, IF**) in this scenario, EA**, IA** and RA** were set to zero.When animal populations were not infected, the disease transmission was between human and sandfly population alone. -) -) ) )