Static environment

If we assume the completely symmetric case with no fitness differences between variants in a static erythropoietic environment, i.e. f_i=f, ∀_i and ∂f/∂t=0, we find the dynamics quickly converging towards one of the possible equilibrium states. System (1)–(3) possesses a disease-free equilibrium

and linear stability analysis shows that this point is unstable as long as

where R can be defined as the pathogen's basic reproductive number. Assuming that (5) holds, we find the system reaching a non-zero state of equilibrium. In the case of a constant invasion profile, i.e. f=K, this is given as

This is interesting as neither the density of infected red blood cells nor the density of free-living merozoites is dependent on the number of merozoites released from an infected red blood cell. If we assume that the invasion profile of free-living merozoites is density dependent, the function f can be defined as

where δ can be taken as the maximum invasion rate in the absence of resource limitation. We assumed that all infected erythrocytes have a negative effect on the invasion success of merozoites of type i. In this case we find two non-trivial fixed-point solutions to system (1)–(3), one of which is density controlled and independent of the immune response:

where n is the number of phenotypes. This means that strong resource limitation can lead to self-limited parasitaemia. Incorporating fitness differences by making f variant specific, for example ordered in a hierarchical manner f₁<f₂< …, does not alter the above results other than leading to structured equilibria in the RBC and antibody population, where phenotype prevalence is in order of their respective fitness. Interestingly, though, free-living merozoites will exhibit no such structure but equilibrate towards the same frequency, despite fitness differences (results not shown here).

Dynamic environment

Rather than assuming a static composition of the red blood cells, we can incorporate a dynamic erythropoietic environment by making f explicitly time-dependent. By taking

with positive constants Δ and c, we make sure that f_i[ges ]0 and that the average availability of red blood cells for invasion stays constant over time, although this is not a necessary requirement to achieve similar results.

For low values of ε the system will quickly converge towards low amplitude oscillations about the fixed point of the equivalent autonomous system (1)–(3), with ∂f/∂t=0. By increasing the frequency at which a merozoite subpopulation is dominated by its originating type, i.e. increasing ε, we can observe sustained oscillations representing sequential selection and amplification of single phenotypes, as shown in Fig. 1. Numerical simulations also reveal a positive correlation between ε and peak heights. That is, clonal selection is more efficient the more clones of a particular variant are produced at each reproduction cycle. In this respect, the clonal phenotypic variation process can enhance the merozoites' survival success by allowing them to adapt to temporal environmental changes. Although the oscillatory behaviour can be attributed to the time-periodic invasion function f, a positive selection process whereby single variants gain an advantage over others is nevertheless apparent. The fluctuating erythropoietic environment thus serves an exogeneous pump to reinforce the natural propensity of the system to oscillate, in a manner analogous to seasonal forcing or demographic stochasticity in simple epidemiological models (as reviewed recently by Levin, Dushoff & Plotkin, 2004).

Static environment

If we assume that the invasion properties of the merozoite are not density dependent, i.e. f=K, we find that for low values of γ, the system will settle into a non-zero steady state with all strains being present at same equilibrium frequency, given as

where η is the number of variants that induce cross-reactive antibodies against i (here η=2). Increasing γ, this steady state looses its stability and is replaced by stable limit cycles that represent sustained oscillations of single merozoite variants. These oscillations, shown in Fig. 2, demonstrate the sequential selection and amplification of variants that have a temporal advantage over other variants as their growth is not restricted by antibodies raised against the preceding type(s). That is, because of the cross-inhibitory relationship between variant types defined by Γ, if we would take type i as the peaking variant at any one time it would be followed by one that is not negatively affected by the immune response launched against i or its preceding types. This scenario leads to a traveling wave of single types.

Fig. 1. Time-series of free-living merozoites, mi, showing a pronounced selection process of single variants for increasing ε. ε=0·4, top, and ε=0·8, bottom. Other parameter values: μm=1 [h−1], μr=1/30 [h−1], μa=5E-3 [h−1], β=1E-4, ρ=20, Δ=120, c=0·2, κ=0·1.

Fig. 2. Depending on the degree of immune pressure, in terms of cross-immunity γ, the system will either settle down onto a stable non-zero equilibrium (top, γ=0·5) or exhibit sustained oscillations of single types (bottom, γ=1). μm=1 [h−1], μr=1/30 [h−1], μa=5E-3 [h−1], β=1E-4, ρ=20, κ=0·1, f=K=10.

The net production rate, ρμ_r, or rather (log (ρ)+1)μ_r, does not have a significant impact on the dynamics, whereas the latter is crucially dependent on the composition of the merozoite subpopulation, ε. Decreasing the frequency of a predominantly produced merozoite phenotype (decreasing ε) will lead to a stable equilibrium, whereas an increase will tend to produce much sharper dynamics in terms of period and amplitude. With a more homogeneous subpopulation (high ε) single types cannot be amplified within the population as all antibodies are present at high level at all time. On the other hand, producing only merozoites of the invading type is favourable in a selection process, although very high peaks could be subject to a trade-off in terms of exhausting the host's resources. The parameter space in Fig. 3 shows how an increase in γ and ε will lead to sustained oscillations of single merozoite types that represent the clonal selection process.

Fig. 3. (γ,ε) – parameter-space showing the region of clonal phenotypic selection. The sequential amplification of single types is promoted by a strong immune selection pressure (high γ) and strongly biased transcription (high ε). μm=1 [h−1], μr=1/30 [h−1], μa=5E-3 [h−1], β=1E-4, ρ=20, κ=0·1, f=K=10.

If we include density-dependent invasion we find that for intermediate or lower degrees of invasion regulation (high values of C), sustained oscillations can still be found in higher regions of cross-protection, whereas low ε and γ will lead to a state of equilibrium. Therefore, provided the invasion success is not too restricted by the number of infected red blood cells, phenotypic variant selection can still be found, as shown in Fig. 4. Adding host-heterogeneity to the density-dependent invasion regulation, in terms of merozoite fitness differences, will lead to similar dynamical behaviour, except that sustained oscillation can occur even in lower regions of γ and peak heights are restricted by a self-regulatory effect.

Fig. 4. A strong dependency on the number of red blood cells already infected will self-limit the disease, independent on the host's immune response (top, C=1E2). For higher values of C phenotypic variation can be observed (bottom, C=1E4). Other parameter values: μm=1 [h−1], μr=1/30 [h−1], μa=5E-3 [h−1], β=1E-4, ρ=20, κ=0·1, δ=10, γ=0·8, ε=0·9.

Dynamic environment

Introducing a dynamical erythropoietic environment we find that the dynamics do not depend that critically on the degree of cross-protection anymore, although an increase in peak heights is still observable with increasing γ. That is, because of the non-autonomous (periodically driven) invasion profile, oscillations are observable even in the very low regions of γ. The weighing factor ε, on the other hand, significantly determines the qualitative behaviour: from low-amplitude oscillations around the associated fixed point to sequential amplification of single types with an increase in ε. In the extreme where ε=1, after an initial peak some of the variants will be driven to extinction by a very strong immune response. Increasing the immune pressure upon the merozoite population (increasing γ) we see more and more variants being driven extinct, and because phenotypes are not being produced by any other types, they are permanently eliminated (see Fig. 5). Under these circumstances the benefit of producing a mixed merozoite subpopulation from a single schizont would be to avoid severe diversity reduction under increased immune pressure.

Fig. 5. Increasing the immune pressure will lead to diversity reduction within the merozoite population, (γ,ε)=(0·2,1), (0·2,1), ( 0·6,1) and (0·8,1) (top left to bottom right). μm=1 [h−1], μr=1/30 [h−1], μa=5E-3 [h−1], β=1E-4, ρ=20, κ=0·1, Δ=120, c=0·2.

A note on variant order

The sequential dominance of single merozoite phenotypes observed seems to follow a particular order. So far, no variant order has been observed within real merozoite populations (Preiser et al. 1999), as the sequential amplification of certain types is still only a hypothesis. In our model, the order of phenotypes that are amplified over time is determined by the non-reciprocal form of cross-protection matrix Γ, and, in the case of a dynamic environment by the invasion function f. Looking at Γ we see that antibodies against type i cross-protect against types i+1 and i+2. Growth of type i+3, for example, is thus not constrained by an immune response against i and can therefore be amplified within the population. (Note, in case i=n, i+1=n+1=:1.) Introducing a time-dependent composition of the red blood cell population promotes oscillations of single variants in order of the distribution defined by f. Growth of variants that are fittest at any one point will thus be promoted. Interestingly though, the order in which variants peak is changed when immune pressure is increased (Fig. 6). Two different selection phenomena are displayed by the system: one which can be attributed to immune evasion (variant order determined by Γ), and the other one to host/environment adaptation (variant order determined by f). We therefore note two driving forces for positive selection that could contribute to a successful evasion and adaptation strategy.

Fig. 6. Low degree of immune pressure will preserve the variant order (arbitrarily chosen 1→2→ … →8→1 … ) imposed by f, γ=0·3 (top). An increase results in an order such that succeeding variants do not interfere with each others growth, determined by cross-protection matrix Γ (γ=0·8, bottom). ε=0·9, μm=1 [h−1], μr=1/30 [h−1], μa=5E-3 [h−1], β=1E-4, ρ=20, κ=0·1, Δ=120, c=0·2.