 Research
 Open Access
 Published:
Refined stratifiedwormburden models that incorporate specific biological features of human and snail hosts provide better estimates of Schistosoma diagnosis, transmission, and control
Parasites & Vectorsvolume 9, Article number: 428 (2016)
Background
Schistosoma parasites sustain a complex transmission process that cycles between a definitive human host, two freeswimming larval stages, and an intermediate snail host. Multiple factors modify their transmission and affect their control, including heterogeneity in host populations and environment, the aggregated distribution of human worm burdens, and features of parasite reproduction and host snail biology. Because these factors serve to enhance local transmission, their inclusion is important in attempting accurate quantitative prediction of the outcomes of schistosomiasis control programs. However, their inclusion raises many mathematical and computational challenges. To address these, we have recently developed a tractable stratified worm burden (SWB) model that occupies an intermediate place between simpler deterministic mean worm burden models and the very computationallyintensive, autonomous agent models.
Methods
To refine the accuracy of model predictions, we modified an earlier version of the SWB by incorporating factors representing essential inhost biology (parasite mating, aggregation, densitydependent fecundity, and random eggrelease) into demographically structured host communities. We also revised the snail component of the transmission model to reflect a saturable form of humantosnail transmission. The new model allowed us to realistically simulate overdispersed eggtest results observed in individuallevel field data. We further developed a Bayesiantype calibration methodology that accounted for model and data uncertainties.
Results
The new model methodology was applied to multiyear, individuallevel field data on S. haematobium infections in coastal Kenya. We successfully derived agespecific estimates of worm burden distributions and worm fecundity and crowding functions for children and adults. Estimates from the new SWB model were compared with those from the older, simpler SWB with some substantial differences noted. We validated our new SWB estimates in prediction of drug treatmentbased control outcomes for a typical Kenyan community.
Conclusions
The new version of the SWB model provides a better tool to predict the outcomes of ongoing schistosomiasis control programs. It reflects parasite features that augment and perpetuate transmission, while it also readily incorporates differences in diagnostic testing and human subpopulation differences in treatment coverage. Once extended to other Schistosoma species and transmission environments, it will provide a useful and efficient tool for planning control and elimination strategies.
Background
Parasitic Schistosoma species pose a significant health burden in many developing countries [1]. Broadbased regional schistosomiasis control and local elimination of parasite transmission have been prioritized in the 2012 London Declaration for Neglected Tropical Diseases (http://unitingtocombatntds.org/resource/londondeclaration) and the recent World Health Organization 2020 Roadmap on Neglected Tropical Diseases (NTDs) [2].
The parasite has a complex ecology in which it cycles between human and snail hosts through intermediate larval stages, in a manner that strongly embeds transmission within atrisk subtropical and tropical ecosystems [3]. A key feature of Schistosoma infections (similar to other metazoan macroparasites) is the highly uneven (heterogeneous) distribution of infection burden among its vertebrate host populations, as evident in both experimental and field data [4–6]. Other heterogeneous factors that play important roles in perpetuation of schistosomiasis include variations in local human demographics and exposure frequencies, patchy transmission habitats, and the seasonality of rainfall and temperature factors that affect snail abundance [7].
Modeling of such systems is a challenging task, but to provide more accurate estimates for current control programs, newer models should account for the influential features of inhost biology, transmission environment, diagnostic uncertainties, and the potential efficiency of different control interventions. Conventional approaches based on mean worm burden (MWB) formulations [4, 8–11] have had several shortcomings in this respect. They have used adhoc assumptions about worm load distributions and appear to have oversimplified some components of the transmission system. As a result, infection rates in such systems and the modeled impacts of treatment tend to be overestimated [12]. Individualbased modeling approaches (e.g. [13]) could potentially address some of these issues, but individualbased models have significant limitations in terms of accessibility and programming requirements, particularly for large populations.
The stratified worm burden (SWB) approach occupies an intermediate place between MWB and individualbased models, and offers many advantages [12, 14]. Among other features, it provides a natural way to account for worm dispersion in demographicallystructured host populations. Our earlier work with SWB [14, 15] was limited in terms of withinhost parasite biology, as we had assumed perfect mating and uniform eggrelease by all worm malefemale pairs, independent of accumulated burden. Here, we refine the earlier SWB model to account for influential components of inhost biology including the aggregated distribution of worm burden, worm mating probabilities, densitydependent worm fecundity, and random features of eggrelease [12, 16]. These newly incorporated biological parameters allow us more realistic simulation of diagnostic eggtests and of environmental egg release and its effect on the force of infection to intermediate host snails. We have developed a new Bayesian calibration procedure that recognizes many uncertainties of modeling infection data based on standard eggcount tests [17–22]. The goal of this calibration approach was to identify the most ‘likely’ model parameter choices consistent with the collected field data. The posterior parameter ensembles for a community could then be used for dynamic simulations of control interventions, incorporating the uncertainties about diagnosis and transmission reflected in the input data. The model/data uncertainties thus yield more robust statistical estimates, including credible ranges for the projected treatment outcomes.
In the current paper, we explain how we have applied the new model and calibration methodology to the field data collected in survey and control studies of S. haematobium in coastal Kenya [23–26]. These data provide a fairly complete demographic coverage of several communities, wellsuited for our analysis. Because we had also used this data set in our previous modelling using a simpler SWB model (i.e. without inhost biology [14, 15]), we can compare the two models to demonstrate how inclusion of these biological factors can result in more accurate projections of posttreatment infection outcomes.
Methods
Description of stratified worm burden system with inhost biology
In the SWB approach [12, 14], the dynamic variables included in the modeling framework are host population strata {h _{ k }(t)}, (Table 1) which are defined by their actual worm burden values. In our formulation, strata {h _{ k }} and transition rates between them are determined by a fixed worm count increment Δw  (Table 2). We think of Δw as mating threshold, so that hosts carrying less than Δw worms (stratum h _{0}) are considered noninfective.
Because low worm burden is very difficult to measure in living humans, there are no accurate estimates of the relevant minimal Δw. Theoretical arguments suggest relatively low values, Δw = 5 [9, 14]. Experimental data from primate infections [27] predicts Δw ≈ 40. Currently, we use an intermediate value of Δw = 10 in our numeric implementation of the SWB system.
A single SWB system describes a homogeneous population, determined by human FOI λ (rate of worm accumulation), worm mortality γ, and host turnover (demographic) rate, μ. Dynamic variables {h _{ k }(t)} obey coupled differential equations with matrix A(λ, μ, γ) and sources {S _{ k }(t)}. The latter account for demographic changes to SWB populations (birth, death, maturation, migration; for details see Additional file 2 and [12, 14]). In some applications (e.g. model calibration), we ask for equilibrium solutions for SWB systems to align with endemic infection levels. For a single SWB system, equilibrium distribution ĥ* = {h _{ k } ^{*} (λ/γ, μ/γ)} depends on two dimensionless parameters (λ/γ, μ/γ) (Additional file 2). For μ = 0 (no population turnover), the sequence {h _{ k } ^{*} } becomes a Poisson distribution with mean \( \overline{w}=\lambda /\gamma \), which is equivalent to the equilibrium worm burden (MWB) of a MacDonaldtype system [8]. For typical demographic turnover (small μ/γ ≪ 1), distribution {h _{ k } ^{*} } values are close to Poisson or negative binomial, with high aggregation (see Fig. 1 and [12]).
Random eggrelease by hosts and SWB communities: simulation of eggtest data
In the updated SWB, two factors determine eggaccumulation by human hosts: the number of fertilized females (mated worm pair count), ϕ, and the worm fecundity factor, ρ. Both depend on worm burden w, and can be estimated as functions of w, or for SWB models, the stratum number, k. The mated worm count depends on worm mating patterns and worm accumulation in hosts [12, 28]. Some commonly used assumptions about mating, i.e. random worm acquisition and monogamous mating, yield a binomiallydistributed sex ratio in each wstratum (w adult worms), hence the mated count,
for the h _{ k }stratum. Worm fecundity ρ(w) (or ρ _{k}) is expected to drop with increased burden due to densitydependent crowding effects. Following Anderson & Medley’s approach [29, 30], we have used the exponential decay function
with maximal value ρ _{0}, and threshold burden, w _{0}, or \( {k}_0=\frac{w_0}{\varDelta w} \) (for h _{ k }) to simulate this phenomenon. Equations (1) and (2) predict mean eggrelease by each host carrying w worms
Observed eggcounts, e.g. diagnostic test data, are typically overdispersed (Fig. 2) and can be highly variable from day to day. A negative binomial (NB) distribution of daily egg output per worm has been proposed [16], and we adopt this premise at the level of individual egg output per host (Fig. 3). Specifically each fertilized female in h _{ k }stratum is assumed to provide a random (NB) daily eggcount in the stool with mean ρ _{ k } and aggregation r. Then eggrelease by a given h _{ k }host (carrying ϕ _{ k } fertilized worms) is also negative binomial, with mean E _{ k } = ρ _{ k } ϕ _{ k } and aggregation r _{ k } = rϕ _{ k }. For the overall SWB community with strata {h _{ k }}, eggtest results are random samples drawn from the mixed NB distribution (Table 3).
In our SWB model development, simulated (random) community eggtests were used extensively for model calibration and later for prediction and analysis of control intervention outcomes.
Unlike the results of individual diagnostic tests, which influence estimates of human infection prevalence, environmental eggrelease by SWB community was considered a deterministic process that accumulates the random contributions of multiple hosts. Hence, in determining humantosnail FOI, random eggrelease (Eq. 4) is replaced by its mean value
and function E represents average host infectivity of SWB community {h _{ k }}.
Details of the local snail population model
The snail dynamics modeling combines snail population biology and infection processes (Table 4). We assume the snail population obeys a logistic growth with maximal reproduction rate, β _{ 0 }, and carrying capacity, K, following models that are commonly used in population biology to account for growth in resourcelimited environments. For modeling Schistosoma transmission within snail populations, three compartments (susceptibleexposedinfected (SEI)) are used, labeled as: x, susceptible; y, prepatent infected; and z, patent/shedding infected, with total population N = x + y + z. It is assumed that shedding snails do not reproduce [31].
Here a successful miracidial invasion (Table 4) is followed by the prepatent stage lasting 1/r days, after which a fraction c of prepatent snails convert to the patent/shedding stage (z), while the remaining fraction (1 − c) may resolve infection to return to the susceptible (x) stage. Snail FOI Λ is determined by human infectivity (E), the human population size, and humansnail contact rate, ω.
Humantosnail force of infection in coupled SWB systems
Humansnail transmission is mediated by two larval stages, cercaria (C) and miracidium (M), which determine snailtohuman FOI (λ) and humanto snail FOI (Λ), respectively. It is convenient to measure all populations (human H, snail N, M, and C) by their densities per unit habitat.
Human and snail FOI are determined by larval equilibria (Eq. (6) of Table 4), but their functional forms require more detailed analysis. For human FOI (rate of worm accumulation) we expect a linear dependence on C, proportional to patent snail prevalence. Hence
for exposure rate ω, patent snail density N z, and (snailtohuman) transmission coefficient α that accounts for intermediate cercaria stage, and the probability of worm establishment in a human host. For snail FOI, many conventional modeling approaches have adopted a similar linear relation Λ = bωH E with (humantosnail) transmission coefficient b, but this relation is questionable. Snail infection is accounted by prevalence variables (y, z), and the functional relation between miracidial density (M) and FOI Λ requires a more careful analysis. We propose, instead, that the link is a nonlinear function that takes into account two processes: (i) multiple possible Minvasions [32] and (ii) sporocyst establishment in susceptible snails [32]. The invasion process likely depends on the average number of miracidia per snail, M/N, and a snail innate resistance level, p = the probability of ejecting an invading miracidium. Multiple biological and environmental factors could contribute to a successful miracidial invasion, and p serves as a crude proxy for their cumulative (mean) effect in the susceptible snail population. Having fixed M/N, we estimate the fraction of successfully invaded snails by 1 − ρ ^{M/N}. The resulting snail FOI is the product of Λ _{0}  the rate of sporocyst development in snails (minimum 1–2 weeks [31, 32]), and the invaded snail fraction:
It is convenient to replace snail resistance p = e ^{− α} by susceptibility parameter, α = 1n(1/p), and using miracidia equilibrium M* = β _{ M } ωH E (Eq. 6) to write Λ as function of human infectivity E.
The key inputs in Λ are sporocyst establishment rate Λ_{0} and the transmission coefficient b = αβ _{ M } (product of “miracidium coefficient” β _{ M } times “snail susceptibility” α). At low transmission intensity (H E ≪ 1) function (Eq. 9) is approximately linear \( \left(\varLambda \approx {\varLambda}_0b\omega \frac{H\;E}{N}\right) \)  the conventional form of snail FOI. But unlike a “linear” Λ function, Eq. (9) would saturate at a maximum level = Λ_{0} for large H E. This has important implications for transmission dynamics and model calibration. In general, one could expect higher values of estimated coefficient b, compared to the linear model case. In dynamic simulations, this would yield higher persistence of transmission and a more rapid rate of human reinfection after mass drug administration (MDA) to the human population.
Thus, the basic inputs needed for running a coupled SWB model system are: (i) the biological SWB parameters of worm fecundity ρ(w) and eggrelease E _{;} (ii) human and snail population densities H, N, along with their demographic (birth, death, migration) parameters; and (iii) the human exposure/ water contact rate ω, and the resulting operative transmission parameters (α, b).
Calibrating the human SWB system
The first input for calibration is a human eggtest data set (~500 cases) from a given community, or a population subgroup (e.g. age group). The goal of this humanside SWB calibration is to find most likely values of biological (fecundity) parameters ρ _{ B } = {ρ _{0}, w _{0}, r} and transmission λ (= “rate of worm accumulation”/”worm mortality”) that are consistent with test data. Parameter λ is proportional to the mean worm burden of the SWB community in its equilibrium (endemic) state, \( \overline{w}=\lambda \varDelta w \).
For each choice (λ, ρ _{ B }), we simulate an ensemble of random eggtests outputs to compare with the observed data. Each simulated eggtest involves two random steps: (1) the random selection of subjects from a population subgroup to be tested and (2) the random eggrelease by each tested host in each stratum. These two steps are combined via a random sampling of a mixed NB distribution (Eq. 4) with prescribed parameters {λ, ρ _{ B } = (ρ _{0}, k _{0}, r)} (see Table 3). By this approach, we generate an SWB equilibrium {h _{ k }(λ)} from λ and then compute the mean eggrelease E _{ k } = ρ _{ k } ϕ _{ k } and aggregation r _{ k } = rϕ _{ k } for each stratum h _{ k }, using the determined biological parameters ρ _{ B }, i.e. \( {\rho}_k={\rho}_0{e}^{k/{k}_0} \) (fecundity) and Φ _{ k } (mating (Eq. 1)). Three inputs {h _{ k }, E _{ k }, r _{ k }} give rise to the mixed NBdistribution D _{ M } (λ, ρ _{ B }) of (Eq. 4). The simulated community eggtest output is thus a random sample of H subjects (the surveyed pool) drawn from distribution \( {D}_M,{\mathsf{E}}_T=\left\{{e}_1;\dots; {e}_H\right\} \) (e _{ i }  eggcount of ith test sample).
For further parameter fitting (calibration), individual counts are binned into eggcount distributions E _{ S } = {c _{0}, c _{1}, …} determined by prescribed sequence of Epartition bins: E _{0} = 0 < E _{1} = 1 < E _{2} < … < E _{ n }
Here, c _{0} counts negatives (uninfected pool), c _{1}  includes the range E _{1} ≤ e < E _{2}, etc. Different binning choices E _{0} < E _{1} < E _{2} < … are possible depending on the range and distribution of test diagnostics. For S. haematobium, 10 ml urine filtration [21] has typical counts 0 ≤ e ≤ 1000, and egg count values are highly overdispersed (Fig. 3). A suitable choice in such context is the logscale E _{ k } = 2^{k}(k = 1, 2, …, 10), such that the eggcount range, ([0–1,000), would split into 10 logscale bins.
The logbin counts of test data and of simulated tests show a typical bimodal pattern, with maximal value c _{0} (uninfected) and another peak between 50 < E < 500 (Fig. 4), depending on mean community burden. The elevated negative egg count category, c _{0}, could be an overestimate due to low test sensitivity for light infections [20, 21], or it could mean low infection prevalence among the tested subject pool.
Simulated egg distribution E _{ S }(λ, ρ _{ B }) (Eq. 10) depends on model parameters, but each output E _{ S } = {c _{0}, c _{1}, …} is random. The same, we expect, should hold for the real test data E _{ D } = {d _{ k }}, due to uncertainties of sampling and diagnostics [17–22]. Our goal is to compare two random samples (E _{ S }; E _{ D }) to assess their “proximity” via a suitable distancefunction, d(E _{ S }; E _{ D }). The task is confounded by randomness of both samples, particularly of simulated E _{ S }. To account for this, we assess the “proximity” between the simulated test ensemble \( \mathsf{E}=\mathsf{E}\left(\uplambda, {\uprho}_{\mathrm{B}}\right)=\left\{{\mathrm{E}}_{\mathrm{S}}\right\} \), and observed data, E _{ D }, instead of comparing individual test samples {E _{ S }}. Such a procedure is cast in a Bayesian framework by asking how likely it is to observe a given data set E _{ D } for a particular (parameter choice) ensemble \( \left(\mathsf{E}=\mathsf{E}\left(\lambda, {\rho}_B\right)\right) \). In other words, we ask what is the probability of observing E _{ D } conditioned on \( \mathsf{E}? \).
A natural answer can be given via the mean and covariance structure of ensemble \( \mathsf{E}\left(\lambda, {\rho}_B\right) \), {Ē(λ, ρ _{ B }), σ(λ, ρ _{ B })}. Specifically, we define the distance (error) function between \( \mathsf{E} \) and E _{ D } as:
Definition (Eq. 11) gives a family of likelihood weights on parameter space (λ, ρ _{ B })
with larger W corresponding to more likely parameter choices, consistent with the data E _{ D } (see Additional file 3 for technical details).
Calibration proceeds by randomly scanning the parameter space (λ, ρ _{ B }) (“uninformed” prior) and generating a meancovariance structure for uniformly sampled 4D hypercube in the parameter space: λ ' < λ < λ "; ρ _{0} ' < ρ _{0} < ρ _{0} "; k _{0} ' < k _{0} < k _{0} "; r ' < r < r ".
Once such a meancovariance testbed ({Ē(λ, ρ _{ B }), σ(λ, ρ _{ B })}) is computed, any specific test data (E _{ D }) will generate a family of likelihood weights (Eq. 12) on the (λ, ρ _{ B }) space and give the empirical posterior distribution \( \mathsf{D}=\mathsf{D}\left(\lambda, {\rho}_B\right) \) consistent with E _{ D }.
Figure 4 illustrates the range of likely egg output simulation runs using a fixed, random parameter choice that makes the Ē _{ S }, E _{ D } distribution fairly close.
The above calibration procedure is applied to the youngest age group (children). The adult calibration requires an additional input, namely a SWB  source term coming from maturing children (see Additional file 2). Thus, the adult parameter space has an additional calibrated parameter  preadult λ _{ C }, along with adult λ _{ A } and biological adult ρ _{ B } = (ρ _{0}, k _{0}, r). To estimate its likelihood weights, we proceed as above by scanning the extended parameter space (λ _{ C }, λ _{ A }, ρ _{ B }) and generating the adult testbed {Ē(λ _{ C }, λ _{ A }, ρ _{ B }), σ(λ _{ C }, λ _{ A }, ρ _{ B })}.
Once the meancovariance testbed is computed, calibration of any specific data set proceeds straightforwardly. Binned eggdata is substituted into the error function (Eq. 11) to get a family of likelihood weights (Eq. 12), resulting in an (empirical) posterior distribution consistent with the data.
Given a community (village) test data, we can generate two (child and adult) posterior distributions on their respective parameter spaces in 2 steps. First, the child group is calibrated to get its posterior and the resulting marginal distribution of child FOI (weights {W (λ _{ C })}). Then the adult posterior is computed using source terms derived from the child SWB strata {h ^{*}(λ _{ C })} evaluated at properly weighted FOI λ _{ C }, i.e. the distribution W (λ _{ C }) estimated earlier. The resulting adult likelihood weights in the extended parameter space are conditioned on the likely child values W (λ _{ C }), as
Calibrating the coupled humansnail system’s transmission coefficients
Having calibrated the human side of SWB dynamics, our next task is to combine human and snail infection data to estimate transmission coefficients (α, b) of coupled humansnail systems. In earlier work [14], we developed such a scheme for simplified SWB  snail systems. The current SWB modelling approach required significant modification, as outlined below.
First we demonstrate calibration for a single human SWB  single snail site, then proceed to more complex, coupled childadult systems. Equilibrium solutions of the snail system (Additional file 4) allow us to relate observed snail infection data to model parameters. Typical data include snail prevalences (prepatent  y, patent  z, total  N) expressed through basic inputs: growth rate β _{ 0 }, mortality \( \nu \), carrying capacity K, as well as snail FOI Λ (which needs to be estimated during calibration). Some unknown (uncertain) parameters, such as the local snail carrying capacity K (assumed here to be stationary), the humansnail contact rate ω, and snail susceptibility α, can be combined into single transmission coefficient. Specifically,

i)
\( A=\alpha \omega =\frac{\lambda /\varDelta w}{N\;z} \)  snailtohuman transmission

ii)
\( c=\frac{\nu }{r}\frac{z}{y} \)  patency conversion fraction (estimated from snail prevalences)

iii)
\( \varLambda =\frac{v\left(v+r\right){y}^{*}}{v\left(c\;r+v\right){y}^{*}} \)  snail FOI, estimated from snail prevalence data, known ν, r and estimated c

iv)
\( B=b\omega =\frac{N^{*}\left(\varLambda \right)}{N\;E}1\mathrm{n}\left(\frac{\varLambda_0}{\varLambda_0\varLambda}\right) \)  humantosnail transmission in snail FOI (Eq. 9).
Here N ^{*} (Λ) is total equilibrium snail density (see Additional file 4), and Λ_{0} is rate of sporocyst establishment in snails (Table 5).
We first estimate (c, Λ)  equations (iiiii), then apply them to transmission coefficients A, B via equations (iiv). Calibration of a mixed SWBsystem, made of several groups (i = 1, 2, …), requires additional assumptions on relative transmission rates B _{ i } /A _{ i }. Namely, B _{ i }/A _{ i } = b/α should be identical for all population groups i = 1, 2, …. To compute b/α we use estimated transmission coefficients \( {A}_i=\frac{\lambda_i}{N\;z} \) and replace human infectivity factor H E in equation (iv) by combined infectivity of all population groups ∑H _{ i } E _{ i }. Then
and humanto snail transmission by each group is given by \( {B}_i=\frac{b}{\alpha }{A}_i \).
Adding MDAbased control to the SWB system
The effect of drug treatment on stratified (SWB) population is to move a treated fraction of stratum h _{ n } (t) to a lowerlevel stratum h _{ m } (t), where m ≈ ε n is determined by the estimated efficacy of drug ε = fraction of adult worms surviving each drug treatment (see [14, 33]). In particular, all strata in the lowest range {h _{ m } : 0 ≤ m < 1/ε} shift to h _{ 0 } (effective clearing), the next interval {h _{ m } : 1/ε ≤ m < 2/ε} would go to h _{1}, etc. In numeric code, each drug treatment is simulated as an instantaneous event, due to the short duration of drug action (days) compared to the slow timescale of transmission dynamics (months to years). Computationally, terminal values of SWB variables at the treatment time t _{0} are reinitialized to new (posttreatment) values, depending on MDA inputs, the treatment coverage fraction (0 < f < 1), the drug efficacy, ε, etc. Each MDA “event” would then reshuffle variables {h _{ m }(t)} according to
The reinitialized system is solved over the prescribed timerange (between two “events”) and the process continues.For more detail, see Fig. 4 of reference [33].
Results
Model calibration
We applied our calibration scheme to infection data collected in the Msambweni region of Coastal Kenya [23–26], where repeated crosssectional surveys were conducted in 12 villages using standard filtration diagnostics (10 ml urine sample test [20, 21]) along with surveys of water contact [26] and local snail infection data [24]. The first rows of Table 6 (children) and Table 7 (adults) summarize basic demographics and epidemiological results of those surveys. For convenience, we ordered villages by their infection prevalence from highest risk (V _{1}) to lowest risk (V _{12}) based on initially observed childhood prevalence values.
For data analysis and model calibration, the total population of each village was split into children (0–20 years) and adults (20+); this choice was partly motivated by distinctive drop of infection about age 20 (Fig. 2). Additional inputs specific for Kenya are listed in Table 5. Two sets of calibrated parameters (Table 8) include agespecific fecundity (ρ, w _{0}, r) and human FOI λ. The final result of calibration was 24 posterior distributions, i.e. the 12 Msambweni villages, with two modeled age groups for each village (Tables 6 and 7).
Posterior distributions and their likelihood weights (Eq. 12) play important roles in our analysis and the control simulations reported below. All statistical outputs (means, correlations, quantiles, etc.) were computed relative to ensemble \( \mathsf{D} \). Thus for MDA simulations, we ran multiple treatment histories using a suite of likely parameter choices (for a given community) based on \( \mathsf{D} \), and assigned each output its respective likelihood weight.
A brief summary of calibration results for the 12 villages (ensemble mean and standard deviation, SD) is given in Table 6 (for children), Table 7 (for adults), and the accompanying Fig. 5 (showing ensemblemean biological parameters). Figure 6 shows the distributions of calibrated parameter for child and adult groups in high and low transmission villages.
Not unexpectedly, fecundity parameters exhibited fairly consistent values across the region with nearly horizontal linear fits for these parameter values across villages (Fig. 5). This suggested that hosts carrying comparable worm burden release similar eggcounts regardless of residing in high vs. low transmission areas. The only significant difference in parameter estimates comes between age groups, with children showing much higher, perworm fecundity than adults (ρ _{ C } ≫ ρ _{ A }). Consistent values of biological parameters across the region allow us to combine the 12 local village distributions of (ρ _{0}, w _{0}, r) into a single posterior ensemble, \( {\mathsf{D}}_B \). This “biological ensemble” \( {\mathsf{D}}_B \) is then used in the subsequent analysis and simulation of coupled humansnail systems. Calibrated child and adult FOIs (λ_{C} ; λ_{A}) show consistently higher values for children than for adults across the region (see Fig. 6).
Estimates for the densitydependent crowding effect
Our calibration results also provide estimates of the worm fecundity ‘crowding’ function ρ (w) of Eq. (2). The data on worm fecundity in human hosts are sparse and difficult to measure in vivo; however, our results give indirect estimates of ρ _{ B } = (ρ _{0}, w _{0}, r) as part of biological ensemble \( {\mathsf{D}}_B \). Ensemble envelopes of fecundity function (based on 1,000 random choices from biological posterior \( {\mathsf{D}}_B \)) are shown in Fig. 7. For these estimates, diagnostic test results (eggs per 10 ml urine) were adjusted for average daily urine release: U _{ C } ≈ 700 ml (for children) and U _{ A } ≈ 1300 (for adults) [34]. The estimated crowding effect was pronounced for both human age groups, with children’s subgroup wormfecundity dropping from a high of > 1,000 eggs/worm/day at very low burden (in children) to fewer than 100 eggs/worm/day at heavy burden. For adults, the estimates were 1,200 eggs/worm/day at low intensity, declining with halfvalue every ≈ 120 worms of burden. Overall patterns were consistent with known data (see, e.g. reference [9], chapter 15, and [29]). An alternative approach to modeling fecundity effects was included in our earlier work using a simpler SWB without inhost biology [14]. There, fecundity was assumed uniform across all strata and the resulting calibrated ρvalues were broadly distributed in the range [0, ρ _{0}]. The two models, simple SWB and the current refined version, differ in their predictions as described below.
Predicting prevalence and intensity curves
The SWB model predicts specific relations between prevalence and intensity based on its simulated eggtest results. Namely, for each parameter choice (λ, ρ _{0}, k _{0}, r), we can take the corresponding SWB equilibrium {h _{ k } (λ)} and compute its mating/fecundity factors ρ _{ k }, ϕ _{ k } (Eqs.1 and 2) in terms of biological parameters ρ _{ B } = (ρ _{0}, k _{0}, r). Then, from mixedNB assumption (Eq. 4) on eggrelease, we get formulae for prevalence P _{ E } and infection intensity M _{ E } as functions of parameters (λ, ρ _{ B }):
Each biological choice P _{ B } gives a particular (parametric) prevalenceintensity curve (Eq.13), by continuously varying λ. Whereas equations (Eq. 13) define an increasing function M _{ E } = f(P _{ E }), there is no simple analytic formula for it. Such curves can be computed and manipulated numerically, however. Fig. 8a compares the envelope of theoretical curves (Eq. 13) based on the biological posterior ensemble calibrated with Kenyan survey data [26]. Most data points for children lie within the 95 % quantile envelope, though we observe some departure (overprediction) at higher prevalence values.
Model validation
First, we generated an ensemble of virtual communities drawn from two calibrated posteriors. Then each group was randomly tested and its binned eggcounts recorded. The resulting ensemble of bincounts were compared to the test data of each group, (see Fig. 9 for children and adult subgroups), which produced reasonable agreement with the data overall. Model prevalence estimates were close to the observed data points (i.e. within the 95 % uncertainty envelopes) for both child and adult groups.
Another validation test involved a functional relation between eggtest prevalence P _{ E }, FOI λ or its inverse function λ(P _{ E }), and predicted mean infection intensity (= mean eggcount). Model equations (Eq. 13) gives both P _{ E }(λ, ρ _{ B }) and M _{ E }(λ, ρ _{ B }) as functions of λ and biological parameters ρ _{ B } = (ρ _{0}, w _{0}, r). To test consistency of such relations with village data across the region, we took a random selection of parameters {ρ _{ B }} drawn from the biological posterior \( {\mathsf{D}}_B \) and for each choice computed M _{ E }(λ, ρ _{ B }) and P _{ E }(λ, ρ _{ B }) curves and their envelopes, shown in Fig. 8a, b. We found the data points lying well within the margins of the predicted mean and prevalencefunctions. [N.B. These functional relations (Eq. 13) illustrated in Fig. 8b are useful for estimating an unknown SWB parameter λ in situations where the available data are incomplete, e.g. when only aggregated, communitylevel egg prevalence data are available instead of individuallevel egg test data].
Projecting the impact of MDAbased control
For this step, we applied our calibrated model to project MDAcontrol effects and posttreatment prevalence in one of the twelve Msambweni villages, Milalani, which was more intensively studied for the impact of MDA. The allages data used for analysis were collected in three villagewide surveys from 2000 (pretreatment baseline), 2003 (2 years posttreatment) and 2009 (six years posttreatment). The communitywide MDA in 2000 had a 79 % treatment coverage whereas the 2003 treatment coverage was only 41 %. The baseline (pretreatment) data served to calibrate the system, while the two posttreatment data sets were used to test model predictions. The coupled humansnail model consisted of two SWB groups (children, adults) linked to a single hypothetical snail site with a snail infection level approximating the average of five known Msambweni snail sites; SEI snail prevalence values in this experiment were taken as {x*, y*, z*} = {.63,.35,.02} based on field observations of bulinid snail PCR positivity for S. haematobium DNA [35] used to determine the number of ‘exposed’ snails and the observed annual frequencies of shedding snails [24, 36] for ‘infectious’ snails.
For our long term MDA prediction, we also scaled in overall population growth based on Kenyan demographics [37] and interseasonal variations of snail density over the 9year study period. An ensemble of 150 likely calibrated parameter choices was drawn from the Milalani posterior and their 9year histories computed. At each time, t, we took dynamic solutions and used equations (Eq. 13) to estimate the corresponding community eggtest results (prevalence and intensity) for child and adult age groups. We also ran this simulation for several values of the snail FOI parameter rate, Λ_{0} (sporocyst establishment rate), ranging from 0.5 to 2 weeks.
The results shown in Fig. 10 correspond to 1/Λ_{0} = 1.5 week. We plotted the ensemble envelope along with predicted mean and compared to the three observed data points for prevalence and intensity, respectively. The data for 2000–2003 fell well within predicted margins. In modeling treatment outcomes after MDA in this test village, our prevalence and intensity estimates proved to be somewhat low by the last study year (Fig. 10). Infection intensity (measured by mean egg test) by year 9 is underpredicted for children and overpredicted for adults. This could be attributed to demographic and/or environmental changes that happened over 6year intervening period, which were not accounted for in the model.
Comparison between simple and advanced SWB models: parameter estimation and control predictions
The key for analysis and calibration of coupled humansnail systems are two forces of infection: snailtohuman λ and humantosnail Λ. Both depend on snail/ human infectivity and the population densities of host and vector. For simple SI (susceptibleinfected) snail system, λ is proportional to infected (patent) snail prevalence λ ∝α y with transmission coefficient a = “mean rate of worm accumulation in host” and infected (patent) snail prevalence y.
In the simple SWB system, snail FOI is taken proportional to mean human ‘infectivity’ = mean eggrelease E by human hosts into the environment, Λ = bE, with transmission coefficient b. Transmission coefficients \( \alpha \), b combine multiple factors and processes (population densities, humansnail contact rates, intermediate larval stages, etc.). Human infectivity function E depends on details of on worm aggregation and inhost biology.
For the earlier, simpler SWB, the calibration procedure [14, 15] employed algebraic relations between test data (prevalence P _{ D } and mean intensity E _{ D }), model functions P _{ W }(λ) = 1 − h _{0}(λ), and MWB \( \overline{w}\left(\lambda \right)=\varDelta w{\displaystyle {\sum}_k{h}_k\left(\lambda \right)} \) expressed through dimensionless human FOI \( \lambda =\frac{\alpha y}{\gamma \varDelta w} \). No distinction was made between worm and eggtest prevalence (P _{ W } = P _{ E }) as eggrelease was assumed to be proportional to MWB \( \left(E=\rho \overline{w}\right) \) with fixed (uniform) fecundity/worm factor ρ to be estimated. That calibration proceeded in 2 steps: from prevalence equation P _{ W }(λ) = P _{ D }, we computed equilibrium λ ^{*} then derived MWB w ^{*} and the fecundity factor ρ
with mating function ϕ(w) = w/2. The transmission coefficients were estimated as
This procedure could be extended to demographically and/or geographically linked SWB systems (see [14, 15]). No uncertainties entered the simple SWB (or its calibration), but explicit algebraic formulae allowed one to relate data noise to parameter estimation.
In comparing the predicted FOI curves (P _{ E }(λ, ρ _{ B }) and P _{ W }(λ)) for the two SWB systems (old and new), only the latter, P _{ W }, exists for the simple SWB, as no distinction was made there between wormbased and eggtestbased prevalence. For the new SWB, Fig. 8b shows quantile envelopes of P _{ E } curves sampled over a range of biological parameters ρ _{ B } (biological posterior). In general, we expect P _{ E }(λ) < P _{ W }(λ), as wormcarrying strata could contribute to zero egg count. For the simpler SWB, the red curve in Fig. 8b shows its significant departure from the median (and surrounding quantile envelopes) of the new SWB that incorporated hostworm biological interactions. We conclude that the old SWB significantly underestimated FOI λ ^{*} and transmission rate a (Eq. 15), while overestimating fecundity factor ρ (Eq. 14).
To explore possible effect of such discrepancies on MDA control prediction, we took the above 9year study for Milalani and compared simple SWB (dashed red line) with new calibrated envelope prediction in Fig. 10. Predicted prevalence of both models is reasonably close (within uncertainty envelopes), but the predicted infection intensity for children is underestimated using the simpler SWB (the dashed red line outside the newly calibrated envelope of likely values).
Discussion
The newer calibrated SWB modeling approach appears well suited for simulating the effects of control interventions, particularly the effects of mass drug therapy. Other data inputs and other control strategies could easily be implemented within our setup, such as targeting treatment to specific populations, and modeling the efficiency of different posttreatment surveillance strategies. In the present analysis, we examined one case of longterm MDA outcomes in a Kenyan study site (population ~2,000) to validate our model in a dynamic setting and assess its predictive accuracy. Predicted relaxation patterns (means, envelopes, and quantiles) were found in good agreement with the observed data. Additional evidence comes from recent work [33] where our model and methodology was applied to program data from the Schistosomiasis Consortium for Operational Research and Evaluation (SCORE) on S. haematobium control from Mozambique. While this Mozambique data covered a wide spectrum of communities and control strategies, it was limited in terms of the number and age range of persons tested. To fill in some of the missing data gaps, we utilized certain parameter estimates from the Kenyan data studied here, assuming that similar population groups in both countries share comparable ‘biological’ infection parameters in terms of worm fecundity and density dependence. These calibrated inputs allowed us to build coupled SWBsnail systems for Mozambique communities and accurately simulate the MDA outcomes of the SCORE project.
Typical diagnostic tests for Schistosoma infection (based on egg counting in stool or urine) exhibit highly uneven distributions in host populations. These heterogeneities exist not only for broader communities, but for specific demographic groups that are assumed to be nearly the same in terms of risk. We conclude that several factors contribute to overdispersion of test results, among them uneven worm load due to varying exposure and/or host susceptibility [6], and irregular (clustered) egg release into human excreta [16]. Conventional populationbased approaches to modeling (mean worm burden (MWB) models [8]) either ignore uneven burden (taking its population mean) or impose adhoc assumptions on worm distribution (e.g. the negative binomial [12]). Either approach has severe limitations (see [12]) that can reduce the utility of the model and the accuracy and robustness of its predictions. While autonomous agent, individualbased model simulations can allow for multiple heterogeneities [15], this alternative type of model has limited capacity in term of population size and program implementation on desktop/laptop platforms.
The SWB approach bridges the gap between these two types of models. It gives a consistent, assumptionfree account of uneven worm load and naturally accommodates essential inhost biology, including worm mating probability and densitydependent reduction in fecundity. While the resulting SWB systems have more variables (depending on stratification), there are only a few model parameters per stratum that need to be calibrated, similar to the requirements for a lowdimensional MWB model. We note that promiscuous mating by female worms would enhance the continuation of transmission following MDA if, as generally seen, only partial elimination of worms is achieved with treatment. We have discussed the issue of mating patterns in some depth in our previous paper [12].
To better mimic the dynamics of humantosnailtohuman parasite transmission, we have revisited and revised conventional approaches to modeling snail population and infection in our coupled SWB modeling systems. Typically, in MWB and related models, the parasite’s shortlived larval stages are not modeled but rather incorporated into an effective ‘force of infection’ (FOI) term, λ, for snailtohuman transmission, and Λ, for humantosnail transmission. The latter, Λ, is often taken be proportional to humantosnail infectivity in terms of cumulative Schistosoma egg release by the local human community, with a fraction of these eggs converting to miracidia to invade susceptible snails (see e.g. [9]). A new, closer look at the miracidium snail invasion process reveals a more likely nonlinear (saturated) form of snail FOI and sporocyst establishment [32]. The resulting estimates of humantosnail transmission could provide an explanation for markedly different reinfection rates in some postMDA communities [38] and the leveraged impact of human inmigration on persistence of transmission [39–41].
To better project future control program outcomes in terms of presentday data, we have developed a Bayesian calibration procedure for our model based on simulating the recognized imperfections of diagnosing infection based on eggcount data [17–22]. Other factors believed to be critical in accurate forecasting of Schistosoma prevalence include human agegroup differences in exposure and susceptibility to infection [42–44], densitydependent (crowding) effects on egg production per worm [30, 45], limitations on mating success among adult schistosomes [4, 46], and the development of antifecundity immunity among older patients [47]. In our current SWB model, which included these factors [12], we found better calibration for the model in terms of projecting treatment and reinfection outcomes. The result of our calibration procedure is a posterior ensemble of likely SWB community values consistent with a given data set. The posterior distribution of parameter values for a given SWB community/group can be used to generate multiple community actualizations (based on (λ, ρ _{ B })  choices) in order to simulate the range of likely outcomes. Each outcome is assigned a significance level determined by its likelihood weight. Thus any data/model uncertainties are propagated into “prediction uncertainty”. In most dynamic simulations, e.g. MDA control, these uncertainties can then be shown as prediction envelopes of possible outcomes.
We have applied the above calibration scheme to a specific data set for communities in coastal Kenya, where only S. haematobium is endemic. While the Kenyan communities differed markedly in terms of their risk and infection levels, we found their agespecific biological parameters ρ _{ B } confined within the same close range, regardless of transmission intensity. This result supported our hypothesis on the constancy of inhost worm biology (mating, fecundity) and its parameterization. It was notable that the estimated crowding function (densitydependent fecundity) was different for child and adult groups, which appears to be in accordance with recent findings about the acquisition of antifecundity immunity in primates and humans [47]. Of note, our calibrated aggregation parameter r (Table 6) is also consistent with estimates of Hubbard et al. [16] for Schistosoma japonicum infection. A limitation of the present paper is its focus on S. haematobium (and its transmission features) in calibration of model predictions. However, work is in progress to repeat calibration and testing for S. mansonicontrol projects in Kenya and Uganda, which should allow comparison of the estimated biological parameters for each species and help to determine if model recalibration is necessary for different species and for different ecological settings.
Conclusions
The SWB provides an efficient, flexible, and viable approach for modeling Schistosoma transmission and control among stratified populations in simple and complex environments. SWB allows for inclusion of inhost biological factors and limitations of diagnostics, and is applicable to a broad range of treatment strategies. Most helpful to program managers, these new features allow us to predict diagnostic eggtest results for modeled SWB population subgroups and for communities atlarge. Where only partial diagnostic data are available, the curves in Fig. 8 can serve to estimate parameters of transmission for program outcomes predictions.
This work is being extended to treatment projections for largescale treatment trials currently implemented in both S. mansoni and S. haematobiumendemic areas. For the near future, as part of the ongoing NTD Modelling Consortium project [48], our refined SWB model will be further validated against a new data and directly compared to the more traditional deterministic model of our consortium partners [49]. Fitted model predictions will be compared for likelihood and precision using two large data sets from recent control programs in subSaharan Africa: (i) the 2003–2006 S. mansoni data from the Ugandan National Schistosomiasis Control Programme in the African Great Lakes region, and (ii) 2010–2015 data from the SCORE/SCI multivillage operational research trial on S. haematobium control in Mozambique, evaluated in our previous paper [33]. Overall, we expect that the methodology developed in the current paper has a broad scope of applications for different Schistosoma species and more generally, for helminth/macroparasites infections where inhost biology plays an important role.
Abbreviations
FOI, force of infection; MDA, mass drug administration; MWB, mean worm burden model; NB, negative binomial distribution; NTD, neglected tropical diseases; PDF, probability density function; SD, standard deviation of the mean; SWB, stratified worm burden model
References
 1.
Hotez PJ, Alvarado M, Basanez MG, Bolliger I, Bourne R, Boussinesq M, et al. The global burden of disease study 2010: interpretation and implications for the neglected tropical diseases. PLoS Negl Trop Dis. 2014;8(7), e2865. doi:10.1371/journal.pntd.0002865.
 2.
Savioli L, Daumiere D. Accelerating work to overcome the global impact of neglected tropical diseases: A roadmap for implementation. Geneva: World Health Organization; 2012.
 3.
Colley DG, Bustinduy AL, Secor WE, King CH. Human schistosomiasis. Lancet. 2014;383:2253–64. doi:10.1016/S01406736(13)619492.
 4.
May RM. Togetherness among schistosomes: Its effects on the dynamics of infection. Math Biosci. 1977;35:301–43.
 5.
Nåsell I, Hirsch WM. The transmission dynamics of schistosomiasis. Commun Pure Appl Math. 1973;26(4):395–453.
 6.
Shaw DJ, Grenfell BT, Dobson AP. Patterns of macroparasite aggregation in wildlife host populations. Parasitology. 1998;117(Pt 6):597–610.
 7.
Bavia ME, Hale LF, Malone JB, Braud DH, Shane SM. Geographic information systems and the environmental risk of schistosomiasis in Bahia, Brazil. Am J Trop Med Hyg. 1999;60(4):566–72.
 8.
MacDonald G. The dynamics of helminth infections, with special reference to schistosomes. Trans R Soc Trop Med Hyg. 1965;59(5):489–506.
 9.
Anderson RM, May RM. Infectious diseases of humans. Dynamics and control. New York: Oxford University Press; 1991.
 10.
Crofton HD. A quantitative approach to parasitism. Parasitology. 1971;62:179–93.
 11.
Bradley DJ. Regulation of parasite populations. A general theory of the epidemiology and control of parasitic infections. Trans R Soc Trop Med Hyg. 1972;66(5):697–708.
 12.
Gurarie D, King CH. Population biology of Schistosoma mating, aggregation, and transmission breakpoints: More reliable model analysis for the endgame in communities at risk. PLoS One. 2014;9(12), e115875. doi:10.1371/journal.pone.0115875.
 13.
Cornell SJ. Modelling stochastic transmission processes in helminth infections. In: Michael E, Spear RC, editors. Modelling parasite transmission and control. New York: Springer; 2010. p. 66–78.
 14.
Gurarie D, King CH, Wang X. A new approach to modelling schistosomiasis transmission based on stratified worm burden. Parasitology. 2010;137(13):1951–65.
 15.
Wang X, Gurarie D, Mungai PL, Muchiri EM, Kitron U, King CH. Projecting the longterm impact of school or communitybased masstreatment interventions for control of Schistosoma infection. PLoS Negl Trop Dis. 2012;6(11), e1903. doi:10.1371/journal.pntd.0001903.
 16.
Hubbard A, Liang S, Maszle D, Qiu D, Gu X, Spear RC. Estimating the distribution of worm burden and egg excretion of Schistosoma japonicum by risk group in Sichuan Province, China. Parasitology. 2002;125(Pt 3):221–31.
 17.
Carabin H, Budke CM, Cowan LD, Willingham 3rd AL, Torgerson PR. Methods for assessing the burden of parasitic zoonoses: echinococcosis and cysticercosis. Trends Parasitol. 2005;21(7):327–33.
 18.
de Vlas SJ, Engels D, Rabello AL, Oostburg BF, Van Lieshout L, Polderman AM, et al. Validation of a chart to estimate true Schistosoma mansoni prevalences from simple egg counts. Parasitology. 1997;114(Pt 2):113–21.
 19.
de Vlas SJ, Gryseels B. Underestimation of Schistosoma mansoni prevalences. Parasitol Today. 1992;8(8):274–7.
 20.
Savioli L, Hatz C, Dixon H, Kisumku UM, Mott KE. Control of morbidity due to Schistosoma haematobium on Pemba Island: egg excretion and hematuria as indicators of infection. Am J Trop Med Hyg. 1990;43:289–95.
 21.
Warren KS, Arap Siongok TK, Hauser HB, Ouma JH, Peters PAS. Quantification of infection with Schistosoma haematobium in relation to epidemiology and selective population chemotherapy. I. Minimal number of daily egg counts in urine necessary to establish intensity of infection. J Infect Dis. 1978;138:849–55.
 22.
Mwinzi PN, Kittur N, Ochola E, Cooper PJ, Campbell Jr CH, King CH, et al. Additional evaluation of the PointofContact Circulating Cathodic Antigen assay for Schistosoma mansoni infection. Front Public Health. 2015;3:48. doi:10.3389/fpubh.2015.00048.
 23.
Clennon JA, Mungai PL, Muchiri EM, King CH, Kitron U. Spatial and temporal variations in local transmission of Schistosoma haematobium in Msambweni, Kenya. Am J Trop Med Hyg. 2006;75(6):1034–41.
 24.
Kariuki HC, Clennon JA, Brady MS, Kitron U, Sturrock RF, Ouma JH, et al. Distribution patterns and cercarial shedding of Bulinus nasutus and other snails in the Msambweni area, Coast Province, Kenya. Am J Trop Med Hyg. 2004;70(4):449–56.
 25.
King CH, Blanton RE, Muchiri EM, Ouma JH, Kariuki HC, Mungai P, et al. Low heritable component of risk for infection intensity and infectionassociated disease in urinary schistosomiasis among Wadigo village populations in Coast Province, Kenya. Am J Trop Med Hyg. 2004;70(1):57–62.
 26.
Muchiri EM, Ouma JH, King CH. Dynamics and control of Schistosoma haematobium transmission in Kenya: an overview of the Msambweni Project. Am J Trop Med Hyg. 1996;55(5 Suppl):127–34.
 27.
Wilson RA, van Dam GJ, Kariuki TM, Farah IO, Deelder AM, Coulson PS. The detection limits for estimates of infection intensity in schistosomiasis mansoni established by a study in nonhuman primates. Int J Parasitol. 2006;36(12):1241–4.
 28.
Nasell I. Mating for schistosomes. J Math Biol. 1978;6(1):21–35.
 29.
Anderson RM, Medley GF. Community control of helminth infections of man by mass and selective chemotherapy. Parasitology. 1985;90(Pt 4):629–60.
 30.
Medley G, Anderson RM. Densitydependent fecundity in Schistosoma mansoni infections in man. Trans R Soc Trop Med Hyg. 1985;79(4):532–4.
 31.
Crews AE, Yoshino TP. Schistosoma mansoni: effect of infection on reproduction and gonadal growth in Biomphalaria glabrata. Exp Parasitol. 1989;68(3):326–34.
 32.
Thiele EA, Minchella DJ. Molecular assessment of trematode coinfection and intraspecific competition in molluscan intermediate hosts. Mol Biochem Parasitol. 2013;187(1):52–9. doi:10.1016/j.molbiopara.2012.12.003.
 33.
Gurarie D, Yoon N, Li E, NdeffoMbah M, Durham D, Phillips AE, et al. Modelling control of Schistosoma haematobium infection: predictions of the longterm impact of mass drug administration in Africa. Parasit Vectors. 2015;8(1):529. doi:10.1186/s1307101511443.
 34.
Hesse A, Classen A, Knoll M, Timmermann F, Vahlensieck W. Dependence of urine composition on the age and sex of healthy subjects. Clin Chim Acta. 1986;160(2):79–86.
 35.
Hamburger J, Hoffman O, Kariuki HC, Muchiri EM, Ouma JH, Koech DK, et al. Largescale, polymerase chain reactionbased surveillance of Schistosoma haematobium DNA in snails from transmission sites in coastal Kenya: A new tool for studying the dynamics of snail infection. Am J Trop Med Hyg. 2004;71:765–73.
 36.
Sturrock RF, Kinyanjui H, Thiongo FW, Tosha S, Ouma JH, King CH, et al. Chemotherapybased control of schistosomiasis haematobia. 3. Snail studies monitoring the effect of chemotherapy on transmission in the Msambweni area, Kenya. Trans R Soc Trop Med Hyg. 1990;84(2):257–61.
 37.
The World Bank. World development indicators. 2015.
 38.
Kahama AI, Vennervald BJ, Kombe Y, Kihara RW, Ndzovu M, Mungai P, et al. Parameters associated with Schistosoma haematobium infection before and after chemotherapy in school children from two villages in the coast province of Kenya. Trop Med Int Health. 1999;4(5):335–40.
 39.
Fenwick A, Jorgensen TA. The effect of a control programme against Schistosoma mansoni on the prevalence and intensity of infection on an irrigated sugar estate in northern Tanzania. Bull World Health Organ. 1972;47(5):579–86.
 40.
Macdonald F, Clarke Vde V, Gaddie P, Atkinson G. Report on a largescale attempt at control of bilharziasis by combined mass treatment and intensive snail control. Cent Afr J Med. 1973;19:22–32.
 41.
Gautret P, Mockenhaupt FP, von Sonnenburg F, Rothe C, Libman M, Van De Winkel K, et al. Local and international implications of schistosomiasis acquired in Corsica, France. Emerg Infect Dis. 2015;21(10):1865–8. doi:10.3201/eid2110.150881.
 42.
Abel L, Demenais F, Prata A, Souza AE, Dessein A. Evidence for the segregation of a major gene in human susceptibility/resistance to infection by Schistosoma mansoni. Am J Hum Genet. 1991;48:959–70.
 43.
Barbour AD. The importance of age and water contact patterns in relation to Schistosoma haematobium infection. Trans R Soc Med. 1985;79:151–3.
 44.
Chandiwana SK, Woolhouse ME. Heterogeneities in water contact patterns and the epidemiology of Schistosoma haematobium. Parasitology. 1991;103(Pt 3):363–70.
 45.
Polman K, De Vlas SJ, Van Lieshout L, Deelder AM, Gryseels B. Evaluation of densitydependent fecundity in human Schistosoma mansoni infections by relating egg counts to circulating antigens through Deming regression. Parasitology. 2001;122(Pt 2):161–7.
 46.
May RM, Woolhouse ME. Biased sex ratios and parasite mating probabilities. Parasitology. 1993;107(Pt 3):287–95.
 47.
Wilson S, Jones FM, van Dam GJ, Corstjens PL, Riveau G, Fitzsimmons CM, et al. Human Schistosoma haematobium antifecundity immunity is dependent on transmission intensity and associated with immunoglobulin G1 to wormderived antigens. J Infect Dis. 2014;210(12):2009–16. doi:10.1093/infdis/jiu374.
 48.
Hollingsworth TD, Adams ER, Anderson RM, Atkins K, Bartsch S, Basanez MG, et al. Quantitative analyses and modelling to support achievement of the 2020 goals for nine neglected tropical diseases. Parasit Vectors. 2015;8:630. doi:10.1186/s1307101512351.
 49.
Anderson RM, Turner HC, Farrell SH, Yang J, Truscott JE. What is required in terms of mass drug administration to interrupt the transmission of schistosome parasites in regions of endemic infection? Parasit Vectors. 2015;8:553. doi:10.1186/s130710151157y.
 50.
Anderson RM, May RM. Herd immunity to helminth infection and implications for parasite control. Nature. 1985;315(6019):493–6.
Acknowledgements
We thank our colleagues R. Alsallaq, A. Galvani, D. Durham, M. Ndeffo Mbah, and E. Gurarie for useful inputs and stimulating discussion.
Funding
This work was funded by the Schistosomiasis Consortium for Operational Research and Evaluation (SCORE), based at the University of Georgia, USA, and by the Children’s Investment Fund Foundation (UK) (“CIFF”) through a grant to the Neglected Tropical Diseases Modelling Consortium at Warwick University, UK. The views, opinions, assumptions or any other information set out in this study are solely those of the authors and should not be attributed to CIFF or any person connected with CIFF. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Availability of data and materials
The data used in this study for model calibration and validation are included as Additional file 5.
Authors’ contributions
DG, CHK, NY, EL conceived of the study, and participated in its design, performance, and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
This project reused anonymized human data from published research surveys [23, 25]. Those surveys were performed with written informed consent of the participants under supervision of the University HospitalsCase Medical Center IRB (protocol 110742), Cleveland, Ohio, and the Kenya Medical Research Institute Ethical Review Committee (protocol nonSSC #087), Nairobi. Kenya.
Author information
Additional files
Additional file 1:
SWB codes15.nb [an example of the SWB model implemented using Mathematica v10.2 software (Wolfram Research, Champaign IL. USA)] (NB 35 kb)
Additional file 2:
Mixed SWB systems and equilibria (DOCX 90 kb)
Additional file 3:
Likelihood estimates for simulated egg test results (DOCX 204 kb)
Additional file 4:
Snail equilibria and calibration (DOCX 64 kb)
Additional file 5:
Infection data from Milalani village crosssectional surveys (XLSX 55 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Mathematical models
 Schistosomiasis control
 Drug therapy
 Disease transmission
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate. Please note that comments may be removed without notice if they are flagged by another user or do not comply with our community guidelines.