### Basic algorithm for modeling dynamics of aging populations

Individuals in a population age with some average duration

*D* before progressing to the next stage of life development or senescence, and not all individuals necessarily take the same amount of time to age. The youngest members of the population might be created at a time-varying rate

*s*. In population modeling, this is a time delay system which can be difficult to solve [

68]. A useful formalism for approximating pathogen population dynamics was adapted by Lloyd from ecology [

69,

70]. The idea is to introduce a set of coupled ordinary different equations (ODEs) that govern the time evolution of a set of of fictitious variables,

*P*
_{1},

*P*
_{2},

*P*
_{3},…

*P*
_{
N
}. The ODEs are chosen so that the rate at which individuals leave the population for the next stage (or death) is approximately a normal function of time

*t* with mean

*D* and variance

*σ*
^{2}=

*D*
^{2}
*N*
^{−1}. Ignoring environmental influences which might remove the population, the ODEs are

In this formalism, the total population is given by the sum of all *P*
_{
n
}, and *s* is the rate (source term) at which new individuals enter the population. The tangible quantities *D* and *σ* set *N*, the total number of components. (If *D*=*σ*, then the system of Equations 1 reduces to a single equation describing exponential decay.) As mentioned above, this formalism has been used to describe within-host populations of the asexual *Plasmodium* cells interacting with various systems of the host [29–31]. For this report, this formalism has been adapted to include gametocytes as well. The populations considered in this report are shown schematically in Figure 1, along with the corresponding values of *D* and *σ*. (All population sizes are stated as number per *μ*
*L* of blood.) Because the full system of differential equations for the modal formalism involves many hundreds of components, it is given in detail in Additional file 2. A qualitative description is given below for each of the major parts of the model.

### Model for asexual parasite population dynamics

As shown in Figure 1, five populations of asexual parasite cells which are morphologically distinct from each other were considered: (1) ring stage, (2) early trophozoites, (3) late trophozoites, (4) schizonts, and (5) merozoites. The dynamics of each of these populations is described by a set of equations similar to Eq. 1 with the addition of terms due to host immune responses described below. (See equation 6 in Additional file 2). The duration and variance of each of the intracellular stages (1-4) were taken to be the same, with duration *D*
_{
Asx
} = 12 hr, but with the variance
varied from simulation to simulation. (The sampling of parameters varied from simulation to simulation as explained and tabulated in Subsection “Sampling the parameter space” below). Since merozoites are short lived in blood, with duration *D*
_{
μ
}=0.1 *hr*[71], just one compartment was used for them, *μ*, and take *σ*
_{
μ
}=*D*
_{
μ
}. There are two sources for the asexual population: (1) the primary release of merozoites from the liver, thus triggering the blood infection phase of the disease, the (2) merozoite that are released from the bursting schizonts, with *p* number of merozoites released on average from each schizont. The simulated infection begins with the primary release of merozoites, an event that takes just a few hours to complete [72]. For every simulated infection, it was assumed that the rate of primary release was 0.002 (*μL*
*hr*)^{−1}, maintained for the one hour. This is equivalent to 10^{4} merozoites being released into an adult human with blood volume 5×10^{6}
*μL*. Whether from primary release or from bursting schizonts, the rate of creation for new trophozoites is *μ*
*ζ*
*V*, where *μ* is the density of merozoites, *ζ* is the binding affinity of merozoites to their target erythrocyte, and *V* is the density of the targeted erythroyctes (reticulocytes only for *P. vivax*, all red blood cells for *P. falciparam*). (Remember that *μ* is also time dependent.) For this report, asexual density *Asx* means the sum of the ring stage, early and late trophozoite, and schizont blood densities. (Although the later stages of the intracellular asexual forms sequester in *P. falciparum* infection, *Asx* was taken as a measure of the full load of the asexual forms on the host.) The parasite days for the asexuals mentioned in the discussion, *PD*
_{
Asx
}, would be the average of *Asx* during infection in host times the duration of infection in host. If
*Asx* ≥ 0.01 *μL*
^{−1}, then it is assumed that a fraction *r*=0.05 of new infected erythrocytes are committed to the cryptic sexual pathway.

Experimental and clinical data put strong constraints on the values of

*p* and

*ζ*. Define the initial reproduction rate

*R*
_{0} as the average number of progeny an individual parasite would have after primary release of merozoites and before host immune responses. One can show that

*R*
_{0} to

*p*,

*ζ*, and

*D*
_{
μ
} are closely related [

29]:

(Here *V*
_{0} is the initial density of vulnerable red blood cells.) Observations of parasite growth in neurosyphilis patients treated with malaria therapy [64, 73] as well as in inoculated volunteers [74] indicate that *R*
_{0}≈15 for both *P. vivax* and *P. falciparum*. Photographs of bursting *P. falciparum* schizonts indicate *p* ≥ 16 [75]. Thus, for this report *p* = 16, *R*
_{0}=15. Equation 1 fixes *ζ* to 3×10^{−5}
*μLhr*
^{−1} for *P. falciparum* and 2.4×10^{−3}
*μlhr*
^{−1} for *P. vivax*.

It is assumed that the four intracellular stages are attacked by the innate immune response of the host, with time-dependent clearance rate *χ*
_{
Inn
}. Host acquired responses in malaria are complicated, (see [76] and [77] for review), but for simplicity, it is assumed for this report that the main antibody response is against the schizont stage with time-dependent clearance rate *χ*
_{
Sch,Ab
}. (Immune dynamics models are explained in more detail in Subsection “Model for immune response dynamics” below.) Also for simplicity, the variation in surface antigens which *P. falciparum* employs for immune evasion was not considered [78].

### Model for population dynamics of sexual forms

As mentioned above, two very different models of gametocytogenesis were considered, the “cryptic sexual” (CS) model and the “non-cryptic sexual” (non-CS) model. As illustrated in Figure 1, in the CS model there is a set of parasite stages that morphologically resemble the asexual stages, but are committed to eventual development into overt gametocytes. The ring stage, early trophozoites, late trophozoites, schizonts, and merozoites each have cryptic sexual counterparts. (See equation 7 in Additional file 2.) An erythrocyte invaded by a cryptic sexual merozoite becomes a gametocyte. For each simulation, it was assumed that the values corresponding to *D*
_{
Asx
}, *σ*
_{
Asx
}, *N*
_{
Asx
}, *D*
_{
μ
}, *p*, and *ζ* are the same as for the asexual populations. Based on studies of gametocyte duration [50, 51]. the immature gametocyte duration *D*
_{
IG
} was set to 216hr for *P. falciparum*, with *σ*
_{
IG
} = 24 hr. For *P. vivax*, *D*
_{
IG
} was set to 72 hr with *σ*
_{
IG
} = 12 hr. For simplicity, the mature gametocyte population was represented with a single compartment,
with exponential decay, *D*
_{
MG
} = *σ*
_{
MG
} = 156hr, based on a study of gametocyte dynamics in malaria therapy patients [56]. (The actual point at which a cryptic sexual form emerges is not yet clear, but the CS model is a extreme formulation of the cryptic sexual thesis.) It was assumed that the cryptic sexual forms are subject to the same innate clearance rate *χ*
_{
Inn
} and antibody clearance rate *χ*
_{
Sc,Ab
}. In the Non-CS model, there are no cryptic sexual stages at all; a certain proportion (*r*=0.05) of erythrocytes invaded by regular merozoites directly become gametocytes once *Asx* ≥ 0.01 *μL*
^{−1} (See equation 8 in Additional file 2.)

Four immune modalities were considered for gametocytes in this report for both *P. falciparum* and *P. vivax*, and for both the CS and Non-CS models: (1) no host immune responses to gametocytes at all, (2) no antibody response to gametocytes, but an innate response with the same clearance rate *χ*
_{
Inn
} as for the asexuals, (3) an innate response to all gametocytes with clearance rate *χ*
_{
Inn
}, and an antibody response to the immature gametocytes but not to mature gametocytes (and different from the one to schizonts) with clearance rate *χ*
_{
IG,Ab
}, and (4) an innate response to all gametocytes with clearance rate *χ*
_{
Inn
}, and an antibody response on the mature but not immature gametocytes, with clearance rate *χ*
_{
MG,A
b
}. The parameters for antibody responses to gametocytes are chosen independently to the response against schizonts. See “Model for immune response dynamics” below.

### Model for red blood cell dynamics

Three populations were used to describe the red blood cells: (1) reticulocytes, the youngest of the erythrocytes, (2) mature red blood cells, and (3) senescent red blood cells ready to be removed by phagocytosis in the spleen, liver or bone marrow [79]. (See equation 9 in Additional file 2.) Based on hematological studies [80], the respective durations of the stages are taken as *D*
_{
Re
} =36 *hr* (with *σ*
_{
Re
} =6 *hr*),), *D*
_{
Ma
} =2796 *hr* (with *σ*
_{
Ma
} =148 *hr*), and *D*
_{
Se
} =48 *hr* (with *σ*
_{
Se
} =12 *hr*). As mentioned above, if *V* is the density of the erythrocyte stage or stages vulnerable to the merozoite invasion, then the total number of erythrocytes loss to infection by merozoites is *ζ*
*μ*
*V*. The basal erythrocyte density for a healthy adult was assumed to be 5×10^{6}
*μL*
^{−1}, and if the total density of uninfected erythrocytes drop to under 60% of this value, it was assumed that the host died.

The rate of production of new reticulocytes from the bone marrow,
, is the source term for the erythrocytes and has its own dynamics. (See equation 10 in Additional file 2). In a healthy human, the erythropoietic system makes new cells at a basal rate
to maintain the red blood cell density at 5×10^{6}
*μ*L
^{−1}, (≈1736 (*μ*L *hr*)^{−1}). In the event of blood loss, the erythropoietic system can boost the rate of RBC production up to
within a couple of days [80]. However, both *P. vivax* and *P. falciparum* can modulate erythropoiesis, causing a dysfunction of erythropoiesis [28] probably through several mechanisms [81–84]. To account for these effects, erythropoietic source dynamics of the model were designed to be controlled by feedback from the dynamics of the RBC populations, so if there were no dyserythropoeisis, the rate of increase
itself is proportional to *λ*
_{
E
S
}×*ζ*
*μ*
*V*. Here
. To account for dyserythropoiesis, an offset *δ*
_{
Dys
}
*ζ*
*μ*
*V* was subtracted from the growth rate in
. For simplicity, it was assumed that the factor *δ*
_{
Dys
} had no dynamics of its own although it varied form from simulation to simulation; see Subsection “Sampling the parameter space” below. A value *δ*
_{
Dys
}=0 means there is no dyserythropoiesis.

### Model for immune response dynamics

As indicated in Figure 1, host immune responses were modeled in a very phenomenological manner: the presence of some stage of the parasite triggers the actuator stage, which is self-amplifying and eventually triggers the attacker component that removes some stage of the parasite. The entire response becomes self-limiting. The equation of dynamics are modifications of Eq. 1 above which incorporate the self-amplification and self-limiting. Although a highly simplified model of real immune processes, it does capture real aspects of human immune reactions in malaria [30], and is similar to models used to describe immune responses against viruses [85]. Up to three immune responses were considered for each simulation: (1) a quickly acting, quickly decaying innate response, (2) a slow-to-start but long lasting antibody response to intracellular asexual parasites, (3) a slow-to-start but long lasting antibody response to gametocytes. (As discussed below, in some simulations gametocytes were invisible to the host immune systems.) Consider the innate response first: the actuator
is present at some background level
when the response is quiet. Various studies indicate that malarias can trigger a rapid and strong cytokine response in the host upon bursting of schizonts, which then decays on a time scale of ≈1 *h*
*r* until triggered again [22, 86–88]. Thus, in models reported here the presence of merozoites triggers the production of more actuator, which in turn triggers the attacker component that clears its targets at rate *χ*
_{
Inn
}. A buildup of
or *χ*
_{
Inn
} then limits the response. The limit on growth of the actuator is set by parameter
. The limit on growth of the attacker is set by maximum clearance rate *χ*
_{
Inn,M
x
}, which was varied from simulation to simulation. The level of merozoites that triggers the model innate response was varied from simulation to simulation; see Subsection “Sampling the parameter space” below. The values of *D* and *σ* used for the actuator and attacker components, as well as the self-amplification and self-limiting terms were chosen to emulate the dynamics seen in the cytokine responses reported in malaria patients [87, 88]. (See equation 11 in Additional file 2.)

The actuator-attack formalism was adapted to model host antibody responses as well, since antibody responses also have feedback [89]. Acquired immunity responses take some days to activate once the triggering antigen is sensed, so a delay stage was included in the antibody models with its own *D* and *σ*. The values of the *D* and *σ* used for actuator and delay stages (as shown in Figure 1) as well as the self-amplification and self-limiting terms were chosen to emulate the T-cell dynamics seen acquired immunity responses [85]. (See equation 12 in Additional file 2.) Every simulation incorporated an antibody response to the schizont stage of the asexuals. As discussed above, some simulations incorporated an additional antibody response to the gametocytes, either against the immature forms or the mature forms. For each antibody response, the density of the targeted parasite stage, and the maximum clearance rate of the response were parameters varied from simulation to simulation. In addition, studies indicate that antibodies to malaria antigens can last from a couple of months to years [90, 91], so *D* for the attacker stage of the antibodies was varied from simulation to simulation as well; see Subsection“Sampling the parameter space” below. For simplicity, the effects of memory B cells were not considered.

### Sampling the parameter space

Several model parameters values were varied from simulation to simulation, either because the parameters vary strongly from patient to patient, or else their values are not known. For a full listing of the parameters that were varied from simulation to simulation, the range of values used, and the equations in which they appear, see Table S1 in Additional file 2. For a given class of sexual pathway, species, and immune structure, the Latin hypercube algorithm [92] was used to sample among the relevant parameters in the following manner: let
be the parameters varied for given model class. The plausible range for each *P*
_{
n
} is divided into ten equal intervals. A 10×*M* matrix M was defined in which the columns consist of a random ordering of the integers 1 through 10, with no integer repeated in a column. The integers are associated with the parameter intervals as follows: integer *k*=*M*
_{
i,n
} labels interval *k* for parameter *P*
_{
n
}. The first simulation uses values of
chosen randomly within the intervals labeled by
. The second simulation uses values chosen randomly within the intervals labeled by
. This procedure is repeated until values in the intervals labeled by
are used. The order of the integers in each column is scrambled to repeat the procedure again. Thus, for a given class of sexual pathway, species, and immune structure, 2000 randomized versions of M were used so that there would be 20000 simulations altogether. The Latin hypercube algorithm was used to attempt uniform sampling of the parameter space for a class of models, although with so many variable parameters the sampling could not be perfect.

### Solving the model equations

To solve the ODE system, a fifth-order Runge-Kutta-Fehlberg algorithm with adaptive step-size control was used for time integration [93, 94]so that the difference between the fourth- and fifth-order solutions for each component of the ODE systems was less than one part in 10^{6}. Since cells are discrete entities the following constraint was enforced: if any density of a cell population indicates less than one cell in a normal blood volume of 5×10^{6}
*μ*
*l*, then all components associated with that population are set to zero. See Additional file 2 for more details on how this was enforced. A simulation stopped when (1) the host died, (2) all parasite stages are cleared, (3) simulated time reached 3 years, or (4), the adaptive time stepping algorithm itself could not converge with the desired precision. With the choice of parameters used, only for *P. vivax* CS models did (4) occur, affecting 39 out of 8×10^{4} simulations.

### Ethical clearance

The study involved no experimental research on humans or animals.