Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Parasitology. Author manuscript; available in PMC 2012 February 3.
Published in final edited form as:
PMCID: PMC3271124

A new approach to modelling schistosomiasis transmission based on stratified worm burden



Multiple factors affect schistosomiasis transmission in distributed meta-population systems including age, behaviour, and environment. The traditional approach to modelling macroparasite transmission often exploits the ‘mean worm burden’ (MWB) formulation for human hosts. However, typical worm distribution in humans is overdispersed, and classic models either ignore this characteristic or make ad hoc assumptions about its pattern (e.g., by assuming a negative binomial distribution). Such over simplifications can give wrong predictions for the impact of control interventions.


We propose a new modelling approach to macro-parasite transmission by stratifying human populations according to worm burden, and replacing MWB dynamics with that of ‘population strata’. We developed proper calibration procedures for such multi-component systems, based on typical epidemiological and demographic field data, and implemented them using Wolfram Mathematica.


Model programming and calibration proved to be straightforward. Our calibrated system provided good agreement with the individual level field data from the Msambweni region of eastern Kenya.


The Stratified Worm Burden (SWB) approach offers many advantages, in that it accounts naturally for overdispersion and accommodates other important factors and measures of human infection and demographics. Future work will apply this model and methodology to evaluate innovative control intervention strategies, including expanded drug treatment programmes proposed by the World Health Organization and its partners.

Keywords: models, theoretical, schistosomiasis, disease transmission, infections, environmental and behavioural risk


Helminth infection typically involves non-linear systems that link, human populations intermediate hosts and multiple environmental factors contributing to efficient transmission in endemic areas. For policymakers, capturing the dynamics of this process has proven challenging in terms of developing predictive models that have quantitative accuracy. A traditional mathematical approach to modelling macro-parasite systems (Macdonald, 1965) exploits a mean worm burden (MWB) formulation for human hosts and a suitable infection prevalence estimate for intermediate invertebrate hosts. However, worm distributions among host populations often exhibit an overdispersed pattern. In the standard modelling approach, overdispersion is either ignored or over-simplified using ad hoc assumptions about its distribution (e.g., the negative binomial distribution). Furthermore, the pattern of infectious burden is assumed to remain static (unchanged) over time. This approach cannot explain overdispersion or justify its ‘stationarity’, and thus leaves many uncertainties about its application to observed epidemiological data from the field and its predictions about the impact of control programmes.

We propose a new approach to modelling schistosomiasis transmission (or any vector-borne macro-parasite), based on a Stratified-Worm-Burden (SWB) system as an alternative to the classical MWB approach. Several papers have studied such or similar systems for macro-parasites in the past (Pugliese,2000; Feng et al. 2001; Rosa and Pugliese, 2002), but most of this work aimed at deriving reduced (MWB-type) equations, or their extended versions that couple MWB to higher moments (variance and/or aggregation) of SWB distribution. Among several examples are papers on (i) regulation and stability of host-parasite populations (Anderson and May, 1978a, b); (ii) aggregation for single parasite species and multi-strain competition, coexistence, invasion (Dobson and Roberts, 1994; Pugliese, 2000; Rosa and Pugliese, 2002) and (iii) parasite strain diversity and evolution of drug resistance (Feng et al. 2001).

In contrast, for the present study, we developed and implemented an extended SWB approach to model, as necessary, the complexity of schistosomiasis transmission for structured host populations in a geographically distributed environment with varying water contact patterns. We take advantage of the extended model systems to consider age-structured communities made of children and adult groups having suitable rates of birth, mortality, and ageing.

Our approach extends the standard Macdonald-type MWB and other prevalence-based models (e.g., Riley et al. 2008). We gain several advantages by combining both infection intensity and prevalence in the model. Primarily, we can naturally account for overdispersed worm burden in host populations (even for homogeneous behavioural risk subgroups) without ad hoc assumptions on distribution patterns commonly used in MWB-type approach such as the negative binomial (Anderson and May, 1991; Alexander et al. 2000). Inclusion of additional features reflecting the distribution of human infection (the prevalence of infected/non-infected, or lightly infected/heavily infected, and mean egg counts (as a measure of infection intensity)), proves particularly helpful for comparison of model outputs with the observed field data for calibration/validation purposes. In using the model, these variables can be derived from data or entered as dynamic variables. In addition, our model accommodates host demographics, so that one can study the effect of multiple factors (both stationary and variable) on transmission and long-term outcomes of infection. These changing variables can include population changes, environmental/behavioural modifications, snail variability, and intermittent control interventions.

Some earlier works (Chan et al. 1996; Medley and Bundy, 1996; Chan and Bundy, 1997) used age-structured MWB to assess infection-associated morbidity and the effects of drug treatment for idealized host populations. The MWB approach, however, severely limited our ability to project outcomes for promising control interventions that target both age and geographical risk factors together (see e.g., Gurarie and King, 2005, 2008; Gurarie and Seto, 2009). The new SWB approach opens new possibilities for modelling the outcomes of control programmes that focus interventions by targeting high burden strata by age or location.


We start with an outline of general SWB-scheme, its implications for burden/prevalence distribution, and the analysis of endemic equilibria. Then we proceed to develop a model calibration scheme that exploits typical demographical and epidemiological infection data to estimate the baseline human-to-snail and snail-to-human transmission rates. These models are built upon each other, starting with a model describing a single homogeneous population, adding an age-structured (adult–child) population, and finally a geographically distributed (human–snail) transmission environment.

Our choice of calibration scheme is partly motivated by the available data from the Msambweni region of eastern Kenya, which were collected in several multi-year studies (Sturrock et al. 1990; Muchiri et al. 1996; Hamburger et al. 2004; Kariuki et al. 2004; Clennon et al. 2006). In a sequel paper we will apply our methodology (SWB model and calibration scheme) to this specific environment, comprised of 10 villages, each one subdivided into adult and child subgroups with 5 shared snail sites. When calibrated and validated, this model will then allow us to make projections on the transmission of vector-borne disease, and to simulate the effects of hypothetical control interventions based on variously targeted, single or combined drug therapy- and vector control-based schemes.

Infection strata

The host population H is divided into burden strata {hi} based on the number of parasites they carry.

Infection level01
Worm burden0< w<w1w1 <w<w2
Population stratah0h1

So each group hi viewed as a local subpopulation carries a parasite range wi < w < wi+1. The dividers {wi} can be chosen arbitrarily, the simplest (extreme) choice being ‘zero-nonzero’ (‘infected’ – ‘non in-fected’) strata with prevalences {h0, h1}. Such partition would result in the standard SIR-type prevalence model. The SIR methodology (h0; h1) however, is too crude for macro-parasite transmission, where force of infection depends primarily on the worm burden of the host population, rather than the infection frequency. On the opposite extreme is continuously distributed worm burden {h(w):0≤w < ∞}, which has limited practical applications outside math-1 ematical analyses. Our stratified worm burden (SWB) models extends the SIR (prevalence-based) modelling to accommodate both the worm burden and infection frequency relevant to population dynamics of macroparasite infections using the Macdonald-type ‘mean burden’ formulation, as explained below.

In what follows, we shall use uniform strata wi = iΔw; (i = 0, 1, 2, …) counted in steps Δw worms. The choice of Δw is related to infectivity threshold wT – the ‘effective’ number below which a population is considered free of infection. It is convenient to take wT a multiple of Δw; thus wT = Δw means h0 = ‘uninfected strata’, wT = 2Δw the first two strata (h0 + h1) = ‘uninfected’, etc. Indeed, host infectivity is proportional to the number of mated worm pairs, which in turn depends on their mating probability [var phi] (w) (Nasell, 1978; see also Anderson-May, 1991). Derivation of such a probability function is based on ‘negative-binomial’ assumption, suggesting a some-what low threshold value wT = 5 (independent of aggregation). On the other hand, experimental/field data for primate infection (Wilson, 2006) predict a substantially higher threshold for egg detection(wT ≈ 40). We assume infectivity threshold to lie somewhere in between 10≤wT≤20, and numerically explore 2 values: wT = 10 and wT = 20. Accordingly, we choose to convert the continuous values to discrete values at 10 worms per step (discretization step, Δw). While the choice of step Δw is not essential for calibration and analysis (as explained below), the choice of threshold wT has significant effect (see Fig. 6).

Fig. 6
Sensitivity analysis of model parameters (snail-to-human transmission rates) {Ai} for children in 10 villages: panel (A) shows heavily infected villages nos 6 and 7, and panel (B) for all others. Two sets (light/dark) of box-and-whisker plots (‘box’ ...

For higher burden levels (w > 10) we assume that host infectivity, as measured by egg counts, is roughly proportional to wi. The choice of uniform strata is not essential from computational standpoint, but has some mathematical advantages. We used 15 strata (i.e. burden levels up to 150 worms/host) in our analysis, which, according to our data set, is sufficient to cover a broad range of human infection, even heavily infected populations.

Fig. 1 shows the diagram of transitions among strata with ‘up-arrows’ → designating the force of human infection λ (determined by infected snail/larval densities, water contact patterns, probability of worm establishment, etc.), and ‘down arrows’, ← indicating burden resolution λi (due to worm mortality or other causes, including drug treatment). Different strata i have specific resolution rates γi that increase with burden level. Additional source terms, Si, represent new recruits to the infection pool (newborns, migrants into the area, and transition from younger to older age groups). Loss rates, μ, combine host mortality as well as out-migration, aging, etc. In some papers μ’s are considered to depend on level i to account for ‘burden-specific’ mortality. We ignore such disease-specific mortality as quite small, i.e., not having impact on population levels or community levels of schistosomiasis (even its more pathogenic S. japonicum species (van der Werf and de Vlas, 2001)). So our μ are purely demographic and considered identical for homogeneous populations.

Fig. 1
Infective strata and transitions of SWB model. Population strata (numbers or fractions of total) {hi} move ‘up’ or ‘down’ the scale at certain rates, due to worm acquisition (λ) and resolution (γi). Each ...

For age-structured populations (e.g., children and adult subgroupings) each group has its own infection strata, but their source terms differ. The children’s source terms come primarily in the lowest (uninfected) strata S0 (as newborns are free of infection). But the adult source, {Si}, will be distributed over all strata and depend on the burden distribution among children by the time they age into adulthood.

For a single (homogeneous) population H=0nhi, variables {hi(t)} obey a system of differential equations (SWB)


We used the following notations (i) transition rate λ is proportional to (per capita) force of infection: λ=AzΔw, made of ‘transmission coefficient’/‘infected snail’ – A, infected snail density – z, and Δw – ‘burden increment’ in the definition of strata {hi} (see Appendix A and Table 1); (ii) resolution rates {γi = γi} are proportional to burden levels {wi} of each strata (γ being natural worm mortality); (iii) host loss rate(s) = μ and new recruit sources = {Si}. Hereweassume identical force of infection and loss for all strata, i.e., a homogeneous community in terms of behaviour and environmental factors. This, however, does not preclude uneven (aggregated) burden distribution within system (1). Fig. 2A shows a few equilibrium distributions of human prevalence {hi(t)} (1) for several values λ contributed by a stationary source S0 at the lowest (uninfected) stratum (a typical pattern for child group). Fig. 2B illustrates dynamic relaxation (to equilibrium) of the coupled SWB-system (1) and snail prevalence y(t) (see Appendix A). Here human force of infection λ = Ay and snail force Λ is proportional to the cumulative worm burden (2) of human population. For dynamic simulation Fig. 2B we initialized the coupled system at a low-level human infection (light bars) – a hypothetical effect of ‘drug clearing’, and let it relax to its natural equilibrium aggregated at higher burden levels (black bars).

Fig. 2
(A) Equilibrium worm burden distribution for different values λ = 0·2 (light grey), 0·6 (medium grey) and 1 (black), with identical human loss rate μ = 0·02 and stationary source at the lowest strata. (B) Time relaxation ...
Table 1
List of notations, variables and parameters

System (1) can be interpreted as cumulative out-come of a stochastic process of ‘worm acquisition-death’ on individual level. Indeed, the RHS of system (1) has a matrix made of ‘subdiagonal’ λ and ‘diagonal’ {λj} that can be viewed as discretized Markov generator of such stochastic process. In this sense the standard Macdonald (1965) MWB gives the ‘statistical mean dynamics’ for such process (as function of time), whereas our SWB – its complete probability distribution function (PDF).

Human system (1) can be coupled in the usual way to snail infection (pre-patent and patent groups), through the force of human infection λ and the corresponding, extra-human, larval densities (Appendix A). Here the force of snail infection Λ will depend on human MWB, similar to the standard Macdonald (1965). But, unlike the Macdonald system, our MWB is not an independent variable but computed from prevalence variables {hi}



In the follow-up analysis and calibration we rescale SWB-system (1) and MWB (2) using non-dimensional (primed) variables and parameters, along with time and source terms. Specifically,

ww=wΔw={0,1,2,}– dimensionless
hihi=hi/H– prevalence strata
SiSi=Si/(ΔwH)– source terms in
prevalence strata

Then the ith equation for system (1) becomes


For the sake of brevity we shall drop ‘prime’ from all rescaled variables.

The SWB system yields several reduced equations based on statistical moments of distribution {hi}. Among them are (i) ‘total population dynamics’ for H (0-moment (sum)),


with S=Hi=0nSi; (ii) Macdonald-type MWB equation for w (1st moment (mean))


Here S¯=i=0niSi designates the cumulative source. Let us note that S¯=0, if new recruits enter system (1) only through its non-infected (zero) compartment (like newborn children). Derivation of (4) exploits a special form of resolution rates γi = γi for i-th strata). So worms in all hosts are assumed to have the same (mean) life-span, independent of the host state (age, burden level etc.). We also ignore all other confounding factors, like density-dependent mortality (due to crowding), or immune regulation of parasites. Such factors can in principle be included in the SWB scheme, but to calibrate such a model would require a more detailed demographic/infection data of host population (e.g. age, behavioural risk, susceptibility etc.), than provided in most community studies of schistosomiasis. Our formulation tries to keep what we believe are the most essential components of the system, suited for the most frequently available data.

Equation (4) looks like the standard Macdonald MWB for dimensionless w(t) with an additional term hn that disappears in the infinite limit systems (n → ∞), used in many theoretical studies. In our case (finite size SWB) this term remains positive, but sufficiently small over a wide range of parameter values, to be safely dropped in the follow-up analysis and calibration.

Equations (3) (total population growth) and (4) (MWB dynamics) represent the first two statistical moments of distribution {hi} driven by SWB system (1). It can be further extended to higher moments: variance, aggregation etc. Such schemes were adopted in some theoretical studies (e.g. Pugliese 2000) as a way to reduce a large (infinite) system (1) to its low-dimensional (mean-field) approximation. In this paper, our goal and approach is different. We utilize reduced (moment) equations (3)–(4) for model calibration. But all dynamic simulations (for prediction-analysis-control) are run with full SWB system (1), using its calibrated parameters. Thus a system of 10 villages in a particular study area (Msambweni), with each population subdivided into child and adult groups will contain 2 × 10 × 15 ‘human’ variables plus a suitable number of ‘snail’ variables proportional to the number of contact sites.

It should be stressed that increased order (dimension) of SWB system (compared to its MWB cousin) posed no substantial mathematical or computational obstacle by itself. Indeed, we can generate such a large, coupled dynamical system using Wolfram Mathematica 7 (Champaign, IL). It also allows efficient ways to solve and manipulate them in various ways for analysis, calibration, prediction and control experiments. Details of these programs are available from the corresponding author.

Analysis of equilibria

The equilibrium solution of system (1), with stationary source S and prescribed λ, depends on 2 parameters μ and λ. The former μ is determined by known demographical or biological factors, the latter λ (proportional to the force of human infection) can vary depending on the transmission environment, and will be calibrated along with other parameters.

The complete SWB system (1) has no exact analytic solutions (even for stationary equilibria), but we can explore a broad range of parameters {λ, μ} numerically. Fig. 2A shows equilibrium burden distributions at several λ. Qualitatively they resemble negative binomial (NB) distributions of decreased aggregation (smaller λ having higher aggregation, hence a larger fraction of low-level infections).

MWB equation (4) can be used either in dimensional form for variable w(t), or non-dimensional w=1nihi. Our calibration scheme exploits 2 epidemiological data values, uninfected human prevalence P = h0, and MWB, w (which is related to measured egg-counts in the field). Both are functions of parameters λ, μ, and source S


Both functions [var phi] and ψ are illustrated in Fig. 3. While there are no analytic formulae, they can be easily computed numerically by interpolation (sampling a range of values λ and μ). These procedures are used extensively in our calibration schemes below.

Fig. 3
(A) Plot of uninfected prevalence function P = [var phi](λ, μ, S) on worm acquisition rate λ for 3 values μ = {0·1γ, 0·2γ, 0·3λ} (top to bottom) and identical stationary source ...

Let us note that the range of uninfected human prevalence from our data (Table 2), 0·2 < P < 1, constrains λ/γ within relatively narrow interval [0, 2·5] (Fig. 3A). In that range, both λ and w obey a nearly linear (Macdonald type) relationship (Fig. 3B),


which is important for our calibration scheme. At large λ > 8 (a strong force of infection) we observe nonlinear departure from simple Macdonald model outcomes (Fig. 3B), due to the extra-term of our MWB-equation (4). But this range (and level of transmission) is well above the likely range of values in our study environment.

Table 2
Human and snail demographics and infection data: ‘egg-count’ refers to arithmetic mean/host (Prevalence, infected fraction of population.)


Data that are typically available for calibrating models of macro-parasite transmission include: (i) Demographics: human population H, and snail density N; (ii) Biological factors: human loss rates (mortality, turnover) – μ, snail mortality – ν;, worm mortality λ; (iii) Epidemiological infection data: infected and patent snail prevalences (y, z), human egg-counts (arithmetic and geometric means) E, G. Depending on the choice of model, all or part of these data can be used for estimating unknown transmission parameters (see, e.g. Riley et al. 2008).

In our calibration scheme below, we use both egg-counts (arithmetic mean) and infection prevalence as 2 independent parameters. For the sake of exposition we precede the discussion of larger, distributed systems, with a few special cases.

The general scheme involves 3 steps:

  1. Estimate human force of infection λ (for suitable populations) from demographic and human prevalence data. This part utilizes prevalence function [var phi](λ, μ, S)
  2. Next estimate snail-to-human transmission coefficient {A} for appropriate population groups, as well as biological parameter ρ (egg-production/worm).
  3. Finally, estimate human-to-snail transmission coefficients {B} (see Table 1 for notations).

The single population case

In this case our calibration scheme will utilize 2 reduced equations resulting from SWB (1): (i) MWB (4) coupled to a snail prevalence equation (as in the standard Macdonald setup), and (ii) an equilibrium zero-prevalence function (5), computed numerically. The former gives a coupled system of 2 differential equations for w(t) and y(t)


with per capita forces of infection


The transmission coefficients A (snail to human) and B (human to snail) depend on several biological and environmental/behavioural factors (Table 1). Let us note that having stratified human population is irrelevant for snail infection, as long as all strata have similar behavioural (water contact) patterns. So the force of snail infection Λ is the same as in the classical Macdonald setup (i.e., proportional to total human burden, Hw)

System (7) has stable endemic equilibrium provided its Basic Reproduction Number (BRN)


Our goal is to estimate transmission coefficients A, B from equilibrium solutions (7) and the available data. The data include egg-count E (related linearly to w) and snail prevalence (y). The relation of E to w is given by another (unknown) parameter egg production rate


Our goal is to calibrate parameters A, B, ρ. We proceed in 2 steps:

  1. solve equation (5) for a given (known) human prevalence P, relative loss rate μ, (demographics), and source S, to find rescaled force of infection λ.
  2. compute A from the first expression of (8), ρ and B equilibrium of equations (7)

The resulting explicit formulae give calibrated parameters in terms of the data


In reality, infection field data, either human (egg-counts E, prevalence P) or snail data (pre-patent and patent prevalences y, z), can be highly inaccurate (de Vlas et al. 1997; Hamburger et al. 2004) due to low sensitivity of the standard tests and observational errors. One way to account for data errors (underestimation) is through a suitable statistical ensemble of data values and the resulting parameters (A, B, ρ). Here explicit formulae (9) prove helpful, as any a priori information of data errors (statistics or distribution) can be translated via (9) into suitable estimates for calibrated parameters. Indeed, we shall explore sensitivity of our calibration scheme to data errors for a specific data set (below).

Multiple risk groups

Next we subdivide human population within a single village (or human–snail site pairing) into several risk groups labelled by index s, that differ in their behaviour (contact rate) hence the resulting force of infection, λs. In the following, those will be children and adult groups. The extended system couples human transmission (either MWB equations for dimensional {ws} or the corresponding SWBs for each risk-group) and a single snail equation (contributed by all groups)


As above, we use equilibrium solutions (10) along with the prevalence equation [var phi](λs, μs, Ss) = Ps, for each group. The calibration problem now asks to estimate transmission coefficients {As, Bs}. The above 2 steps yield (i) the rescaled force of infection for each group {λs}, (ii) the corresponding A-coefficients


To get B-coefficients we observe that both sets (A’s and B’s) are proportional to contact rates {ηs}, assumed to be the only risk factor,


Here a is worm establishment rate/per single infected snail, assumed to be identical for all groups, b – snail infection rate/per single worm × contact. Then we get B-coefficients


expressed through group specific egg-production/worm rates


Equations (11)–(13)–(14) extend the above calibration scheme (9) to heterogeneous populations with multiple risk groups.

Child–adult system

Next we further specify the above ‘multi-risk’ setup for separate child–adult groups, still within a single human–snail pairing context. These two human age-related subgroups differ in their risk factors (contact behaviour) and infection levels. This case, however, differs from the previous one, in that the two groups are also coupled by a demographic birth–aging process, whereby infection is carried from the younger to the older group. The population dynamics forvariables Hc; Ha (total child/adult populations) is given by a Leslie-type differential system


The demographic parameters include (i) mean per capita birth rate β, (ii) specific child and adult mortalities μc, μa, (iii) transition (aging) rate τ.

Note that population dynamics (15) follow from the extended child–adult SWB system, the same way that the previous simple population model (3) follows from (1). The source term for child SWB, Sc = βHa, enters through the uninfected (zero) stratum (new-borns), whereas the adult source, Sa, is equal to children’s SWB distribution ((Eq.B.1) in Appendix B).

The real population and transmission dynamics are more involved. Clearly, different ages have different SWB-distributions, gradually changing from early childhood (aggregated at low burden) to older age-group patterns shifted towards higher w (Fig. 2). The adult source would then come from the older children’s distribution (15–20 years), rather than the mean distribution over all ages (5–20). In principle, both systems (demographic and SWB) can be extended to accommodate refined population structure (e.g. 5-year bins). But here we confine ourselves to the coarse child–adult partition, as the only available epidemiological data for our specific environment has only those two groups.

The corresponding child–adult Macdonald-type MWB equations account for demographics


As in the previous (multiple risk) case, we utilize zero-prevalence to compute forces of infection,


Specifically, we first find children’s force of infection λc with known prevalence Pc and source Sc (in terms of birth rate), then compute the resulting equilibrium-prevalence distribution h^c={hic}i=0n, substitute it to get adult source Sa={τhicHc/Ha}i=0n and get λa from the second equation in (17). Having estimated both λ’s we proceed as in the previous section to get transmission coefficients {Ac; Aa; Bc; Ba; ρc; ρa} for two groups (see Appendix B for details).

Distributed system

Our next goal is to extend the above models and calibration scheme for geographically distributed systems made up of several different human sites (villages), labelled by i = 1, 2, …, n, and connected snail (water), sites labelled by j = 1, 2, …, m. Mathematically, calibration of such dynamical system is an ill-posed (heavily under-determined) problem as far as the available human and snail infection data. Indeed, the n + m equations of a typical ‘MWB + snail prevalence’ system (Eq.C.1) have altogether (2n × m) unknown transmission coefficients {Aij; Bij}, if all connections are allowed (Woolhouse, 1991; Woolhouse et al. 1991; Gurarie and King, 2005; Riley et al. 2008). Different reduction schemes can be used to constrain the resulting parameter set. Thus Gurarie and King (2005) and Gurarie and Seto (2009) proposed a version of a landscape ‘gravity’ model, whereby geographical distances (or other hurdles) control human water contact rates. A suitable parameterization of this hurdle-function yields a drastic reduction of the parameter set.

Other procedures are possible depending on the available data. For example, data collected from our specific environment (Msambweni region of eastern Kenya (Sturrock et al. 1990; Muchiri et al. 1996; Hamburger et al. 2004; Kariuki et al. 2004; Clennon et al. 2006)) has additional information in terms of relative contact rates from various human sites to water (snail) sites. Here we propose a distributed calibration scheme designed specifically for such type of data. Each contact pair (i, j) contributes to human infection in village i and snail infection at site j through transmission coefficients Aij; Bij both pro-portional to absolute contacts per person from village i to water site j, Tiηij with Ti being the group-specific total contacts to all water sites and ηij relative contact (jηij=1). Our basic hypothesis is that each of them is a product of a suitable local (biological/environmental) factor (A0 times ‘absolute contact rate’),


It means that humans at site i will pick up infection from snail site j in proportion to their ‘absolute contact’ × ‘the site’s snail infection level’. By the same pattern, snails at j will get infected by humans from different sites i in proportion to their absolute contact times human infection level (egg release). Since the establishment of miracidia in snails involves 2 stages – hatching of miracidia and the miracidium’s successful invasion into a snail body – which depends heavily on the local environmental conditions, we set the local environmental factor B0j as site specific. This hypothesis clearly reduces the number of unknown parameters to a manageable set {Ai; Bj}, consistent with the number of MWB equations.

Accordingly, for age-structured populations, the MWB system ((Eq.C.2) in Appendix C) has the following unknown parameters to be estimated: (i) local transmission coefficients for child-adult groups {Aic;Aia}; (ii) local transmission coefficients for snail sites {Bj}; (iii) egg production/release rate for children and adults {ρc; ρa}, the latter two are biological factors independent of location.

The calibration scheme proceeds similarly to the previous case. Namely, we first use local prevalence and loss rate for children and adults to get their local rescaled forces of infection, {λic;λia}, as in Table 3. Then we write the proper distributed versions of the calibration formula (11), (14)

Table 3
Calibrated transmission coefficients and egg production rates based on the individual level field data for 2 choices of discretization step (Δω), and infection threshold ωT = 10.

Symbol s = {c, a} indicates child/adult groups, and denominator of A gives a weighted snail infection density (with relative ‘contact weights’). While all A-coefficients are local (i.e., depend on i), the egg production term, ρs, should be biological and site independent. So, based on (19) we expect the numerator of the right hand site of ρs expression to be proportional to the denominator at all sites, with slope ρs. The real data show departures from simple linear relations, so we approximate parameters ρc; ρa with their ‘best-fit’ slopes. Having estimated {Ais} and {ρs}, we compute coefficients {Bj} from the snail equilibrium equations ((Eq.C.4) in Appendix C).


Calibration results for a specific environment

In this section we apply SWB methodology to a specific environment of the Msambweni region of Eastern Kenya sampled from 10 villages and 5 contact snail sites observed in year 2000, shown in Fig. 4 (Sturrock et al. 1990; Muchiri et al. 1996; Hamburger et al. 2004; Kariuki et al. 2004; Clennon et al. 2006). Infection prevalence and egg-counts (arithmetic means) were observed using standard filtration of urine samples (10 ml) from the resident populations (Table 2). Basic biological parameters were included from well-known estimates of snail and worm mortality, as well as basic life table estimates (among other parameters in Table 1). Field data (Table 2) for the models included (i) demographic data: human populations (divided into child/adult groups) and snail densities; (ii) infection data: mean egg count for children in each village and density of susceptible, infected and shedding adult snails in all water sites; and (iii) behavioural data: relative contact distribution for each village to adjacent snail sites.

Fig. 4
Partial area map of the Msambweni region (left), and its schematic view (right) with human sites (black dots), snail sites (grey dots) and established ‘human-to-snail’ links based on observed human contact activity.

We first demonstrate calibration results of the individual level field data, utilizing formulae and procedures of the previous methods sections (Model Calibration). As above we run it in 3 steps:

  1. Estimate ‘local forces’ of human infection in 10 villages (first children, then adults) from the human data (prevalence and demographics), by solving equations at each site. Here we utilize the known (‘universal’) childhood prevalence function [var phi](λ, μ, Sc) (Fig. 3) and the derived adult prevalence function [var phi](λ, μ, Sa) with source Sa from the children’s group.
  2. Next we estimate village-specific transmission coefficients Ai(i = 1, 2, …, 10) for child/adult groups via formulae (19), and forces of infection λi (output of step (i)). We also compute the ‘bestfit’ estimates for biological parameters ρc; ρa (egg-production/worm).
  3. Finally, we estimate snail transmission coefficients Bj(j = 1, 2, …, 5) from equation (Eq.C.4) (Appendix C).

Table 3 shows estimated snail to human transmission coefficients {Ac, Aa} for child and adult groups in 10 villages, along with the best fit egg-production rate estimates {ρc, ρa} for the two groups, and human to snail transmission rates {Bj: j = 1, …, 5}. We also looked at the effect of worm discretization step Δw on estimated parameters. To this end we run 2 calibrations of the individual level field data, with steps Δw = 10 and Δw = 5, and fixed infectivity threshold wT = 10 (Table 3). The 2nd case gave overall higher estimates of transmission rates {Ai}, but the children group was less sensitive (< 15%) than adults. Since children carry the bulk of infection and transmission, our scheme seems overall insensitive to Δw.

Model validation

We ran a simple self-consistency (validation) test for our calibrated system by comparing the predicted equilibria (computed from the model) with the individual level field data. An exact match is not expected as parameters of the distributed system are no longer exact solutions of equilibrium equations. Indeed, the complete set of equations (MWB + human prevalence + snail prevalence) is over-determined now, so we use best-fit approximation to obtain the estimates of Table 3. The comparison test shows good agreement (Fig. 5).

Fig. 5
Comparison between individual level field data (dark disk) and equilibrium solutions of the calibrated model (grey square); the x-axis labels the ‘village’ or snail site and y-axis – ‘mean egg count’ or ‘snail ...

Sensitivity analysis

The infection data based on egg and snail counts is notoriously inaccurate partly due to imperfect diagnostic tools. Another major uncertainty in our model is the infectivity threshold wT. To examine such uncertainties on calibrated parameters we took 2 threshold values wT = 10 and 20, and in each case run an ensemble of 200 data choices, randomly selected within suitable ranges: ±10% of the individual level field data, for egg counts; ±5% – for human prevalences, and ±20% – for snail infection. The error margins are somewhat arbitrary, but they reflect the overall level of uncertainty. For instance, individual egg-counts can be much ‘noisier’ than 10%, but the community mean would reduce their variance byfactor 1/n – community size. Fig. 6 shows the results of sensitivity analysis for some calibrated parameters. The overall estimates remain fairly stable (error <5%). We plan to use the calibrated ensemble of parameters to provide error bars for predicted control outcomes in the future work on control analysis with SWB.


We propose a new approach to modelling schistosomiasis transmission (and other vector-borne macroparasites), based on a Stratified Worm Burden (SWB) formulation. In this approach, a community of hosts is subdivided into burden strata in terms of their infection levels. The transitions among strata depend on (i) force of infection (determined by the transmission environment and human risk factors); (ii) worm attrition (natural or drug-induced); and (iii) human demographics (ageing, mortality, growth sources). Such stratified systems offer many advantages compared to the conventional MWB approach (Macdonald, 1965), or SIR-type prevalence systems. By naturally encoding both types of information, these models allow new (refined) calibration procedures with the available data. Furthermore, they can also include demographics and/or other measured (e.g., behavioural) risk factors for structured human populations (age, occupation, etc.). Finally, SWB naturally explain highly aggregated burden distributions in host populations (typical of Schistosoma and other helminth infections (Anderson and May, 1991)) without imposing ad hoc distribution patterns.

We have demonstrated the SWB approach and model calibration procedures in several settings: the single homogeneous population, a population with multiple risk/age groups (child–adult systems), and a more complex area with geographically distributed human–snail populations and contact patterns. In each case, the modelling combines demographics (population levels, births, deaths, and ageing), human–snail infection levels (egg-counts and prevalences), and for distributed systems, a minimum of environmental data (snail habitats, human contact patterns, etc.).

In simple cases such dynamical (differential) systems can be analysed mathematically to determine their endemic equilibrium levels in terms of transmission parameters such as the Basic Reproduction Number, R0. In all cases they are amenable to efficient numeric implementation using desktop software, despite a large number of variables and parameters. For instance, a complete list of variables for a distributed human–snail system would typically include number of population sites × population age/risk groups × burden strata + all snail variables. Here, we were able to implement it using accessible commercial desktop software.

The SWB methodology has many potential applications. In the present paper, we outlined calibration results for a specific environment (the Msambweni region of eastern Kenya). We also conducted a simple validation test and sensitivity analysis of our calibration scheme to uncertainties in the data (random errors). In subsequent work we shall provide further details and applications of our methodology inexploring the predicted impact of different control strategies such as the currently advocated targeted treatment of school-aged children (WHO, 2006) as compared to community-based delivery to children and adults (Gabrielli et al. 2006).

One might wonder how our extended model and results compare to the standard (Macdonald) MWB approach. There is no direct way to compare the two systems, as MWB by itself has no human prevalence information. Other significant differences can arise in implementation of drug treatment in the two systems. The drug-based clearing of infection can enter an MWB-system only through increased worm mortality, whereas in SWB we implement it by shifting burden strata from ‘high’ to ‘low’ to ‘lowest’ levels as hosts get cleared. The significant difference has to do with the slow process of natural worm mortality versus fast worm clearing by drug. We shall discuss this aspect in more detail in the forthcoming work on control interventions.

While our approach allows more detailed parasite distribution in host populations, it leaves out many other important within-host factors, such as density-dependent survival and egg-production and immune regulation of parasite levels and establishment. Some of these factors (immunity) are age-dependent, and would require more detailed age-structured (doubly stratified) SWB. In the current analysis, we had limited parasitological data (child–adult), and structured our model accordingly. Hence age-dependent factors (immunological etc) enter our system through different transmission and egg-production rates. Table 3 shows marked differences between two groups, and we attribute this fact partly to a protective role of immunity that develops with age.

It is, however, somewhat easier to accommodate density-dependent factors (worm survival and host infectivity) within SWB approach. Resolution rates {γm} (for m-th strata) related to worm mortality, can be made to increase faster (than linearly) in ‘m’ to account for crowding mortality. One can also make infectivity (egg-production) of m-th strata a nonlinear function of m in terms of reduced fertility. While there is no simple way to accommodate such features in MWB, we have demonstrated how SWB provides a natural representation of such macroparasite dynamics.


The authors thank the referees for their pointed questions and stimulating comments that helped us to improve the results and exposition of the work.


This project was supported in part by National Institutes of Health Research Grant R01 TW008067 funded by the Fogarty International Center. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


A. Basic Macdonald system for multiple human risk groups

All snail and larvae populations {N; c, m} are measured by [density]/[unit snail habitat]. Our basic assumptions include

  1. Snails are minimally mobile and relatively attached to their local habitat, so we ignore hydrologic (or other) transport of snails/larvae.
  2. There is an approximately uniform snail density throughout each habitat (their local numbers are part of the collected field data).
  3. All human–snail contacts and transmission (via larval densities) take place locally. So the force of human-to-snail infection (miracidial production term) is proportional to human population density/unit habitat, rather than total human population.

The basic Macdonald setup (Macdonald, 1965) for single habitat modified to include multiple risk groups (labelled by s), consists of (i) mean worm burden (MWB) and egg counts {ws, Es}; (ii) susceptible-infected-patent snail densities (or prevalences) {x, y, z}; (iii) cercarial and miracidial densities {c, m}:


The miracidial production term is given by Qs = πMηsHsEs – proportional to the population density (per unit habitat).

When specified for child–adult groups (the human part of the MWB system) takes the form


Here we allow for aging, i.e. transition of worm burden from ‘child’ to ‘adult’ compartment at the age of 20. Hence we estimated loss-rate (τ = 0·05) in the child-equation, balanced by a gain term, τHcHa, in the adult-equation, weighted in proportion to population afraction Hc/Ha. In this formulation we assumed stationary populations, but the equations can beeasily adjusted for variable (time dependent) case. The notations are listed in Table 1.

The complete system (Eq.A.1) can be reduced by eliminating larval equations, due to their ‘fast’ death rates (μC; μM) compared to other ‘slow’ processes. So we replace (c, m) with their quasi-equilibrium values



The resulting human–snail system has MWB-type dw human equations: dwsdt=AsNy(γ+(...))ws (for s-group), or specifying it for a ‘child-adult’ setup we get


The snail part assuming stationary (seasonal mean) total population N, consist of 3 prevalence variables: x-susceptible, y-infected (pre-patent), z-patent, that obey a system of equations


The snail force of infection Λ (contributed by all contact groups {s}) is given by


The transmission coefficients {As; Bs} (‘snail-to-human’ and ‘human-to-snail’) are shown in Table 1.

B. Child-adult SWB system

Differential equations for SWB system for children and adults


The condition for stationary population in demographic equations (15) is


We estimated τ and β from the demographic data of the Msambweni region (Sturrock et al., 1990; Muchiri et al., 1996; Hamburger et al., 2004; Kariuki et al., 2004; Clennon et al., 2006). The child-adult death rates were computed from


They depend on local populations, hence we make them site specific.

Having estimated rescaled forces of infection λc and λa, we get the transmission coefficients and egg release rate from the equilibrium solution of system (16)


C. Distributed system

MWB equations for distributed human–snail meta-population system:


MWB equations for an age-structured (child–adult) human population with implified transmission coefficients:


Here we use local patent snail densities {Zj} instead of infected snail prevalences of earlier special cases (formulae (9)–(11)–(13)). Since patency conversion rates of infected snails {ψj = Zj/Yj} are local and vary across snail sites (Table 2) we need to account of their contribution to human infection separately. Another modification is the force of infection for snail in (Eq.C.2).


Here {Uj} are relative sizes of snail habitats (based on the number of contact sites of Table 2). We divide ‘human population’ by ‘contact size’ (Hi/Uj), since humans distribute their infection in inverse proportions to contact size of each habitat. Also the mean egg count Ei is used in (Eq.C.3) instead of the mean worm burden wi in the snail force of infection Λ, as different age groups may have different egg-release rates. The transmission coefficients {Ai; Bj} in formula (Eq.C.3) are somewhat different from earlier cases (Table 1). Namely,



where Tis represents the total contact rate for the risk group s in village i to all water sites. Having estimated Aic, Aia, ρc and ρa, we compute Bj from the snail equation (2) at site j, using the relation between A’s and B’s equation (18):



  • Alexander N, Moyeed R, Stander J. Spatial modelling of individual-level parasite counts using the negative binomial distribution. Biostatistics. 2000;1:453–463. [PubMed]
  • Anderson RM, May RM. Regulation and stability of host-parasite population interactions. I. Journal of Animal Ecology. 1978a;47:219–247.
  • Anderson RM, May RM. Regulation and stability of host-parasite population interactions. II. Journal of Animal Ecology. 1978b;47:249–267.
  • Anderson RM, May RM. Dynamics and Control. Oxford University Press; New York, USA: 1991. Infectious Diseases of Humans.
  • Chan MS, Bundy DA. Modelling the dynamic effects of community chemotherapy on patterns of morbidity due to Schistosoma mansoni. Transactions of the Royal Society of Tropical Medicine and Hygiene. 1997;91:216–220. [PubMed]
  • Chan MS, Guyatt HL, Bundy DA, Medley GF. Dynamic models of schistosomiasis morbidity. American Journal of Tropical Medicine and Hygiene. 1996;55:52–62. [PubMed]
  • Clennon JA, Mungai PL, Muchiri EM, King CH, Kitron U. Spatial and temporal variations in local transmission of Schistosoma haematobium in Msambweni, Kenya. American Journal of Tropical Medicine and Hygiene. 2006;75:1034–1041. [PubMed]
  • Dobson A, Roberts M. The population dynamics of parasitic helminth communities. Parasitology. 1994;109(Suppl.):S97–S108. [PubMed]
  • Feng Z, Curtis J, Minchella DJ. The influence of drug treatment on the maintenance of schistosome genetic diversity. Journal of Mathematical Biology. 2001;43:52–68. [PubMed]
  • Gabrielli AF, Toure S, Sellin B, Sellin E, Ky C, Ouedraogo H, Yaogho M, Wilson MD, Thompson H, Sanou S, Fenwick A. A combined school- and community-based campaign targeting all school-age children of Burkina Faso against schistosomiasis and soil-transmitted helminthiasis: performance, financial costs and implications for sustainability. Acta Tropica. 2006;99:234–242. [PubMed]
  • Gurarie D, King CH. 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. 2005;130:49–65. [PubMed]
  • Gurarie D, King CH. Age- and risk-targeted control of schistosomiasis-associated morbidity among children and adult age groups. The Open Tropical Medicine Journal. 2008;1:21–30.
  • Gurarie D, Seto EY. Connectivity sustains disease transmission in environments with low potential for endemicity: modelling schistosomiasis with hydrologic and social connectivities. Journal of the Royal Society Interface. 2009;6:495–508. doi: 10.1098/ rsif.2008.0265. [PMC free article] [PubMed]
  • Hamburger J, Hoffman O, Kariuki HC, Muchiri EM, Ouma JH, Koech DK, Sturrock RF, King CH. Large-scale, polymerase chain reaction-based surveillance of Schistosoma haematobium DNA in snails from transmission sites in coastal Kenya: A new tool for studying the dynamics of snail infection. American Journal of Tropical Medicine and Hygiene. 2004;71:765–773. [PubMed]
  • Kariuki HC, Clennon JA, Brady MS, Kitron U, Sturrock RF, Ouma JH, Ndzovu ST, Mungai P, Hoffman O, Hamburger J, Pellegrini C, Muchiri EM, King CH. Distribution patterns and cercarial shedding of Bulinus nasutus and other snails in the Msambweni area, Coast Province, Kenya. American Journal of Tropical Medicine and Hygiene. 2004;70:449–456. [PubMed]
  • Macdonald G. The dynamics of helminth infections, with special reference to schistosomes. Transactions of the Royal Society of Tropical Medicine and Hygiene. 1965;59:489–506. [PubMed]
  • Medley GF, Bundy DA. Dynamic modeling of epidemiologic patterns of schistosomiasis morbidity. American Journal of Tropical Medicine and Hygiene. 1996;55:149–158. [PubMed]
  • Muchiri EM, Ouma JH, King CH. Dynamics and control of Schistosoma haematobium transmission in Kenya: an overview of the Msambweni Project. American Journal of Tropical Medicine and Hygiene. 1996;55:127–134. [PubMed]
  • Nasell I. Mating for schistosomes. Journal of Mathematical Biology. 1978;6:21–35. [PubMed]
  • Pugliese A. Coexistence of macroparasites without direct interactions. Theoretical Population Biology. 2000;57:145–165. doi: S0040-5809(99)91443-0 [pii] [PubMed]
  • Riley S, Carabin H, Belisle P, Joseph L, Tallo V, Balolong E, Willingham AL, Fernandez TJ, Gonzales RO, Olveda R, Mcgarvey ST. Multi-host transmission dynamics of Schistosoma japonicum in Samar province, the Philippines. PLoS Med. 2008;5:e18. doi: 10.1371/journal.pmed.0050018. [PMC free article] [PubMed]
  • Rosa R, Pugliese A. Aggregation, stability, and oscillations in different models for host-macroparasite interactions. Theoretical Population Biology. 2002;61:319–334. doi: 10.1006/tpbi.2002.1575. [PubMed]
  • Sturrock RF, Kinyanjui H, Thiongo FW, Tosha S, Ouma JH, King CH, Koech D, Siongok TK, Mahmoud AA. Chemotherapy-based control of schistosomiasis haematobia. 3. Snail studies monitoring the effect of chemotherapy on transmission in the Msambweni area, Kenya. Transactions of the Royal Society of Tropical Medicine and Hygiene. 1990;84:257–261. [PubMed]
  • Van Der Werf MJ, De Vlas SJ. Report for WHO Parasitic Diseases and Vector Contol. Erasmus University; Rotterdam, The Netherlands: 2001. Morbidity and Infection with Schistosomes or Soil-Transmitted Helminths; pp. 1–103.
  • de Vlas SJ, Engels D, Rabello AL, Oostburg BF, Van Lieshout L, Polderman AM, Van Oortmarssen GJ, Habbema JD, Gryseels B. Validation of a chart to estimate true Schistosoma mansoni prevalences from simple egg counts. Parasitology. 1997;114:113–121. [PubMed]
  • Wilson RA, Dam GJ, Kariuki TM, Farah IO, Deedler AM, Coulson PS. The detection limits for estimates of infection intensity in schistosomiasis mansoni established by a study in non-human primates. International Journal for Parasitology. 2006;36:1241–1244. [PubMed]
  • World Health Organization . Preventive Chemotherapy in Human Helminthiasis: Coordinated Use of Anthelminthic Drugs in Control Interventions: a Manual for Health Professionals and Programme Managers. World Health Organization Press; Geneva, Switzerland: 2006.
  • Woolhouse ME. On the application of mathematical models of schistosome transmission dynamics. I. Natural transmission. Acta Tropica. 1991;49:241–270. [PubMed]
  • Woolhouse ME, Watts CH, Chandiwana SK. Heterogeneities in transmission rates and the epidemiology of schistosome infection. Proceedings of the Royal Society of London B. 1991;245:109–114. [PubMed]