INTRODUCTION
Gaining insights into the complexity of the interplay between antigens and host defences in malaria infection is both experimentally and ethically difficult when dealing with the human host and mainly rests on formal approaches. The current mathematical models of the dynamics of malaria humoral immunity share the common view that each antigenic variant triggers the synthesis of both a variant-specific antibody and antibodies able to cross-react with, either all variants (variant-transcending immunity) (Antia, Nowak and Anderson, 1996; Gamain, Miller and Baruch, 2001; Molineaux et al. 2001; Paget-McNicol et al. 2002) or only with a limited subset of the variant population (Recker et al. 2004). These models differ with respect to the number and the nature of the dynamic parameters chosen to account for the evolution with time of antigen and antibody populations: variable growth rates of the antigenic variants (Molineaux et al. 2001), variable switching rates of the PfEMP1-encoding var genes (Paget-McNicol et al. 2002), long-lived or transient immune responses elicited by ‘major’ and ‘minor’ epitopes (Recker et al. 2004). This multiplicity of hypotheses partly stems from the fact that the primary processes underlying the formation of antigen variants or accounting for their ability to induce variant-specific and variant-transcending (or partially cross-reactive) immune responses are unknown.
We propose that replication slippage (slipped strand mispairing taking place in replicating DNA), likely to occur within regular, low-complexity regions of the plasmodium genome, is one possible mechanism likely to stand behind the ability of the parasite to show a broad antigen polymorphism and to escape host immune defences.
Some years ago it was postulated that the variability of Plasmodium falciparum antigenic response might simply arise from aberrant but stable heteroduplexes (staggered double strands or bulge-containing duplexes), which may form upon strand slippage in replicating DNA (replication slippage) in the highly repetitive, immunogenic sequences when the DNA helix opens up and closes back during DNA replication, and expressed as polymorphic proteins by the parasite in its sporozoite and merozoite stages (Dodin, 1986). Strand slippage is recognized to be a major mutation mechanism within repetitive sequences mostly causing insertion or deletion of repeated units (Levinson and Gutman, 1987; Dover, 1995; Sinden, 2001).
We reasoned that, in addition to those regular sequences, any other nucleotide sequence, not necessarily repetitive, encoding antigenic proteins or lying close to long stretches of repeated nucleotide motifs (low-complexity genic or non-genic regions (Gardner et al. 1998, 2002a,b) should experience a variegation-like process (Saveliev et al. 2003) resulting in modulation of the expression of neighbouring genes, like switching expression of members in var multi-gene family (Paget-McNicol et al. 2002), or to mutation of genes themselves which would be subsequently expressed as new antigenic proteins. The rifin and stevor genes should also be affected, as they also stand close to the highly repetitive subtelomeric regions of P. falciparum at both ends of chromosomes 2 and 3 (Gardner et al. 1998, 2002a,b; Bowman et al. 1999; Hall et al. 2002), thus sensing the neighbouring perturbation. Telomeres are low-complexity structures in the sense they are constituted by long stretches of repeated motifs. As examples, the first 2000 bases of chromosome 2 show the 7-base tandem repeat of the left telomer. The sequence in the range 3000–9000 shows motifs that are repeated at low frequencies (every 135 and 687 bases respectively). The sequence within rep-20 shows the typical 21 base repeats together with another long regular pattern (Gardner et al. 1998, 2002a,b; Bowman et al. 1999; Hall et al. 2002). Interestingly, in this instance, up to 70% of Watson-Crick pairs would be conserved in the telomere-rep-20 region after the first slippage steps (Figure S1 in Supplementary Information) with the consequence that the corresponding heterodimers are likely to retain a significant stability. Increasing the length of rep-20 DNA through slippage would result in enhanced stability of the heteroduplexes (for a given staggering step, the longer the repeated region the larger the number of Watson-Crick pairs in the duplex) and hence, would increase the probability of switching the expression or mutating the proximal var and rifin genes. The current view of the role of telomeres in inducing the polymorphism in rifin and var genes is that they cluster into heterologous chromosome structures where active recombination takes place (Scherf, Figueiredo and Freitas-Junior, 2001). In addition or parallel to this mechanism, one may also consider slippage within a single chromosome as a general process behind the antigen response polymorphism, and taking place in the telomeres, as well as within any other low-complexity regions (see Concluding remarks).
Fig. 1. Time-evolution of erythrocyte populations with variants from one antigen in a system of 5 antigens, 5 antibodies and 5 plasmocytes. Dynamic parameters (see also text): α=1·38 per day βi,i=0·001; βi,i±1=0·00065; βi,i±2=0·000036; βi,i±3=0·000022; βi,i±4=0·00001; in per day.per antibody; γ=0·0000001 per day.per antigen; Δ=0·035 per day; ε=0·01 per day. (A) Time evolution of the of the iRBC population in host; initial populations of infected erythrocyte for the naive host infected by 1 sporozoite: 17200, 11200, 6100, 3800, 1800. (B) Time evolution of iRBC populations per μl of blood (total blood volume: 5l); The parasitaemia is detectable from 10 iRBC per μl. Initial populations as in (A). (C) Time evolution of the iRBC populations per μl in a system with no cross-reacting antibodies.: α=1·38 per day, βi,i±1=βi,i±2=βi,i±3=βi,i±4=0; γ=0·0000001 per day.per antigen; Δ=0·035 per day; ε=0·01 per day. Initial populations as in (A). (D) Time evolution of the global parasitaemia; blue graph (a), cross-reacting antibodies; green graph (b), no cross-reacting antibodies. (E) Time evolution of iRBC with initial populations respectively: 36000, 1000, 1000, 1000, 1000.
Fig. 2. Time evolution of the global RBC population of 2 variants from 2 antigens. Graphs (a): each iRBC harbours 1 variant of each antigen. Graphs (b): each iRBC has 1 variant of either antigen. Dynamic parameters: antigen 1: α=1·38 per day βi,i=0·001; βi,i±1=0·00065; in per day.per antibody; antigen 2: α=1·38 per day βi,i=0·003; βi,i±1=0·0015; in per day.per antibody. γ=0·0000001 per day.per antigen; Δ=0·035 per day; ε=0·01 per day. The initial populations of infected RBC are denoted for each antigen A and variant V as AiVj. (A) Antigen initial populations: A1V1=19000; A1V2=1000; A2V1=19000; A2V2=1000. (B) A1V1=19000; A1V2=12000; A2V1=6000; A2V2=3000. (C) A1V1=6000; A1V2=3000; A2V1=19000; A2V2=12000. (D) A1V1=12000; A1V2=8000; A2V1=12000; A2V2=8000.
Fig. S1. Percentage of Watson-Crick pairs in the rep-20 region of P. falciparum chromosome 2 conserved after strand slippage (base position: 12000–20000 on chromosome sequence).
Other DNA regions, distant from the telomeres, and encoding proteins that become inserted into the membrane of infected erythrocytes, are also potential hot spots for intense slippage and broad antigen polymorphism. For example slippage within the Erythrocyte Membrane Protein 3 gene (PfEMP3) located at position 91318–98838 on P. falciparum chromosome 2 or within EMP1 (between bases 50387–62872) on chromosome 6, retains more than 90% of the Watson-Crick pairs of the original sequences (Figure S2 in Supplementary Information)
Fig. S2. Percentage of Watson-Crick pairs in the PfEMP3 gene located at base position 91318–98838 on chromosome 2 nucleotide sequence after strand slippage.
The course of malaria invasion
In order to introduce relevant dynamic parameters into the mathematical model we shall first survey the time-course of plasmodium invasion of its host.
The onset of malaria, in its asymptomatic phase starts with the inoculation of the parasite, (Plasmodium falciparum in the most severe form of malaria in humans) during mosquito bite. The parasite enters the liver as a sporozoite and develops into merozoites, 1 sporozoite giving rise, within 12 h to 40000 merozoites (Frevert, 2004) which, after 6 days, are released into the bloodstream and invade the red blood cells (RBC). Multiplication of the merozoites in RBCs, yielding 16 new merozoites for every 1 entering the erythrocyte in 12 h (Margos et al. 2004), leads to the collapse of the infected RBC after 2 days. The subsequent release of parasites coincides with the manifestation of the fever symptoms. Most of these merozoites re-invade new RBCs. A small proportion evolves into gametocytes that will enter the mosquito vector at the occasion of a bite. The gametocytes transform within the gut of the mosquito into sporozoites that eventually migrate to the salivary glands of the insect ready for a new round of infection.
In the naïve individual (non-immunized against sporozoites) the entering parasite reaches the liver rapidly enough to prevent any significant build-up of the immune response (Frevert, 2004). In order to prevent liver invasion, vaccines directed towards sporozoites have been developed using engineered peptide antigens from the circumsporozoite protein (CS protein) (Alonso et al. 2004; Tongren et al. 2004). In the liver, the sporozoite experiences the reaction of the host defence. This response is mainly cellular and is mediated by CD8+ T-cells activated by highly polymorphic antigens from the CS protein (CSP) (and the TRAP protein which is expressed, though to a much lower concentration than CSP, in P.yoelli) (Tsuji and Zavala, 2003; Frevert, 2004). The cellular immune reaction in the liver stage is weak due to the portal vein tolerance, a mechanism that involves rapid decay of the antigen-activated CD8+ T-cells and consequently offers only limited protection against progression of the disease (Frevert, 2004).
Once in the blood-stream, the merozoites rapidly colonize the RBCs. The life-time of the free merozoite is short (10 min) and the immune system is likely not to be able to mount a response towards the free parasite. In the multiplication step within the RBCs (1 merozoite gives rise to 16–20) the parasite inserts proteins encoded by var (PfEMP1), rifin and other genes (see above and the conclusion) into the erythrocyte membrane. In the case of the var immunity, each erythrocyte in the infected RBC population harbours a specific PfEMP1 variant at its surface (Chen et al. 1998; Taylor et al. 2000). The broad repertoire of antigen variants subsequently triggers a complex humoral (mediated by antibodies) immune response.
MATERIALS AND METHODS
Materials
Plasmodium falciparum chromosome 2, 3 and 6 sequences: chromosome 2 (GenBank AE001362·2); chromosome 3 (GenBank and EMBL AL844502); chromosome 6 (GenBank NC004327).
Method
The hypotheses
We have developed a simple dynamical model (see equations below) to account for the humoral immunity conferred by activated B-cell-producing antibodies (plasmocytes) and based on the following hypotheses. (1) Strand slippage takes place during stages where the parasite multiplies (within-host, in the hepatic and erythrocyte asexual phase, and within-mosquito in the gametocyte evolution to sporozoites). We shall suppose that the slippage events leading to either mutation or gene expression modulation are essentially occurring during liver-stage multiplication and that, parasites re-infecting new RBCs will retain the same variant repertoire along the whole blood stage. This latter view allows the development of a simple formal mathematical model which would otherwise be complex when including changes in the variant repertoire after each erythrocyte invasion cycle every 2 days. (2) Strand slippage occurs at repeated, low-complexity, regions of the chromosomes which stand either within genes involved in the immune response or in their close vicinity. Slippage would generate 2 kinds of effects: (i) direct mutations within genes harbouring inner low complexity sequences and (ii) mutation of out-of gene control elements that would affect the gene expression profile (the switching of var-genes might enter this latter category).
The polymorphic epitopes arising from replication slippage, stemming either from mutations or from modulating the expression of antigen variants genes, will not be distinguished in the formal mathematical model.
In the within-host phase, the invading parasite multiplies and transforms in the liver and RBCs giving rise to a progeny where individuals have conserved ancestor antigen genes (their double-stranded DNA is more stable than any other heteroduplex arrangement) together with lower proportions of modified genes that would give rise to new antigenic determinants. These epitopes are expected to present significant degrees of similarity, as deriving from not-too-different sequences, so that the antibodies elicited by activated B-lymphocytes (plasmocytes) should cross-react with all parent determinants (or only a subgroup of them), though the strongest response would occur with its cognate variant. As an example, variants of PfEMP1 exposed at the surface of infected RBC are constituted by various combinations of the same building blocks, the DBL units where each variant triggers the synthesis of antibodies able to cross-react with any other, though with different efficacy.
The mathematical model
Saul, (Saul, 1998) and Gravenor, (Gravenor and Lloyd, 1998) have questioned the relevance of describing discontinuous phenomena (like the instantaneous collapse of the erythrocyte 2 days after infection) by the smooth solutions of differential equations. Using ‘constant rates to describe life-span is an attractive starting point’ (Gravenor and Lloyd, 1998) that, besides being appealing, may also lead to a realistic description of biological systems. For example, the actual evolution with time of the number of infected RBCs (iRBCs) can be described by a discontinuous step diagram where the iRBC population is constant for 2 days, falls instantaneously to zero and jumps back to a stable level after the released merozoites have quickly reinvaded new RBCs. The time between the burst of iRBCs and reinvasion is short (the life-time of free merozoites is 10 min) so that the period T of the collapse-reinvasion events would be roughly 2 days (2 days plus the time for reinvasion). When each infecting merozoite gives rise to 16 reinvading individuals, the iRBC population after a number, n, of periods is 16n (2 periods, 4 days, after invasion by 1 merozoite the iRBC population will be 256). Obviously the same result is obtained when the growth of the iRBC population is now described by a smooth exponential function with a time-constant ln(16)/2=1·38 per day. Hence, exponential functions (and the differential equations these functions are solutions for) are acceptable means for representing time evolution of iRBCs.
A system of 3 sets of differential equations is used to describe the dynamics of the evolution with time of variants from a single antigen and host defences (see below). When several antigens, each expressing several variants at the surface of iRBC, are considered, new interaction terms, stemming from the fact that each iRBC induces cross- reacting variants from the whole antigen population, must be introduced into the equation system (see the example below). For the sake of clarity these terms were not introduced in the formal expressions.
In the equations, the yi represents the population in RBCs harbouring on their surface the variant epitope i in a given family; xi is the population of antibody i directed against this epitope. The population xj will also cross-react with the erythrocytes harbouring epitope subpopulation yi of the same family (second term in equation (1)). In equation (2) the first term indicates that antibody xj derives from antigen yj and its population is proportional to that of activated lymphocytes (plasmocytes), zj, it originates from. The variable xj decays through a non-specific process (second term in equation (2). Equation (3) shows the instant variation of the plasmocyte populations with time. The parameter εk>0 represents the balance of growth and decay rates of plasmocytes. For εk>0, the infected RBC population will eventually go to zero with time; εk=0 is associated to a steady state that would result in persistently infected RBC population.
It is assumed that variants yi, yi+1, yi+2, yi+n are progressively more structurally distant and that the reactivity of the antibody, xi, elicited by yi towards yi, yi+1, yi+2, yi+n steadily declines. Also, the reactivities of xi towards upstream and downstream antigens yi-1 and yi+1 will be regarded as identical.
The choice of numerical values of the parameters
The dynamic parameters
The numerical values of the dynamic parameters have been chosen to match experimental observation.
The parameters αi=α=1·38 per day, is the growth rate of the infected RBC population (see discussion above).
The βi,j rate constants have been taken in a range of values from previous investigations (0·001 per day.per infected RBC) (Antia et al. 1996; Recker et al. 2004). We reasoned that the rate parameters βi,j (where i≠j) should be, somehow, related to the probability of formation of the diverse variants, which would be itself linked to the probability of occurrence of the genes they originate from. Hence the various βi,j (i≠j) could be evaluated, with respect to βi,i, from the percentage of Watson-Crick pairs computed from a formal slippage of DNA strands. In a nucleotide sequence of length L (expressed in nucleotides), with a number Nn,wc of Watson/Crick pairs in the n-slipped duplex, the ratio Nn =Nn,wc/L is linked to the relative stability energy of the duplex. For a fully paired double strand (n=0, no slippage), N0 is 1. Now, when 1 strand is slipped with respect to the other by n bases, Nn becomes less than, or equal, to 1. It is then assumed that the stabilization energy of the heteroduplex is simply linked to its number of Watson-Crick pairs. This is certainly a rough view since, (i) no mismatch energy penalties are introduced into the model, (ii) AT and GC pairs are considered as energetically equivalent and, (iii) the expressed protein reactivity is supposed to follow the order of stability of their cognate genes. Consequently, and within the scope of this hypothesis, the ratio of the population Pn of one specific staggered duplex with respect to that of the fully paired sequence, P0, would obey a Boltzmann distribution. The stabilization energies of duplexes being supposedly proportional to their number of Watson-Crick pairs, the relative stabilization energy of heteroduplex n, is (Nn−N0)/Nn:=(Nn−1)/Nn (where N0=1 as seen above) with, consequently, Pn/P0=exp[(Nn−1)/Nn].
For example, when 70% of the Watson-Crick pairs are retained in a DNA sequence after slippage, (Nn=0·7), the population of staggered duplexes should be 65% of that of the fully paired duplex. This value has been arbitrarily chosen to define, for i=i±1, the value of the dynamic parameters used in equations, β i,j=0·65 βi,i and βi,j=0·00065 per day.per RBC with βi,i=0·001 per day per amount (Antia et al. 1996; Recker et al. 2004).
The parameters γj=γ in equation (2) reflect the probability for the antigen to select the appropriate individual in a broad repertoire of B-lymphocytes. This probability is weak as antigens are not likely to directly activate B-cells without help from cytokine-stimulated auxiliary T-lymphocytes. In order to bring the parasitaemia predicted by the model to levels matching malaria field observations (where the disease can be detected from 10 infected RBC per μl of blood, Paget-McNicol et al. 2002), γ has been taken to be 1×10−7 per day.per antigen.
The second term of the right member in equations (2) represents the slow, spontaneous decay of circulating antibodies. A 3-week half-life time would correspond to a decay rate δj of 0·035 per day when assuming an exponential decrease.
As stated above, εk=ε in equations (3) represents the balance between close growth and decay rates of plasmocytes. Values of ε ranging from 0 to 0·1 per day have been used in the model.
The initial population values
The simulations have been performed considering 1 sporozoite entering the liver and giving rise to 40000 merozoites and the calculations have been limited to systems comprising 3–5 variants per antigen reacting with 3–5 antibodies elicited by plasmocytes. The initial RBC populations harbouring the diverse variant epitopes were computed either following the Boltzmann distribution described above (Fig. 1A–D), or considering that 1 population is dominant (see Fig. 1E). In the naïve individual the initial antibody populations are zero and each antibody stems from 1 given plasmocyte. The case of a patient exposed to a second parasite infection, where the antibodies and the specific memory B-cell from the naïve case would have only limited action of the new variants has not been considered.
The system of differential equations has been solved numerically using the ODE 45 routine in MATLABR..
RESULTS AND DISCUSSION
Firstly, we have simulated the humoral blood-stage immunity elicited by variants from a single antigen. Fig. 1, illustrates the example of 5 antigen variants eliciting 5 antibodies through induction of 5 plasmocytes, (i, j, k=1, 5 in equations). The general dynamic pattern of plasmodium immunity, including periodical bursts and cessation of infection, is satisfactorily accounted for by the model. Interestingly, the general course of the parasitaemia evolution with time is robust with regard to varying the antigen-specific dynamic parameters like β i,j, γ or ε (see for example the case where 2 classes of antigens with their specific parameters were considered, Fig. 2). This prediction is consistent with the common observation that the malaria infection, though multiform in intensity and duration, displays conserved characteristic features (Smith et al. 2004).
Interestingly, Fig. 1 shows that antibody cross-reactivity leads to constant, or even increased populations of infected RBCs during the course of the disease (Fig. 1B, D) whereas in the absence of cross-reaction, the intensities of the bursts of infection steadily decrease with time (Fig. 1C, D). Antibody cross-reactivity also induces a shift in time of the parasitaemia peaks and a narrowing of the period between bursts with respect to the uncoupled case. Patients would hence suffer a more chronically sustained parasitaemia, though to a lower level, than in the absence of antibody cross-reactivity (see on Fig. 1D, the first 100 days after the start of infection). These predictions are in line with Recker's model (Recker et al. 2004). We have also simulated the dynamics of the humoral immunity induced by variants of 2 different antigens. Antibodies elicited by variants in a given antigen family will, in general, not recognize variants from another antigen class. Rapid growth of combinatorial complexity is expected when increasing the number of antigens and variants (k antigens and n variants per antigen would lead to nk different infected erythrocyte families). Fig. 2 shows the evolution with time of 4 infected RBC populations (each harbouring 2 variants from 2 antigens). The time-profiles of the infected erythrocyte populations, and the global course of parasitaemia (total population of infected RBCs) are significantly dependent on the initial populations of individual variants for the same overall initial 40000 infecting merozoites.
Fig. 2 also compare the time-profiles of RBC populations harbouring 2 variants with populations bearing only 1 variant of either antigen and suggest that the parasitaemia is more severe in terms of duration and intensity when only 1 antigen is expressed, in line with the expectation, that multiplication of the number of targets for the immune system is detrimental to the parasite since antibodies directed against only 1 surface antigen variant will suffice to neutralize the infected RBC). An efficient survival strategy for plasmodium would hence be to express only 1 antigen.
Concluding remarks
The approach we have presented complements previous mathematical models with a plausible hypothesis regarding the origin of antigen variability in plasmodium parasites. The observation that ‘low complexity‘ proteins, abundant in the P. falciparum genome (Gardner et al. 1998, 2002a,b; Bowman et al. 1999) are highly polymorphic is consistent with the view that strand slippage in repetitive DNA sequences is a general perturbation mechanism. This view provides an alternative to the hypothesis that the large variability of regions close to repetitive sequences (like subtelomeric regions in P. falciparum chromosomes) stems from active exchange between these regions in several chromosomes (Scherf et al. 2001). Assessing local complexity in a genomic sequence by means of dedicated algorithms (Wootton, 1994; Dodin et al. 2000; Wan and Wootton, 2000) should permit the unravelling of new families of antigen variants. The access to fully sequenced P. falciparum chromosomes makes it possible to undertake a genome-wide survey of repetitive motifs that would confer variability to genes, either directly within the gene, or by a variegation mechanism. In this respect, it is worth noting that the gene of the erythrocyte membrane protein-3 (located far from telomeres at position 91318–98838 of chromosome 2 of P. falciparum), (see Fig. S2 in Supplementary Information) is highly repetitive and may also give rise to mutated proteins that will be inserted into the blood cell surface or that will variegate the neighbouring histidine-rich protein (109564–110580). Recognized antigens, like the liver stage antigen-3 located at position 796750–801584 on chromosome 2 and the Circumsporozoite protein (see above) at position 217997–219190 on chromosome 3 do present stretches of repeated motifs (Gardner et al. 2002a,b), as also the 41 kB antigen (281907–287045, chromosome 2). Similarly, the EMP1 gene at position 50387–62872 on chromosome 6 may also experience strong stabilization of slipped heteroduplexes.
Several proteins of unknown function, like a putative binding protein (337558–342966, chromosome 2) are also plausible candidates whenever accessible to the host immune machinery.
