Description of LYMFASIM

LYMFASIM is a stochastic microsimulation model for the epidemiology of lymphatic filariasis in human populations (Habbema et al. 1996; Plaisier et al. 1998). The model simulates the life-histories of human individuals (birth, acquisition and loss of parasites, death) and individual parasites (maturation, mating, production of Mf, death). Together, the simulated persons constitute the population of an endemic village or area. A detailed description and mathematical formulation of the model has been given in an earlier publication (Plaisier et al. 1998). Here we restrict ourselves to a brief description of the model and the factors that are directly relevant to the effects of vector control. Of particular importance are the regulation of parasite density in both the vector and the human host.

Transmission and parasite dynamics

A graphical representation of the model is given in Fig. 1. In this figure, the monthly force-of-infection (foi) indicates the number of parasites that enter the human host in a month and the proportion that develops successfully into adult worms, sr. The force-of-infection varies between individuals and over time; its calculation is explained below.

Fig. 1. Schematic representation of LYMFASIM. Immune-regulation is optional in LYMFASIM. The shaded boxes are entomological variables; these variables do not vary between individuals. The unshaded boxes are human variables; the index i indicates that the variable may differ between individuals. Along the arrows, the symbols of relevant parameters are shown.

In the case of a constant force-of-infection, the expected equilibrium worm-load (M, number of mature worms) in a person is found by multiplying the force-of-infection for this person times the average reproductive life-span (i.e. total life-span minus duration of immature stage) of an adult parasite. The total life-span of the worm is assumed to vary between parasites, and is described by a Weibull distribution with mean T_l and shape-parameter α_Tl. Estimates for the life-span of Onchocerca volvulus (the filarial nematode species causing human onchocerciasis), indicated less than exponential variation (α_Tl>1) and hence we have fixed the value of α_Tl to 2·0 (Plaisier et al. 1991). The duration Ti of the immature stage is considered to be 8 months (World Health Organization, 1992). The parasite life-span not only determines the equilibrium worm-load, but also the rate of worm-mortality and thereby the rate at which the worm-load declines in the case of interruption of transmission.

Female adult worms produce microfilariae (Mf), provided that the human host harbours at least one adult male parasite, assuming a totally polygamous system in W. bancrofti. The Mf production is equal to r₀ per month per 20 μl blood per female parasite in the absence of an anti-fecundity immune response, but is reduced when the human hosts develop such a response (see below). The simulated true Mf-density, m, in a person is expressed in terms of the average number per 20 μl peripheral blood taken for diagnosis and is updated monthly using the number of Mf produced by each female worm per month. Mf-mortality is governed by a monthly survival fraction s=0·9 for the Mf (Plaisier et al. 1999). The variability (between blood samples within a host) in the actual (discrete) number of Mf counted in the smear is described by a negative binomial distribution with a mean equal to the true Mf-density in an individual and a parameter of dispersion k_m. Overdispersion will be smaller (k_m larger) when a larger volume of blood is examined (due to increased sensitivity). Due to intra-host and observer variability in Mf counts, false Mf-negative cases (count=0) can occur.

Based on experimental data, the relation between the Mf-density m in a human and the number of L3 that will develop in Culex quinquefasciatus mosquitoes, the principal vector of W. bancrofti in Pondicherry, feeding on such a person (L3 from bloodmeal) is described by the following hyperbolic function (Subramanian et al. 1998),

with parameter values in Table 1.

This relationship saturates at φ/ζ at high human Mf densities and has an initial slope of φ. Because of this saturation, the development of the parasite in the vector is one of the density regulation mechanisms in the transmission of the parasite.

The number of L3-stage larvae released per mosquito bite (L3) depends on this relationship between Mf-density in human and L3 developing in a mosquito, and also on a number of mosquito-related factors, such as the survival of the mosquitoes between the uptake of microfilariae and the development to the L3-stage under natural conditions, the fraction of mosquitoes that is potentially infectious (i.e. taking into account that some mosquitoes never had a bloodmeal before), and the probability that a L3-larva will actually be released during the act of biting. These mosquito-related factors have been combined in the factor v. Since v and sr are linear multiplication factors in the same sequence of calculations, we decided to arbitrarily fix the proportion v at 0·1, and estimate sr. The average number of infective larvae L3 released per mosquito-bite is calculated as a population average, by weighting each person with his or her relative exposure.

An individual's relative exposure to bites depends on his/her age, but there is also inter-individual variability. We assume the following relation between age and exposure: at birth a person has a relative exposure of E₀, and thereafter it increases linearly until age a_max at which a maximum exposure is reached, which remains at this level for the remainder of life. The variation in mosquito bites between individuals is captured by a personal ‘exposure index’. This exposure index is assumed to be a life-long characteristic of a person. Its value is randomly selected from a gamma probability distribution with mean=1 and shape-parameter α_E. This gamma-distribution allows for persons to have low or very low relative exposure, but it does not allow for zero exposure. We therefore consider an additional parameter, the fraction f₀ of persons that is never exposed to the bites of mosquitoes. We assume that males and females are equally exposed to mosquito bites.

The monthly transmission potential (mtp) is defined as the number of incoming L3 larvae per person per month, which varies between individuals and over time. The transmission potential is calculated as the product of the average monthly biting rate (mbr, number of mosquito-bites per month for an adult person), the relative exposure to bites of this person, and the average number of infective L3 released per mosquito-bite into a human host. Only a fraction of the larvae that entered the human body will survive the larval stages and develop into mature adult worms. This brings us back to the monthly force-of-infection, which depends on the monthly transmission potential, on the proportion (sr) of inoculated larvae that will survive the L3 and L4 stages, and on the individual's level of immunity to L3-larvae, which may vary between 0 (no immunity) and 1 (full immunity, no larva will survive).

Immune-regulation of parasite numbers

In LYMFASIM we assume two mechanisms for the working of the immune system on the dynamics of the parasite: anti-L3 immunity and anti-fecundity immunity. Anti-L3 immunity is triggered by exposure to L3-antigens and reduces the success of inoculated L3-larvae to mature in the human body. This mechanism is proposed on the basis of work by Day et al. (1991a,b) who found, among people followed for one year, an increase in antibodies to the L3 surface mainly in subjects aged 20 years and older i.e. subjects with the longest history of L3-inoculation. Beuria et al. (1995) also found an age-specific increase in the presence of antibodies and further concluded that antibody levels were highly variable between individuals. Further, a recent study showed that the prevalence of antibodies to L3 surface antigens was higher among amicrofilaraemic persons with or without antigenaemia than in subjects with microfilaraemia (Helmy et al. 2000). Several other epidemiological studies also provide indirect evidence for the possible role of acquired immunity in regulating filarial infections (Vanamail et al. 1989; Das et al. 1990; Bundy & Medley, 1992; Michael & Bundy, 1998, Michael et al. 2001a). However, the above field observations (Beuria et al. 1995; Day et al. 1991a,b) corroborate the evidence from laboratory studies on cat-Brugia pahangi (Denham et al. 1972, 1983; Grenfell, Michael & Denham, 1991; Michael et al. 1998; Devaney & Osborne, 2000) and jird-Acanthocheilonema viteae (Eisenbeiss, Apfel & Meyer, 1994; Bleiss et al. 2002) models that immunity acts against re-infection.

Anti-fecundity immunity reflects that prolonged presence of adult parasites may cause a breakdown in tolerance to the parasites, resulting in clearance of microfilariae and progress of disease (Maizels & Lawrence, 1991). Whether and to what extent the adult worms or the microfilariae are the target of this response is not yet clear. In the model we assume that the immune response causes a reduction in Mf production.

The modelling of these two types of immunity is similar (see Fig. 1), and is analogous to Woolhouse's (1992) ‘larval antigens, anti-larval response (LL)’ and ‘adult worm antigens, anti-egg response (AE)’ models. In the following, those parameters referring to the anti L3-immunity and the anti-fecundity models are denoted by, respectively, attaching a suffix l or w to the corresponding symbols. Cumulative ‘experience’ (H) of, respectively, L3 infection (H_l) and adult worm infection (H_w) determines the level of immunity, R_l and R_w. Loss of experience is governed by THl and THw, the half-life (in years) of experience of infection in the absence of boosting. The factor γ (‘strength of immunity’) translates the experience of infection into an immune response (γ_l and γ_w). The immune responsiveness levels R_l and R_w vary between individuals according to a gamma-probability distribution with mean 1 and shape-parameters α_l and α_w. A list and definitions of the model variables and parameter values for which we use external sources (observations, experiments, and literature), or for which we simply fixed the value within plausible ranges is given in Table 1. Table 2 summarizes the parameters estimated from fitting the models to the Pondicherry data.

Table 2. Parameters of LYMFASIM describing the transmission dynamics of Wuchereria bancrofti in humans and their estimated values arising from the fit of models with and without immunity (Units are in years unless otherwise stated. The sign ‘—’ denotes parameters that are not included in a particular model. Values in parentheses are 95% confidence boundaries for the duration of the immunological memory, success ratio and the estimates for the strength of the immune-response corresponding to lower and upper boundaries of the duration of immunological memory.)

Model quantification

The population of Pondicherry in 1981 is simulated by quantifying the life-table and human fertility from statistics for that year (Registrar General of India & Census Commissioner, 1981). The values for the monthly biting rate (mbr, see Fig. 2) during the vector management programme were estimated from fortnightly collection of human landing mosquitoes in one site in Pondicherry (Ramaiah et al. 1992). The mbr was used to assess the seasonal effect on vector population and also to monitor the impact of integrated vector management. Entomological observations indicated that the vector management programme has achieved a large reduction in transmission but did not achieve a total interruption (Ramaiah et al. 1992): within 2 years the annual infective biting rate was reduced by 86% and in 4 years by 94%; the average annual infective biting rate during the programme period was 45, compared to 228 prior to its start (80% reduction). Assuming that the observed pre-control monthly biting rate is representative for the situation in Pondicherry prior to the year 1981, we fixed the monthly biting rate at 2200 per adult person for the period before the start of vector management and after its cessation.

Fig. 2. Observed monthly biting rate in Pondicherry over the period 1980–1986.

Simulations are always started 150 years before 1981 in order to ensure an equilibrium age-composition of the human population and a dynamic equilibrium for the parasite population. The two types of immune response are considered in separate models. The parameters for the anti-L3 immunity are estimated by assuming that there is no anti-fecundity immunity, and vice versa. Adding the possibility of no immune-regulation, we have three versions of the full LYMFASIM model: anti-L3 immunity model, anti-fecundity immunity model, and no-immunity model.

Data

Epidemiological data are from the five-year Integrated Vector Management programme in Pondicherry. Surveys were carried out right before and after the completion of the programme (in 1981 and 1986). Details of sampling design and parasitological data collection are given by Rajagopalan et al. (1989) and Subramanian et al. (1989b). Mf-counts in 20 μl blood smears for both 1981 and 1986 are available for a cohort of 5071 persons. To enable a comparison of simulation results with the observations, the longitudinal data are represented as age-specific cross-tabulations of the Mf-count in 1981 versus the Mf-count in 1986 (Table 3). Data on overall Mf prevalence in 1989 and 1992 (Manoharan et al. 1997) are used for a first validation of the model.

Goodness-of-fit

Simulation results from the three models are compared with data for each of the cells in Table 3. The agreement between observed and simulated data is assessed by the following statistic,

with: O_aij: Observed no. of persons in age-class a (3–7, 8–10, etc.) of whom the Mf-count in 1981 was in class i (0, 1–5, or >5 Mf per smear) and the Mf-count in 1986 was in class j. E_aij: See O_aij, for the simulated population. C_a: O_a/E_a, with O_a total observed and E_a total simulated no. of persons in age-class a.

In some age-classes, cells with i and j combinations have been merged to ensure that they comprise at least 5 observed individuals. The factor (1+C_a) in the denominator accounts for the stochastic variation in the simulated population (i.e. the ‘expected’ number is derived from a finite simulated population; with increasing simulation size, C_a approaches zero).

A P-value for the goodness-of-fit is calculated by assuming that X² follows a χ²-distribution with D.F.=42 for models with anti-L3 or anti-fecundity immunity, and D.F.=44 for the model without immunity. The number of degrees of freedom is derived from the number of cells in Table 3 (72), minus the number of cells combined with other cells to ensure a minimum of 5 persons in each (combined) cell (12), minus the number of age-groups (8), minus the number of parameters to be estimated on the basis of the data (10 for the immunity models and 8 for the model without immunity). P-values >0·05 are taken to indicate a satisfactory agreement between estimations and observed data.

Due to the stochastic nature of the various processes involved in the model, the simulation output will be subject to random variation and will only represent an estimate of the true outcomes of the model. As a compromise between random variation and computing time for each version of the model (no immunity, anti-L3 or anti-fecundity immunity) a maximum of 1500 simulation runs was carried out and the total number of individuals per simulation run is approximately 50000.

As a result of variability in simulation output the standard estimation procedures (e.g. maximum likelihood estimation) are not applicable. Instead parameters in Table 2 are estimated by minimizing X² in Eqn. (2) through a downhill-simplex routine (Nelder & Mead, 1965). For the immunity models, a 95% CI was determined for the immunological memory (parameters THl and THw) and for the success ratio (parameter sr) following the method of Plaisier et al. (1995). Starting from the best-fitting values of these parameters, alternative lower and higher values are tested and the other parameters re-estimated. Those values that result in a X²-difference of approximately 3·84 (95th percentile of a χ²-distribution with D.F.=1) are considered to be the CI-boundaries.