Population Biology of Schistosoma Mating, Aggregation, and Transmission Breakpoints: More Reliable Model Analysis for the End-Game in Communities at Risk

Mathematical modeling is widely used for predictive analysis of control options for infectious agents. Challenging problems arise for modeling host-parasite systems having complex life-cycles and transmission environments. Macroparasites, like Schistosoma, inhabit highly fragmented habitats that shape their reproductive success and distribution. Overdispersion and mating success are important factors to consider in modeling control options for such systems. Simpler models based on mean worm burden (MWB) formulations do not take these into account and overestimate transmission. Proposed MWB revisions have employed prescribed distributions and mating factor corrections to derive modified MWB models that have qualitatively different equilibria, including ‘breakpoints’ below which the parasite goes to extinction, suggesting the possibility of elimination via long-term mass-treatment control. Despite common use, no one has attempted to validate the scope and hypotheses underlying such MWB approaches. We conducted a systematic analysis of both the classical MWB and more recent “stratified worm burden” (SWB) modeling that accounts for mating and reproductive hurdles (Allee effect). Our analysis reveals some similarities, including breakpoints, between MWB and SWB, but also significant differences between the two types of model. We show the classic MWB has inherent inconsistencies, and propose SWB as a reliable alternative for projection of long-term control outcomes.


Introduction
In the last decade, greater recognition of the sub-clinical, but disabling effects of schistosomiasis has led to a new awareness of the importance of preventing Schistosoma infection and reinfection among populations at risk [1]. Because of a better understanding of the long-term consequences of the chronic inflammation triggered by anti-Schistosoma immunity [2,3], it is no longer considered sufficient to provide just 'morbidity control' (via suppression of infection intensity), as has been advocated in the past [4,5]. Instead, it has become a priority to find practical means to interrupt transmission and provide local elimination of infection wherever possible. This objective has been outlined in the 2012 London Declaration for Neglected Tropical Diseases (NTDs, http://unitingtocombatntds. org/resource/london-declaration) and in the World Health Organization's 2020 Roadmap on NTDs [6], and codified in World Health Assembly resolution 65.21. These initiatives now seek tools to aid the goal of local elimination of schistosomiasis. As a result, it seems appropriate to re-evaluate existing transmission models of dioecious macroparasites (like Schistosoma) for their usefulness in implementing effective policy in areas that will experience declining human and intermediate host infection prevalence under the pressure of current infection-and transmission-control interventions.
In 1965, the dynamic modeling of MacDonald [7] prompted hopes there could be an ecological 'breakpoint' in Schistosoma spp. transmission if the local numbers of intermediate host snails could be reduced by 90% or more [8][9][10]. This local extinction was projected as a consequence of the obligate need for sexual reproduction by the parasite (dioecy) within the human host; if male and female worms could not combine within the same host, then transmission was effectively ended. MacDonald's analysis projected that there would be special leverage in obtaining reductions in transmission by interventions that limit snail-to-human transmission [7]. Such reductions were seen to be potentially achievable by existing modalities of chemical mollusciciding and/or snail habitat destruction. Effectively, the breakpoints described in MacDonald's work [7] and subsequent studies [11,12], are points or regions within the transmission parameter space below which existing worm burden is unsuccessful in maintaining transmission, and parasite numbers ultimately go to zero without further intervention.
Questions remain: do breakpoints exist in real world settings?, and does this phenomenon have relevance, i.e., yield a practical benefit in the context of community-wide Schistosoma control campaigns? In practice, where prolonged and extensive reductions of snail numbers were achieved [8,[13][14][15][16], Schistosoma prevalence often dropped, but transmission was not interrupted, suggesting flaws in the basic assumptions in the MacDonald model [7,17], and additional complicating factors with regard to parasite mating patterns [12,18]. In addition, heterogeneities in water habitat distribution and in human water contact behavior [19,20], were suspected to contribute strongly to the persistence of Schistosoma within suitable ecosystems.
Bradley and May [18] point out that the typically observed 'clumping' of high worm burdens among a small fraction of human hosts (overdispersion with aggregation) could lead, overall, to a lower stability of transmission. However, with greater aggregation, if male and female worms are transmitted in roughly equal numbers to each human host (or if worm mating is promiscuous and not monogamous) [21], then the projected breakpoint phenomenon might not prove relevant, because females are increasingly likely to be successfully mated in such scenarios. As a consequence, egg output, and hence human-to-snail transmission, will persist, albeit at lower levels. They note, however, that if the acquisition of male worms and of female worms is aggregated in separate fashion for each sex (as might occur with very low worm burdens in a low-transmission environment) then breakpoints are more likely to be relevant, and spontaneous failure of parasite transmission is more likely to occur.
To evaluate the likely relevance of the breakpoint phenomenon in current control efforts, the present study compares the projections of two established modeling approaches to Schistosoma transmission: i) the modified MacDonaldtype Mean Worm Burden (MWB) model proposed by May and colleagues [12,18] and by Nåsell and colleagues [11,17] utilizing their assumed negative binomial (NB) or Poisson distributions of worm burden; and ii) a distributionfree Stratified Worm Burden (SWB) approach we have previously developed [22,23], and now modify to include the effects of parasite mating probability (see Table 1 for a list of abbreviations used in this paper).
Of special interest is each model system's projections (and their limitations) when parasite burdens get very low. In brief, we find that the classical MWB approach and its extensions have shortcomings that limit their usefulness in projecting elimination in various transmission settings (see S1 Appendix in S1 Text,). In prior work [23], we have advocated the use of SWB systems as an improved alternative to the more analytically tractable (but potentially less realistic) MWB systems [24]. While SWB are high-dimensional models (depending on the number of strata), they can now be efficiently implemented, simulated, and studied numerically.
Our current approach allows a more detailed account of density-dependent mating factors-we are now able to include the so-called Allee effect (common for many species [25]), in which growth and subsequent mating success are significantly impaired when population numbers get very low in a given withinhost environmental 'patch'. For schistosomes, we propose that the Allee effect obtains when female worms fail to mature in the absence of sufficient males [26], leading to disproportionately lower transmission despite persistent (albeit low) mean worm burden in human hosts. Additionally, we indicate how the SWB approach can be extended to accommodate highly influential geographic and agerelated demographic factors [19,20].
Increased k gives a more aggregated (clumped) distribution for h n w,k ð Þ f g; in the limiting case where k~?, the NB distribution becomes the Poisson distribution.
Overdispersed parasite burden has been noted in many wildlife populations, and the NB has been proposed as the best model for this phenomenon [27][28][29]. In most cases, the apparent aggregation factor was relatively low, but found to vary widely among different species and locations. Despite its resemblance to empirical data, there is no biological dictate for choosing the NB distribution based on first principles. A multitude of biological, environmental and other factors can affect parasite distributions within definitive hosts, and the only justifiable pattern derived from the underlying principles (random worm acquisition), has been the Poisson case advocated by Nåsell & Hirsch [17].
The NB assumption has been widely used in modeling studies of macroparasite transmission [7,11,12,[29][30][31]. In these works, a prescribed distribution of worm burden (NB or Poisson), has been used to reduce a large (infinite-dimensional) stratified system of h n f g to a low-dimensional (MacDonald-type) ''moment system'' for the MWB variable w t ð Þ, and/or higher moments (reviewed in S2 Appendix in S1 Text). The reduced models, unlike the SWB [23], can be analyzed mathematically [24]. On the other hand, the SWB approach requires no a priori assumptions on distribution or aggregation of strata h n f g. Both SWB values arise naturally from the underlying processes of worm acquisition and loss. In future, our very practical interest will be to apply the calibrated SWB systems to demographically-or geographically-structured populations resembling problem areas that confront elimination program planners [32].
For such a complex, extended community, each population group can be represented by its own SWB system (Fig. 1), and these separate systems coupled via suitable source parameters. In Fig. 2, SWB equilibrium distributions are compared to the standard NB/Poisson case. The simplest 'single population' SWB system ( Fig. 1) has three main parameters: l, human FOI, balanced by worm mortality, c, and demographic loss, m. Its equilibrium distribution h n l',m' ð Þ f g depends on rescaled values l'~l=c and m'~m=c, the former l' having dimension  Comparison of worm distribution patterns for a negative binomial-based MWB system vs. a SWB system. Here, all SWB distributions (right-hand panels) are produced from an uninfected source S n f g. (a) NB-MWB with fixed mean w~10 and increasing k (k~? is the limiting Poisson case); (b) Equilibrium SWB distribution for uninfected source with FOI (mean) l~10, and varying demographic parameter m (see equation (13)); m plays the role of aggregation k for NB, with small m corresponding to large (infinite) k. Panels (c2d) Poisson distributions with different means, m i (left panel) vs. SWB distributions with different l i~mi (right panel) exhibit striking similarity. Panels (e)2(f) compare infectious prevalence V for the two models as a function of NB-MWB w or SWB l. of [worm burden] like the MacDonald-MWB term w, while m' is dimensionless. The resulting SWB solutions depend on demographic source terms (equations (13) below).
For better comparison to the basic NB/Poisson distributed models, we can use a simplified SWB system with only an uninfected source, i.e., a source term adding to the uninfected strata, h 0 , only (S 0 w0, S n~0 for n.0). Such sources would obtain in closed (isolated) populations, and be relevant to the youngest age-group (a newborn source) in the studied population. In general, equilibrium SWB distributions have no analytic formulae, so most results below are computed numerically. The only analytically tractable case corresponds to the limiting (degenerate) system, m~0, S~0 (no population turnover and sources). Here equilibrium solution gives the standard Poisson distribution h n l,0 ð Þ~e {l l n n! , consistent with Nåsell & Hirsch [17].

Adding mating patterns as functions in the extended MWB and SWB
A simple way to account for mating behavior of adult worms is via a pairing function h i, j ð Þ5the number of mated (egg-shedding) females for i -males, jfemales. Several studies have looked the effect of mating on snail infection in the MWB type models [11,12]. Examples for possible types of mating included: The mating patterns will affect transmission by determining the number of mated females, their effective egg production, and, as a consequence, the resulting force of snail infection, L. Specifically, for a given sex distribution P i, jjn ð Þ in the n-th stratum h n (izj~n), the expected number of fertilized females is given by Hence their mean (or total) number in host population The force of snail infection L~B : Y is proportional to Y, with the transmission rate/worm, B, dependent on multiple factors including egg production/release per female, intermediate larval survival in the transmission environment, and human host behavior.
To compute Y and L, one needs some assumptions on worm acquisition, the sex ratio distribution P i, jjn ð Þ, the mating pattern h(above), and other fecundity/ fitness parameters. May [12] distinguished two types of worm acquisition: (i) Togetherness whereby both sexes come through the same accumulation process with equal probability (50.5), where the sex ratio obeys a standard binomial: (ii) Separateness whereby each sex comes from its own (independent) accumulation process, hence For the present, field data on Schistosoma infection in rats [29] suggest that togetherness is the proper mating pattern in the wild. Using equation (4) leads to a closed form expression for the mated female count with m~n=2 ½ -''integer part'' of n=2 (see [12] and S3 Appendix in S1 Text for details).
Function y n f g can now enter our SWB formulation of L (as equation (3)). It was used with prescribed (NB, Poisson) distribution patterns, h n f g in earlier works [7,12] to derive a suitable mating function 0vw w ð Þv1. The latter measures the expected fraction (probability) of mated females per host, and the total mated (female) count expressed as: In special cases of distribution h n f g (uniform burden, Poisson, or NB) the mated female worm count Y, and mating function w can be computed in closed analytic form (see Table 2, Fig. 3, and, for derivation, S3 Appendix in S1 Text).
A parameter commonly used to define the efficacy of transmission control is the infectious prevalence, V, defined as the population fraction that carries at least one mated couple (different from the commonly used infection prevalence based on worm count). To estimate V we note that the probability of ''zero couples'' in the n-th strata is 1{2 1{n . Hence for SWB formulation, For MacDonald-MWB type systems with prescribed distribution (NB, Poisson), May [12] derived explicit formulae for V (Table 2). Therefore, in our revised estimation of community transmission, we may now include the function V in calibrating model parameters based on diagnostic egg count data from control programs.

Numeric implementation
Our analysis of MWB and SWB systems combines analytic tools with a substantial amount of numeric simulations. The latter applies to equilibria and parameter Table 2. Mating function for specific distribution patterns; 2 F 1 a, b, cjz ð Þis the hypergeometric function, I m z ð Þ -modified Bessel functions of order m.

Distribution Mating function
Prevalence of mated couples (host infectivity) V (abbreviation: NB, negative binomial, see S3 Appendix in S1 Text for details).
doi:10.1371/journal.pone.0115875.t002  space analysis on the one hand, and to dynamic simulations for prediction/control on the other. All numeric codes and procedures were implemented and run in Wolfram Mathematica 9, with differential equation solvers that offer eventcontrol tools adapted for simulation and analysis of control interventions. A basic version of our Mathematica program notebook (nb) for this analysis is posted online as S1 Workbook, to this paper.

The Macdonald MWB system
This simpler model of Schistosoma transmission for a single population has two variables: w t ð Þthe MWB of host population, and patent (or infectious) snail prevalence, y t ð Þ. When including a mating factor, these variables obey a coupled differential system: The mating function 0vw w ð Þƒ1 depends on underlying assumptions on parasite distribution (NB, Poisson, etc.) and, in many cases, it can be computed explicitly (see Table 2, Fig. 3, and S3 Appendix in S1 Text).
Transmission rates A and B lump together multiple biological, environmental, and human behavioral factors and reflect the success of intermediate larval stages. In particular, A is proportional to snail population density (N), A~a N (with per capita rate5a), while B is proportional to total human population (H), B~b H. Equation system (9) assumes stationary values for (H, N), but it can be easily extended to non-stationary cases (e.g., changing human demographics, or seasonal variability of snail population and transmission). One can also include the population turnover (rate m) in equations (9), by changing worm loss term c?czm in the w-equation. The basic model can be further extended to various heterogeneous settings (age-structured and/or spatially distributed communities). Such modified MWB systems have been used extensively in the prediction/control analysis (see, e.g., [20,22,33]).

The extended MacDonald MWB systems: moment equations and dynamic aggregation
A serious drawback of using the NB assumption in MacDonald-MWB systems is the uncertainty about (or evident variability) of the aggregation parameter k across time, age groups, and communities. This applies to different geographic environments and, more importantly, to the same system subjected to dynamic changes (e.g., by drug treatment) [29]. So, fixing k, as estimated from observed data, in the mating function w w, k ð Þ of equations (9) leads to inconsistency, as shown in drug treatment simulation studies [34][35][36]. To accommodate changing k values, some workers have proposed making k a dynamic variable, e.g., k~k 0 zk 1 w -a linear function of w, and then estimating coefficients k 0 , k 1 ð Þfrom field data [35,36]. In general, one could expect an increase of k with w (a higher average burden implies higher aggregation), but, in reality, the relationship may not be linear.
A more consistent way to introduce dynamic aggregation within the NBframework is via moment equations derived from the underlying SWB, namely 1 st moment (MWB)w t ð Þ~P ?
have human FOI, l~A y, decay rates (c -worm mortality, m -population turnover), and the external (demographic) sources S w , S u ð Þ derived from the underlying SWB source S n f g (see S2 Appendix in S1 Text). The moment system (10) for variables w t ð Þ, u t ð Þ f gcan be coupled to snail equation as in (9) via mating function w w, k ð Þ. Then the NB assumption on host strata gives an algebraic formula for aggregation k expressed through variables w, u ð Þ as So the snail equation in (9) turns into (see [37,38], and S2 Appendix in S1 Text).

The SWB system
The stratified worm burden (SWB) approach has been used before in theoretical studies (e.g., [33,37,38]), mostly to derive the reduced (MWB-type) ''moment'' equations, like (9) or (10)-(12), above. The basic dynamic variables of the SWBsystem are population strata h n t ð Þ : n~0,1,::: Fig. 1). They obey a differential Þh 0 zc 1 h 1 zS 0 ; ::::: :::::: Here l is per capita force of human infection (rate of worm accretion), mpopulation turnover (combined mortality, aging, migration), c k~k c -resolution rate for k-th strata (proportional to worm mortality, c), and S k f g -demographic source term. The latter represent infections brought into a given population group from outside. Thus the youngest age-group has only a newly born (uninfected) source S 0 w0 (proportional to the birth rate), while all other S n~0 . The older groups have their sources coming from younger groups, while in-and outmigration can also create additional sources and sinks for geographically coupled populations. The human part of the system (13) is coupled to the snail equation by two FOI factors: the human, l~A : y, and the snail, L~B : Y, i.e., proportional to the mated worm count given by function Y of equation (3).
The worm strata in the SWB setup are defined by a worm increment Dw~1,2,::: per stratum, so h n consists of hosts carrying n Dwƒwv nz1 ð ÞDw worms. In theoretical studies, fine-grain strata (Dw~1) are commonly used, but for practical modeling applications, larger increments are more appropriate (see [23,32]).
To extend the basic SWB setup [23,32] by including parasite mating in the force of snail infection, L, we can use May's estimate (equation (6)) of the expected number of mated couples y n (see [12] and S3 Appendix in S1 Text), but these estimates employ the optimal (combinatorial) male-female pairing count. Not all such couples are likely to be realized at low parasite densities, when the maturational (trophic) effects of male-female worm pairing may go missed [26]. To account for a low-density mating hurdle (the Allee effect [25,26]), we augment May's factors y n f g with an additional density-dependent success rates, 1{q n , where parameter 0vqv1 -measures the probability of mating failure. This factor would predict lower mating success in low-n strata, but approach a higher mating success rate of 1 at higher n. The resulting count of mated pairs takes the form: The force of snail infection is now expressed as a combination of variables h m t ð Þ f g with weights y m, q ð Þ, Similar to MacDonald-MWB system's equations (9), the coupled SWB-snail system consists of the human part (equation (13)) with FOI l~A y, and snail equation Both human and snail equations need some modification when worm increment Dww1. Here FOI of equation (15) is changed into while the human FOI of equation (13) is changed from l~A y (at Dw~1) to l~A y=Dw.
It is important to elaborate the parallels and differences between the two types of models, and how this influences their predictions for transmission control. Table 3 summarizes the system components and formulae for the MWB and SWB modeling approaches. The key inputs for both cases are: m -host turnover, cworm mortality, n -snail mortality. Some can be estimated from published studies (e.g., c~:25/year, n~4-5/year, (see Table 4)), while others involve known demographics in endemic areas (e.g., m~1=15 year for children, etc.). The most important parameters are transmission rates A, B, which must be estimated from observed human and snail infection data. Results

MWB equilibria and breakpoints
The simplest MacDonald-MWB system (equations (9)), without a mating component (i.e., w~1) has a stable-unstable pair of equilibria (infection-free + endemic state), provided that the Basic Reproductive Number (BRN, also known as R 0 ) For R 0 v1 it goes to elimination (stable zero equilibrium). The BRN R 0 is made of two factors that measure relative input of snail-to-human transmission (T SH~A =c), and human-to-snail transmission (T HS~B =n). The former, T SH , can be viewed as the maximal MWB-level (for a given transmission system) attained at the highest snail prevalence, y Ã~1 . As mentioned, to account for population turnover, formula (18) should be modified by changing c?czm. It is important to note that, in this model, any system with R 0 w1 cannot go to elimination, as any positive infection level (no matter how small) is predicted to eventually bring it back to the stable endemic state.
Addition of a mating function, 0vw w ð Þv1, vanishing at w50, has profound effect on equilibria and the dynamics of MWB equation (9) by turning it into a bistable system, [12,39]. Now the transition from ''stable zero'' to a bistable (endemic) regime requires higher transmission rates. Specifically, there exists a critical value, R c~Rc k,a ð Þw1, depending on aggregation, k, and snail-to-human transmission a~A=c, such that R 0 wR c is bistable (endemic), while R 0 vR c goes to elimination. Function R c k,a ð Þ has no simple algebraic form like equation (18), but we can explore it numerically.
The analysis exploits the reduced snail equation obtained from equilibrated worm burden w Ã~A y c of system (9), so roots of F give equilibrium values y Ã . Function, F y ð Þ, has a typical S-shaped pattern over 0ƒyƒ1, with either a single root y~0, or triple root 0vy b vy Ã (zero,  Panel (a) shows 3 critical (bifurcation) curves R c k, a ð Þ in the k, R 0 ð Þplane for a~:5; 1; 4. The region below each curve R 0 vR c k, a ð Þ has stable zero (elimination), whereas the region above R 0 wR c k, a ð Þ has a bistable (triple equilibrium) state. Panel (b) shows the same ''stable/unstable'' critical partition in the (a,b) parameter plane for three aggregation values :05vkv10. Panel (c) shows functions F y; R 0 , a, k ð Þof equation (19) for fixed R 0~6 :54 (5R c 1,1 ð Þ), and k~:4; 1; 1:5 (above, at, and below critical R c ) corresponding to the three marked points on panel (a). Panel (d) shows functions F y; R 0 , a, k ð Þat fixed k~1 and 3 values R 0 (above, at, below 6.54) marked on panel (a)). Panel (e) shows the phase plane of a bistable MWB system with three marked equilibria. Two ''separatrices'' at the breakpoint (in the middle) divide the phase plane into two attractor regions: ''zero'' -on the left and ''endemic'' -on the right. Panel (f) shows the bifurcation diagram of the MWB system for three types of equilibria, y Ã R 0 , a, k ð Þ , where the upper branch is stable ''endemic'', the low is stable ''zero'', the middle (gray) are breakpoints. The dashed curve corresponds to critical (bifurcation) values R c a, k ð Þ. Here k510, and the four curves arise from four different values of a in the range 1ƒaƒ6 (increased a pushes bifurcation curves to the left). Roots of F y ð Þ depend on three parameters: a~A c , R 0 , k , or a, b, k ð Þ. As these parameters change, the system undergoes a transition from a stable ''zero'' state (single equilibrium) to a triple equilibrium state, illustrated in bifurcation diagrams of Fig. 4 (f).
(a) shows k, R 0 ð Þparameter plane separated by the critical R c k, a ð Þ into ''infection-free'' range (below each curve) and ''endemic/bistable'' range (above it). In particular, values k~1, a~1, correspond to an R c~6 :54. Panel (b) of Fig. 4 shows a similar ''zero-endemic'' partition in the parameter space (a~A=c, b~B=n) for different levels of aggregation, k. As in panel (a), the respective infection-free ranges lie below the marked curves (bv 2 a R c k, a ð Þ), and the endemic ranges above them. For comparison we also show in (b) the ''infectionfree'' range ( a b 2 ƒ1) of a simple MacDonald-MWB system without mating (shaded). As shown, the effect of obligate mating (function w w, k ð Þ) is to raise the thresholds for endemicity, R c k, a ð Þw1, and higher k (clumping) gives higher R c values.
These results indicate that, under similar environmental conditions, sustained transmission in the ''mated'' MWB system is less likely that in the ''simple'' MWB, hence it would be easier to eradicate. Among different ''k -systems'' (NB vs. Poisson) higher clumping k makes the endemic state less tenable, and eradication potentially easier. Panels 4(c) and 4(d) illustrate profiles of function F y ð Þ above and below bifurcation value R c for several levels of aggregation (k), and BRN (R 0 ) values. Panel 4(e) shows a typical (w, y) phase portrait for a bistable MacDonald system with three equilibria, and schematic trajectories (arrows). The saddle-type breakpoint equilibrium (in the middle) has two ''separatrices'' (stable orbits) that divide the phase plane into two attractor regions: one solution driven to zero (elimination) on the left, and those relaxing to the endemic state on the right. Another effect of mating function w w, k ð Þ displayed in panels 4c, d, and f is a finite jump of endemic equilibria w Ã , y Ã ð Þas R 0 moves past the bifurcation value (R 0 wR c ), typical for many bistable systems. In contrast, a simple (w~1) MWB endemic equilibrium undergoes a gradual transition (y Ã~1 {1=R 0 ) with R 0 . This feature is related to hysteresis, whereby a gradual change of model transmission rates could, at a specific point, bring about a significant jump of endemic levels (an outbreak), rather than slow change.
In sum, sustained infection in any MacDonald-type system can be interrupted (brought to local extinction) by reducing transmission rates A, B (or R 0~A : B=n : c) below critical levels. But a typical intervention (MDA or snail control) does not affect the core transmission environment reflected by (A, B). Mathematically, a simple ''no-mating'' MacDonald system with BRN.1 cannot be brought to extinction, even after many control steps.
If valid, the breakpoint phenomena demonstrated for the extended MWB with mating in Fig. 4 could have significant implications for Schistosoma control [7,12,40]. The impact of control interventions are explored for MacDonald and SWB systems below. The extended MWB system (10)- (12) can be analyzed similarly to the basic MWB case. Equilibria w Ã , u Ã , k Ã ð Þof system (10) can be computed in terms of rescaled rates l, m ð Þ (relative to c), and demographic sources S w , S u ð Þcontributed by birth, aging, and migration (for details see S2 Appendix's formulae (32)- (33) in S1 Text). Depending on population type (e.g., a young cohort), S w , S u ð Þcould be zero or nonzero.
Of special note, our analysis revealed substantial differences in predictions between the zero and nonzero source conditions. In the ''zero'' source case, equations (32) from S2 Appendix in S1 Text, give equilibrium value k Ã~1 z2=m, independent of the transmission intensity l, and indicate a k greater than 1, whereas observed aggregation values are often,1. The ''nonzero'' (positive) source case give a more complicated expression for k Ã that can be studied numerically. The problem arises when function k Ã l, m, S w , S u ð Þturns negative (non-physical) in certain regions of the l, m ð Þparameter space. Such regions always exist, as long as S w , S u w0. While demographic parameter m is typically fixed, FOI l(proportional to snail infection y) could undergo big changes due to interventions (MDA) that could take it into an ''unphysical'' domain. The only way to maintain strictly positive k Ã throughout the l, m ð Þ plane is to have zero source terms (S w~Su~0 ). The unphysical l, m ð Þ-regions create problems for dynamic simulations in situations after MDA when l changes abruptly, if computed results fall into unreal/impossible ranges.

SWB equilibria and breakpoints
A complete SWB system consists of an infinite set of variables h n t ð Þ f g (n~0,1,2,:::), but for practical applications and computation, we truncate it at a finite (maximal) burden level N (0ƒnƒN). The choice of N has minor effect on the system's behavior and outputs, as long as FOI l (or MWB w<l) remain small relative to N. Another practical consideration concerns SWB-granularity, defined by using worm increments Dw §1. In some applications, it might be advantageous to reduce the number of variables (system dimensionality) by coarse-graining, from N (Dw~1) to N=Dw ð Þ. Overall, an increased step Dw would lower FOI, L, but when done consistently, its effect could be minimized (see Fig. 5).
Another SWB input -mating hurdle 0vqv1, related to the Allee phenomenon, has a more pronounced effect on model projections, particularly at low burden (small n). One way to assess it is through an estimated ''mated count per worm'' in the n-th strata, with mating factor y n, q ð Þ given by (14). As n increases, w n approaches its maximal value, 1/2, at a q-dependent relaxation rate. Fig. 59s left hand panels show the differential q effect, whereby an increasing hurdle factor from 0.1 to 0.95  (20) for fine-grain system Dw~1 (gray line) vs. coarse-grained Dw~5 (black dots).
Step Dw has only a minor effect on fraction w n (i.e., low sensitivity), but the force of infection L l, q, Dw ð Þ(column (b)) shows more sensitivity to Dw (higher Dw predicts stronger force, particularly at low l). Overall, the mating effect on L is significant -all curves in (b) depart from the simple linear relation L~B l (the thin, straight line), and increased hurdle factor q also has a significant effect in lowering L. would significantly slow the system's relaxation rate by 1/2. The immediate effect of the reduced mating capacity at low n is a significant reduction of FOI, L (right hand panels of Fig. 5). As explained below, this behavior is primarily responsible for the breakpoint phenomena in SWB systems.
Turning to equilibrium analysis of the SWB system (13), under prescribed FOI, population turnover, and population sources, we use rescaled values l, m, S ð Þover worm mortality, c. As mentioned, no analytic solutions for (13) exist except the limiting case m~0, S50 (stationary host population without turnover). Here h Ã n l,0 ð Þ È É becomes the principal (zero) eigenvector of the Frobenius-type matrix A of (13) which gives h Ã n l, 0 ð Þ~e {l l n n! -a Poisson distribution. We expect h Ã n l, m ð Þ È É for small m could be approximated by the limiting Poisson h Ã n l ð Þ. Fig. 2 shows numeric simulation of for the young-age group with a stationary uninfected source. Panels (a-b) compare NB/Poisson distributions with fixed mean w~10, to SWB-distributions h Ã n l, m ð Þ È É with fixed l~10, across a range of NB aggregation values k. They produce similar patterns, whereby clumping increases as k?? (for NB), and m?0 (for SWB). In that sense, the dimensionless SWB-turnover time 1=m plays the same role as NB-aggregation k. Note that realistic demographic values are relatively small m~:1{:2 (based on 20-40 year human life span [41], vs. worm 1=c54 years [33]). As expected from the m~0 case, they look like Poisson distribution results with mean m<l; Fig. 2 (c-d) demonstrate this for m~:1. We observe closely matched Poisson cases (c) and ''small m'' SWB cases (d), whereas large m cases (b) correspond to increased NB aggregation k. We conclude that in the absence of other confounding factors, a simple (homogeneous) SWB host system with an uninfected population source would attain a Poisson-like equilibrium state with w<l.
The last panel, 2(f), shows infection zero-prevalence V l, m ð Þ for several values, m. Once again, we note parallels with the corresponding MacDonald NB prevalence functions for increased k (panel 2(e)).
Turning to fully coupled SWB + snail systems, equilibria can be computed from the reduced snail equation, described by function F SW y ð Þ F SW y ð Þ~B n : X n §1 y n Dw, q Dw which plays a similar role to Macdonald function F y,::: ð Þ, (equation 19).
In the Macdonald case, F y ð Þ undergoes a transition from ''stable zero'' to a ''bistable'' (breakpoint) state, depending on model parameters, and zero is always a stable equilibrium.
Our analysis of SWB (S4 Appendix in S1 Text) reveals a more complicated picture. Specifically, we identified three cases: (i) single stable zero (eradication); (ii) double stable/unstable pair (''zero + endemic''), as in the Macdonald MWB system without mating (w~1); and (iii) a bistable (zero-breakpoint-endemic) case, like the mated extended-MWB case. The outcome depends on transmission rates A, B ð Þ and the mating hurdle 0ƒqv1. Unlike the Macdonald case (see Fig. 4(b)), the parameter space is now divided into 3 regions. A condition for a stable zero equilibrium (y Ã~0 ) is negative slope dF SW dy 0 ð Þv0. The slope can by computed in terms of parameters A,B, q, and the worm-step Dw used in SWB formulation (see S4 Appendix in S1 Text), namely. d dy where R 0 is the standard Macdonald BRN (18) adjusted for population turnover (c?czm) and y n -mated fraction (2). In the simulations described below, we used step Dw~2. Condition (22) is similar to stability of the ''zero'' equilibrium (R 0 v1) for a simple Macdonald system (w~1). In that sense, the SWB system occupies an intermediate place between two Macdonald types: the simple ''no mating'' (w~1), and the ''breakpoint'' (wv1) containing system. A more challenging task was to identify stable endemic regions in the A, B, q ð Þparameter space of SWB. Unlike Macdonald F y ð Þ, the SWB F SW y ð Þ has no simple algebraic form, so we computed it numerically (see S4 Appendix in S1 Text). The results are shown in Fig. 6: the shaded region in the A,B plane marks the bistable (breakpoint) parameter values, above this shaded area, the system is ''stable endemic'', below the area, it goes to ''stable zero'' (eradication). We observe that having an increased q would expand the breakpoint region, and shift it up in the (A, B) plane. Qualitatively, the stable endemic regions of the Macdonald system in Fig. 4(b) and those of the SWB (Fig. 6) look similar. Fig. 7. demonstrates F SW y ð Þ for fixed rates A~6, B~5:8 and several q (above, in, and below the breakpoint region) to observe all three equilibrium patterns. Their profiles resemble MacDonald-MWB functions F y ð Þ of Fig. 4 (panels 4c and 4d) within the breakpoint parameter ranges (dashed and light gray curves), but they look qualitatively similar to the simple (no-mating) MacDonald function F y ð Þ (black curve for q~:95) above the shaded region in Fig. 6(c). Panel 7(b) shows the resulting endemic and breakpoint equilibrium distributions h Ã n È É in the breakpoint case (a) q~:95, the former being more aggregated with higher MWB value.
Overall, a significant finding of this SWB modeling is that the breakpoint regions occupy a relatively small fraction of (A, B, q) space. So, randomly chosen A, B would most likely result in a simple Macdonald-type dichotomy: either ''stable zero'' or ''stable endemic''. In the real world, however, situations A, B might be related, and the breakpoint could actually play a more significant role in determining Schistosoma persistence. This prompted us to explore their dynamic implications.

Dynamic responses and projected effects of drug treatment
To reflect the impact of control programs in endemic Schistosoma transmission settings, it is important to compare the respective dynamic responses of the MWB and SWB models, including their endemic equilibria, relaxation patterns, and the long term impact of control interventions. While there are some major differences in their structure, the two models have some similarities in terms of a comparable set of parameters: A, B, aggregation k, for MacDonald; and A, B, and the mating hurdle, q, for SWB. One way to compare the two types of model is to calibrate them with an identical data set and conduct numeric simulations of control. (The calibration procedures are outlined in S5 Appendix in S1 Text).
Drug treatment with praziquantel clears a sizable fraction of adult worms (up to 90-95%), and thereby reduces worm burden in treated populations. There are two essential parameters of mass drug administration (MDA): i) drug efficacy human population fraction covered by treatment 0vf v1. Other important factors in program efficacy are the frequency (timing) and the number of MDA sessions. Mathematically, MDA is implemented differently in the two different types of model.
Let us note that our data set was chosen in a special way, to produce a breakpoint-type SWB system. Fig. 8 shows two reduced equilibrium functions: MacDonald's F y ð Þ of equation (19) (shown in gray), and the SWB-function F SW y ð Þ of equation (21) (in black). Both exhibit breakpoints near y50, but SWB has higher value y Ã B than MacDonald's (see Table 5). It suggests that SWB infection would be easier (faster) brought to elimination compared to MacDonald case.
Here worm mortality for the treated group undergoes an abrupt change at the treatment time t 0 d t{t 0 ð Þ. The combined force of snail infection by the two groups is given by  Table 4 doi:10.1371/journal.pone.0115875.g008 Table 5. Calibrated model parameters from the data inputs listed in Table 4.
, represented by Dirac delta-function Another way to implement it numerically is to terminate solution (23) at t~t 0 , and reinitialize using MWB variable w t ð Þ~1{f ð Þ : u t ð Þzf : v t ð Þ, as In practice, this means that population is randomly divided into treated/ untreated fractions in each session and no prior treatment data are incorporated. Such schemes can be repeated any number of times (t 0 vt 1 v:::) with prescribed frequency, and prescribed treatment fractions (f 0 ; f 1 ; :::). Furthermore, these schemes could be augmented with additional features of monitoring and control, e.g., control termination after infection levels are brought below a specified level (the natural choice would be a 'breakpoint'). Fig. 9 compares a hypothetical MDA control program, having 70% coverage and a 90% cure rate, for two calibrated MacDonald systems: simple w w ð Þ~1  (36) in S3 Appendix in S1 Text). Both systems are initialized at their endemic states. The former scenario indicates that the region would require an indefinite series of treatments to maintain control, while the latter suggests that MDA would bring the system to eradication after 4 sessions, when wv:13, drops below the breakpoint value.

SWB drug control
The effect of drug treatment on the SWB model is to move a treated fraction from higher strata h n t ð Þ to lower-level strata h m t m : 0ƒmv1= f would shift to h 0 (complete clearing), the next higher range h m : 1= f The corresponding MDA reinitialization event then becomes h m t 0 ð Þ h 1 j t 0~: :: However, the snail equation (force of infection L) doesn't change its functional form as variables h n t ð Þ f g get reshuffled. To examine the effect of control on long term SWB histories we took the calibrated SWB system with parameters of Table 4 and ran the same four-year control strategy as for MWB described above. The results of simulation are compared to the calibrated MacDonald-MWB system in Fig. 10(a). Overall, SWB simulations predict more efficient reduction of worm burden and prevalence in the four-session control program. Both systems predict elimination below their breakpoints.
The concerning feature of SWB prediction is its high sensitivity to the Allee parameter q. A slight change from q~:9 (breakpoint case I) to q~:87 (saddlenode case II), has important long-term implications, shown in Fig. 10(b). In both case I and case II, MWB can be brought to relatively low levels after 4 sessions, but case I goes to extinction while case II relaxes back to pretreatment endemic levels -i.e., a finite eradication time for case I vs. the requirement for indefinitely sustained effort for case II.

Conclusions
Uneven parasite burden and sex distribution have a significant impact on infection levels and sustainability of transmission. Indeed, both depend on ð Þ , with m< n determined by the drug   efficacy [23]. In particular, all strata having h g  ƒmv2= gwould go to h 1 etc.      fertilized female count, and uneven parasite loads create hurdles for worm mating that reduce resulting egg production, particularly in low-level worm burden host strata. May and colleagues [12,18] have addressed these issues in the context of Macdonald MWB formulation, utilizing an assumption about infection distribution patterns (e.g., NB, Poisson) to facilitate solutions to their model. This analysis produced a modified Macdonald-MWB system that includes a ''mating factor''. The mating factor makes profound changes in mathematical structure and equilibria of Macdonald-type system, creating a breakpoint between its ''zero'' and endemic levels. Hence, a modified Macdonald system (with mating), Fig. 10. Control of the calibrated SWB system with the same MDA strategy as for previous figure. The left column (a) compares SWB outputs (mean worm load, human and snail prevalence) in the solid curves, compared with the MWB model outputs shown in the dashed curves. SWB predicts faster (more efficient) reduction of all infection outcomes. Panel (b) compares the SWB system of column (a) (solid) with an MDAperturbed system where mating hurdle was changed to q~:87 (dashed). The perturbed system has no breakpoint (which has sensitive dependence on q in the SWB system) and after an initial four-treatment reduction, infection gradually relaxes to the pre-MDA endemic state doi:10.1371/journal.pone.0115875.g010 unlike its simple cousin, predicts elimination after a finite number of control interventions.
Consideration of the MWB work raises several questions: how reasonable is the NB assumption?; should NB aggregation parameters be fixed or subject to change?; how can these models be reconciled with underlying host demographics?; and how reliable are assumptions about the effects of interventions (MDA) on model dynamics and parameters? In place of the MWB approach, we believe that a proper way to account for mating and uneven distribution is through a stratified worm burden approach (SWB) to model development. Mathematically, worm strata can be viewed as population level distribution function (PDF) of the underlying stochastic process of worm acquisition/loss. Some important conclusions of the MWB-SWB comparative analysis: i) In the modified MacDonald-MWB systems, the NB aggregation parameter k should not be treated as fixed, but should be treated as a dynamic variable. When such systems are calibrated based on equilibrium relations, the outcomes are highly sensitive to k. Furthermore, k-values can undergo significant changes after intervention. ii) A consistent dynamic formulation (based on the underlying SWB) requires an extension of the Macdonald system with either k as an additional variable, or an equivalent ''second-moment'' variable (where the MWB w represents its ''first moment''). iii) However, this extended (1st+2nd moment) system has inherent inconsistencies when coupled to host demographics. The consistent formulation is possible only within a limited context -when the sole source of host population enters the uninfected (''zero-level'') SWB strata, e.g., as newborn children. So one can provide an extended Macdonald system for a children's age category (with a newborn source), but not for a corresponding adult group or for any coupled ''child-adult'' systems. By contrast, the SWB system is free of such limitations and can be set for any age-or locationstructured populations. iv) In exploring the breakpoint phenomena for both MWB (fixed k) and SWB systems, in the former case (MWB) a breakpoint comes automatically in any setting with a positive endemic level, i.e., with sufficiently high transmission rates (a, b); in the latter case (SWB), zero and endemic equilibria are not rigidly connected, so there are three possible outcomes: ''zero-breakpointendemic'', ''zero-endemic'', and ''zero''. The existence of breakpoint depends on an additional (low-density) mating constraint. Here we have accounted for such an effect by a single parameter q -the probability of mating failure per single adult. We showed that in the SWB system, the breakpoint phenomena depends strongly on q and we have explored the bounds of breakpoint regions in the (a, b, q)-parameter space. v) In both cases, (MacDonald and SWB), we can explore the effect of drug treatment (MDA) with different parameter values (coverage fraction, drug efficacy, treatment frequency/year). In all cases we tested, a breakpoint was shown to bring elimination after a finite number of interventions whenever infection levels fell below the breakpoint value.