Spatio-temporal modelling of Leishmania infantum infection among domestic dogs: a simulation study and sensitivity analysis applied to rural Brazil

Background The parasite Leishmania infantum causes zoonotic visceral leishmaniasis (VL), a potentially fatal vector-borne disease of canids and humans. Zoonotic VL poses a significant risk to public health, with regions of Latin America being particularly afflicted by the disease. Leishmania infantum parasites are transmitted between hosts during blood-feeding by infected female phlebotomine sand flies. With a principal reservoir host of L. infantum being domestic dogs, limiting prevalence in this reservoir may result in a reduced risk of infection for the human population. To this end, a primary focus of research efforts has been to understand disease transmission dynamics among dogs. One way this can be achieved is through the use of mathematical models. Methods We have developed a stochastic, spatial, individual-based mechanistic model of L. infantum transmission in domestic dogs. The model framework was applied to a rural Brazilian village setting with parameter values informed by fieldwork and laboratory data. To ensure household and sand fly populations were realistic, we statistically fitted distributions for these entities to existing survey data. To identify the model parameters of highest importance, we performed a stochastic parameter sensitivity analysis of the prevalence of infection among dogs to the model parameters. Results We computed parametric distributions for the number of humans and animals per household and a non-parametric temporal profile for sand fly abundance. The stochastic parameter sensitivity analysis determined prevalence of L. infantum infection in dogs to be most strongly affected by the sand fly associated parameters and the proportion of immigrant dogs already infected with L. infantum parasites. Conclusions Establishing the model parameters with the highest sensitivity of average L. infantum infection prevalence in dogs to their variation helps motivate future data collection efforts focusing on these elements. Moreover, the proposed mechanistic modelling framework provides a foundation that can be expanded to explore spatial patterns of zoonotic VL in humans and to assess spatially targeted interventions. Electronic supplementary material The online version of this article (10.1186/s13071-019-3430-y) contains supplementary material, which is available to authorized users.


Household-level host distributions
The number of each type of host at each household was assigned in each model run by sampling from distributions of host numbers per household. We obtained these distributions by fitting to survey data from the Marajó region collected in July and August of 2010 at 140 households across seven villages; via a questionnaire, data were collected on the number of adults, adolescents, and children resident in the home, as well as the number of dogs and chickens kept at the home [1].
We fit a Poisson distribution to the data for each host; we also fit a negative binomial distribution when the sample variance was not less than the sample mean. The Poisson distribution was initially chosen and fitted because these data are count data; the negative binomial distribution is also a discrete distribution and is considered an alternative to the Poisson distribution for data where the Poisson assumption of the mean being equal to the variance is not appropriate. Distributions were fitted using maximum likelihood estimation via the poissfit and fitdist functions from the Matlab R version R2016b Statistics and Machine Learning Toolbox. We compared the fitted Poisson and negative binomial distributions using the Akaike information criterion (AIC) [2]. The AIC estimates the quality of each fitted distribution relative to the other. Specifically, the difference in AIC values between the Poisson and negative binomial models (∆AIC) indicates the level of support for the Poisson model fit relative to the negative binomial model fit. As a general rule, a larger difference indicates that the Poisson model is less plausible compared to the negative binomial model. Notably, a ∆AIC≤ 2 corresponds to the Poisson model having substantial support relative to the negative binomial model, while a ∆AIC> 10 indicates no support for the Poisson model [3].
The results of fitting the Poisson and negative binomial distributions (Supplementary Figure 1) determined that where both distributions were fit the negative binomial distribution was preferred as all ∆AIC were greater than 2 (Supplementary Table 1). Visual inspection provided additional support for this outcome (Supplementary Figure 1). Figure 1: Distributions of the number of hosts per household. Data (bars), best fit Poisson distributions (blue solid line) and negative binomial distributions (red dashed line, considered when the sample variance was greater than the sample mean) fitted using maximum likelihood estimation for the number of adults and adolescents, children, dogs, and chickens resident at households in Marajó.

Sandfly abundance
Sand fly trapping data from villages in Marajó were used to obtain realistic estimates of the abundance of sand flies, L h , at households. For each household h, L h comprised of two parts: a constant initial estimate K h 1 1−ζ and a seasonal scaling v(t).
Data on the abundance of female sand flies, specifically the vector species Lutzomyia longipalpis, were available from a previous study of 180 households in fifteen villages on Marajó island where sand fly numbers were surveyed using CDC light-traps [4]. The trap-count abundance, K h , was sampled from these data. With traps only capturing a proportion of female sand flies expected at households, a proportion of the female sand fly population, ζ, remained unobserved. Accounting for this inconsistency necessitated the scaling of the trap-count abundance by a factor of 1 1−ζ .
Sand fly populations have been observed to exhibit temporal dependencies. To incorporate seasonality into the model, at the beginning of each time step we applied a time-dependent scaling factor, v(t), to all initial abundance estimates. To produce the scaling factor v(t), a smooth trend line was fitted, via a Lowess smoother, to the mean number of female Lutzomyia longipalpis trapped over an eight month period across eight different households in the village of Boa Vista, Marajó [5]. The curve was extrapolated over the remaining four months of the year for which no data were available. This highlighted a peak in sand fly abundance during January, at the transition from the dry to wet season (Supplementary Figure  2). Expected vector abundance then dropped and attained its minimum level in May and June, coinciding with the end of the wet season. Normalising the curve between 0 and 1 gave our seasonal scaling factor v(t). Similar temporal patterns were observed in the data split by the eight households (Supplementary Figure 3) and split by location within household (Supplementary Figure 4).
Amalgamating the initial abundance estimate and seasonal scaling components results in the following seasonally-scaled sand fly abundance at household h at time t, Supplementary

Proportion of infectious sand flies
The proportion of infectious vectors at household h was comprised of a time-independent background level of prevalence, φ, constant across all households plus an additional proportion dependent on the number of infectious dogs in the neighbourhood of household h. The contribution from each type of infectious dog (high and low infectiousness) was computed separately. We matched the radius r defining this neighbourhood with the maximum sand fly travel distance (taken as 300m at the baseline [6], see Table 1 in main text).
Initially, we computed the proportion of biomass of infectious dogs of type x within radius r of household h. This proportion of biomass was subsequently modified by a linear weighting function to account for a reduction in impact with increasing distance from household h: where H h (r) is defined as the set of households within distance r of household h, d(h, k) is the distance between households h and k, and N h,s is the number of dogs of infectiousness type s at household h.
Using these weighted biomass computations for infectious dogs, the proportion of sand flies that were infectious at household h on day t was computed as: where φ is the constant background level of prevalence, and m high and m low are upper bounds on the proportion of infectious sand flies obtained when the only hosts present were high infectiousness dogs or low infectiousness dogs respectively. Under an assumption that 80% of transmission from dogs to sand flies is caused by high infectiousness dogs, with the remaining 20% of total transmission events contributed by infected dogs with low infectiousness [7], the explicit calculations for m low and m high were as follows: withπ low andπ high denoting the proportion of infectious dogs that have low and high infectiousness, respectively, and m avg corresponding to the proportion of infectious sand flies obtained when the only hosts present are infectious dogs, obtained by averaging over both high and low infectiousness dogs.

Sensitivity coefficients
In a stochastic modelling framework, like the one developed in this study, outputs do not take a unique value. Instead, they take a range of values with a given probability, defined by a probability density function f . Therefore, to calculate a stochastic sensitivity coefficient for each parameter we followed the procedure outlined in Damiani et al. [8]. In brief, this technique evaluates the sensitivity coefficient Υ u p of the output variable of interest u with respect to each parameter p, Υ u p = Ωp Ωu ∂f (u(p)) ∂p f (u(p)) du dp, where Ωu is the domain of integration of u. Due to the computational demands of evaluating the density function for the entire parameter space, the integrals in Equation 1 were calculated on a finite domain.
The probability density function f (u(p)) and the partial derivatives ∂f (u(p)) ∂p were estimated using nonparametric kernel methods using simulation outputs from the model.