Introduction

Onchocerciasis (“river blindness”) is a neglected tropical disease (NTD) caused by the filarial parasite Onchocerca volvulus. It is transmitted through bites of infected black flies of the genus Simulium. These flies typically breed in fast-flowing waters, because the high oxygen content of this environment is required for larval development [1]. Infected humans may experience visual impairment that can progress to total blindness. Onchocerciasis can also cause severe dermatitis, skin depigmentation, subcutaneous nodules, epilepsy, and excess mortality [2]. The social consequences of onchocerciasis can be devastating [3].

Côte d’Ivoire is endemic for several NTDs [4]. Efforts to reduce the onchocerciasis burden started in 1975 with vector control activities that were coordinated by the Onchocerciasis Control Programme (OCP) in West Africa [5,6]. Aerial spraying of insecticides was carried out mainly in savannah areas in the northern and central parts of the country where blinding onchocerciasis was prevalent. Since 1992, onchocerciasis control has been largely based on the strategy of administering ivermectin to eligible populations in endemic communities [6], although vector control continued beyond 1992 in some river valleys. Ivermectin distribution was initially piloted by nongovernmental organizations. Between 1995 and 2016 the strategy of community directed treatment with ivermectin (CDTi) was adopted, with initial villages receiving this intervention in 1996. This was provided by the Ministry of Health in collaboration with OCP and the African Programme for Onchocerciasis Control (APOC) that supported ivermectin distribution in countries of Africa that were not covered by OCP. CDTi uses local volunteers to distribute ivermectin in endemic communities. This approach enabled control efforts to be extended to areas that were previously excluded from the aerial spraying programme. APOC coordinated CDTi in some 19 countries that provided an estimated 1 billion ivermectin treatments and prevented some 2 million cases of blindness [7]. APOC ended in 2015, and control efforts for onchocerciasis and other NTDs in the AFRO region of the World Health Organization (WHO) are now coordinated by the Expanded Special Project for Elimination of Neglected Tropical Diseases (ESPEN) [8].

Thus, Côte d’Ivoire has had 45 years of intervention that started with vector control in 1975–1991 that was mostly replaced by ivermectin MDA/CDTi since 1992 [4]. Years of civil unrest (2002 to 2007) interrupted public health interventions [9], and this was especially true in the rebel-held northern regions of the country. However, as onchocerciasis has been eliminated in four countries in the Americas and from several foci in Mali [10], Senegal [10,11], Nigeria [12], Sudan [13,14], and Uganda [15,16], there is increasing evidence that the interventions for onchocerciasis can lead to elimination in some areas and that elimination targets set by the WHO are feasible. It is important, therefore, to identify countries and regions where elimination may be achieved using the current interventions and regions where modified strategies may be needed.

Baseline endemicity is a key determinant that affects the feasibility of onchocerciasis elimination and the time required for treatment with ivermectin [17]. Community microfilarial load (CMFL) is widely regarded as an important measure for categorizing onchocerciasis endemicity [18]. This is a measure of the intensity of infection that is calculated as the geometric mean number of MF per skin snip in adults aged ≥20 years within a community. During early stage of intervention, CMFL is considered a robust measure for determining the true epidemiological situation within an endemic population. However, MF (skin snip) prevalence data is commonly used for classifying endemicity, as only meso- and hyperendemic areas are considered to have a high risk for blinding disease.

Previous modelling studies suggested that either annual or biannual MDA, depending on endemicity, coupled with high levels of therapeutic coverage, should be adequate to achieve the elimination threshold suggested by APOC [19]. Maps that classify infection levels could be useful for identifying areas with the highest risk of infection. They can also identify areas that might be more susceptible to recrudescence of infection following local elimination, either from nearby infection hotspots or from more distant foci within the same transmission zone.

Recent advances in disease prediction provide methods for producing continuous maps with high resolution spatial scales using various modelling approaches [20–24]. These geostatistical models are able to predict infection across a large geographical space using a suite of potential disease drivers such as remotely-sensed climate and environmental data together with relevant socio-demographic data to improve model predictions. Furthermore, continuous maps have been published for a variety of helminth infections [24,25] including vector-borne parasites [20,26] such as lymphatic filariasis [27–30]. O’Hanlon et al recently presented a geostatistical map of the pre-control prevalence of onchocerciasis in OCP countries in West Africa [25]. Although that study was a valuable contribution, the ground-truth data used to build the models only considered data from savannah areas, and that excluded important endemic areas in Côte d’Ivoire. Model predictions were extrapolated to the forested southern parts of Côte d’Ivoire where there were no associated survey data. Uncertainty in model predictions increase as one moves farther away from ground-truth data. In addition, data based on MF prevalence (skin snip) alone was used to determine onchocerciasis endemicity. Modeling based on CMFL should provide a more complete picture of the disease burden. Finally, as maps are usually built with high-dimensional satellite data, these data are usually non-linear. Modelling within a machine learning framework can efficiently handle complex relationships between predictor and response variables [22,31].

Therefore, to delineate remaining infection hotspots in Côte d’Ivoire, we modelled the CMFL and MF prevalence prior to and during the use of ivermectin MDA/CDTi. We used a trained quantile regression forest (QRF) model to: i) predict the MF prevalence and intensity of onchocerciasis in unsampled locations at 5 × 5km spatial resolution; ii) identify important socio-demographic and other factors that correlate with onchocerciasis endemicity; and iii) delineate areas at risk of ongoing transmission and potential infection hotspots.

Methods

Study areas and onchocerciasis survey data

Côte d’Ivoire is divided into four main ecological zones: savannah in the northern region, pre-forest in the central region, forest in the southern region, and mountainous areas in the west-central region. Blinding onchocerciasis predominantly occurs in the northern region, and early intervention (OCPs use of aerial spraying of vector breeding sites) was focused in eight districts in this region from 1975 to 1991. Following the commencement of mass treatment with ivermectin from 1992, intervention was extended to villages in 53 districts located in the pre-forest and forest regions.

Methods used in the onchocerciasis epidemiological surveys have been previously described [4]. Briefly, proximity to rivers (≤ 5 km) and population size (less than 2000 people) were key determinants that informed the selection of survey villages. Two to four villages were selected for surveys per district depending on accessibility and district size. Selection of survey sites was mainly based on accessibility, although expert knowledge from local health staff and prior onchocerciasis endemicity data were also used to spread survey locations to cover as much geographical space as logistically possible.

For each survey, parasitological examinations included collection of a single skin snip from each posterior iliac crest using a Holth-type corneoscleral punch according to WHO protocols. In sentinel villages where more than one survey was conducted, we chose the infection estimates from the latest survey conducted. The MF load for each person was defined as the arithmetic mean MF count per mg skin.

Figs 1 and 2 show the geographical distribution and endemicity level of survey villages, grouped by period.

Fig 1

Location of survey sites in Côte d’Ivoire.

Plots on the left and right are MF prevalence surveys from 1975 to 1991, and 1992 to 2016, respectively.

Fig 2

Location of survey sites in Côte d’Ivoire.

Plots on the left and right are CMFL surveys from 1975 to 1991, and 1992 to 2016, respectively.

Socio-demographic, climatic and environmental data

Interpolated climate and remote sensing data used in combination with statistical models to map the distribution of diseases as well as other health indices has accelerated greatly in the past decade. These data layers are created by collating large amounts of point-level data from weather stations globally. Smoothing algorithms are used to produce continuous maps. We downloaded climate variables related to temperature from the WorldClim database [32]. Here climate layers are presented as smooth maps of mean monthly climate data obtained for the period 1950–2000 from thousands of weather stations globally (data from ~15,000 weather stations were used to estimate the minimum and maximum temperature variables).

Access to communities was a key determinant for selecting survey sites for onchocerciasis surveys. In order to account for this in our analysis, we processed data from the WorldPop repository [33,34] that measures remoteness and proximity to human settlements. Variables such as travel time to nearest large settlement, proximity to major roads and night-time lights were considered. Furthermore, gridded population density estimates, Euclidean distance to permanent inland water bodies and the ocean, terrain slope and altitude were all downloaded from the WorldPop repository. Vegetation cover types (according to the United Nations land cover classification system) were extracted from the GlobCover project at the European Space Agency [35]. Here maps are derived by an automatic and regionally-tuned classification of a 300m medium resolution imaging spectrometer (MERIS) sensor on the ENVISAT satellite mission.

Because onchocerciasis is known to be endemic in agricultural communities, we considered covariates such as designated croplands and areas with established irrigation infrastructure. These data layers were downloaded from the Global Map of Irrigation Areas of Food and Agricultural Organization (FAO) [36]. Here, a digital map of irrigation areas was developed by collating over 10,000 sub-national irrigation records from census surveys and reports available to the FAO and World Bank to create geo-spatial layer for irrigation density (defined as continuous grid-cells that are equipped for irrigation) [37]. Finally, data for household wealth and maternal education were downloaded from the Socioeconomic Data and Applications Center, Columbia University [38], and data on housing type were obtained from the Malaria Atlas Project, University of Oxford [39]. Here housing type was categorized as either being built with finished materials (e.g. cement, bricks, or tiles) or built with natural or unfinished materials (e.g. earth, sand, or palm flooring) [23]. These variables were considered as proxies for wealth.

Changes in the covariates over time were assumed to be negligible for this analysis. For example, the WorldClim dataset uses interpolated data averaged over a 50-year period. Although this dataset pre-dates some of the epidemiological surveys in our study, we assumed that changes that may have occurred in this covariate have been negligible and that it accurately reflects climatic conditions during the study period. Other covariates considered in our model were modelled estimates, and we treated any temporal trends that may exist as negligible. Table 1 shows a list of the covariates and their sources.

All input grids were resampled to a common spatial resolution using the nearest-neighbor algorithm [40] and then clipped to align to the geographical boundaries of Côte d’Ivoire. Raster manipulation and processing were done using the raster package in R [41].

Table 1

Environmental variables used in analysis and their sources.

Variables	Source
Distance to permanent inland water bodies (rivers and streams)	WorldPop [33]
Distance to coastline
Distance to major roads
Night-time lights
Population density
Slope of terrain
Altitude
Travel time to cities	Malaria Atlas Project [23,39]
Housing type
Cropland	Food and Agriculture Organization of the United Nations [36,42]
Irrigation
Vegetation cover	European Space Agency [35]
Household wealth	Socioeconomic Data and Applications Center [38]
Maternal education
Average minimum temperature	WorldClim [32]
Average maximum temperature

Results

Variogram analysis

Fig 3 shows that there is significant spatial autocorrelation in the observed MF prevalence data, although spatial autocorrelation starts to decay beyond 250 km. In contrast, there was limited spatial autocorrelation for CMFL, even at shorter distances.

Fig 3

Variogram plot showing the spatial autocorrelation in observed onchocerciasis infection data.

A. MF prevalence B. CMFL. The empirical variogram is represented by the black dots; the theoretical variogram is represented by the solid black line.

Variable importance analysis

Fig 4 is a plot of the percentage increment in mean square error computed for the final suite of variables included in the trained QRF model. The most important predictors were proximity to permanent rivers and streams, proximity to the coast, altitude and household wealth.

Fig 4

Variable importance for onchocerciasis infection in the trained quantile regression forest (QRF) model.

Marginal effects plots of covariates included in the model

Fig 5 indicate that the probability of onchocerciasis infections decreases with increasing distance from rivers and streams and the coast. Also, higher household wealth and better housing type (houses built with modern materials) were positively correlated with onchocerciasis infection. Travel time to major cities had little effect on onchocerciasis infection.

Fig 5

Marginal effects plots for covariates included in the QRF model.

The Y-axis is the response (probability of onchocerciasis infection) and the X-axis is the covariate values.

MF prevalence and CMFL circa 1991 and 2016

The overall mean observed MF prevalence decreased from 43.2% circa 1991 to 7.6% circa 2016, and the observed mean CMFL decreased from 13.8 MF/snip circa 1991 to 0.6 MF/snip circa 2016. The relatively greater decrease in CMFL means that ivermectin had a greater effect on infection intensity than on prevalence.

Maps presented in Fig 6 show predicted MF prevalence circa 1991 and circa 2016. These maps indicate that prior to mass treatment with ivermectin, all regions in Côte d’Ivoire were endemic for onchocerciasis, although areas in the north were hypo-endemic and south were generally meso- to hyper-endemic. Following the use of ivermectin, endemicity levels decreased significantly throughout the country. However, infection persists in focal areas in the south and central regions. The overall predicted mean MF prevalence decreased from 51.8% circa 1991 to 3.9% circa 2016. The mean Pearson’s correlation between observed and predicted MF prevalence was 0.687 for model trained with surveys conducted from 1975 to 1991and 0.632 for model trained with surveys conducted from 1992 to 2016.

Fig 6

Predicted mean MF prevalence.

Plot on the left and right are predictions for onchocerciasis endemicity for circa 1991 and circa 2016, respectively.

Predicted CMFL data shown in Fig 7 also indicate that onchocerciasis infection was widespread in Côte d’Ivoire in the pre-ivermectin, although values were higher in the northern and southern regions of the country. The overall predicted mean CMFL was 10.1 MF/snip circa 1991. This decreased dramatically following the use of ivermectin to 0.1 MF/snip circa 2016, when CMFL predictions show no infection in most pixels in the map, although ongoing infections were predicted in areas in the south and west-central regions. The mean Pearson’s correlation between observed and predicted CMFL in surveys conducted prior to ivermectin MDA/CDTi was 0.664, and 0.583 for surveys conducted during CDTi.

Fig 7

Predicted mean CMFL.

Plot on the left and right are predictions for onchocerciasis endemicity for circa 1991 and circa 2016, respectively.

The RMSE and R² values for the trained QRF model were 1.14 and 0.62 respectively for the model trained with onchocerciasis data from 1975 to 1991. The RMSE and R² values for the model trained with infection data collected from 1992 to 2016 was 1.28 and 0.57, respectively.

Discussion

In this study, we have produced continuous maps of onchocerciasis MF prevalence and CMFL in Côte d’Ivoire circa 1991 and 2016. We have also identified important socio-economic and environmental variables that correlate with onchocerciasis infection. Our maps illustrate significant reductions in both MF prevalence and CMFL that occurred after the introduction of ivermectin in 1992. We have also identified potential remaining infection hotspots, mainly in the southern region of Côte d’Ivoire. We believe that the data presented in this study are useful for understanding changes in the spatial distribution of onchocerciasis in Côte d’Ivoire over time.

The predictive accuracy of machine learning models is usually assessed by exploring the ability of models to correctly predict results for an independent dataset [51]. As there was no independent dataset available, we trained the QRF model on a random sample of 70% of the total data and predicted results for the held-out 30%. Predictive accuracy was presented as R², which is the percentage of variation explained by the covariates in our model. The model trained with data collected before and during the ivermectin era had R² values of 62% and 57%, respectively. Although these values are reasonably high compared to those presented in other modelling exercises [28], we believe that our model would have performed even better if we had been able to incorporate more detailed information on onchocerciasis morbidity and interventions (vector control and ivermectin treatment) as predictors in our model. This is also shown in our uncertainty maps, presented as standard deviations around per-grid square predictions (S1 Fig). Where standard deviations are high, confidence in model predictions is lower than in locations where the standard deviation is low. Here the contrasting uncertainty maps in the two periods studied may reflect missing intervention data; vector control (mostly in the northern and central regions from 1975 to 1991) and ivermectin treatment in the southern region that started in 1992.

Migration of humans and the drifting of black flies by monsoon winds from far distant hotspots with ongoing infection can also lead to resurgence of infection in areas previously cleared [52–55]. Our variogram analysis shows that the spatial autocorrelation in the observed MF prevalence data has a range of ~250 km (Fig 3A). This is consistent with previous findings [25] and with reports that black flies can travel and infect humans several hundred kilometers from their breeding sites [55–57], although biting density drops to 10% of its highest levels in areas that are more than 5 km from breeding sites [58]. However, if infections levels in humans are reduced to elimination thresholds due to effective intervention measures, resurgence is unlikely regardless of local increases in vector populations.

Although machine learning algorithms are increasingly being used for spatial modelling, they sometimes fail to account for the spatial structure of the outcome of interest. In order to correct for this, in addition to using spatially referenced predictors, we included the geographical coordinates of survey villages in our dataset as covariates in our final model. This adds spatial structure and improves model predictions, as previously reported [49]. We extended the methods of Eneanya et al. [27] to carry out a robust, step-wise selection procedure to identify the best suite of uncorrelated explanatory variables. Multicollinearity often arises in statistical models, and it can lead to unstable estimates of the variance of regression coefficients [43].

Remotely-sensed covariates that are associated with the modeled outcome aid in defining the natural geographical limits of the prediction. This improves the ability of the model to explain the variability in predicted outcomes and to account for further spatial structure during the modelling process. We considered land areas designated as cropland, the presence of irrigation infrastructure, type of housing, and household wealth as covariates for infection. These have not been used as predictors in previous spatial models for onchocerciasis. However, we view these as key determinants of the potential distribution of onchocerciasis for the following reasons; i) Agricultural communities are historically known to be highly endemic for onchocerciasis. Farms in these communities are often situated near rivers that provide water for crops. As black flies bite outdoors, farmers and others working outdoors have an increased risk of infection. ii) Wealthier households are more likely to have land for farming.

Our model predicts meso- to hyperendemic hotspots in areas of northern Côte d’Ivoire prior to ivermectin treatment. Previous predictions for these areas [25] were generally higher than ours. In their work, O’Hanlon et al. only considered villages in the OCP region that were intervention naïve, whereas villages in our dataset for this area in Côte d’Ivoire had received multiple years of intervention by vector control. Our model also predicted high-level endemicity in the southern part of the country, and those predictions are consistent with results from studies such as that of Dadzie et al. [59]. The authors documented high onchocerciasis prevalence in 11 first-line villages in the lower areas near the rivers Bandama and Comoe. Similarly, known hyperendemicity in the Cavally River valley in the western region [60] was also captured in our model.

Our study confirmed that onchocerciasis infection is negatively associated with increasing distance from rivers. Elevation and proximity to breeding sites were reported to be important predictors of onchocerciasis endemicity and severity in Cameroon [61] and Venezuela [62]. Temperature, an important predictor in our model, is linked to development of the black fly vector and to development of parasite larvae in the vector [63]. Proximity to major cities was the least important variable in our model, perhaps due to the extensive geographical coverage of the survey dataset used to construct the models.

Our study predicted high onchocerciasis prevalence for much of southern Côte d’Ivoire. The native vegetation type in this area is rain forest, and high prevalences have been observed for forest-type onchocerciasis in Côte d’Ivoire [64] and elsewhere [65]. Furthermore, data presented by Adjami et al. suggest that the savanna-dwelling species S. damnosum s.s. appear to thrive in areas along the middle Bandama river leading to an expansion of savanna type O. volvulus [66]. This may be as a result of climatic and anthropogenic changes that have resulted in large scale deforestation in West Africa.

Ground-truth data for villages to the east and west of Abidjan indicated meso- or hyperendemic. A recent meta-analysis of onchocerciasis-induced epilepsy in West Africa recorded cases in the south of Côte d’Ivoire [67], in keeping with our model predictions for this area. As there are no ground-truth data for greater Abidjan city, our model made predictions based on these nearby villages and corresponding environmental limits of the covariates. This explains why our maps predict that areas in greater Abidjan were hyperendemic prior to treatment with ivermectin. Therefore, our model may be missing some important covariates that are peculiar to the commercial capital city of Abidjan. As it is known that onchocerciasis is not endemic in Abidjan, we have excluded Abidjan from our final prediction maps. However, it might be useful to perform field surveys to verify the absence of onchocerciasis infection and onchocerciasis-induced epilepsy in the Greater Abidjan area.

Data used to build our predictive spatial models were from surveys conducted by the National Onchocerciasis Control Programme. One limitation of these data is that surveys were not conducted at random; thus, villages surveyed may be skewed towards areas with known high endemicity. In addition, the original programme areas were selected using entomological knowledge about vector breeding sites and information on the distribution of clinical onchocerciasis in villages. Despite these limitations, data from the National OCP comprises an extensive source of standardized parasitological survey data that cover large areas within Côte d’Ivoire.

In conclusion, despite the disruption of CDTi due to civil unrest (especially between 2002 and 2007), our results clearly show significant reductions in onchocerciasis prevalence in Côte d’Ivoire after the scale up of mass treatment with ivermectin. Our identification of potential foci with ongoing infection may help control programmes target intervention with ivermectin (with or without vector control) or more frequent ivermectin distribution in areas where infection persists despite adequate ivermectin coverage. Such focused efforts are likely to be very important in the late stages of the country’s onchocerciasis elimination programme. Although this study focused on Côte d’Ivoire, this approach may be useful for identifying endemic areas and targeting interventions to eliminate onchocerciasis in other African countries.