A Comparative Assessment of Visceral Leishmaniasis Burden in Two Eco-epidemiologically Different Countries, India and Sudan

The two hyper–endemic regions for Visceral Leishmaniasis (VL) in the world are located in India and Sudan. These two countries account for more than half of the world’s VL burden. The regional risk factors associated with VL vary drastically per region. A mathematical model of VL transmission dynamics is introduced and parametrized to quantify risk of VL infection in India and Sudan via a careful analysis of VL prevalence level and the control reproductive number, , a metric often used to characterize the degree of endemicity. Parameters, associated with VL-epidemiology for India and Sudan, are estimated using data from health departmental reports, clinical trials, field studies, and surveys in order to assess potential differences between the hyper–endemic regions of India and Sudan. The estimated value of reproduction number for India is found to be 60% higher than that of Sudan ( and ). It is observed that the is most sensitive to the average biting rate and vector-human transmission rates irrespective of regional differences. The treatment rate is found to be the most sensitive parameter to VL prevalence in humans for both India and Sudan. Although the unexplained higher incidence of VL in India needs to be carefully monitored during long-term empirical follow-up, the risk factors associated with vectors are identified as more critical to dynamics of VL than factors related to humans through this modeling study. Author Summary The Visceral Leishmaniasis (VL) is a neglected tropical disease, primarily endemic in five countries, with India and Sudan having the highest burden. The risk factors associated with VL are either unknown in some regions or vary drastically among empirical studies. In this study, we collect VL-related data from multiple sources for the two different countries, India and Sudan, and use techniques from mathematical modeling to understand factors that may be critical in the spread and control of VL. The results suggest that the risk factors associated with disease progression are important in explaining high VL prevalence in both the countries. However, the likelihood of disease outbreak in India is much higher than that in Sudan and the probability of transmission between human and sandfly populations vary significantly between the two. The results have implications towards VL elimination and may require a review of current control priorities.


Introduction
1.  The Leishmania donovani transmission cycle is anthroponotic and takes place from human to human The model system is given as:  Disease-induced mortality is not included because, due to institutionalized treatment, deaths from VL are negligible. For simplicity, the human population is assumed to be constant. Λ h denotes the recruitment rate into the susceptible population, and µ h denotes the per-capita death rate. Because N h approaches Λ h µ h when t approaches ∞, we assume, without loss of generality, that N h = Λ h µ h [16]. A susceptible individual acquires the L. Donovani parasite following an effective contact with an infectious sandfly. The rate λ vh , the force of infection on humans, is given by where the right-hand expression (Equation 3) is given by the product of the per-vector daily biting rate of sandflies (b), the VL infection transmission probability, given a bite from an infected sandfly to human (β vh ), the average number of sandflies per humans m v:h , and the proportion of infectious sandflies in the vector population ( Iv /Nv). It is assumed that all newly VL-infected humans go through an asymptomatic (symptomless) stage (A h ). After an asymptomatic period of several months, humans develop clinical symptoms at the per capita rate φ h , moving to the infectious class I h . During the infectious period, humans will seek VL treatment at the per capita rate θ h , proper treatment leads to recovery at the per capita rate γ h (recovered individuals gain lifelong immunity). Newly emerging adult female sandflies are recruited into the susceptible population at rate Λ v and die at the per-capita rate µ v . The sandfly population is assumed constant. A susceptible sandfly is infected following an effective contact with infectious humans at the per capita rate λ hv (force of infection on sandflies). The rate λ hv is given by where the right-hand side is the product of: the per vector daily biting rate (b); the probability that 139 susceptible sandflies acquire the Leishmania parasite while feeding on a VL-infected individuals (β hv ); 140 the proportion of VL infectious humans in the human population ( I h /N h ). It is also assumed that the Interactions between vector biting behavior and uneven pathogen transmission potential between hosts 146 may lead to difficulty in controlling infection. How vector species respond to availability of hosts is 147 highly variable and has fostered considerable interest among vector borne disease modelers for decades.

148
The proportion of blood-meals taken by vectors from the host species of interest is generally assumed to 149 increase directly with increasing human availability and changing levels of vector density. Hence, vector 150 biting can play a significant role in the transmission process [88]. The biting rate of sandflies is typically 151 a function of ambient air temperatures, humidity, wind speed, vector density and local habitat. There 152 remains a couple of challenges in effectively using biting rates, namely, which is a proper functional 153 response to capture biting rates in the model and how to measure it precisely from the field data.

154
To effectively use models to make reasonable definitions, models must be carefully parameterized and 155 validated with epidemiological and entomological data. On the other hand, researchers have modeled 156 biting rate in different ways but realistically the biting rate may vary according to the abundance of 157 hosts and to vector preference [87]. In this study, we suggest alternative forms of transmission terms 158 as well as use distinct data sets to estimate parameters of the two different terms (vector-to-host and 159 host-to-vector terms).

161
This section first provides a careful derivation of incidence rates expression as a function of landing and 162 biting rates and then use landing rate data to estimate the transmission probabilities from sandflies to 163 humans (β vh ) and humans to sandflies (β hv ).

164
The human incidence rate (Equation (3)) is a function of the average rate of interactions between vectors and humans, which in turn is directly proportional to the proportion of infectious sandflies Let b denote the average number of bites per sandfly per unit time and ρ the average number of bites received per human per unit time. Assuming that all sandfly bites are to humans only, we must have that the total number of bites made by all sandflies per unit of time (bN * v ) equals the total number of bites received by all human hosts per unit of time (ρN * h ). Thus, we have that The assumption that ρ is constant is customary in the literature although there are some studies where the host vector ratio is assumed not constant over time [85]. We further assume that the average number of bites received by a human per unit time is proportional to the number of sandflies landing on an individual per unit time, that is, ρ ∝ .
Hence, the total number of effective landings on all humans from all sandflies per unit time is where the first (second) term of (7) accounts for the total number of effective landings on all humans from 165 all susceptible (infected) sandflies per unit time. That is, the total effective landing/feeding of vectors on 166 humans is a function of the total vector population, which includes both susceptible and infected vectors.

167
In epidemiology, of importance are only the two cases when landing occurs from a susceptible sandfly If β vh is the per-person transmission efficiency (that is, probability that infection is successfully transmitted from vector to human given an infected bite), then the rate at which VL is transmitted to humans is using Equations (5) and (6). 174 Similarly, we can derive the infection rate in the vector population generated by infected humans. If¯ accounts for the average number of times a sandfly lands on humans per unit time, then the total number of effective landings by all sandflies on all humans is It should be noted that the total N h =¯ N v and that, while accounting for new incidences in sandflies, we are interested in landings occurring from susceptible sandflies on infected humans only. Hence, the term¯ is the one that plays a role in accounting for new sandflies incidences, while the remaining terms aren't. If we let β hv be the per-person transmission efficiency from human to vector (i.e., transmission probability per bite on infectious humans that leads to infection in a susceptible sandfly), then the total number of sandflies who acquire infection while effectively landing on infected humans per unit time is using Equations (5) and (6) as the number of secondary infections caused by a single infective introduced in a primarily susceptible 187 population (i.e., N ≈ S 0 ) but in the presence of interventions [12,20,53,84].

188
Using our model, we compute R C using the next generation operator approach [12,13,84], a process that requires the computation of the matrix of new infection terms, F, and the matrix of transition between compartments, V. The R C is the spectral radius of the next generation matrix, ρ FV −1 (see section B.1 for derivation), in the presence of treatment program (where the treatment rate is θ h ) and is given by where is the landing rate on a human, b is the biting rate per sandfly, and β vh the number of infections in humans generated by one infected vector. The expression is the average number of new cases vectors generated by one infected human and bβ hv µv represent average number of new cases in humans produced by one infected vector. Hence, R C (θ h ) is given by the geometric mean of two sub "reproduction" numbers The nocturnal activities of various sandfly species start at around 6:00 pm -9:00 pm, peaks between 239 the hours of 11:00 pm -1:00 am and ends between the hours of 3:00 pm and 6:00 am. A rapid rise 240 to a maximum pick and then a sharp decline observed in data from various field studies suggest the 241 probability distribution for biting and landing rate would best be fitted with a triangular distribution.  . Collected data of P. Argentines was first averaged over seasons and then fitted to the triangular distribution to estimate parameters of the distribution representing landing rate. The mean and 95% Confidence Interval for landing rate distribution are given in Table 2.  Figure 4. Collected data of P. Orientalis was first averaged over months and then fitted to the triangular distribution to estimate parameters of the distribution representing landing rate. The mean and 95% Confidence Interval for landing rate distribution are given in Table 2. for which we can obtain data with relatively high certainty and mathematical methods to estimate the parameters representing transmission probabilities required in our model. Two novel approaches, that uses endemic prevalence from the model, were designed to estimate the transmission probabilities as described in the next two sub-sections. Note, the unique endemic equilibrium of the model is stable (as shown in the Section 3) and is given by Since VL is endemic in both India and Sudan, we use these expression to obtain prevalences and thereby 250 use them to estimate transmission probabilities (i.e. β vh , β vh ). We assume and hence, the host and vector populations becomes constant. The prevalences in humans and sandflies 252 populations are given by The equations (16a)  and (15b) for the prevalences can be rewritten in terms of R C as follows: Isolating β vh and β hv from (17a) and (17b) we obtain, The estimated distributions using the this approach are given in Figure 5.   of parameters on which R C depends) results in N samples for R C in each realization.

290
After 10 realizations, the mean (µ) of point estimates of R C , standard error (σ) of R C , and the 291 probabilities that R C estimates fall below and above the threshold value one for India and Sudan were 292 collected ( Statistical analysis on the differences in the means of R C for India and Sudan was carried out using  that R C for India is indeed significantly higher than that the one for Sudan using model generated R C , 305 re-scaled by population size, we proceed to determine what are the parameters that if modified generates 306 the larger change in R C .

307
The PRCCs were calculated for each country in order to quantify sensitivity of model parameters on 308 the R C estimates. We observe the sign and the magnitude of the PRCC values for each parameter above 309 the line y = ±0.3 for each respective country.

On endemic prevalences (P
Parameter uncertainty and sensitivity analyses are also performed on the Prevalence of the infected 312 populations (P A h , P I h , and P Iv ). These analysis are used to assess which of the same eight input are most significant to estimating endemic prevalences.

314
As described in the Section 4.2.1 on R C , the similar sensitivity and uncertainty analysis procedure 315 was carried out on the endemic prevalences for both the countries. However, higher number of samples 316 (50, 000) for each parameter were obtained. The first 10, 000 sample-sets (out of the 50, 000) that resulted 317 in R C > 1 (condition for existence of the endemic prevalence) were eventually used in the analysis. This 318 is because that endemic equilibrium only exists and stable when R C > 1 . The estimated distribution of R C I from uncertainty analysis, is shown in Figure 7a. The mean estimate 322 of R C I for India is found to be 2.05 with a standard deviation of 1.09. The sensitivity analysis of R C I 323 provides the ranking of parameters based on their influence on R C I (Figure 7e). In decreasing order of 324 influence, the parameter ranking was θ h , being the most sensitive parameter, followed by b, , β vh , β hv , 325 and the least sensitive parameters are φ h followed by µ v .

327
The estimated distributions of prevalence (P A h , P I h , and P Iv ) are shown in Figure 7b-7d. The mean 328 estimate of P A h was found to be 0.0045 with a standard deviation of 0.0019. The parameter φ h was 329 found to be the most influential parameter on the prevalence of asymptomatic, P A h . The remaining 330 parameters in descending order of magnitude of PRCC were, θ h , , β vh , and β hv , with µ v and µ h being 331 least sensitive parameters to P A h . The sensitivity analysis performed on P I h reveal that the treatment 332 rate, θ h is the most influential parameter for changing disease prevalence. The mean estimates of vector 333 prevalence were found to be 0.0526 with a with a standard deviation of 0.0432. From our sensitivity 334 analysis of P Iv we observe in Table 7 and Figure 7h, that there are four most influential parameters.

346
For the asymptomatic prevalence, P A h , we estimated a mean prevalence of 0.0024 with a standard 347 deviation of 0.0018. Results of uncertainty analysis is shown for Sudan in Figure 9b. From Table 8 and 348 Figure 9f, we observe that the prevalence of asymptomatic population is negatively correlated but most 349 sensitive to φ h , followed by the parameters β hv , θ h , b, β vh , and . The natural death rates, µ v , and  Table 8 and displayed in Figure 9g shows that the treatment rate of infectious  Figure 9. Results for Sudan: Uncertainty of the Reproduction Number (Subfigure 9a) and the Prevalence (Subfigures 9b -9d) of Asymtomatics, Infectious Humans and Infectious Sandflies, respectively. Tornado plot showing partial rank correlation coefficients (PRCCs) of the Reproduction Number (Subfigure 9e) and the Prevalence (Subfigures 9f -9h) of Asymtomatics, Infectious Humans and Infectious Sandflies, respectively. deviation of 0.0088. Our analysis identified the parameters sensitivity to changes in P Iv (Table 8 and 359 Figure 9h). The result shows that the treatment rate, θ h is the most dominant parameter followed by b, 360 β vh , and µ v . The less influential parameters on P Iv are µ h , , φ h , and β vh .

362
Parameter estimates were obtained either from the literature or estimated from field data, and were used 363 for an evaluation of country-specific risks. The risk was quantified to study differences and similarities in 364 VL disease burden in India and Sudan. In this section, we conduct comparative (between two countries) 365 assessment by studying impact of change in parameter estimations on VL disease burden in these two 366 countries when risk is measured either in terms of R C or prevalence of infection. The assessment was 367 based on uncertainty and sensitivity analyses..  Table 6 for model parameters (h)  Sudan (t-test with H 0 : µ(R C S ) = µ(R C I ) against H 1 : µ(R C S ) = µ(R C I ) where the µ represents mean 373 of R C I and R C S ). The analysis suggested rejection of null hypothesis (Table 9), that is, the obtained 374 point estimates of R C between India and Sudan are different. We also performed Kolmogorov-Smirnov 375 test between empirical distributions of R C for the two countries and found that empirical distributions 376 are not the same.  Figure 11. (a) The comparison between estimated distributions of R C for India and Sudan. The box plot compares the mean(•), median, minimum, and maximum of R C estimates for both countries. It is found that the gamma, is a best-fitted distribution for the samples from the uncertainty analysis. Table 2 Figure 14a -14b). Sensitivity analysis shows that P I h is most sensitive to θ h , , b, β vh , and β hv and 402 least sensitive to µ h , µ v and φ h for both countries (Table 12 and Figure 14c). The treatment rate, the 403 first most sensitive parameter, and µ v , µ h and φ h in same decreasing order of influence, are common 404 parameters for both countries. For India, parameters ranking in descending order is , β vh , β hv and b 405 whereas for Sudan the order of parameters is β vh , b , β hv , and .  were most sensitive to the prevalence of infection in sand flies, P Iv , for both countries (see Table 13 and 411 Figure 15c). Between the two parameters, β hv (b) is relatively more sensitive for India (Sudan). The  Table 13. A comparison of the partial rank correlation coefficients for input parameters of the output value (P Iv ). Where (*) denotes p < 0.01. for India and Sudan.   Table 3. Summary of estimates of the transmission probabilities, β hv and β vh , using the two approaches with mean and ranges for other parameters (Table 2) for India and Sudan were fixed.   Table 9. Statistical estimates of quantities, R C , P A h , P I h , and P Iv , for VL in Sudan and India using the 2 sample t-test and two-sample Kolmogorov-Smirnov test. All analysis were found to be significant, i.e. p < 0.05. Table 10. A comparison of the partial rank correlation coefficients for input parameters of the output value (R C ), where (*) denotes p < 0.01. for India and Sudan.

415
The regional risk factors associated with VL are complex and ambiguous. In face of this uncertainty 416 systematic evaluation of ongoing VL control programs is essential but it remains challenging, as appro- proof of infection, hence in this study we used multiple data sets to obtained ranges of the parameters.

429
This is the first study to best of our knowledge that review and make use of extensive collection of 430 available data on epidemiological and ecological parameters to understand the dynamics of the Visceral 431 Leishmaniasis (VL) and identify risk factors in India and Sudan using mathematical modeling approach.

432
The study compares and contrasts quantities from two nations where the disease is endemic and spread analysis on the R C also showed that there were eight parameters (see Table 2) that should be taken vectors. This campaign reduced the number of VL cases during this period (1962-1963,) showing no new 452 prorated cases. It was observed that soon after the DTT spraying campaign stopped the number of VL 453 cases were elevated to higher epidemic levels [5]. High treatment rate is also found to be a critical factor 454 in impacting the dynamics of VL but primarily in India. However, we assumed effective treatment for all periods and therefore may not be a representative of the country. However, the study clearly identifies 464 the type of data that are relevant and needs to be collected for thoroughly understanding of VL dynamics.

465
In our future research, we plan to provide elaborate analytical methods for the estimation of partially 466 observed data (usually temporal incidence data) for the two developing countries.   but not yet infectious individuals move into the asymptomatic population (sub-clinical infection, exposed 757 to VL but not yet infectious), who may exit the system through natural death or through progress to 758 clinical VL. The change in A h population is The asymptomatic can then progress to a VL clinical symptoms stage (I h ) at the rate φ h : where θ h is the per-capita treatment rate and µ h is the per-capita departure rate.
The population of recovered individuals from VL (R h ) is increased following successful treatment, leading 764 to permanent immunity into the R h class (at the rate γ h ). The population is decreased by natural death 765 and is given by The population of new female sandflies (S v ) is increased by an adult recruitment rate (λ v ) and decrease is described by The population of infected female sandflies is generated at the per-capita rate λ hv and diminished by the 771 natural death rate µ v . Thus,

774
For simplification, we let G 1 = φ h + µ h , G 2 = θ h + µ h and G 3 = γ h + µ h . Considering the infected subpopulations I h (t), A h (t), and I v (t), we let F be the rate of new infections into the infected compartments and V be the rate of exit of humans into infected compartments: of rates of inflow of new infections in each compartment and V = V + + V − is the vector of rates 776 transfer rates of individuals into and out of the infective compartments by all other processes. Taking   777 the Jacobian matrix of each vector with respect to each of the infectious classes and evaluating at Computing FV −1 , we obtain Taking the spectral radius of the next generation matrix operator, ρ FV −1 , gives  Proof. Because this model is of epidemiological relevance, we first show that the region Ω is positively 788 invariant in R 7 + , with respect to the system (1) and (2). It is easy to see  (1) is

B.2 Positivity and Boundedness of Solutions
Hence, on applying a (comparison) theorem from Birkhoff 793 and Rota ( [10]) on differential inequality, we get Similarly, let (S v , I v ) ∈ R 2 + be the solution with non-negative initial solution. Taking the time 796 derivative along the sum of all solutions curves of model (2) gives By differential inequality theorem in [10], we find where N v (0) represents the initial sandfly population at the initial phase of the disease. As t → ∞, the 799 inequality becomes Proof. Linearization at DFE gives The characteristic polynomial of the Jacobian matrix J(E 0 ) is given by where Thus by Routh-Hurwitz criteria, E 0 is locally asymptotically stable for R C < 1 and is unstable for asymptotically stable in Ω whenever R C < 1 and unstable if R C > 1.

821
Proof. Consider a candidate Lyapunov function defined in Ω, where the constants L i , i = 1...5 are taken to be The function L is positive definite, in the sense that it vanishes only at the disease-free equilibrium while 823 otherwise it is positive in Ω. Moreover, taking the time derivative of the function in (25) along solutions 824 of system 1-2 and then substituting the expression for the derivatives, gives Substituting the L i constants in equation 26 and then grouping and collecting terms, giveṡ The first two terms are negative, as the arithmetic mean is greater than or equal to the geometrical mean.

827
However, the third term is negative for values of R C < 1. Therefore, by Lyapunov-LaSalle asymptotic 828 stability [52], the disease-free equilibrium E 0 is globally asymptotically stable if R C < 1 for all t > 0.

830
As a result of no disease deaths, observeD in Figure ??, the existence of a DFE and an Endemic Equilib-831 rium (EE) that depends on R C . In this section, we show the local and global stability of the EE when 832 R * C become 1.

835
Proof. The EE of the Model system equations 1-2 is given by E * . The Jacobian matrix at EE gives by It's characteristic polynomial is given by We observe that the characteristic polynomial P (λ) can be factored to roots λ = −µ h , −µ v , −G 3 and 838 P (λ) = λ 4 + h 3 λ 3 + h 2 λ 2 + h 1 λ + h 0 . Applying the Routh-Hurwitz conditions: h i > 0, (i = 0, ..., 4), hold when R C > 1. Thus, the endemic equilibrium, E * , is locally asymptotically stable because all 841 eigenvalues of the septic polynomial have all negative real parts for R C > 1.

845
Proof. Consider a candidate Lyapunov function defined in Ω, where the constants L i , i = 1...5 are given by Taking the time derivative of the Lyapunov function in (28) along solutions of system 1-2 and then substituting the expression for the derivatives giveṡ Substituting the L i in 29 and performing some algebra giveṡ The first two terms in parenthesis and the remaining expression are negative, as the arithmetic mean 846 is greater than or equal to the geometrical mean. Therefore, by LaSalle's Invariable Principle [52], the 847 endemic equilibrium point E * is globally asymptotically stable in Ω for R 0 > 1 for all t > 0.

849
After extensive searching of the literature, annual reports, and census data, ecological and epidemiological 850 parameter ranges for the respective human and sandfly populations in India and Sudan were gathered 851 and estimated. See Table 2 for a summary of these estimates. India Parameter estimates were generated by solving for β vh in our R C expression and then pairing samples of known values in Table 2

India
Year