Modelling the elimination of river blindness using long-term epidemiological and programmatic data from Mali and Senegal.

The onchocerciasis transmission models EPIONCHO and ONCHOSIM have been independently developed and used to explore the feasibility of eliminating onchocerciasis from Africa with mass (annual or biannual) distribution of ivermectin within the timeframes proposed by the World Health Organization (WHO) and endorsed by the 2012 London Declaration on Neglected Tropical Diseases (i.e. by 2020/2025). Based on the findings of our previous model comparison, we implemented technical refinements and tested the projections of EPIONCHO and ONCHOSIM against long-term epidemiological data from two West African transmission foci in Mali and Senegal where the observed prevalence of infection was brought to zero circa 2007-2009 after 15-17 years of mass ivermectin treatment. We simulated these interventions using programmatic information on the frequency and coverage of mass treatments and trained the model projections using longitudinal parasitological data from 27 communities, evaluating the projected outcome of elimination (local parasite extinction) or resurgence. We found that EPIONCHO and ONCHOSIM captured adequately the epidemiological trends during mass treatment but that resurgence, while never predicted by ONCHOSIM, was predicted by EPIONCHO in some communities with the highest (inferred) vector biting rates and associated pre-intervention endemicities. Resurgence can be extremely protracted such that low (microfilarial) prevalence between 1% and 5% can be maintained for 3-5 years before manifesting more prominently. We highlight that post-treatment and post-elimination surveillance protocols must be implemented for long enough and with high enough sensitivity to detect possible residual latent infections potentially indicative of resurgence. We also discuss uncertainty and differences between EPIONCHO and ONCHOSIM projections, the potential importance of vector control in high-transmission settings as a complementary intervention strategy, and the short remaining timeline for African countries to be ready to stop treatment safely and begin surveillance in order to meet the impending 2020/2025 elimination targets.

Introduction

Human onchocerciasis or river blindness is caused by the filarial nematode Onchocerca volvulus and is earmarked for elimination by the World Health Organization (WHO) as articulated by the 2012 Roadmap (WHO, 2012) and the London Declaration on Neglected Tropical Diseases (2012). The principal strategy to achieve elimination is mass drug administration (MDA) with ivermectin. Ivermectin kills the skin-dwelling microfilariae (mf) that are the progeny of adult O. volvulus and are infectious to biting blackfly species vectors. Ivermectin may also kill and/or sterilize adult worms (Gardon et al., 2002). Multiple rounds of mass treatment are effective in lowering the prevalence and intensity of onchocerciasis and – if given for long enough at high enough coverage – can lead to the interruption of transmission and elimination of the infection. In Latin America, (biannual) MDA has successfully eliminated onchocerciasis from Colombia (West et al., 2012), Ecuador (Lovato et al., 2014), northern Venezuela (Convit et al., 2013) and Mexico (Rodríguez-Pérez et al., 2015), with Guatemala awaiting certification. Good progress towards elimination has also been made in Africa (Tekle et al., 2016) which bears 99% of the onchocerciasis burden, with notable successes in regions of Mali, Senegal (Diawara et al., 2009, Traore et al., 2012), Nigeria (Tekle et al., 2012), Sudan (Higazi et al., 2013) and eastern Uganda (Katabarwa et al., 2014) but also conspicuous regions of ongoing transmission despite years of intervention in Ghana (Lamberton et al., 2015), Cameroon (Katabarwa et al., 2013a, Wanji et al., 2015, Eisenbarth et al., 2016) and northwestern Uganda (Katabarwa et al., 2013b), and evidence of recrudescence in Burkina Faso (Koala et al., 2017).

EPIONCHO and ONCHOSIM are mathematical transmission models of human onchocerciasis that have been used to evaluate the effectiveness of different intervention strategies in reaching the 2020/2025 elimination targets for onchocerciasis (WHO, 2012, APOC, 2012). Both models have been used to predict the time to eliminate onchocerciasis under annual or biannual MDA (Coffeng et al., 2014a, Turner et al., 2014), and have been used in collaboration with the African Programme for Onchocerciais Control (APOC) to evaluate how progress towards elimination can be accelerated by implementing alternative treatment strategies (WHO, 2015). ONCHOSIM has been used to inform the Onchocerciasis Control Programme in West Africa (OCP) and the African Programme for Onchocerciasis Control on expected trends of infection during vector control (Plaisier et al., 1991), mass ivermectin treatment (Winnen et al., 2002, Coffeng et al., 2014a, Tekle et al., 2016) or their combination (Plaisier et al., 1997), the health impact of the interventions (Coffeng et al., 2014b), and the expected time to elimination (Kim et al., 2015).

More recently, the developers of EPIONCHO and ONCHOSIM have started to compare and refine their models, with the objective of reaching consensus on the feasibility of and time horizons for elimination, and optimum interventions for achieving the WHO goals in Africa (Stolk et al., 2015). The first formal comparison, which attempted to ‘dock’ the two models by making parameter assumptions as equivalent as possible, revealed some important differences between the models’ structural assumptions and resulting outputs. Both models demonstrated a benefit of biannual versus annual treatment but differed in their predicted time to elimination, particularly when using transmission breakpoints as endpoints rather than operational (microfilarial) prevalence thresholds such as those provisionally proposed by APOC (2010). One of the conclusions of this first comparison was that the use of a single prevalence threshold across different endemicity (transmission intensity) settings may not be appropriate and needs further investigation (Stolk et al., 2015).

The epidemiology of onchocerciasis has changed as a result of the introduction of ivermectin MDA, which means that EPIONCHO and ONCHOSIM are now being used in transmission contexts far from those in which they were originally developed and parameterized (Basáñez and Boussinesq, 1999, Plaisier et al., 1991, Plaisier et al., 1995). Hence, to maintain confidence in projections, these models must be tested and validated against data collected before, during and after MDA, reflecting the spectrum of transmission conditions, from endemic infection to transmission interruption and elimination (Basáñez et al., 2012a, Basáñez et al., 2012b). Given the prolonged nature of onchocerciasis interventions, capturing and recording such data is challenging and resource-intensive and would probably have been impossible without the logistical and technical support of the OCP and APOC.

Here we test the EPIONCHO and ONCHOSIM models against epidemiological data collected over two decades from 27 sentinel communities in two onchocerciasis foci in Mali and Senegal. In these foci, MDA with ivermectin alone has successfully brought the prevalence of infective larvae (in blackflies) and of mf (in humans) – detectable by skin snips – to zero, indicative of elimination (Diawara et al., 2009, Traore et al., 2012). In particular, prevalence of skin mf reached the previously proposed operational thresholds set by APOC of <5% in all surveyed villages and <1% in 90% of those surveyed (APOC, 2010). We use programmatic data on the frequency and coverage of MDA with ivermectin to simulate these interventions from endemic baseline, circa 1987, through to cessation of treatment, circa 2006, and then beyond, projecting forwards to 2020. We train the model projections using subsets of epidemiological data comprising community-specific estimates of microfilarial prevalence collected throughout the intervention. We determine whether the models predict sustained elimination or resurgence of infection in the post-treatment and post-elimination periods (WHO, 2016) and explore how these predictions change when using increasing amounts of epidemiological training data, from pre-intervention data only to multiple longitudinal data points per community.

Models and methods

EPIONCHO

EPIONCHO is a deterministic onchocerciasis transmission model that uses partial differential equations to describe changes (with respect to time and host age) in mean number of fertile and non-fertile female adult worms per host, mean number of mf per milligram (mg) of skin in humans and mean number of larvae per simuliid (blackfly) vector. Briefly, the model is based on a prototype presented by Basáñez and Boussinesq (1999), and extended to include age and sex structure of the human population (Filipe et al., 2005); the temporal dynamics of mf following ivermectin treatment (Basáñez et al., 2008), and increased programmatic realism related to patterns of treatment coverage and systematic non-adherence (Turner et al., 2013). The model allows for age- and sex-specific patterns of exposure to blackfly bites (Filipe et al., 2005) and variation in adherence to treatment. The latter is modelled by partitioning the population into four groups, namely, a full adherence group that takes treatment every round; two semi-adherent groups that take treatment every other round alternately, and a systematically non-adherent group that never takes treatment. Treatment with ivermectin is assumed to kill 98–99% mf, temporarily sterilize adult female worms and cumulatively reduce their capacity to produce mf (Basáñez et al., 2008, Turner et al., 2013). EPIONCHO has been used to address public health policy questions including an economic evaluation of implementing community-directed treatment with ivermectin (CDTI) biannually compared to annually (Turner et al., 2014); the potential epidemiological impact of MDA with moxidectin (Turner et al., 2015a), and the future deployment of an onchocerciasis vaccine (Turner et al., 2015b).

The latest version of EPIONCHO used here incorporates a number of refinements from the version described in Turner et al. (2013) and reviewed by Basáñez et al. (2016). Some of these are motivated by the results of our recent comparison with ONCHOSIM (Stolk et al., 2015), and others by our improved understanding of the population biology of O. volvulus. Modifications include: (1) the introduction of a latent period for the development of O. volvulus within the blackfly vector, implemented by dividing the larval state L into L1, L2 and L3 compartments (Basáñez et al., 2007; Basáñez et al., 2009); (2) a reparameterization of the microfilarial ‘uptake’ curve (describing parasite establishment within the vector) to reflect that an overdispersed distribution of parasites among human hosts affects the overall severity of density dependence (Churcher et al., 2005); (3) a new functional relationship between the prevalence and intensity of mf that reflects the underlying distribution of adult worms, incorporates the sensitivity of skin snips (Bottomley et al., 2016) and better captures the initial ‘bounce back’ of microfilarial prevalence during the first few rounds of ivermectin MDA (Plaisier et al., 1995); and (4) a more realistic distribution of adult parasite survival times (previously exponential), implemented by introducing multiple nominal ‘age’ compartments for adult O. volvulus and parameterised using published estimates of life expectancy as well as the ‘tail’ of the distribution of survival times (Plaisier et al., 1991). Full details of these refinements and the complete description of this latest version of EPIONCHO are given in the Supplementary information, 1.1 EPIONCHO and a list of parameter definitions, descriptions and values in Supplementary information, Table S1.

ONCHOSIM

ONCHOSIM is a long-established microsimulation model of onchocerciasis transmission (Habbema et al., 1996, Plaisier et al., 1990). Briefly, the model considers a human population, typically consisting of approximately 400 individuals, with a demographic composition that changes stochastically over time. The life-histories of individual male and female adult parasites and the populations of mf within individual human hosts are tracked and transmission is mediated by a population of blackflies (the parasite-within-the-vector component is treated deterministically). The survival times of adult worms follow a Weibull distribution. The biting rate (number of bites per person per unit time) is seasonal, varies among hosts according to an individual-specific exposure variable (assigned at birth), and changes systematically with host age and sex in a manner similar to that of EPIONCHO. Treatment in ONCHOSIM is modelled in a broadly similar manner to EPIONCHO, albeit (100%) mf are killed instantaneously in ONCHOSIM rather than over 1–2 months as in EPIONCHO and the rates at which female worms recover their fertility are subtly different between the two models (Plaisier et al., 1995; Basáñez et al., 2008, Coffeng et al., 2014a, Basáñez et al., 2016). Human hosts, excluding children aged <5 years and a fraction of women of reproductive age, participate in each round of MDA with a probability that depends on age and sex, and a parameter governing their lifelong adherence to MDA. ONCHOSIM has been used to inform and support OCP (Plaisier et al., 1997, Winnen et al., 2002) and APOC (Coffeng et al., 2014a, Coffeng et al., 2014b) intervention policies and more recently has been used to evaluate the progress of APOC projects targeting elimination (Tekle et al., 2016). Technical details of the ONCHOSIM model can be found in Supplementary information, 1.2 ONCHOSIM, Coffeng et al. (2014a) and Basáñez et al. (2016). The ONCHOSIM parameter definitions, descriptions and values are given in Supplementary information, Table S2.

Epidemiological and programmatic data

We use long-term parasitological and programmatic data from the River Bakoye and River Gambia foci in Mali and Senegal where onchocerciasis has been declared eliminated after 15–17 years of, respectively, annual or biannual MDA (Diawara et al., 2009, Traore et al., 2012) (Fig. 1). We do not use data from the cross-border River Faleme focus. The predominant blackfly vector species in both foci is Simulium sirbanum (a savannah member of the S. damnosum s.l. complex; Boakye et al., 1998). Sentinel communities were initially selected by the OCP circa 1987 in areas of the Western Extension where vector control had never been implemented and therefore provided a unique opportunity to evaluate the impact of ivermectin as the sole measure of control. Communities were routinely monitored by the OCP from 1987 to 2002 and from 2006 by a team studying the feasibility of eliminating onchocerciasis with MDA alone (Diawara et al., 2009).

Fig. 1

Endemic microfilarial prevalence (before intervention) in selected sentinel communities in (a) the River Bakoye focus, Mali and (b) the River Gambia focus, Senegal circa 1987. Endemicity category is defined by a microfilarial prevalence in people aged ≥5 years of <40% (hypoendemic), 40%–59% (mesoendemic) and ≥60% (hyperendemic) and coloured sequentially from yellow to red. Error bars denote 95% confidence intervals. Community names are anonymised.

The epidemiological data comprise community-specific sequences of microfilarial prevalence estimates. These estimates are based on the detection, by skin snips, of mf in all those aged 1 year and older using standard OCP protocols (the observation of at least one microfilaria in 2 skins snips incubated in distilled water for at least 30 min, but for 24 h if negative after this time) and standardized by age and sex using the OCP reference population (Moreau et al., 1978). Because young children tend to be underrepresented in community skin-snip surveys, we assumed that the data better represented the prevalence in those aged ≥5 years rather than ≥1year, thereby harmonizing the sampled demographic group with usual EPIONCHO and ONCHOSIM outputs (Stolk et al., 2015). We excluded from our analysis communities where only pre- or post-intervention data were available (Table 1).

Table 1

Summary of epidemiological data on community prevalence of microfilariae.

Focus, country	Number of communities	Endemic microfilarial prevalence before intervention (range)	Sequence length of microfilarial prevalence training data (range)
River Bakoye, Mali	13	31%–70%	3–5
River Gambia, Senegal	14	36%–82%	4–11

The programmatic data comprise estimates, within a range, of the overall therapeutic coverage (proportion of the total population treated) with ivermectin within each river basin for the periods 1988–1991, 1992–1996 and 1998–2006 (Table 2). Ivermectin was first distributed by OCP mobile teams before CDTI – which empowered communities to appoint local community-drug distributors who took responsibility for distributing ivermectin – was introduced in 1997, following the inception of APOC in 1995 (Boatin, 2008). There is a notable drop in coverage in the first year of CDTI (Table 2) associated with logistical challenges in changing the distribution system (Diawara et al., 2009).

Table 2

Summary of programmatic data on community coverage with ivermectin.

Focus, country	Time period	Treatment frequency	Treatment coverage
River Bakoye, Mali	1989–1991	Annual	59%–62%
	1992–1996	Annual	75%–78%
	1997	Annual	∼38%b
	1998–2006	Annual	73%–83%
River Gambia, Senegal	1988–1991	Annual/biannuala	64%–69%
	1992–1996	Biannual	76%–77%
	1997	Biannual	∼38%b
	1998–2006	Biannual	77%–81%

Annual 1988–1990; biannual thereafter.

First year of community-directed treatment with ivermectin, coverage was low (Diawara et al., 2009) and here is assumed to be half that of the preceding period.

Simulating interventions and estimating community-specific biting rates

We simulated MDA interventions using the (focus-specific) programmatic data on frequency and coverage summarized in Table 2. Both EPIONCHO and ONCHOSIM are initialized using the annual biting rate (ABR), defined here as the average number of blackfly bites received per person per year. Hence, to run community-specific simulations, we required an estimate of the ABR in each community. We estimated ABRs by a maximum likelihood approach (Supplementary information, 1.4 Estimating community-specific biting rates) using subsets of the microfilarial prevalence time series from each community. We refer to these subsets as training datasets, each comprising a variable sequence length of longitudinal prevalence estimates from each community. The smallest training dataset comprised pre-intervention data only and the largest included all data from all communities. We iteratively re-estimated ABRs by incrementally increasing the length of the longitudinal training datasets, evaluating how increasing amounts of training data altered and influenced the model projections. Data on the exact time of treatment (unlike the time of epidemiological surveys) were not available and hence, for the purposes of estimating ABRs, we assumed that surveys were undertaken exactly 11 or 5 months after the last treatment for, respectively, the annual and biannual distributions in the River Bakoye and River Gambia foci (Supplementary information, Table S3).

Parametric uncertainty

The estimation of ABRs depends on its relationship with infection which, in turn, depends on structural assumptions and parameter values of the model. In ONCHOSIM, we estimated ABRs, accounting for seasonality in transmission, while all other parameters were fixed at their default values (Coffeng et al., 2014a, Basáñez et al., 2016). In EPIONCHO, we explicitly considered that there is uncertainty in some essential parameters: multiple unique parameter sets are all concordant with an independent dataset comprising coupled pre-intervention community-level prevalence and ABR estimates from Cameroon (Renz and Wenk, 1987, Basáñez and Boussinesq, 1999) (Fig. 2a). We identified these parameter sets using a sampling importance resampling (SIR) approach (Gambhir et al., 2015), in which we sampled from the posterior distribution of parameter values given uniform priors defined by the ranges of published estimates (Supplementary information, Table S1). Details of this approach are given in Supplementary information, 1.3 Parametric uncertainty.

Fig. 2

Modelled relationship between annual biting rate (ABR) of blackfly vectors and endemic microfilarial prevalence using (a) EPIONCHO and (b) ONCHOSIM. In panel (a), the coupled ABR-prevalence data are from 9 communities in northern Cameroon (Basáñez and Boussinesq, 1999) where each ABR was measured as an average from multiple years and locations within and around each community, weighted by the proportion of time community residents spent at these locations (Renz and Wenk, 1987). Each thin line corresponds to an EPIONCHO parameter set identified by the sampling importance resampling (SIR) procedure described in Supplementary information, 1.3 Parametric uncertainty. These are coloured sequentially from yellow to red in accordance with endemicity category as defined by a microfilarial prevalence in people aged ≥5 years of <40% (hypoendemic), 40%–59% (mesoendemic) and ≥60% (hyperendemic). The thick black line corresponds to the parameter set that achieved the highest likelihood. Inset is the same graph with the x-axis transformed to a logarithmic scale. In panel (b), the thin lines correspond to stochastic realizations of the ONCHOSIM default parameter set (Coffeng et al., 2014a), coloured sequentially according to endemicity category, and the thick black line is the median of 500 simulations. The coupled ABR-prevalence data are not shown in panel (b) because ONCHOSIM has not been re-fitted to these data.

Selection of a maximum a posteriori parameter set

The SIR approach applied to EPIONCHO involves defining a prior distribution for discrete parameter sets, thereby supporting a wide range of model projections for an estimated ABR. We calculated the posterior probability of each community projection as the likelihood of the combined training dataset from all communities (i.e. a single combined dataset for each focus) weighted by the parameter set (prior) probability calculated from the SIR procedure. For a more direct comparison with ONCHOSIM – in which we only consider a single default parameter set – we selected the parameter set that maximised the posterior probability, akin to a maximum a posteriori (MAP) approach. Details are given in Supplementary information, 1.5 Selecting maximum a posteriori (MAP) parameter sets. We highlight that, unlike the estimation of ABRs, undertaken on a community basis, identification of MAP parameter sets was done at the regional level (i.e. including both foci); fundamental biological and population dynamics parameters were assumed constant among communities and between foci (but ABRs were permitted to vary).

Elimination criterion

We evaluated the outcome of each community-based simulation as sustained elimination or resurgence according to the MAP parameter sets for EPIONCHO and the default parameter set for ONCHOSIM. We defined sustained elimination as local parasite extinction, indicated by the parasite population tending terminally to zero. If elimination is not sustained, resurgence occurs. In EPIONCHO, sustained elimination arises by crossing the transmission breakpoint. This arises because of the necessity for female O. volvulus to mate to produce mf (modelled by the so-called mating probability (May, 1977), Supplementary information, 1.1 EPIONCHO). The population density below this threshold becomes unsustainably low and tends to extinction (evaluated numerically as when the mean number of adult worms immediately after the last treatment is greater than the corresponding value 20 years later, i.e. indicative of terminal decline). The mating probability is governed by the assumed sex ratio, the parasite mating system (assumed to be polygamous, whereby one male can mate all females within a single host; Schulz-Key and Karam, 1986, Hildebrandt et al., 2012) and the degree of (adult) worm overdispersion among hosts which is included as an uncertain parameter in our SIR approach (Supplementary information, 1.3 Parametric uncertainty). In ONCHOSIM the breakpoint (also determined by the assumed polygamous mating system and the degree of worm overdispersion) operates in combination with chance elimination of the parasite population caused by demographic stochasticity in the human host population, a non-negligible phenomenon because of the relatively small simulated human population size of 400–440 individuals (see Supplementary information, Table S2). Neither EPIONCHO nor ONCHOSIM consider possible reintroduction of infection (by blackflies or humans) from proximate foci where transmission maybe ongoing.

Coverage and adherence

We repeated the simulations, re-estimating ABRs (and for EPIONCHO, identifying MAP parameter sets) and evaluating projections of sustained elimination or resurgence using different assumed values of the estimated treatment coverage, which was given as within the ranges shown in Table 2. Specifically, we ran repeated simulations using the lower, median and upper estimates of coverage within each timeframe (Table 2). The degree of non-adherence to MDA was unknown and was set to 5% for both models, in accordance with previous simulations (Stolk et al., 2015).

Results

Infection dynamics towards elimination

We present in Fig. 3 the observed microfilarial prevalence and the EPIONCHO- and ONCHOSIM-modelled dynamics in the River Bakoye focus in Mali and in the River Gambia focus in Senegal. The simulations shown are those with estimated community-specific ABRs that maximise the likelihood of the community training datasets using pre-intervention microfilarial prevalence data only (Fig. 3a, b, e and f) or the complete longitudinal training dataset sequence from each community (Fig. 3c, d, g and h; the additional fits to each sequence length of longitudinal data are shown in the Supplementary information, Figs. S1 and S2; the maximum likelihoods of the data given the estimated community-specific ABRs are presented in Supplementary information, Table S4). We present EPIONCHO projections for MAP parameter sets for a viz-a-viz comparison with the fixed parameter set used by ONCHOSIM. In the River Bakoye focus, under annual ivermectin, EPIONCHO predicts a slightly faster initial decline in prevalence than ONCHOSIM. This is not apparent in the River Gambia focus under biannual treatment. Nevertheless, both models adequately capture the long-term broad trends in the data from pre-intervention endemic prevalence, through the treatment phase towards elimination. It is noteworthy that, in general, the estimated community-specific ABRs (and the corresponding model-derived pre-intervention endemicity) decrease with increasing longitudinal training dataset sequences. This is particularly apparent for the EPIONCHO projections in the River Bakoye focus where, when using the complete training dataset, no communities are predicted as being hyperendemic (i.e., no pre-intervention microfilarial prevalence ≥60%, Fig. 3c) despite five communities being putatively classified as such (although only one statistically significantly so) using pre-intervention data only (Fig. 1).

Fig. 3

Observed and modelled dynamics of microfilarial prevalence in 13 communities from the River Bakoye focus, Mali (panels a–d) and in 14 communities from the River Gambia focus, Senegal (panels e–h). Panels on the left (a, c, e and g) and on the right (b, d, f and h) show EPIONCHO and ONCHOSIM projections, respectively. The thin lines correspond to community-specific simulations using maximum likelihood estimates of the community-specific annual biting rates (ABRs) and either the maximum a posteriori (MAP) parameter set (EPIONCHO) or the default parameter set (ONCHOSIM). The estimated ABRs and MAP parameter sets are derived using the pre-intervention microfilarial prevalence data only (a, b, e and f) or using the complete longitudinal sequence for each community (c, d, g and h). For brevity, the dynamics predicted from estimates using one or two interim time points are not shown. For ONCHOSIM there are many stochastic projections for each community projection; for EPIONCHO there is a single deterministic projection for each community, corresponding to the MAP parameter set. The thick solid lines show the median dynamics by endemicity category as categorised by a model-derived pre-intervention microfilarial prevalence in people aged ≥5 years of <40% (hypoendemic), 40%–59% (mesoendemic) and ≥60% (hyperendemic) and coloured sequentially from yellow to red. In panel c the estimated ABRs (from EPIONCHO, using the complete longitudinal data sequences) indicated all communities in the River Bakoye focus were either hypoendemic or mesoendemic. Panel insets show the period between 2010 and 2020 using a transformed y-axis for a better visual appraisal of the model projections compared to the data close to zero.

Elimination or resurgence

We present in Fig. 4 the number and percentage of communities where the EPIONCHO-predicted outcome of MDA was sustained elimination (local parasite extinction) using incrementally increasing subsets of the longitudinal training datasets (Table 1) and different assumed levels of treatment coverage within the range of reported estimates (Table 2). The percentage of (stochastic) ONCHOSIM projections resulting in sustained elimination was greater than 99% for all communities, irrespective of the length of the training datasets. By contrast, EPIONCHO projections were variable. Longer training datasets resulted in more communities with a prediction of sustained elimination in the River Bakoye focus. In the River Gambia focus the length of the training data made little difference to the outcome of sustained elimination or resurgence. It is noteworthy that in some instances, there are more communities with a prediction of elimination when lower levels of treatment coverage were assumed (Fig. 4b). This reflects parametric uncertainty in the model projections since different MAP parameter sets could be selected when estimating community-specific ABRs using different training data subsets and different assumed levels of treatment coverage.

Fig. 4

Number and percentage of community-specific EPIONCHO projections resulting in elimination in (a) the River Bakoye focus, Mali and (b) the River Gambia focus, Senegal. Each grid element indicates the number and percentage of communities with a projected outcome of sustained elimination when training the model projection using different sequence lengths of the community-specific microfilarial prevalence estimates (x-axis, see Table 1) and assuming different levels of coverage (y-axis, see Table 2). The outcomes are coloured sequentially from red (50% of communities with predicted elimination) to green (100% of communities with predicted elimination).

Generally, in EPIONCHO the outcome of sustained elimination is sensitive to the model-derived endemicity, as driven by the estimated community-specific ABRs. Fig. 5 shows that – regardless of the length of the training dataset or the assumed level of coverage – sustained elimination was rarely predicted in hyperendemic communities (where the pre-intervention endemic microfilarial prevalence was estimated to be ≥60%; a slightly higher percentage of simulations with pre-intervention hyperendemicity occurred in the River Gambia focus, Senegal compared with the River Bakoye focus, Mali because of the biannual treatment strategy). When elimination is predicted, the rate of resurgence is variable but sometimes very protracted, depending on community endemicity (and the corresponding ABR). We illustrate this in Fig. 6 by showing the EPIONCHO-predicted microfilarial prevalence levels in 2007, 2009 and 2012 – one, three and five years after the last round of treatment in 2006 – among communities where either sustained elimination (microfilarial prevalence tending to zero) or resurgence was indicated. Many of the simulations predict a microfilarial prevalence between 1% and 5% at these times in communities where infection is (slowly) resurgent.

Fig. 5

The percentage of communities with an EPIONCHO-projected outcome of sustained elimination in (a) 13 communities in the River Bakoye focus, Mali and (b) 14 communities in the River Gambia focus, Senegal. Endemicity is categorized according to the model-derived pre-intervention microfilarial prevalence estimates (using the estimated community-specific annual biting rates) in people aged ≥5 years of <40% (hypoendemic), 40%–59% (mesoendemi) and ≥60% (hyperendemic) and coloured sequentially from yellow to red. Each bar indicates the combined results from estimating the community annual biting rate (and the subsequent pre-intervention endemicity) using different sequence lengths of the community training datasets (Table 1) and different assumed levels of treatment coverage (Table 2). Hence, in any simulation where the pre-intervention endemicity was estimated as hyperendemic, the chance of elimination was low, irrespective of the annual (a) or biannual (b) treatment strategy.

Fig. 6

EPIONCHO predictions of microfilarial prevalence in 2007, 2009 and 2012, one, three and five years after the last round of treatment in 2006, among communities in (a–c) the River Bakoye focus, Mali and (d–f) the River Gambia focus, Senegal. Endemicity is categorized according to model-derived pre-intervention microfilarial prevalence estimates (using the estimated community-specific annual biting rates) in people aged ≥5 years of <40% (hypoendemic), 40%–59% (mesoendemic) and ≥60% (hyperendemic) and coloured sequentially from yellow to red. Each category includes the combined results from estimating the community annual biting rate (ABR) and the subsequent pre-intervention endemicity using different sequence lengths of the community training datasets (Table 1) and different assumed levels of treatment coverage (Table 2). The circular data points in each panel indicate simulations where resurgence is the predicted outcome (only in red hyperendemic categories, see Fig. 5) and the triangular data points (all close to 0%) are those corresponding to sustained elimination. The violin plots illustrate the density of the microfilarial prevalence predictions. The horizontal dashed lines represent nominal (operational) microfilarial prevalence thresholds of 1% and 5% (APOC, 2010).

Discussion

We have tested the EPIONCHO and ONCHOSIM onchocerciasis transmission models against long-term epidemiological data from two foci in Mali and Senegal where the (detectable) prevalence of O. volvulus mf has been brought to zero by, respectively, annual and biannual MDA with ivermectin alone. Both models capture adequately community infection dynamics through MDA towards elimination. Elimination in all 27 communities is the ubiquitous outcome predicted by ONCHOSIM, irrespective of inferred vector biting rates and uncertain coverage. In EPIONCHO, elimination is less certain, being dependent on the inferred community-specific blackfly biting rates (ABRs) and corresponding (model-derived) levels of infection endemicity. In the more highly endemic communities in the River Gambia focus, Senegal – with the highest vector biting rates – EPIONCHO indicates a substantive risk of infection resurgence, despite the biannual MDA strategy. Resurgence dynamics can be very protracted such that the microfilarial prevalence often remains between 1% and 5% for at least five years after the last treatment, before rising more substantively.

The model-predicted protracted nature of possible resurgence has important implications for the design of elimination-verification protocols, post-treatment surveillance (PTS) and post-elimination surveillance (PES). Notwithstanding the new WHO guidelines (WHO, 2016) for verifying elimination – which recommend serological surveillance of children (possibly in combination with skin snipping) and entomological monitoring during at least three years of post-treatment surveillance (PTS) (Golden et al., 2016) – hitherto in Africa, monitoring, evaluation and surveillance have still mostly been based on skin snipping to detect mf. APOC set operational prevalence thresholds for stopping treatment and starting surveillance (<5% microfilarial prevalence in all surveyed villages and <1% in 90% of those surveyed) that should be met for a three-year PTS period to confirm elimination (APOC, 2010). Indeed, the progress of 54 APOC project areas (comprising clusters of 5–20 villages) targeting these thresholds has recently been evaluated by comparing observed microfilarial prevalence versus ONCHOSIM predictions (Tekle et al., 2016), with 46/54 (85%) judged on track or exceeding prediction in their progress towards PTS. The protracted nature of some of the resurgence dynamics indicated by EPIONCHO, initially characterized by latent, difficult-to-detect infections levels before more marked increases, suggests that the proposed 3–5 years of PTS (APOC, 2010, WHO, 2016) must be extremely robust – and include optimized and efficient sampling methodologies – to be effective in identifying the early signs of infection resurgence. This also means that for African countries to confirm elimination by 2020 (WHO, 2012) or 2025 (APOC, 2012) they will need to be in a position to stop treatment safely either imminently or at the latest by 2022. For many African countries, the former is unlikely. Meeting the latter 2022 stopping-treatment date will require careful monitoring and evaluation of progress.

Critically, the resurgence dynamics projected by EPIONCHO are intimately linked with the (poor) sensitivity of skin snipping at low infection levels (Bottomley et al., 2016); in near-elimination settings many infections are likely to be missed by skin snipping alone (the distribution of mf in the skin is aggregated by contrast to the random Poisson distribution assumed by ONCHOSIM, a difference which also helps to explain the more rapid initial decline in prevalence projected by EPIONCHO compared with ONCHOSIM). That is, EPIONCHO assumes more often than ONCHOSIM that mf counts are false-negative. Thus, application of the newly recommended and more sensitive (serological) diagnostics, might enable earlier detection of resurgence. Indeed, an epidemiological survey conducted in 2014 in the River Gambia focus, including some but not all of the communities included in this analysis, found a 2.5% seroprevalence among children aged 5–9 years (Wilson et al., 2016), above the 0.1% threshold (upper 95% confidence interval) in children under 10 years recommended by the WHO for verification of elimination (WHO, 2016). Hence, by this criterion, elimination cannot be verified in the River Gambia focus and it is probable that low-level transmission continues. Resurgence – based on parasitological (mf) and entomological data – has also been reported in Burkina Faso in a focus previously thought to have eliminated the infection (Koala et al., 2017).

Notwithstanding the superior sensitivity of serological techniques, how seroprevalence profiles in children manifest in near-elimination settings and, more specifically, how the seroprevalence threshold (WHO, 2016) relates to transmission breakpoints and the likelihood of elimination, is a subject of ongoing research. A seroprevalence sub-model was recently included in ONCHOSIM (Lont et al., 2017) and is under development for EPIONCHO, in order to explore these issues more thoroughly. Early analysis suggests that it is unlikely that a single seroprevalence threshold (as an indicator of sustained elimination) will be valid across diverse transmission settings (Lont et al., 2017). This would concord with our previous analysis which indicated that in hyperendemic settings the proposed microfilarial prevalence thresholds may be too high – and hence pose a risk of resurgence if applied to stop treatment – but probably too low in lower meso- and hypoendemic settings (Stolk et al., 2015). Nevertheless, prevalence thresholds (both serological and parasitological) are extremely useful operational tools for decision-support in near-elimination settings if they can be appropriately constructed – and informed using mathematical models – to reflect the diversity of transmission settings and the vagaries of different diagnostic tools.

The greater circumspection on sustained elimination (and possible resurgence) indicated by EPIONCHO compared with ONCHOSIM reflects fundamental technical differences between the two models. One important difference is that the highly non-linear relationship between microfilarial prevalence and ABR modelled by EPIONCHO – reflecting additional density-dependent biological processes that are not included in ONCHOSIM (Basáñez et al., 2016) – supports a wide range of biting rates (Fig. 2a) and resultant transmission intensities. For example, EPIONCHO predicts an ABR of approximately 20,000 to simulate 80% microfilarial prevalence in the population aged ≥5 years. In ONCHOSIM, the corresponding ABR is approximately 14,000. Thus, EPIONCHO requires higher biting rates to simulate hyperendemicity, implying a higher risk of resurgence if some residual infection remains after intervention. Hence, the epidemiological data from some communities are associated (by EPIONCHO) with conditions highly propitious to transmission and a parasite population resilient to even long and intensive interventions. Resurgence is a probable projected outcome in many putatively hyperendemic villages, especially in the River Gambia focus in Senegal, where inferred ABRs tended to be higher than in the River Bakoye focus. In ONCHOSIM, the much less non-linear ABR-prevalence relationship (Fig. 2b) ensures more homogeneity among inferred community transmission conditions, explaining why sustained elimination is the ubiquitous modelled outcome for both foci. Furthermore, ONCHOSIM assumes less efficient transmission than EPIONCHO at low intensities (due to density-dependent parasite establishment; Basáñez et al., 2002) and a decline in mf production with increasing worm age, both making resurgence less likely. A future priority is to understand better the discrepancy between the ABR-prevalence relationship modelled by ONCHOSIM and the data collected by Renz and Wenk (1987) (where ABRs were measured as averages from multiple years and locations within and around communities, weighted by the proportion of time community residents spent at these locations) and to revisit the density dependencies included within both models.

Another reason that resurgence is more likely in EPIONCHO compared to ONCHOSIM is the large population assumption inherent to deterministic modelling approaches. Chance elimination by stochastic demographic events like the death of (all) infected hosts is not considered; elimination can only be achieved and sustained by crossing the transmission breakpoint, the theoretical population density below which the parasite population cannot sustain itself and declines terminally towards extinction. This threshold – which depends on vector abundance and the resulting ABR – is thought to be very low for onchocerciasis (Basáñez et al., 2009), and human helminthiases more generally (Anderson and May, 1991) This is because of the assumed polygamous mating system (a broadly-applied assumption in helminth transmission models, excepting schistosomes which are thought to form monogamous mating pairs; May, 1977) and the degree of (adult) worm overdispersion among hosts. A less efficient mating system and a more randomly (less aggregated) distribution of worms among hosts would make elimination by both models easier to reach. Additionally, elimination in ONCHOSIM is also influenced by chance (stochastic) demographic events.

Typically, ONCHOSIM models a relatively small (400–440) human population such that stochastic events may contribute appreciably to the elimination outcome, especially when infection prevalence is low and because parasites are aggregated among hosts; the death of the few individuals who harbour the bulk of the parasite population can lead to an extinction event. However, for the long and intensive interventions (biannual treatment frequency for the River Gambia focus in Senegal) considered here, even increasing the modelled population size to 800 individuals did not alter the projected outcome of elimination across all 27 communities. However, future analyses should be conducted to understand better the interplay between how assumptions on population size might affect the minimum required duration or frequency of treatment with ivermectin and to update current estimated timeframes for onchocerciasis elimination among different endemicity settings (Stolk et al., 2015).

Although the individual rural communities affected by onchocerciasis tend to be relatively small, whether the large-population assumption of EPIONCHO or the explicitly small (or ‘intermediate’ 800) population size modelled by ONCHOSIM is most appropriate will crucially depend on the degree to which communities are coupled by the movement of blackflies and people. Neither EPIONCHO nor ONCHOSIM consider the spatial coupling of proximate populations and especially the possible reintroduction of infections from foci with ongoing transmission. The degree to which individual communities are exposed to such foci will depend on – in addition to blackfly and human movement – the geographical coverage of an intervention over the wider area and, in the case of MDA, the treatment coverage in nearby communities. These factors, and the more general role of spatial transmission processes on the dynamics, control and elimination of NTDs are poorly understood, constituting a major gap in our understanding of how best to reach and sustain elimination.

An important limitation of both models’ projections is that entomological data from these foci – which between 2008 and 2010 were strongly indicative of transmission interruption (Diawara et al., 2009, Traore et al., 2012) – were not used to train the models. This is because entomological surveillance was conducted away from the communities at a few capture points, in a manner similar to OCP collection protocols during vector control. In principle, these data could be – and ultimately should be – used to inform model projections, but without detailed spatially-explicit models, this would require assumptions that flies from a limited number of collection points (4 in the River Bakoye and 3 in the River Gambia foci respectively; Diawara et al., 2009) were representative of infectivity rates among flies biting in communities. Linking spatially-distinct epidemiological and entomological data (where they exist) will be an important future challenge to refining our modelling projections.

A further technical limitation of the modelling approach is that statistical discrepancies between model predictions and data – given a particular parameter set, and assumed level of coverage and adherence – were attributed to variations in community-specific ABRs (which were estimated as part of the training procedure). While such variations in ABR are entirely plausible among geographically distinct populations, it is also possible that community-specific coverage (and adherence rates) were outside of the considered intervals (or the assumed point adherence rate) in some communities. It is further plausible that underlying population processes that were assumed to be constant between foci in this analysis are somewhat geographically variable (e.g. responses of O. volvulus to ivermectin in Ghana are distinct among circumscribed foci, possibly driven by geographically distinct parasite genetic backgrounds; Frempong et al., 2016). Hence, some model-data deviations could be alternatively ascribed to variations in coverage, adherence or more fundamental population processes rather than only to ABRs. Hence, without more precision (at the individual community level) in the programmatic information, it is not possible to conclude that more epidemiological data always give a better ABR estimate, the key determinant of local transmission conditions. Indeed, a variety of other factors such as ecological or environmental changes that might systematically have changed ABRs over the two decades of intervention also impede our ability to estimate accurately local transmission conditions.

The results of this analysis together with our previous work (Stolk et al., 2015, Coffeng et al., 2014a, Turner et al., 2014) suggest that whether the 2020/2025 onchocerciasis elimination goals will be met critically depends on when an intervention began; how intensively and effectively it was implemented (e.g. annual or biannual MDA with ivermectin; levels of coverage) and on the local transmission conditions. In areas where the prevailing ecological conditions are highly propitious to transmission – including many foci in West Africa where the initial endemic prevalence was very high (O’Hanlon et al., 2016) and vector abundance has returned to original levels following OCP vector control (Lamberton et al., 2014) – elimination by 2020 or 2025 may be unfeasible with ivermectin alone. In such places (and elsewhere, such as where suboptimal responses to ivermectin are persistent, Frempong et al., 2016), localized low-cost vector control should be considered as a complement to MDA alongside other so-called alternative treatment strategies (WHO, 2015) including anti-wolbachial therapies (Walker et al., 2015) to accelerate progress towards the 2020/2025 goals. Emerging new treatment options (Kuesel, 2016), particularly moxidectin which is similar but more efficacious than ivermectin (Awadzi et al., 2014), present promising alternatives (Turner et al., 2015b). Where elimination appears to have been achieved, it is vital that appropriate and robust PTS protocols are implemented for confirmation. Due to the poor sensitivity of the skin snip diagnostic, 3- to 5-year PTS periods of microfilarial prevalence may not be long enough to offer sufficient confidence that resurgence has not occurred. Models must now be used to explore how resurgence manifests in seroprevalence and entomological data collected using the latest and newly recommended diagnostic tools which will become the cornerstone of future onchocerciasis surveillance in Africa.

Author contributions

MW, WAS and MGB conceived the study design, undertook the modelling, the data analysis and interpretation of the results. MAD and CB helped to refine EPIONCHO and assisted with the data analysis. LD and MOT led the collection of the original epidemiological field study and assisted with the interpretation of the study design and the programmatic and epidemiological data. MW, WAS and MGB drafted the manuscript which was read, revised for critical intellectual content, modified and then approved by all authors.

Funding

We acknowledge funding of the NTD Modelling Consortium by the Bill and Melinda Gates Foundation in partnership with the Task Force for Global Health. MW and MGB also received financial support from the UNICEF/UNDP/World Bank WHO Special Programme for Research and Training in Tropical Diseases (TDR) and the Wellcome Trust (grant no. 092677/Z/10/Z). MAD would like to acknowledge the Medical Research Council-Doctoral Training Programme at Imperial College London. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The views, opinions, assumptions or any other information set out in this article are solely those of the authors.

Conflicts of interest

None.