Correction
1 Apr 2015: The PLOS ONE Staff (2015) Correction: Population Biology of Schistosoma Mating, Aggregation, and Transmission Breakpoints: More Reliable Model Analysis for the End-Game in Communities at Risk. PLOS ONE 10(4): e0123592. https://doi.org/10.1371/journal.pone.0123592 View correction
Figures
Abstract
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.
Citation: Gurarie D, King CH (2014) Population Biology of Schistosoma Mating, Aggregation, and Transmission Breakpoints: More Reliable Model Analysis for the End-Game in Communities at Risk. PLoS ONE 9(12): e115875. https://doi.org/10.1371/journal.pone.0115875
Editor: Ulrike G. Munderloh, University of Minnesota, United States of America
Received: August 13, 2014; Accepted: December 2, 2014; Published: December 30, 2014
Copyright: © 2014 Gurarie, King. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. All relevant data are within the paper.
Funding: D. Gurarie was supported by the sabbatical fellowship at NIMBioS (University of Tennessee) during Academic Year 2012–13. This work is also supported by the Schistosomiasis Consortium for Operational Research and Evaluation (SCORE) funded by the University of Georgia Research Foundation through a grant from the Bill and Melinda Gates Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
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]–[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]–[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 MacDonald-type 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 distribution-free 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 within-host 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 age-related demographic factors [19], [20].
Methods
Modeling worm aggregation in MWB and SWB
Different approaches have been used to describe Schistosoma worm distributions, ranging from a nearly uniform host burden (whereby each host carries approximately worms- the total parasite load W distributed evenly over host population – H), to types of over-dispersed distribution such as the negative binomial (NB) distribution or the Poisson distribution. The former, NB, is defined by two parameters: aggregation, k, and mean w, where probability
and
(1)
Increased k gives a more aggregated (clumped) distribution for ; in the limiting case where
, 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]–[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]–[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 to a low-dimensional (MacDonald-type) “moment system” for the MWB variable
, 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
. 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: , human FOI, balanced by worm mortality,
, and demographic loss,
. Its equilibrium distribution
depends on rescaled values
and
, the former
having dimension of [worm burden] like the MacDonald-MWB term w, while
is dimensionless. The resulting SWB solutions depend on demographic source terms (equations (13) below).
The SWB model includes population strata , sources
, population turnover/loss rates
, the force of infection
( = worm accretion rate), and worm clearing rates
(
is worm mortality).
Here, all SWB distributions (right-hand panels) are produced from an uninfected source . (a) NB-MWB with fixed mean
and increasing k (
is the limiting Poisson case); (b) Equilibrium SWB distribution for uninfected source with FOI (mean)
, and varying demographic parameter
(see equation (13));
plays the role of aggregation k for NB, with small
corresponding to large (infinite) k. Panels (c−d) Poisson distributions with different means,
(left panel) vs. SWB distributions with different
(right panel) exhibit striking similarity. Panels (e)−(f) compare infectious prevalence
for the two models as a function of NB-MWB w or SWB
.
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, , only (
,
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,
,
(no population turnover and sources). Here equilibrium solution gives the standard Poisson distribution
, 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 = the number of mated (egg-shedding) females for i - males, j - females. 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, . Specifically, for a given sex distribution
in the n-th stratum
(
), the expected number of fertilized females is given by
Hence their mean (or total) number in host population
The force of snail infection is proportional to
, 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 and
, one needs some assumptions on worm acquisition, the sex ratio distribution
, the mating pattern
(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 ( = 0.5), where the sex ratio obeys a standard binomial:(4)
(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
(6)with
- “integer part” of
(see [12] and S3 Appendix in S1 Text for details).
Function can now enter our SWB formulation of
(as equation (3)). It was used with prescribed (NB, Poisson) distribution patterns,
in earlier works [7], [12] to derive a suitable mating function
. 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 (uniform burden, Poisson, or NB) the mated female worm count
, and mating function
can be computed in closed analytic form (see Table 2, Fig. 3, and, for derivation, S3 Appendix in S1 Text).
Shown are results for (Poisson). Note that a higher degree of clumping lowers
and reduces transmission potential of the system.
A parameter commonly used to define the efficacy of transmission control is the infectious prevalence, , 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
we note that the probability of “zero couples” in the n-th strata is
. Hence for SWB formulation,
For MacDonald-MWB type systems with prescribed distribution (NB, Poisson), May [12] derived explicit formulae for (Table 2). Therefore, in our revised estimation of community transmission, we may now include the function
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 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 event-control 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.
Transmission Models that Are Compared
The Macdonald MWB system
This simpler model of Schistosoma transmission for a single population has two variables: - the MWB of host population, and patent (or infectious) snail prevalence,
. When including a mating factor, these variables obey a coupled differential system:
(9)
The mating function 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), (with per capita rate = a), while B is proportional to total human population (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
) in equations (9), by changing worm loss term
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 of equations (9) leads to inconsistency, as shown in drug treatment simulation studies [34]–[36]. To accommodate changing k values, some workers have proposed making k a dynamic variable, e.g.,
- a linear function of w, and then estimating coefficients
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 NB-framework is via moment equations derived from the underlying SWB, namely 1st moment (MWB) - ; 2nd moment -
, etc. The resulting two-moment systems
(10)have human FOI,
, decay rates (
- worm mortality,
- population turnover), and the external (demographic) sources
derived from the underlying SWB source
(see S2 Appendix in S1 Text). The moment system (10) for variables
can be coupled to snail equation as in (9) via mating function
. Then the NB assumption on host strata gives an algebraic formula for aggregation k expressed through variables
as
(11)So the snail equation in (9) turns into
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 SWB-system are population strata with total population
(schematic diagram shown in Fig. 1). They obey a differential equation system ([37]–[39])
(13)
Here is per capita force of human infection (rate of worm accretion),
- population turnover (combined mortality, aging, migration),
- resolution rate for k-th strata (proportional to worm mortality,
), and
- 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
(proportional to the birth rate), while all other
. The older groups have their sources coming from younger groups, while in- and out-migration 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,
, and the snail,
, i.e., proportional to the mated worm count given by function
of equation (3).
The worm strata in the SWB setup are defined by a worm increment per stratum, so
consists of hosts carrying
worms. In theoretical studies, fine-grain strata (
) 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, , we can use May’s estimate (equation (6)) of the expected number of mated couples
(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
with an additional density-dependent success rates,
, where parameter
- 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 with weights
,
Similar to MacDonald-MWB system’s equations (9), the coupled SWB-snail system consists of the human part (equation (13)) with FOI , and snail equation
(16)
Both human and snail equations need some modification when worm increment . Here FOI of equation (15) is changed into
while the human FOI of equation (13) is changed from (at
) to
.
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: - host turnover,
- worm mortality,
- snail mortality. Some can be estimated from published studies (e.g.,
/year,
–5/year, (see Table 4)), while others involve known demographics in endemic areas (e.g.,
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., ) has a stable-unstable pair of equilibria (infection-free + endemic state), provided that the Basic Reproductive Number (BRN, also known as R0)
(18)
For it goes to elimination (stable zero equilibrium). The BRN
is made of two factors that measure relative input of snail-to-human transmission (
), and human-to-snail transmission (
). The former,
, can be viewed as the maximal MWB-level (for a given transmission system) attained at the highest snail prevalence,
. As mentioned, to account for population turnover, formula (18) should be modified by changing
. It is important to note that, in this model, any system with
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, , vanishing at w = 0, 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,
, depending on aggregation, k, and snail-to-human transmission
, such that
is bistable (endemic), while
goes to elimination. Function
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 of system (9), so roots of F give equilibrium values
. Function,
, has a typical S-shaped pattern over
, with either a single root
, or triple root
(zero, breakpoint, endemic), provided
is sufficiently large,
- critical (bifurcation) value. These features are demonstrated in Fig. 4.
Panel (a) shows 3 critical (bifurcation) curves in the
plane for
. The region below each curve
has stable zero (elimination), whereas the region above
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
. Panel (c) shows functions
of equation (19) for fixed
( =
), and
(above, at, and below critical
) corresponding to the three marked points on panel (a). Panel (d) shows functions
at fixed
and 3 values
(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,
, 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
. Here k = 10, and the four curves arise from four different values of a in the range
(increased a pushes bifurcation curves to the left).
Roots of depend on three parameters:
, or
. 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 - parameter plane separated by the critical
into “infection-free” range (below each curve) and “endemic/bistable” range (above it). In particular, values
,
, correspond to an
. Panel (b) of Fig. 4 shows a similar “zero-endemic” partition in the parameter space (
,
) for different levels of aggregation, k. As in panel (a), the respective infection-free ranges lie below the marked curves (
), and the endemic ranges above them. For comparison we also show in (b) the “infection-free” range (
) of a simple MacDonald-MWB system without mating (shaded). As shown, the effect of obligate mating (function
) is to raise the thresholds for endemicity,
, and higher k (clumping) gives higher
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 above and below bifurcation value
for several levels of aggregation (k), and BRN (
) 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
displayed in panels 4c, d, and f is a finite jump of endemic equilibria
as R0 moves past the bifurcation value (
), typical for many bistable systems. In contrast, a simple (
) MWB endemic equilibrium undergoes a gradual transition (
) with
. 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 (or
) 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 of system (10) can be computed in terms of rescaled rates
(relative to
), and demographic sources
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),
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 , independent of the transmission intensity
, 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
that can be studied numerically. The problem arises when function
turns negative (non-physical) in certain regions of the
- parameter space. Such regions always exist, as long as
. While demographic parameter
is typically fixed, FOI
(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
throughout the
plane is to have zero source terms (
). The unphysical
-regions create problems for dynamic simulations in situations after MDA when
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 (
), but for practical applications and computation, we truncate it at a finite (maximal) burden level N (
). The choice of N has minor effect on the system’s behavior and outputs, as long as FOI
(or MWB
) remain small relative to N. Another practical consideration concerns SWB-granularity, defined by using worm increments
. In some applications, it might be advantageous to reduce the number of variables (system dimensionality) by coarse-graining, from N (
) to
. Overall, an increased step
would lower FOI,
, but when done consistently, its effect could be minimized (see Fig. 5).
Left-hand panels show the distribution of mated pairs in host strata, and right-hand panels show the force of snail infection . The three plots in each column correspond to different choices of q. Column (a) compares mated fraction
of equation (20) for fine-grain system
(gray line) vs. coarse-grained
(black dots). Step
has only a minor effect on fraction
(i.e., low sensitivity), but the force of infection
(column (b)) shows more sensitivity to
(higher
predicts stronger force, particularly at low
). Overall, the mating effect on
is significant - all curves in (b) depart from the simple linear relation
(the thin, straight line), and increased hurdle factor q also has a significant effect in lowering Λ.
Another SWB input –mating hurdle , 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,
(20)
with mating factor given by (14). As n increases,
approaches its maximal value, 1/2, at a q-dependent relaxation rate. Fig. 5′s left hand panels show the differential q effect, whereby an increasing hurdle factor from 0.1 to 0.95 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,
(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 over worm mortality,
. As mentioned, no analytic solutions for (13) exist except the limiting case
, S = 0 (stationary host population without turnover). Here
becomes the principal (zero) eigenvector of the Frobenius-type matrix A of (13)
which gives
- a Poisson distribution. We expect
for small
could be approximated by the limiting Poisson
.
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 , to SWB-distributions
with fixed
, across a range of NB aggregation values k. They produce similar patterns, whereby clumping increases as
(for NB), and
(for SWB). In that sense, the dimensionless SWB-turnover time
plays the same role as NB-aggregation k. Note that realistic demographic values are relatively small
(based on 20–40 year human life span [41], vs. worm
= 4 years [33]). As expected from the
case, they look like Poisson distribution results with mean
; Fig. 2 (c–d) demonstrate this for
. We observe closely matched Poisson cases (c) and “small
” SWB cases (d), whereas large
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
.
The last panel, 2(f), shows infection zero-prevalence for several values,
. 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
(21)which plays a similar role to Macdonald function
, (equation 19).
In the Macdonald case, 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 (); and (iii) a bistable (zero-breakpoint-endemic) case, like the mated extended-MWB case. The outcome depends on transmission rates
and the mating hurdle
. Unlike the Macdonald case (see Fig. 4(b)), the parameter space is now divided into 3 regions. A condition for a stable zero equilibrium (
) is negative slope
. The slope can by computed in terms of parameters A,B, q, and the worm-step
used in SWB formulation (see S4 Appendix in S1 Text), namely.
(22)where
is the standard Macdonald BRN (18) adjusted for population turnover (
) and
- mated fraction (2). In the simulations described below, we used step
. Condition (22) is similar to stability of the “zero” equilibrium (
) for a simple Macdonald system (
). In that sense, the SWB system occupies an intermediate place between two Macdonald types: the simple “no mating” (
), and the “breakpoint” (
) containing system.
A more challenging task was to identify stable endemic regions in the -parameter space of SWB. Unlike Macdonald
, the SWB
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.
Shown are values for () in A, B - parameter space at different values of mating hurdle q. The shaded region in each plot marks the bistable (zero-breakpoint-endemic) range; all values above it are saddle-nodes (“unstable zero + stable endemic”); below it is the region of “stable zero” (eradication).
Fig. 7. demonstrates for fixed rates
,
and several q (above, in, and below the breakpoint region) to observe all three equilibrium patterns. Their profiles resemble MacDonald-MWB functions
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
(black curve for
) above the shaded region in Fig. 6(c). Panel 7(b) shows the resulting endemic and breakpoint equilibrium distributions
in the breakpoint case (a)
, the former being more aggregated with higher MWB value.
Panel (a) shows equilibrium snail function, , with transmission rates
,
, exhibiting three equilibrium patterns: saddle-node (zero - endemic), bistable (breakpoint) and stable “zero”, for increasing mating hurdle q. Panel (b) shows SWB equilibrium distributions for the stable (endemic) equilibrium and unstable breakpoint for the middle curve
of panel (a).
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 - in terms of the fraction of surviving worms (
), and ii) the human population fraction covered by treatment
. 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 of equation (19) (shown in gray), and the SWB-function
of equation (21) (in black). Both exhibit breakpoints near y = 0, but SWB has higher value
than MacDonald’s (see Table 5). It suggests that SWB infection would be easier (faster) brought to elimination compared to MacDonald case.
Shown are : MacDonald equation (19) (gray), and SWB equation (21) (black), for calibrated model parameters of Table 4
MDA for Macdonald-type systems
For the Macdonald-MWB system, one divides population into treated+ untreated groups, with burden - treated,
- untreated) and sets up an extended version of system (9) for variables
Here worm mortality for the treated group undergoes an abrupt change at the treatment time ,
, represented by Dirac delta-function
. The combined force of snail infection by the two groups is given by
Another way to implement it numerically is to terminate solution (23) at , and reinitialize using MWB variable
, 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 () with prescribed frequency, and prescribed treatment fractions (
). 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 (panel (a)) vs. a modified NB-MacDonald with mating function
(panel (b), see equation (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
, drops below the breakpoint value.
Panels indicate results for (a) a simple Macdonald MWB model without mating; (b) a MWB model with NB distribution and an included mating function (equation (36), S3 Appendix in S1 Text) having breakpoint level (shaded range, lower left hand panel). On both left hand plots, solid black is the overall community MWB,
; gray is the untreated group
, and the dashed line is the treated group
. The upper right hand panel indicates yearly snail prevalence of infection without (solid line) or with (dashed line) inclusion of the mating factor, and hence the breakpoint, in the model. The lower right hand panel indicates expected human prevalence with MDA treatment in the breakpoint setting.
SWB drug control
The effect of drug treatment on the SWB model is to move a treated fraction from higher strata to lower-level strata
, with
determined by the drug efficacy
[23]. In particular, all strata having
would shift to
(complete clearing), the next higher range
would go to
etc. The corresponding MDA reinitialization event then becomes
However, the snail equation (force of infection ) doesn’t change its functional form as variables
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 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 MDA-perturbed system where mating hurdle was changed to (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
The concerning feature of SWB prediction is its high sensitivity to the Allee parameter q. A slight change from (breakpoint case I) to
(saddle-node 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 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), 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:
- 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.
- 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”).
- 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 location-structured populations.
- 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-breakpoint-endemic”, “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.
- 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.
Our study of mating, breakpoints, and infection persistence raises several issues on the role and meaning of BRN in transmission models. BRN by itself can predict a transition from “no infection” to “endemic state” (above/below critical value), but has no direct links to “breakpoints”. Mathematically the “breakpoint phenomena” amounts to a two-parameter space analysis. It has clear implications for parasite elimination: MDA, by itself, won’t necessarily affect BRN (transmission environment), but by driving the system below an inherent breakpoint, we might achieve the requisite transition (bifurcation) to zero endemicity. Future work with extended control program datasets and improved SWB methodology will shed the light on the probable existence of breakpoints for elimination. In extending our published work [23], we plan to apply the newly extended SWB methodology to model structured host populations and distributed human -snail environments [19], [20]. We believe that the modified SWB approach will provide more accurate and reliable prediction than the conventional MacDonald MWB-based methods.
Supporting Information
S1 Text.
S1–S5 Appendices. Extended description and formulae detailing: S1 Appendix) The background on the problems of Schistosoma transmission model formulation; S2 Appendix) Moment equations for SWB, and the extended MWB model; S3 Appendix) Mating function in Macdonald and SWB systems; S4 Appendix) Equilibria of a coupled human SWB - snail system; and S5 Appendix) Model calibration.
https://doi.org/10.1371/journal.pone.0115875.s001
(DOCX)
S1 Workbook.
Mathematica software notebook file (.nb) with programming for MWB and SWB transmission models.
https://doi.org/10.1371/journal.pone.0115875.s002
(NB)
Author Contributions
Conceived and designed the experiments: DG CHK. Performed the experiments: DG CHK. Analyzed the data: DG CHK. Contributed reagents/materials/analysis tools: DG CHK. Contributed to the writing of the manuscript: DG CHK.
References
- 1. King CH, Sturrock RF, Kariuki HC, Hamburger J (2006) Transmission control for schistosomiasis - why it matters now. Trends Parasitol 22:575–582.
- 2. Colley DG, Bustinduy AL, Secor WE, King CH (2014) Human schistosomiasis. Lancet 383:2253–2264.
- 3.
Bustinduy AL, King CH (2009) Parasitic Helminths. In: Fratamico PM, Smith JL, Brogden KAeditors. Post-Infectious Sequelae and Long-Term Consequences of Infectious Diseases. Washington: American Society for Microbiology Press. 291–329.
- 4. Warren KS (1982) Selective primary health care: strategies for control of disease in the developing world. I. Schistosomiasis. Rev Infect Dis 4:715–726.
- 5. Jordan P, Christie JD, Unrau GO (1980) Schistosomiasis transmission with particular reference to possible ecological and biological methods of control. A review. Acta Trop 37:95–135.
- 6.
Savioli L, Daumiere D (2012) Accelerating Work to Overcome the Global Impact of Neglected Tropical Diseases: A Roadmap for Implementation. Geneva: World Health Organization.
- 7. MacDonald G (1965) The dynamics of helminth infections, with special reference to schistosomes. Trans R Soc Trop Med Hyg 59:489–506.
- 8. Coura-Filho P, Mendes NM, de Souza CP, Pereira JP (1992) The prolonged use of niclosamide as a molluscicide for the control of Schistosoma mansoni. Rev Inst Med Trop Sao Paulo 34:427–431.
- 9. Duke BO, Moore PJ (1976) The use of a molluscicide, in conjunction with chemotherapy, to control Schistosoma haematobium at the Barombi Lake foci in Cameroon. I. The attack on the snail hosts, using N-tritylmorpholine, and the effect on transmission from snail to man. Tropenmed Parasitol 27:297–313.
- 10. Fenwick A (1970) The development of snail control methods on an irrigated sugar-cane estate in northern Tanzania. Bull World Health Organ 42:589–596.
- 11. Nåsell I (1978) Mating models for schistosomes. Journal of Mathematical Biology 6:21–35.
- 12. May RM (1977) Togetherness among schistosomes: its effects on the dynamics of the infection. Mathematical Biosciences 35:301–343.
- 13. Coura Filho P, Rocha RS, de Lima ECMF, Katz N (1992) A municipal level approach to the management of schistosomiasis control in Peri-Peri, MG, Brazil. Rev Inst Med Trop Sao Paulo 34:543–548.
- 14. Fenwick A, Jorgensen TA (1972) 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 47:579–586.
- 15. Barnish G, Jordan P, Bartholomew RK, Grist E (1982) Routine focal mollusciciding after chemotherapy to control Schistosoma mansoni in Cul de Sac valley, Saint Lucia. Trans R Soc Trop Med Hyg 76:602–609.
- 16.
Clarke VV, Blair DM, Weber MC (1973) The Mayfield experiment: an attempt to eliminate bilharziasis from a small community by repeated treatment of the infected people combined with intensive snail control. Cent Afr J Med: 14–22.
- 17. Nåsell I, Hirsch WM (1973) The transmission dynamics of schistosomiasis. Communications on Pure and Applied Mathematics 26:395–453.
- 18. Bradley DJ, May RM (1978) Consequences of helminth aggregation for the dynamics of schistosomiasis. Trans R Soc Trop Med Hyg 72:262–273.
- 19. Woolhouse ME, Chandiwana SK (1989) Spatial and temporal heterogeneity in the population dynamics of Bulinus globosus and Biomphalaria pfeifferi and in the epidemiology of their infection with schistosomes. Parasitology 98:21–34.
- 20. Woolhouse ME, Etard JF, Dietz K, Ndhlovu PD, Chandiwana SK (1998) Heterogeneities in schistosome transmission dynamics and control. Parasitology 117:475–482.
- 21. May RM, Woolhouse ME (1993) Biased sex ratios and parasite mating probabilities. Parasitology 107:287–295.
- 22. Gurarie D, King CH (2005) Heterogeneous model of schistosomiasis transmission and long-term control: the combined influence of spatial variation and age-dependent factors on optimal allocation of drug therapy. Parasitology 130:49–65.
- 23. Gurarie D, King CH, Wang X (2010) A new approach to modelling schistosomiasis transmission based on stratified worm burden. Parasitology 137:1951–1965.
- 24. Barbour AD (1996) Modeling the transmission of schistosomiasis: an introductory view. Am J Trop Med Hyg 55:135–143.
- 25. Thomas JD (1973) Schistosomiasis and control of molluscan hosts of human schistosomes with particular reference to possible self-regulatory mechanisms. Adv Parasitol 11:307–394.
- 26. Ribeiro-Paes JT, Rodrigues V (1997) Sex determination and female reproductive development in the genus Schistosoma: a review. Rev Inst Med Trop Sao Paulo 39:337–344.
- 27. Dobson A, Roberts M (1994) The population dynamics of parasitic helminth communities. Parasitology 109 Suppl: S97–108
- 28. Shaw DJ, Grenfell BT, Dobson AP (1998) Patterns of macroparasite aggregation in wildlife host populations. Parasitology 117 (Pt 6):597–610.
- 29. Morand S, Pointier J-P, Borel G, Théron A (1993) Pairing probability of schistosomes related to their distribution among the host population. Ecology 74:2444–2449.
- 30. Anderson RM, May RM (1978) Regulation and stability of host-parasite population interactions. I. Journal of Animal Ecology 47:219–247.
- 31. Anderson RM, May RM (1978) Regulation and stability of host-parasite population interactions. II. Journal of Animal Ecology 47:249–267.
- 32. Wang X, Gurarie D, Mungai PL, Muchiri EM, Kitron U, et al. (2012) Projecting the long-term impact of school-or community-based mass-treatment interventions for control of Schistosoma infection. PLoS Neglected Tropical Diseases 6:e1903.
- 33.
Anderson RM, May RM (1991) Infectious Diseases of Humans. Oxford: Oxford University Press. 757 p.
- 34. Sturrock RF (1973) Field studies on the transmission of Schistosoma mansoni and on the bionomics of its intermediate host, Biomphalaria glabrata, on St. Lucia, West Indies. Int J Parasitol 3:175–194.
- 35. Chan MS, Guyatt HL, Bundy DA, Booth M, Fulford AJ, et al. (1995) The development of an age structured model for schistosomiasis transmission dynamics and control and its validation for Schistosoma mansoni. Epidemiol Infect 115:325–344.
- 36. Chan MS, Guyatt HL, Bundy DAP, Medley GF (1994) The development and validation of an age-structured model for the evaluation of disease control strategies for intestinal helminths. Parasitology 109:389–396.
- 37.
Cornell SJ (2010) Modelling stochastic transmission processes in helminth infections. In: Michael E, Spear RC, editors. Modelling Parasite Transmission and Control: Springer. 66–78.
- 38. Pugliese A (2000) Coexistence of macroparasites without direct interactions. Theor Popul Biol 57:145–165.
- 39. Prichard RK, Basáñez MG, Boatin BA, McCarthy JS, García HH, et al. (2012) A research agenda for helminth diseases of humans: intervention for control and elimination. PLoS Neglected Tropical Diseases 6:e1549.
- 40. Churcher TS, Ferguson NM, Basáñez MG (2005) Density dependence and overdispersion in the transmission of helminth parasites. Parasitology 131:121–132.
- 41. Warren KS, Mahmoud AAF, Cummings P, Murphy DJ, Houser HB (1974) Schistosomiasis mansoni in Yemeni in California: duration of infection, presence of disease, therapeutic management. Am J Trop Med Hyg 23:902–909.