The timing of investment into reproduction is a key determinant of lifetime reproductive success (fitness). Many organisms exhibit plastic, i.e., environmentally-responsive, investment strategies, raising the questions of what environmental cues trigger responses and why organisms have evolved to respond to those particular cues. For malaria parasites (Plasmodium spp.), investment into the production of specialized transmission stages (versus stages that replicate asexually within the host) is synonymous with reproductive investment and also plastic, responding to host- and parasite-derived factors. Previous theory has identified optimal plastic transmission investment strategies for the rodent malaria parasite, Plasmodium chabaudi, as a function of the time since infection, implicitly assuming that parasites have perfect information about the within-host environment and how it is changing. We extend that theory to ask which cue(s) should parasites use? Put another way, which cue(s) maximize parasite fitness, quantified as host infectiousness during acute infection? Our results show that sensing a parasite-associated cue, e.g., the abundance of infected red blood cells or transmission stages, allows parasites to achieve fitness approaching that of the optimal time-varying strategy, but only when parasites perceive the cue non-linearly, responding more sensitively to changes at low densities. However, no single cue can recreate the best time-varying strategy or allow parasites to adopt terminal investment as the infection ends, a classic expectation for reproductive investment. Sensing two cues—log-transformed infected and uninfected red blood cell abundance—enables parasites to accurately track the progression of the infection, permits terminal investment, and recovers the fitness of the optimal time-varying investment strategy. Importantly, parasites that detect two cues more efficiently exploit hosts, resulting in higher virulence compared with those sensing only one cue. However, parasites sensing two cues also experience larger fitness declines in the face of environmental and developmental fluctuations. Collectively, our results suggest that sensing non-redundant cues enables more optimal transmission investment but trades off against robustness in the face of environmental and developmental noise.

Fig 1. Schematic of P. chabaudi within-host infection dynamics and fitness optimization.

(A) Flow diagram of model tracking the dynamics of red blood cells (R; RBCs), RBCs infected with asexual parasites (I), RBCs infected with sexual— but not yet mature — parasites (), asexually-committed merozoites (M), sexually-committed merozoites (), and gametocytes (G). (B) Temporal dynamics of cue perception. Cue perception occurs one day before the production of asexual/sexual merozoites, thereby introducing a delay between cue perception and the consequences of cue perception (i.e., production of more asexual infected RBCs or gametocytes). A further delay occurs in the sexual cycle since gametocytes of P. chabaudi take two days to mature [38]. Dashed circles indicate the bursting/maturation of a particular cohort of iRBCs. (C) Modeling reaction norm evolution. We model the relationship between the perceived cue and transmission investment (i.e., the reaction norm) as a flexible spline (black lines). Specifically, a cubic spline can specify any third order polynomial at a given x value as a linear combination of the columns in its basis matrix [39] (exemplified by the red, green, and blue dashed lines; combined they produce the solid, black line). For local optimization, an arbitrary starting spline is picked (left panel) and an optimization algorithm is used to adjust the relative weights of the basis functions until a fitness maximum is achieved (going from left to right). The small coloured arrows depict how the weight of a given basis function has changed from the previous panel.

https://doi.org/10.1371/journal.pbio.3003081.g001

Cue perception mediates parasite life history and fitness

Based on previous work modeling P. chabaudi transmission investment as a function of time [23], we expect the following as the best strategy over a 20-day infection: parasites should delay investment early on, display a rapid increase in investment when sufficient gametocytes can be produced to ensure transmission, subtly reduce investment in the face of declining resources, and ultimately rise to the highest levels of investment shortly before the end of infection (roughly, “terminal investment” [7,23]). Although the model we employ here has a few key differences from that work (see Materials and Methods), we find a similar optimal transmission investment pattern as in [23] when we assume that investment is a function of time (Fig 2A). When parasites instead respond to their environment and are constrained to perceive a single within-host cue, these same patterns of transmission investment cannot be fully recapitulated. Consequently, parasite fitness is lower in every case than when the cue is time (Fig 2A, left panel).

The environmental cues that result in the highest parasite fitness include logged parasite-centric cues (e.g., log₁₀ densities of sexually-committed iRBCs). Perceiving these cues is sufficient to produce the initial delay in transmission investment, but the delay is longer than parasites sensing time as the cue, and there is no increase in investment near the end of infection (Fig 2A, right panel). Given that parasite densities rise non-linearly at the start of infections (both simulated and real [11,43–46]), sensing logged rather than raw densities results in reaction norms that are more sensitive to those early changes (Fig 2B), enabling parasites to increase transmission investment after the densities have reached a rough threshold. In contrast, sensing raw parasite densities results in reaction norms that essentially ignore changes at the low parasite densities characteristic of early infections (Fig 2B). In principle, raw and log-transformed densities encode the same information, so parasites sensing either cue should have the same optimized fitness. In practice, the log transformation amplifies the influence of low-density changes, allowing smooth spline curves to model threshold-like behaviour on a raw scale. One could also capture abrupt changes at a specific threshold by the inclusion of knots. However, optimizing the number and placement of knots is computationally challenging.

The density of host resources (RBCs; raw or logged) is also a poor cue for initial stages of infection, since they change very little early on and, therefore, do not allow parasites to engage in delayed transmission investment (Fig 2A). Overall, our results broadly suggest that early transmission delay necessitates sensitivity to changes at low parasite density.

Dual cue perception is necessary for terminal investment

We found that parasites sensing a single cue cannot reproduce the expected optimal transmission investment strategy, particularly the rapid increase in transmission investment during late infection (i.e., terminal investment). Achieving terminal investment presumably requires sensing a cue that provides a unique signal at the later stages of infection. Inspecting infection dynamics produced by parasites that sense time revealed that only RBC density provides such a distinctive late-stage signal (Fig 3A). Given the poor performance of this cue in early stages of infection, these results suggest that simultaneous perception of RBC and a parasite-centric cue could provide infection stage-specific information that enables both early reproductive restraint and terminal investment. Consistent with this expectation, in a rodent malaria experiment, including variables that reflect parasite density and host resources (specifically, replication rate and changes in RBC density) in statistical models explained substantially more variation in transmission investment than the parasite-centric variable alone [13].

To investigate the fitness consequences of sensing two cues, we optimized the transmission investment strategy of parasites sensing all pairwise combinations of log-transformed cues, using two-dimensional spline surfaces and focusing on combinations of log-transformed cues given their higher performance in the single cue optimizations. Our results suggest that sensing dual cue combinations improves parasite fitness relative to sensing a single cue (Fig 3B), with parasites sensing logged RBC and logged iRBC (asexual or total) exhibiting notably higher fitness (Fig 3B). This higher fitness is conferred by transmission investment patterns that are close to the time-varying ideal, with an initial delay in transmission investment and a late infection (nearly) terminal investment before returning to a low transmission investment (Figs 3C and S2). Increasing the spline flexibility of the single cue models (to match the degrees of freedom in dual cue models) could not reproduce this same pattern (Fig 3C), confirming that dual cue perception specifically enables terminal investment and its associated fitness benefits (S3 Fig).

To better understand the optimal patterns of transmission investment, we next analyzed the reaction norm of the dual cue combination that conferred the highest fitness: logged asexual iRBC and logged RBC. The pattern of transmission investment is defined by a looping journey through the reaction norm surface (Fig 3D), confirming that the cue combination traces a unique path over the course of an infection. We categorized the reaction norm into four distinct sections/stages (Fig 3D). Stage I corresponds to the initial phase of infection, characterized by low transmission investment, low parasite density, and high RBC density (Fig 3C and 3D). As the infection progresses, the reaction norm enters stage II, characterized by asexual iRBC production and RBC decline. During this stage, the optimal transmission investment increases considerably, followed by subtle declines as the iRBC population continues to increase (Fig 3C and 3D). The transition from a low to high transmission investment observed as stage I transitions into stage II constitutes delayed transmission investment. Following peak parasite density, the reaction norm enters stage III, which is marked by decreasing iRBC population and gradual replenishment of RBCs, signaling imminent infection clearance. At this stage, transmission investment increases to peak levels (i.e., terminal investment). Finally, as the infection nears clearance and the reaction norm enters stage IV, transmission investment sharply drops as the host heads towards recovery (Fig 3C and 3D). We note that this drop in predicted transmission investment is likely a consequence of the constraints of the spline surface, and different strategies during this stage would have negligible impact on fitness [23]. Overall, our results suggest that each cue in this combination plays a crucial role in defining progression through infection stages I to III, with asexual iRBC density exhibiting the largest variation in early and late infection (stages I and III), and logged RBC density exhibiting the largest variation during mid-infection (Stage II; Fig 3D). This temporal staggering of iRBC and RBC dynamics (i.e., their hysteretic relationship) ensures that these infection stages can be tracked. As a result, dual cue perception provides parasites with the ability to “state differentiate," to accurately monitor the progress of the infection and adjust transmission investment accordingly.

Intriguingly, the dual cue combinations that are most informative for parasites making life history decisions (top two in Fig 3B) are also the most informative for predicting patient outcomes in clinical malaria infections [47]. In both cases, two slightly out-of-phase within-host measures can uniquely define a patient’s (or simulated host’s) position on the path from health through sickness to recovery. Plotting this trajectory through disease maps (i.e., the two measures plotted against each other, as in Fig 3D) produces a looping, or hysteretic, pattern, with each infection stage occupying a unique coordinate in this two-dimensional space [47]. To further validate that forming the characteristic looping pattern is a unique feature of our identified best cue combinations, we constructed disease maps using experimental data from laboratory mice infected by a single strain of P. chabaudi [11,43–46,48] and simulated data from our models where parasite senses time. Both experimental and simulated data indicate that raw or logged RBC with logged iRBC (asexual or total in the model; only total is available in the experimental data) are the only cue combinations that generate fully separated loops (e.g., Fig 4A), while other cue combinations do not (e.g., Fig 4B and4C; full data shown in Figs S4 and S5). When early and late infection stages are “condensed" in space, simulated parasites fail to state differentiate and do not delay transmission investment or terminally invest (S2 Fig). In summary, our results demonstrate that optimal transmission investment requires parasites to unambiguously state differentiate, track infection stages, and dynamically adjust investment in response to changing conditions. Our data suggest that perceiving hysteretic cues—of which a combination of RBC and logged iRBC provides a clean delineation of different infection stages—is a viable method of attaining state differentiation.

Dual cue perception exacerbates the fitness loss incurred by environmental and developmental stochasticity

The benefits of plasticity depend on the quality of information organisms obtain from environmental cues. Adaptive phenotypic plasticity is contingent on proper phenotype–environment matching, where the perceived cue is a reliable predictor of the optimal phenotype [49–53]. Theory suggests that a weak correlation between cue(s) and the optimal phenotype (i.e., low cue reliability) could disfavour plasticity relative to other strategies (i.e., local adaptation [50,51,54], canalization [53], bet-hedging [54]). For malaria parasites, there may be unpredictable variation in the environment (e.g., RBC availability [42,55,56], immunity [42,57–59]) or in traits governing parasite development within the host [42]. Such environmental and developmental stochasticity could preclude evolution towards the strategies we predict above. Although our model allows parasites to sense components of their state, and their resulting level of transmission investment will feedback to influence that state, the results above assume that all else is equal between hosts. Specifically, the optimal reaction norms identified above are associated with a particular set of model parameters (the median estimates from posterior distributions in [42]; we refer to this as the deterministic model). This implicitly assumes that hosts are all the same or that parasites are moving between hosts sufficiently quickly that they adapt to an “average” environment.

To investigate how biologically relevant variation affects cue performance, we built a semi-stochastic model based on the work of Kamiya et al. [42], which employed a hierarchical Bayesian method to quantify inter-host variation in parasite burst size, erythropoiesis (RBC replenishment), and immune strengths. For all log-transformed single cues and dual cues, we performed 1000 20-day simulations, holding the reaction norm at the optimum from the relevant deterministic simulation, and simultaneously varying the RBC replenishment rate (ρ), parasite burst size (β), immune activation strengths ( and ) and immunity half-lives ( and ), sampling parameter values from the reported posterior distributions [42] (S6 Fig). Each choice of parameter set will influence the within-host infection dynamics, potentially altering the trajectory of the perceived cue(s), resulting in different, expressed patterns of transmission investment and, thus, fitness. Here we are asking, if parasites are “optimized” to the average host, does environmental and developmental variation change inferences about which cues are best? We quantified the effect of variation on cue performance using the geometric mean fitness, commonly used to assess fitness in fluctuating environments [60,61].

We find that the geometric mean fitness of parasites responding to a given cue, in the face of this variation, is largely uncoupled from fitness in the deterministic model (i.e., the environment in which the strategy was optimized; Fig 5A), emphasizing that the “best” reaction norm for a certain environment does not always confer the best phenotype in another. Interestingly, all dual cues have a lower geometric mean fitness than the best-performing single cue (I log, Fig 5A). Furthermore, this result emphasizes the challenge of using time as an environmental proxy, since its geometric mean fitness in the semi-stochastic model is among the bottom half of cues. Allowing only one parameter to vary similarly tends to reduce the fitness of the best dual cue strategies disproportionately (S7 Fig). These results indicate a potential tradeoff between fine-tuning fitness in a given environment and susceptibility towards environmental and developmental stochasticity.

While sensing two cues provides more information about “where" an infection is, it could also increase the amount of perceived noise and exacerbate the negative effects of environmental/developmental stochasticity. To test whether increased noise perception reduces the fitness of dual cue-sensing parasites, we compared the trajectories of the most “robust" single cue (I log; Fig 5B left panel) and the dual cue combination with the highest deterministic model fitness (R log & I log; Fig 5B right panel). As expected, the variation in asexual iRBC density across our simulations is much lower than that of both asexual iRBC and RBC (Fig 5B top panels). While parasites sensing both types of cues delay transmission, those sensing logged RBC and iRBC display highly variable transmission investment dynamics (Fig 5B bottom right). In contrast, the transmission investment pattern of parasites sensing only logged asexual iRBC resembles that of the deterministic model (Fig 5B bottom left). These results suggest that sensing certain dual cue combinations can reduce parasite fitness by increasing noise perception that induces suboptimal transmission investment.

Dual cue perception leads to higher virulence

Given that cue perception mediates infection dynamics, it could also influence the severity of infections, or virulence. We estimated virulence using the methods described in Torres et al. [47], where infection dynamics are mapped in disease space, with the y-axis quantifying host health (RBC density) and the x-axis quantifying parasite load (total iRBC density). A more virulent infection is expected to produce a higher parasite load while consuming more host resources, tracing a larger loop through disease space. To evaluate the effects of cue perception on virulence, we generated 30-day disease maps of infections, allowing sufficient time for the infection to fully resolve in our model. We use the area enclosed by each disease map as a proxy for virulence, with larger areas indicating higher virulence. We note that we are using “virulence” as synonymous with infection severity and not host mortality, thus we assume no link between virulence and the duration of infection. Recent modeling suggests that resource limitation imposes a more important constraint on the evolution of transmission investment in P. chabaudi than does host mortality due to infection [62].

Parasites perceiving only a single cue trace smaller disease maps than parasites that adopt the ideal pattern of transmission investment (by sensing time; Fig 6A). Additionally, the area of the disease map is negatively correlated with parasite fitness across all single cues (Fig 6B), suggesting that if parasites are constrained to perceive a single within-host cue, strategies that result in lower virulence would be favoured. Parasites sensing the top two dual cue combinations produce disease maps that are similar to the map traced by parasites adopting the time-varying strategy (Fig 6C) and result in higher virulence relative to parasites sensing a single cue (Fig 6D; compare x-axis values of points in B and D). Overall, these findings suggest that if parasites have limited information or have evolved to respond to only one cue (e.g., because of physiological constraints or high environmental/developmental variation), virulence will be lower, while if parasites are favoured to more finely track infection stages through dual cue perception, they will exploit host resources to a greater extent and infections will be more virulent.

Within-host model of P. chabaudi infection dynamics

We modified and extended a previously published mathematical model describing the within-host infection dynamics of the rodent malaria parasite, P. chabaudi [23] (Fig 1). The original model consists of ordinary and delay differential equations representing the interactions between the parasite, host red blood cells (RBCs), and host immunity over a 20-day infection. There are three key differences in our model from that previous work. First, we include the activity of two forms of immunity to better capture host defense mechanisms: a targeted immune response that only kills infected RBCs (iRBCs) and an indiscriminate response that kills any RBC, as in [42]. Second, we revised the model to account for the observation that sexually-committed parasites do not engage in gametocytogenesis until the next infection cycle (Next Cycle Conversion; [40]). This necessitated separately tracking the production of sexual merozoites () and asexual merozoites (M) from infected cells and consequently introduced a delay between sexual commitment and the outcome of sexual commitment (i.e., production of an iRBC that will produce a gametocyte, ). Third, we modeled transmission investment as a function of one or two within-host cues, instead of time as in [23], to better capture the mechanism of phenotypic plasticity, since there is no evidence malaria parasites track time across multiple generations. We further describe the implications of these changes on the model below.

We depict an immunologically naive host given the availability of experimentally derived parameters. Throughout an infection, we assume indiscriminate immunity eliminates all RBCs (R, I, I_g) whereas targeted immunity eliminates only iRBC (I, I_g). As in [42], we define N and W as the proportion of cells cleared per day by indiscriminate and targeted responses, respectively. The magnitude of this activity is dependent on the total abundance of iRBC ():

(1)(2)

where and denote the immune activation strengths, and represent the immune activity half-lives, and ι represents the theoretical maximum total iRBC density, which determines the threshold at which immunity saturates. Experimental data from mouse infections indicate that indiscriminate immunity is quickly activated but short-lived, whereas targeted immunity has a longer activation period but a longer half-life [42]. Because N and W have been defined as proportions of cells cleared, they can be rewritten to reflect daily rates of clearance. For example, , where x_N is equal to the rate of clearance on a given day, and is defined as , with some rearranging of terms. A similar expression provides the daily clearance rate for targeted killing, .

Mimicking an experimental P. chabaudi infection, the host receives an initial dose of asexually-committed iRBC (I₀) to start the infection. We assume no parasites in this initial dose commit to gametocytogenesis and all produce asexual merozoites (M). Therefore, before the end of the first asexual development period (, where α is the development period of one day for this species), the dynamics of asexual iRBC (I) in the host are determined solely by the invasion of RBC by asexual merozoites, natural RBC death, maturation of the initial cohort of infected cells, and immune-mediated killing such that

(3)

The first term, , represents the production of iRBC through asexual merozoite invasion of RBCs, which occurs at a rate of ε per day, given contact. Infected RBCs die at a background rate of μ, while –ln(1–N(t)) and –ln(1–W(t)) describe their rate of death due to indiscriminate and targeted immunity, respectively. The infection-age structure of the initial cohort of asexual iRBCs, denoted as Beta(sp,sp), is modeled using a Beta distribution with a shape parameter sp = 1, allowing iRBC maturation to occur with a uniform probability throughout the development period [90]. The uniform probability of maturation for the initial cohort is computationally efficient and provided extremely similar results to more narrow age distributions in a previous modeling study [23]. The rate of asexual iRBC maturation depends on an exponential survival function S, which describes the proportion of iRBC that survives until maturation. S is defined by a cumulative hazard function that is calculated by integrating the total rate of death across the development period:

(4)

During the period , asexual merozoite production is described by

(5)

where the first term represents the production of asexually-committed merozoites by asexual iRBC that have survived the duration of the developmental cycle and burst to produce β merozoites each. Asexual merozoites are lost by natural death at a rate of and are removed by their invasion of susceptible RBC.

After all the initially injected iRBC have matured or died (), the next cohort of asexual iRBC starts to reach maturity. Transmission investment is defined as the proportion, c, of this new generation that commits to gametocytogenesis, resulting in the production of β sexually-committed merozoites (M_g) per cell. These sexual merozoites invade susceptible RBCs to produce sexual iRBCs (I_g). We assume that sexually-committed merozoites have the same death rate () and rate of RBC invasion (ε) as asexual merozoites. The remaining proportion (1–c) of iRBC produces asexual merozoites and continues the asexual replication cycle. Altogether, when , the dynamics of iRBC and merozoites are as follows:

(6)(7)(8)

The parameter governing transmission investment (c) represents the proportion of iRBCs that produce sexually-committed merozoites and is further defined in the next section.

Before the first cohort of sexual iRBC matures (), the dynamics of sexual iRBC is governed solely by the rate of RBC invasion by sexually-committed merozoites, natural RBC death, and immunity:

(9)

When the first cohort of sexual iRBC matures (), each sexual iRBC produces a single gametocyte, G, such that

(10)(11)

where is the background death rate of gametocytes. The survival probability of sexual iRBCs during their development cycle is defined by the exponential survival function S_g, which is similar to that of asexual iRBC (S), except the cumulative hazard function is integrated over the length of sexual iRBC development, :

(12)

Throughout the infection, the density of uninfected RBCs, R, is affected by erythropoiesis, natural death, indiscriminate immunity, and merozoite invasion, such that:

(13)

In the absence of parasites, RBCs are replenished at a baseline rate, μ, that maintains a homeostatic equilibrium RBC concentration (R₀). During infection, the rate of erythropoiesis increases, and the extent of this upregulation is proportional to the deviation of RBC concentration from equilibrium [91], where ρ represents the maximum proportion of the RBC deficit recovered per day. RBCs are eliminated at a natural death rate of μ (which balances new production in the absence of infection), removed through indiscriminate immunity (–ln(1–N(t))), or converted to iRBC via invasion by merozoites.

We parameterized our model based on prior work on malaria intrahost dynamics [23,42] (Table 2). Because parameter values in Kamiya et al. [42] depend on initial dosage (I₀), we estimated our starting dose using the same method, which corresponds to approximately 10⁴-10⁵ iRBCs. Preliminary investigations indicate that varying the initial iRBC dose did not lead to quantitative differences in parasite relative fitness. Consequently, we parameterized and ρ using the median posterior estimates from [42] assuming an initial dosage of 10⁴ iRBCs.

Modeling transmission investment

In our model, transmission investment (c) is determined by a transformed cue variable, f_t, perceived α days ago by a parasite infecting a cell, to account for NCC (). Following the approach taken by Greischar et al. [23], we modeled the transmission investment reaction norm as the complementary log-log basis-spline with three degrees of freedom (B₃). Briefly, this method models the reaction norm as a linear combination of basis functions, generating a continuous, smooth cubic curve. The relative contribution of each basis function (also known as the model coefficient) defines the overall shape of the reaction norm (see Fig 1C). We further modify the basis-spline by applying a complementary log-log transformation to constrain the transmission investment value to be between zero and one. This method allows us to use numerical optimization to search for the set of model parameters (i.e., reaction norm shape) that maximizes parasite fitness. The shape of the reaction norm is constrained by several factors, including the inherent smoothness of cubic splines and the degrees of freedom. Cubic splines were previously shown to be sufficient to model optimal transmission strategies with respect to time, and further spline complexities led to negligible improvements in fitness [23]. For our models, increasing the flexibility of cubic splines also did not alter the shape of the optimal reaction norm for both log-transformed parasite and RBC density (Fig 3B). Because cubic splines are fundamentally smooth, the resultant reaction norms may have difficulty representing sharp transitions. For that reason, we used untransformed and log-transformed cue values to build our models and capture greater variation in putative reaction norms. Log transformation is used in lieu of other operations because it magnifies small changes when parasite density is low, and because parasites in vivo exhibit near-exponential growth in the initial phase of infection [10,96,97].

For numerical optimization, the perceived cue must be within the bounds of the reaction norm. To achieve this, the perceived cue variable governing transmission investment is transformed via a linear combination of Heaviside step functions (H) from the R package CRONE [98] that constrains the cue value to be between zero and the maximum value (max) defined in Table 3. In summary, the transmission investment function is defined as

(14)

where the Heaviside-transformed environmental cue is given by

(15)

For the best time-varying reaction norm, we defined the perceived cue as the time since infection, which has a range of 0-20 days and a step size of 0.001 (same as the simulation time step). In this case,

(16)

To incorporate dual cue perception, we used the mgcv package [99] to model transmission investment as the response variable of a generalized additive model. This model incorporates a full tensor product smooth of two environmental cues with nine degrees of freedom, represented as follows:

(17)

where transmission investment is governed by two Heaviside-transformed cues, f_t1 and f_t2, perceived α days ago. The transmission investment surface (i.e., 2-dimensional reaction norm) is generated by taking the tensor product of the basis-splines of two cues (B_3,1 and B_3,2). The shape of this spline surface is defined by the model coefficients .

Quantifying parasite fitness

To quantify parasite fitness, we calculated the cumulative transmission potential of different transmission investment strategies, following the method described by Greischar et al. [23]. For a given time step, transmission potential, , is defined as the probability that a mosquito vector becomes infected after biting an infected host and has been shown empirically to be a sigmoid function of gametocyte density in this parasite species [10],

(18)

The cumulative transmission potential, or overall fitness of the parasite, is calculated by integrating transmission potential over the entire infection duration of δ days:

(19)

This fitness metric assumes that the mosquito population is at a constant abundance and continuously biting, for the duration of the infection. Although this is a simplification of the evolutionary pressures faced by malaria parasites [36], it enables us to focus our investigation at the within-host level, where transmission investment “decisions" are being made and environmental cues are being perceived by parasites. In our numerical simulations, we estimated fitness with a step size of 0.001 days unless stated otherwise.

Cue choice and perception

Transmission investment has been shown to vary in response to within-host competition [11,12] and RBC age structure and availability [13,17–19,24], lysophosphatidylcholine level (a key component of the RBC plasma membrane) and exosome-like particles [100,101]. Whether parasites are responding to these cues directly, or if they are instead a proxy for changes in some other within-host environmental factor is unclear. To capture some of the diversity of putative cues, we defined transmission investment as a function of one of ten different host and parasite-centric measures (Table 3). Additionally, we implemented dual cue models, allowing parasites to simultaneously perceive two different cues. As transmission investment in P. falciparum varies non-linearly with cue strength (e.g., P. falciparum in response to lysophosphatidylcholine concentration [100]), and within-host dynamics are known to vary only when infectious dose changes by orders of magnitude [30], we introduced non-linear cue perception, enabling parasites to perceive cues either unaltered or log₁₀-transformed. We excluded merozoite density as a cue due to its short half-life and the rapid pace of RBC invasion [102]. We also did not include measurements of immunity as cues, given that our model does not track individual components of immunity (rather, it tracks immune killing, phenomenologically), nor do we have sufficient empirical evidence detailing the relationship between different immunity components and transmission investment.

For all cues, we adopt a standardized range that spans the values typically observed when infecting immunologically naive mice with P. chabaudi [11], given by the minimum and maximum values in Table 3. All values that exceed the maximum bounds are “perceived" by the parasite as the maximum bound via the modified Heaviside transformation function.

Fitness optimization

Acute P. chabaudi infection tends to resolve within 20 days [25]. For all models, we therefore confined our simulations and optimizations to a 20-day infection period using the parameters detailed in Table 2 and using the deSolve package in R [103]. To ensure the convergence and continuity of output dynamics, we simulated the models with a small time step size of 0.001 days using the vode solver [104]. For the dual cue model, we used a larger time step of 0.01 days (which does not greatly alter parasite dynamics but improves computational speed immensely).

To investigate how transmission investment strategies affect parasite fitness, we performed numerical optimization of the basis-spline (single cue) or generalized additive model coefficients (dual cue) to maximize the cumulative transmission potential. We optimized the models using two methods to increase the chances of obtaining the global fitness optima. First, we performed local optimization following the methodology outlined by Greischar et al. [23] with the same arbitrary initial values of 0.5 for all spline coefficients, but using the more efficient parallelized L-BFGS-B algorithm (Limited-Memory Broyden–Fletcher–Goldfarb–Shanno algorithm-B) implemented in the OptimParallel package [105]. Second, we optimized each model using a hybrid approach to avoid convergence to a local fitness maximum. This technique involved initial global optimization using the Differential Evolution (DE) algorithm from the DEoptim package [106] (maximum iterations = 500, step tolerance = 50, differential weighting factor = 0.8, cross-over constant = 0.5 (single cue) or 0.9 (dual cue), population size = 40 (single cue) or 90 (dual cue), time step = 0.01 (single cue) or 0.05 (dual cue)). DE is an efficient and robust algorithm with comparable/better performance relative to other stochastic optimization algorithms [107]. Following optimization by the Differential Evolution algorithm, we further optimized the objective function using the L-BFGS-B algorithm to pinpoint a probable fitness maximum. For the differential Evolution algorithm, we used a parameter bound of (–5,–500,–500,–500) and (5,500,5000,500) for the single cue model and (–10,–500,–1000,–1000,–250,–500,–1000,–500,–250) and (10,500,1000,1000,250,500,1000,500,250) for the dual cue model. These bounds effectively encompass the majority of basis-spline flexibility. From the two optimization results, we selected the parameter sets that resulted in the highest fitness.

For dual cue models, we implemented an additional optimization method that relies on the single cue models to inform the initial starting point. For dual cue models, we reasoned that the best single cue strategy out of the two provides a theoretical baseline fitness. In other words, the worst a dual-cue strategy should do is as good as the best single cue (out of the two cues), while being unresponsive to the second cue. As such, we applied the L-BFGS-B algorithm to derive a dual cue strategy that “resembles" the best single cue strategy, specifically by minimizing an objective function with two terms: (i) the summed absolute difference between the dual and single cue reaction norms, and (ii) the total standard deviation in transmission investment across the worse cue. Inclusion of the second term ensures that transmission investment is not affected by variations in the worse cue. Together, the objective function is represented by:

(20)

We first discretized the cue ranges into 500 evenly spaced values (, better performing cue; , worse performing cue). For each cue value with the index i, we calculate the deviation of the dual cue model using the absolute difference between , the dual cue transmission investment averaged over the worse cue, and , the corresponding transmission investment in the single cue model. The objective function also quantified the extent to which the worse cue affects transmission investment, specifically by summing the standard deviation across the worse cue value for each cue value indexed i.

Following this initial optimization, the resultant parameter is set as the starting point for an additional round of local optimization using the L-BFGS-B algorithm to maximize parasite fitness (same parameters as mentioned above).

As a further validation step, we compared the fitness of each putative single cue optimal strategy to the fitness obtained by 1000 randomly generated strategies. These random strategies were generated within the predefined bounds of (–5,–500,–500,–500) and (5,500,5000,500), established previously to encompass much of the spline space. In no case was the putative best strategy outperformed by a random one (S1 Fig).

Semi-stochastic infection model

In order to explore the impact of individual infection-level variation on putative optimal transmission investment strategies, we incorporated the following random variables in the within-host model: maximum proportion of RBC deficit recovered per day (ρ), parasite burst size (β), immunity activation strengths ( and ), and immunity half-lives ( and ). For each iteration, random variables were drawn from 167 Markov chains (selected at regular intervals from a pool of 16,000 chains) used to estimate the posterior distribution in [42]. For the maximum proportion of RBC recovered per day and parasite burst size, the posterior distributions (originally centered at zero) are re-centered at the deterministic model values of 0.263 and 5.721, respectively. Since the paper from which we draw parameters explored the impact of dose [42], we only considered estimates pertaining to hosts infected with an initial dose of 10⁴ parasites (constituting around 1000 parameter values), which corresponds to our choice for this model parameter (Table 2).

For each optimal transmission investment strategy (i.e., the best strategy for each cue or dual cue combination), we performed 20-day simulations with (i) all parameters allowed to vary from their deterministic default and (ii) only a single parameter varying. Here, we used a time step of 0.01 days, which does not change fitness substantially relative to performing the simulation with a time step of 0.001 days.

Disease mapping of laboratory mice infections

To explore the dynamics of host and parasite-centric cues in an in vivo setting, we curated data from publicly available sources of P. chabaudi infection dynamics in mice [11,43–46,48]. The curated dataset includes experiments that met the following criteria: a) single parasite strain infection, b) no administration of anti-malarial drugs, and c) availability of data on iRBC, RBC, and gametocyte densities. We then visualized two-dimensional disease maps (e.g., iRBC vs. RBC path plots) using the R package ggplot2 [108].

Quantifying parasite virulence

We used the area encompassed within a disease map (where the x-axis displays total iRBC density and the y-axis displays uninfected RBC density, as in [47]) as a proxy for parasite virulence. To calculate this area, we first used the chull function [109] to define the convex hull, the smallest set of coordinates that encloses the points within the infection trajectory. We then used the areapl function from the splancs package [110] to quantify the area of the convex hull.

S1 Fig. Validation of transmission investment optimization.

We conducted 1000 infection simulations for parasites sensing different cues, with each simulation using randomly generated transmission investment strategies. The fitness values of these simulations were compared to the fitness of parasites adopting the optimal transmission investment strategy. Each histogram represents the difference between the optimized fitness and the fitness obtained from the randomly generated strategy. The fitness of the random strategies is always lower.The data and code needed to reproduce S1 Fig can be found in DOI: 10.5281/zenodo.17114049, plos-bio_R-notebook.Rmd file (l296-318) and S6 Data.

https://doi.org/10.1371/journal.pbio.3003081.s001

(TIFF)

S2 Fig. Transmission investment dynamics produced by parasites sensing different dual cue combinations.

Transmission investment is plotted for all possible dual cue combinations, with heatmaps ranked in descending order (from top to bottom) based on their associated cumulative transmission potential. The data needed to reproduce S2 Fig can be found in S7 Data.

https://doi.org/10.1371/journal.pbio.3003081.s002

(TIFF)

S3 Fig. Terminal investment increases the fitness of dual cue-perceiving parasites.

The cumulative transmission potential of parasites sensing time, the best performing dual cue combination (R log & I log), or the individual single cues (I log, R log) are indicated with a black cross, blue square, orange circle, and yellow triangle, respectively. The cumulative transmission potential is normalized with respect to the fitness of the best time-varying strategy. Note that time, the dual cue combination, and I log all permit an initial delay in transmission investment, while R log does not; the fitness advantage of this delayed strategy can be clearly viewed on these plots as fitness accrues faster. Only time and the dual cue combination lead to terminal investment, leading to additional fitness gains from day 18 onward. The data needed to reproduce S3 Fig can be found in S8 Data.

https://doi.org/10.1371/journal.pbio.3003081.s003

(TIFF)

S4 Fig. Logged RBC and logged iRBC form a looping trajectory in silico.

Graphs display the two-dimensional reaction norm and infection trajectories for parasite sensing time with nine degrees of freedom. The trajectories depict the dynamics of the cues (indicated on the x- and y-axes) along with the corresponding transmission investment (represented by the colour of the line). Note that the dynamics of gametocytes and sexual iRBC are hard to display completely given that these densities start off at zero. The distance between each hollow circle represents one day. The data needed to reproduce S4 Fig can be found in S3 Data.

https://doi.org/10.1371/journal.pbio.3003081.s004

(TIFF)

S5 Fig. RBC/logged RBC and logged iRBC form a looping trajectory in vivo.

To identify combinations of within-host variables that exhibit a looping relationship, we plotted all possible permutations of experimentally derived P. chabaudi variables. Each line represents the dynamics of a single strain of infecting parasite. The line colour indicates the progression of time. The data needed to reproduce S4 Fig can be found in S3 Data.

https://doi.org/10.1371/journal.pbio.3003081.s005

(TIFF)

S6 Fig. Posterior distribution of parameter values used for the semi-stochastic model.

The histograms display the distribution of 1000 parameter values drawn from 167 Markov chains derived by Kamiya et al. [42]. The parameter values used in the deterministic model are indicated by the dotted line. The data needed to reproduce S6 Fig can be found in S9 Data.

https://doi.org/10.1371/journal.pbio.3003081.s006

(TIFF)

S7 Fig. Perceiving the optimal dual cue combinations increases parasite susceptibility to environmental and developmental stochasticity.

We conducted 1000 simulations, each lasting 20 days, in which parasites adopted the optimal transmission investment strategy (from the deterministic model) using a semi-stochastic model where only one parameter was randomly drawn from the posterior distributions (see S6 Fig). For each parameter that is varied (indicated in the facet panels), we plot the geometric mean of the fitness on the y-axis and the fitness from the deterministic model on the x-axis. All fitness values are normalized with respect to the best time-varying strategy. The dashed line represents a 1:1 relationship between the x-axis and the y-axis; the line is not visible on the plot for burst size, since variation in this parameter has the greatest negative influence on fitness. The data needed to reproduce S7 Fig can be found in S10 Data.

https://doi.org/10.1371/journal.pbio.3003081.s007

(TIFF)