Modelling the distribution and transmission intensity of lymphatic filariasis in sub-Saharan Africa prior to scaling up interventions: integrated use of geostatistical and mathematical modelling

Data sources and data selection

Data on the prevalence of LF infection were identified from searches of the formal and informal literature and direct communication with LF control programmes. Details of the search strategies, inclusion criteria, data abstraction and geolocation procedures are provided by Cano et al. [12]. In brief, only population-based survey data based on random sampling were included, whereas data from non-random data, including surveys from hospitals, prisons, mental institutions or military facilities, were excluded. In the current analysis, only surveys conducted prior to the implementation of countrywide, population-based MDA were included. Infection prevalence was defined as either (i) the proportion of surveyed individuals with detectable microfilaraemia or (ii) the proportion of surveyed individuals with detectable antigenaemia. A table with a brief description of the surveys compiled and used in this work is provided as supplemental information (see Additional file 1).

Microfilaraemia (mf) was estimated from examination of thick blood films collected during night blood surveys to detect the presence of microfilariae using microscopy [14, 15]. It is generally assumed that the sensitivity of methods for detecting microfilaraemia is influenced by the volume of blood sampled and previous authors have adjusted prevalence estimates according to blood volume [16]. We therefore investigated this issue further but found that although estimates of infection prevalence varied according to blood volume, we were unable to derive consistent adjustment factors and therefore did not adjust prevalence estimates according to blood volume (see Additional file 1). It is also known that concentration and counting chamber methods have greater sensitivity [17] and therefore where estimates were derived using both concentration methods and thick blood smears (n = 10 surveys), we used only the thick smear data.

Wuchereria bancrofti antigenaemia was typically estimated using an immuno-chromatographic card test (ICT) [14, 18]. These tests are more sensitive than mf detection and can be conducted on blood collected at any time of the day and therefore since 2000 have been the diagnostic method of choice for mapping the distribution of LF caused by W. bancrofti. Recent work has highlighted potential cross-reactivity of the ICT test with Loa loa [19], therefore we excluded ICT-based surveys conducted in areas of L. loa transmission (n = 314), as defined by Zouré et al. [20].

Data derived from parasitological blood examinations during day time (n = 71), ELISA (n = 70), clinical examinations (n = 24) and other molecular or antibody-based tests (i.e. PCR, IFI, new rapid test) (n = 17) were excluded because of the lack of comparable data. Finally, we excluded data from Egypt since this country is no longer considered endemic for LF and has a different transmission ecology from SSA [21, 22]. Figure 1 summarises the survey data selection and Fig. 2 shows the geographical distribution of survey data by diagnostic method.

Covariate variables

There are a number of climatic, environmental and demographic factors that influence the transmission and epidemiology of LF, including temperature, humidity, elevation and population density [12, 23–25]. We generated a suite of climatic and environmental covariates from surfaces interpolated from meteorological stations or remote sensing imagery, including estimates of mean, minimum and maximum temperature and precipitation at 1 km² resolution [26], averaged long-term estimates of land surface temperature [27], elevation at 1 km² resolution [28], and aridity index [28], averaged enhanced vegetation index for the period 2000 to 2012 [27], and land cover data [29]. Estimates of population density were obtained from the Gridded Population of the World [30] and the United Nations Environment Programme [31], which we used to classify areas as urban, peri-urban or rural areas, based on the assumption that urban extents (UE) have a population densities ≥1000 persons/km², peri-urban >250 persons/km² within a 15 km distance from the UE edge, and rural <250 persons/km² and/or >15 km from the UE edge. From these datasets of population density, population growth rate for the period 1960 to 2010, as a measure of the change rate in population over this period, was also calculated. A detailed description of each covariate and source is provided in the supplemental file (Additional file 2: Table S3). Each covariate grid was resampled to a 10 km spatial resolution using bilinear interpolation for continuous surface and a majority approach for categorical data [32]. Covariate data were extracted to survey locations using ArcGIS Desktop v10.2 (Environmental Systems Research Institute Inc., Redlands CA, USA).

Vector distribution maps

LF is transmitted mainly by mosquito species belonging to the Anopheles, Culex, and Mansonia genera and to lesser extent Aedes, Coquillettidia and Ochlerotatus genera [33, 34]. Anopheles mosquitoes are the main vectors of LF through much of west and central Africa and inland east Africa, whereas Culex species are the main vectors in coastal east and southern Africa [33]. Studies have shown that the survival, parasite uptake and development of infective L3 stages and overall transmission potential varies by vector species [35–37]. In order to capture these species differences, we developed maps of the distribution of each dominant LF vector species: Anopheles, Culex and Mansonia.

Maps of the distribution of mosquitoes belonging to the An. gambiae complex and An. funestus complex were obtained from the Malaria Atlas Project project (http://www.map.ox.ac.uk/) [38, 39] and a binary map displaying the distribution of each complex was created. Maps of Culex and Mansonia mosquitoes were obtained from the VectorMap project (http://www.vectormap.org/), which collates collection records of major vector insects and uses maximum entropy ecological niche modelling to develop occurrence maps [40, 41]. The maps present the probability of occurrence of Cx. quinquefasciatus, Cx. pipiens, Cx. univittatus, and Mansonia Africana, and we used a 90 % probability threshold to indicate dominance of a particular species according to a set of presence records obtained from the Global Biodiversity Information Facility database [42]. Each of the species-specific binary maps were combined to produce gridded maps of mosquito distribution by genus which, in turn, were combined into a single map of the different LF vector species at a 5 km spatial resolution (Fig. 3). The developed maps of Anopheles, Culex and Mansonia genera were subsequently incorporated as covariates in the goestatistical and mathematical modelling.

Geostatistical modelling approach

The prevalence of mf and prevalence of antigenaemia were modelled separately due to a poor correlation between outcomes and the inability to predict one from another, as detailed elsewhere [43]. Bayesian geostatistical models, which included both fixed and random effects, were used to predict the spatial distribution of each outcome. Fixed effects quantify the effects of the covariates on LF infections, whereas random effects account for unexplained spatial variation whose structure can eventually help identify anomalous areas of high or low risk that can be further investigated.

where Κ _υ is the modified Bessel function of second kind and order υ > 0. The integer value of υ determines the mean square differentiability of the underlying process. υ is usually fixed since typically it is poorly identified in applications. σ ² is the marginal variance and k > 0 is a scaling parameter related to the range ρ, the distance at which the S(x _i ) and S(x _j ) are almost independent. In particular, we used the empirically derived definition \( \rho = \sqrt{8\upupsilon}/k \), with ρ corresponding to the distance at which the spatial correlation is approximately 0.1. We also include in the model a spatially unstructured random effect, u _i , that captures the effects of unmeasured characteristics that affectall members in the survey. u _i are assumed to be independent zero-mean Gaussian distributed with precision τ _u .

where |a| <1, and ξ(x _i , 1) follows the stationary distribution of a first-order autoregressive process, namely N(0, σ _w ² /(1 − a ²)). W(x _i , t) is a spatial correlation term. Each w(x _i , t) follows a zero-mean Gaussian distribution. The w(x _i , t) is temporally independent but spatially dependent at each time t, with Matérn covariance function. u _it denotes the unstructured random effect and are assumed to be independent zero-mean Gaussian distributed with precision τ _u .

Variable selection and model development

where \( l\left(\widehat{\theta}\right) \) is the maximum log-likelihood function and k is the number of parameters. After that, we explored further simplification of the model by backward elimination of selected variables until it was no longer possible to reduce AIC by elimination of any of the remaining variables.

The age group (either below or above 15 years of age) of individuals studied on the prevalence surveys was considered at every step of the modelling process. In addition, population density was considered for the periods at which the surveys were undertaken. Since ICT-based data were available for surveys conducted from 1990 onwards, estimates of population density at 2000 were used. For mf data, in order to take into account temporal changes in population density, gridded maps of population density at 1960, 1980 and 2000 were used accordingly with the period of survey (1950–1969, 1970–1989, 1990–onwards).

The selection of covariates at this stage in the overall model-fitting process ignores spatial and temporal correlation, and as a result is likely to over-state the statistical significance of covariate effects. In the final stage of the model-fitting process we re-assess the covariate effects and their significance within a spatial mixed model for the ICT data and a spatiotemporal mixed model for the Mf data that take into account the spatial and spatiotemporal correlation, respectively.

Model validation

We assessed the predictive ability of the model using a leave-one-out cross-validation procedure. In this approach, a single observation is retained as the validation data, and the spatial model is fitted to the remaining data. Then, the observation in the validation data is predicted using the fitted model. The validation data needs to spatially represent the whole region where the prevalence is predicted. Therefore, instead of repeating the cross-validation procedure using each observation once as the validation data, we used each of the locations of a spatially representative sample of the prediction surface. To obtain a valid data set, 20 % of the observations were sampled without replacement where each observation had a probability of selection proportional to the area of the Thiessen polygon surrounding its location, that is, the area closest to the location relative to the surrounding points.

where m is the number of observations in the validation data, p*(x_i) is the predicted prevalence, and p(x_i) is the observed prevalence at location x_i, i = 1, …, m.

Implementation and spatial prediction

where G is the number of vertices in the triangulation, ψ _g (⋅) are piecewise polynomial basis functions on each triangle, and {S _g } are zero-mean Gaussian distributed weights. The INLA and SPDE approaches are easily applied thanks to their implementation in the R-INLA software package [5].

We assigned a flat improper prior for the intercept, \( {\beta}_0\sim N\left(0,{\tau}_{\beta_0}^{-1}\right) \) with \( {\tau}_{\beta_0}=0 \), and independent vague Gaussian priors with fixed precision for all other components of the fixed effects, β _i  ∼ N(0, 1/0.001), i = 1, … p. The smoothness parameter ν was considered fixed to 1 implying a continuous domain Markov field. In the spatial model, the vector of weights S = (S _1,..., S _G )′ is assigned a Gaussian distribution, S ∼ N(0, Q ^− 1) where Q is a sparse precision matrix depending on the Matérn covariance function parameter κ and variance σ ². In the spatiotemporal model, for ξ = (ξ ₁, …, ξ _n ) ', we use a Gaussian prior with zero mean and precision matrix depending on the autocorrelation hyperparameter a, k and σ _w ² . To complete the models, we assign τ _u a vague Gamma prior.

The model output provided the posterior distribution for each of the model parameters. We summarized these distributions to obtain the posterior mean and 95 % CI of the fixed effects and hyperparameters. Predictions of infection prevalence were provided at a spatial resolution of 10 km. Maps showing the predicted LF prevalence as well as its uncertainty were generated by summarizing the posterior distributions of the prevalence obtained at each of the prediction locations. Specifically, the predicted prevalence was represented as the posterior mean, and its uncertainty was represented using the 95 % CI. Predictions were standardized to provide estimates for entire communities (children and adults combined) and for the period 2000-onwards.

Mathematical modelling of the intensity of LF transmission

The prevalence of LF infection provides a useful metric to guide the planning of control, but provides limited insight to the dynamics of transmission. Instead the intensity of infection transmission is best quantified by the basic reproduction number, R₀, which for macroparasites such as LF can be defined as the average number of female offspring per adult female worm surviving to reproduction (in the absence of density-dependence i.e. phenomena such as host immunity or worm mating probability, which accelerate or curtail the production of parasite life stages) [48]. Existing analytic expressions for R₀ for LF, and helminths with similar natural history such as onchocerciasis, are defined in terms of average worm burden or microfilarial load in a community [49–51], yet the majority of LF studies only present data on the prevalence of infection. In the framework published by Gambhir et al. a simple, differential equation transmission model of LF is used to estimate R₀ for a given prevalence value [6]. This model assumes various density-dependent functions which alter the rates of transition between state variables in the model (for example, the rate of parasite establishment and rate of parasite fecundity). The model parameter prior distributions were defined using data from the literature and parameter posteriors were found by fitting to baseline mf prevalence data from low, moderate and high transmission settings in Tanzania [52] – settings which are representative of much of SSA.

For a setting with a particular endemic prevalence, there is an underlying bite rate and force of infection which leads to this endemic equilibrium. This parameter contributes to R₀ (as explained in detail below), but can only be calculated indirectly through the method described in detail by Gambhir et al. [53] and outlined here.

where λ, α and β are the ‘immigration’ rates of each of the life stages (adult worms, microfilariae, and L3 infective larvae respectively; λ is a composite parameter, incorporating the mosquito-human biting rate which is estimated using the methodology describe above; μ, δ and σ are the death rates of each of the life stages respectively; and f ₁(I), f ₂(W), f ₃(W), f ₄(M) represent the modifying density-dependent functions acting on each of the respective parasite life-stage intensities (W, the number of adult worms per definitive host, M, the mf load per host, refer to parasite life stages within the definitive host population, whereas the L3 infective larval load per mosquito (larvae develop through L1 and L2 stages but only become infective once they reach the L3 stage), L, refers to the parasite life-stage within the vector host. The host immunity variable, I, increases in magnitude at a rate equal to the adult worm burden but can also decay over time, allowing the model to capture the waning of immunity). Note that for R _0, host immunity levels and worm burdens are assumed to be zero as this expression calculates the number of offspring in a wholly susceptible population. The above mathematical model was initially fitted separately for settings where either Anopheles or Culex species are the predominant vector, based on the develop vector distribution maps (Fig. 3). However, the relationship between the mf prevalence and R₀ is practically identical for Anopheles or Culex species, as the model assumes similar mosquito life expectancy for each species, and therefore we present a single relationship between prevalence of mf and R₀. Finally, the threshold mf prevalence value above which LF transmission can persist, which is higher than the prevalence at which R ₀ = 1 due to the density dependencies assumed in the modelling framework [53], was calculated.

Collated survey data

Vector distribution maps

Figure 3 presents maps of the predicted distributions of Anopheles, Culex, and Mansonia species. Anopheles mosquitoes of gambiae and funestus complexes are widely distributed across sub-Saharan Africa (Fig. 3a). In contrast, M. africana and Culex mosquitoes show a more limited and distinct distribution: Culex mosquitoes occur in eastern Africa, east coast of Madagascar and in restricted areas of west Africa (Fig. 3b), whereas M. africana occurs in west and middle Africa, and coastal areas of east Africa and Madagascar (Fig. 3c). A wider presence of Culex mosquitoes is predicted in Nigeria, north-east of Cameroon and north of Angola. Areas where the distribution of each of the major LF mosquito vectors occur across west Africa and coastal areas of east Africa (Fig. 3d).

Spatial prediction of LF prevalence

In the spatial model for ICT prevalence, the model selection process excluded land surface temperature, slope and distance to water bodies, while in the spatiotemporal model for mf prevalence no covariates were excluded during variable selection, as based on the AIC criterion (Additional file 2: Table S5). The predicted prevalence of LF, based on mf and ICT-based data, are presented in Fig. 4a and b, respectively, along with estimates of 2.5 and 97.5 % quantiles. Overall, predicted mf prevalence is lower than predicted ICT prevalence across sub-Saharan Africa, not exceeding the 5 % threshold in most endemic areas. Only a few small pockets of mf prevalence higher than 30–40 % are predicted in east Africa, at the north of coastal Tanzania and the southeast coast of Madagascar. Broader areas of high mf prevalence are predicted in west Africa; south of Mali (Sikasso region), large central areas of Benin and north-west of Ghana (bordering with Benin) and at the south of Abuja, the federal capital of Nigeria.

The ICT-based model indicate large areas of high prevalence (above 40 %) are predicted along the coast of east Africa, from southern Kenya to central regions of Mozambique, and the east coast of Madagascar. High prevalences of antigenaemia are predicted across large areas of west Africa and pockets at the north-west of the Democratic Republic of Congo (DRC) and central Congo. In general, lower (<5 %) estimates of ICT prevalence are predicted in inland areas of east Africa (Ethiopia, Rwanda, Burundi, central and south of Uganda) and middle Africa (southern and central DRC).

Predictive distribution of R₀

The relationship between the prevalence of mf and R₀ is presented in Fig. 5. Figure 6 shows a predicted map of R₀ based on mf prevalence. The estimated R₀ values vary from 2.7 to 30 across sub-Saharan Africa. The areas of low and moderate intensity of LF transmission (R₀ value lower than 10) predominate in all endemic areas and the highest intensity is predicted in restricted areas of west Africa (Nigeria, Ghana and Burkina Faso), along the coast of east Africa, from north-east of Kenya to southern Mozambique and larger areas in central and eastern Madagascar.

Strengths and challenges in LF mapping

The presented work has a number of strengths and limitations. Our use of Bayesian model-based geostatistics provides a flexible and robust approach for spatial modelling that takes into account the spatial dependence structure of the data and provides a formal expression of uncertainty in the prediction estimates [67]. Traditionally, Bayesian predictive inference has been implemented via MCMC methods which make inference tractable for complex models but present problems in terms of convergence and computational time [68]. We overcome these issues by using the INLA approach [46] which is a computationally effective alternative to MCMC designed for latent Gaussian models. By using a combination of analytical approximation and numerical integration, INLA produces fast and accurate approximations to posterior distributions. The analysis of spatial point data is possible by combining INLA and SPDE approaches [47], whereby a continuously indexed Gaussian field is represented as a discretely indexed Gaussian Markov random field (GMRF) using a basis function representation defined on a triangulation of the domain of interest.

Intensity of LF transmission is influenced by a complexity of factors, yet our models incorporated only environmental and demographic drivers of LF transmission in pre-intervention scenarios. The models do not incorporate socioeconomic status [69–72] or the coverage of LF interventions, including MDA and vector control, which will influence transmission, especially once interventions have been scaled up. The role of vector control measures in modifying patterns of LF transmission has recently been highlighted in The Gambia where the large-scale distribution of bed nets for malaria control has, in part, contributed to the elimination of LF in the country [73]. In coastal Kenya, large-scale distribution of ITNs for malaria control has been proposed to sustain reduction of LF infection levels, even after MDA is interrupted [74]. An increasing number of studies have demonstrated the impact of repeated rounds of albendazole and ivermectin alter patterns of transmission [75] and over time dramatically reduce transmission [76, 77]. A methodological challenge for incorporating intervention-related factors into future models is that much of the available data are presented at an area-level. Bayesian hierarchical models have been used to address this problem by providing a natural way to combine data from different sources taking into account the different uncertainties [78–80]. However, such models have three important limitations. First, they rely on MCMC methods for Bayesian inference which are computationally intensive and may become unfeasible for large datasets. Second, few approaches for temporally misaligned data exist and they do not account for correlation in spatial random effects over time. Third, correlation studies use the predictions from several models as covariates in regression models and as such predictions contain error as the predicted values do not equal the true exposures. These important issues indicate that there is considerable scope to develop new model-based methods for the analysis of misaligned data.

The vector distribution maps are also not without their limitations. First, different methodologies were used for different genera: maps for Culex and Mansonia were constructed by using maximum entropy ecological niche (MaxEnt) modelling [40]; and boosted regression tree (BRT) modelling was used to modelling environmental suitability for Anopheles mosquitoes [38]. Second, there is a lack of data for some areas of central and northeast Africa, introducing imprecision into maps.

Finally, the presented map of R0 is based on a mathematical model of transmission dynamics [53]. The relationship between prevalence and R0 is non-linear and there are very high values of R0 for moderate to high levels of prevalence, suggesting great challenges for control. However, it should be remembered that these values are for a population when there is no immunity, unlike endemic communities, and that we do not yet have a reliable estimate for how long this immunity will last. Therefore the high values of R0 may not represent such a large challenge for control in the short term. For recrudescence or re-emergence however, they may represent a higher risk. The model uses prevalence to estimate transmission intensity because of the availability of prevalence data. Alternative measures of transmission intensity include annual transmission potential and annual biting rates [81], but these are very rarely measured. Ideally, the model would have been parameterised with data from a range of transmission settings, but there are a lack of detailed data and we recognise that transmission dynamics will vary in settings other than Tanzania from where data were used to parameterise the model.

Mapping LF to guide control and surveillance

The developed map of R₀ can be used to identify areas of higher transmission intensity where additional control efforts may be required to achieve the elimination of LF infection (i.e. implementing vector control or extending MDA rounds). Our model predicts a marked geographical heterogeneity in the intensity of transmission in SSA. Current MDA strategies do not account for this spatial heterogeneity, and therefore coverage targets and implementation schemes are common for all endemic areas [3]. The precise level and duration of treatments to achieve LF elimination in different endemic regions remains unknown [82], such that it is difficult to decide when to stop ongoing MDA interventions [51]. Indeed, although some model predictions suggest that LF transmission can be interrupted by annual MDA alone, it is not clear that this can be achieved everywhere with 4–6 yearly rounds of MDA [83].

In recent years, some experts have suggested that current control programmes based on standard MDA schemes should give rise to more tailored control interventions adapted to the local specificities of LF transmission [6, 14, 84]. Our maps provide programme managers with a thorough picture of the initial level of LF endemicity and potential geographical variation in the intensity of transmission prior to the implementation of large-scale interventions. The R₀ maps also provide insight into (i) the risk of persistent transmission despite repeated MDA and highlight that pockets of residual transmission may occur in settings where overall transmission has declined [77], and (ii) the possibility of recrudescence, especially where transmission is highly heterogeneous [85]. Careful use of the developed maps can guide decision makers and ultimately programme managers first to define the most appropriate control strategy and adapt it to local conditions of transmission and second to identify areas which may require special attention during monitoring and surveillance stages. How to integrate an understanding of the spatial heterogeneity of LF into the design of optimal LF surveillance is the subject of ongoing work.