\begin{matrix}c \\ ◯ \end{matrix} Revista de Matemática: Teoría y Aplicaciones 2020 27(1) : 221–239 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v27i1. 39973

modeling post-kala-azar dermal leishmaniasis as an infection reservoir

for visceral leishmaniasis

leishmaniasis dérmica post-kala-azar

como reservorio de infección para la leishmaniasis visceral

Andrea Calderon^{∗}  Ryan Landrith^{†}  Nhan Le^{‡}
Ileana Muñoz^{§}  Christopher M. Kribs^{¶}

Received: 16/May/2019; Revised: 6/Jun/2019; Accepted: 13/Sep/2019

∗University of Texas at Arlington, Department of Biology, Arlington TX, United States. E-Mail: andrea.calderon@mavs.uta.edu

†Misma dirección que/Same address as: A Calderon. E-Mail: ryan.landrith@mavs.uta.edu ‡University of Texas at Arlington, Department of Computer Science and Engineering, Arling-

ton TX, United States. E-Mail: nhan.le@mavs.uta.edu

§University of Texas at Arlington, Department of Mathematics, Arlington TX, United States. E-Mail: ileana.munoz@mavs.uta.edu

¶University of Texas at Arlington, Departments of Mathematics and Curriculum & Instruction, Arlington TX, United States. E-Mail: kribs@uta.edu

222 a. calderon – r. landrith – n. le – i. muñoz – c. kribs

Abstract

Visceral Leishmaniasis (VL) is a potentially fatal disease caused by the protozoan parasite Leishmania donovani. This disease is a health problem for the very poor because it results in thousands of deaths and illnesses every year. Some countries, such as India and Bangladesh, have started programs to reduce the occurrences of VL by focusing on early diagnosis and complete treatment of VL. Post-Kala-azar Dermal Leishmaniasis (PKDL) is a cutaneous manifestation of Leishmaniasis that can occur following the incomplete treatment of VL. Diagnosis and treatment of PKDL are limited in affected regions, and PKDL has been identified as a possible reservoir for infection. This study develops a mathematical model of the relationship between the level of PKDL treatment and the incidences of VL during a given period. The results indicate a nearly linear relationship between PKDL treatment rates and the percent reduction of VL incidences. With the current treatments available and considering achievable levels of treatment, the model predicts that up to 20% of VL cases could be prevented by treating new PKDL cases. Hypothetical combined treatment initiatives including bed nets and insecticide spraying are also considered. Results suggest that the population of individuals with PKDL is certainly a significant factor in the transmission of L. donovani infection, with treatment of new cases particularly important.

Keywords: disease reservoir; experimental treatment; sandfly vector.

Resumen

La leishmaniasis visceral (VL) es una enfermedad potencialmente mortal causada por el parásito protozoario Leishmania donovani. Esta enfermedad es un problema de salud para los muy pobres porque da lugar a miles de muertes e infecciones cada año. Algunos países, como la India y Bangladesh, han comenzado programas para reducir las ocurrencias de VL centrándose en el diagnóstico precoz y el tratamiento completo de la VL. La leishmaniasis dérmica post-kala-azar (PKDL) es una manifestación cutánea de la leishmaniasis que puede ocurrir después del tratamiento incompleto de la VL. El diagnóstico y el tratamiento de PKDL son limitados en las regiones afectadas, y PKDL se ha identificado como un posible reservorio para la infección. Este estudio desarrolla un modelo matemático de la relación entre el nivel de tratamiento de PKDL y las incidencias de VL durante un periodo determinado. Los resultados indican una relación casi lineal entre las tasas de tratamiento de PKDL y la reducción porcentual de incidencias de VL. Con los tratamientos actualmente disponibles y teniendo en cuenta los niveles alcanzables de tratamiento, el modelo predice que hasta 20% de los casos de VL podrían prevenirse mediante el tratamiento de nuevos casos PKDL. También se tienen en cuenta iniciativas de tratamiento combinadas hipotéticas, como mosquiteros y pulverización de insecticidas. Los resultados sugieren que la población de individuos con PKDL es sin duda un factor importante en la transmisión de L. donovani, con el tratamiento de nuevos casos particularmente importantes.

Palabras clave: reservorio de enfermedad; tratamiento experimental; vector flebótomo.

Mathematics Subject Classification: 92C37.

1 Introduction

Visceral Leishmaniasis (VL), also known as Kala-azar, is a fatal disease caused by the Leishmania donovani protozoan parasite and is characterized by fever, weight loss, hepatosplenomegaly, and pancytopenia [19]. The parasite is transmitted by the bite of the vector Phlebotomus argentipes, commonly known as the sandfly. VL is a health problem among the very poor; 90% of VL cases arise in Bangladesh, Nepal, Sudan, and India, where the disease is dire in the eastern and more rural parts of the country, causing the deaths of thousands and the severe sickening of hundreds of thousands every year [19]. Post-Kala-azar Dermal Leishmaniasis (PKDL) is a cutaneous manifestation of leishmaniasis following treatment of VL, characterized by skin lesions and nodules or papules which are often to be found on the face [13]. PKDL is not a life threatening disease, and the treatment of this disease is considered to be a burden by many of the affected. For example, in Bangladesh, current treatment guidelines call for 120 intramuscular injections of sodium stibogluconate. Hence, many patients remain undiagnosed and untreated [23].

In India, humans are considered the primary reservoir for the L. donovani parasite due to the high population density making the transmission from human to human via sandflies (anthroponosis) common. For the purposes of this study, animals are not considered a significant reservoir in India [22]. Similar to India, a major reservoir in Sudan for the parasite is in humans; however, the findings of infected flies in the uninhabited Dinder Park strongly suggest the presence of a reservoir other than man. L. donovani has been detected in dogs and some other animals including livestock [26]. As a result of extended exposure time and reduced attention relative to VL, PKDL-infected humans have been posited as a reservoir for sandfly infection.

VL has been the focus of concerted treatment efforts due to its severity. In 2005, the Indian, Nepalese, and Bengali governments initiated a plan to reduce the occurrence of VL to less than 0.01% by 2015. They planned to do this by focusing on early diagnosis and complete treatment (treatment-related control strategies) and spraying insecticides in homes (vector-related control strategy) [23]. Four drugs are available to treat VL. These drugs include pentavalent antimonials, which have been the “first-line” treatment for 70 years, but are said to be toxic and accompanied by failure rates due to drug resistance; miltefosine, which is the first oral treatment against VL, but is said to lead to resistance because of its long half-life; amphotericin B, used in conventional and liposomal formulations, yet too expensive and complex to be used on a large scale; and paromomycin (PMM), which is currently being tested in a Phase IV trial in India [23].

Treatment of PKDL, in contrast, is less developed. Diagnosis of PKDL is based on a history of VL, distribution and appearance of lesions, and by parasitological confirmation when the diagnosis is doubtful [8]. In Sudan, some patients that were diagnosed with PKDL underwent treatment with intravenous sodium stibogluconate [13]. The result of the study showed a complete disappearance of any indication of PKDL. Treatment with intravenous sodium stibogluconate varied among patients due to their reaction or prior diseases. Some patients were treated as planned for 30 days and were healed, some patients needed to be treated with ketoconazole for 30 days and then put back on treatment with sodium stibogluconate before they were healed, and for some patients, even after being returned to regular treatment received a higher dose of the treatment [13]. Another study was reported with a 26 year old Ethiopian man (patient 1) and a 42 year old Ethiopian man (patient 2). Both patients were given treatment options: intravenous sodium stibogluconate or oral miltefosine. The patients chose to try miltefosine. Patient 1 was cured after six months of treatment with side effects during the elapsed time. Patient 2 was cured after three months of treatment with no experienced side effects [2]. This was the first reported use of miltefosine in the treatment of PKDL; according to Belay et al. “Miltefosine appears to be a promising treatment for PKDL, and its use in this context merits further investigation” [2, p. 226]. These treatments of PKDL can serve to reduce the reservoir for the VL disease.

In this study, we will use a dynamical system model to investigate the question: What is the relationship between VL incidence during a given period and the level of PKDL treatment during the same period? We assume that the infection reaches an endemic state and define the proportion of PKDL cases treated as the ratio of PKDL cases treated to the incidences of PKDL over the time period.

2 Mathematical model

The model’s dynamics are based on the epidemiologically pertinent stages of the disease. Diagnostically the stages of VL can be separated by corresponding levels of parasitemia [23]. VL has a SIRS infection cycle structure in hosts, since recovery confers temporary immunity, with some recovered hosts spontaneously developing PKDL. Spontaneous development of PKDL by asymptomatic individuals was neglected from the model because the reported incidences are insignificantly low. Vectors undergo an SEI cycle with a significant latent period during which the parasite reproduces in the gut. Due to elapsed time for the parasite to grow in the gut and be transmitted, in addition to the SIR model, a vector category will include an exposed (E) vector category. The model is compartmentalized into five human classes and three classes of sandflies listed in Table 1 related to L. donovani infection as summarized in the flow chart (Figure 1). The flow chart reflects the natural cycle of the disease with the arrows between each state variable representing per capita rates of transition between categories.

Figure 1: Flow diagram for human compartments.

Figure 2: Flow diagram for vector compartments.

Table 1: Human and sandfly state variables.

Variable Definition
S  Susceptible humans
I_{A}  Asymptomatic VL-infected humans
I_{S}  Symptomatic VL-infected humans
R  Recovered humans
P  PKDL-infected humans
S_{v}  Susceptible vectors
E_{v}  Exposed vectors
I_{v}  Infectious vectors

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 221–239, Jan–Jun 2020

The model is a system of ordinary differential equations with one equation representing each class or state variable. Each transition rate from the flow chart represents a term in the respective equations. Though India is experiencing roughly exponential human population growth, birth rates into the susceptible classes are represented as constant rates [4]. This simplifying assumption is made since the time period to be simulated is short and to allow the system to reach an endemic equilibrium at which the efficacy of PKDL treatment can be evaluated. A synthetic birth rate also allows maintenance of a realistic sandfly to human density ratio and provides conditions consistent with vector dependent transmission as indicated by the literature [16]. The infection rates were developed assuming that the transmission of the disease is dependent on the density of the sandfly population and not the human population [16]. Thus the terms are a product of the contact rate (c), probability of infection per contact (β_{v}, σ_{A}, σ_{P} , σ_{S}), density of sandflies involved for the respective transmission (susceptible or infected), and the proportion of humans in the state variable involved. The rate of PKDL treatment (ψ) is a parameter varied in the numerical simulations. The other simple per capita rates (ω, θ, γ, ϕ, ρ) are derived as the inverse of the average time spent in a state before the respective transition or as the product of the inverse of the average time spent in a state regardless of the transition out and the fraction of individuals in that state who will make that particular transition. Table 2 summarizes the model parameters and their estimated values (for further details see Appendix A).

The resulting system is given below. The last three equations (\frac{dT}{dt} , \frac{dZ}{dt} , \frac{dV}{dt} ) are used to tally the number of cases of PKDL treated (T), the new incidences of PKDL (Z), and the new symptomatic incidences of VL (V) respectively. Totaling the number of new symptomatic VL cases and the proportion of new PKDL cases treated over the simulated time period allows for evaluation of a relationship between the cases treated and symptomatic VL incidences prevented as presented in the results section.

\frac{dS}{dt} =Λ_{H} + ωR − β_{v}cS \frac{I_{v}}{N_{H}} − ξS,  (1)
\frac{dI_{A}}{dt} =β_{v}cS \frac{I_{v}}{N_{H}} − (θ + γ + ξ)I_{A},  (2)
\frac{dI_{S}}{dt} =γI_{A} − (ϕ + δ + ξ)I_{S},  (3)
\frac{dR}{dt} =θI_{A} + ϕI_{S} + ψP − (ρ + ω + ψ)R,  (4)

\frac{dP}{dt} =ρR − (ψ + ξ)P,  (5)
\frac{dS_{v}}{dt} =Λ_{v} − c(σ_{A}I_{A} + σ_{S}I_{S} + σ_{P} P)\frac{S_{v}}{N_{H}} − μS_{v},  (6)
\frac{dE_{v}}{dt} =c(σ_{A}I_{A} + σ_{S}I_{S} + σ_{P} P)\frac{S_{v}}{N_{H}} − (k + μ)E_{v},  (7)
\frac{dI_{v}}{dt} =kE_{v} − μI_{v},  (8)
\frac{dT}{dt} =ψP,  (9)
\frac{dZ}{dt} =ρR,  (10)
\frac{dV}{dt} =γI_{A}.  (11)

Table 2: Model parameters and estimated values from the literature.

Parm. Definition  Unit Est. Value Source
Λ_{H} Human reproductive rate  \frac{humans}{sq.km. day} 0.0309  [15]
ξ Natural death rate for humans  1/day  0.0079  [15]
N_{H} Population density of Bihar, India  \frac{humans}{sq.km.} 1102.39  [10]
γ Progression rate from asymptomatic  1/day  0.0006  [24]

to symptomatic

θ Spontaneous recovery from asymptomatic 1/day  0.0139 [1, 11, 23]
δ Death rate due to VL (with treatment)  1/day  0.13  [20]
ϕ Rate of recovery with treatment for VL  1/day  0.0306  [23, 25]
ρ Rate of development to PKDL  1/day  0.0028  [9]
ψ PKDL treatment rate  1/day  0.981  [6]
ω Rate of immunity loss  1/day  0.0135  [23]
Λ_{v} Reproductive rate of vectors  \frac{vec.}{sq.km. day} 195.6003 [12, 21]
μ Natural death rate of vectors  1/day  0.0714  [22]
N_{V} Density of female sandflies  \frac{vec.}{sq.km.}  987.63 [12, 21, 23]
c Biting rate of flies  1/day  0.25  [14]
σ_{A} Probability of sandfly infection from  0.01458  [23]

biting asymptomatic VL-infected humans

σ_{S} Probability of sandfly infection from 1 Assumed

biting symptomatic VL-infected humans

σ_{P} Probability of sandfly infection from 1 Assumed

biting PKDL-infected humans

k Rate of sandflies becoming infectious 1/day 0.2 [18]

after exposure

β_{v} Infectivity of vectors  1/day  1  Assumed
Ratio of sandflies to humans  527:100  [23]

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 221–239, Jan–Jun 2020

3 Analysis

The control reproduction number R_{c} provides a measure of the persistence of L. donovani infection or lack thereof, with a theoretical PKDL treatment rate in the model. By evaluating at the disease-free equilibrium it can be determined whether transmission of the infection will be sustained when the number of infected individuals is close to zero. This condition is equivalent to the infection persisting at some endemic level in the population. When the value of R_{c} is 1 or greater the infection will persist in the population.

At the disease-free equilibrium all classes are empty except for the wholly uninfected groups S and S_{v}. Here equations (1) and (6) can be solved to find the equilibrium population densities, S^{∗} = N_{H}^{∗} = Λ_{H}/ξ and S_{v}^{∗} = N_{v}^{∗} = Λ_{v}/μ.

To compute R_{c} we use standard next-generation methods, which identify the control reproduction number as the dominant eigenvalue (spectral radius, or largest growth factor) of the next-generation matrix, which we calculate (see Appendix B) as

[ \begin{matrix}0 & 0 & 0 & \frac{cβ_{v}}{μ} \\ \frac{γ}{γ+θ+ξ} & 0 & 0 & 0 \\ \frac{θ}{γ+θ+ξ} \frac{ρ}{ρ+ω+ξ} & \frac{ϕ}{ϕ+δ+ξ} \frac{ρ}{ρ+ω+ξ} & 0 & 0 \\ \frac{k}{k+μ} \frac{cσ_{A}}{γ+θ+ξ} \frac{Λ_{v}/μ}{Λ_{H}/ξ} & \frac{k}{k+μ} \frac{cσ_{S}}{ϕ+δ+ξ} \frac{Λ_{v}/μ}{Λ_{H}/ξ} & \frac{k}{k+μ} \frac{cσ_{P}}{ξ+ϕ \frac{ω+ξ}{ρ+ω+ξ}} \frac{Λ_{v}/μ}{Λ_{H}/ξ} & 0 \end{matrix} ] .

Although we cannot write a closed-form expression for R_{c} as the root of a quartic equation, this matrix can be used to compute R_{c} during numerical analysis.

When R_{c} < 1, the disease-free equilibrium is locally asymptotically stable. Numerical explorations indicate that when R_{c} > 1 a unique endemic equilibrium exists and is stable.

3.1 Numerical methods

Simulations and graphics were produced in MATLAB version 7.0.4.365 using the built in ordinary differential equation solver ode45. Initial conditions were calculated using the prevalence for each category based on parasitemia level and population density given in Table 2. The treatment rate of PKDL varies from none to 0.0333(day)^{−1}, a rate which corresponds to an average time of seeking treatment and becoming non-infectious of 30 days. This range was chosen based on the treatment duration of PKDL which lasts 30 days so scenarios are considered in which the average time to receive treatment ranges from never (infinity) to 30 days. We realize that patients may be treated to the point of non-infectivity sooner than 30 days; however, we assume 30 days to be a reasonable lower bound considering any delay in receiving treatment.

Preliminary results indicated that values of the state variables reached within one tenth of one percent of their equilibrium values within two years of simulation. All simulations presented in the results were run for two years. This condition functions well under the assumption of a constant population density since the density of Bihar is not likely to increase significantly over two years.

Scenarios with hypothetical sandfly control measures were also simulated. These simulations were performed at ψ = 0.025/day which corresponds to an average time of treatment until noninfectiousness of 40 days. This value showed the typical dynamics of the treatment rates considered and is still an optimistic estimate, but not unreasonably so. For these simulations two parameters were varied: the contact rate (c) and the death rate of sandflies (μ). These represent varying levels of implementation of bed nets and indoor residual insecticide spraying respectively.

One measure of the level of treatment presented is the proportion of new incidences of PKDL treated over the two year simulation. This reflects a hypothetical treatment initiative where individuals who develop PKDL either seek treatment or are reexamined at regular intervals for PKDL following treatment for VL.

4 Results

Under the parameter estimates given in Table 2, R_{c} > 1, so numerical work explored the ability of various control measures (primarily PKDL treatment) to reduce the spread of L. donovani infection. Two series of simulations are presented. In the first series (Figures 3 through 7), the simulation is run with the treatment rate of PKDL (ψ) varying each simulation from none to 0.0333(day)^{−1}. In the second series (Figure 8) the rate of treatment, ψ = 0.025/day, is held constant, and the parameters representing the biting rate of sandflies (c) and their death rate(μ) are varied for each simulation instead.

Figures 4, 5 and 7 represent the level of treatment via the proportion of new incidences of PKDL treated over the two year simulation. These reflect the impact of the hypothetical treatment initiative. Figures 3 and 6 represent the treatment rate as it is in the equations of the model, the per capita rate ψ.

Figure 3 graphs the treatment rate ψ against the percent reduction of total VL cases. Though this curve may appear roughly linear at these values, the change in concavity is not numerical error, and for a larger range of treatment rates the curve assumes a more pronounced sigmoid shape.

Figure 3: Percent reduction in VL cases as a function of PKDL treatment rate (in 1/day).

Figure 4: Percent reduction in VL cases as a function of proportion of new PKDL cases treated.

Figure 4 shows the proportion of new PKDL cases treated plotted against the percent of total symptomatic VL incidences averted over the two years (referred to as percent reduction) for each treatment level. This curve would level off at extremely high treatment rates; however, in the range of reasonable treatment rates considered, the relationship is roughly exponential.

Figure 5: Proportion of new PKDL cases treated as a function of PKDL treatment rate (in 1/day).

Figure 6: Control reproduction number as a function of PKDL treatment rate (in 1/day).

Figure 5 plots the proportion of new incidences of PKDL cases that were treated into remission during the two year period against the control reproductive number R_{c}, the value of which actually depends on the value of ψ for each simulation (through the next-generation matrix, q.v.). Similarly, Figure 6 shows the treatment rate ψ graphed against the control reproduction number R_{c}.

Figure 7: Control reproduction number as a function of proportion of new PKDL cases treated.

Figure 8: Effect of sandfly control measures on control reproduction number.

Figure 7 demonstrates the effect of the treatment rate on the proportion of new PKDL cases treated.

Figure 8 illustrates the effect of varying rates of sandfly control measures on the value of RC. Hypothetical controls were simulated deviating from the natural values of contact rate (c) and the death rate of sandflies (μ). These values correspond to the average amount of time between bites per sandfly (feeding cycle duration) which is normally 4 days and the average lifespan of sandflies which is naturally 14 days.

In the first series of simulations, a treatment rate of ψ = 0.08/day yielded R_{c} = 0.9999; however, this level of treatment corresponds to an average time of treatment to the point of non-infectivity of 12.5 days, which is not possible given current treatments considered for PKDL.

The simulations also show that based on our parameter estimates the maximum percent reduction possible, assuming treatment is sought immediately after developing symptoms and takes 30 days to complete, is 20.2% (Figure 5). Likewise, with PKDL treatment alone the most optimistic reduction of R_{c} was to a value of 1.18.

In the second series, the simulations of sandfly control measures (Figure 8) showed that as little as a 24.7% increase in the feeding cycle duration (from 4 to 5 days) and a 13.2% reduction in the average lifespan of a sandfly (from 14 to 12.14 days) could check the spread of L. donovani infection (bringing R_{c} < 1).

5 Discussion

Conclusions

Figure 3 illustrates an important result of the model, showing that with the current treatments available and considering only achievable levels of treatment, the impact of treating more new PKDL cases on symptomatic VL cases does not diminish at higher proportions of cases treated. Likewise, the graph of R_{c} vs. the proportion of new PKDL cases treated (Figure 6) indicates that effect of treatment on the spread of the disease within the range considered will continue to increase within reasonable treatment ranges.

Conversely, Figures 4 and 5 show that though the effect of PKDL treatment on symptomatic VL incidence increases at the higher proportions of new PKDL cases treated, there is a diminishing return of this proportion with increased treatment rate ψ indicating how high the actual treatment rate would have to be to capture a large proportion of new PKDL cases. The opposite concavities of Figures 6 and 7 display the same effect. The effect of an increasing proportion of cases treated on R_{c} produces increasingly higher reduction of R_{c} whereas the effect of higher treatment rates diminishes.

These results indicate that public health initiatives aimed at preventing or decreasing the spread of L. donovani will be more effective when focusing on new cases by both raising awareness about the epidemiological effect of seeking treatment and keeping track of individuals who received treatment for VL.

The series of sandfly control simulations indicated that PKDL treatments in tandem with achievable vector control initiatives may produce even more pronounced effects on the transmission of L. donovani, even to the point of reducing the spread of the disease into decline.

Limitations

The particular values reported may not accurately reflect the actual results if these treatments were implemented the large number of parameters estimated from independent sources certainly produced some numerical inaccuracies in the simulations, though these inaccuracies should not extend into the general behavior of the model assuming these estimate are close to the actual values. However, the results certainly produce meaningful insight into the relationship between PKDL treatment and symptomatic VL incidence and corroborate the conclusion that high levels of PKDL treatment will significantly impact the rate of L. donovani infection. Many assumptions may have distorted the accuracy of the model. It should be noted that among other limiting assumptions this model is only applicable to regions where L. donovani infection has no substantial nonhuman reservoir, such as the endemic areas of India where human population density is high.

Further study

This study suggests that there is a significant epidemiological impact on L. donovani infection in humans by treating PKDL, but more research is required to determine the precise level of treatment necessary. Further study is required to fit models like this to actual data involving varying treatment levels of PKDL to determine both goodness of fit and fit the model to real rates of impact. Clinical studies to determine how infectious PKDL patients are may improve understanding of how significant a reservoir for transmission PKDL is.

Acknowledgements

This research was supported by an NSF UBM-Institutional grant, DUE#0827136, as part of the UTTER program at UT Arlington.

Appendix A. Parameter estimates

The initial values for the five human classes were taken from [23], in which the authors tested a population for VL using two different tests: PCR (polymerase chain reaction) and DAT (direct agglutination test). The PCR test is used to replicate a small segment of DNA into larger amounts; this is then used to determine whether or not the bacteria are present. The test will result positive if the bacteria are present. The DAT test is used to determine the presence of antibodies to a specific antigen. If the DAT test results positive the antibodies are present and if it results negative the antibodies are not present.

For susceptible human hosts, denoted in Figure 1 as S, the value was taken from the percentage of cases in which both tests resulted negative, in this case 76%. The proportion of asymptotic infected individuals, denoted in Figure 1 as I_{A}, was taken from the percentage of individuals who tested positive for PCR without having symptoms, in this case 11.985%. The value for symptomatic individuals, denoted in Figure 1 as I_{S}, was taken from the percentage of individuals who tested positive for PCR or DAT while showing symptoms, in this case 0.015%. The value for recovered individuals, denoted in Figure 1 as R, was taken from the cases that resulted PCR negative and DAT positive. The PCR test with negative results determines the current infection status of the individual, whereas the DAT test resulting positive indicates prior exposure. This means that the individual does not currently have visceral leishmaniasis but did at one point. The value for the prevalence of recovered individuals who develop PKDL, denoted in Figure 1as P, was taken by simply adding the total prevalence from individuals treated with the first line treatment to the total percentage of prevalence from individuals treated with the second line treatment, which amounted to 6%.

The human population density in Bihar, India, which is denoted in Figure 1 as N_{H}, was calculated with information from [10]. We divided the total population of Bihar, India by the area in square kilometers. The estimated density of female sandflies, denoted as N_{V}, was then calculated as follows. We multiplied the total human population density by the ratio of vectors per person found in [23]. This value gave us the vector population density in Bihar, India, which we then multiplied by the proportion of male to female sandflies given in [21]. We then divided the female population density by vector lifespan [22] to estimate the vector reproductive rate, denoted in Figure 2 as Λ_{V}.

The value for susceptible vectors, denoted in Figure 2 as S_{V}, was taken from [17] by subtracting the total percentage of vectors already infected with VL from 100%, which yielded 96.6%. The value for the prevalence of exposed sandflies, which is denoted in Figure 2 as E_{V}, was also calculated with information from [17]. We subtracted the prevalence of infectious vectors from the total percentage of vectors infected with VL, which amounted to 2.9%. The value for the prevalence of infectious vectors, which is denoted in Figure 2 as I_{V}, was taken from [3] as 0.5%.

The human death rate, denoted in Figure 1 as ξ, was calculated using information using the death rate per 1,000 individuals given in [15]. The human birth rate, denoted in Figure 1 as Λ_{H}, was also estimated from [15], from the birth rate per 1,000 individuals.

The progression rate from asymptomatic cases to symptomatic cases, which is denoted in Figure 1 as γ, was calculated from [24], which gave the rate of cases that progressed from asymptomatic to symptomatic was given in terms of 1,000 people per months. In order to obtain the rate in one person per day we divided this by 1,000, multiplied it by the number of months in a year, and then divided it by the total number of days in one year.

The rate of spontaneous recovery from asymptomatic VL, denoted in Figure 1 as θ, is calculated with information from [23]. First we added the sojourn time in the early asymptomatic stage to the sojourn time in the late asymptomatic stage to get the total time an individual spends in the asymptomatic stage. We then multiplied the inverse of the time spent in the asymptomatic stage by the fraction of cases that spontaneously recover from asymptomatic VL.

The VL recovery rate due to treatment, denoted in Figure 1 by ϕ, was calculated with information from [25]. We multiplied the fraction of individuals who respond to treatment against Visceral Leishmaniasis by the inverse of the duration of treatment. The death rate due to Visceral Leishmaniasis when treated, denoted in Figure 1 as δ, was taken from [20].

The PKDL development rate, denoted in Figure 1 as ρ, was calculated with information from [9], based on the rate given in terms of per 180 days. The PKDL recovery rate due to treatment, denoted in Figure 1 as ψ, was taken from [6].

The rate of loss of immunity, which is denoted in Figure 1 as ω, was calculated from [23] as the inverse of the mean duration of the period when the DAT tests are positive.

The probability of the vector becoming infected after biting an asymptomatic VL-infected human, denoted in Figure 2 as σ_{A}, was calculated with information from [23], as a weighted average of (a) the probability that a susceptible vector becomes infected when feeding on a human host in the early asymptomatic stage and (b) the probability that a susceptible vector becomes infected when feeding on a human host in the late asymptomatic stage, where the respective weights are the amount of time spent in the corresponding asymptomatic stage divided by the total time spent in the asymptomatic stage. The value of the relative infectivity of symptomatic VL-infected humans, denoted in Figure 2 as σ_{S}, was taken from [23] and is assumed to be 1. The probability of vectors becoming infected after biting a PKDL-infected human, denoted in Figure 2 as σ_{P}, was taken from [23]. This value is assumed to be 1.

Appendix B. R_{c} derivation

The control reproductive number R_{c} can be derived using standard next-gener-ation methods, e.g., [7]. We begin by decomposing the state vector (1–8) into three parts: uninfected classes X = {S, S_{V} }, noninfectious infected classes Y = {R, E_{V} } (we must include R due to the spontaneous development of PKDL from R), and infectious classes Z = {I_{A}, I_{S}, P, I_{v}}. We then compute the submatrix A of the Jacobian of the system involving the infectious variables Z, with the classes in Y replaced by their equilibrium expressions (obtained by setting dY/dt = 0 and solving for Y in terms of X and Z) and the classes in X and Z evaluated at the disease-free equilibrium (DFE):

A = [ \begin{matrix}\frac{∂}{∂I_{A}} (\frac{dI_{A}}{dt} ) & \frac{∂}{∂I_{S}} (\frac{dI_{A}}{dt} ) & \frac{∂}{∂P} (\frac{dI_{A}}{dt} ) & \frac{∂}{∂I_{V}} (\frac{dI_{A}}{dt} ) \\ \frac{∂}{∂I_{A}} (\frac{dI_{S}}{dt} ) & \frac{∂}{∂I_{S}} (\frac{dI_{S}}{dt} ) & \frac{∂}{∂P} (\frac{dI_{S}}{dt} ) & \frac{∂}{∂I_{V}} (\frac{dI_{S}}{dt} ) \\ \frac{∂}{∂I_{A}} (\frac{dP}{dt} ) & \frac{∂}{∂I_{S}} (\frac{dP}{dt} ) & \frac{∂}{∂P} (\frac{dP}{dt} ) & \frac{∂}{∂I_{V}} (\frac{dP}{dt} ) \\ \frac{∂}{∂I_{A}} (\frac{dI_{V}}{dt} ) & \frac{∂}{∂I_{S}} (\frac{dI_{V}}{dt} ) & \frac{∂}{∂P} (\frac{dI_{V}}{dt} ) & \frac{∂}{∂I_{V}} (\frac{dI_{V}}{dt} ) \end{matrix} ]_{DF E}

= [ \begin{matrix}−(ξ + γ + θ) & 0 & 0 & cβ_{V} \\ γ & −(ξ + δ + ϕ) & 0 & 0 \\ \frac{θρ}{ξ+ρ+ω} & \frac{ϕρ}{ξ+ρ+ω} & − (ξ + \frac{ψ(ξ+ω)}{ξ+ρ+ω}) & 0 \\ \frac{kcξΛ_{V} σ_{A}}{Λ_{H}μ(k+μ)} & \frac{kcξΛ_{V} σ_{S}}{Λ_{H}μ(k+μ)} & \frac{kcξΛ_{V} σ_{P}}{Λ_{H}μ(k+μ)} & −μ \end{matrix} ] .

We decompose A = M − D where M is nonnegative and D is diagonal:

M = [ \begin{matrix}0 & 0 & 0 & cβ_{v} \\ γ & 0 & 0 & 0 \\ \frac{θρ}{ξ+ρ+ω} & \frac{ϕρ}{ξ+ρ+ω} & 0 & 0 \\ \frac{kcξΛ_{V} σ_{A}}{Λ_{H}μ(k+μ)} & \frac{kcξΛ_{V} σ_{S}}{Λ_{H}μ(k+μ)} & \frac{kcξΛ_{V} σ_{P}}{Λ_{H}μ(k+μ)} & 0 \end{matrix} ] ,

D = [ \begin{matrix}(ξ + γ + θ) & 0 & 0 & 0 \\ 0 & (ξ + δ + ϕ) & 0 & 0 \\ 0 & 0 & (ξ + \frac{ψ(ξ+ω)}{ξ+ρ+ω}) & 0 \\ 0 & 0 & 0 & μ \end{matrix} ] .

Then R_{c} is given as the spectral radius (dominant eigenvalue, or maximal growth ratio) of the next-generation matrix

MD^{−1} =[ \begin{matrix}0 & 0 & 0 & \frac{cβ_{v}}{μ} \\ \frac{γ}{γ+θ+ξ} & 0 & 0 & 0 \\ \frac{θ}{γ+θ+ξ} \frac{ρ}{ρ+ω+ξ} & \frac{ϕ}{ϕ+δ+ξ} \frac{ρ}{ρ+ω+ξ} & 0 & 0 \\ \frac{k}{k+μ} \frac{cσ_{A}}{γ+θ+ξ} \frac{Λ_{v}/μ}{Λ_{H}/ξ} & \frac{k}{k+μ} \frac{cσ_{S}}{ϕ+δ+ξ} \frac{Λ_{v}/μ}{Λ_{H}/ξ} & \frac{k}{k+μ} \frac{cσ_{P}}{ξ+ϕ \frac{ω+ξ}{ρ+ω+ξ}} \frac{Λ_{v}/μ}{Λ_{H}/ξ} & 0 \end{matrix} ] .

References

[1] Ad Hoc Committee on Health Research Relating to Future Intervention Options, Investing in Health Research and Development. World Health Organization, Geneva, 1996 (Document TDR/Gen/96.1). Available at:

    https:// apps. who. int/ iris/ handle/ 10665/ 63024

[2] A.D. Belay, Y.Y. Asafa, J.J. Mesure, R.N. Davidson, Successful miltefosine treatment of post-kala-azar dermal leishmaniasis occurring during antiretroviral therapy, Ann. Trop. Med. Parasitol. 100(2006), no. 3, 223–227.

doi: 10.1179/ 136485906X91440

[3] N.R. Bhattarai, M.L. Das, S. Rijal, G. van der Auwera, A. Picado, B. Khanal, L. Roy, …, J.C. Dujardin, Natural infection of Phlebotomus argentipes with Leishmania and other trypanosomatids in a visceral leishmaniasis endemic region of Nepal, Transactions of the Royal Society of Tropical Medicine and Hygiene 103(2009), no. 11, 1087–1092. doi: 10.

1016/ j.trstmh. 2009. 03.008

[4] Central Intelligence Agency, The world factbook, https:// www.

cia. gov/ library/ publications/  the-world-factbook/
index. html

[5] F. Chappuis, S. Sundar, A. Hailu, H. Ghalib, S. Rijal, R.W. Peeling, J. Alvar, M. Boelaert, Visceral leishmaniasis: what are the needs for diagnosis, treatment and control?, Nature Reviews Microbiology 5(2007), no. 11,

873–882. doi: 10.1038/ nrmicro1748

[6] S. Collin, R. Davidson, K. Ritmeijer, K. Keus, Y. Melaku, S. Kipngetich, C. Davies, Conflict and kala-azar: determinants of adverse outcomes of kala-azar among patients in southern Sudan, Clinical Infectious Diseases 38(2004), no. 5, 612–619. doi: 10.1086/ 381203

[7] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R_{0} in models for infectious diseases in heterogeneous populations, Journal of Mathematical Biology 28(1990), no. 4, 365–382. doi: 10.1007/ BF00178324

[8] S. Ganguly, N.K. Das, J.N. Barbhuiya, M. Chatterjee, Post-kala-azar dermal leishmaniasis — an overview, International Journal of Dermatology 49(2010), no. 8, 921–931. doi: 10.1111/ j.1365-4632. 2010.

04558. x

[9] S. Gasim, A.M. Elhassan, A. Kharazmi, E.A.G. Khalil, A. Ismail, T.G. Theander, The development of post-kala-azar dermal leishmaniasis (PKDL) is associated with acquisition of Leishmania reactivity by peripheral blood mononuclear cells (PBMC), Clinical & Experimental Immunology 119(2000), no. 3, 523–529. doi: 10.1046/ j.1365-2249. 2000.

01163. x

[10] Government of India Ministry of Home Affairs, Office of Registrar General

& Census Commissioner, India, http:// www. censusindia.  gov.

in/2011-common/  censusdataonline.  html

[11] A. Hailu, M. Graiccia, P.A. Kager, Visceral leishmaniasis in Aba-Roba, south-western Ethiopia: prevalence and incidence of active and subclinical infections, Annals of Tropical Medicine and Parasitology 103(2009), no. 8,

659–670. doi: 10.1179/ 000349809X12554106963555

[12] J.G. Harre, K.M. Dorsey, K.L. Armstrong, J.R. Burge, K.E. Kinnamon, Comparative fecundity and survival rates of Phlebotomus papatasi sandflies membrane fed on blood from eight mammal species, Medical and Veterinary Entomology 15(2001), no. 2, 189–196. doi: 10.1046/ j.

1365-2915. 2001. 00278. x

[13] A.M. Hassan, H.W. Ghalib, E.E. Zijlstra, I.A. Eltoum, M. Satti, M.S. Ali, H.M.A. Ali, Post kala-azar dermal leishmaniasis in the Sudan: clinical features, pathology and treatment, Trans R Soc Trop Med Hyg 86(1992),

no. 3, 245–248. doi: 10.1016/ 0035-9203(92)90294-m

[14] A.K. Hati, S. Sur, N. De, H.N. Dwivedi, J. Bhattacharyya, H. Mukherjee, G. Chandra, Longitudinal study on distribution of phlebotomus argentipes sandflies at different heights in cattleshed, Indian Journal of Medical Research 93 (1991), 388–390.

[15] K. Nayar, A. Kumar, Health analysis – Kerala and Bihar: a comparison, Yojana 49(2005), no. 7, 9–13.

[16] B. Ostyn, V. Vanlerberghe, A. Picado, D.S. Dinesh, S. Sundar, F. Chappuis, S. Rijal, …, C. Davies, Vector control by insecticide-treated nets in the fight against visceral leishmaniasis in the Indian subcontinent, what is the evidence?, Tropical Medicine & International Health 13(2008), no. 8,

1073–1085. doi: 10.1111/ j.1365-3156. 2008. 02110. x

[17] K. Pandey, S. Pant, H. Kanbara, M. Shuaibu, A. Malik, B. Pandey, O. Kaneko, T. Yanagi, Molecular detection of Leishmania parasites from whole bodies of sandflies collected in Nepal, Parasitology Research 103(2008), no. 2, 293–297. doi: 10.1007/ s00436-008-0967-7

[18] D.L. Sacks, P.V. Perkins, Development of infective stage Leishmania promastigotes within phlebotomine sand flies, The American Journal of Tropical Medicine and Hygiene 34(1985), no. 3, 456–459. doi: 10.4269/

ajtmh. 1985. 34.456

[19] P. Salotra, R. Singh, Challenges in the diagnosis of post kala-azar dermal leishmaniasis, Indian J Med Res 123(2006), no. 3, 295–310.

[20] J. Seaman, D. Pryce, H.E. Sondorp, A. Moody, A.D. Bryceson, R.N. Davidson, Epidemic visceral leishmaniasis in Sudan: a randomized trial of aminosdidine plus sodium stibogluconate versus sodium stibogluconate alone, The Journal of Infectious Diseases 168(1993), no. 3, 715–720. doi:

10.1093/ infdis/ 168. 3.715

[21] N.S. Singh, D.P. Singh, Seasonal occurrence of phlebotominae sand flies (Phlebotominae: Diptera) and its correlation with kala-azar in eastern Uttar Pradesh, India, The Southeast Asian Journal of Tropical Medicine and Public Health 40(2009), no. 3, 458–462.

[22] R. Srinivasan, K.N. Panicker, Laboratory observations on the biology of the phlebotomid sandfly, Phlebotomus papatasi (Scopoli, 1786), Southeast Asian J Trop Med Public Health 24(1993), no. 3, 536–539.

[23] A. Stauch, R. Sarkar, A. Picado, B. Ostyn, S. Sundar, S. Rijal, M. Boelaert, …, H.P. Duerr, Visceral leishmaniasis in the Indian subcontinent: modelling epidemiology and control, PLoS Neglected Tropical Diseases

5(2011), no. 11, 1–12. doi: 10.1371/ journal. pntd. 0001405

[24] R. Topno, V. Das, A. Ranjan, K. Pandey, D. Singh, N. Kuman, N. Siddiqui, …, P. Das, Asymptomatic infection with visceral leishmaniasis in a disease-endemic area in Bihar, India, The American Journal of Tropical Medicine and Hygiene 83(2010), no. 3, 502–506. doi: 10.4269/

ajtmh. 2010. 09-0345

[25] J. van Griensven, M. Boelaert, Combination therapy for visceral leishmaniasis, Lancet 377(2011), no. 9764, 443–444. doi: 10.1016/

S0140-6736(10)62237-4

[26] E.E. Zijlstra, E.G. Khalil, P.A. Kager, A.M. El-Hassan, Post-kala-azar dermal leishmaniasis in the Sudan: clinical presentation and differential diagnosis, British Journal of Dermatology 143(2000), no. 1, 136–143. doi:

10.1046/ j.1365-2133. 2000. 03603. x