Dynamics of immune response and drug resistance in malaria infection

Abstract

Background

Malaria parasites that concurrently infect a host compete on the basis of their intrinsic growth rates and by stimulating cross-reactive immune responses that inhibit each others' growth. If the phenotypes also show different drug sensitivities ('sensitive' vs. 'resistant' strains), drug treatment can change their joint dynamics and the long-term outcome of the infection: most obviously, persistent drug pressure can permit the more resistant, but otherwise competitively-inferior, strains to dominate.

Methods

Here a mathematical model is developed to analyse how these and more subtle effects of antimalarial drug use are modulated by immune response, repeated re-inoculation of parasites, drug pharmacokinetic parameters, dose and treatment frequency.

Results

The model quantifies possible effects of single and multiple (periodic) treatment on the outcome of parasite competition. In the absence of further inoculation, the dosage and/or treatment frequency required for complete clearance can be estimated. With persistent superinfection, time-average parasite densities can be derived in terms of the basic immune-regulating parameters, the drug efficacy and treatment regimen.

Conclusion

The functional relations in the model are applicable to a wide range of conditions and transmission environments, allowing predictions to be made on both the individual and the community levels, and, in particular, transitions from drug-sensitive to drug-resistant parasite dominance to be projected on both levels.

Background

Since the 1920s it has been clear that interactions of immune-system and drug dynamics are critical to the elimination or persistence of drug-resistant malaria parasites in individual patients, but details and dynamic principles remain to be established [1, 2]. Later, in the 1980s, it was found that an isolate from a malaria infection is typically a mixture of parasite genotypes, some of which differ in antigenic profile, growth rate and drug response, and show complex population dynamics in mixed cultures [3, 4]. Here a mathematical model is developed to investigate the dynamics of parasite phenotypes in a malaria-infected host, with respect to critical interactions between their immune-mediated competition, relative drug sensitivities and persistent superinfection. The results characterize conditions under which a malaria infection dominated by relatively drug-sensitive parasites transitions to one dominated by relatively drug-resistant parasites.

Antimalarial drugs taken prophylactically or early in infection limit the development of immune responses protective against clinical attacks in subsequent infections [5, 6]. In endemic regions, antimalarial drug treatment is usually meant to cure a symptomatic malaria infection, thus, usually, a host with incomplete clinical immunity. Hosts are repeatedly re-infected and superinfected with mixtures of parasites. Within an infection, in the absence of drug, the dynamics of the competing parasites are determined in large part by their differences and similarities in stimulating and succumbing to the host's immune responses. The frequency with which drug clears parasites, including drug-resistant parasites, increases with host age, a surrogate for exposure and acquired immunity [7, 8]. Stronger, more diverse antibody responses are associated with greater success in antimalarial drug treatment [9, 10]. Thus, for instance, all else equal, both a higher prevalence of drug-resistant parasites and a higher frequency of drug treatment would be expected in younger age groups.

Understanding the aggregate effects of immunity can help to guide the introduction of new antimalarial drugs, explain discrepancies between in-vitro and in-vivo tests in drug-resistance surveys [11], and, in general, inform population-based intervention strategies. Indeed, these community patterns summarize very broad effects of immunity across many infections, so their interpretation and application can be greatly improved by a more detailed understanding of immune-drug-parasite interactions within individual infections.

Most studies of intra-host parasite-drug dynamics focus on the clearing effect of drug (the key medical objective) or on optimal strategies for treatment (dosage, frequency). These issues become less relevant under persistent superinfection in endemic areas, however. The focus here is more on longer-term outcomes of treatment (or prophylaxis), and their effects on competition among parasites. For instance, when joined to a proper model of parasite transmission, the crude estimates derived in relating the mean parasitaemia of phenotypes in a 'typical' host to drug dose and frequency should allow community-level analysis of the spread and prevalence of drug -resistant parasites.

The model developed here builds on earlier models of the within-host dynamics of malaria [12, 13] to address critical interactions between immune response and drug resistance. In the baseline model, the relative cross-reactivities, immune efficiencies and growth rates of competing parasites determine their dynamics and joint equilibria. Further developments examine the effects of repeated re-inoculations of parasites, then incorporate differences in drug sensitivity between the competing parasites, with the costs of resistance summarized by differences in intrinsic growth rate, to examine how dose, decay rate and frequency of drug administration affect competitive dynamics.

Methods

Model of intra-host competition mediated by immunity

The basic model of intra-host competition of multiple parasite phenotypes is based on the models of cross-reactive 'stimulating-clearing' immunity formulated in previous papers [12, 13]. Briefly, it describes interactions in which parasite phenotypes stimulate the production of non-specific and specific immune effectors, which clear (kill) them at rates that depend on the strength, specificity and quantity of effector.

Mathematically the model is a coupled system of differential equations for parasite densities x, y, and immune-effector densities: non-specific I and specific J,K for x,y-species respectively (see [12, 13]),

$\begin{array}{l}\dot{x}=\left[a_X-f\left(\frac{x}{x_p}+\frac{y}{y_p}\right)-c_{n}I-c_X\left(J+\eta K\right)\right]x+S_X;\\ \dot{y}=\left[a_Y-f\left(\frac{x}{x_p}+\frac{y}{y_p}\right)-c_{n}I-c_Y\left(\xi J+K\right)\right]y+S_Y;\left(1\right)\\ \dot{I}={\sigma }_n\left(x+y\right)-{\mu }_{n}I;\dot{J}={\sigma }_{X}x-{\mu }_{s}J;\dot{K}={\sigma }_{Y}y-{\mu }_{s}K;\end{array}$

In the follow up analysis of drug treatment x will represent a resistant strain, while y -a sensitive one. This is similar to the basic model of [12] but augmented with two additional features:

(i) a 'fever regulator' f(z) = F₀φ _q (z), described by a sigmoid function φ _q (z) = z^q/(1 + z^q), switched 'on' and 'off by the combined parasite density (over their pyrogenic thresholds x _p , y _p ), with maximal pyrogenic removal rate F₀.

(ii) super-infection (inoculation) sources S _X ;S _Y , which can be either stationary (convenient for qualitative analysis and numeric simulation), or stochastic with a prescribed 'spacing distribution' between inoculations.

The fever plays a limited role in the dynamics of (1). At the initial stage of infection it serves to arrest the parasite growth about the pyrogenic threshold, while immune effectors build up, but then becomes negligible below that level. Therefore it is omitted in the analyses of subfebrile equilibria, but is included in the dynamic simulations.

The immune-stimulation terms (in the I,J,K – equations) are taken here as simple linear functions, so that the effectors are stimulated in proportion to parasite densities. Other parameters in system (1) are similar to [12, 13]:

(i) Parasite growth rates: a _X ;a _Y

(ii) Stimulation coefficients, for nonspecific and specific effectors: σ _n ;σ _X ;σ _Y

(iii) Clearing coefficients, for nonspecific and specific effectors: c _n ;c _X ;c _Y

(iv) Decay/deactivation rates, for nonspecific and specific effectors: μ _n ;μ _s

(v) Cross-reactivities, in terms of the relative clearing rates of the x-effector on y and the y-effector on x: 0 ≤ ξ, η < 1

System (1) can be rescaled to a nondimensional form which yields the essential (non-dimensional) parameters. The most important of those are cross-reactivities and the immune efficiencies (non-specific e _n and specific e _x , e _y ) introduced in [13]. The immune efficiencies combine the appropriate 'stimulation' and 'clearing' coefficients against the product of 'immune decay' and 'parasite growth' into a single nondimensional parameter.

System (1) can be considered a 'multi-prey, multi-predator' model of parasites and immune effectors, or alternatively as 'immune-mediated' competition between parasites x and y. The outcome of such competition depends on phenotypic traits with respect to the host immune response, as expressed through the essential, non-dimensional parameters above (cross-reactivities and immune efficiencies), as well as relative parasite growth rate α = a _Y /a _X . Note, for instance, that different parasite species may differ in replication rates (a _X ≠ a _Y ) and show low cross-reactivity, while different "strains" of a species (e.g. drug-sensitive or drug-resistant/tolerant) may be more strongly cross-reactive

Basic pharmaco-dynamics of drug treatment

The next development is a simple pharmaco-dynamic model for a single parasite phenotype. The degree to which drug treatment reduces the growth rate of a given parasite depends on both the pharmacokinetic characteristics of the drug and the sensitivity of the parasite to the drug. Key assumptions are that the phenotypic cost of drug resistance is a lower intrinsic growth rate in the absence of drug, and that the drug-induced reduction can be approximated from the basic drug parameters and the frequency of (periodic) treatment. The effect of 'heavy treatment' may be to remove the more sensitive strain completely, or, for repeated treatment, with superinfection, to tip the effective relative growth rate α in favour of x (below its critical levels α_1,2, explained in the following sections), releasing the more resistant strain from competitive constraints.

The pharmacodynamic model adopts the basic premises of [14]:

(i) after intake, drug concentration decays from its initial dose d₀ at an exponential rate d(t) = d₀e^-βt(for simplicity, blood and tissue decay are not distinguished);

(ii) the drug removes parasites through a 'clearing function' Bφ _p (d/d _S ), where B denotes the maximal clearing rate, d _S the parasite 'sensitivity threshold' (often denoted d₅₀, as it clears parasite at 50% maximal rate), and φ _p (z) – a sigmoid function φ = φ _p (z) = z^p/(1 + z^p) with a suitable Hill exponent p.

A large p in φ means that the clearing rate remains relatively steady ≈ B over a wide range of drug concentrations, provided that d stays above d₅₀. For the numeric simulations below, the exponent p = 3, close to its estimated value for mefloquine [15].

The resulting dynamic treatment model [14] with constant parasite growth and a source (inoculation) S, takes the form

$$\dot{x}$$

(1)

= a[1 - bφ(D(t))]x + S(t)     (2)

where D(t) = d(t)/d _S is the dimensionless drug concentration, and b = B/a the relative parasite clearing rate. Two special cases are:

(i) D(t) = $D_0{e}^{-\beta \left(t-t_0\right)}$ – for a single dose administered at t₀, written in terms of a dimensionless initial dose D₀ = d₀/d _S ≫ 1;

(ii) a periodic function (equation 5; presented below) for multiple applications with period T.

Equation (2) contains the essential pharmaco-kinetics of drug (the relative parasite clearing rate b and 'clearing pattern' φ), as well as its pharmaco-dynamics (initial dose and the frequency of treatment encoded in D(t)).

Drug treatment will enter system (1) or its rescaled version (3) as additional attrition terms: b_X/Yφ(D_X/Y(t)) in the x, y equations. Two phenotypes are assumed to have different drug-sensitivities, with x the more tolerant/resistant and y the more sensitive strain. Hence the different clearing coefficients, b _X <b _Y (to indicate higher maximal removal rate of y), and different drug-thresholds, expressed through their relative initial dosages, D _X <D _Y (y being cleared at lower concentrations of drug, compared to x). Thus strain y has a higher natural growth rate, α = a _Y /a _x > 1, and is competitively superior in the absence of drug, but, as shown below, drug treatment can offset this advantage.

Results

Analysis of stationary equilibria

The analysis at this stage is meant to outline the long-term effect of parasite-immune interactions. In particular, it attempts to predict the resulting parasite densities (long term 'outcome of competition') through the basic host parameters (growth rates, immune efficiencies, cross-reactivities), and does not yet involve drug treatment.

The rescaled system (1) with drug treatment represented by time-dependent concentration D = D(t) is given by a coupled differential system

$\begin{array}{l}\dot{x}=a_X\left[1-f\left(x+\frac{y}{y_p}\right)-b_X\phi \left(D\right)-e_{n}I-e_X\left(J+\eta K\right)\right]x+S_1;\\ \dot{y}=a_X\left[\alpha -f\left(x+\frac{y}{y_p}\right)-b_Y\phi \left(D\right)-e_{n}I-e_Y\left(\xi J+K\right)\right]y+S_2;\left(3\right)\\ \dot{I}={\mu }_n\left[\left(x+y\right)-I\right];\dot{J}={\mu }_s\left(x-J\right);\dot{K}={\mu }_s\left(y-K\right);\end{array}$

It depends on several dimensionless parameters: fever efficiency $f_0=\frac{F_0}{a_X}$, cross-reactivities 0 ≤ ξ, η < 1, and immune efficiencies: $e_n=\frac{c_n{\sigma }_n}{a_X{\mu }_n}{x}_p$, $e_X=\frac{c_X{\sigma }_X}{a_X{\mu }_s}{x}_p,e_Y=\frac{c_Y{\sigma }_Y}{a_X{\mu }_s}{x}_p$. To facilitate analysis we drop the fever term, assuming 'subfebrile equilibria' x*, y* ≪ 1 that result from (relatively) high immune efficiencies e _n ;e_X,Y≫ 1. The equilibrium equations for (3) without treatment or fever, are then reduced to the classical Volterra-Lotka system

$\begin{array}{c}F\left(x,y\right)=\left[1-\left(e_X+e_n\right)x-\left(\eta e_X+e_n\right)y\right]x+S_X=0\\ G\left(x,y\right)=\left[\alpha -\left(\xi e_Y+e_n\right)x-\left(e_Y+e_n\right)y\right]y+S_Y=0\end{array}\left(4\right)$

In the absence of sources (S_X,Y= 0), it has three types of equilibria (Figure 1). Two of them, 'x-domination' (x₁, 0) or 'y-domination' (0, y₂), are given by coordinate intercepts of the 'linear factors' in the F-equation: x₁ = 1/(e _X + e _n ); y₁ = 1/(ηe _X + e _n ), and the G-equation x₂ = α/(ξe _Y + e _n ); y₂ = α/(e _Y + e _n ). The third possibility is a state of 'coexistence' (x*, y*). The stability and qualitative behaviour of equilibria depend on the relative position of two null-clines (4). In our case, the equations y₁ = y₂ and x₁ = x₂ give two critical (bifurcation) values of the relative growth-rate (fitness parameter) α = a _Y /a _X , namely

Figure 1

Phase-plane of dynamical system (4). (a) Phase-plane views of 3 regions of system (4) in terms of relative growth rate α : y-domination (equilibrium y₂, left); coexistence (equilibrium I, middle) and x-domination (equilibrium x₁, right). Stable equilibria are marked in black. In all 3 cases, adding sources (inoculations) produces a stable coexistence equilibrium (intersection of asymptotic hyperbolae). The two solid lines are x-nullclines F(x, y) = 0, the two dashed lines are y-nullclines G(x, y) = 0. (b) Equilibria (4) for stationary sources S_X,Y(x(ε) – solid, y(ε) – dashed), with the parameters of Table 1, in 3 different ranges of α, as functions of EIR, 0 <ε < 100. Two sets of curves on each plot correspond to different partitions of EIR among strains: (i) p _X = p _Y = .5 (outer curves), (ii) p _X = .9; p _Y = .1 (inner 'x-biased' curves). The two solid lines are x-strains, and the dashed lines are y-strains.

${\alpha }_1=\frac{e_Y+e_n}{\eta e_X+e_n}>{\alpha }_2=\frac{\xi e_Y+e_n}{e_X+e_n}\left(5\right)$

The three fitness regions include:

$\begin{array}{l}\alpha >{\alpha }_1\text{-}\text{'}\text{y-domination}\text{'}\text{;}\\ {\alpha }_{\text{1}}>\alpha >{\alpha }_2\text{-coexistence;}\left(\text{6}\right)\\ {\alpha }_2>\alpha \text{-'x-domination'}\end{array}$

Here 'domination' in the absence of inoculation sources, means complete removal of the competitor. The familiar phase-plots of such competition are illustrated in Figure 1.

The presence of stationary sources of infection S _X ; S _Y > 0 (no matter how small) will shift all three types of equilibria to the upper-right quadrant into a stable 'coexistence state', defined by two asymptotic hyperbolae (solid-x and dashed-y, Figure 1). Now 'domination' does not mean complete removal, but a 'dominant density', e.g. y*/x* ≫ 1. The use of α as the 'control parameter' is relevant for qualitative analysis of drug intervention (below). Indeed, given two parasites with different drug-sensitivities, drug treatment effectively lowers their growth rates, and thus tips the outcome of competition in favour of the resistant strain.

This qualitative analysis makes two simplifying assumptions – a reduced 2D model in place of the full 5D (3), and the omission of fever. The latter has minor consequences for high immune efficiencies. The dimensional reduction (5D to 2D) maintains the equilibrium values x*, y*, but changes their stability types somewhat, from '2D stable nodes' (real negative eigenvalues of the Jacobian matrix), to '5D spiral sinks' (complex negative eigenvalues), (see [13] for details).

Dynamic relaxation

Next, the dynamics of this baseline model are examined by numerically computing solutions of rescaled system (3) with and without inoculation sources, using the parameter values described in Table 1. Some parameter values are based on earlier work, e.g. the relative pyrogenic clearing rate [18] and immune loss rates [12]. The essential 'uncertain parameters' are the immune efficiencies and cross-reactivities. Because the immune efficiencies control equilibria as ≈ 1/(e _n + e_X,Y) we chose {e _n ,e_X,Y} sufficiently large to keep equilibria below pyrogenic levels, as expected for a typical 'asymptomatic parasitaemia'. The exact choice of e – values makes little qualitative difference, provided equilibria stay below the 'pyrogenic level' x _p , but could affect the pharmacokinetic parameters (doses) for parasite clearing, as explained in the next section.

Table 1 Basic parameters for dynamic model (3)

The role of cross-reactivities in the equilibrium analysis above is to narrow (ξ, η ≈ 0) or widen (ξ, η ≫ 0) the α - range of coexistence. Intermediate values are chosen to indicate a relative genetic (antigenic) proximity of the competing strains.

Under this choice of parameters the predicted thresholds for the relative growth rate α (5)-(6) become α₁ = 2.17; α₂ = .54. Also introduced is a lower cut-off for parasite densities, x _c = 10^-11 (in dimensionless units), assuming a pyrogenic threshold x₀ = 10⁴/μl, corresponding to a density of 10^-7/μl (i.e. fewer than one parasite for the entire blood volume of an adult). If either of the densities x,y falls below x _c it is set to zero – i.e. complete clearing.

The inoculation sources are considered either stationary, of strength EIR × "injected density s₀", or 'random,' with exponentially distributed waiting times and mean spacing EIR ('entomological inoculation rate,' i.e. the frequency of infectious bites by mosquitoes). The 'injected density' is estimated in terms of the parasites (primary merozoites) released from the liver, assuming a mean of 10 parasites (sporozoites) per mosquito inoculum, with each sporozoite developing into 30,000 primary merozoites, thus a total 300,000 per 4.5 liters, or 0.67/ml. This gives a dimensionless value, s₀ = 6.7·10^-5. The infection source is distributed among the two phenotypes in different proportions meant to reflect parasite prevalences in the community: p _X + p _Y = 1. Dynamic simulations with stochastic S use relatively low EIR ε = .3/day, due to computational constraints, but the use of equivalent stationary (mean) sources allows sampling of a broad range of EIR. Figure 1(b) demonstrates the effect of stationary EIR on x,y-equilibria. It remains marginal for a wide range of EIR, but high values O(10 – 100) can bring about a significant shift.

Figure 2 shows computed solutions of system (3) with the parameters of Table 1 in two cases: (a) 'y-dominant' α > α₁; (b) x,y coexistence. Solid curves on each of the upper panels show x, and dashed curves y; the same marking is used for their specific (SS) effectors on the lower panels, while the non-specific (NS) is marked by thin curves. In each case three solutions are compared, marked in shades of gray: black corresponds to the unperturbed system (i.e. no inoculation), dark gray has steady sources S _X = p _X εs₀; S _Y = p _Y εs₀, and light gray has random sources of the same proportions. In case (a) inoculation has a marginal effect on the 'dominant' species, as its equilibrium level O(10^-1) stays far above the source level O(10^-4), but the 'losing side' gains in strength at a level comparable to the effective stationary source (Figure 2(a): 2 gray 'x-curves'). Indeed, the complete clearance of x by day 30 in plot (a) (black 'x-curve') is changed into the 'quasi-equilibrated' value x* ≈ 10^-2y*. Note also that random inoculation paths closely follow the 'mean inoculation' curves. This observation is used to replace (the computationally extensive) 'stochastic source' with its stationary mean.

Figure 2

Dynamic histories. Three long-term dynamic histories for 2 cases of Figure 1 over the time range 0 <t < 250 days (horizontal axis). The top plot in each panel shows dimensionless densities of two strains (x(t) – solid, y(t) – dashed), relative to pyrogenic levels x _p = y _p = 1; the bottom plot shows their immune effectors, I(t) – thin, J(t) – solid thick, K(t) – dashed, weighted by their respective efficiencies e _n ;e _X ;e _Y . For each panel, the shading (black, heavy gray, light gray) corresponds to the case of 'no-source', 'stationary source', 'stochastic source'. Note that 'y-curves' (dashed) are essentially the same in all 3 cases, leaving y unaffected, while x is strongly affected by the 'stationary' and 'random' sources.

For coexistence, in Figure 2(b) the equilibrium values show little change, but the amplitude and phases of 'damped oscillations' are shifted during the initial stage. Overall the dynamic simulations confirm the earlier conclusion that inoculation (stationary or random) will change the outcome of competition to a coexistence pattern, though the two equilibria may differ by orders of magnitude.

Single drug treatment

The simple growth-treatment model (2) for a single parasite without immune response allows exact analytic solution, expressed through the multiplier function $X\left(t\right)=a\left[t-\frac{b}{\beta p}\ln \left(\frac{1+D_0^p}{1+D_0^p{e}^{-p\beta t}}\right)\right]$,

$x\left(t\right)=e^{-X\left(t\right)}\left[x_0+{\displaystyle {\int }_0^t{e}^{X\left(t-\tau \right)}S\left(\tau \right)d\tau }\right]\left(7\right)$

This form allows one to estimate the drug efficacy (maximal burden reduction) and the time required to reach it. Assuming zero inoculation (S = 0) in (2) we get

$\begin{array}{l}{t}_{\min }=\frac{1}{\beta }ln\left[{\left(b-1\right)}^{1/p}{D}_0\right];\hfill \\ X_{\min }=\frac{a}{\beta p}\left\{\ln \left[\frac{b^b}{{\left(b-1\right)}^{b-1}}\right]+p\ln D_0-b\ln \left(1+D_0^p\right)\right\}\hfill \end{array}\left(8\right)$

in terms of its relative clearing b > 1, Hill exponent p, drug decay-rate β, and (relative) initial dose D₀. The (dimensionless) clearing level for malaria is set at x _c = 10^-11 (relative to the pyrogenic x _p = 1), i.e. less than one parasite in the entire 4.5-liter blood volume of an adult. Figure 3 shows (upper left) the maximal parasite removal as a function of the relative initial dose D₀ for three hypothetical drugs with decay rates β = 0.2/day, 0.1/day, and 0.05/day; (upper right) the corresponding clearing time t_min; (lower left) the resulting solutions x(t) for several initial doses D₀. The dashed line at the bottom is the hypothetical (relative) clearing level for malaria parasites. In particular, the decay rate 0.05/day (close to mefloquine) would require an initial dose d₀ = 2d₅₀ to kill the parasites by day 18. The lower right plot compares the effect of a drug with β = 0.05/day on two hypothetical strains: 'sensitive' (thin line) with clearing rate b _Y = 3 and relative initial dose D _Y = 4, and 'resistant' (thick line) with b _X = 1.5, D _X = 2. So both parameters (parasite clearing rates and sensitivity thresholds) differ by factor 2. Note that the different relative D's do not imply different actual doses, as both strains inhabit the same 'treated host', but rather refer to different sensitivity thresholds d₅₀ for x and y.

Figure 3

Pharmaco-dynamics of the simple 'growth-treatment' model. The upper-left panel shows maximal parasite reduction (8) as a function of (relative) initial dose D₀ for fixed clearing efficiency b = 3 and several decay rates β. The upper right panel shows the corresponding time to attain 'maximal clearing' (lowest parasite density). The lower left panels plot exact solutions x(t) for fixed b = 3, three choices D₀ = 2,3,4, and two decay rates β = .2 (thin curves) and β = .1 (thick curves). Only one of them (β = .1; D₀ = 4) leads to complete clearing by day 12. The lower right panel compares treatment outcomes of a sensitive strain (thin) and resistant strain (thick) at the decay rate β = .05, with clearing rates b _X = .5b _Y ; D _X = .5D _Y .

This method also works in the case of zero parasite growth (a = 0), when equation (2) is solved for initial value x₀. The corresponding solution x(t) (7) would then stabilize at its minimal value, $x_{\min }/x_0={\left(1+D_0^p\right)}^{-b/\beta p}$. So the 'clearing dose' to bring parasitaemia down from x₀ below x _c is estimated by

$D_0\ge {\left[{\left(x_0/x_c\right)}^{\beta p/b}-1\right]}^{1/p}\left(9\right)$

While immune regulation complicates the dynamics of drug treatment beyond the simple model (2), it still affords some clues for possible treatment strategies. Indeed, sustained immune levels (I*; J* > 0) in the 'mixed' system (3) or a simpler 'single-strain' version would effectively lower the parasite growth rate from its natural value a to a' = a(1 - c _n I* - c _X J*). This would automatically raise the clearing drug-efficiency for such an 'immune-competent' host from its 'natural' value b to b' = B/a', hence to higher clearing in a shorter time, as illustrated in Figure 3 (lower left panel). This result supports earlier findings [9, 10] on the positive impact of acquired immunity on antimalarial therapies.

Unlike this 'simple growth' (2), however, dynamic immune regulation (with or without treatment) renders system (3) unsolvable, so it is analysed using numeric simulations. The analysis of treatment for mixed phenotypes starts with a single drug dose applied at the peak parasitaemia, close to the pyrogenic threshold, and considers two phenotypes with the parameters of Figure 4(a), where y dominates and, in the absence of inoculation, drives x to extinction. As above we assume different 'drug clearing rates' and sensitivity thresholds for x,y: b _X = 1.5;b _Y = 3; d _X = 2d _Y (so the x- strain exhibits substantially higher resistance than y). The drug is introduced by day 10, when the 'fast' y reaches its peak parasitaemia. The outcomes are shown in Figure 4 for two different values of the relative initial dose D _Y = d₀/d _Y . For D _Y = 1.5 (left plot) both x and y are brought to their minima by day 20 and 28 respectively, with x taking a temporal lead through day 45. Then the competitive dominance of y is restored, and, as in Figure 4(a), leads to the demise of x by day 70. Increasing the dose to D _Y = 1.85 is sufficient to eliminate y, making x the sole survivor, by day 27. Note that D _Y = 1.5 or 1.85 exceeds the predicted value D₀ = 1.37, if one applies a 'neutral-growth' model (2) initiated at the pyrogenic level x₀ = 1. This example demonstrates important differences between the 'balanced states' of nonlinear interacting parasite-immune-drug systems and the simple linear-growth system (2) often used in pharmacokinetic estimates of dosage.

Figure 4

Single treatment. The effect of single treatment (shaded region at the bottom) on the x-y competing pair of Figure 2(a) with different initial dose D _Y = D₀ = 1.5 vs. 1.85, over an 80-day time span. The upper plots show (as above) the rescaled x,y densities (x-solid, y-dashed), the lower plots the corresponding immune effectors (J,K – solid/dashed, I – thin line).

Thus the above has shown how a high single-treatment dose can lead to the dominance of an otherwise less-fit drug-resistant/tolerant phenotype, in the absence of inoculation. Clearly, the introduction of a y-inoculum following such treatment would restore 'y-domination,' as predicted in Figure 4(a). Successful malaria prophylaxis or treatment in endemically-exposed hosts typically requires repeated doses.

Repeated periodic treatment

Now the effect of multiple, periodically-spaced drug doses at time intervals T is explored. The drug concentration and its 'clearing function' become periodic as well, namely,

$D\left(t\right)=D_0\frac{e^{-\beta \text{Mod}\left(t,T\right)}}{1-e^{-\beta T}}\left(10\right)$

where Mod(t,T) designates a periodic linear function t on interval [0,T]. Turning to the clearing function φ(D(t)), observe that a high Hill exponent of φ(z) allows φ(D(t)) to be approximated by a periodic step-function, taking value 1 on interval [0, T _M ] and 0 on the complimentary range [T _M ,T], with T _M ≈ ln D₀/β. Hence as the first-order approximation the periodic φ(D(t)) can be replaced by its mean value,

$\overline{\phi }=\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}\phi }\left(D_0{e}^{-\beta t}\right)dt\approx \frac{\ln D_0}{\beta T}\left(11\right)$

where D₀ is either D _X or D _Y . The long-term effect of such treatment is to effectively lower the parasites' growth-rates to a _Y - b _Y $\overline{\phi }$ and a _X - b _X $\overline{\phi }$. That in turn can modify (reduce) the basic relative growth-rate parameter $\alpha =\frac{a_Y-b_Y\overline{\phi }}{a_X-b_X\overline{\phi }}$ used in the qualitative analysis of equilibria (see Figure 5). Based on the comparison with the stationary (mean) case, it is to be expected that α <α₁ should change y-domination to a (periodically-modulated) coexistence pattern, while α <α₂ would bring about 'x- domination'.

Figure 5

Effect of periodic treatment on inter-strain competition. The outcomes vary from long periods of 'y-domination' (top), to a coexistence pattern (middle) to persistent 'x-domination' (bottom). The upper plot has relatively short windows of 'x-domination', the lower one has persistent 'dominant x' by 1 or 2 orders of magnitude over y. The treatment cover is shown as gray areas; curve marking and axes are the same as Figure 4.

Therefore approximate formula (11) is applied to get crude estimates of the treatment frequency for the establishment of the resistant strain. Assuming for simplicity b _X = 0 (fully-resistant strain), the equation is obtained for two 'critical' treatment frequencies/periods, T _C (for coexistence) and T _R (for 'x-domination'),

$T_C=\frac{b\ln D_0}{\beta \left(\alpha -{\alpha }_1\right)}>T_R=\frac{b\ln D_0}{\beta \left(\alpha -{\alpha }_2\right)};\left(12\right)$

in terms of pharmacokinetic parameters β and $b=\frac{B}{a_X}$, initial dose D₀, and relative growth rate $\alpha =\frac{a_Y}{a_X}$. As in section 1 the sensitive strain y is assumed competitively superior in the absence of drug, i.e. α > α₁, and we expect the following long-term outcomes:

(i) treatment frequency T > T _C will maintain y-domination;

(ii) T _C > T > T _R will bring about coexistence, with alternating x,y phases;

(iii) T <T _R will lead to dominance by the resistant x-strain.

To test these predictions the y-dominant case of Figure 2(a) is subjected to random inoculation with EIR = 0.3/day, distributed among 2 strains in proportions p _X = .6; p _Y = .4 (60% of inoculations are x-type, 40% are y). Periodic treatment (10) is then applied to such a host, and a range of values of period T sampled, using the pharmacokinetic parameters of the two strains as in the previous section, b _X = 1.5;b _Y = 3; d _x = 2d _Y .

Figure 5 shows the resulting outcomes, ranging from y-dominance at period T = 120 days, to alternating x or y dominance (T = 45 days; 'periodic coexistence'), to x-dominance at T = 20. The inoculations, though barely visible in plots, still play a role here. Indeed, without them the process would terminate (at T = 120) after the first two cycles. Namely, the first cycle would drop 'treated y' to near-extinction, after which it would rebound, at a lower level, and drive x to extinction (Figure 5, left), but a second treatment would y altogether. These numeric results confirm the earlier qualitative predictions on the 'dominance – coexistence' transition, but the 'critical' periods T _C ;T _R differ from the simple estimates (12).

The effect of treatment on mean parasite densities is demonstrated by taking a 'sensitive + fully resistant' pair (with a different set of parameters) and examining the period-average values of parasitaemia, $\overline{x}\left(T\right)=\frac{1}{T}{\displaystyle {\int }_{nT}^{\left(n+1\right)T}x\left(t\right)dt}$, $\overline{y}\left(T\right)=\frac{1}{T}{\displaystyle {\int }_{nT}^{\left(n+1\right)T}y\left(t\right)dt}$ and their dependence on T. The relative initial dose 2.5 <D _Y < 6 is also varied. In each plot in Figure 6, observe that high treatment-frequency (short T) drives y to extinction and leaves the dominant x at (or close to) its equilibrium value. But an increased initial dose D _Y = D₀ extends the frequency range for x-domination, from T ≈ 42 days at D₀ = 3.5, to T = 60 at D₀ = 6. It can be shown that the competition model (3) with periodic coefficients (due e.g. to 'periodic treatment') has a stable 'periodic equilibrium' – a counterpart of the stable 'stationary' equilibrium.

Figure 6

Time average densities. Periodically treated mean values $\overline{x}$, $\overline{y}$ (thick/thin) as functions of period 35 <T < 120 for increased initial dose 1.5 <D _Y < 3.5. Increased D _Y and frequency 1/T facilitate the onset of drug-resistance (dominant $\overline{x}$), e.g. the transition period T _R = 30 for D _Y = 1.5, is shifted towards T _R = 49 for D _Y = 3.5. Dashed lines indicate untreated equilibria for y (thin) and x (thick). The parameter α = 2.06 > α₁ assures 'y-domination' for the untreated case.

The above shows that persistent treatment creates a new 'effective environment' that could significantly alter the time-averaged distribution of phenotype densities within a host. Provided that such 'over-treated' hosts make up a sizable fraction of the population and that transmission to mosquitoes is in some way proportional to parasite densities, this could further promote the spread of drug resistance through the community.

Stationary treatment

The above has demonstrated some critical effects of periodic drug treatment and superinfection on immune-mediated competition. Many applications require a more 'exact' account of the relationships among all these factors, however. Specifically, this section examines the effects of (i) drug efficiencies b _X ≪ b _Y , (ii) treatment intensity, as measured by the (period average) factor $\overline{\phi }$ (11), and (iii) superinfection sources {S _X , S _Y }, through analysis of an (approximate) stationary system obtained by averaging the sources and the treatment regimen. The resulting stationary model allows all three factors to be incorporated in a simple and efficient way.

The 'averaging' proceeds by replacing the stochastic inoculation source by its steady (mean) value (S _X , S _Y ) = (p _X , p _Y ) $\overline{S}$, where $\overline{S}$ is the product of EIR and 'mean inocula' (merozoite release s₀), and (p _X , p _Y ) – the relative fractions of the two phenotypes. By the same pattern we replace the 'periodic drug function' φ [D(t)] with its mean value (treatment intensity) (11), $\overline{\phi }$ = $\overline{\phi }$(T,D₀); 0 ≤ $\overline{\phi }$ ≤ 1, considered as a function of the initial dose D _Y = D₀ and treatment period T. The resulting 'average' equilibrium system is similar to (4),

$\begin{array}{l}x\left[\left(\overline{\phi }-\frac{1}{b_X}\right)+\frac{1}{b_X}\left(\frac{x}{m_1}+\frac{y}{n_1}\right)\right]=S_1/b_X;\\ y\left[\left(\overline{\phi }-\frac{1}{b_Y}\right)+\frac{1}{b_Y}\left(\frac{x}{m_2}+\frac{y}{n_2}\right)\right]=S_2/b_Y\left(13\right)\end{array}$

But somewhat different rescaling and notations are used for the 'x,y'-intercepts (designated here by {m _i , n _i }), and for the stationary (mean) sources (S₁,S₂) = (αp _X , p _Y )S of relative strength S = EIR × s₀/a _Y , where α = a _Y /a _x – the above 'fitness' parameter. All three factors (intensity $\overline{\phi }$, clearing rates b _x ≪ b _Y , and sources {S _i }) enter (13) explicitly. The resulting 'stationary equilibria' of (13) are expected to approximate the stable periodic (or quasi-periodic/stochastic) equilibria of the original system.

The exact solution of the 4-th order algebraic system (13) is given by (grossly cumbersome) functions: x = x* ($\overline{\phi }$, S,b _X ,b _Y ) and y = y* (φ,S,b _X ,b _Y ). (Note that these can be brought into a more manageable 'analytic form' under simplifying assumptions, e.g. a 'highly sensitive' y-strain.) These can easily be manipulated, however (for numeric and graphic purposes), by any symbolic algebra package (here Mathematica 5). Having four independent parameters in x*, y* some of them can be fixed (e.g. clearing rates b _Y ;b _X – for the drug-sensitive/resistant phenotypes), to examine the effect of the treatment intensity 0 <$\overline{\phi }$ < 1, and/or the source strength S.

Figure 7(a) shows equilibria x*,y* as functions of treatment intensity 0 <$\overline{\phi }$ < 1, and Figure 7(b) does the same in terms of the source strength 0 <S < 3. In the upper plots (a) we fix b _Y = 10 (for the drug-sensitive strain), and take three dispersed values of the source strength S. The lower plots (b) show a similar exercise for x*,y* as functions of source strength S for three treatment intensities. Three solid curves on each panel correspond to three choices of b _X : highly drug-resistant b _X = .25 (top solid curve), moderate b _X = 1 (middle curve), and relatively sensitive b _X = 4 (bottom curve). The overall effect of increased intensity $\overline{\phi }$ is to bring down y* (by several orders of magnitude), while maintaining near-stable x*. In the upper panel (a), x becomes dominant at relatively low treatment levels; this threshold $\overline{\phi }$ _c increases with source strength (0 <S < 3), but remains confined within a narrow range .1 <$\overline{\phi }$ _c < .2. In the lower panel (b), equilibria increase with the source strength S (as expected), and, at φ = .5, the x-strain becomes fully dominant. Note that in both panels widely divergent drug sensitivities b _X of the x-strain have only marginal effects on the behaviour of y-equilibria (either its $\overline{\phi }$ – or S-dependence).

Figure 7

Equilibria of stationary model. Three plots of panel (a) show equilibria of system (13) (x-solid, y-dashed) as functions of relative treatment intensity 0 <φ < 1, for selected values of the superinfection source S. Panel (b) does the same in terms of the variable 0 <S < 3 (source strength) for 3 selected treatment intensities. Three solid curves on each plot (in descending order) correspond to 3 values of clearing rate for resistant strain: b _X = 0.25 (highly resistant), 1 (medium); 4 (low).

This suggests that treatment intensity is more crucial for the onset of resistance than inoculation frequency. The increased source by itself (at all treatment levels) would raise equilibria but diminish their relative difference: x*/y* ≈ 1 at high S.

Discussion

The above has examined how immune-mediated competition between parasites is perturbed by a drug to which the competing parasites are differentially resistant, and how drug dose, drug timing, and inoculations of new parasites affect these interactions. For clarity, the examples focused on the situation in which, in the absence of drug, the drug-sensitive strain is numerically dominant in immune-mediated competition, but the model spans the full range of possibilities. That is, while the outcome of immune-mediated competition is either that one phenotype dominates, or the phenotypes coexist, further possibilities arise when phenotypes of differing drug sensitivity are subjected to treatment and superinfection. Figure 7 illustrates possible outcomes for a range of sensitivities, drug efficacies, treatment intensities, and superinfection frequencies, and shows the potential for comparing model predictions across the range of possibilities that may arise in empirical studies.

Actual frequencies of parasite inoculation vary with the size of the vector mosquito population, and, like the composition of inocula, with the prevalence and characteristics of parasites and immunity in the human population. The baseline results reflect an inoculum that initiates an infection, i.e. after any previous infections have cleared. The parasites in superinfecting inocula are in general greatly outnumbered by those in an ongoing infection, particularly when the ongoing infection is at or near its peak. In the model, parasite inoculation enters either as a random (stochastic) source, or as its stationary (mean) value. The latter allowed exploration of a broad range of source strengths (EIR): overall, increased strength drives both equilibrium densities (or their 'period means') up, but increased intensity and efficacy of treatment eventually tips the outcome from the strain that is more 'fit' in terms of a host immune response to the one that is more 'fit' in terms of a drug.

What is transmitted from an infected human to a mosquito, and onward, poses another complex set of questions, but our results here seem to support the hypothesis that immune-mediated interactions can shape the spread of drug resistance, even if the phenotypic traits are not linked genetically [19, 20]. The fitness cost of drug resistance, considered here simply in terms of replication rate, is likely to be multifactorial in a population of hosts heterogeneous with respect to infection histories and immune profiles [21]. It would be interesting, in future work, to explicitly consider the genetics of parasite drug resistance, with respect to both origin and spread, and the common use of anti-malarial drugs to "treat" non-malarial fevers or other symptoms.

Antigenic variation can also play a role in parasite competition. From the standpoint of this model the only potentially relevant effect is a change of immune stimulation/clearing, since parasites' intrinsic replication or drug sensitivities should not be affected. A simple way to accommodate antigenic variation (in lieu of more complicated 'multi-strain/multi-clone' models [22, 23]) is to make the specific immune clearing function to change in time c(t). A drop in c can be thought of as resulting from a new variant of 'low cross-reactivity' to prior antibodies, taking over and growing into a dominant strain. As the clone keeps replicating, specific antibodies develop, and function c(t) gradually 'relaxes' to its normal (relatively high) value. Such a 'random drop + relaxation' form of c(t) was proposed in [13].

Time dependent c(t) would make the model non-stationary, the same way as the 'variable inoculation source' or 'non-stationary treatment'. In each case the strategy here was to 'average variable coefficients' over time. Applying the same 'averaging methodology' to variable clearing rates would replace them with somewhat lower 'mean values' {$\overline{c}$ _X ; $\overline{c}$ _Y }, and the new 'effective' c would enter formulae and analyses. Indeed, lower {e _X ;e _Y } (proportional to {$\overline{c}$ _X ;$\overline{c}$ _Y }) would change (decrease) 'fitness thresholds' {α₁,α₂}, and drive up equilibrium levels of x,y (as both are cleared less 'efficiently' on average).

Thus the implication of antigenic variation is that degrees of cross-reactivity between parasites may fluctuate during the course of an infection, and its net effect would be to (effectively) lower the immune efficiencies, and change the related equilibria, the fitness ranges of 'domination and coexistence' etc. The major points of interest here have to do with the relationship between the parasite entities, however, which would shift if antigenic variation differs between them in some systematic, biased way, which would then be expressed in terms of relative competitive advantage. Thus most of our results would maintain their qualitative form, but some quantitative changes on the predicted outcomes of treatment would be expected (Figures 6, 7). It might prove interesting, in future work, to explicitly consider the extent to which antigenic variation gives rise to fevers which give rise to drug-taking, for instance.

Apparently no previous work has examined these critical within-host interactions between parasite, immune-response and drug dynamics in a malaria infection, but several models have touched on important aspects of these analyses or on closely-related issues. Davis and Martin [24] compared several simple descriptive models of parasite clearance dynamics during curative drug treatment. Based on the results of pharmacokinetic-pharmacodynamic models, Hoshen et al. [25] argued that well-timed follow-up doses might eliminate even resistant phenotypes, and Simpson et al. [15] that, to prevent resistance, larger doses should be deployed as standard at the first introduction. While these models did not consider host immune responses or parasite replication, Austin et al. [14] noted that "most drugs act best against replicating pathogens in combination with effective immunological responses," joined a pharmacokinetic model to a simple model of pathogen-immune dynamics [26], and derived drug doses and frequencies necessary to reduce the average lifespan of infected RBCs below a critical threshold. Their model did not consider drug- or immune-mediated competition between strains, however.

Gatton et al. [27] modeled the risk of a drug-resistant mutant arising during an infection, in terms of the rate at which a single parasite genotype switches between antigenic variants, the response rate of the corresponding specific antibody, and a simple time-of-treatment model. They found mutants most likely to arise in hosts lacking any such specific antibodies, in hosts treated before antibody response to a switch of variants, and, in recent work [28], in hosts taking drugs with long half-lives or in subcurative doses. Hastings [29] recently developed parasite-population-genetic models that represent competing forces behind shifts in equilibrium gene frequencies in terms of within-host averages and proportions; though these models do not yet encompass considerations of immune responses or other dynamic factors, he emphasized that "intense competition between separate malaria clones co-infecting the same human can generate complex dynamics," and that "the dynamics underlying the evolution of antimalarial resistance may therefore be much more complex than previously realized." This appears irrefutable, hence the model above was developed to address these "complex dynamics."