Mathematical modelling of the use of insecticide-treated nets for elimination of visceral leishmaniasis in Bihar, India

Visceral leishmaniasis (VL) is a deadly neglected tropical disease caused by a parasite Leishmania donovani and spread by female sand flies Phlebotomus argentipes. There is conflicting evidence regarding the role of insecticide-treated nets (ITNs) on the prevention of VL. Numerous studies demonstrated the effectiveness of ITNs. However, KalaNet, a large trial in Nepal and India did not support those findings. The purpose of this paper is to gain insight into the situation by mathematical modelling. We expand a mathematical model of VL transmission based on the KalaNet trial and incorporate the use of ITNs explicitly into the model. One of the major contributions of this work is that we calibrate the model based on the available epidemiological data, generally independent of the KalaNet trial. We validate the model on data collected during the KalaNet trial. We conclude that in order to eliminate VL, the ITN usage would have to stay above 96%. This is higher than the 91% ITNs use at the end of the trial which may explain why the trial did not show a positive effect from ITNs. At the same time, our model indicates that asymptomatic individuals play a crucial role in VL transmission.

Visceral leishmaniasis (VL) is a deadly neglected tropical disease caused by a parasite Leishmania donovani and spread by female sand flies Phlebotomus argentipes. There is conflicting evidence regarding the role of insecticide-treated nets (ITNs) on the prevention of VL. Numerous studies demonstrated the effectiveness of ITNs. However, KalaNet, a large trial in Nepal and India did not support those findings. The purpose of this paper is to gain insight into the situation by mathematical modelling. We expand a mathematical model of VL transmission based on the KalaNet trial and incorporate the use of ITNs explicitly into the model. One of the major contributions of this work is that we calibrate the model based on the available epidemiological data, generally independent of the KalaNet trial. We validate the model on data collected during the KalaNet trial. We conclude that in order to eliminate VL, the ITN usage would have to stay above 96%. This is higher than the 91% ITNs use at the end of the trial which may explain why the trial did not show a positive effect from ITNs. At the same time, our model indicates that asymptomatic individuals play a crucial role in VL transmission. poorest did not [42]. At the same time, even in relatively poor areas, the use of ITNs was already demonstrated to be very cost-effective [43,44].
A cluster-wide distribution of ITNs during the KalaNet trial reduced the vector density [45] but did not reduce the risk of L. donovani infection or clinical VL [46]. As a consequence, since 2010, ITNs are not part of the standard government VL control programme in India. Other trials in Bangladesh showed that ITNs may reduce the VL incidence rate [36,47]. The use of ITNs was also recommended for PKDL and VL-HIV patients [9]. This apparent contradiction raises the question about the role that ITNs may play in controlling VL [46]. The purpose of this paper is to gain insight into the situation by mathematical modelling.

Mathematical models of visceral leishmaniasis
There are many mathematical models of VL dynamics, see for example [48][49][50][51] for recent reviews. A model of VL transmission at the district-level of Bihar to estimate the basic reproduction number was created in [2]. Different vector control measures for VL elimination were considered in [52]. A multistate Markov model of VL was developed in [53]; an individual-based, stochastic, compartmental model of a temperature-driven sand fly population was developed in [54], which also simulated the effects of the use of drugs administered to cattle on the vector control. Chapman et al. [55] developed methods to analyse longitudinal spatial incidence data on VL and PKDL. A set of three age-structured model variants based on [56], each with individuals from a different disease stage being the main contributors to transmission: asymptomatic individuals, previously immune individuals in whom infection has reactivated, and individuals with PKDL, was created in [57]. The cost-effectiveness of different drug treatments was studied in [58,59].
A comprehensive model of VL for the Indian subcontinent to fit data collected from the KalaNet trial was developed in [56]. Their model extends the Susceptible-Infected-Recovered structure for the human population by segmenting it into five distinct stages according to an individual's infection status determined by the results of three diagnostic markers: (i) a polymerase chain reaction (PCR) test, the earliest infection marker able to pick up the presence of antigens [60], (ii) a direct agglutination test (DAT) [39] which measures the antibody response, and (iii) the leishmanin skin test (LST), also called Montenegro test [61] which detects the cellular immunity [62]. The model incorporates the role and significance of asymptomatic cases on VL transmission, which remains uncertain to this day [63]. The model also includes two lines of treatment of symptomatic VL cases, a possible treatment failure, relapse into PKDL and PKDL treatment. Most of the model parameters were found by fitting to data from the KalaNet trial using maximum-likelihood optimization method.

Game theory and disease prevention
From the behavioural perspective, disease prevention, such as ITN usage, produces public goods (such as herd protection and potentially the elimination-achieving zero cases-of the disease) that are nonrivalrous and non-exclusive [64]. When a sufficient proportion of the population is protected, then the slightest cost associated with using protection will outweigh the risk from infection [65]. Since the individuals often act in a way that maximizes their self-interests rather than the interests of the entire group [66,67], disease prevention is prone to free-riding. The 'free-riders', incorporated in our model as individuals that do not use ITNs, avoid the costs of prevention while they benefit from the preventive actions of others.
Residents in VL endemic areas seem to be reasonably aware of the role of ITNs in the prevention of VL and other vector-borne diseases [44]. We will assume that all individuals are provided with the same information such as VL incidence rates, treatment costs and ITN coverage, and use the information in the same and rational way to assess costs and risks [65]. In our setting, the game theoretical framework assumes that people will sleep under ITNs not to protect others, but to protect themselves. Specifically, people will not start sleeping under ITNs once they learn they have VL, but they may start sleeping under ITNs once they learn that someone in their community has VL or PKDL.

. Main objectives
In this paper, we will study an individual's choice to use ITNs to protect themselves against sand fly bites. We are interested to see why the KalaNet trial found that the ITNs were not effective in reducing VL incidence. We would also like to see if ITNs, possibly together with other control measures, could assist with VL elimination. This paper is organized as follows. We will adapt a compartmental ODE model of VL originally introduced in [56] and extend it by adding a compartment for untreated PKDL cases and by a parameter describing the ITN use. We then derive explicit formulae for the disease-free and endemic equilibria of the ODE system. We calibrate the model based on data found in literature. We then use maximum-likelihood method to estimate the probability of VL transmission from asymptomatic cases to sand flies. We validate our model by demonstrating it fits the KalaNet trial data. We show that asymptomatic individuals play a crucial role in VL transmission and without the asymptomatic transmissions, VL would be eliminated. We derive conditions for the VL elimination (no cases of VL), perform the game theory analysis and compare the levels of ITN use needed for the elimination with the Nash equilibria levels. We show that, even if ITNs are not provided for free, it is in individuals' self-interests to use ITNs more than the current coverage. Moreover, if all individuals used ITNs optimally, VL would be very close to being eliminated; in fact, the disease would become eliminated as a public health problem. We include detailed mathematical analysis as well as the Policy-Relevant Items for Reporting Models in Epidemiology of Neglected Tropical Diseases (PRIME-NTD) Summary Table [ [56]. The model is shown in figure 1.   Figure 1. Scheme of the ODE model for VL transmission; based on [56]. Majority of the human population is in the asymptomatic VL cycle: individuals are born as susceptible (S) and if bitten by an infected sand fly (I F )-signified by the red dotted arrow-they progress through early asymptomatic (I P ), late asymptomatic (I D ), recovering asymptomatic (R D ) and recovered (R C ). Recovered individuals lose immunity and eventually become susceptible again. Sand flies are born susceptible (S F ). If they bite an individual in any of the red or brown compartments-signified by the red curly brace-they become exposed (E F ) and eventually infectious (I F ). All human compartments are colour-coded based on their PCR, DAT and LST tests as shown on the left; informally, the intensity of the colour corresponds to the intensity of the infection which increases from left to right. The solid black arrows represent transitions between compartments.

S I
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 Humans are born susceptible (S) at rate Λ. They test negative in all PCR, DAT and LST tests. At some point, they may be bitten and infected by an infectious sand fly. For simplicity, we assume that ITNs offer a perfect protection and thus the force of infection is where p is the percentage of the population that uses ITNs, β −1 is the duration of the sand fly's feeding cycle, i F is the probability a human gets infected by a bite of an infectious sand fly, and I F is the number of infected sand flies. The infection starts by an early asymptomatic stage (I P ). Individuals in I P are potentially infectious to sand flies, but do not exhibit symptoms. They are PCR-positive, DAT-negative and LST-negative. The stage lasts for a period g À1 P and patients then progress to a late asymptomatic stage (I D ). They are still potentially infectious to sand flies, and PCR-positive, DAT-positive and LST-negative. The stage lasts for a period g À1 D . The vast majority (a fraction of f D ) recovers without ever showing any symptoms. At first, the recovering individuals are PCR-negative, DAT-positive and LST-negative. We will denote those R D . Eventually, after time r À1 D , the recovering individuals become DAT-negative and LSTpositive. We will denote those by R C . Both R D and R C individuals are immune to re-infection. However, the individuals in R C lose immunity at rate r C and become susceptible once again.
Only a small fraction, f S , of late asymptomatic patients (in I D ) develop symptoms of an acute VL and progress to a compartment I S . Since the symptoms are persistent, individuals typically seek diagnosis and treatment. The time from the onset of the symptoms to receiving a proper diagnosis and the beginning of the treatment is denoted g À1 S . The first-line treatment (I T1 ) lasts for a time t À1 1 . However, with probability s 1 , the treatment fails and individuals have to go through a second line of treatment (I T2 ) which lasts for a time t À1 2 . The most common first-line drug used in Bihar, India is miltefosine and the most common second-line drug is amphotericin B [83].
The first and the second line of treatments are successful with probabilities s 3 and s 5 , respectively. The individuals recover and move to the R T stage. The R T individuals are immune to re-infection, PCRnegative, DAT-positive and LST-negative. After time r À1 T , they will move to R C . With probability s 2 (s 4 for the second line of the treatment), the treated individuals appear to be recovering, but their infection becomes dormant. These patients, denoted R L , are PCR-negative, DAT-positive, LSTnegative. They are not infectious to the sand flies, but unless they die of natural causes, they will relapse and develop PKDL in time r À1 L . A fraction p TL of PKDL cases seek a relatively long and expensive treatment that lasts a time t À1 3 . We will denote those patients by I TL . The remaining PKDL cases, denoted by I UL , remain untreated. They spontaneously recover at rate g L . The cases that recover move to R T .
All individuals in I D , I S , I T1 , I T2 , I TL and I UL are PCR-positive, DAT-positive and LST-negative. All individuals in I X compartments for X ∈ {P, D, S, T1, T2, TL, UL} are potentially infectious to sand flies. The probability that a fly gets the infection while biting such an individual is denoted by i X .
All individuals die naturally at rate μ. The symptomatic VL cases without treatment, (I S ), have a higher mortality and die at additional rate m K , i.e. at the rate m S ¼ m þ m K . The VL treatments are aggressive, and the drug toxicity may cause permanent, irreversible damage (although a newer drug miltefosine is far less toxic than previous SSG [84]). A fraction d T1 (or d T2 ) of the treated patients die from the first (or second) line of treatment. The mortality rates in I T1 and I T2 are thus The sand flies follow the standard Susceptible-Exposed-Infected dynamics. All sand flies are born susceptible (S F ) at rate of m F N F . They become exposed (E F ) after feeding on an infectious human (in I P , I D , I S , I T1 , I T2 , I TL , I UL ) at rate where β is the bite rate and i X is the probability that an individual in stage I X infects a sand fly. Note that, unlike for the force of infection λ in (2.1), we do not consider the factor (1 − p) because we assume that only individuals who do not use any protection can get infected (and become infectious to sand flies). The duration of the latent stage is denoted by s À1 F . After that period, the sand flies become infectious (I F ). Since the mortality rate of sand flies, m F , is relatively high, we do not consider any recovered stage.

Summary of the notation
The scheme of the VL dynamics is shown in figure 1.

5
The compartments are denoted S (susceptible), I (infected and potentially infectious), R (recovering) and E (exposed). Subscript F denotes sand flies. For the human compartments, we use the subscript P for PCR-positive, D for DAT-positive, C for cellular immunity (LST-positive), S for symptomatic, T for treated, U for untreated and L for PKDL. The notation is summarized in table 1. Parameters are summarized in tables 2-5. The parameter values are estimated in §4 based on empirical evidence and data from the literature. Greek letters stand for the rates. Lower case roman letters stand for probabilities/proportions of the populations.
2.3. Differences between our model and the original model in [56] Aside from slight differences in the notation, we made the following changes and additions to the original model from [56].
(i) We explicitly added a parameter p signifying the level of protection against sand fly's bites by using ITNs. (ii) We separated PKDL cases (I HL in [56]) into treated (I TL ) and untreated (I UL ) cases to better account for the different duration of the stages and infectivity of the cases.   i T1 probability that an individual in I T1 infects a feeding sand fly 0 [0,1] [24] i T2 probability that an individual in I T2 infects a feeding sand fly 0 [0,1] [24] i TL probability that an individual in I TL infects a feeding sand fly 0 [0,1] [24] i UL probability that an individual in I UL infects a feeding sand fly 0.  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 (iii) We assume that the human population size is constant with a birth rate Λ (as opposed to α H N H considered in [56]). This change makes the ODE system less sensitive to changes in the birth rate and death rate. In the original model of [56], even a small change in α H , μ or m K would destabilize the system and result in either exponential growth or exponential extinction of the population. In our model, a (reasonable) change in Λ, μ or m K will cause the system to converge to an equilibrium with a potentially different population size, but the population will not go extinct nor grow above any bound. We can also solve for the equilibria of the dynamics. (iv) We keep the exponential growth of sand flies, denoted a F N F in [56], but we stipulate that a F ¼ m F in order to keep the population size constant. (v) We consider the death rate of treated individuals to be Specifically, we do not include the death rate of untreated individuals (m K ) and we consider the rates different with different treatment.

Numerical implementation
We coded the model in Matlab (v. 2020a with optimization toolbox) and made the code available in the electronic supplementary material. Following the best practices highlighted in [102], we included several tests to ensure correctness of the code. Specifically, we checked that analytical and graphical solutions (which were coded independently) yield the same results. We closely tracked the code execution to make sure the code runs as expected and values are passed correctly from function to function. Throughout the work on this manuscript, the formal analysis and numerical code were developed side by side and checked against each other. Independent cross-checks on a graphical calculator were also performed.

Analysis
The detailed analysis is shown in appendix A. Here we present only the summary. There are two possible equilibria: (1) the disease-free equilibrium with all humans and sand flies susceptible in S 0 = Λ/μ and S 0 F ¼ n F (L=m), and (2) the endemic equilibrium given by The explicit formulae are given in appendix A.3. When the ITNs usage is p, the effective reproduction number, R e (p), is the average number of new infections caused by a single-infected individual [103,104]. The formula for R e (p) can be derived using the next-generation matrix method [105] but also directly as shown below.
Assume that there is a single-infected person in compartment I P . As derived above, that person spends an expected time T Comp in each of the compartments Comp [ {I P , I D , I S , I T1 , I T2 , I TL , I UL }; the  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 formulae for T Comp are given in (A 33)-(A 43). During the time T Comp , the infectious individuals expose sand flies at rate bi X N F ¼ bi X n F (L=m) where X ∈ {P, D, S, T1, T2, TL, UL}. Each of the exposed sand flies becomes infectious with the probability s F =(m F þ s F ). Each sand fly stays infectious for the time m À1 F and during that time it infects humans at rate (1 À p)bi F N ¼ (1 À p)bi F (L=m). Putting it all together yields where It follows from (A 59) that the endemic equilibrium exists only if R e (p) . 1. While we did not perform a formal stability analysis, the ODE system-although large-is quite standard and similar to one considered in [106]. Consequently, the disease-free equilibrium is globally asymptotically stable if R e (p) 1, and unstable if R e (p) . 1. The endemic equilibrium is locally asymptotically stable if R e (p) . 1.

Model calibration 4.1. Sand fly parameters
The number of sand flies per human is n F ¼ 3 with range 2.1-3.5 ( [7], electronic supplementary material, table S15). This is in line with the India KalaNet site: 938 sand flies (94.2% of which were P. argentipes) from 325 households [45]. We note that this is different to the estimated value n F ¼ 5:27 with the range of 3.45-9.90 used by Stauch et al. [56].
The duration of the feeding cycle is β −1 = 4 days as in [56] with the range 2-8 days [85]. Sand flies normally take one blood meal per oviposition cycle [86,107,108] and the cycle usually takes 4 days [85] although infectious sand flies seem to feed more often [109,110]. We note that other models such as [3,111,112] used β −1 = 10 days based on [86], which estimated β −1 between 6 and 30 days based on experimental results of [107], which seemed to let the flies P. ariasi oviposit for six or more days. While the literature varies on the actual length, it all agrees on the fact that the flies bite only once per oviposition cycle.
The sojourn time in the latent stage E F is s À1 F ¼ 5 days, with the range 3-7 days [87]. The life expectancy of sand flies is m À1 F ¼ 14 days with the range 6-31 days [88].

Human parameters
The birth rate in rural Bihar, India is 27.7 births per year per 1000 people [89]. For the simulation purposes, we will assume the range of Λ to be [0.0105/12, 0.0333/12] per month. The life expectancy in Bihar is about 67.8 years [90]. For simulation purposes, we will use the range 50-90 years.
If the symptomatic VL patients remain untreated, their life expectancy, m À1 K is between 2 and 3 years [2]. We note that we did not find any data confirming this estimate. In [56], the authors assumed m À1 K ¼ 5 months. In ( [4], p. 14, table 6), it is reported that a significant number of patients die within a month from the diagnosis; yet it is not quite clear how long it took to be diagnosed. The report itself states that the mortality is associated with many factors, including late health care seeking. Consequently, while we did not find data supporting [2], we will use m À1 K ¼ 30 months with the range 5-36 months.
The sojourn time of early asymptomatic stage, g À1 P , is between four and six months [91]. This value is based on unpublished laboratory experiments on L. donovani-infected grivet monkeys (Cercopithecus aethiops) in which the seroconversion (DAT-to DAT+) takes place four to six months after infection [91]. The value g À1 P ¼ 5 months differs significantly from 60 days used in [56]. However, since [56] derived their values by fitting data to their model, we opted to use the value directly measured in the experiments. Also, the models considered in [113] use 150-202 days. Finally, a multi-state Markov model of VL developed by Chapman et al. [53] predicts the duration to be 147 days (95% CI 130-166).
The sojourn time of the late asymptomatic stage, g À1 D , is about one month. In [91], the authors were able to determine the time from the seroconversion to the development of active VL in 28 incidents of VL during their study. While they do not provide specific times, they report that most patients (23, i.e. 82.0%) royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 developed active VL in two to three months after seroconversion, while three (11.0%) took less than two months and the remaining two (7.0%) took more than three months. This means that within three months, all but two patients developed the symptoms, i.e. g D ¼ À ln (2=28)=3, i.e. g À1 D % 1:13 months (figure 2). We will consider the range to be 0.5-4 months.
The time between onset of symptoms to the start of treatment, g À1 S , is one month, with the range 0.5-2.5 months [92]. According to [92], there are 30 days (with the range 17-73 days) between the onset of symptoms and the diagnosis, and then 1 day (with the range 0.5-3 days) between the diagnosis and the start of treatment. Also, Sundar et al. [100] reports the median diagnosis time (symptoms to diagnosis) as five weeks with the range three to eight weeks. Finally, WHO [4] reports the median delay from onset of fever till diagnosis and then diagnosis to treatment is 59 days. In [56], g À1 S is 1 day because, as a condition of the KalaNet trial, it took the patients 1 day to get treated once diagnosed (as in [92]). However, the proper meaning of g S is the rate at which individuals leave the compartment I S (conditional on not dying). Since it takes a long time to get diagnosed [92,100], g À1 S is more in the order of a month than in the order of a day.
We estimate that the immunity in R C lasts for r À1 C ¼ 33 months, ranging from 10 to 38 months. An unexpected loss of LST reactivity in 150 (30%) of 499 cases over 1 year and in 127 (51%) of 246 cases over 2 years was observed in [62]. This indicates an exponential decay at rate r C ¼ À ln ((499 À 150)=499) ¼ 0:3575 per year. At this rate, there should be 120 = 246 exp(−0.3575 × 2) of LST-positive cases after 2 years which is in close agreement with observed data from [62]. We also used Matlab to fit an exponential decay function to data from [62]. It resulted in r À1 C ¼ 0:36176 À1 years with the range from 0.321 −1 to 0:4025 À1 years, see figure 2. We note that [91] also observed a loss of LST positivity in the range of l = 15−40% within several six-month-long periods corresponding to r À1 C ¼ À6= ln (1 À l) between about 12 and 37 months. They attributed such a great variability to the variability of VL incidence. Data from [62] seem more reliable and more suited for the estimate as they track the same individuals for a longer period. Finally, our estimates are significantly more than 307 days that [56] used to fit the KalaNet trial data. However, as seen later in the validation, our model with these parameters still fits the KalaNet reasonably well, and since this parameter value fits data from [62], we will use r À1 C ¼ 33 months. We will estimate r À1 D by six months based on data from [93]. This means that the DAT-positivity is vanishing at the rate À ln (50=375) ¼ 2:015 per year, i.e. r À1 D ¼ 5:96 % 6 months. We did not find data for the range, but estimate it between four and eight months.
Since there is no significant difference between compartments R T and R D , we assume, as in [56], that r À1 T ¼ r À1 D . We will assume that the baseline value of p, the fraction of population using ITNs, is 0.7 as in [38]. However, as p is variable, we will be interested in what values of p yield the elimination of VL.
We set f S , the fraction of I D that develop symptoms to be 0.035 with the range from 0.01 to 0.15. There is a wide range of values found in the literature for this parameter. For example, 0.0142 in [113], 0.013 in [114] and 0.147 in [53]. Out of 375 DAT-positive individuals in Bihar and Nepal followed by Ostyn et al. [93], seven developed an acute case of VL, meaning that f S ¼ 7=375 ¼ 0:019. Rahman et al. [21] surveyed N = 22 699 individuals for KA and PKDL in the past 6 years and determined that 813 respondents had  [91] fitted to y = 28 e −0.8845t which yields r À1 D ¼ 1:13 months. (b) Data from [62] fitted to y = 1.0008 e −0.36176t which yields an estimate for r À1 C as 0:36176 À1 years, i.e. 33 months. (c) Data from [22] fitted to y = 1.0441 e −0.21625t which yields g À1 L ¼ 0:21625 À1 years, i.e. 55.5 months.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 KA and eight had PKDL with no history of KA. From these findings, we calculate f S as follows. If we disregard any symptomatic infections and assume that individuals only undergo asymptomatic cycles S ! I P ! I D ! R D ! R C ! S, we will get that g P I P ¼ r C R C (the influx to and the outflow from a compartment is the same within and across all compartments). Since g P I P is the (monthly) incidence rate of late asymptomatic cases, we know that in t years there would be 12tr C R C late asymptomatic cases. While [21] does not provide the value of R C , we can use R C =N ¼ 0:35 for the prevalence of LST-positive cases from [62]. This gives us the estimate that over the period of 6 years, the cohort of Because [21,22] used a much larger sample size than [93] also notes that incidences were larger in India than in Nepal, we will adopt f S ¼ 0:035 and consider a range from 0.01 to 0.15. We estimate the fraction of asymptomatic patients who become dormant to be f L ¼ 0:00055 with the range [0.0001, 0.001]. This will be done in a similar fashion as the estimate for f S based on data from [21,22] We note that authors of [56] provide an estimate f L ¼ 0:0001 based on [21] although it seems that they wanted to use f L ¼ 8=22 699 ¼ 0:00035.

Estimates for visceral leishmaniasis transmission probabilities
Unlike [56], we will assume that individuals receiving treatment cannot infect the sand flies, i.e. i T1 ¼ 0, i T2 ¼ 0 and i TL ¼ 0. The treated individuals are not in their regular environments but rather at the treatment facility. Consequently, they are probably better protected from sand fly bites. Moreover, even if they are bitten and infect the fly, they are too far from their home for this event to meaningfully contribute to VL transmission in their home/village. This assumption is in agreement with [24], which reported that after receiving treatment, the VL patients are probably less infectious to the sand flies resulting in less transmission to family members during their 18-month follow-up. Similar findings were presented in [94].
To estimate the remaining transmission probabilities from humans to sand flies, we use the recent findings from [94], see also [115]. In [94], the authors found that 42 (54.5%) or 60 (77.9%) of 77 patients with active visceral leishmaniasis (as assessed by microscopy or quantitative PCR, respectively) and 11 (42.3%) or 23 (88.4%) of 26 patients with active PKDL transmitted parasites to sand flies. The results for PKDL transmission are slightly higher than those observed in [95] where only 27 out of 47 PKDL patients transmitted the disease. The difference can probably be attributed to the facts that [94] used more sand flies per patient (30-35 females and 10-12 males versus 20-25 females and 5-10 males) and for a longer time (30 min versus 15 min) than [95]. While the number 88.4% is large and it may seem that the parasite transmission from untreated PKDL patient to the vector is almost certain, it does not automatically mean that every bite results in the transmission. Unfortunately, neither [94] nor [95] report exact numbers of infected sand flies per patient and categorize patients only to those that did not transmit the parasite at all or those that transmitted it to at least one fly. Based on ( [94], electronic supplementary material, figure 1B), the average number of infected sand flies was around 1.5 while the average number of dissected blood-fed sand flies was around 15 per patient. This would give i S % 0:1. For the lack of better data, we adopt i UL % 0:1 as royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 well. We demonstrate in §5 that the exact values of i S and i UL have only a small effect on the overall results of the model.
The role of asymptomatic individuals in VL transmission was also studied in [94]. They found that none of 187 asymptomatic VL patients transmitted the parasite to a sand fly. However, based on ( [94], table 2), the vast majority of their patients tested DAT-positive and PCR-negative, i.e. they are categorized as recovering asymptomatic (R D ) by Stauch et al. [56] and our model, and as such already assumed to not be infectious. Despite an unfortunate difference in terminology, there is thus no factual difference between results of [94] and the assumptions of our model. None of the patients studied in [94] were PCR-positive and DAT-negative (early asymptomatic, I P ). Only 11 patients were DAT-positive and PCR-positive (late asymptomatic, I D ). Given the results presented in §5, this sample size may still be too small to draw conclusions about the roles of these early and late asymptomatic individuals in VL transmission. Consequently, in §5 we estimated the values of i P and i D differently by fitting the model to KalaNet data.

Treatment parameters
We will assume that the first line of treatment fails with the probability s 1 = 0.06 [96] and that the range is [0.01, 0.12]. Miltefosine is a common and effective first-line treatment in Bihar, India, with a clinical efficacy of 94% [96]. Studies like [116] give s 1 = 0.03 and lower, while other studies such as [117] show that the treatment failures by common treatments can be as high as s 1 = 0.12.
The fraction of KA patients who appear to recover under KA first-line treatment but will develop PKDL (conditional on surviving treatment, 1 À d T1 ) is s 2 = 0.063 [117]. This is based on the presentation 'Cohort observational study to estimate the cumulative incidence of PKDL in VL patients treated with three regimens in Bihar (Presenter: Suman Rijal)' that reported on a cohort of 1622 KA patients treated between 2012 and 2015 and followed up in 2016 and 2017 to determine the occurrence of PKDL. The cumulative incidence of PKDL was 6.3%, with a PKDL rate of 4.8%, 5.7% and 9.2% depending on the treatment. Since Stauch et al. [56] reported s 2 = 0.03, we will use [0.02, 0.1] for the range.
We will assume the probability to die during the treatment (as a result of the treatment) to be about d T1 ¼ d T2 ¼ 0:04. This is based on ( [99], table 6) which reported the death incidence of ( presumably treated) KA cases as 45.75 per 1000 (compared with a regular mortality rate of 2.9 per 1000 in the state of Bihar for the age group 15-59). The mortality varies greatly by age group from about 0.02 (5-14 years) to about 0.13 (60+ years). We will thus assume the range to be [0.02, 0.13]. We note that [56] assumed f T = 0.05. Overall, the numbers seem in line with other studies; Ahasan et al. [118] reported 27 deaths out of 553 patients (4.8%) treated with sodium antimony gluconate and, cumulatively, nine studies discussed in [118] report 51 deaths out of 1819 patients (2.8%).
The time from the apparent cure from KA until relapse to PKDL is assumed to be r À1 L ¼ 21 months, with the range 16-26 months [21]. This is in line with 19 months reported by Islam et al. [22].
The duration of the first line of treatment is t À1 1 ¼ 1 month, with the range [0.6, 1] months [96]. The duration varies by the treatment type, but the most common first-line drug used in Bihar, India is miltefosine [83].
Similarly, the duration of the second line of treatment is t À1 2 ¼ 1 month, with the range [0.6, 1] months [96]. The most common second-line drug used in Bihar, India is amphotericin B [83].
The duration of the PKDL treatment is t À1 3 ¼ 6 months with range one to seven months. There are three different kinds of treatments: sodium stibogluconate (SSG), miltefosine (MF) and amphotericin B. The SSG treatment takes six months (six 30-days-long cycles [21]), the MF treatment takes approximately three to four months [97,98], and amphotericin B treatment takes three weeks [22]. While amphotericin B appears to be the fastest treatment, it does not guarantee that PKDL will be cured within that time period. In fact Islam et al. [22] followed 30 PKDL cases treated with amphotericin B, and while 27 reported an improvement in four months, only one case was completely cured. We do not have additional data but it seems that even in the case of this treatment, the duration of PKDL is six months.
The duration of untreated PKDL, g À1 L , will be estimated as follows. In [22], the authors followed 98 PKDL patients that never received treatment and provides estimated resolution rates: 8% within 1 year of onset, 34% within 2 years and 67% within 5 years. Fitting the exponential decay to these data by Matlab (see figure 2) yields g À1 L ¼ 55:5 months with the range of 32 months to 16 years. The fraction of PKDL patients that receive treatment is estimated as p TL ¼ 0:5. This estimate is based on [22], which reports that out of 185 PKDL cases, 98 did not seek the treatment.

The cost parameters
Once infected, all individuals have to go through I P and I D . Those stages are asymptomatic and consequently with no associated costs. The costs appear only if and when an individual experiences symptoms (in I S ), is being treated (I T1 , I T2 and I TL ) or is recovering after the treatment (in R T and R L ). The costs are denoted by C I S , C IT1 , C IT2 , C ITL , C RT and C RL and the actual values are estimated below.
The Indian government provides free care to VL patients; however, hoping for a quick cure of seemingly minor illness, patients prefer to access private providers, which contributes to high out-ofpocket expenditures [119]. Even though VL is often misdiagnosed in private hospitals, patients continue to access private care until they can no longer afford it [119].
The cost of being in compartment I S is estimated as C I S ¼ 19 USD [100], and we will assume the range [11,26]. The average monthly income in Bihar, India in 2010 was 38 USD per month [100]. It takes about a month to get a diagnosis [92]. During that time, the patient loses 2.14 weeks (about half a month) of work, i.e. the cost is about 19 USD.
The cost of getting the first line of treatment is C IT1 ¼ 146 USD [100], and we will assume the range [120,170]. This includes 127 USD of direct medical costs for one month of treatment [100], and the indirect cost of lost work. For every month of illness, half a month of work is lost (19 USD loss) [100]. This is in line with 100 USD direct medical costs used by Hasker et al. [83].
The cost of getting the second line of treatment is C IT2 ¼ C IT1 ¼ 146 USD [100]. Both lines of treatment last approximately one month, so the cost of the individual second treatment will be about the same.
The cost of getting treatment for PKDL is assumed to be C ITL ¼ 349:00 USD [101]. We will assume the range [290, 410]. The direct cost of SSG treatment for PKDL costs 179 USD [101]. The treatment lasts about six months and so the loss of productivity during this illness is 170 USD.
The cost of recovering from treatment without a dormant infection is C RT ¼ 57 USD [100], and we will assume the range [60,110]. During recovery, a patient can only work 3.21 weeks/month on average, instead of the normal 4.29, i.e. losing 25% of the wage per month [100]. Since the average monthly wage is 38 USD and R T lasts approximately r À1 T ¼ 6 months, the total loss is about 0:25 Á 38 Á 6 ¼ 57 USD. The cost of recovering from treatment, but with a dormant infection is C RL ¼ C RT ¼ 57 USD [100]. We assume that even though individuals in this stage do have a dormant infection, the time and the cost needed to recover from treatment will be about the same.
The overall cost of VL infection is thus given by C VL ¼ P(I P ! I S )C I S þ P(I P ! I T1 )C IT1 þ P(I P ! I T2 )C IT2 þ P(I P ! I TL )C ITL þ P(I P ! R T jI T1 or I T2 or I TL )C RT þ P(I P ! R L jI T1 or I T2 )C RL , (4:4) where the probabilities P(I P ! C) (and P(I P ! CjC 0 )) of getting to a compartment C (through a compartment C 0 ) when currently at a compartment I P are given by P(I P ! I S ) ¼ T ID f S g D , ( 4 :5) P(I P ! I T1 ) ¼ T I S g S , ( 4 :6) P(I P ! I T2 ) ¼ T IT1 s 1 t 1 , ( 4 :7) P(I P ! I TL ) ¼ T RL p TL r L , ( 4 :8) P(I P ! R T jI T1 or I T2 or I TL ) ¼ T IT1 s 3 t 1 þ T IT2 s 5 t 2 þ T ITL t 3 (4:9) and P(I P ! R L jI T1 or I T2 ) ¼ T IT1 s 2 t 1 þ T IT2 s 4 t 2 : (4:10) The cost of ITN is C ITN ¼ 3:62 USD [36].

The role of asymptomatic individuals in visceral leishmaniasis transmissions
The transmission probabilities i P , i D for asymptomatic individuals, and i F for the transmission from sand flies to humans were estimated using the built-in Matlab function to fit the model to KalaNet trial data ( [56], table 1) as described in appendix B. The probability that a susceptible human becomes infected by an infected sand fly was estimated as i F % 1 (figure 3). This high transmission probability is consistent with the fact that the parasite can modify the sand fly's feeding apparatus so that the fly feeds more persistently and releases more parasites [109,120].
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 The estimates for i P and i D depend on p, the ITN coverage during KalaNet trial as well as the transmission probabilities from untreated symptomatic individuals, i S and i UL . However, figure 3 demonstrates that, for reasonable values of p, i S and i UL , the estimates for i P and i UL are fairly stable.
For small p ≈ 0, we get i D % 0:02 and i P % 0:004 which is consistent with [56]. However, for p ≈ 0.7, we get i D % 0:05 and i P % 0:01. If p grows even more, the values of i D and i P grow as well. Note that the estimates for i P and i D are indeed increasing in p-to achieve the same disease prevalence as measured in KalaNet trial when the population uses higher level of protection p, the disease must be more transmittable.
At the same time, as the probability of transmission from an untreated symptomatic individuals increases, both i P and i D decrease, i.e. the asymptomatic individuals become less important if the symptomatic individuals are more likely to transmit the disease. Nevertheless, even if symptomatic individuals transmit the disease with 100% probability, the role of late asymptomatic individuals is still not negligible.
When i S and i UL are around 0.1 as recently measured by Singh et al. [94], we get i P ¼ 0:01 and i D ¼ 0:05. This means that late asymptomatic individuals (PCR-positive, DAT-positive and LST-negative) are roughly 50% as important to VL transmission as untreated symptomatic VL and PKDL cases.

Model validation
The model is validated using KalaNet trial data ( [56], table 1). Our model gives the following prevalences (the KalaNet data are in the parenthesis): 0.7599 (0.76) for S þ R C , 0.0979 (0.1) for I P , 0.0221 (0.02) for i D , 0.1173 (0.12) for R D and 0.0108 (0.005) for I F .
To validate the model on another dataset, one would have to potentially update the parameter values to properly reflect the time and location of the experiment where the data came from. Then, we can use formulae for endemic equilibrium from appendix A.3 to obtain distribution of population across different model compartments.
We note that while not impossible, it is hard to make KalaNet trial data, our model and the model from [56]

Minimal insecticide-treated net coverage needed for visceral leishmaniasis elimination
To find the ITN usage level necessary to achieve the complete elimination of VL, we need to find the smallest p 0 ∈ [0, 1] such that when p ≥ p 0 , then R e (p) 1. It follows from (3.2) that R e (p) is the effective reproduction number when nobody is using the protection. Consequently, , if R e (0) . 1, 0, otherwise:

<
: For the parameter values as specified in tables 2-4, p 0 = 0.95963, i.e. one needs just under 96% ITN coverage for a complete VL elimination.

Optimal voluntary use of insecticide-treated nets
In this section, we will find the optimal proportion of the use of ITNs. We are looking for Nash equilibrium, i.e. a proportion that, when adopted by the population, no individual has an incentive to deviate from their choice. To find p NE , a Nash equilibrium value of p, we have to solve where C ITN is the cost protection, bi F I Ã F =(m þ bi F I Ã F ) is the probability of getting infected by an infected sand fly, and C VL is the expected cost one pays after such an event. Note that (5.3) is an equation for p because I Ã F is a function of p. Figure 5 illustrates a graphical solution of (5.3). After algebraic manipulations shown in appendix C, we get , ( 5 :4) where I Ã F is given by (C5). It follows that p NE % p 0 . For our parameter values, we have p 0 = 0.95963, p NE ¼ 0:95956 with R e ( p NE ) ¼ 1:0018. It means that the disease can be almost eliminated by the optimal voluntary use of ITNs alone ( provided broken ITNs are replenished and people get perfect information about incidence rates, ITNs usage and various costs). In fact, with the voluntary use of ITNs, the disease would become eliminated as a public health problem.

Sensitivity analysis
The sensitivity of the outcomes (ITN use for disease elimination, p 0 , and optimal ITN use, p NE ) on different parameter values is displayed in table 6 and shown in figures 8 and 9. Since p NE % p 0 , only the sensitivity of p 0 is shown. It follows that p 0 (and p NE ) are not overly sensitive to any parameter. The sensitivity index is at most 0.0846 (for μ or Λ) or −0.0849 (for β −1 ). It is the second highest for m À1 F and closely followed by n F and i F . Sensitivity to other parameters is between −0.011 and 0.02.
It follows that increasing the time between the bites and the reduction of the lifespan or the number of sand flies are the most promising control measure apart from using ITNs.

Discussion
In this section, we will mostly discuss the model limitations and possible extensions.
The model of sand fly dynamics could be improved in at least four possible ways. First, one could incorporate a seasonal dynamics by allowing for a variable birth rate. With this change, the analysis would become more involved as we could not focus on equilibria any more. Second, it is known that infectious sand flies tend to feed more often than non-infectious sand flies [109,110]. One may thus need to consider two different bite rates: one for susceptible sand flies (this will play a role in l F , i.e. humans infecting sand flies) and another one for sand flies infecting humans (this will play a role in λ, i.e. sand flies infecting humans).
Third, we did not consider any ITN-related mortality of the sand flies. This is in agreement with [122], although slight reduction of sand fly prevalence was observed by Picado et al. [45]. Fourth, we also do not consider increased mortality of infected sand flies [110,123]. To incorporate any of these last two improvements, one would have to consider m F ¼ m F (p) for an appropriate increasing function m F (p). Furthermore, the total mortality of the infected sand flies would be m F (p) þ m I (p) where m I (p) represents the additional mortality caused by the infection and the fact that infected flies bite more often, i.e. could suffer from an increased ITN-related mortality as well. To make these additions work, we would have to assume the birth rate to be L F ; otherwise, the sand fly population would not have a non-trivial finite equilibrium.
We also did not consider that ITNs offer only imperfect protection. First, smaller mesh size and insecticide treatment increase the net efficacy, but a small number of sand flies were still found inside ITNs [124]. Furthermore, as noted in [46], recent entomological findings in India indicate that L. donovani vectors are more exophilic and exophagic than previously reported [18,122]. If P. argentipes bite people outdoors (e.g. in the early evening when and where ITNs are not deployed), ITNs will have a limited impact on L. donovani transmission. To account for these, one would have to introduce a parameter e, the entomological efficacy of ITNs against the vector. The factor (1 − p) in (2.1) for the force of infection in humans would then change to (1 − pe) = (1 − p) + p(1 − e). Here 1 − p corresponds to an individual not using an ITN and p(1 − e) corresponds to the fact that ITN was used but did not protect (because the fly got inside anyway, or the bite occurred outside). Further changes would have to be made in (2.2) in the calculation of the force of infection for the sand flies. We would have to properly account for two distinct scenarios depending whether the fly bites an infected individual royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 who does not use ITNs, or whether it bites an individual despite them using ITNs. The first scenario accounts for (1 − p)/(1 − pe) cases, the second one for the remaining p(1 − e)/(1 − pe) cases. Moreover, in the second scenario, we would need to add one more factor (1 − e) to stand for a second 'failure' of the ITN.
Recent results of [55] suggest that incorporating spatial structure into the model would greatly increase its realism. Table 6. The sensitivity index SI y of a variable y on a parameter x was calculated as (x/y) · (∂y/∂x) [121]. The numbers were rounded to the three decimal places. Parameters are as specified in tables 2-4. The sensitivity index −0.5 means that a 1% increase of a parameter value x will result in the 0.5% decrease of the variable y. Also, there are several limitations of the game-theoretical framework, particularly the assumption that all individuals are provided with the same information and use it in the same (and rational) way to assess costs and risks [65]. First, ITN ownership is associated with wealth, and the cost of ITNs is often prohibitive for the poor [42]. Moreover, VL susceptibility is bound to reduce with wealth, as malnutrition and certain interior wall types are VL risk factors [125,126]. One would have to account for this by considering a non-homogeneous population. Second, VL has a long incubation period during which individuals are asymptomatic yet may be contagious. This makes the risk assessment prone to errors as people may think that there is no or very low risk of contracting VL (and thus stop using ITNs) while in reality someone in their community may already be asymptomatic. The risk assessment is further complicated by the fact that even professionals often misdiagnose acute symptomatic VL cases [119] and by the social and cultural stigma associated with VL [20].
For simplicity, we focused our analysis on equilibria of the ODE system. This approach had several limitations and disadvantages. As already mentioned above, it precludes us from incorporating seasonal sand fly dynamics. Similarly, the model cannot properly capture dynamical changes in pricing or availability of ITNs or drugs as was recently happening due to COVID-19 [127]. Most importantly, to understand and model the final stages of VL elimination, one should consider not only the equilibria but also how long it takes to reach an equilibrium. One also needs to explicitly model the interactions of the potential dynamical ITNs coverage and the disease dynamics.
Unlike vaccines that offer a long-lasting protection after just one decision, ITNs have to be used repeatedly. In principle, individuals decide every night whether or not to use ITNs. This dynamic decision process is not captured in our model. However, we hypothesize that coupled with long incubation periods of VL and especially PKDL, the seemingly optimal behaviour of using ITNs only if acute VL and PKDL cases reach above certain threshold, could lead to periodic spikes and drops of VL and PKDL cases. This is illustrated in figure 6. We can see that soon after the ITN coverage is relaxed, the early asymptomatic, I P , and symptomatic VL cases are on the rise, followed by a rise of asymptomatic recovered cases, R D , and only a very small and slow rise of royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 late asymptomatic cases, I D . Once the population starts to use more ITNs again, I P declines almost immediately, soon followed by a decline of I D and symptomatic VL cases and a trailing decline of asymptomatic recovered cases. At the same time, even as ITN coverage is on the rise, there is a slow but steady increase of PKDL cases that act as a reservoir of infections and can restart the cycle again. Overall, the importance of understanding of the consequences of dynamical decisions warrant future studies.

Conclusion
We expanded the mathematical model of VL transmission [56] to understand the role of ITN usage in VL dynamics and possible elimination. We calibrated the model based on the data found in the literature. We validated the model on KalaNet data. We concluded that in order to completely eliminate VL, the ITN usage would have to be about 96% or more. At such an ITN coverage, the effective reproduction number is less than 1 and the disease-free equilibrium is stable. Our work seems to be the first explicitly giving the ITN coverage threshold for VL elimination. Our work, specifically figure 5, may also explain why there was no significant effect of ITN use on VL incidence during the KalaNet trial-the risk of VL infection (which correlates with the mean cost function) is essentially a step function with a jump just around the disease elimination threshold. At the end of the KalaNet trial, the use of ITNs was very high: 91% of the individuals slept more than 80% of the nights under treated nets [46]. Yet, even with such high ITN use, our model predicts that the disease is not eliminated and the risk of infection almost does not change. The difference between 91% and 96% does not seem large, but note that 91% coverage leaves 9% of the population unprotected. This is more than double the unprotected population if the ITNs use is 96% or more.
From the policy-making perspective, while our model demonstrates that the elimination of VL by ITN use is possible in theory, it is unlikely in practice; the ITN coverage of 96% required for complete VL elimination is unrealistically high. At the same time, as seen from figure 6, a decrease of ITNs coverage can lead to sharp and dramatic increases of VL cases. Instead of abandoning the use of ITNs completely after the unfortunate findings of KalaNet trial, we recommend that the ITN use should be combined with other intervention methods. Figure 6 also carries a warning that regardless of what method is used to bring down VL incidence, the methods should stay in place long after VL seems eliminated-otherwise the disease can quickly resurface back.
Our model, specifically the sensitivity analysis illustrated in figure 8, suggests that increasing β −1 , the interval between two sand fly bites, from 4 days to 6 or more days can reduce the ITN coverage needed for VL elimination from 96 to 90% or even lower (which falls below 91% observed during the KalaNet trial). The increase of β −1 could be achieved by giving the sand flies other biting opportunities. While controversial, this intervention method has been studied for other vector-borne neglected tropical diseases such as Chagas disease [129]. The impact of the presence of cattle on VL incidence has already been studied separately. The cattle are associated with increased VL risk in some studies and decreased risk in others [38,130], reflecting the complexity of the effects on sand fly abundance, aggregation, feeding behaviour and VL incidence [131]. An increase in blood meal opportunities could lead to increase in the sand fly population size, possibly erasing the positive effect of increased biting interval. Consequently, an effective elimination strategy should include multiple interventions at the same time-(i) using ITNs to prevent insect bites as much as possible, (ii) using cattle and/or other animals (that cannot serve as hosts for L. donovani) to divert vector bites away from humans [132], and (iii) using insecticide residual spraying [7], or (iv) destroying breeding sites [52] to keep the vector population under control.
Another important finding of our model is that asymptomatic individuals, specifically those that test PCR-positive and DAT-positive, play a crucial role in VL transmission. Our model predicts that those individuals are about 50% as infectious as the untreated symptomatic individuals; although we admit that such an estimate seems about twice as high as assumptions made in other studies, such as [55,56]. Based on our model, without the asymptomatic transmissions, VL would already be eliminated. The apparent difference between our results and recent experimental results of [94] stems only from a different use of the word asymptomatic. Majority of cases from [94] tested DATpositive but PCR negative, i.e. were considered as asymptomatic recovered and assumed noninfectious in our model. There is thus no factual difference between results of [94] and the assumptions of our model. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 We believe that empirical research needs to be done to understand the role of asymptomatic patients in VL transmissions-specifically those that tested DAT-positive and PCR-positive and were infected roughly four to six months in the past. As shown in figure 6 as well as in the sample of cases studied in [94], these individuals may be a relatively elusive part of the population; yet they may be an important missing piece in understanding VL dynamics.
Data accessibility. The Matlab code supporting this article has been uploaded as part of the electronic supplementary material.
Authors' contributions. A.K.F., J.A.W.: conceptualization, formal analysis, investigation and writing: original draft. C.P.G: conceptualization, formal analysis, investigation, software and writing: original draft. Y.L. supervision, writing: original draft. J.R. writing: original draft, writing: review and editing, software, methodology, supervision, conceptualization, resources. D.T. writing: original draft, writing: review and editing, methodology, supervision, conceptualization, formal analysis, validation, funding acquisition. The order of the authors was determined as follows: The first three authors are undergraduate students in alphabetical order. The last three authors are professors in alphabetical order.
Competing interests. There are no competing interests. Funding. A.K.F., C.P.G. and J.A.W. were supported by the VCU REU program in mathematics funded by the National Security Agency (https://www.nsa.gov/what-we-do/research/math-sciences-program/proposal-guidelines/) grant no. H98230-20-1-0011 awarded to D.T. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. The publication fees were covered by the Virginia Tech's Open Access Subvention Fund.

A.1. Equilibria of the ODE system
The equilibria of the dynamics (A 1)-(A 15) are found by solving the following system of algebraic equations: ) The system (A 16)-(A 30) will be solved as follows, corresponding to starting at compartment S and going downstream. By (A 17), and since l ¼ (1 À p)bi F I F , Solving for S Ã yields : By (2.2), Thus, : From (A30), From (A29), royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 Thus, and : Once I Ã F is calculated from (A 59), we can calculate E Ã F , S Ã F and S Ã as and : Finally, using above, we can also find all the human compartments as below and

Appendix B. Numerical estimates of transmission probabilities
The transmission probabilities i P , i D for asymptomatic individuals, and i F for sand flies were estimated using the built-in function fmincon from Matlab optimization toolbox. We minimized the difference between the model predicted prevalences and the prevalences measured during KalaNet trial ( [56], royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 table 1). Depending on the initial guess for the optimum, the function fmincon sometimes returns a local optimum rather than global optimum. To make sure we seeded properly and indeed found global optima, we initialized the search with 1000 random seeds and collected all results. We collected all the local minima and searched for the global optima from those. The results are shown in figure 7. It follows that the global optima found by this complicated procedure, i P ¼ 0:0110, i D ¼ 0:0487 and i F ¼ 1 is close enough to local optima i P ¼ 0:0111, i D ¼ 0:0481 and i F ¼ 1 found by initiating the searches for i P and i D at 0 and for i F at 1 (figures 8 and 9),

Appendix C. Finding Nash equilibrium
To find a Nash equilibrium value of p, we have to solve (5.3), i.e. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 Recall that, by (A 59), and in particular I Ã F is decreasing in p and thus bi F I Ã F =(m þ bi F I Ã F ) is decreasing in p. The maximum value of I Ã F is thus (1 À (1=R e (0))) bi F T Cycle (N Ã F =R e (0)) þ (m F =s F ) þ 1 (C 3) attained for p = 0, and consequently (5.3) has a solution only if It follows from (5.3) that (when (C 4) is true) : Thus, after algebraic manipulations of (C 2), we get that p NE is given by where I Ã F is given by (C 5).    royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201960 Table 7. Policy-Relevant Items for Reporting Models in Epidemiology of Neglected Tropical Diseases (PRIME-NTD) Summary Table as specified in [82].
principle What has been done to satisfy the principle?
Where in the manuscript is this described?

stakeholder engagement
We did not directly engage the stakeholder. However, we based our model on [56]  We detailed procedures of how we extracted parameter values that were not found directly. §4.

communicating uncertainty
We performed sensitivity analysis based on [121]. §5.5, table 6, figures 8 and 9 and the Matlab code.

Testable model outcomes
We predicted that asymptomatic individuals (PCR+, DAT+) can transmit VL and their infectiousness is at about 50% of the level of untreated symptomatic individuals. §5.1 and 5.2.
We gave formula for the minimal ITN coverage needed to achieve complete elimination of VL and predicted that 96% ITN coverage is needed. §5.3.
We also predicted that voluntary use of ITN should yield elimination of VL as public concern. §5.4 with mathematical details in appendix C.