Basic reproduction matrix and endemic equilibria of the extended model

By analogy with the single-site, two dimensional case (1.1), equilibria of system (1.8) depend on a matrix analog the basic reproduction number R₀ (1.2), namely, an L×L (L – number of snail sites) basic reproduction matrix,

Its j, k entry measures a contribution of k-th snail site to schistosome reproduction at the j-th site, through distributed human intermediaries. We compute equilibria of (1.8), by first solving the ŵ-equations and then substituting it in the ŷ-equation. The offshoot is a nonlinear system for vector variable ŷ, and a linear operator relating ŵ to ŷ

The introduction of regular chemotherapy will change the worm attrition matrix Γ in (1.11)–(1.12) to Γ+Θ, and thus make R, ŷ, ŵ and the resulting community mean burden,

functions of the treatment frequencies

.

System (1.8) is known to have a nonzero stable (endemic) equilibrium (1.12), provided the largest eigenvalue λ₁(R)>1. So λ₁ plays the role of the ‘basic reproduction number’ for equation system (1.8).

To proceed further we make another simplifying assumption.

Independence hypothesis: we assume the basic environmental parameters are independent of age/behavioural factors. Mathematically, it would turn all essential variables and parameters into products of age-dependent and environmental terms. Specifically, we assume:

1. identical age structure at each locale, with population fractions
   
   (age distribution), and
   
   (site distribution). So the population of a-th age stratum at i-th site is H_i,a=Hh_a*h_i,
2. a universal water contact frequency pattern defined for various age groups {ω_a} (independent of location), multiplied by a hurdle factor {χ_ij} to account for the geographical distance between the i (human) and j (snail) sites. In other words, each age stratum at site i redistributes its maximal contact rate ω_a among different snail sites {j}, according to hurdle factors
   
   .

Thus our model environment has ‘product-type’ populations and contact patterns

For our numerical example (discussed below), we assume hurdle factors {χ_ij} depend on geographical distances {d_ij} between human and snail sites via an ‘inverse square law’ (Haynes & Fotheringham, 1984). Let us note, however, that any other realistic hurdle pattern can easily be accommodated in our scheme.

While the demographic component (a) of (1.14) seems reasonable based on observed population structure in rural endemic households and villages (King et al. 2004), the validity of contact rate estimates (b) remains open (Muchiri, 1991). Further work, based on the available field data will test its (approximate) validity, find proper form of geographical hurdles, and possibly better spatial and socioeconomic models, in order to pattern a more realistic transmission environment (Woolhouse et al. 1998).

Under the independence hypothesis (1.14), the reproduction matrix R (1.11), the resulting worm burden distribution ŵ, and the community mean burden

, will all factor into products of environmental and age-dependent terms. Precisely,

is a product of the basic reproduction scalar (B.R.S.)

and a simplified ‘environmental’ matrix, R⁰, with entries

The former, ρ, encodes age-behavioural data, the latter accounts for the environmental factors (hurdles and spatial distribution of human and snail populations).

The resulting equilibrium snail infection prevalences and human worm burden (1.12) will depend on parameter ρ, through the infection potential. Namely, the (i,a) human stratum has burden

where û=(u_a) designates a universal age pattern of burden distribution

while coefficient

represents (environmental) infection potential at site i due to snail prevalence distribution (at various sites j=1,…, L) via interconnecting hurdle factors. Hence, we find the community mean burden over the entire region

factors into a product of age-dependent terms

and the global infection potential

Rescaled model and its parameters

While certain variables and parameters used in system (1.8) are readily available from observational data (el Kholy et al. 1989; Clennon et al. 2004; Kariuki et al. 2004) and/or experiment, some others (e.g. probabilities of worm establishment {α_a} or contamination rates {β_a}) are uncertain, and can be estimated only by proxy data or through predicted outcomes. In general, given a multi-parameter system like (1.8), one tries to simplify (non-dimensionalize) it by rescaling variables and identifying the essential parameters.

We shall do it for system (1.8) assuming independence hypothesis. To this end we rescale all vector/matrix quantities by their ‘typical/maximal’ value, designated by ‘bar’ and the corresponding normalized dimensionless quantities (patterns) designated by ‘star’,

The analysis reveals two essential parameters of (1.8), namely a typical B.R.S. and typical worm burden,

We also introduce the dimensionless age pattern (of worm distribution)

and rescaled infection potential

All other parameters and equilibria of (1.8) are expressed through (1.24)–(1.26), namely,

and equilibrium burden distributions w_i,a=w₀f_i*(ρ)u_a* (by ‘age’ and ‘site’), hence

all w₀ – multiples (1.24) of infection potential functions of ρ.

The age pattern û* represents a universal worm distribution by age observable in any locale (or over the region), under the independence hypothesis. Fig. 3 shows û* for the 12 age-bin strata. Here we used contact rates {ω_a} of Chan et al. (1997), (tabulated in Table 1), a typical worm attrition

, and exponentially decreasing population fractions

. Qualitatively our worm burden pattern û* resembles the observed w-distributions (see, for example, Butterworth et al. (1996)). In endemic communities they typically peak between ages 10 and 20, then fall off (for older ages) by factor 10 or more in adult life (Butterworth et al. 1996). Such patterns cannot be explained by changing contact rates alone (for older groups). The additional protection, represented here by α*, could be acquired either through natural aging process (e.g. skin thickening or changing patterns of water exposure), or represent a function of the life-long worm exposure or egg-burden experience (acquired immunity). Here we shall consider only the former case (age-protection). The analysis of acquired immunity in a heterogeneous environment is more involved and we shall pursue it in a separate study. We demonstrate the effect of age-protection on burden patterns ŵ or û, by comparing the ‘unprotected’ distribution vs an ‘age-protected’ one (grey vs black on the right panel of Fig. 3). As expected, factors {0<α_a^*[les ]1} drive down the worm burden distribution for older groups. More significantly, the ratio of maximal-to-minimal burdens takes on a higher (more realistic) value in the latter case.

Fig. 3. The effect of an age-related protection factor {αa*} (left panel) that limits the rate of new infection on the mean worm burden age-distribution {wa} for fixed population fractions {ha*}, and contact rates {ωa*} of Table 1. In the right panel, grey bars indicate age-specific mean worm burden in the setting without age-related protection factors vs the setting where age-related protection is included (black bars).

Equations (1.20)–(1.23), or rescaled ones (1.25)–(1.27), yield a manageable formulation for analysis and therapy control. Indeed, control parameters {τ} will modify (augment) the worm attrition matrix Γ→Γ+Θ, and thus enter explicitly into scalar factors ρ,σ of (1.21). The infection potential, f(ρ), is constructed via numeric interpolation of functions (1.20), for computed snail infection prevalences ŷ(ρ) of (1.12).

Evaluation of control strategies based on drug therapy

We shall consider specifically, three cases of therapy-based control. (1) Age targeted strategy where all locales are treated identically, so all treatment matrices are equal {Θ_i=Θ; i=1,…, M}. (2) The general case that allows for site-specific and age-specific treatment at high risk locales. (3) An example of site-specific treatment of high risk locales (sites 4, 10, 11 on Fig. 2) treated in isolation.

Case 1

We apply the above ‘mean burden’ formulae (1.21)–(1.23) with teatment-augmented attrition matrix Γ+Θ. The corresponding scalars ρ(1.27) and σ(1.22) will depend now on treatment frequencies {τ_a}, hence the resulting community mean burden is

(τ₁,…,τ_P). We compute numerically equilibrium snail infection prevalences {y_j(ρ)} for a range of B.R.S. values:

, and find the corresponding (normalized) infection potential f*(ρ) of (1.20). Fig. 4 shows function f(ρ) for the model environment of Fig. 2, while Fig. 5 demonstrates site distribution of equilibrium snail infection prevalences {y_j*(ρ)} and worm burdens {w_i*(ρ)} for several selected values of the basic reproductive scalar, ρ.

Fig. 4. Infection potential f(ρ) (1.20), as a function of the basic reproduction scalar ρ, for model environment of Fig. 2.

Fig. 5. Site distribution of relative infection levels at equilibrium for different model locations, based on area transmission potential: mean worm burden (left column), and snail infection prevalence (right column) at selected values of normalized reproduction scalar ρλ1=3 (top panel); 5·5 (middle panel); and 8 (bottom panel).

Case 2

One can target chemotherapy at both high-risk population subgroups and/or locations. The location risk factor clearly exhibits itself in the elevated infection levels found at i=4, 10, 11, 12 (Fig. 5).

A general approach would require an array of treatment matrices {Θ₁,…, Θ_M} adopted for individual sites. Therefore, the basic reproduction matrix (1.11) would depend on several treatment frequency vectors

. The increased number of control parameters, however, is only part of the problem. More importantly, mathematical formulation of site-specific control problem becomes more involved. We provide details in Appendix 1, and here just state the results. The basic reproduction matrix will depend now on several parameters: R=ρR⁰ (r₁,…), namely the previous B.R.S., ρ, and additional (site-specific) parameters

. One can think of ρ as B.R.S. of untreated populations, while

measures the long-term effect of focal treatment at site {i}. Hence equilibrium snail-site infection prevalences become multivariable functions ŷ(ρ; …, r_i,…), and the corresponding community worm burden takes on the form

Here, σ corresponds to untreated value (1.22),

gives the effect of local treatment, and

The treatment frequencies

enter such

through variables r_i, s_i, and we can formulate the corresponding control/optimization problem similar to Case 1. We do not pursue it here for two reasons (i) extra variables

add a substantial computational burden, (ii) in future programming development we do not foresee that this strategy wins compared to the age targeted (‘uniform’) teatment scheme of Case 1. Rather than the full-scale optimization, we shall estimate the efficacy of site-targeted control at high-risk locales, comparing them to untreated populations (see Case 3, below).

Case 3 (site-specific chemotherapy)

We want to examine the effect of selective treatment at high-risk locations. Our environment (Fig. 2) suggests that sites 4, 10, 11, and 12, close to a high-density snail site I, are high-risk locations, and Fig. 4 (untreated equilibria) corroborates such conclusion. For analysis of site-specific treatment, we choose sites i=4;10;11, hence an array of treatment matrices

According to (1.28) the community mean burden will take on the form

As discussed above for Case 2, the extra variables {r_i,s_i} represent ratios of treated/untreated values of B.R.S., ρ and σ. Clearly, a treatment can only decrease ρ_i and σ_i, according to (1.16) and (1.22), so variables r,s change between 0 (extreme treatment) and 1 (untreated). They go down to zero with increasing frequencies,

. Thus, we get a crude estimate of the efficacy of focal treatment (and its effect on untreated populations) by comparing the two extremes.

Synthesis of chemotherapy strategies: analysis and optimization

Schistosomiasis can produce both short-term (acute) effects and long-term (chronic) illness, exemplified by hepatomegaly and Symmers's fibrosis. Here we consider only its short-term effects, and use worm burden as a crude proxy of acute morbidity, as correlated with the egg count and consequent illness. For other models of morbidity and its long-term (chronic) effects, we refer to (Medley & Bundy, 1996; Chan & Bundy, 1997). In particular, Chan et al. (1997) proposed, based on their analysis, a mass-treatment strategy comprised of 4 year treatment cycles administered over a 1 to 21 year span.

Here, we propose a similar strategy targeted at high-risk/exposure groups but, instead of a fixed 4-year cycle, we shall adjust treatment frequency {τ_a} for each age group, so as to minimize the overall mean worm burden of the entire community.

The important constraint for mass chemotherapy in many developing countries is limited available resources. In asking for the most efficient distribution of treatment efforts among various risk groups, we aim to attain the minimal community worm burden within a given cost constraint. We shall formulate the corresponding constrained optimization problem, assuming as above that the long-term effect of time-dependent (periodic) chemotherapy amounts to changing the natural worm attrition matrix Γ by the treatment matrix Γ+Θ, Θ=diag(τ_a).

The optimal control strategy then asks to distribute the available teatment resources among affected human groups, so as to attain the lowest community burden

(τ₁,…,τ_P), subject to non-fixed (per capita) cost constraint C(τ₁,…,τ_P)=C₀, and the natural time constraints {0<τ_a<τ_max}.

Although a treatment programme carries both fixed and non-fixed expenses, the non-fixed cost of treatment is clearly proportional to its frequency and the size of the targeted populations. Hence, the per capita non-fixed cost function in terms of population fractions, {h_a*} becomes

where c₀ denotes per capita cost of each treatment dispensed, and C₀ is the amount actually allocted for control on a per capita basis. Here we used simple considerations of cost that result in a linear constraint (1.30). Other considerations, e.g. fixed costs or economies of scale, could produce nonlinear

, but, computationally, a particular form of C makes no difference in our approach to numeric solution.

For our analysis, we developed and implemented a numerical scheme based on Mathematica software package that computes equilibria ŷ(ρ) for a range of values 1/λ₁(R⁰)=ρ_min<ρ<ρ_max, and runs optimization routine of the resulting mean burden

=σf(ρ).

In our model environment (Fig. 2), we stratified the population into 12 age-bins (with step Δa=5 years), and adopted the population data, contact patterns ω, and contamination rates β (Table 1), by the continuous model of Chan & Bundy (1997). Chan & Bundy suggested functional form:

,

with properly adjusted coefficients A_i, B_i. Their estimates are derived, initially, from data of the long-term study of S. mansoni transmission in St Lucia (Jordan, 1985), and validated against data from field studies of S. mansoni control in Kenya (Butterworth et al. 1991).

Based on these inputs, we computed and compared the age burden distribution (averaged over the entire community) without treatment and with optimal treatment at a fixed level of cost allocation, C₀/c₀=0·2, or 20 of the cost of total population coverage, and two distinct values of (untreated) B.R.S.: low values ρ=5ρ_min (left panel of Fig. 6), and relatively high value ρ=20ρ_min (right panel). Here ρ_min=1/λ₁(R⁰) is the minimal (bifurcation) value of B.R.S. required for the establishment of endemic infection equilibria. Although relatively low in the allocated cost ratio, such age-based optimized treatment would still allow a significant reduction of worm burden across all ages, but most significantly for the younger groups. The corresponding optimal treatment intervals {τ_a} for the first 6 age bins are shown in Table 2 for low and high B.R.S., ρ=5ρ_min and 20ρ_min. The optimization scheme here was applied only to the 6 youngest groups (below 30 years old), leaving older groups untreated.

Fig. 6. Comparison of age distributed equilibrium worm burdens {wa} wihout treatment (grey), and with optimal frequency age-specific treatment (black), at allocated/full treatment-cost ratio: C0/c0=0·2, for two values of basic reproduction scalar (B.R.S.), ρλ1=5 (low), and ρ=20ρmin (high).

Table 2. Optimal treatment periods (in years) for the first 6 age bins, at 20 allocation of full population coverage, at two values of the area-wide basic reproduction scalar (B.R.S., (1.16) ρ=5ρmin and 20ρmin. The last column shows the residual fraction of the community mean worm burden, , after implementation of drug-based control for those two settings.)

The next plot (Fig. 7) shows the optimally reduced community worm burden

of the entire population (solid black), and the corresponding treatment frequencies {τ_a} of the first 6 age groups (dashed curves of increasing dashing). We plot them over a range of the relative cost allocation parameter ξ=C₀/c₀, namely, 0·15<ξ<0·5 – for low B.R.S. ρ=5ρ_min, and 0·2<ξ<0·95 for high B.R.S., ρ=20ρ_min.

Fig. 7. Effect of optimized treatment on community mean worm burden (solid black). The optimal therapy frequencies {τa} for the 6 youngest age groups (dashed curves), are shown as functions of available funding, in terms of the ratio (allocated cost/cost of total population coverage) ξ=C0/c0, shown on the ordinate. In the setting where only 50 or fewer can be covered, then annual treatment of 6–15 year age groups (τ=1), combined with treatment every 4–10 years for younger and older age groups (τ=0·25 to 0·1), provides the maximal possible reduction in community worm burdens.

Let us note, that the optimization procedure must be terminated when the community mean burden

falls to a low (near eradication) level. So the former case (low B.R.S.) is terminated at about ξ=0·5, which means 50 population coverage could practically eradicate the infection. In the latter case (high B.R.S.) it would take a substantially higher value ξ>0·9. In either case, we see a significant reduction of

at moderate values of ξ.

Next we turn to control strategies targeting high risk sites, such as Site 4 on Fig. 2, having the greatest mean untreated worm burden (w-level, Fig. 5). Note that extreme intervention at selected locales will set a lower bound (on worm burden) for any other control strategy, including mixed age- and site-specific optimal control. Fig. 8 compares the community mean burdens

, as functions of the natural (pre-treatment) B.R.S. ρ, in two cases: (i) untreated population (thick line), which corresponds to r_i=s_i=1, and (ii) ‘extremely treated’ site 4 (thin curve), i.e. r₄=s₄=0 in (1.29). The overall reduction is marginal, about 5–10 of the total mean burden, in curing the same proportion of the human population. The impact is reduced in the setting of increased ρ, and does not appear to offer any special leverage for effecting control.

Fig. 8. Community mean worm burden as a function of the basic reproduction scalar ρ0 in two cases: (a) untreated population (thick curve), and (b) location-specific treatment with village site 4 infections removed through eradication τa→∞ (thin curve). The difference between the two curves shows the best possible gain in overall area worm burden (about 10 reduction) attainable through the targeted treatment of site 4 alone. The dashed curve shows the relative overall gain in area worm burden between the two cases. It falls rapidly with ρ, so the targeted treatment of site 4 alone gives <10 reduction at moderate or higher values of reproduction scalar ρ (>30).

Fig. 9 compares site burden distribution

in two cases: ‘extreme control’ at sites i=4, 10, 11 (light gray) vs. ‘untreated populations’ (dark), for two values of the B.R.S. ρ. Again the gain at all sites (excepting 4, 10, 11) is marginal (far less than 50 reduction in worm burden), and becomes even less pronounced at larger ρ.

Fig. 9. Site distribution of worm burden for pre-treatment human populations (light bars), and the effect of extreme treatment for sites 4, 10 and 11 on the remaining untreated sites (dark bars). Left panel indicates overall impact under low transmission potential (ρ=7ρmin). Right panel indicates the lesser impact for untreated sites under conditions of high transmission potential (ρ=20ρmin).

We propose a plausible explanation for apparent inefficiency of the site-targeted control in our setup. This has to do with an overall connectivity of the snail-human system, that links (with nonzero contact probability η_ia,j>0) any human site i and snail site j. By selectively suppressing sites i=4, 10, 11, but leaving other transmission channels open, the system can effectively restore the pre-treated high levels of infection. For comparison, the age control strategy, in particular, an optimal one applied across the entire region, appears far more efficient.