Modelling exposure heterogeneity and density dependence in onchocerciasis using a novel individual-based transmission model, EPIONCHO-IBM: Implications for elimination and data needs

Background Density dependence in helminth establishment and heterogeneity in exposure to infection are known to drive resilience to interventions based on mass drug administration (MDA). However, the interaction between these processes is poorly understood. We developed a novel individual-based model for onchocerciasis transmission, EPIONCHO-IBM, which accounts for both processes. We fit the model to pre-intervention epidemiological data and explore parasite dynamics during MDA with ivermectin. Methodology/Principal findings Density dependence and heterogeneity in exposure to blackfly (vector) bites were estimated by fitting the model to matched pre-intervention microfilarial prevalence, microfilarial intensity and vector biting rate data from savannah areas of Cameroon and Côte d’Ivoire/Burkina Faso using Latin hypercube sampling. Transmission dynamics during 25 years of annual and biannual ivermectin MDA were investigated. Density dependence in parasite establishment within humans was estimated for different levels of (fixed) exposure heterogeneity to understand how parametric uncertainty may influence treatment dynamics. Stronger overdispersion in exposure to blackfly bites results in the estimation of stronger density-dependent parasite establishment within humans, consequently increasing resilience to MDA. For all levels of exposure heterogeneity tested, the model predicts a departure from the functional forms for density dependence assumed in the deterministic version of the model. Conclusions/Significance This is the first, stochastic model of onchocerciasis, that accounts for and estimates density-dependent parasite establishment in humans alongside exposure heterogeneity. Capturing the interaction between these processes is fundamental to our understanding of resilience to MDA interventions. Given that uncertainty in these processes results in very different treatment dynamics, collecting data on exposure heterogeneity would be essential for improving model predictions during MDA. We discuss possible ways in which such data may be collected as well as the importance of better understanding the effects of immunological responses on establishing parasites prior to and during ivermectin treatment.


Methodology/Principal findings
Density dependence and heterogeneity in exposure to blackfly (vector) bites were estimated by fitting the model to matched pre-intervention microfilarial prevalence, microfilarial intensity and vector biting rate data from savannah areas of Cameroon and Cô te d'Ivoire/Burkina Faso using Latin hypercube sampling. Transmission dynamics during 25 years of annual and biannual ivermectin MDA were investigated. Density dependence in parasite establishment within humans was estimated for different levels of (fixed) exposure heterogeneity to understand how parametric uncertainty may influence treatment dynamics. Stronger overdispersion in exposure to blackfly bites results in the estimation of stronger density-dependent parasite establishment within humans, consequently increasing resilience to MDA

Introduction
The World Health Organization (WHO)'s roadmap on neglected tropical diseases [1] has earmarked onchocerciasis for elimination by 2020 in selected African countries, and the Joint Action Forum (JAF) of the WHO African Programme for Onchocerciasis Control (APOC) proposed elimination in 80% of endemic countries by 2025 [2]. As onchocerciasis programmes based on mass drug administration (MDA) of ivermectin transition from morbidity control to parasite elimination [3], the usefulness of mathematical models will rest on our ability to identify and understand processes that may make parasite populations resilient to MDA and able to persist at low prevalence [4]. Density-dependent processes acting on various parts of parasite lifecycles are well recognised as an important aspect of helminth transmission dynamics, stabilising parasite populations and contributing to their resilience-their capacity to recover-during (and after) control interventions [5,6,7,8,9]. Positive or facilitative density dependencies (e.g. the mating probability in dioecious, separate sexes species) limit transmission at low parasite population densities and create so-called transmission breakpoints [10]. Negative or constraining densitydependent processes limit transmission at high parasite population densities and, as they are relaxed during intervention, enhance transmission at low population densities [6,7,11]. It follows that density-dependent parasite establishment may have important implications for the resilience of onchocerciasis to interventions [12,13,14]. Since parasites are typically overdispersed among hosts (i.e. the variance is substantially greater than the mean worm load), the severity of density dependence differs among individuals and, therefore, the net populationlevel effect is altered by the degree of parasite overdispersion [10,15]. A description of the biological motivation for accounting for a variety of density-dependent processes in onchocerciasis transmission models can be found in [6], [9], [16] and [17]. There are two representations of density dependence in the establishment of the parasite life stage in humans. The first was introduced by Dietz [18] and was later used in the deterministic precursor of the EPIONCHO model (see [22] for a review) and is also adapted for use (for individual human hosts) in this paper. The proportion of parasites that establish in humans is considered to be a function of the annual transmission potential (ATP = L3 × ABR, the number of L3 larvae potentially received per person per year, where L3 is the mean number of L3 larvae in the fly population and ABR is the annual biting rate (the average number of bites received per person per year)). This decreasing proportion of establishing parasites with increasing ATP has implications for parasite resilience during mass drug administration (Fig 1A). The red points indicate how the proportion of establishing L3 larvae may increase after treatment as the ATP declines. This increase in the proportion of establishing parasites constrains the ability of treatment to reduce transmission because-although the number of parasites in the fly population is reduced-a higher proportion of these parasites establish. Duerr et al. (2006) [9] presented evidence for density-dependent establishment using individual-level data. The equivalent output of Dietz's function can be reproduced using data from the 14 villages with paired nodulectomy and ATP information in the OCP database, and the model in [9] (Fig 1B). Similarly, to Fig 1A, we see that the proportion of establishing adult parasites decreases as the ATP increases.
In Onchocerca volvulus, density-dependent processes affect the proportion of microfilariae (mf, the progeny of adult worms) establishing in blackfly (Simulium spp.) vectors [16], the survival of infected blackflies [17], and the establishment of juvenile larvae as adult worms in humans [18,19,20]. A large amount of work has considered the impact of these density dependencies on the transmission dynamics of onchocerciasis using population-based deterministic models [6,21,22,23]. However, these approaches have not permitted investigation of the interactive effects of density dependence with heterogeneity in individuals' exposure to blackfly bites.
Individuals may differ in their exposure to blackfly bites due to their attractiveness to flies [24], occupation [25,26] or age and sex [21,27]. Variation in exposure to vector bites has received theoretical and empirical attention in the context of lymphatic filariasis [28,29], schistosomiasis [30] and malaria [31,32] but is less well studied for onchocerciasis (although see [33]). In lymphatic filariasis, high levels of exposure heterogeneity are associated with increased resilience to MDA, allowing parasite persistence at low prevalence [28]. Moreover, it has been shown that identical prevalence values can be produced by different combinations of vector to host ratios (indicative of population average exposure) and levels of heterogeneity; high exposure heterogeneity reduces parasite prevalence [28] and therefore, a higher vector to host ratio is required to achieve a given prevalence. Thus, it has been suggested that prevalence alone should not determine expectations on intervention success [28]; a population with high prevalence and high exposure heterogeneity will be more resilient to treatment than a population with the same prevalence, but lower heterogeneity in exposure.
Here, we consider how the interaction between overdispersion in exposure to vector bites and density-dependent establishment of O. volvulus in human hosts affects the parasite's population dynamics. We present a novel individual-based stochastic onchocerciasis transmission model, EPIONCHO-IBM (an analogue of the population-based deterministic EPIONCHO model (see [22] for a recent review)), which we parameterise using historical pre-intervention epidemiological data from savannah settings in Africa to capture the interaction between exposure heterogeneity and density-dependent parasite establishment within humans. We use the model to explore the impact of these interactive processes on O. volvulus population dynamics during MDA with ivermectin. We discuss how uncertainty in key parameters influences our projections and highlight data requirements for improving the accuracy of modelbased predictions. In addition to individual-based models (IBMs) being suited to modelling heterogeneity in exposure (and potentially also in susceptibility) to infection, they are particularly useful in context of public health questions (for example, allowing the consideration of sensitivity and specificity to diagnostic tests; accounting for screen and treat protocols, among others). A broader aim of this work is, therefore, to develop a tool with which a wide range of questions (not within the scope of this paper) may be answered that is not possible with its deterministic precursor (EPIONCHO [22]).  where NER i = k i /a i ; k i is the number of nodules palpated in individual i at age a i ; Z is the number of people across all villages entering the analysis for which nodulectomy (and paired ATP) data were available for 14 villages in the Onchocerciasis Control Programme in West Africa (OCP) database; n is the population size of the village under consideration; n a is the number of people aged a, and z a is the number of people aged a across all villages. The incidence of adult worms (parasite establishment rate, PER) in a village is then, PER ¼ m p NER, where p is the proportion of nodules accessible to palpation, and m is the average number of adult female worms in a nodule. We set m (which was multiplied by 2 to find the total number of adult worms, assuming a balanced sex ratio) and p to values proposed by [34].

Methods
There are three principal methodological components of this work: i) development of a novel individual-based onchocerciasis transmission model (EPIONCHO-IBM); ii) parameterisation of heterogeneity in exposure to blackfly bites and density-dependent establishment of worms within humans, and iii) simulations to investigate the role of these processes and their interactions on the dynamics of O. volvulus during MDA with ivermectin.

EPIONCHO-IBM
EPIONCHO-IBM is a stochastic, individual-based analogue of a previously developed population-based model, EPIONCHO (see [22] for a recent review and [23] for the latest refinements to the deterministic version). EPIONCHO-IBM follows each human in a closed population, keeping track of the number of infecting adult O. volvulus (of each sex and reproductive status) and microfilariae. The presence of both male and female worms is required for the production of microfilariae, assuming a completely polygamous mating system [10,34]. The model accounts for age-and sex-dependent exposure of humans to blackfly bites, as in the deterministic version of the model [21], whilst additionally incorporating individual-level variation in exposure. An individual-specific exposure factor is assigned at birth and is drawn from a gamma distribution, where k E and β E are the shape and rate parameters, respectively. We assume always that k E = β E , such that the mean exposure in the population is unity, permitting the ABR (the average number of bites per person per year), to be distributed among the host population.

Parasite life-history traits
We assume senescence in parasite longevity and fecundity based on existing data [35]. Mortality rates of adult worms and microfilariae are assumed to increase with age, according to a Weibull model fitted, in the case of microfilariae, to data presented by [36] (S1 Text, Text A, Formal description of EPIONCHO-IBM, Fig A). Following [37], parasite fecundity decreases with age (Fig B in S1 Text). Figures showing the dependency of parasite life history traits (mortality and fecundity) on age, and the resulting fit of the model to temporal (declining) trends of community microfilarial load (CMFL, as defined in [38]) since the inception of vector control in Burkina Faso villages of the Onchocerciasis Control Programme in West Africa (OCP), are included in S1 Text ( Fig C).

Parasite population regulation
Density-dependent processes are assumed to act on three stages of the O. volvulus lifecycle, namely: establishment of larvae within the vector; parasite-induced mortality of the vector, and establishment of adult worms within the human. Herein, reference to 'density dependence parameters' refers to the density dependence in the establishment of the parasite in humans unless otherwise specified, and the parameters determining other density-dependent processes remain fixed [23]. After accounting for individual-level variation, as well as age-and sex-specific exposure, the density-dependent establishment of adult worms in individual i is given by, where τ H is the delay between L3 larvae entering the host and establishing as adult worms to account for the duration of development into L4 stages and L5 (juvenile adults), a (i) is the age of individual i, O T (a (i) ) is the total exposure for individual i, δ H0 is the proportion of L3 larvae developing to the adult stage within the human host per bite when ATP(t) ! 0, δ H1 is the proportion of L3 larvae developing to the adult stage within the human host per bite when ATP(t) ! 1 and C H is the severity of transmission intensity-dependent parasite establishment within humans. The model accounts for a latent period in the development of the parasite in the vector by including L1, L2 and L3 stages [23,39], based on data for African (Cameroon) settings [40]. The dynamics of the parasite within the vector are modelled deterministically at a fly population level, as in the deterministic EPIONCHO model [22,23]. This population-based modelling of the vector population is also found in other individual-based models for vector borne diseases [28,41]. Treatment with ivermectin is assumed to have a large but finite microfilaricidal effect which decreases with time since treatment (following the dynamics presented in [42]). Ivermectin temporarily sterilises some female worms while making others permanently infertile [43]. A complete description and mathematical definition of the model, together with parameters values, is given in S1 Text (Text A: Formal description of EPIONCHO-IBM, Tables A to H, Figs A to D).

Parasitological data and model parameterisation
The degree of exposure heterogeneity and density dependence in parasite establishment in humans were estimated by fitting the model simultaneously to pre-intervention microfilarial prevalence and intensity (mean number of microfilariae per mg of skin) data and their corresponding annual biting rate (ABR, the number of bites per person per year) from savannah settings in [44,45] (Northern Cameroon) and [46] (Burkina Faso and Côte d'Ivoire). Although the model explicitly considers hypothetical individuals (which then allow the calculation of microfilarial prevalence and mean microfilarial infection intensity in the population), no individual-based data on exposure were used for its parameterisation as these data are not yet available at suitable scales. Instead, population-based epidemiological data, which were collected from different locations and timepoints, were used to parameterise the model regarding the processes of interest (i.e. exposure and density dependence within humans). For other processes, the model was parameterised based on epidemiological setting-independent data (e.g. density dependence within blackflies; effects of ivermectin on parasite life stages, etc.), as described in [22, 23] and S1 Text.
Parameter sets for density dependence were generated using Latin hypercube sampling (LHS, n = 100), with each set simulated for a range of k E values (0.2, 0.3, 0.4). (S1 Text, Text B: Uncertainty and sensitivity analysis using the Latin hypercube parameter sets, Fig E, illustrates the results for k E = 0.3 in the form of box-and-whiskers plots for microfilarial prevalence and intensity.) The estimation of the density dependence parameters separately for each value of k E allows the simulation of MDA for each level of exposure heterogeneity, and consequently exploration of how uncertainty in this parameter may influence the outcome of intervention programmes. The fit produced by each value of k E (with the corresponding density dependence parameters estimated for that value of k E ) can then be compared.
We calculated the sum of squared residuals (the discrepancy between the modelled microfilarial prevalence and intensity and the observed data) as a measure of the goodness-of-fit of each parameter set. Residuals were normalised between 0 and 1 to allow the prevalence and intensity to influence the fitting equally. We calculated the partial rank correlation coefficients for δ H0 , δ H1 and c H (following [47]) to quantify how each parameter influenced the pre-intervention model predictions (Fig F in S1 Text).
The microfilarial prevalence and intensity data used to fit the model [44,45,46] were complemented with additional data from [9] (OCP, microfilarial prevalence only) and other epidemiological settings with simuliid species without cibarial armatures and vector competence characteristics similar to those of Siumulium damnosum sensu stricto (s.s.)/S. sirbanum [7], namely Amazonian focus communities with transmission by S. guianense sensu lato (s.l.) [48], and Ecuadorean communities with transmission by S. exiguum s.l. [49,50] to validate the model fit. For microfilarial prevalence data, binomial 95% confidence intervals (95% CIs) were calculated using the Clopper-Pearson method [51]. Individual host microfilarial intensity data were only available for [44] and these were used for calculating 95% CIs around the arithmetic means using bootstrapping [52].

Population dynamics during mass drug administration (MDA) with ivermectin
We used the best fit parameter sets (for each value of k E ) which gave the lowest sum of squared residuals to project the dynamics of microfilarial prevalence (from two skin snips, assuming that a Holth corneoscleral punch-derived skin snip weighs, on average, 2 mg [22]) during MDA with ivermectin. Since ivermectin has been distributed in savannah areas of northern Cameroon since the late 1980's we generated predictions for annual (yearly) treatment for 25 years. For comparative purposes, we also simulated biannual (6-monthly) treatment for the same programme duration. We made predictions for populations with baseline microfilarial prevalence values � 30% (hypoendemic), 50% (mesoendemic), 60% and 70% (hyperendemic). The choice of baseline conditions was motivated by the frequency of those that were observed in the OCP prior to the commencement of control [53]. We assume that 1% of the population are non-compliant (they never receive treatment), and that the mean treatment probability in any treatment round is 0.8 (i.e. there is 80% therapeutic coverage of the total population). An overview of how compliance structure influences resilience to treatment in a generalized setting can be found in [54]. Fig 2 shows the relationship between pre-control microfilarial prevalence, arithmetic mean microfilarial intensity and ABR, using the parameters obtained by LHS that achieved the minimum squared residuals compared with the observed data from [44,45,46] for each value of k E . A value of k E = 0.3 and the associated density dependence parameters (δ H0 = 0.186, δ H1 = 0.003, c H = 0.005) gave the lowest squared residuals.

Endemic microfilarial prevalence and intensity
When using k E = 0.3, EPIONCHO-IBM captures well not only the pre-intervention observed microfilarial prevalence and intensity across the biting rates used for comparing the model fit to data, but also the other data sources used for validation (i.e. excluded from the least squares calculation) [9,48,49,50]. A value of k E = 0.2 predicted well the microfilarial intensity but tended to under predict microfilarial prevalence; conversely, a value of k E = 0.4 underpredicted intensity but was better able to capture microfilarial prevalence across all ABRs than k E = 0.2.

Threshold biting rate
Increasing the level of exposure heterogeneity (i.e. decreasing the value of k E ) decreases the threshold biting rate (TBR, the minimal ABR necessary for endemic onchocerciasis, i.e. for which the basic reproduction number of the parasite, R 0 , is equal to 1 [16]), giving TBR values ranging from approximately 97 for k E = 0.2 to 429 for k E = 0.4, assuming that the proportion of bloodmeals taken on humans is 0.63 for S. damnosum s.s./S. sirbanum [50] (Fig 3).

Population and individual-level parasite establishment within humans
We see also that the population mean density dependence in the stochastic model developed here departs from that found in the simpler deterministic framework without exposure heterogeneity (when using the estimated density dependence parameters from the stochastic model) (Fig 4A). In the stochastic setting, density-dependent establishment of the parasite in an individual depends on their exposure to fly bites (Eq 2), which is determined by a gamma-distributed individual level exposure, as well as age-and sex-dependent exposure. This creates a distribution of P H curves in the host population (Fig 4B), which can vary depending on k E (as well as on the host sex ratio and age distribution). Since (for the k E values tested), most individuals in a population have low exposure (and therefore weak density-dependent establishment relative to the annual transmission potential), there is a higher mean proportion of establishing parasites in the human population for higher annual transmission potentials than in the deterministic model.

Dynamics of microfilarial prevalence under ivermectin treatment
Treatment dynamics indicate that resilience to MDA is markedly higher for k E = 0.2 than for k E = 0.3 and 0.4. This difference between the different levels of exposure heterogeneity tended to increase as the pre-intervention microfilarial prevalence increased but decreased with an increase in treatment frequency (Fig 5). Although this difference in resilience between the different values of k E is in part due to the resulting exposure heterogeneity, it is also associated with different estimated strengths of density dependence (Fig 4). More heterogenous exposure leads to the estimation of stronger density dependence, principally due to a higher proportion of parasites establishing as ATP decreases. Since increasing exposure heterogeneity (reducing k E ) reduces the microfilarial prevalence, a higher proportion of parasites establishing (than when exposure is less overdispersed) is required to reach a given observed prevalence and intensity. Simulations assuming either a fixed level of exposure heterogeneity, and different strengths of density dependence, or a fixed level of density dependence and different levels of exposure heterogeneity are shown in S1 Text, Text B: Uncertainty and sensitivity analysis using the Latin hypercube parameter sets, Figs G and H. Pre-intervention Onchocerca volvulus microfilarial prevalence and intensity vs. the annual biting rate of simuliid vectors. The predicted microfilarial prevalence (percent) and microfilarial intensity (mean no. of microfilariae, mf, per mg of skin) (from 2 skin snips) for the annual vector biting rates reported in the combined epidemiological dataset (i.e. fitting and validation data, solid colour circles), using the estimated parameters, are represented by solid lines. The best-fit parameter values of the density-dependent parasite establishment within humans were δ H0 = 0.186, δ H1 = 0.003 and c H = 0.005 for exposure heterogeneity parameter k E = 0.3 (thick red line); δ H0 = 0.385, δ H1 = 0.003, c H = 0.008 for k E = 0.2 (thin light blue line), and δ H0 = 0.118, δ H1 = 0.002, c H = 0.004 for k E = 0.4 (thin dark blue line). The EPIONCHO-IBM predictions are based on a host population size of 500 and 100 runs for 80 years to reach endemic equilibrium (human demography was simulated to equilibrium before simulating the epidemiology to equilibrium). The error bars are binomial (Clopper-Pearson) 95% confidence intervals for prevalence and bootstrapped 95% CIs for intensity (for which the raw individual microfilarial intensity data were available [44]). The main vectors in each setting are Simulium damnosum s.s./S. sirbanum (African savannah in Benin, Burkina Faso, Cameroon, Côte d'Ivoire, Ghana, Guinea, Mali, Togo); S. guianense s.l. (Amazonian focus in southern Venezuela) and S. exiguum s.l. (Cayapas focus in Ecuador), all species are without armed cibaria and have vector competence features similar to those of S. damnosum s.s. [7]. Fitting data are from [44,45,46]; validation data are from [9,48,49,50]. The prevalence and intensity data points for an annual biting rate of 81,000 were excluded from the fitting of the model as according to [45]

Discussion
We have developed a novel stochastic individual-based onchocerciasis transmission model, EPIONCHO-IBM, based on the deterministic population-based analogue, EPIONCHO. EPIONCHO-IBM captures explicitly the interactive effects of heterogeneity among humans in exposure to blackfly vector bites and density-dependent processes that operate fundamentally at the individual host level and that are key determinants of the resilience of parasite populations to intervention. Modelling this interaction is, therefore, essential for predicting the likely impact of interventions on parasite population dynamics and ultimately the feasibility of elimination. This is the first study to explore this effect in the context of onchocerciasis using Exposure heterogeneity and density dependence in onchocerciasis transmission epidemiological data to parameterise and validate these interactive population processes. In doing so, we demonstrate that data collection on exposure heterogeneity may make an important contribution to reducing uncertainty when modelling parasite dynamics during MDA.
The range of exposure heterogeneity levels explored (k E = 0.2 to 0.4) fits to the variation in pre-control microfilarial prevalence and intensity with ABR, encompassing most of the data points collated (including fitting and validation datasets), with k E = 0.3 providing the best overall fit. A value of k E = 0.2 tended to underestimate prevalence but fitted better low prevalence values, whilst k E = 0.4 tended to overestimate prevalence but fitted better high prevalence. The results suggest that exposure heterogeneity may decrease as the annual biting rate increases (i.e. k E increases with increasing ABR).
This decrease in heterogeneity with increasing ABR may represent that, in a small population of blackflies, the physical limit on how many times a fly can bite constrains the number of people that may be bitten. This would result in a high heterogeneity in exposure to bites. Conversely, a large population of flies is capable of biting many people, leading to lower levels of heterogeneity in exposure to bites. If such a relationship between exposure heterogeneity and ABR were supported by other data sources and incorporated into the model, it may have substantive implications for the predicted population dynamics during MDA. In particular, the degree of heterogeneity in exposure would be less in highly endemic settings with high prevalence and high ABRs. Consequently, the resilience to intervention in these settings would be reduced and the impact of MDA enhanced. The projected parasite population dynamics during treatment with ivermectin would likely be somewhat more closely aligned with the corresponding deterministic predictions of EPIONCHO, which incorporates a functional relationship between the degree of adult worm overdispersion and ABR [23]. The assumption of a fixed level of exposure heterogeneity (either across the wide range of annual biting rates collated, as in Fig 2, or for single modelled endemic communities, as in Fig 5) is unavoidable due to data availability. More likely, there will be a distribution of k E values depending on the specific ecological and sociological context of a geographical area. Even if such a distribution allowed k E to vary between 0.2 and 0.4, we might expect substantial spatial variation in population responses to treatment for a given baseline microfilarial prevalence.
The threshold biting rate decreased with increasing levels of exposure heterogeneity which is consistent with similar work on lymphatic filariasis [28]. Furthermore, the threshold biting rate for the best fit parameters was similar to previous deterministic versions of the model assuming the same proportion of blood meals taken on humans (380 in [6], 300 in [19]), but lower than a different model based on the same data (730 [20]).
A key distinction between EPIONCHO-IBM and its deterministic counterpart is that in the former, overdispersion is generated mechanistically via heterogeneity in exposure and density dependence operating at the level of the individual host. In EPIONCHO, the effects of parasite overdispersion on the severity of density dependence are modelled empirically [15,21,22,23]. More generally, worm overdispersion in the host population is an input in deterministic

Fig 5. Microfilarial prevalence dynamics during 25 years of annual (a-d) and biannual (e-h) mass drug administration (MDA) with ivermectin for various levels of endemicity, heterogeneity in exposure to vector bites (measured by k E ) and the associated density-dependent establishment parameters.
The baseline microfilarial prevalence (30%, 50%, 60%, 70%, indicative of hypo-, meso-, hyper-, and high hyperendemicity, respectively) is modelled by increasing the annual biting rate. At each treatment round, 80% of the population receive treatment, excluding children aged <5 years who are ineligible to receive ivermectin. We assume 1% of the total population are systematic non-adherers (i.e. they never take treatment during the programme duration). The EPIONCHO-IBM predictions are based on a host population size of 500, 1000 runs and with 80 years to reach endemic equilibrium before initiating MDA. Parameters for the thick red lines (k E = 0.3), thin light blue lines (k E = 0.2) and thin dark blue lines (k E = 0.4) are as in Fig 2. https://doi.org/10.1371/journal.pntd.0007557.g005 Exposure heterogeneity and density dependence in onchocerciasis transmission population-based models, and an output of individual-based stochastic models (calculated by fitting a negative binomial distribution to predicted worm burdens), depending on exposure heterogeneity and potentially density dependent processes. This is a crucial difference and implies individual-based models are not simply approximations of their deterministic counterparts. How density-dependent processes influence the relationship between exposure heterogeneity and parasite overdispersion (of adult worms and microfilariae) in EPIONCHO-IBM is an important question which remains to be addressed and will be presented elsewhere.
The phenomenological nature of density-dependent establishment of the parasite in the human host is an important limitation of this work. The use of population-level data to simultaneously estimate both density dependence and exposure heterogeneity allows the parameters involved to counteract, reducing identifiability and diminishing mechanistic interpretation. Since direct estimation of density-dependent parasite establishment within humans is not feasible experimentally (with the only available observational study on parasite establishment rates being that of [9]), data collection on heterogeneity in exposure to fly bites will be an important step in better resolving density dependence. This has been discussed previously [26] but available data remain limited. Since blackflies have very specific environmental requirements, with breeding sites varying in distance to human settlements, the estimation of individual-level variation in exposure to bites poses a substantial challenge.
Heterogeneity in exposure to mosquito bites has been estimated by comparing DNA extracted from (indoor resting) mosquito bloodmeals with DNA from individuals living in study households [32]. Since Anopheles mosquitoes breed in standing water, which may be littered throughout villages, this is a viable method for finding mosquitoes which have fed on study participants. In the context of onchocerciasis, this process may be impractical because blackflies do not bite and rest indoors, and outdoor resting sites are very difficult to find, so studies investigating bloodmeal origin rely on host-dependent and host-independent sampling methods [55]. The former (human landing catches) would require a large number of participants to find enough flies which that could be matched to the individuals on which they previously fed, and the latter (using oviposition traps in nearby breeding sites) may not sample sufficient flies which would have fed on the residents of the nearby communities [55]. Irvine et al. [29] have suggested there may be spatial variation in exposure heterogeneity to mosquito bites, adding an additional layer of complexity and uncertainty in obtaining generalizable estimates. Perhaps more suitable for estimating exposure to blackfly bites is the method which has recently been applied in an empirical study of Leishmania infantum transmission among dogs in Brazil [56]. Levels of (IgG) antibody responses to sandfly saliva were found to vary with the intensity of transmission, declining after periods of low transmission, as well as to vary greatly between dogs, correlating with the intensity of transmission experienced by individual dogs. Studies such as these require the development of anti-saliva antibody assays, which are currently lacking in simuliids (partly because their colonisation in the lab is very difficult). If these assays could be developed in the context of onchocerciasis, based initially on titres of antibody responses to crude blackfly salivary antigens, they could be used in population samples across all age and sex groups and may provide useful data for estimating heterogeneity in individuallevel exposure. Such a study may also provide data to increase the accuracy of assumptions regarding age-and sex-dependent exposure. Currently this is estimated by fitting models to data on age-and sex-specific profiles of microfilarial prevalence and intensity [21]. However, these modelled age and sex exposure profiles are determined not by exposure to vectors and parasite infective stages alone, but also by processes involving parasite development and parasite fecundity within the host, as they rely on (downstream) microfilarial data. Reducing the uncertainty in age-and sex-dependent exposure patterns is particularly important for optimal selection of age (and potentially sex) groups for (serological) assessment of exposure to parasite antigens in foci thought to be nearing elimination [57]. Ideally, assays for anti-blackfly saliva could be combined with assays for exposure to O. volvulus to investigate both exposure to vector bites and parasite antigens (e.g. Ov16 [58]).
In addition to exposure, susceptibility to infection may also vary between individual hosts. Work on schistosomiasis has shown that there is marked variability in susceptibility to infection [59,60], and that accounting for both exposure and susceptibility can better capture Schistosoma transmission to snails than accounting for exposure alone [61]. It is possible, that what we term exposure heterogeneity, does also account for individual variation in susceptibility, given the way it is estimated. However, as there is a lack of data on individual heterogeneity in exposure and susceptibility to infection, making this distinction explicit is unlikely to lead to better parameterisation or improved accuracy of the model. In addition, since density-dependent establishment of adult worms might represent immunological processes (which interact with exposure at the individual level in the model), it is not necessarily clear how variation in susceptibility (which is likely determined by genetic and immunological processes) should be accounted for. Therefore, although we do not account for variation in susceptibility explicitly, the proportion of incoming parasites which establish in an individual depends on their exposure, implying a relationship between exposure and susceptibility.
The density-dependent within-human parasite establishment function used in EPIONCHO and EPIONCHO-IBM, derived from the work in [18], is a phenomenological representation of a possible immune-mediated response to infective (L3) larvae driven by the intensity of exposure to these parasite stages. In areas of low ATP (or in individuals with little exposure), low levels of exposure to L3 larvae may be responsible for a poor development of protective immune responses, leading to high parasite establishment rates. Conversely, in areas of high ATP (or in individuals highly exposed), stronger immune responses against incoming worms may decrease parasite establishment rates. The operation of anti-L3 responses in putative immunes living in areas of hyperendemic onchocerciasis has been shown by [62,63]. In the Teladorsagia circumcincta-Scottish Blackface sheep parasite-host system, priming of the immune system early in the season by exposure to (ingestion of) infective L3 larvae in pasture, reduces parasite establishment and growth and, therefore, faecal egg counts in a density-dependent fashion, not unlike our own parasite establishment function (see Figs 2 and 3 of [64]).
Density-dependent processes may also (or instead) act on parasite fecundity [65]. Here we assume for simplicity, that one worm of each sex is required for reproduction, i.e. that one male worm is sufficient for all adult female worms to be fertilised and produce mf [10,34], and that increasing female worm density does not reduce the per capita fecundity rate [19]. However, if in addition to density-dependent parasite establishment, the fecundity rate was to decrease as the density of adult parasites increases, we might expect additional resilience during treatment, since the fecundity rate per worm would increase as MDA reduces the number of parasites in the population. Of interest, previous work indicates that, when investigated separately, density-dependent fecundity contributes less to the rate of bounce back following treatment than density-dependent adult worm establishment [8].
The EPIONCHO-IBM projections indicate that there is a disproportionate increase in resilience with increasing pre-intervention prevalence for k E = 0.2. That is, as pre-intervention prevalence increases, particularly from low hyperendemic (microfilarial prevalence � 60%) to high hyperendemic (microfilarial prevalence � 70%) settings, the O. volvulus population becomes disproportionally more resilient to MDA. This can be partly understood in the context of the strongly nonlinear prevalence vs. ABR relationship, in which prevalence begins to saturate at ABR values exceeding 7,000 (daily biting rates > 20/person,~70% mf prevalence). This phenomenon would also be found for lower levels of exposure heterogeneity, as the prevalence vs. ABR relationship begins to saturate, albeit at higher ABRs. Another contributory factor to resilience to MDA-based interventions in our model is that our parasite establishment function is not affected by the potential interaction between microfilaricidal treatment and the immune response. Therefore, as treatment progresses and the ATP declines, the parasite establishment rate inevitably increases. However, ivermectin-facilitated immunity in onchocerciasis has been reported [66,67]. Human immunological studies have demonstrated that filarial parasites induce a state of hypo-responsiveness in the host that is associated with the presence of circulating mf (patent infection), and O. volvulus is no exception [68]. The reversal of this mf-associated immune-suppression, following clearance of skin mf due to ivermectin treatment, may contribute to controlling O. volvulus infection [66]. In the O. ochengicattle system, animals treated with ivermectin and exposed to blackfly bites under natural transmission conditions did not develop patent infections whilst treated but recovered their susceptibility and acquired infection (at rates higher than untreated counterparts) once treatment stopped [69]. Similar observations have been made in a large pharmaco-epidemiological study of Dirofilaria immitis in dogs [70]. A better understanding of how ivermectin-facilitated immunity could impact the establishment of incoming worms in human onchocerciasis would greatly improve the modelling of parasite establishment rates during the implementation and after cessation of MDA programmes.
The microsimulation model for onchocerciasis ONCHOSIM [22,33,41,71], which also uses a gamma distribution to model individual-level variation in exposure, uses parameter k E as 1 or 3.5, giving substantially less exposure heterogeneity than in EPIONCHO-IBM (k E = 0.3). It has been shown that EPIONCHO (the deterministic analogue of EPIONCHO-IBM) predicts more resilience to intervention than ONCHOSIM under a range of treatment scenarios (using k E = 3.5 for ONCHOSIM) [13,14]. Although our treatment dynamics were simulated under high (and likely unrealistic) levels of therapeutic coverage and adherence, the programme duration (25 years) was motivated by that in the Vina valley of northern Cameroon (from which the data used to estimate the parameters investigated in this paper originated [44,45]). In Cameroon, there was no interruption of transmission after 15 [72,73], 17 [74], 18 [75] or 25 years [76] of annual ivermectin MDA in some communities. Although this may be due to a number of factors including lower than reported coverage, decreased ivermectin efficacy [77], movement of infected individuals between foci and spatial variation in exposure heterogeneity, it may indicate that the more pessimistic predictions of EPIONCHO-IBM for hyperendemic populations, are not necessarily out of touch with observed trends. It is noteworthy that EPIONCHO-IBM can mimic the behaviour of ONCHOSIM (i.e. less pessimistic treatment dynamics) by removing density-dependent parasite establishment within humans (absent in ONCHOSIM) and increasing k E . This process-and its interaction with heterogeneity in exposure-is a critical determinant of the different population dynamics predicted by the two models (to be formally discussed elsewhere).
EPIONCHO-IBM has a wider scope of application than its population-based predecessor, albeit at the cost of reduced tractability (a general drawback of many individual-based stochastic models). A particular advantage is the capacity to model individual-based interventions. For example, individuals co-infected with O. volvulus and Loa Loa are at risk of severe adverse events (e.g. encephalitis) following treatment with ivermectin if harbouring heavy loiasis microfilaraemia [78]. Consequently "test and treat" (testing for O. volvulus) or "test and not treat" (testing for L. loa) strategies have been proposed and trialled for control in areas coendemic with onchocerciasis and loiasis [79]. EPIONCHO-IBM, as an individual-based model, can simulate such strategies, incorporating performance features (i.e. sensitivity, specificity) of a variety of diagnostics and detailed information on observed screening, coverage and adherence patterns among demographic groups. It follows that co-infection dynamics may also be modelled, although the within-host interaction between O. volvulus and L. loa and helminth species more generally is poorly understood (but see [80]).

Conclusions
In conclusion, we have developed a novel individual-based stochastic onchocerciasis transmission model, EPIONCHO-IBM, based on the well-established deterministic analogue, EPIONCHO. We have used EPIONCHO-IBM to better understand how density-dependent processes-in particular the density-dependent establishment of newly acquired O. volvulus parasites-and heterogeneity in individual human exposure shape both the relationships between microfilarial prevalence, microfilarial intensity and ABR, and the resilience of onchocerciasis to MDA with ivermectin. In future, EPIONCHO-IBM will be used to model the control and elimination of onchocerciasis using current and alternative interventions, including the use of anti-Wolbachia therapies [81], moxidectin [82], ground-based vector control [83] and new macrofilaricidal therapies [84]. Our work also highlights the importance and uncertainty in the key and interactive population processes of density dependence and heterogeneity in exposure to blackfly vectors. Novel approaches for generating data on exposure heterogeneity and anti-L3 immunity during and after the cessation of ivermectin MDA programmes would be particularly valuable in helping to resolve outstanding uncertainty on their relative importance to the population dynamics of O. volvulus. Combined with more statistically advanced approaches for estimating the parameters of individual-based models [85], and fitting such models to longitudinal epidemiological trends, such data would greatly enhance the reliability and accuracy of onchocerciasis modelling projections.

Acknowledgments
We are grateful to Alfonso Renz for providing data, and Gaby Harrow for data entry. We thank Simon O'Hanlon, Wilma Stolk and Sake de Vlas for helpful discussions during model development.

Data/code availability
The model code is available at https://github.com/jonathanhamley. All data used for the analyses are contained in the figures.