Understanding heterogeneities in mosquito-bite exposure and infection distributions for the elimination of lymphatic filariasis.

It is well known that individuals in the same community can be exposed to a highly variable number of mosquito bites. This heterogeneity in bite exposure has consequences for the control of vector-borne diseases because a few people may be contributing significantly to transmission. However, very few studies measure sources of heterogeneity in a way which is relevant to decision-making. We investigate the relationship between two classic measures of heterogeneity, spatial and individual, within the context of lymphatic filariasis, a parasitic mosquito-borne disease. Using infection and mosquito-bite data for five villages in Papua New Guinea, we measure biting characteristics to model what impact bed-nets have had on control of the disease. We combine this analysis with geospatial modelling to understand the spatial relationship between disease indicators and nightly mosquito bites. We found a weak association between biting and infection heterogeneity within villages. The introduction of bed-nets increased biting heterogeneity, but the reduction in mean biting more than compensated for this, by reducing prevalence closer to elimination thresholds. Nightly biting was explained by a spatial heterogeneity model, while parasite load was better explained by an individual heterogeneity model. Spatial and individual heterogeneity are qualitatively different with profoundly different policy implications.

Introduction

Heterogeneities in disease transmission play an important role in the epidemiology of vector-borne diseases and influence opportunities for control. The heterogeneous exposure to mosquito bites can drive vector-borne disease hotspots [1], and is a crucial factor in the optimal design of disease control intervention [2]. The degree of exposure heterogeneity can be as important as mean transmission rates in driving patterns of disease [3], but methods for measuring this heterogeneity vary and are rarely compared in the same setting [2,4,5]. It is also unclear how best to evaluate the heterogeneity of exposure within individuals in order to inform modelling and policy [6].

There are multiple levels of heterogeneity that contribute to the aggregation patterns of disease observed within a community: spatial heterogeneity, which is largely governed by ecological variation and environmental conditions [7], and individual heterogeneity, which is governed by many factors such as socioeconomic, behavioural and physiological variation among hosts [8,9]. In the context of heterogeneous exposure to mosquito bites, spatial heterogeneity may be due to landscape, rainfall, breeding site productivity, insecticide use, household size or urbanicity [1,10-13]. Individual attractiveness to mosquitoes will differ by sex, age, size and variability in human odours [14-16]. Spatial variation exists at multiple scales [9,13,17]—with different transmission dynamics between neighbouring villages and even between households [18]. These are rarely studied in the same place or through multiple measures. The transmission of lymphatic filariasis (LF), a mosquito-borne helminth infection, provides the opportunity to explore heterogeneity in bite exposure as well as in parasite burden. LF affects over 120 million people worldwide but is currently targeted for elimination. Our aim is to evaluate the multiple sources of heterogeneity which could undermine the LF elimination campaign.

Global efforts to eliminate LF through the mass distribution of anti-helminthic drugs have resulted in a large-scale reduction of prevalence [19]. However, there are numerous challenges to achieving LF elimination targets using community-wide treatments. Top-down, uniform strategies which aim for a specific intervention coverage or duration are unlikely to achieve elimination without appreciation for the significant heterogeneity driving transmission and extinction dynamics [20,21]. For LF, the target of less than 1% microfilaria (mf) prevalence set by WHO as a mark of success gives poor confidence in the probability of elimination [20]. While the recommended strategy may be sufficient in some areas [22], other areas can require many more rounds [23]. The true threshold prevalence below which transmission cannot be sustained depends on competence of the dominant vector, vector biting rates and microfilaria intensity [24]. Failure to break transmission would require community-wide mass drug administration (MDA) for the duration of the adult worm's lifespan, or direct testing and treatment, both of which may be prohibitively costly for a scaled-down LF programme. For successful elimination we require clear targets of the MDA coverage and duration needed to break transmission.

Vector control can increase the likelihood that an elimination campaign of recommended coverage and duration will achieve local elimination. The breakpoint prevalence of a vector-borne parasitic disease, below which transmission cannot be sustained, is dependent on vector biting density [25-26], so vector control will help to raise the threshold microfilaria prevalence. Supplementing MDA with vector control was recommended in countries where the burden is the heaviest [27], and evidence is mounting that vector control should be an essential component of the global elimination strategy [28]. In addition to reducing vector-borne disease transmission, vector-based interventions may also influence the spatial patterns of exposure and risk. For preventive chemotherapy vector-borne diseases such as LF, onchocerciasis and schistosomiasis, the success of community-wide coverage will be influenced by the degree of aggregation. For example, higher intervention coverage will be required in communities with highly aggregated bite risk to ensure appropriate coverage of hotspots [26]. If aggregation in biting differs significantly from village to village, a uniform strategy may underestimate the coverage required to break transmission across the implementation area (see table 1).

Table 1.

Policy consequences for different types of heterogeneity.

Heterogeneity	High	Low
Spatial	different villages may have drastically different prevalence, one cannot be compared with the other	use of sentinel sites can be justified; reduction in one village comparable to reduction in another with the same intervention
Individual	small group of individuals highly burdened and disproportionately contributing towards ongoing infection; targeted treatment may be necessary	low variation in individuals implies blanket coverage would be effective; no small subset of population driving disease implying systematic non-adherence less of an issue

Statistical models can be used in order to determine both spatial and individual heterogeneity within a count distribution [17] (figure 1). When individual heterogeneity is high, the count distribution is heavy-tailed and an individual's parasite count can be far from the mean. When the individual heterogeneity is low, the count distribution has a variance similar to the mean. When both spatial and individual heterogeneities are high, a highly over-dispersed distribution is produced with a greater than expected number of zeros observed compared with when the distribution is more spatially homogeneous. The result of the aggregation observed under high heterogeneity implies that an intervention that does not obtain good geographical coverage and population coverage may not be able to achieve targets in reduction or elimination. These differences in the type of heterogeneity have a profound impact on the control and elimination of parasitic disease; table 1 outlines the policy implications for each of these scenarios.

Figure 1.

Teasing apart different types of heterogeneity. Size of houses represents relative risk in space and size of people represents relative risk in individuals. A Gaussian process is used to simulate the mean rate (e.g. biting rate) across space, with both high (left-hand side) and low (right-hand side) variance. Compounding this is the variance around the mean at each spatial location, which is referred to as intrinsic heterogeneity. Example probability distributions with a mean of 10 and high and low heterogeneity are shown across the middle. Example outcomes for the four cases are given in the bottom row. How count data is aggregated and whether there is heterogeneity among individuals (individual) and/or among space leads to qualitatively different forms of count distributions. Policy implications for each of these situations are described in table 1. (Online version in colour.)

Our understanding of the sources of heterogeneity within a vector-borne disease transmission system is crucial for control and elimination because high heterogeneity is often associated with a higher basic reproduction number (R0) and a hard-to-reach threshold at which elimination can be achieved [21,29]. However, the effect of heterogeneity on disease prevalence will depend on numerous factors, including the transmission dynamics of the parasite or pathogen. Malaria parasites are cyclopropagative in the mosquito vector, while filarial worms are cyclodevelopmental with sexual reproduction occurring in the vertebrate host. Malaria transmission is highly efficient and one infective mosquito could successfully transmit malaria to multiple people. Filariasis transmission on the other hand is inefficient, requiring continuous high exposure to the infective stage larvae (up to 15 500 bites in one setting [30]) for a patent infection. Filariasis transmission models show that heterogeneous exposure results in a higher disease prevalence at low mean biting rates compared with homogeneous exposure, but this relationship changes at higher biting densities [26]. The threshold biting rate, leading to a non-zero endemic equilibrium, is significantly lower with heterogeneous biting [26]. In other words, heterogeneous exposure can sustain transmission at a comparatively lower prevalence, making it more difficult to break transmission with community-based interventions. Universal coverage of community-based interventions in a heterogeneous system may be inefficient, even leading to greater heterogeneity, and not protecting the high-risk households [31]. The spread of infection to the broader community from these households is a threat to elimination programmes and may require the integration of targeted interventions [32]. Properly implemented targeted control can result in impacts up to 4-fold higher than untargeted control [5,31].

Study aims

For elimination programmes to succeed, we must achieve the appropriate coverage, continuity and combination of interventions to break transmission and prevent resurgence. However, the approach and target coverage will depend on the aggregation of exposure and disease. It is therefore imperative to understand the impacts of heterogeneity on disease breakpoints to better tailor interventions and elimination campaigns. The aims of this study are twofold: (1) to determine what drives the heterogeneity in LF prevalence and intensity; and (2) to determine how aggregated biting patterns are influenced by vector control and the implications for LF elimination. The first aim compares the spatial relationships between breeding sites, anopheline biting rates, and infection prevalence and intensity to determine whether heterogeneity in disease status is driven by heterogeneity in either spatially dependent bite exposure or through individual variation. The second aim considers heterogeneity on a village scale by quantifying spatial aggregation of mosquito biting in five neighbouring villages before and after bed-net distribution. The fitted village biting heterogeneities are then used to parametrize an individual-based transmission model to estimate the implications of bed-net introduction on the sustainability of ongoing transmission. More broadly, this study aims to evaluate which of our standard measures and analyses of heterogeneity are most appropriate to evaluate heterogeneities which are relevant for infectious-disease control.

Methods

In order to understand the causes and effects of heterogeneity on the prevalence of LF and its underlying intensity we consider two approaches to analyse the heterogeneity of risk and infection. The first approach considers the non-spatial heterogeneity in bites and mf count by fitting these distributions by village to an over-dispersed distribution and measuring the amount of overdispersion for each fit. These fitted distributions for biting density are then applied to an individual-based model of LF transmission in order to understand how vector control impacts the ongoing transmission of LF.

The second approach considers how these indicators vary spatially and what the spatial association is between them in order to understand whether the heterogeneity in disease status or intensity is driven by spatial heterogeneity, individual heterogeneity or both. This was done by first fitting a model of individual and spatial variation to each disease outcome (mf prevalence, mf intensity and antigenic prevalence) in turn. A combined model was then used where the spatial variation is dependent on the biting density.

Study sites

Five villages in the East Sepik province of Papua New Guinea have been the focus of extensive research into filariasis epidemiology and transmission [20,33-34]. These villages received annual MDA from 1993 through 1998, with no further interventions until long-lasting insecticidal nets (LLINs) were distributed in August 2009. Self-reported LLIN use ranged 75–90% [35].

Infection prevalence

Antigen prevalence and microfilaria prevalence were measured in these communities in 2008 as part of the post-MDA evaluation [35]. This was done by BinaxNow filariasis antigen test and by microscopic evaluation of 1 ml filtered venous blood, collected at night (21.00–03.00). The age and sex of participants were recorded as well as the time of blood collection. The GPS coordinates of all households were recorded.

Mosquito collection

Mosquitoes were collected monthly by the human landing catch method from July 2007 through July 2010 as described by Reimer et al. [35]. Villages were divided into four quadrants and houses were chosen from each quadrant every month for even sampling across the village. Mosquitoes were collected in the front of the house from July 2007 through July 2010. Quarterly collection continued in Nanaha and Yauatong through December 2011. The total collection effort ranged between 40 and 48 collection nights. All host-seeking anophelines were included in the density summary. Anopheles punctulatus comprised the majority of mosquitoes collected; additional members of the An. punctulatus group included An. koliensis, An. hinesorum, An. farauti 4 and An. farauti s. s. [36].

All temporary and permanent breeding sites were geolocated. These breeding sites were further categorized as confirmed or potential depending on the presence of anopheline larvae at the time of the survey.

Non-spatial modelling

To determine the heterogeneity at each village before and after the distribution of bed-nets, a negative-binomial distribution was fitted to both the mf count and bite data, parametrized by the mean m, and the heterogeneity parameter k. Here, a smaller k indicates a more over-dispersed distribution and a higher k indicates a less dispersed or more Poisson-like distribution. For a count n, the probability distribution is defined asThe negative-binomial distribution was fitted to village-level count data using a maximum-likelihood approach. For the nightly mosquito catches, the data were stratified before and after LLINs were distributed as a further measure of the impact of vector control on heterogeneity.

Infection transmission model

In order to understand the fitted heterogeneity k and mean biting density before and after the introduction of bed-nets in the context of disease transmission, the results were compared with an established model of LF transmission, TRANSFIL [26]. The model is a multi-scale stochastic simulation of individuals with worm burden, microfilaraemia and other demographic parameters relating to age and risk of exposure. Humans are modelled individually, with their own male and female worm burden. The density of mf in the peripheral blood is also modelled for each individual and is dependent on the number of female worms. The total mf density in the population contributes towards the instantaneous density of L3 larvae in the human-biting mosquito population. This density combined with the mosquito biting rate and an intrinsic factor that varies between individuals determines the probability of an individual being infected with a new adult worm. See [26] for a full model description.

Spatial modelling

In order to determine how much spatial variation and individual variation contribute towards differences in disease status between individuals a number of geospatial models were implemented. These models take into account distance to breeding sites, anopheline biting rates, and infection prevalence and intensity. The first group of models compare the measured disease statuses dependent on a random spatially varying risk. The second group of models combine together the biting density and disease status, by assuming that spatial variation in status is determined by the biting density alone.

Gaussian process. The underlying spatial variability (in mosquito bites, mf intensity, parasitaemia and antigenaemia) is modelled using a Gaussian process S(x) for each spatial location x. A Gaussian process describes the spatial relationship between different spatial locations and is defined as, given a set of locations {x}, the probability of observing the set {S(x)} is a multivariate Gaussian probability with zero mean and a defined covariance function. For flexibility and computational reasons, a Matérn covariance function was used, which is defined aswhere K is a modified Bessel function of the second kind and order ν > 0. The Matérn covariance function has three free parameters which control the marginal variance, the distance of spatial correlation and the sharpness of the function. These parameters can be combined to give the distance at which there is less than 10% correlation between two points, which is given as . This is referred to as the practical range.

Modelling distance to breeding sites. A second set of data includes the spatial locations of sites that may contribute to mosquito breeding. These sites include pig houses, creeks, gardens and garden houses. These data and the household bite data do not match up in terms of their geolocations. In order to circumvent this problem, we may instead use the minimum distance to a breeding site as a covariate to inform the models. All breeding sites are assumed to be equivalent whether they contained anopheline larvae or not, as there are only a limited number of sites, so as to increase power of the covariate.

Random walk latent model for breeding site distance. The relationship between minimum distance to breeding sites and the number of bites was found to be nonlinear (see the electronic supplementary material). In order to capture the full complexity observed in this relationship a more general functional form of the distance d was used, i.e. f(d). For this functional form, distances were split into a discrete lattice of points. Each lattice point k then has a corresponding coefficient x related to the nightly bites through the log intensity in the negative binomial. The assumed model was a random walk of length one, i.e. x∼N(x, σ2). The fitted function f(d) then returns the coefficient x corresponding to the nearest lattice point to the actual distance d.

Infection status and mosquito catch models. Both the infection status of each individual and the mosquito nightly bites were modelled separately using a generalized linear model including fixed effects for each observed variable (electronic supplementary material, figure S1). The general form of the model is, given a set of observations {y} at spatial locations {x}, the outcomes y have the distributionFor a general random variable X, with underlying mean m. The mean is also a random variable with the structurewhere f is a link function transforming the mean from the positive numbers to the entire real line, z is the jth covariate for the ith data-point and β is the regression coefficient for the jth covariate, which also includes an intercept. S(x) is the Gaussian process as previously defined.

For both the mosquito bite model and the mf count model, the random variable X that the observations are drawn from is assumed to be negative binomial with aggregation (heterogeneity) parameter k. The microfilaraemia and antigenaemia models have observations that are either positive (1) or negative (0), and hence the observations are assumed to be drawn from a Bernoulli random trial, i.e.The link function for the negative-binomial models was taken to be log and the link function for the Bernoulli models was taken to be logit.

The fixed effects for each model considered for the infection status models were the distance to breeding site, age of the individual and the sex of the individual. For the mosquito bite model, the covariates included were the presence or absence of bed-nets and the distance to breeding site.

Combined model. The models for the infection status were generalized to have this Gaussian process derived from the mosquito catch data, rather than including an independent spatial Gaussian process. The two-level hierarchical structure incorporates both the disease status for an individual Y and the bites at household locations B. Both observations are drawn from their own random variables X and Y , with underlying means m and n, respectively. The random variables are connected through a linear model structure of these two means with their covariates combined with a GP fitted to the transformed bite means. This GP is related to the infection status mean through the coefficient η, which measures the dependency of the infection status mean on the underlying spatial distribution of bites conditioned on the fixed effects of the bites. Mathematically, the model is defined asHere the underlying Gaussian process is assumed to capture the distribution of bites, and is fitted to both the bites and infection status simultaneously. η gives the strength of the dependency on the underlying spatial structure of bites on infection status and u is a random effect with variance σ2 used to capture the variation observed in the infection status that is not captured by the bite data.

Model fitting. The model fitting was performed using the R-INLA package [37]. This implements an integrated nested Laplace approximation (INLA) method, which is a faster alternative to Markov chain Monte Carlo (MCMC) for certain classes of models. The package also approximates the Gaussian process as a Gaussian Markov random field (GMRF), which approximates the continuous space used in GP, by a discretization of space using a triangulation based on the spatial location of the data points [38-39]. The spatial model fitting also provides an estimate of the underlying mean and variation across space.

Model comparison. In order to systematically compare the spatial and non-spatial variants of the mosquito nightly catches (bites) and mf, the Akaike information criterion (AIC) was used to assess which model produces a better fit to the data [40]. This is defined aswhere L is the maximum likelihood of the model and k is the number of model parameters. Two AIC were compared by taking the difference of the two.

Results

The LF antigen prevalence, microfilaria prevalence and intensity were collected from all consenting community members (n = 1046 individuals). Nightly biting data were collected from 170 households with 2180 sampling nights total (figure 2).

Figure 2.

Heterogeneity data used in the study. (a) Nightly bite total by village. (b) Distribution of mf count by village. (c) Spatial distribution of bites with colours on a log scale (distance approx. 13 × 6 km). (d) Spatial distribution of mf intensity with colours on a log scale (distance approx. 13 × 6 km). The spatial data indicates Yauatong is a hotspot for biting, and Albulum and Yauatong are hotspots for the presence and intensity of mf. Grey values in (c,d) indicate zero values for the nightly bites and mf concentration, respectively. (Online version in colour.)

Heterogeneity within villages

A negative-binomial distribution was fitted independently to bite counts in each village before and after the introduction of bed-nets (LLIN) (figure 3a). There is significant variation in the heterogeneity between villages pre-LLIN. A reduction in k, corresponding to an increase in heterogeneity, is observed across most villages. Where this reduction is significant is in Albulum, Nananha, Ngahmbule and Yauatong. For Peneng, there is observed a reduction in the maximum-likelihood estimate of k, although the confidence intervals of the estimate overlap.

Figure 3.

Heterogeneity of disease and mosquito bites at village and spatial levels. (a) Comparison of the heterogeneity as measured from the negative-binomial distribution before and after bed-nets. (b) Bites pre-LLIN compared with heterogeneity in mf count among individuals. The maximum-likelihood estimates for each are given as points with 95% CI given as error bars. (c,d) Spatial fits of hierarchical model for (c) mf count and (d) nightly bite rate.

To determine whether aggregated biting patterns were associated with an aggregated parasite population, we compared the pre-LLIN bite rate heterogeneity with the mf count heterogeneity (figure 3b). There was observed a large amount of variation in the mf count heterogeneity, with Nanaha and Nghambule less than 0.0125, and Peneng, Yauatong and Albulum with heterogeneity greater than 0.03. The heterogeneity in the mf counts is significantly greater than the bites in all cases. There is a positive relationship between the two heterogeneities, although the correlation is extremely weak (correlation coefficient 0.012).

Impact on elimination

The change in the vector-to-host ratio and the heterogeneity in bites after the intervention of bed-nets was explored using the stochastic model of LF transmission TRANSFIL [26]. The number of rounds to 1% microfilaraemia (which is used as an assessment for halting MDA [41]) and the prevalence at baseline before the start of any intervention were calculated across a range of bite heterogeneity and vector-to-host ratio values (figure 4). The threshold at which transmission is broken and infections are no longer sustained in the population was also calculated from these simulations. For increased heterogeneity, a smaller vector-to-host ratio, and therefore mean monthly bite rate, can sustain transmission. However, for decreased heterogeneity the vector-to-host ratio required to sustain infection increases. The effect of bed-nets can be seen to both reduce the vector-to-host ratio as well as increase the heterogeneity of bites, although for the villages in the study, the reduction in the vector-to-host ratio more than offsets the increased heterogeneity. Figure 4b highlights the number of rounds required to pre-TAS without any prior intervention. For high heterogeneity many more rounds would be required than for the equivalent bite rate at smaller heterogeneity. The impact of bed-nets can clearly be seen to rapidly reduce the number of rounds required in each village. The range in predicted rounds between villages is also large; this is, however, reduced by the introduction of bed-nets.

Figure 4.

Comparison between the theoretically predicted (a) prevalence at baseline and (b) number of rounds until reaching pre-TAS levels for varying heterogeneity and vector-to-host ratio. The red and yellow dots represent the fitted bite data before and after bed-nets, respectively. The theoretical threshold for the break in transmission is shown as a red dotted line. LLINs caused a reduction in biting density and an increase in heterogeneity, which is associated with fewer rounds of MDA to cross the threshold.

With a weak but positive association between heterogeneous biting and heterogeneous infection, we sought to determine whether this pattern could be interpreted due to spatial variation. Spatial heterogeneity was therefore explored for both the bite distribution and distribution of mf count, microfilaraemia and antigenaemia. The fixed effects considered for each of the infection status spatial models were sex of the individual, age of the individual and the minimum distance to breeding site. For the spatial mosquito catch model, both the presence of LLIN and distance to breeding site were considered as fixed effects. There was found to be no significant seasonal trend in bites, and hence month at which bite survey was conducted was not included. There was also found to be no significant trend in the time at which bleeds were taken, and hence this was also not included.

For the infection status models, age was found to be statistically significant in all cases (p < 0.001), although with a small effect in all cases. Sex of the individual was found to only be marginally significant (p = 0.15, 0.05, 0.02 for microfilaraemia, antigenaemia and mf count, respectively), with males having an increased risk of microfilaraemia, antigenaemia and mf count. The distance to breeding site was not found to be significant for any of the cases.

For the spatial mosquito catch model the presence of LLIN was found to be statistically significant, with the presence of bed-nets decreasing the coefficient by 76%. The distance to breeding site was both not statistically significant and had a very small effect compared with the intercept. The overall calculated k for the infection status mf count model was estimated at 0.05, with a standard deviation of 0.0043 and the k for the bites model was estimated to be 0.73 with a standard deviation of 0.035. The overall heterogeneity in both cases was, therefore, broadly in keeping with the estimates of the villages separately, where the mf count heterogeneity varied between 0.05 and 0.01 and the heterogeneity of bites for the villages where the mf surveys were conducted was between 0.9 and 0.3.

The fitted spatial model also provides an estimate of the mean intensity for infection status and bites across space (figure 3c,d). The estimated bite rate intensity is distributed around Yauatong, with a maximum bite rate of around 60 (figure 3d). This decreases smoothly to zero out towards Ngahmbule to the southeast and Peneng to the northwest. By contrast, both the prevalence of antigen and mf have the highest intensity around Peneng in the northwest, with a smooth decrease down towards Ngahmbule. Both prevalence spatial patterns exhibit different underlying intensity. There is a significant increase in antigenaemia around Yauatong, with a similar, but less pronounced increased risk of microfilaraemia. The uncertainty in prevalence of mf is high for Peneng, Albulum and Yauatong in the northwest, and small in the southwest (near Yauatong), where the mean prevalence is around 40%. The mf count mean is more varied than for the other distributions, with high-intensity areas around the southern part of Yauatong, Albulum and the southeast perimeter of Peneng. Ngahmbule and Nanaha have the lowest mf count matching with the lowest villages for antigenaemia and bite rate. Peneng has high levels of antigenaemia and mf count; however, it is the lowest for bites.

In order to understand the difference in spatial scale between infection status and bites, the fitted covariance structure from the mf count spatial model was compared with the structure from the bite model (electronic supplementary material, figure S1). The practical range for the mf spatial correlation was smaller than the range for the bites (0.014° to 0.008° or ≈1.5 km to ≈0.9 km). The difference in the marginal (non-spatial) variance is also great, with the variance for the mf field σ2 = 38 and the variance for the bites σ2 = 5. There was found to be no significant difference in the fitted κ (sharpness of covariance) between the mf and bite count.

The change in AIC between the spatial and non-spatial model for mf was 1308.26, whereas for the bites model was −352.03 indicating the bite distribution is better explained by a spatial model and mf count is better explained by a non-spatial model.

Combined model. In the final analysis, the bite data and infection status data were combined to produce a bite-count-dependent spatial field that is also used to predict the distribution of mf count, antigenaemia and microfilaraemia separately. The fixed effects used in the first spatial analysis were kept for the combined model. All fixed effects were found to have a similar strength and significance as in the separate models. The coefficient that described the strength of the bite rate spatial field on the outcome of the disease status spatial distribution was also calculated. These coefficients were found to be significant for mf count, microfilaraemia and antigenaemia. The largest dependency was for mf count (1.71 (1.38, 2.06)), with antigenaemia the second strongest (1.17 (1.13, 1.21)) and microfilaraemia the weakest (1.00 (0.96, 1.00)). A spatially independent random-effects term was also included in the model to account for all variation not already accounted for by the fixed effects or the spatial distribution produced by the bites model. These were found to be negligible in all cases.

Discussion

The primary aim of the study was to determine whether heterogeneous biting activity drives the observed heterogeneity in LF prevalence and intensity. We observed a much more complex picture than previously expected, with heterogeneity being driven by both spatial biting patterns and individual processes. The secondary aim was to determine how aggregated biting patterns are influenced by vector control and determine the implications for LF elimination. Combining a statistical and modelling approach, we demonstrate that vector control increases heterogeneity while also uniformly reducing biting. This resulted in decreased variability in the predicted number of years required to achieve elimination between neighbouring villages in Papua New Guinea.

Heterogeneity poses numerous challenges to global elimination programmes that rely on broad-scale mapping to inform distribution of community-wide interventions. Heterogeneities in exposure and infection are well-known drivers of persistent disease transmission. Diseases such as LF have complex ecological interactions that can lead to threshold behaviour, where sustained transmission is dependent on biting density or parasite load [24]. The basic reproduction number R0 is expected to be greater under heterogeneous biting [42]. The probability of re-introduction of a disease in a fully susceptible population can also have a nonlinear relationship with the basic reproductive number of the infection [43]. The potential for transmission can be dependent on heterogeneous exposure (some people bitten by mosquitoes more than others), poor mixing (non-random contacts between hosts and mosquitoes) and finite population sizes (each host can contribute at most one new infection towards the population total) [44].

Heterogeneities in infection can complicate disease surveillance programmes since public health infection mapping is usually performed at village/town level, while interventions are often implemented at a broader administrative level. Data aggregated by population can hide the true patterns, which are more apparent when data are considered in a spatially explicit fashion. It is therefore imperative to understand the relationship between the underlying heterogeneity for these scales and how this heterogeneity impacts the efficacy of interventions [17]. Spatial and individual heterogeneity should be considered in order to ensure implementation policy is appropriate to local transmission and epidemiology (table 1).

Individual heterogeneity of infection in a population reduces the likelihood that community-wide interventions are protecting the highly exposed. As a result there is a lower threshold mf prevalence required to break transmission and a greater likelihood that high-density infections in a few individuals can seed new infections (figure 1). Spatial heterogeneity can also conflate an elimination campaign, as there may be regions of high disease burden adjacent to regions with low rates of transmission. This poses a significant challenge when sentinel and spot-check sites are used to determine the prevalence for an entire region for the purposes of implementation [45]. This may lead to limited resources being wasted on MDA distribution in villages that have lower than threshold prevalence, while possibly missing areas that will require a longer duration of MDA to break transmission.

We observed a strong spatial correlation between biting density and antigenaemia, which captures current or prior presence of adult worms, including amicrofilaraemic infections. While biting density was also significantly associated with microfilaria intensity, this association was weaker than for either antigen or microfilaria prevalence, indicating a weak relationship between high exposure to mosquito bites and intensity of infection. Although the distance to breeding site was not statistically significant for infection status or bites, there is a stronger mean effect from the bite counts than from the infection status. The distance to breeding site would naturally be more associated with biting density, whereas infection status is more strongly dependent on other factors, hence the weaker regression effect. As current microfilaria intensity is the result of fecund adult worms that have established after years of exposure to infective bites, long-term changes in the mosquito population may not be taken into account from the current distribution of breeding sites. Previous studies in these study villages have shown that high density infections are associated with reduced strain variation, and are not necessarily due to multiple adult worms [46]. Host immunity or W. bancrofti strain fecundity probably play a greater role in the intensity of microfilariae.

In our study heterogeneity in bite exposure varied substantially from village to village before the LLIN distribution, and this was associated with wide-ranging predictions on the number of rounds of MDA required to break transmission. Heterogeneity increased significantly after the introduction of LLINs in all villages except the one village with very low pre-LLIN biting rates, resulting in a very similar heterogeneity parameter across the five villages. The greater heterogeneity observed post-LLIN in all communities is associated with a transmission threshold at a lower mean biting rate. In this particular transmission system, the reduction in vector density caused by the LLIN distribution compensated for this change in threshold biting rate. However, it does highlight the extreme importance of considering heterogeneity in elimination strategies, because changes in heterogeneity cause elimination targets to move. Globally there are 54 countries engaging in preventive chemotherapy for the elimination of LF [47], and each of these countries will need to decide when to stop MDA and switch to long-term surveillance. That decision will be made based on the available evidence that microfilaria prevalence has fallen below 1% in sentinel villages, but there is a risk that the minimum duration of MDA will differ significantly between neighbouring villages.

Statistical models can help shed light on the complex factors that contribute towards heterogeneity in disease-burden and transmission [48]. While heterogeneity in the vector population can lead to a particular aggregated exposure, there is an assumption that aggregated exposure leads to an aggregated burden of infection. In transmission models for macroparasites, there is usually an explicit distribution of risk across individuals in the same community, but it is often parametrized against the resulting distribution of infection burden rather than vector data [26,49,50]. Here we have demonstrated that factors including individual and spatial heterogeneity can all contribute towards the perceived variation in the aggregated distribution (figure 1). Although here bites are measured at the household as opposed to the individual-level, these results suggest that modelling should more explicitly take into account other aspects that lead to the final distribution, such as strength of infection-blocking immunity. However, we acknowledge that variation in infection burden has been much more frequently measured than variation in vector biting rates, due to understandable practical challenges, and therefore this may be the best way to proceed in the absence of more vector data. In other vector-borne disease, heterogeneity in exposure is rarely explicitly included in transmission models, despite a number of measurements and theoretical studies highlighting its importance [1,3,5,7]. Our study once again demonstrates the likely impact of these heterogeneities, and the need for more epidemiological and entomological studies performed at the same time and in the same place. While these data are challenging to interpret, larger studies would allow us to identify the right correlates of current transmission rates and the likely impact of control.

Conclusion

Understanding sources of heterogeneity is important both in disease modelling and ultimately in the control and elimination of a disease. We have comprehensively demonstrated here that individual and spatial heterogeneity can impact disease prevalence and intensity in different ways, and has direct implications to policy. There are many logistical and financial challenges to sustaining long-term MDA campaigns in a setting like Papua New Guinea, where communities are hard to reach and departments of health have competing priorities. The risks of resurgence if programmes fail to break transmission thresholds would compromise the gains already made by global elimination efforts. Therefore, knowledge of the degree of heterogeneity is necessary to understand where transmission thresholds lie, and understanding the sources of heterogeneity is essential to designing and delivering interventions with the greatest chance of success.