## Figures

## Abstract

Progress towards controlling and eliminating parasitic worms, including schistosomiasis, onchocerciasis, and lymphatic filariasis, is advancing rapidly as national governments, multinational NGOs, and pharmaceutical companies launch collaborative chemotherapeutic control campaigns. Critical questions remain regarding the potential for achieving elimination of these infections, and analytical methods can help to quickly estimate progress towards—and the probability of achieving—elimination over specific timeframes. Here, we propose the *effective reproduction number*, *R*_{eff}, as a proxy of elimination potential for sexually reproducing worms that are subject to poor mating success at very low abundance (positive density dependence, or Allee effects). *R*_{eff} is the number of parasites produced by a single reproductive parasite at a given stage in the transmission cycle, over the parasite’s lifetime—it is the generalized form of the more familiar basic reproduction number, *R*_{0}, which only applies at the beginning of an epidemic—and it can be estimated in a ‘model-free’ manner by an estimator (**‘**ε’). We introduce ε, demonstrate its estimation using simulated data, and discuss how it may be used in planning and evaluation of ongoing elimination efforts for a range of parasitic diseases.

## Author summary

Critical questions remain regarding the potential for achieving elimination of helminth (worm) infections, and methods are needed to quickly estimate the probability of achieving elimination over specific timeframes. Here, we make use of methods in mathematical epidemiology to propose an estimator of elimination potential, and then apply this estimator of *R*_{eff}−the number of parasites produced by a single reproductive parasite at a given stage in the transmission cycle, over the parasite’s lifetime—as a so-called ‘elimination feasibility coefficient,’ or **‘**ε.’ We compare this metric ε to commonly used metrics of population viability for the parasitic disease schistosomiasis. The value of this estimator is that it allows us to determine whether an elimination program is succeeding and how much effort might be required to bring it to completion.

**Citation: **Arakala A, Hoover CM, Marshall JM, Sokolow SH, De Leo GA, Rohr JR, et al. (2018) Estimating the elimination feasibility in the 'end game' of control efforts for parasites subjected to regular mass drug administration: Methods and their application to schistosomiasis. PLoS Negl Trop Dis 12(11):
e0006794.
https://doi.org/10.1371/journal.pntd.0006794

**Editor: **Andy P. Dobson, Princeton University, UNITED STATES

**Received: **September 6, 2017; **Accepted: **August 27, 2018; **Published: ** November 12, 2018

**Copyright: ** © 2018 Arakala et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability: **Data for the model are available in the main text and SI of the article; these are parameter values that allow the computer models to be built and run. All code used to conduct these analyses can be found at https://github.com/cmhoove14/EliminationFeasibility

**Funding: **CMH and JVR and JRR were supported by National Institutes of Health grant #R01TW010286 and CMH and JVR also by #R01AI125842, and by the National Science Foundation Water, Sustainability and Climate grants 1360330 and 1646708, and by the University of California Multicampus Research Programs and Initiatives award MRP-17-446315. AA and MG were supported by National Institutes of Health grant #R01TW010286. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests: ** The authors have declared that no competing interests exist.

## Introduction

The World Health Organization (WHO) Neglected Tropical Disease (NTD) Roadmap [1] advocates targeting several parasitic diseases, including schistosomiasis, onchocerciasis, and lymphatic filariasis (LF), for elimination or control. Control measures generally rely on mass drug administration (MDA) and environmental interventions that interrupt the parasites’ life cycles [2]. The WHO has set goals to locally eliminate (e.g. malaria from the Greater Mekong Region by 2030 [3]) or globally eradicate (e.g. yaws by 2020 [4]) certain NTDs in the coming years, and elimination may also be within reach for onchocerciasis in Guatemala [5]; Chagas disease in parts of Central and South America [5]; lymphatic filariasis in Brazil, Haiti and The Dominican Republic [5]; trachoma in Mexico [5]; and schistosomiasis in China [6].

However, assessing the control effort required to reach elimination in these settings remains a major challenge. The standard mathematical epidemiological approach relies on fitting differential equation models to epidemiological data, and examining model parameters and associated dynamical features. Such models are simple enough to be rapidly solved and yet rich enough to encode many of the mechanisms that give rise to epidemics, e.g., person-to-person contact processes, incubation times, and infectious periods [7]. However, these models require substantial knowledge of the epidemiological system so that the essential biological processes can be properly mathematically formulated. For this reason, non-parametric approaches (i.e., approaches that do not require model parameters to be fitted to data) that analyze empirical data, such as time series of infection intensity or disease prevalence, have also been proposed to answer policy-relevant questions, such as determining when a disease is close to elimination [8].

A powerful way to use mathematical epidemiology in systems in which the full model cannot be easily specified is to use quantities from well-formulated models to generate model-free analyses of commonly collected data. We refer to analyses that are ‘model-free’ as those that do not require the construction and parameterization of a mechanistic (or simulation) model.

Here, the *effective reproduction number*, *R*_{eff}, is presented as one such candidate quantity for assessing the elimination feasibility of parasitic diseases. Because morbidity associated with parasitic infections is often related to infection intensity (i.e., the number or *burden* of adult parasites or worms harbored by human hosts) rather than the presence or absence of infection in a particular individual [9], the effective reproduction number is defined as the number of parasites produced by a single reproductive parasite at a given stage in the transmission cycle, over the parasite’s lifetime. The effective reproduction number is the generalized form of the more familiar basic reproduction number, *R*_{0}, which only applies at the beginning of an epidemic. *R*_{eff} is a function of parasite population density (it should correctly be referred to as *R*_{eff}(W) but we use the more compact form *R*_{eff} throughout this article) and thus changes as the parasite density changes over time [10]. When *R*_{eff} = 1, the parasite population is at an equilibrium, where it is expected to be perfectly replaced over the course of its lifetime. Furthermore, *R*_{eff} > 1 indicates a parasite population expected to increase in size while a population with *R*_{eff} < 1 would be expected to decrease.

The *R*_{eff} profile (Fig 1) captures the relationship between *R*_{eff} and parasite population density, providing a summary of system stability across a range of parasite burdens within the human host population. Under the influence of both positive density dependence (PDD; e.g. the mating probability of dioecious—two gendered—parasites) and negative density dependence (NDD) (see [10]), the *R*_{eff} profile intersects *R*_{eff} = 1 twice and has a humped shape. Further examples of PDDs include parasite suppression of the host’s immune response [11], and of NDDs include host immunity and parasite crowding effects acting upon parasite establishment or fecundity [12]. An *endemic equilibrium* (Fig 1, *W*_{eq}) occurs at high parasite burdens where parasite population growth is restricted by NDD (e.g., because of crowding, competition, or host immunity). The *breakpoint* (Fig 1, *W*_{bp}) occurs at a low parasite burden where strong PDD, also known as an Allee effect [10,13], facilitates parasite elimination due to reductions in parasite reproduction (e.g., due to mate limitation) [10]. The *R*_{eff} profile can be used to estimate key features of a parasite population such as the rate of rebound following treatment, the expected parasite burden at the endemic equilibrium, or the expected breakpoint parasite population size [10]. Even though the mating probability PDD is often ignored in traditional parasitic models, there are a number of studies that have included this PDD, and breakpoints, in their theoretical analyses [10,14–18].

An *R*_{eff} profile with equilibria noted at both low (unstable; *W*_{bp}) and high (stable; *W*_{eq}) parasite burdens. Between these equilibria lies a point of parasite burden (*W*_{peak}) where intense transmission is expected (transmission being proportional to the value of *R*_{eff}) due to minimized influence of restrictive negative density dependence and facilitation of transmission through amplified positive density dependence. The critical influence of positive density dependence on expected values of *R*_{eff} at low worm burdens (solid line) and its deviation from expected values in the absence of positive density dependence (dashed line), where *R*_{eff} → *R*_{0}, is also shown.

Below the breakpoint (Fig 1, *W*_{bp}), *R*_{eff} < 1 and the effects of PDD are expected to drive the parasite population to elimination. The breakpoint therefore provides a logical target for intervention campaigns that decrease parasite burden, e.g., through the use of MDA, as reducing the parasite population below its breakpoint should result in elimination even without further resources devoted to intervention. Detecting the influence of PDD on parasitic disease systems and summarizing its potential influence on the elimination feasibility of sexually reproducing parasites is therefore of great value to public health practitioners developing campaigns to eliminate parasitic diseases from human populations. We hypothesize that parasite populations influenced by PDD exhibit distinct dynamics that are detectable in epidemiological datasets measuring parasite burden or prevalence through the course of an intervention campaign, and that straightforward analysis of these datasets can shed light on the elimination feasibility of parasite populations subjected to regular interventions.

To demonstrate this, the *R*_{eff} profile of the dioecious macroparasitic disease schistosomiasis is estimated, and shown to approximate the rate at which infection intensity rebounds following MDA, a quantity referred to as the Bounce Back Rate, *BBR*. Furthermore, it is shown that data collected from longitudinal intervention campaigns (e.g., annual MDA campaigns conducted in endemic schistosomiasis communities) can generate time-series of the *BBR* that can be used to develop an estimator that quantifies the elimination feasibility of the parasite population. The estimator is evaluated by comparing it to the results of commonly used population stability analyses. This estimator is presented as a simple, model-free quantity which can be used to determine the elimination feasibility of a parasitic disease subjected to regular control efforts, and more stringent and systematic collection of data that could be used to derive it is suggested.

## Methods

A mathematical model of schistosomiasis transmission, adapted from the compartmental model developed by Anderson and May [19] and from other previous work [20,21] (which focused upon *S*. *haemotobium*) is developed to describe parasite burden, measured as mean worm burden per person, *W*, in the human host population, and three disease states of the intermediate host snail population: susceptible *S*, exposed (pre-patent, *E*) and infectious (patent, *I*). The model, expressed below as Eqs 1–4, incorporates NDD via the snail population carrying capacity, *C*, acquired host immunity, *ρ*, and parasite crowding, *γ*. A PDD is implemented via the mating probability of the obligatory dioecious *Schistosoma* spp. parasites, 𝜙. The mating probability represents the probability that an individual female worm will be successfully mated, and therefore able to produce viable eggs. It is quantified as a function of the mean worm burden per person, *W*, and its distribution amongst human hosts assuming a negative binomial distribution with aggregation parameter, *κ* [22].

### Mathematical model of schistosomiasis transmission

The basic schistosomiasis model is given by the following system of differential equations:
(1)
(2)
(3)
(4)
where the density dependent parameters, *ϕ*, *γ*, and *ρ*, are estimated as:
(5)
(6)
(7)
and *N* = *S* + *E* + *I* is the total snail population.

Parameter and state variable descriptions and values used in the model are shown in Table 1. The transmission parameters, *β* and *λ*, representing man-to-snail and snail-to-man transmission, respectively, are estimated by fitting the full model to epidemiological data collected in a community in the Senegal river basin as described elsewhere [20,21]. The aggregation parameter of the negative binomial distribution, *κ*, is also estimated directly from the data (S1 Text). The degree of parasite aggregation may differ for different stages of the parasite life cycle, but even complex simulation models generally use only a single aggregation parameter [25]; we therefore use a single parameter for the parasite negative binomial distribution to allow simple interpretation of our modelled results.

#### The effective reproduction number and bounce back rate.

Setting state variable time derivatives to zero in Eqs 1–3, and substituting into Eq 4, we obtain the rate of change of worm burden, expressed in terms of *R*_{eff} as:
(8)
where *R*_{eff} is given by (details in SI):
(9)
where , and . In Eq (8), *R*_{eff} performs the same role as the basic reproduction number *R*_{0} for the parasite dynamics close to the disease-free equilibrium (DFE) (and reduces to *R*_{0} at the DFE). We define the Bounce Back Rate, *BBR(W)* (henceforth, *BBR*), as the rate of change of parasite burden per adult parasite, expressed as:
(10)
and therefore:
(11)
The linear relationship between *R*_{eff} and *BBR*, implied by Eq 11, indicates that estimation of *BBR* can inform the estimation of *R*_{eff}, thus establishing a link between the model-based *R*_{eff} and the value of *BBR* measured through field study (see below).

#### Modeling positive density dependence.

To demonstrate the influence of PDD on *R*_{eff} and *BBR*, two versions of the model are simulated across a range of worm burdens, *W*, to generate and compare *R*_{eff} profiles: (1) the *PDD model* is run with *κ* = 0.08 [20] as estimated from the epidemiological data; and (2) the *PDD-free model* is run with *ϕ* = 1 implying all worms are successfully mated and therefore PDD has no influence on model dynamics.

### Estimating BBR from data

While the *BBR* can be estimated from the model (Eq 11), its value lies in its model-free estimation from empirical data. Consider a schistosomiasis intervention program in which an infected human population is treated with annual rounds of MDA. Previous to MDA in year (or other time frame) *i*, mean worm burden, *W*_{prei}, in the human population is assessed, while mean worm burden in the human population following MDA, *W*_{posti}, is estimated based on the coverage of treatment and the efficacy of the administered anthelminthic drug. Therefore, in what follows, *W*_{posti} is the worm burden in the population following MDA in year *i*, and *W*_{prei+1} is the worm burden in the human population before MDA in the following year *i+1*. Using the negative binomial aggregation parameter, *κ*, the prevalence before, , and following, , MDA can also be estimated via simple conversion of *W* to the probability that an individual is infected. The *BBR* in year *i* is then estimated using these longitudinal infection or prevalence data as:
(12)

This is the finite-difference estimator of the continuous time *BBR* (Eq 11), which is linearly related to *R*_{eff}, as previously established. Based on this relationship, it is next demonstrated that the underlying dynamics of disease transmission driven by PDD are detectable in empirical datasets by calculating the expression in Eq 12.

#### Model fitting and data generation.

The model expressed as Eqs 1–4 was fitted to epidemiological data detailed elsewhere [20]. Briefly, the model is used to simulate mean worm burden in treated and untreated segments of the population which is then used to estimate the maximum likelihood transmission parameter set using the *optim* function in R [26]. A profile likelihood approach is then used to estimate the 95% confidence interval of the transmission parameters that are included in subsequent simulations. In order to focus on system dynamics that demonstrate the influence of PDD (i.e., low intensity transmission), 100 low transmission intensity parameter sets within the 95% CI were used to simulate transmission over the course of 20 years with annual MDA at 80% coverage and 94% drug efficacy—corresponding to the egg reduction ratio for *S*. *haematobium* estimated in [27]—using both the *PDD model* and the *PDD-free model*. From these simulated datasets, *W*_{prei} and *W*_{posti} were calculated annually, allowing for the estimation of *BBR* over the course of a feasible, long-term MDA campaign, using Eq 12. Following 20 rounds of annual MDA, 40 years of intervention-free transmission are simulated to demonstrate the influence of PDD on the sustainability of worm burden reductions after an intervention campaign is halted.

### The elimination feasibility coefficient ‘ε’: Development and evaluation

From Eq 11 it can be seen that *BBR* should decrease and eventually become negative as *R*_{eff} approaches and passes the breakpoint over the course of long-term MDA. Because this behavior will only occur in the presence of strong PDD, longitudinal analysis of *BBR* can provide information on the presence and strength of PDD and the elimination feasibility of the disease. Specifically, adopting the simplest possible relationship, if *BBR* is related to time using a simple linear model:
(13)

The parameter ε is expected to be positive, or not significantly different from zero, for the *PDD-free model*, while it should be negative for the model accounting for mating probability and PDD. The coefficient, ε fit to longitudinal *BBR* data—termed the *elimination feasibility coefficient*—is estimated using the *PDD* and *PDD-free model*s run with each of the 100 parameter sets explained above. The ability of ε to identify disease transmission dynamics that are influenced by strong PDD—which are expected to be more susceptible to elimination from MDA-based intervention—is determined by estimating ε for each model for each year 3 < *i* < 20 in the 20 years of simulated annual MDA. In addition, for each set of estimates of ε derived from the 100 parameter sets and for each year where 3 < *i* < 20, one-sided *t*-tests are used to identify significantly negative values of ε in data derived from the *PDD* and *PDD-free* models. An estimate of ε significantly less than 0 provides evidence of the influence of PDD, as it implies *BBR* is consistently decreasing with successive treatments.

The ability of ε to predict the feasibility of elimination in parasite populations influenced by PDD was evaluated by comparing ε to the probability of extinction, *P*(*e*)—or strictly-speaking *local* extinction, for parasite control programs, though we will not use the word local in what follows—a common estimate of population viability in ecology and conservation biology [28]. Derived from stochastic models of population dynamics, *P*(*e*) explicitly incorporates demographic stochasticity, known to have a strong influence on population dynamics [28]. Simulating transmission and annual MDA with a stochastic model across a range of transmission intensities 1.3x10^{-4} ≤ λ ≤ 3.0x10^{-4}, corresponding to the upper and lower bounds of the 95% CI of *λ*, and parasite aggregation (0 ≤ κ ≤ 2) therefore provides an estimate of population viability in a variety of transmission scenarios. Stochastic versions of the *PDD* and *PDD-free model*s were developed using the Gillespie Stochastic Simulation Algorithm [29,30] with the *adaptivetau* package in R [31]. Additionally, random observation noise was added to simulated measurements of worm burden by resampling the estimated mean worm burden with the R function *rnbinom* to represent imperfect sampling and diagnosis associated with estimates of parasite burden in human host populations [15]. These same resampled worm burden data were then also used to estimate prevalence.

Simulated instances (*n* = 1,000) of longitudinal MDA intervention campaigns, as described in section 2.2.1, were generated for each of 2500 parameter sets drawn from the ranges of λ and κ above and holding all other parameters to values shown in Table 1. For each parameter set, the probability of extinction, *P*(*e*), was estimated as the proportion of model instances that result in elimination (i.e., induced extinction of the parasite population) out of all 1,000 instances. In each simulation, ε was estimated as in Eq 13 and the mean of all 1,000 estimates of ε for each parameter set was examined with respect to its correlation with *P*(*e*) estimates using a simple linear model with logit-transformed estimates of *P*(*e*) (thus transformed due to the natural [0,1] range of *P*(*e*)):
(14)
Agreement between ε and *P*(*e*) implies that the model-free elimination feasibility coefficient is a reliable estimate of how likely it is the that a parasite population will be eliminated over the course of a routine MDA campaign. In summary, the higher the positive value of the elimination potential, the further from elimination the system currently lies. As the estimator becomes negative, the system is moving towards elimination and growing negative values indicate progress in this direction. In all modelled parasite prevalence and intensity trajectories that follow, we display the output values from the model unadjusted for diagnostic sensitivity or specificity.

## Results

### Influence of positive density dependence on profiles of R_{eff} and BBR

The *R*_{eff} profile (Fig 1) indicates that, without PDD, *R*_{eff} → *R*_{0} as *W* → 0 (Fig 1, dashed line). The disease-free equilibrium (*W* = 0) is unstable for the *PDD-free model*, while it is always stable for the *PDD model* (Fig 2B, dashed lines). As anticipated, the *PDD model* produces a *BBR* profile with a humped shape, intersecting a line at *BBR* = 0 twice: at the breakpoint, *W*_{bp}, and the endemic equilibrium, *W*_{eq} (Fig 2A, solid line). This closely resembles the *R*_{eff} profile shown in Fig 1, which can indeed be recovered from the *BBR* profile through transformation according to Eq 10. In such a system, as *W* → 0, *R*_{eff} → 0 (Fig 1A, solid line) and < 0 while *W* < *W*_{bp} and *BBR* < 0 (Fig 2A & 2B).

The *BBR* profile for the *PDD* (a, solid line) and *PDD-free* (a, dashed line) *models*. Note the symmetry with the *R*_{eff} profile shown in Fig 1, including the location of key equilibria at the breakpoint (*W*_{bp}) and endemic equilibrium (*W*_{eq}) where *R*_{eff} and *BBR* = 0. Also shown is the rate of change of worm burden, (b), across log-transformed worm burden, *W*, for the *PDD* (solid line) and *PDD-free* (dashed line) models. At low values of *W* (b, inset), the models diverge as transmission in the *PDD model* (solid lines) is restricted by reduced mating probability leading to < 0 below the breakpoint, *W*_{bp}.

### Elimination feasibility coefficient from generated data

A consequence of the global stability of the endemic equilibrium of the *PDD-free model* is that the parasite density will always rebound after MDA is halted [19]. In the *PDD-free model*, remains above 0 even with repeated MDA driving *W* → 0 (Fig 3A, red line), with the potential to rebound to pre-treatment infection levels. The low value of *W* after 20 rounds of treatment in the *PDD-free model* could mistakenly be interpreted as successful elimination, yet rebound occurs after release of treatment (Fig 3A, red line). By contrast, the *PDD model* exhibits no such rebound, since repeated MDA suppresses *W* below *W*_{bp} (Fig 3A, black line), restricting transmission in the absence of exogenous parasites added to the system.

(a) Worm burden profiles derived from the *PDD model* (black line) and *PDD-free model* (red line) during 20 rounds of simulated annual MDA followed by 40 years with no intervention (the worm burden trajectory from the *PDD-free* model is shifted two months along the axis to improve clarity). In the presence of PDD, *W* remains at 0 even after releasing MDA as it has been suppressed below *W*_{bp} (a), while a lack of PDD allows *W* to rebound back to pre-MDA levels once MDA stops. (b) Mean *BBR* values of 100 model runs across 20 annual treatment rounds with error bars indicating standard deviation. (c) The elimination feasibility coefficient, ε, for the *PDD model* (black bars) and the *PDD-free model* (red bars) across each round of MDA. Error bars represent standard deviation.

With respect to bounce back rate, the *PDD model* yields a steady decline over the period of MDA and beyond caused by a decreasing mating probability as *W* → 0, and resulting in the anticipated decline in *BBR*. As *W* falls below *W*_{bp}, a negative estimate of *BBR* is also observed (Fig 3B). The *PDD-free model* shows no such decline in *BBR* with successive MDA, as the mating probability remains high irrespective of reductions in *W*, and in fact *BBR* rises with each round of MDA as the effects of NDD diminish as *W* → 0 (Fig 3B).

The slope of the *BBR* profile provides an estimate of ε (Eq 12), the elimination feasibility coefficient and, as anticipated, the *PDD model* yields a *BBR* profile with a value of ε *<* 0, while the *PDD- free model* shows ε *>* 0. When ε *<* 0, elimination is feasible using successive MDA treatments that drive *W* below *W*_{bp}, while ε *>* 0 implies greater resistance to achieving elimination.

Trajectories of prevalence rather than worm burden through the course of annual MDA intervention reveals similar dynamics that yield sustainable elimination in the *PDD model* and rebound to endemic prevalence in the *PDD-free model* (Fig 4A). Prevalence based estimation of *BBR* and ε reveals two key differences from the worm burden based estimation, however. The first is that ε is initially both positive and relatively large in magnitude (Fig 4C) as a result of relatively minor changes in the overall prevalence initially induced by MDA. This leads to ε becoming negative only after 7 years of annual MDA. However, as opposed to *BBR* estimates derived from worm burden data, every estimate of the prevalence-based *BBR* is < 0 (Fig 4B) implying that, at least in the early stages, negative estimates of the prevalence-based *BBR* may serve as a better indicator of progress towards elimination than prevalence-based estimates of ε.

(a) Prevalence profiles derived from the *PDD model* (black line) and *PDD-free model* (red line) during 20 rounds of simulated annual MDA followed by 40 years with no intervention (the prevalence trajectory from the *PDD-free* model is shifted two months along the axis to improve clarity). In the presence of PDD, prevalence remains at 0% even after releasing MDA (a), while a lack of PDD allows the prevalence to rebound back to pre-MDA levels once MDA stops. (b) Mean prevalence-based *BBR* values of 100 model runs across 20 annual treatment rounds with error bars indicating standard deviation. (c) The prevalence-based elimination feasibility coefficient, ε, for the *PDD model* (black bars) and the *PDD-free model* (red bars) across each round of MDA. Error bars represent standard deviation and the dashed line between years 7 and 8 indicates the timepoint at which ε becomes significantly (p<<0) negative.

### Probability of extinction and elimination feasibility coefficient correlations

The probability of extinction, *P*(*e*), was analyzed in a stochastic modeling framework across a range of realistic transmission intensities and degrees of aggregation of parasites among the human host population, which moderates the strength of PDD. A broad range of transmission scenarios were encapsulated in the tested parameter sets as *P*(*e*) ranged from 0 to 1. Fig 5 shows the mean ε for each parameter set and its corresponding value of *P*(*e*), each derived from 1,000 simulations of the stochastic model. The elimination feasibility coefficient was strongly correlated with logit-transformed *P*(*e*) (Eq 14, *R*^{2} = 0.94). Additionally, the proposed threshold of ε < 0 as an indicator of feasible elimination appears to be corroborated by corresponding estimates of the probability of elimination with *P*(*e*) *~ 0*.*25* when ε = 0 and *P*(*e*) → 1 as ε decreases.

Correlation between mean ε and *P*(*e*) for each parameter set derived from 1,000 simulations of the stochastic model. The model-free elimination feasibility coefficient is a strong predictor of the model-predicted probability of extinction even in the presence of demographic stochasticity and observation noise (R^{2} = 0.94; probit-transformed *P*(*e*)).

### Elimination feasibility coefficient sensitivity

At low transmission intensities, the elimination feasibility coefficient remains negative in the presence of any PDD (i.e. κ > 0; Fig 6). Examining the surface resulting from estimation of ε across a range of transmission intensities and degrees of PDD illustrates that, as transmission intensity increases, ε increases and remains greater than 0 irrespective of the degree of PDD. At these high transmission intensities, single rounds of annual MDA—as were simulated in this experiment—are insufficient to make sustainable reductions in transmission.

Sensitivity analysis showing how the elimination feasibility coefficient (ε) changes across transmission intensity, *λ*, and degree of PDD, *κ*. The surface shows a positive relationship between *λ* and ε and a negative relationship between *κ* and ε. At low transmission intensities, ε remains negative as long as κ ≠ 0 whereas high transmission intensities mostly eliminate the influence of PDD on ε, implying that annual MDA is insufficient to successfully interrupt transmission.

## Discussion

Achieving local elimination of global infectious diseases, particularly parasitic worms that cause a variety of NTDs, has become a major focus of the global health community. The availability of chemotherapeutic drugs is at an historic high due to drug donations from large pharmaceutical companies, and to a commitment from multinational NGOs, such as The Bill and Melinda Gates Foundation, to combat NTDs [32]. Yet measuring and tracking progress towards elimination remains challenging [33,34].

Here, a model-free estimator of elimination potential for parasites with sexual reproduction and positive density dependence was presented, based on longitudinal measurements of worm burden or disease prevalence in human populations treated with regular MDA. The estimator, ε, is based on observations of the bounce back rate, *BBR*, of parasitic infections in the host population following MDA, which are shown to be strongly associated with the effective reproduction number, *R*_{eff}. Negative values of ε derived from worm burden and/or negative values of prevalence-based *BBR* indicate consistent reductions in transmission that are expected to result in elimination if conditions such as MDA frequency and coverage hold.

Short-term estimates of *BBR* (measurable following each round of MDA) and long-term estimates of ε *(*measurable following three or more years of well-documented MDA) can provide a robust estimate of control status and progress towards elimination. Values of *BBR* derived from worm burden data that approach or reach 0 may indicate reductions in transmission that result in little to no reinfection in the measured time period. Furthermore, negative *BBR* following MDA may indicate a parasite population that has been driven below its breakpoint and is heading towards elimination due to reduced mating probabilities and below-replacement reproduction (i.e., *R*_{eff} < 1). Negative values of *BBR* derived from prevalence data may indicate successful reduction in transmission that can be coupled with longer term prevalence-based estimation of ε to indicate progress towards elimination. Meanwhile, high values of *BBR* (> 0.3) may indicate ineffective interruption of transmission and high rates of transmission (*R*_{eff} > 2).

Since most parasites that infect humans and have sexual reproduction exhibit PDD [19], ε estimates should tend to be negative as *W* → 0 in the presence of effective interventions. The elimination feasibility coefficient can therefore serve as a measure of progress and success in an intervention campaign, with negative estimates of *ε* implying a campaign is achieving sustained success in transmission interruption. On the other hand, positive values of ε could be indicative of a campaign failing to induce sustainable reductions in transmission due to unmeasured, untreated, or otherwise unaccounted for demographic groups [15], environmental conditions such as temperature and precipitation [35], and connectivity between transmission sites [36].

Sensitivity analysis of the model revealed the possibility that high, but reasonable, levels of transmission can overwhelm efforts to induce significant reductions in parasite burden using annual MDA alone. Positive estimates of ε and/or high *BBR* values could indicate such a situation, signaling especially challenging conditions for achieving elimination by MDA alone and possibly necessitating additional sanitation and/or education and/or snail control based interventions. Simultaneous estimation of *BBR* and ε can provide public health practitioners with a quantitative framework to assess the status of an NTD intervention as it transitions from morbidity control to transmission interruption to elimination.

Coupling estimates of ε and *BBR* may thus be valuable to infer disease dynamics over the course of an intervention campaign. Importantly, however, estimation of ε is reliant on consistent surveys of parasite burden (that is, measures of infection intensity rather than simple prevalence) in the human host population over multiple years. Indeed, the viability of any quantitative measure of the progress of control efforts hinges upon the data that are collected before, during, and after the campaign. Here, documentation of infection intensity in the host population prior to treatment, accurate assessment of the treatment coverage, and common assumptions regarding treatment efficacy are all necessary to derive values of *W*_{pre} and *W*_{pos} or *PREV*_{pre} and *PREV*_{pos} that are used to quantify *BBR* and ε. Such data are available from long-standing, well-documented MDA campaigns, yet rarely are these data made available in public databases [14]. It goes without saying that the collection and reporting of data on parasite burden (infection intensity) during MDA campaigns remains a critical priority, particularly as new methods become available to analyze and interpret such data with respect to estimating intervention efficacy. While infection intensity data may be preferred, we acknowledge the great difficulty in collecting parasite or egg counts from infected hosts, and demonstrate how similar methods can be adapted to more commonly and easily collected prevalence data [37].

In the present work, we examined performance of the elimination feasibility coefficient in predicting MDA success in the presence of observation noise, and under the influence of demographic stochasticity. Estimates of ε were compared to the probability of extinction, *P*(*e*), a well-established measure of population viability in conservation biology and population ecology which has recently been investigated in relation to elimination strategies in the parasitic disease Onchocerciasis [38]. A strong correlation between *P*(*e*) and ε was found, including a meaningful threshold at ε = 0, where there is approximately a 20% chance of elimination under the simulated MDA intervention scenario, and chances increase as ε decreases. These results suggest that the magnitude of the elimination feasibility coefficient may be a valuable indicator of the probability of extinction (i.e. the success of elimination efforts).

Future research should explicitly incorporate known sources of heterogeneity that may influence elimination feasibility, beyond those explored here. These should include seasonality, connected transmission sites, untreated demographic groups, and other sources of heterogeneity that can limit the impact of intervention campaigns and threaten their success. Methods to account for these sources of variability and uncertainty have been suggested, such as recent work to develop an accurate predictor of parasite extinction based upon a measure of the variance of disease prevalence over time [39]. Still other work [40] has sought to derive prevalence-based indicators of the elimination probability for STHs using an individually-based stochastic model, finding that end-point prevalence serves as the most reliable indicator of elimination probability. The work represented here adds to these analyses by proposing a potential mechanism underlying the observed effects, and providing additional approaches for estimating elimination probability. As a whole, this body of work demonstrates the utility of considering phenomena such as critical slowing down around population thresholds that is extensively investigated in other fields, but underappreciated for its applications to disease elimination. Critical slowing down is a phenomenon encountered in dynamical systems that are close to phase transitions, initially investigated in physics but carried over to ecology by Holling [41]. Since phase transitions occur in the vicinity of unstable equilibria, where the rates of change of system variables approach zero, we naturally expect this slowing down and we argue here that this is exploitable information in the context of disease elimination.

A key remaining question that deserves further theoretical and applied exploration is how to most effectively increase the breakpoint population size of sexually reproducing parasites subject to PDD, and thus increase the population’s susceptibility to elimination, especially when MDA alone appears insufficient to reach elimination due to very high transmission rates. If the breakpoint is very low, then immigration of a single infected individual from outside the treated area could result in the system suddenly returning to the trajectory towards the locally stable endemic equilibrium. Interventions that alter the degree of PDD over time by targeting high risk and/or heavily infected individuals may successfully increase the breakpoint, and interventions such as snail control, sanitation infrastructure, and education campaigns that reduce exposure may also increase the breakpoint by reducing the population’s overall capacity for transmission. The elimination feasibility coefficient, ε, represents a model-free summary estimate of the probability of elimination that can be used to evaluate ‘end game’ control programs aiming for parasite elimination.

## Supporting information

### S1 Text. Supporting information text file containing mathematical details of the compartmental model used in the main text and the derivation of the effective reproduction number.

https://doi.org/10.1371/journal.pntd.0006794.s001

(PDF)

## References

- 1.
World Health Organization. Department of Control of Neglected Tropical Diseases (2012) Accelerating work to overcome the global impact of neglected tropical diseases: a roadmap for implementation. World Health Organization.
- 2. Chidambaram JD, Lee DC, Porco TC, Lietman TM (2005) Mass antibiotics for trachoma and the Allee effect. The Lancet Infectious Diseases 5: 194–196. pmid:15792732
- 3.
World Health Organization. Regional Office for the Western Pacific (2015) Strategy for malaria elimination in the Greater Mekong Subregion: 2015–2030. Manila: WHO Regional Office for the Western Pacific.
- 4. Kazura JW (2015) Yaws Eradication—A Goal Finally within Reach. New England Journal of Medicine 372: 693–695. pmid:25693007
- 5. Pan American Health Organization, World Health Organization. Regional Office for the Americas (2016) Neglected infectious diseases in the Americas: Success stories and innovation to reach the neediest. Washington, D.C., United States
- 6.
Pan American Health Organization, World Health Organization. Regional Office for the Americas (2014) PAHO/WHO Schistosomiasis Regional Meeting: Defining a road map toward verification of elimination of schistosomiasis transmission in Latin America and the Caribbean by 2020. Pan American Health Organization, World Health Organization.
- 7.
Keeling MJ, Rohani P (2008) Modeling Infectious Diseases in Humans and Animals: Princeton University Press.
- 8. O'Regan SM, Lillie JW, Drake JM (2016) Leading indicators of mosquito-borne disease elimination. Theor Ecol 9: 269–286. pmid:27512522
- 9. Hotez PJ, Alvarado M, Basanez MG, Bolliger I, Bourne R, et al. (2014) The global burden of disease study 2010: interpretation and implications for the neglected tropical diseases. PLoS Negl Trop Dis 8: e2865. pmid:25058013
- 10. Gambhir M, Singh BK, Michael E (2015) The Allee effect and elimination of neglected tropical diseases: a mathematical modelling study. Adv Parasitol 87: 1–31. pmid:25765192
- 11. Gambhir M, Michael E (2008) Complex ecological dynamics and eradicability of the vector borne macroparasitic disease, lymphatic filariasis. PLoS One 3: e2874. pmid:18716676
- 12. Churcher TS, Filipe JA, Basanez MG (2006) Density dependence and the control of helminth parasites. J Anim Ecol 75: 1313–1320. pmid:17032363
- 13. Courchamp F, Clutton-Brock T, Grenfell B (1999) Inverse density dependence and the Allee effect. Trends in Ecology & Evolution 14: 405–410.
- 14. Anderson R, Turner H, Farrell S, Yang J, Truscott J (2015) What is required in terms of mass drug administration to interrupt the transmission of schistosome parasites in regions of endemic infection? Parasites & Vectors 8: 553.
- 15. Truscott JE, Gurarie D, Alsallaq R, Toor J, Yoon N, et al. (2017) A comparison of two mathematical models of the impact of mass drug administration on the transmission and control of schistosomiasis. Epidemics 18: 29–37. pmid:28279453
- 16. Ciddio M, Mari L, Gatto M, Rinaldo A, Casagrandi R (2015) The temporal patterns of disease severity and prevalence in schistosomiasis. Chaos 25: 036405. pmid:25833443
- 17. Truscott JE, Turner HC, Farrell SH, Anderson RM (2016) Soil-Transmitted Helminths: Mathematical Models of Transmission, the Impact of Mass Drug Administration and Transmission Elimination Criteria. Adv Parasitol 94: 133–198. pmid:27756454
- 18. Anderson RM, Turner HC, Farrell SH, Truscott JE (2016) Studies of the Transmission Dynamics, Mathematical Model Development and the Control of Schistosome Parasites by Mass Drug Administration in Human Communities. Adv Parasitol 94: 199–246. pmid:27756455
- 19.
Anderson RM, May RM (1992) Infectious Diseases of Humans: Dynamics and Control: OUP Oxford.
- 20. Sokolow SH, Huttinger E, Jouanard N, Hsieh MH, Lafferty KD, et al. (2015) Reduced transmission of human schistosomiasis after restoration of a native river prawn that preys on the snail intermediate host. Proceedings of the National Academy of Sciences 112: 9650–9655.
- 21. Halstead NT, Hoover CM, Arakala A, Civitello DJ, De Leo GA, et al. (forthcoming) Agrochemical pollution increases risk of human exposure to schistosome parasites. Nature Communications.
- 22. May RM (1977) Togetherness among Schistosomes: its effects on the dynamics of the infection. Mathematical Biosciences 35: 301–343.
- 23.
Anderson RM (1980) The Dynamics and Control of Direct Life Cycle Helminth Parasites. In: Barigozzi C, editor. Vito Volterra Symposium on Mathematical Models in Biology. Lecture Notes in Biomathematics: Springer, Berline, Heidelberg. pp. 278–322.
- 24. Woolhouse ME (1992) A theoretical framework for the immunoepidemiology of helminth infection. Parasite Immunol 14: 563–578. pmid:1470476
- 25. Turner HC, Truscott JE, Bettis AA, Farrell SH, Deol AK, et al. (2017) Evaluating the variation in the projected benefit of community-wide mass treatment for schistosomiasis: Implications for future economic evaluations. Parasit Vectors 10: 213. pmid:28454578
- 26.
Team RC (2018) R: A language and environment for statistical computing Vienna, Austria: R Foundation for Statistical Computing.
- 27. Zwang J, Olliaro PL (2014) Clinical efficacy and tolerability of praziquantel for intestinal and urinary schistosomiasis-a meta-analysis of comparative and non-comparative clinical trials. PLoS Negl Trop Dis 8: e3286. pmid:25412105
- 28. Melbourne BA, Hastings A (2008) Extinction risk depends strongly on factors contributing to stochasticity. Nature 454: 100–103. pmid:18596809
- 29. Gillespie DT (1976) A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics 22: 403–434.
- 30. Cao Y, Gillespie DT, Petzold LR (2007) Adaptive explicit-implicit tau-leaping method with automatic tau selection. J Chem Phys 126: 224101. pmid:17581038
- 31.
Johnson P (2016) adaptivetau: Tau-Leaping Stochastic Simulation. R package version 2.2–1. https://CRAN.R-project.org/package=adaptivetau.
- 32. World Health Organization. Department of Control of Neglected Tropical Diseases (2017) Global resolve to end neglected tropical diseases amid unprecedented progress.
- 33. Spear RC, Seto EYW, Carlton E, Liang S, Remais JV, et al. (2011) The challenge of effective surveillance in moving from low transmission to elimination of schistosomiasis in China. International journal for parasitology 41: 1243–1247. pmid:21920366
- 34. Tchuem Tchuenté L-A, Rollinson D, Stothard JR, Molyneux D (2017) Moving from control to elimination of schistosomiasis in sub-Saharan Africa: time to change and adapt strategies. Infectious Diseases of Poverty 6: 42. pmid:28219412
- 35. Remais J, Hubbard A, Zisong WU, Spear RC (2007) Weather-driven dynamics of an intermediate host: mechanistic and statistical population modelling of Oncomelania hupensis. Journal of Applied Ecology 44: 781–791.
- 36. Gurarie D, Seto EYW (2009) Connectivity sustains disease transmission in environments with low potential for endemicity: modelling schistosomiasis with hydrologic and social connectivities. Journal of the Royal Society Interface 6: 495–508.
- 37. Utzinger J, Becker SL, van Lieshout L, van Dam GJ, Knopp S (2015) New diagnostic tools in schistosomiasis. Clin Microbiol Infect 21: 529–542. pmid:25843503
- 38. Coffeng LE, Stolk WA, Hoerauf A, Habbema D, Bakker R, et al. (2014) Elimination of African onchocerciasis: modeling the impact of increasing the frequency of ivermectin mass treatment. PLoS One 9: e115886. pmid:25545677
- 39. O’Regan SM, Drake JM (2013) Theory of early warning signals of disease emergence and leading indicators of elimination. Theor Ecol 6: 333–357.
- 40. Truscott JE, Werkman M, Wright JE, Farrell SH, Sarkar R, et al. (2017) Identifying optimal threshold statistics for elimination of hookworm using a stochastic simulation model. Parasit Vectors 10: 321. pmid:28666452
- 41. Holling CS (1973) Resilience and Stability of Ecological Systems. Annual Review of Ecology and Systematics 4: 1–23.