Hostname: page-component-74656f9884-rdgjk Total loading time: 0 Render date: 2026-11-08T05:50:17.344Z Has data issue: false hasContentIssue false
  • English
  • Français

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

Published online by Cambridge University Press:  13 July 2010

D. GURARIE ,
C. H. KING  and
X. WANG
D. GURARIE*
Affiliation:
Department of Mathematics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, Ohio 44106 USA
C. H. KING
Affiliation:
Center for Global Health and Diseases, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, Ohio 44106 USA
X. WANG
Affiliation:
Department of Mathematics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, Ohio 44106 USA
*
*Corresponding author: Department of Mathematics, 220 Yost Hall, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, Ohio 44106-7058, USA. Tel: 001 216 368 2857. Fax: 001 216 368 5163. E-mail: dxg5@cwru.edu
Rights & Permissions [Opens in a new window]

Summary

Background/Objective. 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 oversimplifications can give wrong predictions for the impact of control interventions. Methods. 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. Results. 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. Conclusion. 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

modelstheoreticalschistosomiasisdisease transmissioninfectionsenvironmental and behavioural risk

Information

Type
Research Article
Information
Parasitology , Volume 137 , Issue 13 , November 2010 , pp. 1951 - 1965
Copyright
Copyright © Cambridge University Press 2010

INTRODUCTION

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, Reference Macdonald1965) 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 oversimplified 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 macroparasite), 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, Reference Pugliese2000; Feng et al. Reference Feng, Curtis and Minchella2001; Rosa and Pugliese, Reference Rosa and Pugliese2002), 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, Reference Anderson and May1978a, Reference Anderson and Mayb); (ii) aggregation for single parasite species and multi-strain competition, coexistence, invasion (Dobson and Roberts, Reference Dobson and Roberts1994; Pugliese, Reference Pugliese2000; Rosa and Pugliese, Reference Rosa and Pugliese2002) and (iii) parasite strain diversity and evolution of drug resistance (Feng et al. Reference Feng, Curtis and Minchella2001).

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. Reference Riley, Carabin, Belisle, Joseph, Tallo, Balolong, Willingham, Fernandez, Gonzales, Olveda and Mcgarvey2008). 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, Reference Anderson and May1991; Alexander et al. Reference Alexander, Moyeed and Stander2000). 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. Reference Chan, Guyatt, Bundy and Medley1996; Medley and Bundy, Reference Medley and Bundy1996; Chan and Bundy, Reference Chan and Bundy1997) 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, Reference Gurarie and King2005, Reference Gurarie and King2008; Gurarie and Seto, Reference Gurarie and Seto2009). 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.

MATERIALS AND METHODS

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. Reference Sturrock, Kinyanjui, Thiongo, Tosha, Ouma, King, Koech, Siongok and Mahmoud1990; Muchiri et al. Reference Muchiri, Ouma and King1996; Hamburger et al. Reference Hamburger, Hoffman, Kariuki, Muchiri, Ouma, Koech, Sturrock and King2004; Kariuki et al. Reference Kariuki, Clennon, Brady, Kitron, Sturrock, Ouma, Ndzovu, Mungai, Hoffman, Hamburger, Pellegrini, Muchiri and King2004; Clennon et al. Reference Clennon, Mungai, Muchiri, King and Kitron2006). 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 {h i} based on the number of parasites they carry.

So each group h i viewed as a local subpopulation carries a parasite range w i<w<w i+1. The dividers {w i} can be chosen arbitrarily, the simplest (extreme) choice being ‘zero-nonzero’ (‘infected’ – ‘non infected’) strata with prevalences {h 0, h 1}. Such partition would result in the standard SIR-type prevalence model. The SIR methodology (h 0; h 1) 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 mathematical 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 w i=iΔw; (i=0, 1, 2, …) counted in steps Δw worms. The choice of Δw is related to infectivity threshold w T – the ‘effective’ number below which a population is considered free of infection. It is convenient to take w T a multiple of Δw; thus w T=Δw means h 0=‘uninfected strata’, w T=2Δw the first two strata (h 0+h 1)=‘uninfected’, etc. Indeed, host infectivity is proportional to the number of mated worm pairs, which in turn depends on their mating probability ϕ(w) (Nasell, Reference Nasell1978; see also Anderson-May, Reference Anderson and May1991). Derivation of such a probability function is based on ‘negative-binomial’ assumption, suggesting a somewhat low threshold value w T=5 (independent of aggregation). On the other hand, experimental/field data for primate infection (Wilson, 2006) predict a substantially higher threshold for egg detection (w T≈40). We assume infectivity threshold to lie somewhere in between 10⩽w T⩽20, and numerically explore 2 values: w T=10 and w T=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 w T has significant effect (see Fig. 6).

For higher burden levels (w>10) we assume that host infectivity, as measured by egg counts, is roughly proportional to w i. 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, S i, 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, Reference Van Der Werf and De Vlas2001)). 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) {h i} move ‘up’ or ‘down’ the scale at certain rates, due to worm acquisition (λ) and resolution (γ i). Each stratum has its source-term S i, and loss rate μ.

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 S 0 (as newborns are free of infection). But the adult source, {S i}, 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=∑0nh i, variables {h i(t)} obey a system of differential equations (SWB)

(1)
$$\openup-3\eqalign{ & {{dh_0 } \over {dt}} = - \left( {\lambda + \mu } \right)h_0 + \gamma _1 h_1 + S_0 ; \cr & ..... \cr & {{dh_i } \over {dt}} = \lambda h_{i - 1} - \left( {\lambda + \gamma _i + \mu } \right)h_i + \gamma _{i + 1} h_{i + 1} + S_i ; \cr & ...... \cr & {{dh_n } \over {dt}} = \lambda h_{n - 1} - \left( {\gamma _n + \mu } \right)h_n + S_n ; \cr} $$

We used the following notations (i) transition rate λ is proportional to (per capita) force of infection: $\lambda = {{Az} \over {{\rm \Delta} w}}$, made of ‘transmission coefficient’/‘infected snail’ – A, infected snail density – z, and Δw – ‘burden increment’ in the definition of strata {h i} (see Appendix A and Table 1); (ii) resolution rates {γ i=γi} are proportional to burden levels {w i} of each strata (γ being natural worm mortality); (iii) host loss rate(s)=μ and new recruit sources={S i}. Here we assume 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 {h i(t)} (1) for several values λ contributed by a stationary source S 0 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 of SWB-system coupled to infected snail prevalence. Left panel shows prevalence distributions of 15 human strata {h i(t)} at t=0 (light grey), t=20 (medium grey) and t=40 (black). For snail prevalence we show the dynamic relaxation pattern (right panel).

Table 1

List of notations, variables and parameters

System (1) can be interpreted as cumulative outcome 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 (Reference Macdonald1965) 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 (Reference Macdonald1965). But, unlike the Macdonald system, our MWB is not an independent variable but computed from prevalence variables {h i}

(2)
$$w = \left( {\rm \Delta} w \right)\sum\limits_1^n {\left( {ih_i /H} \right)} .$$

Rescaling

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,

Then the i th equation for system (1) becomes

$${{dh'_i } \over {dt}} = \lambda h'_{i-1} - \left( {\lambda + \gamma i + \mu } \right)h'_i + \gamma \left( {i + 1} \right)h'_{i + 1} + S_i' ;$$

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 {h i}. Among them are (i) ‘total population dynamics’ for H (0-moment (sum)),

(3)
$${{dH} \over {dt}} = S-\mu H$$

with S=H∑i=0nS i; (ii) Macdonald-type MWB equation for w (1st moment (mean))

(4)
$${{dw} \over {dt}} = \lambda \left( {1 - h_n } \right) - \left( {\gamma + \mu } \right)w + \overline S.$$

Here $\overline S = \sum\nolimits_{i = 0}^n {iS_i } $ designates the cumulative source. Let us note that $\overline 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 h n 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 {h i} 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 Reference Pugliese2000) 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=∑1nih i. Our calibration scheme exploits 2 epidemiological data values, uninfected human prevalence P=h 0, and MWB, w (which is related to measured egg-counts in the field). Both are functions of parameters λ, μ, and source S

(5)
$$\eqalign{ & P = \varphi \left( {\lambda, \mu, S} \right) \vskip6pt\cr & w = \sum\nolimits_1^n {ih_i } = \psi \left( {\lambda ,\mu ,S} \right) .\cr} $$

Both functions ϕ 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=ϕ(λ, μ, S) on worm acquisition rate λ for 3 values μ={0·1γ, 0·2γ, 0·3γ} (top to bottom) and identical stationary source (S) at the lowest strata; (B) Plot of MWB-function w=ψ(λ, μ, S) on λ for the same values of μ as (A) (top to bottom) with dashed lines corresponding to a would-be ‘linear relation’ (slope) ${{\partial w} \over {\partial \lambda }} \approx {1 \over {\gamma + \mu }}$ at λ=0 (appropriate for Macdonald MWB).

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),

(6)
$$\lambda \approx \left( {\gamma + \mu } \right)w$$

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.)

MODEL CALIBRATION

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. Reference Riley, Carabin, Belisle, Joseph, Tallo, Balolong, Willingham, Fernandez, Gonzales, Olveda and Mcgarvey2008).

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. (i) Estimate human force of infection λ (for suitable populations) from demographic and human prevalence data. This part utilizes prevalence function ϕ(λ, μ, S)

  2. (ii) Next estimate snail-to-human transmission coefficient {A} for appropriate population groups, as well as biological parameter ρ (egg-production/worm).

  3. (iii) 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)

(7)
$$\eqalign{ & {{dw} \over {dt}} = \lambda {\rm \Delta} w - \left( {\gamma + \mu } \right)w; \cr & {{dy} \over {dt}} = {\rm \Lambda} \left( {1 - y} \right) - \nu y; \cr} $$

with per capita forces of infection

(8)
$$\eqalign{ & \lambda {\rm \Delta} w = ANy \cr & {\rm \Lambda} = BHw .\cr} $$

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)

$$ {R_0 = {{ANBH} \over {\left( {\gamma + \mu } \right)\nu }} \gt 1.} $$

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

$$\rho = {{{\rm Egg \ release}} \over {{\rm worm} \times {\rm day}}}.$$

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

  1. (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. (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

(9)
$$\matrix{ {A = \displaystyle{{\lambda ({\rm \Delta} w)} \over {Ny}};} & {\rho = \displaystyle{{\left( {\gamma + \mu } \right)E} \over {\lambda ({\rm \Delta} w)}};} & {B = {\displaystyle{\rho \over {HE}}{{\nu y} \over {\left( {1 - y} \right)}}}}. \cr} $$

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. Reference de Vlas, Engels, Rabello, Oostburg, Van Lieshout, Polderman, Van Oortmarssen, Habbema and Gryseels1997; Hamburger et al. Reference Hamburger, Hoffman, Kariuki, Muchiri, Ouma, Koech, Sturrock and King2004) 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 {w s} or the corresponding SWBs for each risk-group) and a single snail equation (contributed by all groups)

(10)
$$\eqalign{ & {{dw_s } \over {dt}} = A_s Ny - \left( {\gamma + \mu _s } \right)w_s ; \cr & {{dy} \over {dt}} = \left( {\sum\limits_s^{} {B_s H_s w_s } } \right)\left( {1 - y} \right) - \nu y. \cr} $$

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

(11)
$$A_s = {{\lambda _s ({\rm \Delta} w)} \over {Ny}}.$$

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,

(12)
$$\eqalign{ & A_s = a\eta _s \cr & B_s = b\eta _s \cr} $$

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

(13)
$$B_s = {{A_s } \over {\sum\nolimits_q^{} {\left( {A_q H_q E_q /\rho _q } \right)} }}{{\nu y} \over {\left( {1 - y} \right)}}$$

expressed through group specific egg-production/worm rates

(14)
$$\rho _s = {{(\gamma + \mu _s )E_s } \over {\lambda _s ({\rm \Delta} w)}}$$

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 for variables H c; H a (total child/adult populations) is given by a Leslie-type differential system

(15)
$$\eqalign{\openup6 & {{dH_c } \over {dt}} = - \left( {\mu _c + \tau } \right)H_c + \beta H_a \vskip14pt\cr & {{dH_a } \over {dt}} = \tau H_c - \mu _a H_a . \cr} $$

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, S c=βH a, enters through the uninfected (zero) stratum (newborns), whereas the adult source, S a, 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

(16)
$$ \eqalign{ & {{dw_c } \over {dt}} = A_c Ny - \left( {\gamma + \mu _c + \tau } \right)w_c \cr & {{dw_a } \over {dt}} = A_a Ny - \left( {\gamma + \mu _a } \right)w_a + \tau {{H_c } \over {H_a }}w_c \cr & {{dy} \over {dt}} = (B_c H_c w_c + B_a H_a w_a )(1 - y) - \nu y .\cr} $$

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

(17)
$$ \openup3\eqalign{ & \varphi \left( {\lambda _c ,\mu _c + \tau ,S^c } \right) = P_c \Rightarrow \lambda _c \cr & \varphi \left( {\lambda _a ,\mu _a ,S^a } \right) = P_a \Rightarrow \lambda _a . \cr} $$

Specifically, we first find children's force of infection λ c with known prevalence P c and source S c (in terms of birth rate), then compute the resulting equilibrium-prevalence distribution $\hat h_c\!=\! \left\{ {h_i^c } \right\}_{i = 0}^n $, substitute it to get adult source S a={τh icH c/H a}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 {A c; A a; B c; B a; ρ 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 {A ij; B ij}, if all connections are allowed (Woolhouse, Reference Woolhouse1991; Woolhouse et al. Reference Woolhouse, Watts and Chandiwana1991; Gurarie and King, Reference Gurarie and King2005; Riley et al. Reference Riley, Carabin, Belisle, Joseph, Tallo, Balolong, Willingham, Fernandez, Gonzales, Olveda and Mcgarvey2008). Different reduction schemes can be used to constrain the resulting parameter set. Thus Gurarie and King (Reference Gurarie and King2005) and Gurarie and Seto (Reference Gurarie and Seto2009) 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. Reference Sturrock, Kinyanjui, Thiongo, Tosha, Ouma, King, Koech, Siongok and Mahmoud1990; Muchiri et al. Reference Muchiri, Ouma and King1996; Hamburger et al. Reference Hamburger, Hoffman, Kariuki, Muchiri, Ouma, Koech, Sturrock and King2004; Kariuki et al. Reference Kariuki, Clennon, Brady, Kitron, Sturrock, Ouma, Ndzovu, Mungai, Hoffman, Hamburger, Pellegrini, Muchiri and King2004; Clennon et al. Reference Clennon, Mungai, Muchiri, King and Kitron2006)) 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 A ij; B ij both proportional to absolute contacts per person from village i to water site j, T iη ij with T i 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 (A 0 times ‘absolute contact rate’),

(18)
$$ \eqalign{ & A_{ij} = A_0 T_i \eta _{ij}^{} = A_i \eta _{ij} ;\quad ( {i = 1,2,\ldots,n} ) \cr & B_{ij} = B_{0j} T_i \eta _{ij}^{} = B_j A_{ij} ;\quad ( {\,j = 1,2,\ldots,m} ) .\cr} $$

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 B 0j as site specific. This hypothesis clearly reduces the number of unknown parameters to a manageable set {A i; B j}, 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 {A ic; A ia}; (ii) local transmission coefficients for snail sites {B j}; (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)

(19)
$$ \eqalign{ & A_i ^s = {{\lambda _i ^s ({\rm \Delta} w)} \over {\sum\nolimits_j^{} {\eta _{ij}^{} Z_j } }}; \quad \rho _c \approx {{(\gamma + \mu _i ^c + \tau )E_i ^c } \over {\lambda _i ^c ({\rm \Delta} w)}}; \cr & \rho _a \approx {{(\gamma + \mu _i ^a )E_i ^a } \over {\tau (H_i ^c /H_i ^a )(E_i ^c /\rho _c ) + \lambda _i ^a ({\rm \Delta} w)}}. \cr}$$
Table 3

Calibrated transmission coefficients and egg production rates based on the individual level field data for 2 choices of discretization step (Δw), and infection threshold w 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 {A is} and {ρ s}, we compute coefficients {B j} from the snail equilibrium equations ((Eq.C.4) in Appendix C).

RESULTS

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. Reference Sturrock, Kinyanjui, Thiongo, Tosha, Ouma, King, Koech, Siongok and Mahmoud1990; Muchiri et al. Reference Muchiri, Ouma and King1996; Hamburger et al. Reference Hamburger, Hoffman, Kariuki, Muchiri, Ouma, Koech, Sturrock and King2004; Kariuki et al. Reference Kariuki, Clennon, Brady, Kitron, Sturrock, Ouma, Ndzovu, Mungai, Hoffman, Hamburger, Pellegrini, Muchiri and King2004; Clennon et al. Reference Clennon, Mungai, Muchiri, King and Kitron2006). 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. (i) 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 ϕ(λ, μ, S c) (Fig. 3) and the derived adult prevalence function ϕ(λ, μ, S a) with source S a from the children's group.

  2. (ii) Next we estimate village-specific transmission coefficients A i(i=1, 2, …, 10) for child/adult groups via formulae (19), and forces of infection λ i (output of step (i)). We also compute the ‘best-fit’ estimates for biological parameters ρ c; ρ a (egg-production/worm).

  3. (iii) Finally, we estimate snail transmission coefficients B j(j=1, 2, …, 5) from equation (Eq.C.4) (Appendix C).

Table 3 shows estimated snail to human transmission coefficients {A c, A a} 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 w T=10 (Table 3). The 2nd case gave overall higher estimates of transmission rates {A i}, 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 density’.

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 w T. To examine such uncertainties on calibrated parameters we took 2 threshold values w T=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 by factor 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.

Fig. 6.

Sensitivity analysis of model parameters (snail-to-human transmission rates) {A i} 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’=0·25–0·75 quantiles; ‘whiskers’=10–90) correspond to infectivity thresholds w T=10 and w T=20. Black dots mark calibrated values of the ‘individual level field data’ in each case.

SUMMARY AND CONCLUSIONS

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, Reference Macdonald1965), 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, Reference Anderson and May1991)) 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, R 0. 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 in exploring 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. Reference Gabrielli, Toure, Sellin, Sellin, Ky, Ouedraogo, Yaogho, Wilson, Thompson, Sanou and Fenwick2006).

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 non-linear 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.

ACKNOWLEDGEMENTS

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

FINANCIAL SUPPORT

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.

APPENDICES

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. (i) Snails are minimally mobile and relatively attached to their local habitat, so we ignore hydrologic (or other) transport of snails/larvae.

  2. (ii) There is an approximately uniform snail density throughout each habitat (their local numbers are part of the collected field data).

  3. (iii) 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, Reference Macdonald1965) for single habitat modified to include multiple risk groups (labelled by s), consists of (i) mean worm burden (MWB) and egg counts {w s, E s}; (ii) susceptible- infected-patent snail densities (or prevalences) {x, y, z}; (iii) cercarial and miracidial densities {c, m}:

(Eq.A.1)
$$\openup-8pt\eqalign{ & {\rm (Human) }\left\{ \openup-2pt\matrix{ \displaystyle{{dw_s } \over {dt}} = a\eta _s c - \left( {\gamma + \cdots} \right)w_s ; \hfill \cr\vskip-12pt\cr E_s = \rho w_s \hfill \cr} \right. \cr & {\rm (Snail) }\left\{ \openup14pt\matrix{ \displaystyle{{dy} \over {dt}} = bm(1 - y) - \nu y; \hfill \cr z = \psi y \hfill \cr} \right. \cr & {\rm (Larvae)}\left\{ \matrix{ \displaystyle{{dm} \over {dt}} = \sum\limits_s {Q_s } - \mu _M m; \hfill \cr \displaystyle{{dc} \over {dt}} = \pi _C Nz - \mu _C c. \hfill \cr} \right. \cr} $$

The miracidial production term is given by Q s=π Mη sH sE s – proportional to the population density (per unit habitat).

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

(Eq.A.2)
\left\{ \matrix{ \displaystyle{{dw_c } \over {dt}} = a\eta _c c - \left( {\gamma + \mu _c + \tau } \right)w_c \hfill \cr\cr \displaystyle{{dw_a } \over {dt}} = a\eta _a c + \tau \displaystyle{{H_c } \over {H_a }}w_c - \left( {\gamma + \mu _a } \right)w_a .\hfill \cr} \right.

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, $\tau {{H_c } \over {H_a }}$, in the adult-equation, weighted in proportion to population fraction H c/H a. In this formulation we assumed stationary populations, but the equations can be easily 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

$$\eqalign{ & c^* = {{\pi _C } \over {\mu _C }}Nz \vskip3pt\cr & m^* = \sum\limits_s {{{\pi _M } \over {\mu _M }}\eta _s H_s E_s } . \cr} $$

The resulting human–snail system has MWB-type human equations: $\displaystyle{{dw_s } \over {dt}} = A_s Ny - ( {\gamma + \cdots}) w_s $ (for s-group), or specifying it for a ‘child–adult’ setup we get

(Eq.A.3)
$$\left\{ \openup4pt\matrix{ \displaystyle{{dw_c } \over {dt}} = A_c Ny - \left( {\gamma + \mu _c + \tau } \right)w_c \hfill \cr\cr \displaystyle{{dw_a } \over {dt}} = A_a Ny + \tau {{H_c } \over {H_a }}w_c - \left( {\gamma + \mu _a } \right)w_a . \hfill \cr} \right. $$

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

(Eq.A.4)
$$\openup3\eqalign{ & {{dy} \over {dt}} = {\rm \Lambda} (1 - y) - \nu y; \cr & x = \left( {1 - y - z} \right); \cr & z = \psi y; \cr} $$

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

(Eq.A.5)
$${\rm \Lambda} = \sum\limits_s {B_s H_s w_s } .$$

The transmission coefficients {A s; B s} (‘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

(Eq.B.1)
$$\eqalign{ & {\rm (Child)}\left\{ \matrix{ \displaystyle{{dh_0^c } \over {dt}} = - \left( {\lambda _c + \mu _c + \tau } \right)h_0^c + \gamma _1 h_1^c + \beta \displaystyle{{H_a } \over {H_c }}; \hfill \vskip9pt\cr .....\hfill\vskip-13pt\cr \hfill \cr \displaystyle{{dh_i^c } \over {dt}} \!=\! \lambda _c h_{i - 1}^c\! -\! \left( {\lambda _c + \gamma _i + \mu _c + \tau } \right)h_i^c + \gamma _{i + 1} h_{i + 1}^c ; \hfill\vskip9pt \cr .....\hfill\vskip-13pt\cr \hfill\vskip-6pt \cr \displaystyle{{dh_n^c } \over {dt}} = \lambda _c h_{n - 1}^c - \left( {\gamma _n + \mu _c + \tau } \right)h_n^c ; \hfill \cr} \right. \eqalign{ {\rm (Adult)}\left\{ \openup6\matrix{ \displaystyle{{dh_0^a } \over {dt}} = - \left( {\lambda _a + \mu _a } \right)h_0^a + \gamma _1 h_1^a + \tau \displaystyle{{H_c } \over {H_a }}h_0^c ; \hfill \cr ..... \hfill\cr\vskip-12pt\hfill\cr \displaystyle{{dh_i^a } \over {dt}} =\lambda _a h_{i - 1}^c\! -\! \left( {\lambda _a + \gamma _i + \mu _a } \right)h_i^c + \gamma _{i + 1} h_{i + 1}^c\cr\tab\quad\hskip18pt + \tau \displaystyle{{H_c } \over {H_a }}h_i^c ; \hfill \cr ...... \hfill\vskip-4pt\cr\hfill\vskip-4pt \cr \displaystyle{{dh_n^a } \over {dt}} = \lambda _a h_{n - 1}^a - \left( {\gamma _n + \mu _a } \right)h_n^a + \tau \displaystyle{{H_c } \over {H_a }}h_n^c. \hfill \cr} \right. \cr} $$

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

$${{H_c } \over {H_a }} = {\beta \over {\mu _c + \tau }} = {{\mu _a } \over \tau }.$$

We estimated τ and β from the demographic data of the Msambweni region (Sturrock et al., Reference Sturrock, Kinyanjui, Thiongo, Tosha, Ouma, King, Koech, Siongok and Mahmoud1990; Muchiri et al., Reference Muchiri, Ouma and King1996; Hamburger et al., Reference Hamburger, Hoffman, Kariuki, Muchiri, Ouma, Koech, Sturrock and King2004; Kariuki et al., Reference Kariuki, Clennon, Brady, Kitron, Sturrock, Ouma, Ndzovu, Mungai, Hoffman, Hamburger, Pellegrini, Muchiri and King2004; Clennon et al., Reference Clennon, Mungai, Muchiri, King and Kitron2006). The child–adult death rates were computed from

$$\mu _c = \beta {{H_a } \over {H_c }} - \tau ; \quad \mu _a = \tau {{H_c } \over {H_a }}.$$

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)

(Eq.B.2)
$$\openup6pt\eqalign{ \openup4pt& A_c = {{\gamma \lambda _c {\rm \Delta} w} \over {Ny}}; \quad A_a = {{\gamma \lambda _a {\rm \Delta} w} \over {Ny}}; \cr & \rho _c = {{\left( {\gamma + \mu _c + \tau } \right)E_c } \over {\gamma \lambda _c ({\rm \Delta} w)}}; \quad \rho _a = {{\left( {\gamma + \mu _a } \right)E_a } \over {\gamma \lambda _c ({\rm \Delta} w) + \tau {{H_c E_c } \over {H_a \rho _c }}}}; \cr & B_c = {{\nu yA_c } \over {\left( {1 - y} \right)\sum\nolimits_s^{} {A_s H_s E_s /\rho _s } }}; \cr & B_a = {{\nu yA_a } \over {\left( {1 - y} \right)\sum\nolimits_s^{} {A_s H_s E_s /\rho _s } }}. \cr} $$

C. Distributed system

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

(Eq.C.1)
$$\left\{ \matrix{ \displaystyle{{dw_i } \over {dt}} = \sum\limits_{\,j = 1}^m {A_{ij}^{} N_j y_j } - (\gamma + \mu _c )w_i^{} ; \quad ( {i = 1,\ldots,n} ) \hfill \cr\hfill\cr \displaystyle{{dy_j } \over {dt}} = \sum\limits_{\,j = 1}^m {B_{ij}^{} H_i w_i (1 - y_j )} - \nu y_j ; \quad ( {\,j = 1,\ldots,m} ) .\hfill \cr} \right. $$

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

(Eq.C.2)
$$\eqalign{ & ({\rm Human})\left\{ \matrix{ \displaystyle{{dw_i^c } \over {dt}} = A_i^c {\rm }\sum\limits_{\,j = 1}^m {\eta _{ij} Z_j } - (\gamma + \mu _c + \tau )w_i^c ; \hfill \cr \displaystyle{{dw_i^a } \over {dt}}\! =\! A_i^a {\rm }\sum\limits_{\,j = 1}^m {\eta _{ij} Z_j } + \tau {{H_i^c } \over {H_i^a }}w_i^c - (\gamma + \mu _a )w_i^a ;\cr\qquad\quad ( {i = 1,\ldots,n} ) \hfill \cr }\right. \cr & ({\rm Worm})\left\{ \openup3\matrix{ E_i^c = \rho _{}^c w_i^c ; \hfill \cr E_i^a = \rho _{}^a w_i^a ; ( {i = 1,\ldots,n} ) \hfill \cr} \right. \cr & ({\rm Snail})\left\{ \openup3\matrix{ \displaystyle{{dY_j } \over {dt}} = {\rm \Lambda} _j \left( {N_j - Y_j } \right)- \nu Y_j ; \vskip6pt\hfill \vskip6pt\cr Z_j = \psi _j Y_j ;\left( {\,j = 1,\ldots\!,m} \right) .\hfill \cr} } \right.} \cr} $$

Here we use local patent snail densities {Z j} instead of infected snail prevalences of earlier special cases (formulae (9)–(11)–(13)). Since patency conversion rates of infected snails {ψ j=Z j/Y j} 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).

(Eq.C.3)
$${\rm \Lambda}_j = B_j \sum\limits_i {\eta _{ij} (A_i^c H_i^c E_i^c + A_i^a H_i^a E_i^a )} /U_j .$$

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

$$\eqalign{ & A_i^s = a\left( {{{\pi _C } \over {\nu _C }}} \right)T_i^s; \cr & B_j = \displaystyle{{b_j \left( {{{\pi _M } \over {\nu _M }}} \right)} \bigg/ {a\left( {{{\pi _C } \over {\nu _C }}} \right)}}, \cr} $$

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

(Eq.C.4)
$$B_j = {{\nu Y_j } \over {(N_j - Y_j )\sum\limits_i {\eta _{ij} (A_i^c H_i^c E_i^c + A_i^a H_i^a E_i^a )/U_j } }} .$$

References

REFERENCES

Alexander, N., Moyeed, R. and Stander, J. (2000). Spatial modelling of individual-level parasite counts using the negative binomial distribution. Biostatistics 1, 453–463.Google Scholar
Anderson, R. M. and May, R. M. (1978 a). Regulation and stability of host-parasite population interactions. I. Journal of Animal Ecology 47, 219–247.Google Scholar
Anderson, R. M. and May, R. M. (1978 b). Regulation and stability of host-parasite population interactions. II. Journal of Animal Ecology 47, 249–267.CrossRefGoogle Scholar
Anderson, R. M. and May, R. M. (1991). Infectious Diseases of Humans. Dynamics and Control. Oxford University Press, New York, USA.CrossRefGoogle Scholar
Chan, M. S. and Bundy, D. A. (1997). 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 91, 216–220.Google Scholar
Chan, M. S., Guyatt, H. L., Bundy, D. A. and Medley, G. F. (1996). Dynamic models of schistosomiasis morbidity. American Journal of Tropical Medicine and Hygiene 55, 52–62.Google Scholar
Clennon, J. A., Mungai, P. L., Muchiri, E. M., King, C. H. and Kitron, U. (2006). Spatial and temporal variations in local transmission of Schistosoma haematobium in Msambweni, Kenya. American Journal of Tropical Medicine and Hygiene 75, 1034–1041.Google Scholar
Dobson, A. and Roberts, M. (1994). The population dynamics of parasitic helminth communities. Parasitology 109, (Suppl.) S97–S108.CrossRefGoogle ScholarPubMed
Feng, Z., Curtis, J. and Minchella, D. J. (2001). The influence of drug treatment on the maintenance of schistosome genetic diversity. Journal of Mathematical Biology 43, 52–68.CrossRefGoogle ScholarPubMed
Gabrielli, A. F., Toure, S., Sellin, B., Sellin, E., Ky, C., Ouedraogo, H., Yaogho, M., Wilson, M. D., Thompson, H., Sanou, S. and Fenwick, A. (2006). 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 99, 234–242.Google Scholar
Gurarie, D. and King, C. H. (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.Google Scholar
Gurarie, D. and King, C. H. (2008). Age- and risk-targeted control of schistosomiasis-associated morbidity among children and adult age groups. The Open Tropical Medicine Journal 1, 21–30.Google Scholar
Gurarie, D. and Seto, E. Y. (2009). Connectivity sustains disease transmission in environments with low potential for endemicity: modelling schistosomiasis with hydrologic and social connectivities. Journal of the Royal Society Interface 6, 495–508. doi: 10.1098/rsif.2008.0265.Google Scholar
Hamburger, J., Hoffman, O., Kariuki, H. C., Muchiri, E. M., Ouma, J. H., Koech, D. K., Sturrock, R. F. and King, C. H. (2004). 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 71, 765–773.Google Scholar
Kariuki, H. C., Clennon, J. A., Brady, M. S., Kitron, U., Sturrock, R. F., Ouma, J. H., Ndzovu, S. T., Mungai, P., Hoffman, O., Hamburger, J., Pellegrini, C., Muchiri, E. M. and King, C. H. (2004). 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 70, 449–456.Google Scholar
Macdonald, G. (1965). The dynamics of helminth infections, with special reference to schistosomes. Transactions of the Royal Society of Tropical Medicine and Hygiene 59, 489–506.Google Scholar
Medley, G. F. and Bundy, D. A. (1996). Dynamic modeling of epidemiologic patterns of schistosomiasis morbidity. American Journal of Tropical Medicine and Hygiene 55, 149–158.CrossRefGoogle ScholarPubMed
Muchiri, E. M., Ouma, J. H. and King, C. H. (1996). Dynamics and control of Schistosoma haematobium transmission in Kenya: an overview of the Msambweni Project. American Journal of Tropical Medicine and Hygiene 55, 127–134.Google Scholar
Nasell, I. (1978). Mating for schistosomes. Journal of Mathematical Biology 6, 21–35.CrossRefGoogle ScholarPubMed
Pugliese, A. (2000). Coexistence of macroparasites without direct interactions. Theoretical Population Biology 57, 145–165. doi: S0040-5809(99)91443-0 [pii].CrossRefGoogle ScholarPubMed
Riley, S., Carabin, H., Belisle, P., Joseph, L., Tallo, V., Balolong, E., Willingham, A. L., Fernandez, T. J., Gonzales, R. O., Olveda, R. and Mcgarvey, S. T. (2008). Multi-host transmission dynamics of Schistosoma japonicum in Samar province, the Philippines. PLoS Med 5, e18. doi: 10.1371/journal.pmed.0050018.Google Scholar
Rosa, R. and Pugliese, A. (2002). Aggregation, stability, and oscillations in different models for host-macroparasite interactions. Theoretical Population Biology 61, 319–334. doi: 10.1006/tpbi.2002.1575CrossRefGoogle ScholarPubMed
Sturrock, R. F., Kinyanjui, H., Thiongo, F. W., Tosha, S., Ouma, J. H., King, C. H., Koech, D., Siongok, T. K. and Mahmoud, A. A. (1990). 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 84, 257–261.CrossRefGoogle ScholarPubMed
Van Der Werf, M. J. and De Vlas, S. J. (2001). Morbidity and Infection with Schistosomes or Soil-Transmitted Helminths. Report for WHO Parasitic Diseases and Vector Contol. pp. 1–103. Erasmus University, Rotterdam, The Netherlands.Google Scholar
de Vlas, S. J., Engels, D., Rabello, A. L., Oostburg, B. F., Van Lieshout, L., Polderman, A. M., Van Oortmarssen, G. J., Habbema, J. D. and Gryseels, B. (1997). Validation of a chart to estimate true Schistosoma mansoni prevalences from simple egg counts. Parasitology 114, 113–121.CrossRefGoogle ScholarPubMed
Wilson, R. A., Dam, G. J., Kariuki, T. M., Farah, I. O., Deedler, A. M. and Coulson, P. S. (2006). The detection limits for estimates of infection intensity in schistosomiasis mansoni established by a study in non-human primates. International Journal for Parasitology 36, 1241–1244.Google Scholar
World Health Organization (2006). 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.Google Scholar
Woolhouse, M. E. (1991). On the application of mathematical models of schistosome transmission dynamics. I. Natural transmission. Acta Tropica 49, 241–270.Google Scholar
Woolhouse, M. E., Watts, C. H. and Chandiwana, S. K. (1991). Heterogeneities in transmission rates and the epidemiology of schistosome infection. Proceedings of the Royal Society of London, B 245, 109–114.Google Scholar
Figure 0

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 stratum has its source-term Si, and loss rate μ.

Figure 1

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 of SWB-system coupled to infected snail prevalence. Left panel shows prevalence distributions of 15 human strata {hi(t)} at t=0 (light grey), t=20 (medium grey) and t=40 (black). For snail prevalence we show the dynamic relaxation pattern (right panel).

Figure 2

Table 1 List of notations, variables and parameters

Figure 3

Fig. 3. (A) Plot of uninfected prevalence function P=ϕ(λ, μ, S) on worm acquisition rate λ for 3 values μ={0·1γ, 0·2γ, 0·3γ} (top to bottom) and identical stationary source (S) at the lowest strata; (B) Plot of MWB-function w=ψ(λ, μ, S) on λ for the same values of μ as (A) (top to bottom) with dashed lines corresponding to a would-be ‘linear relation’ (slope) ${{\partial w} \over {\partial \lambda }} \approx {1 \over {\gamma + \mu }}$ at λ=0 (appropriate for Macdonald MWB).

Figure 4

Table 2 Human and snail demographics and infection data: ‘egg-count’ refers to arithmetic mean/host

(Prevalence, infected fraction of population.)
Figure 5

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

Figure 6

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.

Figure 7

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 density’.

Figure 8

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’=0·25–0·75 quantiles; ‘whiskers’=10–90) correspond to infectivity thresholds wT=10 and wT=20. Black dots mark calibrated values of the ‘individual level field data’ in each case.