Complexity and critical thresholds in the dynamics of visceral leishmaniasis

We study a general multi-host model of visceral leishmaniasis including both humans and animals, and where host and vector characteristics are captured via host competence along with vector biting preference. Additionally, the model accounts for spatial heterogeneity in human population and heterogeneity in biting behaviour of sandflies. We then use parameters for visceral leishmaniasis in the Indian subcontinent as an example and demonstrate that the model exhibits backward bifurcation, i.e. it has a human infection and a sandfly population threshold, characterized by a bi-stable region. These thresholds shift as a function of host competence, host population size, vector feeding preference, spatial heterogeneity, biting heterogeneity and control efforts. In particular, if control is applied through human treatment a new and lower human infection threshold is created, making elimination difficult to achieve, before eventually the human infection threshold no longer exists, making it impossible to control the disease by only reducing the infection levels below a certain threshold. A better strategy would be to reduce the human infection below a certain threshold potentially by early diagnosis, control animal population levels and keep the vector population under check. Spatial heterogeneity in human populations lowers the overall thresholds as a result of weak migration between patches.

We study a general multi-host model of visceral leishmaniasis including both humans and animals, and where host and vector characteristics are captured via host competence along with vector biting preference. Additionally, the model accounts for spatial heterogeneity in human population and heterogeneity in biting behaviour of sandflies. We then use parameters for visceral leishmaniasis in the Indian subcontinent as an example and demonstrate that the model exhibits backward bifurcation, i.e. it has a human infection and a sandfly population threshold, characterized by a bi-stable region. These thresholds shift as a function of host competence, host population size, vector feeding preference, spatial heterogeneity, biting heterogeneity and control efforts. In particular, if control is applied through human treatment a new and lower human infection threshold is created, making elimination difficult to achieve, before eventually the human infection threshold no longer exists, making it impossible to control the disease by only reducing the infection levels below a certain threshold. A better strategy would be to reduce the human infection below a certain threshold potentially by early diagnosis, control animal population levels and keep the vector population under check. Spatial heterogeneity in human populations lowers the overall thresholds as a result of weak migration between patches. disease-free and endemic states, i.e. their ability to persist under perturbation in epidemiological and ecological conditions, to the changing threshold conditions. In the next section, we describe the hierarchical multi-host model, followed by calculation of system thresholds, their resilience and sensitivity to parameter variations.

Material and methods
2.1. The model 2.1. 1

. Model outline
We assume that the VL transmission cycle consists of humans, two non-human reservoirs and sandflies. This assumption encapsulates the generality of a multi-reservoir-human-sandfly transmission because it captures the effects of competing reservoirs in terms of biting preferences, host competences, lifespan and relative population levels, which is otherwise not possible when only one reservoir is present in addition to humans. Further justification for considering two reservoirs follows from the fact that goats and cows are the most probable reservoirs living near human dwellings especially in India [31]. Additionally, the movement of humans and sandflies between cities/villages/towns (called patches) can also affect the transmission dynamics [37].
We model the transmission dynamics by assuming that humans, sandflies and other non-human reservoirs live on n-patches and only humans can move between different patches. Sandfly movement is ignored as their flight range is limited (approx. 300 m) [37]. Although human movement happens on two levels (i) short-term commutes (affecting transmissibility based on visitation to other patches) where their patch identity is preserved [38] and (ii) long-term movement which includes migration to other patches (affecting population dynamics) [38][39][40][41], here we assume only short-term movement between patches. The susceptible humans (S Hi ) of patch i get bitten by infectious sandfly vectors (I vj ) while visiting a patch j, remain asympotomatic (E Hi ) for a period of time before becoming infectious I Hi ). Untreated infectious humans die because of infection or recover due to immunity, or get hospitalized (T Hi ) and receive treatment. Treated humans either recover or get PKDL (K Hi ) due to treatment failure. Humans with PKDL recover at a certain rate. Recovered humans are classified into (R Hi ). Recovered humans can lose immunity and join the susceptible class again. The two non-human reservoir populations are susceptible (S Ai , S Bi ) and become infectious (I Ai , I Bi ) by bites of infectious sandflies (ignoring the asymptomatic compartments due to lack of sufficient data) in patch i, lose immunity and join their respective susceptible classes again (SIS). The susceptible sandfly vectors in patch i (S vi ) bite an infectious/asymptomatic/PKDL human (either resident or visiting patch i) or an infectious member of non-human reservoir, become latent (E vi ) and finally infectious (I vi ). The resulting sandfly-human-non-human interaction network and its corresponding flow chart is shown in figure 1. The equations describing the full n-patch spatial model are as follows: Humans   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200904 In equation (2.5) index l indicates human (denoted by H ), and non-human reservoirs (denoted by A and B), indices i and j are used for space discretization indicating different patches, Δ ij is a notation used for the elements of matrix ℙ used to denote residence times of non-human host in different patches (for example, Δ ij = 0 indicates that non-human hosts do not move between different patches, p ij are the elements of residence times matrix ℙ capturing the strength of connectivity between patches i and j due to human movement). The symbols L B j, m Bj are used to represent recruitment and the natural death rate, respectively, of human, non-human hosts and sandfly vectors in the patch j. The symbol b is used to represent transmission probability in humans, non-humans and vectors. The symbols k ▪ represent patch specific within host biting heterogeneity (more on this in subsections 2.1.2 and 2.1.3). Additional human parameters in patch j are: 1/σ Hj is the latency period, τ Hj is the natural recovery rate, d Hj is the disease-induced death rate, γ Hj is the rate of seeking treatment/hospitalization in the infectious stage, r Tj is the rate of recovery/discharge following treatment/hospitalization, η j is the fraction of treated exhibiting PKDL, r Kj is the rate of recovery from PKDL stage. The parameter δ lj is the loss of immunity for human and non-human hosts in patch j, 1/σ vj is the latency period of sandfly vectors in patch j, N lj is the population of human and non-human hosts in patch j. The symbols C ▪ represent host competences, and α ▪ represent vector biting preferences for each host.
In an isolated patch, biting heterogeneity in sandfly vectors on hosts may arise due to a limited flight range within the patch. Uniform homogeneous biting is recovered in equation (2.5) (also see electronic supplementary material, Document) for large values of k j , the biting heterogeneity parameter in patch j. However, as k j gets closer to zero, the degree of heterogeneity in biting increases [42,43]. On the other hand, sandfly vectors may have different feeding preferences for humans and reservoirs A and B. The 'vector feeding index' α lj assesses the proportion of blood meals from host-l in relation to the proportional abundance of host-l in the host community [27], in general, it could be a complicated function of host population proportions [44]; however, here we take them as constant parameters. Additionally how competent are humans and reservoirs in transmitting the infection to susceptible sandflies is governed by host competence parameter, C lj , defined as the ability of host-l to successfully transmit the virus/parasite to the biting vector [45]. For simplicity, we assume host competency and vector biting preference do not change from patch to patch, i.e. α lj = α l and C lj = C l for all patches j. Note, in our model we allow for hosts to be dead-end (C l = 0) as well as competitive hosts (0 < C l ≤ 1). A detailed list of parameters and their interpretation is given in 1. different types of hosts (e.g. human, vector and animal reservoirs), 2. distinct reservoirs based on host competence, 3. spatial population distribution and structure, and 4. sandfly differential biting preferences.
Apart from human and sandfly vectors, we have considered two types of host reservoirs potentially representing animals. The efficiency of reservoirs to successfully transfer pathogens during a bite may vary depending on host competence. We assume that vector may have differential biting preference towards different host species and individuals in the population. The model captures discrete space via distinct population structure and movement patterns. Many mathematical models in the literature primarily assume homogeneous biting patterns between hosts. In our model, we consider differential biting preferences towards different hosts, hence, capturing heterogeneity in biting rates.

The force of infection (λ ▪ ) in the model
The force of infection in equation (2.5), is defined as the per capita rate at which susceptible hosts/vectors contract the infection [46]. The force of infection on humans, reservoirs and sandfly vectors is derived by combining the short-term movement (Lagrangian movement [38]), within host sandfly biting sandfly biting heterogeneity in patch j unitless variable p ij elements of residence times matrix for humans between patch i and j unitless variable D ij elements of residence times matrix for non-human reservoirs between patch i and j and reservoir B (l = B) unitless  [42,43], vector feeding indices {α lj } [27] and host competences {C lj }. Feeding preference assess the proportion of blood meals sandfly vectors have from host-l, whereas the biting heterogeneity k j indicates that all sandflies do not bite individuals in the host population equally and also that within similar host types all members may not receive the same number of sandfly bites. First, consider a modelling scenario with an isolated patch. Sandfly biting heterogeneity in a patch has two components: first, biting preferences for reservoirs and humans vary (because of sandfly feeding behaviour). Second, the members of the same host (humans or each of the two reservoirs) get varying number of bites (because of host individual's characteristics). Most models in the literature assume that the sandfly-host encounters are well mixed, i.e. any given host (human or non-human reservoir) receives the same number of sandfly bites. In general, this assumption is unrealistic, as a given sandfly may only have the opportunity to bite a limited number of hosts (e.g. some humans may receive more bites than others) in the population. Following the formalism developed for dengue [42,43,47], we assume that the number of VL transmission causing bites on mth host are Poisson distributed with mean θ m . If all host members receive the same bites on average, then θ m is the same for all the hosts in the population, implying the standard homogeneous transmission model. However, the biting means themselves may have a probability distribution over the population and under the assumption that the means θ m follow a Gamma distribution (as mean bites are greater than zero and can take a large range of values) with shape parameter k j and a rate parameter s c (scale parameter 1/s c ). Consequently, the marginal probability distribution of effective risky bites follows a negative binomial with mean k j =s c ¼ ðbb l a l = P l a l N lj ÞI vi (using composition of functions formula, Gamma distribution of Poisson distribution). With a few steps of calculation using the negative binomial [42,43,47], the risk of susceptibles becoming infectious is 1 À 1 þ ðbb l a l =k j P l a l N lj ÞI vi À Á k j (details are given in the electronic supplementary material, Document). Using the relationship between risk and rate [48] we obtain the rate of infection of susceptible hosts as Through a similar argument, the rate of infection for susceptible sandflies after biting infectious humans/reservoirs is In these expressions k j characterizes the level of heterogeneity of bites (susceptible humans/reservoirs by infectious sandfly vectors and infectious humans/reservoirs by susceptible vectors) [42,43]. Finally, incorporating the movement of humans between n-patches, equation (2.5) is obtained.
In this study, we study a one-patch model in relative detail and present preliminary results for a twopatch model separately. In the one-patch model, we increase its complexity as follows: first, switching off the parameters related to treatment and PKDL and other reservoirs we obtain a base model (BM; considered as our first model); second, including human treatment, we get the treatment model (BMT); third, with both treatment and PKDL parameters, the BMT-PKDL model is obtained. Finally, we include reservoirs in the transmission cycle. We will work our way upwards systematically adding complexity as described and later we will study a two-patch model with only humans and sandflies in the transmission cycle. In the subsequent sections, we will assume homogeneous biting and obtain system thresholds for simplicity of algebraic calculations. We come back to discussing the implications of heterogeneous biting on system thresholds in the Results and Discussion and conclusion sections.

Analysis
In this section, we assume homogeneous biting (k j is large so that the logarithmic term in equation (2.5) reduces to the standard mass action force of infection [46]) and obtain threshold conditions for the model.

Thresholds and multiple equilibria in one-patch model
In our model R 0 acts as one of the thresholds and it can be translated to a threshold biting or threshold sandfly population. This means that VL will invade the population only if R 0 > 1 [49]; however, if the disease-induced death rate is sufficiently strong then it can invade and persist even at a lower value of R 0 , i.e. R 0 ≥ R c ≤ 1, provided that the initial proportion of the infectious population is above the infection threshold. The infection threshold serves as an additional threshold and occurs due to a phenomenon known as BB. That this phenomenon is possible in VL is highlighted from the broad spectrum of parameter values estimated in previous studies on VL in the Indian subcontinent. In particular, one such study reported extra deaths caused by VL [14,20] which we find responsible to generate multiple equilibria and therefore multiple thresholds via BB.
The system of equations (2.1)-(2.4) exhibit two stable and one unstable infectious state under a restrictive parameter regime if the disease-induced death is sufficiently strong. The unstable infectious royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200904 states act as a threshold infection prevalence, above (below) which the VL invades (goes extinct in) the population. Another threshold arises out of a minimum sandfly population: below this sandfly number VL cannot persist in the population regardless of the initial prevalence of VL. Next we obtain expressions for these thresholds and discuss the regimes of their validity.
Assuming homogeneous biting, the model equations (2.1)-(2.4) have a disease-free equilibrium (DFE) given by The basic reproduction ratio R 0 evaluated at the DFE using the next generation matrix approach is The above expression for R 0 is true for the general one-patch model. The corresponding expression for the BM is recovered by setting A,B ¼ 0. In this VL system, a threshold R 0 is translatable to a sandfly population threshold. In general, the VL will invade the population only if R 0 > 1 [49]; however, if the disease-induced death rate d H is sufficiently strong then it can invade and persist even at a lower threshold, i.e. R 0 ≥ R c ≤ 1, provided that initial infection prevalence is above the infection threshold due to BB [50]. This invariably generates two thresholds in the system giving rise to a characteristic bi-stable region R c ≤ R 0 ≤ 1 so that the DFE coexists with a stable and an unstable endemic equilibrium: the value of R c gives the threshold vector population (N Sc ) while the unstable endemic equilibrium forms the human infection prevalence threshold (I Hc ).
The endemic equilibria, and hence the human infection prevalence threshold (I Hc ), are easily obtained in terms of the force of infection on humans λ H which satisfies a sextic polynomial equation The steps to obtain the coefficients {A i } are given in the electronic supplementary material, Document. The endemic equilibria , force of infection on reservoir A and reservoir B are related to the force of infection on humans as Humans In our model, the multiple thresholds are created (due to BB) if the disease-induced death rate parameter d H exceeds a certain value such that @l 1 =@R 0 j R0¼1 , 0 in equation (2.8) [46], assuming d A = d B = 0. Note that the condition for BB is never satisfied if d H = 0 as well. Given the multi-host nature of our model (and equation (2.8)) a closed-form relationship between the disease-induced death rate and other parameters for the existence of multiple thresholds to occur is cumbersome, therefore, we numerically confirmed that this criterion is satisfied for our choice of d H .

System resilience and integral stability (d IS )
We quantify the resilience of disease-free/endemic equilibria in the system using the integral stability index. Integral stability (IS) [33], 0 ≤ d IS ≤ 1, is a measure of the resilience of individual equilibria/ attractors in a multi-states-stable dynamical system. Resilience is defined as the capacity of a system to recover in the face of perturbations. Since in BB the endemic state coexists with a DFE, we characterize their resilience using the integral stability. For a detailed mathematical definition of IS see electronic supplementary material, Document.

Sensitivity of system thresholds
The system thresholds are sensitive to variations in parameters of the system. We use the formalism of partial rank correlation coefficient (PRCC) [51][52][53] to characterize the sensitivity of thresholds to variations in systems parameters. To calculate sensitivity to model parameters each parameter is assigned a Gaussian distribution with their respective mean as in table 1 and standard deviation to be 1/30 of the mean (so that neither of the thresholds are eliminated and therefore facilitating their sensitivity analysis). A more realistic assumption on the distribution of parameters will require the distribution of parameters to be calibrated against field data (which we do not attempt here).

Thresholds in two-patch model
To study the effect of spatial heterogeneity we consider the BM in two spatially isolated locations. Coupling between the patches is assumed to be only via human movement between them, making the simplifying assumption that sandflies do not move between patches. The threshold R 0S for the two-patch system is given by where R 0h21 , R 0h12 and R sqrt are described in the electronic supplementary material, Document, R 2 01 and R 2 02 are thresholds for isolated patches. The expression for R 0S shows that system can be above the global threshold even if threshold condition in isolated patch may not be satisfied. The threshold infection prevalence (discussed in the Results section) is obtained numerically using Newton-Raphson root finding method.

Results
In the one-patch BM with homogeneous biting (large k), i.e. no reservoir A and reservoir B (and no treatment in humans) the sextic polynomial, equation (2.8), reduces to the quadratic equation For a sufficiently large disease-induced death rate in humans [14,20] BB, with the characteristic two system thresholds, is observed: the first threshold is the human infection prevalence threshold (calculated using equations (2.10)-(2.12)) and the second threshold is the threshold vector population related to R c as obtained from the condition when A 2 1 À 4A 2 A 0 ¼ 0: here the function f (Λ l , μ l , d l , Λ v , μ v , …) on the right-hand side of equation (  To investigate the effects of biting preference, we consider the following cases: (i) when all hosts are equally competitive, i.e. C H = C A = C B = 1, and sandflies prefer to bite them equally α H = α A = α B = 1, (ii) all hosts are equally competitive, i.e. C H = C A = C B = 1 but α H = 1, α A = 0.6, α B = 0.1, (iii) cows are deadend hosts, i.e. C H = C A = 1, C B = 0 but sandflies prefer to bite humans and goats more as compared with cows α H = 1, α A = 0.6, α B = 0.1, and (iv) cows and goats are dead ends and biting preference is the same as in case (iii). In cases (i)-(iii) presence of reservoirs increases the bi-stable region, by first reducing the sandfly threshold (N Sc ) and increasing the human infection prevalence threshold, while in case (iv) the threshold (N Sc ) is larger and the human infection prevalence threshold is smaller relative to the BM, thereby reducing the bi-stable region, see figure 3a.
Another interesting feature that multiple hosts highlight is the changing resilience of the DFE/ endemic state as the reservoirs are added to the system. The integral stability (IS) d IS shows that if bistable region increases (decreases) then the resilience of the endemic state increases (decreases). The increased (decreased) resilience can be explained by the observation that when bi-stable region is large (small) a smaller (larger) number of sandflies are needed to stabilize the endemic state. In addition, a smaller number of initial human infections can push the system to the endemic state for a given number of sandflies (above the sandfly threshold for the BM). In the case where other reservoirs are dead-end hosts, these two factors combine to show an increase in resilience (d IS ) of DFE, as shown in figure 4.
Controlling the spread of infection requires looking for mechanisms to stabilize the DFE, and when BB is present the traditional way of ensuring R 0 < 1 for stabilization of DFE does not work, as control in such circumstances requires keeping track of the two thresholds: 1. the minimum vector population threshold N Sc 2. the threshold infection prevalence level in humans as defined by the unstable endemic equilibrium.
If either of these two thresholds are crossed from above, i.e. decreasing sandflies population to below N Sc and reduce human infection prevalence than the infection threshold, DFE is stabilized. When humans are subjected to treatment at a constant rate, both these thresholds shift such that the human infection threshold is reduced while increasing the threshold vector population N Sc for the disease to emerge. Humans with PKDL also transmit the disease resulting in an endemic state intermediate to BM and the model with a completely efficient treatment, as shown in figure 2a. Although having less efficient treatment is better than no treatment at all. This implies that, assuming the vector population does not increase, treatment can eliminate the disease; however, if the vector population is allowed to grow then the disease can remain endemic. Thereby the goal of reaching the human infection threshold (or the threshold prevalence) via treatment alone would become ineffective as it cascades to zero with increasing treatment rates, as shown in figure 5. Similar shifts in thresholds are observed when treatment of human hosts or vector control is applied in presence of other  Figure 4. Integral stability of the DFE and endemic equilibria (indicated by arrows) for the base model (black and red) and multi-host models (blue and green). In the multi-host model, dead-end reservoirs were considered to generate this diagram, demonstrating that the DFE is more resilient when dead-end hosts are introduced. Based on this diagram and figure 3, it is easy to see that the endemic state is more resilient when competitive reservoirs are considered.  Figure 5. The human infection prevalence threshold as a function treatment rate γ H (in the absence of reservoirs). It shows that a cascade of prevalence threshold values are created as the treatment rate is increased, before eventually the prevalence threshold drops to zero as shown in (a). This leaves behind only the sandfly population threshold as shown in (b), equivalently only the threshold R 0 = 1 (c). royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200904 reservoirs; however, the multi-host system relatively amplifies (dilutes) the effect of control strategies when dead-end (competitive) reservoirs are present, making VL control either easier or more difficult relative to when humans are the only hosts.
To study the effects of biting heterogeneity on system thresholds we take the BM, and allow an increasingly heterogeneous sandfly biting behaviour by changing the parameter k. As heterogeneity increases in the system the sandfly threshold increases, implying that large number of flies are required for infection to emerge; however, the relatively large sandfly threshold also lowers the threshold human infection prevalence. The behaviour of thresholds follow similar trend when host complexity is increased by adding more disease stages as discussed above for the homogeneous mixing case; however, in comparison the prevalence level is reduced. The thresholds shift in two ways. First, the sandfly threshold is increased and the infection threshold is larger as compared with the BM. Second, at the new sandfly threshold the infection threshold is lower in comparison with the BM. A bifurcation diagram indicating these trends is shown in figure 2b.
The thresholds obtained for the BM are sensitive to changes in parameters. We evaluate the sensitivity to model parameters corresponding to the BM such as sandfly biting rate on humans, PKDL and other parameters using PRCC. Detailed tornado plots showing PRCC for threshold sandfly population (N Sc ) and human infection prevalence threshold (I Hc ) are shown in figure 6a,b. The PRCC was obtained by assuming Gaussian distribution for varied parameters, giving rise to a distribution of threshold (N Sc , I Hc ) as shown in figure 6c-e. A positive/negative PRCC for a parameter means that the threshold increases/decreases with increasing/decreasing the parameter. The behaviour of PRCC shows that sandfly threshold (N Sc ) is strongly but negatively sensitive to biting rate (b), biting heterogeneity parameter (k), the factor of asymptomatic/latent moving to infectious class ( f ), transmission parameters from fly to human (vice versa) (β H , β v ), natural mortality rate of humans (μ H ) and proportion of infectious bites from PKDL cases (ρ 2 ). On the other hand it is strongly but positively sensitive to sandfly mortality (μ v ), human birth rate (Λ H ) and disease-induced death rate (d H ). Similarly, the human infection prevalence threshold is positively sensitive to the factor of asymptomatic/latent moving to infectious class ( f ) and loss of immunity (δ H ). The biting parameter (b) has a strong negative effect on this threshold as expected (if more bites, then infection can emerge at a lower infection prevalence). To control VL it is, therefore, desirable to enhance the parameters increasing the sandfly population threshold as well as increasing the human infection threshold for an optimal control of visceral leishmaniasis: as an example, PRCC indicates that the number of bites have to be decreased (so that a large number of sandflies are required for disease to emerge).
The movement of humans between neighbouring villages/towns can also affect the local thresholds. To study the effect of human movement between otherwise isolated locations we assume the BM in two patches. The coupling between the patches is due to short-term human movements, where humans stay for most time p ii in their home patch and spend a fraction p ij in the other patch, thus the fraction of time spent is the coupling strength (see equations (2.2)-(2.6) and their description). The spatial effect is revealed as soon as asymmetry in parameters of the models in two locations are introduced, as an example we chose the sandfly birth rate and human birth rate as asymmetric parameters and considered two different values of the coupling ( p ii , p ij ). Following this scheme, a bifurcation diagram demonstrating (2.3) effects of spatial coupling and parameter heterogeneity is shown in figure 7. In the first instance, we assume that birth rate of sandflies in patch one is same as that in a one-patch model discussed above, while in patch two it is 0.081 times that in patch one. We find that in presence of human movement, patch one exhibits higher sandfly threshold and higher infection prevalence threshold, whereas patch two has a lower sandfly threshold and lower infection prevalence threshold as compared with the isolated one-patch model. As coupling strength is increased the dynamics in the two patches is as follows: the sandfly threshold in patch one increases while infection threshold decreases, and in patch two the sandfly threshold decreases and so does the threshold infection prevalence (figure 7b). However, the effect of spatial coupling is to lower the sandfly threshold and infection prevalence threshold in the two patches when compared with an isolated patch. Similar effects are observed when asymmetry in human birth rate (and hence population) is also considered, as shown in figure 7c.

Discussion and conclusion
In this paper, we investigated the dynamics of visceral leishmaniasis via a dynamical system model that may exhibit multiple thresholds due to BB phenomenon in a certain parameter regime. BBs belong to a royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200904 class of disease systems where, in addition to the basic reproduction number, R 0 (equivalently the threshold levels of vector bites/vector population), initial size of infected individuals in hosts act as a threshold. Such systems have been observed in both micro-and macro-parasitic systems such as, tuberculosis (TB), dengue, lymphatic filariasis and onchocerciasis [1,[4][5][6]12,54,55]. Although the multithreshold systems have been explored both in the field [55] as well as in theoretical literature [6], the theme which it has demonstrated consistently is that control is not trivial in systems exhibiting multiple thresholds. Our current study investigates dynamical regimes where VL might exhibit multiple thresholds; thus, a potential cause for presenting challenges in control of VL in the Indian subcontinent. This model exhibits multiple thresholds due to BB phenomenon driven by nonlinear relationship between human mortality (disease-induced and natural mortality rates) and sandfly royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200904 transmission and demographic parameters (see electronic supplementary material, Document). Human infection prevalence spontaneously increases when the initial prevalence crosses a threshold level; however, there are two points to note: (i) new infections are generated by transmission from other population (sandfly population) even though disease-induced death in humans acts like a leakage in the system and (ii) increment in infection does not occur in human but in sandfly population, which then feeds back into human population. On the other hand, several reports have indicated the continued persistence of VL despite low death rates and no reported incidence in neighbouring regions [56,57], potentially as a result of higher number of initial infected arrivals in the region as compared with its neighbours. The initial infection prevalence has been recently understood to be a driving factor behind the emergence of chikungunya in South America [58,59], where despite international travel of people from affected regions for decades, the disease could not become endemic until a threshold of imported infections was reached in 2013. Importantly, the imported infection could not be accounted for in calculations of R 0 , which would have been less than one because the disease was not present in the region in the first place. This draws similarities with BBs vis-à-vis the threshold prevalence and highlights the importance of investigating initial infection prevalence among other factors. We began with a BM wherein only humans interact with sandflies without the influence of any control measures. Using the parameters from literature for India, where a sufficient amount of disease-induced death rate for humans was reported [14,20], we find that this BM exhibits BB. This presents a scenario where the system exhibits two thresholds: (i) human infection and (ii) sandfly  4)). (a) Bifurcation diagram is the same as base model in isolated patches when each patch has identical demographic parameters for humans and sandflies, regardless of coupling strength. (b) Bifurcation diagram for different coupling strengths (0.8 and 0.9 as indicated by arrows) when sandfly birth rate (and subsequently sandfly population) is higher in patch one than patch two. (c) Bifurcation diagrams in the two patches when human birth rate is higher in patch one than patch two. (d ) Bifurcation diagram for different coupling strengths (0.8 and 0.9 as indicated by arrows) in the two patches when human and sandfly birth rates are greater in patch one than patch two. These diagrams show that spatial heterogeneity in combination with parameter asymmetry can lower the thresholds as compared with an isolated system and the degree to which this happens depends on coupling strength.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200904 population density thresholds. Introducing treatment of infection in humans can help cross the human infection threshold, beyond which the disease is eventually eliminated; however, if sandfly population is left unchecked and continues to grow then human infection threshold presents itself as a cascading function of treatment, before eventually becoming zero, thereby making only a complete elimination of infection, as a rather impractical and only treatment goal. On the other hand, if the sandfly population remains the same, then the disease-free equilibrium can be stabilized when a sufficient rate of treating humans is employed. The scenario of increasing sandfly thresholds in conjunction with decreasing infection prevalence threshold can be equated to cascading thresholds [60], which is defined as the 'tendency of the crossing of one threshold to induce the crossing of other thresholds' [60]; here the changing human infection threshold also changes the sandfly population threshold with increasing treatment. The cascading thresholds are known to increase the resilience of one of the states [60]; here the resilience of the DFE is increased.
When more reservoirs, such as goats and cows, are available for sandflies to bite then scenarios for amplification (dilution) and subsequent increase (decrease) of the thresholds are possible; however, in all the possible scenarios, a shift in thresholds alters the region of stability of the endemic state. This also implies that multiple hosts either increase or decrease the resilience/persistence of DFE and endemic equilibria. For example, adding dead-end reservoirs to the BM, assuming equal biting preference for humans/reservoirs, lowers the equilibrium endemic state and in tandem the human infection threshold. Simultaneously, since some bites from sandflies are taken up by the dead-end hosts, a higher number of sandflies are required for the endemic state to emerge. Similarly, if the additional reservoirs are competitive for disease transmission, then a higher equilibrium endemic state-human infection prevalence threshold emerges, while requiring a lower number of sandflies. This follows from the fact that the competitive reservoirs act as a buffer for the disease, therefore, a smaller number of sandflies are required for endemic state to emerge. Additionally, the increasing and decreasing stability region of endemic state also follows if biting preferences are increased/decreased. Amplification/dilution of infection in different parameter regimes have also been studied in the context of host diversity although with density-dependent biting preference and fixed total population of all hosts [61]; this is the first time that an amplification and dilution phenomenon has been studied in the context of thresholds due to species diversity (here we assume population of individual reservoirs is a constant which is more realistic for visceral leishmaniasis).
Heterogeneity in sandfly biting behaviour implies that most bites come from a limited number of individuals. The interplay between sandfly population and host infection generates qualitative changes in respective thresholds: first, the force of infection is lowered, therefore, higher number of infectious hosts are required for sandflies to pick up sufficient infection, thereby increasing the infection invasion threshold for a given number of sandflies. Second, a new sandfly threshold is created which is higher than the corresponding homogeneous case. The higher sandfly threshold requires less infectious hosts for the endemic state to emerge.
The shifting thresholds indicate that control measures, such as treatment in humans and vector control used in combination with a control over reservoir population, easily stabilize the DFE. Thus, having multiple reservoirs can be beneficial if their population relative to humans can be managed given that the additional reservoirs are dead end. On the other hand, control will be much more difficult if the reservoirs are competitive and readily available for sandflies to bite. Spatial heterogeneity and connectivity via human movement change thresholds in non-trivial ways. Our simple model shows that infection can persist for a lower sandfly population and human infection prevalence threshold as compared with the geographically isolated case due to weak migration. Thus, control efforts focused to bring down infection cases to less than 1 in 10 000 at the block level (in Bihar, India) will have to take into account of the fact that human movement lowers the local thresholds. Sensitivity analysis using PRCC is able to identify parameters affecting the sandfly population threshold and human infection prevalence threshold. In particular, our results on PRCC show that sandfly threshold (N Sc ) decreases with an increase in biting rate and homogeneous mixing, the factor of asymptomatic/latent moving to infectious class, transmission parameters from sandfly to human (vice versa), natural mortality rate of humans, and proportion of infectious bites from PKDL cases (ρ 2 ). On the other hand, it is strongly but positively sensitive to sandfly mortality, human birth rate and disease-induced death rate. Similarly, the human infection prevalence threshold increases with the factor of asymptomatic/latent moving to infectious class, loss of immunity and heterogeneous biting pattern of sandflies, and it decreases with increasing disease-induced death rate, biting rate and treatment. Therefore, for control purposes it is desirable to reduce the contribution of PKDL cases by increasing the efficiency of treatment, decreasing the biting rates by covering the exposed skin and royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 200904 reducing disease-induced deaths in addition to vector control. Our results also show that an optimal combination of treatment and vector control will need to be evaluated as the control programmes move forward.
This work highlights that despite control efforts, importation of cases can keep the VL transmission going due to feedback coming from continued movement from high endemic regions to low endemic regions. This is particularly relevant in areas like Muzaffarpur in Bihar India, where recently spatial patterns in VL transmission were studied [56,57]. Our results are based on simplistic assumptions on host competence and vector biting behaviour in the transmission cycle when multiple hosts are present. However, host competence is governed by the biology of the animals involved, while biting preference in general could be a nonlinear function of reservoir availability/population [44]. In these later circumstances analysing the equations with the corresponding functional dependencies will be complicated, we propose that the amplification and dilution behaviour presented in this work will be relevant. Consideration of a general vector-feeding-function based on host availability could be a project for future investigations. Although our two-patch model considered here captures the intuitive idea that lower thresholds can arise due to spatial heterogeneity, investigating a more realistic model by combining short-term [38] and long-term migration will be the goal for our future investigation.