How to quantify developmental synchrony in malaria parasites

Abstract

Malaria infections represent an iconic example of developmental synchrony, where periodic fevers can result when the population of parasites develops synchronously within host red blood cells. The level of synchrony appears to vary across individual hosts and across parasite species and strains, variation that—once quantified—can illuminate the ecological and evolutionary drivers of synchrony. Yet current approaches for quantifying synchrony in parasites are either biased by population dynamics or unsuitable when population growth rates vary through time, features ubiquitous to parasite populations in vitro and in vivo. Here we develop an approach to estimate synchrony that accounts for population dynamics, including changing population growth rates, and validate it with simulated time series data encompassing a range of synchrony levels in two different host-parasite systems: malaria infections of mice and human malaria parasites in vitro. This new method accurately quantifies developmental synchrony from per capita growth rates using obtainable abundance data even with realistic sampling noise, without the need to sort parasites into developmental stages. Our approach enables variability in developmental schedules to be disentangled from even extreme variation in population dynamics, providing a comparative metric of developmental synchrony.

1 Introduction

Malaria parasites replicate in the blood of their hosts, and when this replicative cycle is synchronized, it can result in the periodic fevers classically associated with malaria infection (Kitchen, 1949; Nerlich et al., 2008). Synchrony indicates the degree to which parasites within infected red blood cells (iRBCs) show similar developmental timing across the population of iRBCs that make up an infection, and it can have consequences far beyond triggering symptoms. Since antimalarial drugs are most effective against certain parasite developmental ages (ter Kuile et al., 1993), parasite populations may be vulnerable—or resilient—to clearance by fast-acting drugs depending on the timing and synchronization of development (White et al., 1992). Theory suggests that pathogenic organisms could evolve different developmental schedules as a form of non-classical drug resistance (Neagu et al., 2018). Such evolutionary changes in synchrony could have knock-on consequences for malaria epidemiology, since synchrony is thought to influence infection severity (Touré-Ndouo et al., 2009) and the odds of onward transmission (Greischar et al., 2014; Schneider et al., 2018). Synchrony also influences capacity to accurately measure population expansion, generating spuriously large parasite multiplication rates and making it challenging to quantify how novel therapeutics (e.g., drugs, vaccines) impact parasites in vivo (Greischar and Childs, 2023). Anticipating synchrony and its potential to evolve requires understanding the degree to which parasites show heritable variation in their developmental schedules, but quantifying that variation has so far been stymied by methodological challenges. Existing methods are incapable of quantifying synchrony in ways that enable comparisons across genetic backgrounds and environments (Greischar et al., 2019).

The difficulty of quantifying synchrony arises in diverse systems (including insect pests, Bjørnstad et al., 2016), but addressing the challenges for malaria parasites—where data are available from artificial culture, experimental animal models, and natural infections—offers scope for making the comparisons needed to resolve the proximate and evolutionary drivers of developmental rhythms. In the context of malaria parasites, synchrony refers to the period of development and multiplication within a red blood cell (RBC)—the “intraerythrocytic development cycle” or IDC—that is completed when the RBC bursts to release merozoites, stages capable of invading RBCs. While periodic fevers are indicative of synchrony, the iRBC abundance required to trigger fevering (the “pyrogenic threshold”) varies across hosts (reviewed in McKenzie et al., 2008), and infections can be synchronized even when hosts show no outward symptoms (Färnert et al., 1997). Fevers are therefore an inconsistent marker of synchrony. In practice, synchrony is nearly always quantified as the percentage of parasites in a particular age range (i.e., a morphologically distinct stage of the IDC; Lambros and Vanderberg, 1979; Deharo et al., 1994, Deharo et al., 1996; Reilly et al., 2007; Touré-Ndouo et al., 2009; Allen and Kirk, 2010; O’Donnell et al., 2011), a metric that suffers from multiple sources of bias. Attempting to distinguish developmental age from morphology is inherently subjective (Ciuffreda et al., 2020) and the duration of a particular morphological stage can vary across genotypes (e.g., the early “ring” stage in Plasmodium falciparum, Reilly et al., 2007).

Distinct from these practical challenges, stage percentages represent a biased estimate of synchrony because parasite age distributions vary with population growth rates (White et al., 1992; Khoury et al., 2014; Greischar et al., 2019). This bias will be far more widespread—impacting infections in vivo and in vitro—since malaria infections exhibit extreme variability in population dynamics, such that parasite numbers within an infection can vary over several orders of magnitude on the timescale of days or weeks (e.g., Miller et al., 1994; Huijben et al., 2010; Wacker et al., 2012). Robust comparisons can be made by focusing on time windows where population dynamics are similar across treatments (e.g., comparing synchrony in typical versus perturbed host feeding rhythms, Prior et al., 2018), but that approach cannot be generalized, e.g., to the case where treatments alter parasite population dynamics. In vitro experiments reflect this quandary: Allen and Kirk (2010) hypothesized that shaking P. falciparum cultures would help maintain synchrony by reducing spatial variation in nutrient availability and concomitant variability in developmental duration, and they tested that hypothesis by checking for a greater percentage of iRBCs in the early part of development in shaken versus static cultures following initial synchronization. However, shaking cultures also accelerates population expansion by reducing merozoites that are wasted invading RBCs that have already been infected, as would be frequent in static cultures (Allen and Kirk, 2010). Thus, a greater percentage of young iRBCs would be expected even if shaking did not reduce variability in developmental duration. Recent efforts to quantify synchrony have focused on using gene expression data to obtain parasite age distributions (Ciuffreda et al., 2020), but the fundamental problem remains: age distributions, no matter how well-resolved, are still influenced by population dynamics in ways that preclude broader comparison of developmental synchrony (Greischar et al., 2019).

The population dynamics that undermine the stage percentage approach also preclude the most commonly used method for quantifying synchrony that does not rely on age distribution data. In some malaria species, iRBCs in the latter half of parasite development sequester in the microvasculature of their hosts where they cannot be easily sampled (e.g., P. falciparum; White et al., 1992), causing periodic fluctuations in iRBC abundance where the period corresponds to the duration of the IDC. The fluctuations can be fit with a model to quantify synchrony. Those models assume log-linear expansion (i.e., a constant rate of parasite multiplication) and fit a sine wave, the amplitude of which serves as a metric for the degree of synchrony (e.g., Simpson et al., 2002; reviewed in Gnangnon et al., 2021). That approach sidesteps the bias inherent to stage percentages but cannot be applied broadly, for example, when sequestration is absent or when parasite multiplication rates change through time. Sequestration is not ubiquitous and cannot take place in vitro or in certain hosts (e.g., rag1^−/− mice, Khoury et al., 2014). While parasites often multiply exponentially initially, that log-linear expansion cannot continue indefinitely as resources become exhausted and immune pressure mounts, so the latter portions of time series must be discarded prior to fitting the model. It is also not obvious how the amplitude of fitted sine waves can be translated into a biologically meaningful and readily comparable measure of developmental synchrony, especially for malaria species that do not sequester and hence do not exhibit oscillations in observable iRBC abundance.

Here we present a new method for quantifying synchrony that treats the age distribution as unobserved data and instead quantifies developmental synchrony by fitting a simple model, a Leslie matrix (Leslie, 1945), to infection time series (specifically, repeated measures of infected RBC densities). This approach estimates the time window over which a cohort of parasites within an infection complete their IDC and burst out of RBCs: a short window, relative to the cycle length, indicates a highly synchronized infection, while a long window indicates asynchrony. Since no existing approach can provide accurate estimates of synchrony from empirical data against which to compare (Greischar et al., 2019), we instead validate this approach by applying it to simulated data (i.e., where the true answer is known). We do this for two malaria species, Plasmodium falciparum and P. chabaudi, that show considerable differences in the duration of the IDC and population dynamics. Our approach can recover differences in synchrony using tractable sampling schedules, and, equally important, can detect similar levels of synchrony despite differences in IDC duration and population dynamics across species. This method relies on iRBC counts sampled multiple times per IDC rather than labor-intensive morphological staging. The focus on stage percentages has left a dearth of relevant data sets for applying this method, but these data are readily obtainable whenever sequestration does not obscure iRBC abundance (e.g., in vitro). Unlike existing metrics, this definition of synchrony enables comparisons across species and environments even when population dynamics also vary, providing a framework for understanding the evolutionary drivers and practical consequences—including the capacity for harm and vulnerability to control—of synchrony.

2 Methods

We develop an approach for measuring synchrony by fitting a model to per capita replication rates calculated from iRBC time series. We estimate synchrony from the age distribution of parasites within iRBCs, similar to past efforts (Ciuffreda et al., 2020) but with an important difference: instead of directly measuring age distributions—a composite measure of synchrony and population dynamics—we fit the initial age distribution with a model that accounts for subsequent changes in age structure caused by population dynamics. We focus on initial ages because changes in the age distribution through time reflect both synchrony, the variation in developmental timing across the population of iRBCs, and population growth (Greischar et al., 2019). The simple model we fit to data to estimate synchrony represents a null hypothesis for how age distributions would change through time in response to population dynamics if synchrony were constant. Other hypotheses can then be tested against the null model, e.g., that synchrony increases or decreases with time, or that there are multiple cohorts with distinct levels of synchrony. Ideally, we would apply this method to real data, comparing its performance against the current gold standard approach. However, in this case the gold standard metric—the percentage of parasites in a particular developmental window—is known to be unacceptably biased (Greischar et al., 2019), so we instead simulate data where the true level of synchrony is known (that is, specified in the mechanistic model) and constant. We use those simulated data both to illustrate the bias inherent to the gold standard approach and to validate our novel approach.

2.1 Simulated time series data

We use a previously described mechanistic model (Greischar et al., 2014) to simulate experimental infections of mice with the rodent malaria parasite, P. chabaudi, and in vitro cultures of human malaria parasites, P. falciparum. These two scenarios exhibit divergent population dynamics: P. chabaudi requires 24 hours for the IDC and RBCs are replenished by the host causing iRBC abundance to rebound from infection-induced decline; P. falciparum requires 48 hours for the IDC and replicates in an artificial culture of RBCs without RBC replenishment (in our simulations). The model tracks the abundance of infected and uninfected RBCs, as well as the number of RBC-invasive merozoites, per microliter (details in Supplementary Methods). We simulate four different levels of synchrony, identical for each of the two scenarios (Figure 1), and sample our simulated infections on empirically feasible schedules. To test low resolution time series, we simulate “even” sampling where infections are censused twice per IDC with equal time between each sample (12 hours apart for P. chabaudi and 24 hours apart for P. falciparum). We also examine “uneven” sampling where more samples are taken per IDC and sampling is concentrated around peak bursting. For uneven sampling, we assume only three samples per IDC are possible for P. chabaudi infections in vivo (gathered in an eight hour window around peak bursting), whereas we simulate sampling P. falciparum in vitro every two hours for the 12-hour period spanning peak bursting (seven samples per IDC, as described by Reilly Ayala et al., 2010). By using distinct sampling schedules for the two species, we also test whether the new approach can cope with practical differences in the timing and frequency with which malaria species can be assessed.

Figure 1

We use simulated iRBC abundance as is and also with realistic sampling noise. We assume that iRBC abundance will be assessed via microscopy, since qPCR estimates parasite genome copies that can only be related to iRBC abundance if the level of synchrony is known (e.g., Cheesman et al., 2003). Microscopy protocols for estimating iRBC abundance via microscopy vary, so we test the consequences of either counting iRBCs until a target total number of RBCs have been observed (binomial error) or counting to a target iRBC number (negative binomial error). Details of the model, parameter values, and simulated sampling schedules and sampling noise can be found in Supplementary Methods. We confirm that simulated data generate biased estimates of synchrony using existing methods (percentage of iRBCs in the early half of development) and then test our new approach for quantifying synchrony.

2.2 A new metric for developmental synchrony

Our metric for synchrony requires estimating the initial iRBC age distribution, while accounting for changes due to population dynamics. We do this by fitting a symmetric Beta distribution as detailed below. A Beta distribution is flexible enough to yield iRBC age distributions ranging from uniform to narrow bell curves (Figure 1). The fitted initial iRBC age distribution can then be used to demarcate the time required for 99% of iRBCs to burst (as in Figure 1A), and dividing that quantity by the duration of the IDC can serve as a measure of developmental synchrony that can be readily compared across species:

We subtract the fraction from one to arrive at a measure of synchrony that scales from zero (completely asynchronous) to one for perfectly synchronized infections, consistent with past efforts to establish a “synchronicity index” (Deharo et al., 1996). The time required for 99% bursting corresponds to the 0.005 and 0.995 quantiles of the symmetric Beta distribution given by the shape parameter s_P, multiplied by the duration of the IDC. Thus, for the same shape parameter, the time required for 99% of bursting would take twice as long for P. falciparum compared with P. chabaudi (Figure 1), since P. falciparum requires twice as long to complete its cycle (reviewed in Paul et al., 2003; Mideo et al., 2013).

By default, the initial age distribution (the symmetric Beta distribution) is centered halfway through development, but that need not be the case in real populations. We therefore fit an additional offset parameter that reflects the initial median age of the cohort, around which the Beta distribution is centered. In this way, we can independently specify the median age and the variability around that median age. In contrast, if the initial age distribution were specified with an asymmetric Beta distribution, the shape parameters would alter the mean and the variance simultaneously. The choice of 99% of the time required for bursting can yield values comparable to the full duration of the IDC for asynchronous infections, and hence a synchrony of zero according to Equation 1. Note that synchrony can be calculated given an age distribution (here specified by one parameter, s_P) and the duration of the IDC. We test our approach on simulated data by comparing known s_P values (i.e., the values used to simulate the time series) with estimated values () obtained by fitting a simple model (next section). We examine the simpler case of synchrony being maintained through time and later discuss how this approach could encompass changing synchrony.

2.3 Projecting population dynamics using Leslie matrices

We next derive a simple model to fit to time series data and thereby estimate the level of synchrony. In formulating that model, we focus on one of the simplest ways of modeling the continual feedbacks between population dynamics and age distributions, a Leslie matrix (Leslie, 1945; reviewed in Hastings, 1997). Given information on the number of new iRBCs generated by each bursting schizont, a Leslie matrix can be used to project the population expansion (or decline) of an age-structured population, illustrated with a hypothetical example of a population with four age classes in Figure 2. In the context of malaria infections, a Leslie matrix can be used to specify how iRBCs containing parasites in different developmental age classes transition to subsequent developmental age classes or burst to generate new iRBCs (Figure 2B).

Figure 2

We assume that the iRBC population is initially distributed across n developmental age classes, that is, a discretization of the continuous age distribution used in the data-generating model. The fraction of the IDC spent within each class is given by 1/n (e.g., if samples were taken, at minimum, 4 hours apart for a species with an IDC lasting 24 hours, then n = 24/4 = 6, and each age class lasts 4 hours). When iRBCs are sampled at evenly spaced intervals throughout the developmental cycle, n is equal to k, the number of samples per cycle (e.g., in Figure 2, n = k = 4). If sampling is uneven, n is equal to the duration of the IDC divided by the minimum time between samples. For example, if sampling occurred three times per cycle (k = 3) so that iRBCs were censused four hours prior to peak bursting, at peak bursting, and four hours after peak bursting, we treat that situation as though iRBCs were counted every four hours (with missing counts for three of the six samples), so that again n = 24/4 = 6.

For a population with n age classes, the Leslie matrix (L) is

where m values indicate the fecundity of each age class and θ values give the probability of individuals in a particular age class surviving to mature into the subsequent age class in the next time point. We assumed that each age class is certain to survive and mature into the subsequent age class at the next time step (θ₁ = θ₂ = θ₃ =… = θ_n−1 = 1). In the simulated data, approximately 2% of iRBCs will not survive development due to background mortality of RBCs (µ, see Supplementary Tables S1, S2). If we are nonetheless able to recover the true duration of bursting, it would suggest the approach is robust to at least small deviations from the assumption of 100% survival of iRBCs during the IDC. Further, background RBC mortality is likely to be the dominant force removing iRBCs in the acute portion of infection, where RBC limitation can largely explain dynamics (Mideo et al., 2008) and immune clearance is thought to be minimal (White et al., 1992). Immunity will certainly be minimal in vitro. Note that immune clearance of merozoites would be accounted for through the calculation of parasite multiplication rates (PMRs), the fold-increase in iRBC abundance over each IDC represented in the time series (defined in Table 1). However, when immune clearance of developing iRBCs is substantial, the assumption of 100% iRBC survival through each IDC may need to be revisited.

Since only the final age class of malaria parasites are capable of bursting out to generate new iRBCs, m₁ = m₂ = m₃ = … = m_n−1 = 0 and m_n= B_t, where B_t is the number of new rings generated by each iRBC completing development (Figure 2A). Note that B_t is not equivalent to the burst size, the average number of merozoites emerging from each bursting schizont, because some emerging merozoites will fail to successfully invade. Once B_t values are obtained by interpolating between PMR values (see below), the Leslie matrix (Equation 2) can be parameterized for each time point at which population projections are needed. Given a vector to indicate the numbers of iRBCs in each age class at time t,

The number of iRBCs in each age class is then summed to obtain the projected total number of iRBCs at each time point. Projecting the total iRBC abundance (Equation 3) enables calculation of the predicted per capita growth rate, for each time point, where

and the _i subscript indicates the developmental age class.

The initial number of iRBCs in each developmental age class, denoted by the vector , can be specified with a symmetric beta distribution, which has the useful property of ranging from uniform (i.e., completely asynchronous, when the shape parameter, s_P, is equal to one) to increasingly narrow bell-shaped curves with increasing s_P values (Figure 1A). For a symmetric beta distribution, the mean falls at 0.5 by default—equivalent to the median parasite age falling exactly halfway through the IDC—but that need not be the case for real data. We therefore fit an “offset” parameter, constrained to vary between zero and one, to allow the mean of the beta distribution (that is, the median parasite age) to deviate from 0.5. We use the cumulative density function of the beta distribution to obtain the initial stage distribution, , for a given number of stage compartments (n). The initial I₀ iRBCs are sorted into n stages according to the beta distribution specified by a given shape parameter, (the quantity to be fitted). The initial number of iRBCs is arbitrary, since we are fitting per capita replication rates (rather than iRBC abundance), so we set I₀ = 100 to avoid errors of precision that might occur with lower iRBC counts.

2.4 Estimating developmental synchrony

We describe the algorithm for estimating synchrony using Equation 1 (illustrated in Figure 3):

• Step 1. Obtain time series of iRBC abundance, sampled multiple times per IDC. Age distribution data are not required.

• Step 2. (a) We first calculate the parasite multiplication rates (PMRs) and then (b) calculate the per capita replication rates (R_e, Table 1).

• Step 3. By interpolating between the observed PMRs, we obtain estimated B_t values needed to project population dynamics for an age-structured population (see Figure 2).

Figure 3

The subsequent steps 4 and 5 are carried out numerous times, using an optimization algorithm (details follow) to determine the starting age distribution that yields predicted values that best match the observed R_e values:

• Step 4. We project iRBC abundance for a large number of starting age distributions using a Leslie matrix (see example in Figure 2). The starting age distribution is specified as a symmetric Beta distribution with a shape parameter that will be fitted () and used to calculate the time required for 99% bursting and the level of synchrony according to Equation 1. Though not used to calculate synchrony, we also fit an offset parameter to allow the median parasite age to deviate from the default of half the IDC duration.

• Step 5. Each projected iRBC abundance is subsampled to obtain values at the same times reported in the actual data. That is, any time points that are missing in the observations are dropped prior to calculating predicted per capita replication rates () for comparison with observed per capita replication rates (R_e).

• Step 6. Identify the initial age distribution associated with the best fitting (Equation 4). Determine the time required for 99% bursting in that initial age distribution by subtracting the 0.5% quantile from the 99.5% quantile of the Beta distribution defined by the best fit . Calculate synchrony according to Equation 1.

2.5 Numerical optimization and estimates of synchrony

We use the Nelder-Mead algorithm to locate the value of (and the offset parameter) that minimizes the sum of the absolute error between the observed and predicted per capita replication rates (R_e and , respectively). The machinery to perform the optimization is readily available from R (R Core Team, 2020; see Supplementary Methods for details and annotated code). For each simulated time series, we run the optimization 1000 times with randomly chosen starting parameters to increase the chances of arriving at a globally optimal best fit. We then retain the best fit shape parameter () and offset parameter from those 1000 optimizations (i.e., the values corresponding to the fit yielding the lowest sum absolute error) and use it to calculate the time for 99% bursting and then synchrony according to Equation 1.

3 Results

3.1 Stage percentages represent biased synchrony estimates

Our simulated data reveal problems associated with the common approach of quantifying synchrony as the percentage of iRBCs in a particular life stage (reviewed in Greischar et al., 2019). Since the morphologically distinguishable “ring” stage occupies roughly the first half of intraerythrocytic development (Cheesman et al., 2003; Reilly et al., 2007), 50% rings would correspond to asynchrony in a population that is neither expanding nor contracting, while 100% rings is typically thought to indicate perfect synchrony (Greischar et al., 2019). The full simulated dynamics of percentage rings (Supplementary Figure S1) already indicate challenges with this metric: only if infections are sampled at the perfect point in time, and with perfect accuracy, can we recover the correct rank order of synchrony. Even then it is difficult to correctly identify the asynchnronous infection because the percentage rings consistently deviates from the null expectation of 50%. Of course, with real data, this resolution is not typically possible, nor is the perfect accuracy of our simulated data. Instead, we plot the data only from sample time points when maximal ring percentages would be expected in a synchronous infection (i.e., the first half of the IDC, for each of the IDCs simulated) and add realistic measurement error (Figure 4). Depending on the sampling schedule, highly synchronized infections can be indistinguishable from the asynchronous expectation or may appear “semi-synchronized” (a bare majority of parasites in the ring stage, e.g., ∼64%; Khoury et al., 2017). Note that simulations of P. falciparum in vitro yield higher percent rings on average (Figure 4B, D), since that population is only expanding, in contrast to the simulations of P. chabaudi, which expand and then decline (Figure 1B, C). The age distribution assessed at a single time point—the equivalent of choosing one of the point estimates on these plots—will frequently fail to detect even substantial rank order differences in synchrony, but synchrony is frequently assessed in this way (e.g., Touré-Ndouo et al., 2009; Ciuffreda et al., 2020).

Figure 4

Khoury et al., 2014). Horizontal lines indicate expected values for perfect synchrony (100%, long dash), asynchrony (50% rings in a static population, dotted line), or a semi-synchronous population (64% rings, following Khoury et al., 2017). The corresponding initial age distribution and iRBC abundance are shown in Figure 1 (details on data generating model in Supplementary Methods, parameter values in Supplementary Tables S1, S2).

3.2 The new approach recovers true differences in synchrony

3.2.1 The effect of sampling timing

Differences in synchrony can be difficult to discern from iRBC trajectories due to large changes in abundance that can occur in a matter of days, even with extremely high resolution time series in the absence of sampling noise (Figure 1B, C). Per capita replication rates make distinct levels of synchrony much clearer (R_e, Supplementary Figure S2A, C). However, per capita replication rates vary with the sampling schedule, and sparse sampling reduces the ability to visually detect differences in synchrony (e.g., twice per IDC; Supplementary Figure S2B, D). We first examined the fitted per capita replication rates compared to simulated P. chabaudi and P. falciparum values assuming no sampling noise. The simple fitted model (colors as in Figure 1) yields predicted values that are very close to the R_e values from simulated data (Supplementary Figures S3, S4).

Model fitted estimates represent the best fit from 1000 randomly chosen starting values, so we re-ran fits 20 times on simulated time series, without sampling noise, to estimate variability in the fitting process itself. When sampling was uneven, the fitted values did not vary, but even sampling twice per IDC yielded a range of estimates for the duration of bursting and hence synchrony (closed squares in Figure 5). For even sampling only, multiple combinations of median parasite age (fitted with the offset parameter) and duration of bursting yielded equally good fits (Supplementary Figure S5). This issue—nonidentifiability—was foreshadowed by the extremely similar per capita replication rates (R_evalues) for all but the asynchronous simulations (Supplementary Figure S4). In contrast, R_e trajectories are distinct for each level of synchronization with uneven sampling (Supplementary Figure S3) and identifiability is not an issue. We chose the timing of uneven sampling so that more samples were taken around the median time of bursting for each cohort, since that will yield the largest changes in iRBC abundance and hence the largest R_e values. If the time of median bursting is known, clustering samples around that time is likely to make it easier to detect differences in synchrony. However, if the timing of bursting is not known, uneven sampling would still be useful to avoid nonidentifiability between bursting duration and median parasite age.

Figure 5

Focusing on the fits to noiseless data, unevenly sampled, we obtain estimates of bursting duration and synchrony extremely close to their true initial values (open circles, Figure 5). The fits to simulated P. chabaudi infections appear to modestly overestimate the duration of bursting, and hence underestimate the true level of synchrony (Figure 5A, C), but that does not represent a problem with the method itself. Rather, these estimates are detecting a modest decline in initial synchrony over the 16 IDCs simulated, due to very slight differences in IDCs within each cohort arising from the exponentially-distributed waiting times in the merozoite stage assumed in the underlying model. High levels of initial synchrony result in discrete pulses of merozoite abundance that become wider through time as synchrony declines (Figure S6). The P. falciparum fits are to shorter time series—4 IDCs—and the synchrony estimates are therefore nearly identical to their initial true values (Figures 5B, D). When we recalculate synchrony using the mean hours for 99% bursting calculated from simulated merozoite abundance (details in Supplementary Code), we find that the new method returns estimates nearly identical to the true value (Figure S7). Despite large differences in iRBC dynamics (Figure 1B, C), the estimates of the duration of bursting (relative to IDC length) and synchrony are extremely similar for P. chabaudi and P. falciparum (Figures 5, S7). Thus, this approach yields synchrony estimates that account for—and can therefore be uncoupled from—the population dynamics that render existing methods ineffective.

3.2.2 The effect of sampling noise

We then tested the new approach on the same underlying time series with simulated sampling noise assuming percent parasitemia was subject to a binomial or negative binomial distribution. For simulated P. chabaudi infections, fitted estimates of synchrony matched well with their true initial values on average, but binomial error in percent parasitemia yielded a wide range of fitted values that could sometimes suggest the incorrect rank order of synchrony (Figure 6A). For negative binomial error, estimated synchrony is clustered tightly around the true value, though the highest levels of synchronization may be difficult to distinguish from one another (Figure 6B). Initial fits to noisy P. falciparum time series using the same target counts as for P. chabaudi (500 for binomial noise and 100 for negative binomial) gave an unacceptably wide range of estimates and poor match to true values (Supplementary Figure S8). Fits to noisy P. falciparum time series show a wider spread of estimated values, since parameter values were chosen to keep percent parasitemia below approx. 10% (Supplementary Table S2). Using larger, but still plausible, target counts (2000 and 400, respectively), we find that average fitted values can deviate from the true rank order when error is binomially distributed (Figure 6C), but the true rank order can usually be recovered when error is negative binomially distributed (Figure 6D). Noisy asynchronous time series present the biggest challenge to accurate estimation of synchrony, since the ratio of signal to noise is lowest. That is, oscillations in per capita replication rates (R_e values) are minimal (panels D and H in Supplementary Figures S3, S4) so that sampling noise can easily masquerade as a higher level of synchrony. Taken together, these fits suggest that sampling to a target iRBC count should make it easier to distinguish the true level of synchrony, and that greater effort (e.g., higher target counts) will be needed for P. falciparum infections or other species where percent parasitemia tends towards lower values.

Figure 6

Reilly et al., 2007), so larger target RBC (iRBC) counts were used for (C, D), following reported experimental protocols (see details in Supplementary Methods). For corresponding plots showing the estimated duration of bursting, see Supplementary Figure S9.

4 Discussion

Developmental synchrony has the potential to alter malaria parasite fitness by influencing replication rates, drug efficacy, and transmission success (Hawking et al., 1968; Hawking, 1970; White et al., 1992; Greischar et al., 2014; Schneider et al., 2018; Owolabi et al., 2021; O’Donnell et al., 2022; reviewed in Mideo et al., 2013; Greischar et al., 2019; Prior et al., 2020). Further, parasites at different developmental ages vary in their susceptibility to host defenses (e.g., γδ T cells; Costa et al., 2011) and antimalarial drugs (Yayon et al., 1983), including the current front-line drug, artemisinin (ter Kuile et al., 1993; Owolabi et al., 2021). The ability to account for underlying differences in population dynamics is central to comparing synchrony across environments, populations and species. Intuitive metrics (e.g., the percentage of iRBCs in a particular developmental stage) fail to account for population dynamics and yield biased estimates of synchrony (Greischar et al., 2019), except where population dynamics do not vary across treatments (Prior et al., 2018). So far the only method that sidesteps the bias generated by population dynamics is fitting a dynamical model to data (White et al., 1992; Simpson et al., 2002; reviewed in Greischar et al., 2019), but these models assume constant population growth rates, preventing comparisons within and across infections. We developed a new approach, validated against simulated data, that can distinguish the level of synchrony while accounting for considerable variation in population growth rates.

The Leslie matrix model presented here describes the dynamics for a population that begins at a particular level of synchrony and whose age distributions are subsequently unperturbed except by replication. That model is nonetheless able to return the correct average level of synchrony when the underlying dynamics deviate from the assumption of fixed IDC duration and synchrony changes through time (Supplementary Figure S7). As is, these methods could be applied whenever the goal is to detect average differences in synchrony, e.g., the problem of determining which treatments maintain higher levels of synchrony (such as shaking in vitro P. falciparum cultures to maintain uniform distribution of nutrients, Allen and Kirk, 2010). When synchrony is hypothesized to be changing through time, whether due to processes internal or external to the organisms in question, the fitted Leslie matrix model represents the null expectation, serving as a starting point to compare against more complicated scenarios. For example, this approach could also be extended to quantifying how synchrony changes over the course of infection by fitting the model to portions of the time series. Each time window would then have a different fitted shape parameter , and by extension, a different estimated duration of bursting and level of synchrony according to Equation 1. If splitting the time series into multiple windows results in a better fit compared to fitting a single shape parameter for the entire infection—reducing the sum absolute error enough to justify the added complexity of the fitted model—that would suggest that the level of synchrony is changing through time. A similar approach could be used to test for deviations from our null assumption of 100% survival across all developmental classes. For example, drugs may alter synchrony (i.e., the age distribution) directly since anti-malarial drugs disproportionately impact parasites in the middle of development when they are most metabolically active (ter Kuile et al., 1993). As a consequence, population growth rates will also be reduced, with indirect impacts on age distributions. A model that incorporates enhanced removal of middle-aged parasites would allow for direct impact of drugs on parasite age distributions, while in our null model, drugs could only impact age distributions via changes to population growth rates. Comparing these models would reveal how drugs influence synchrony.

Our results have implications for the type of data and the frequency and duration of sampling needed to quantify synchrony in other systems, pathogenic and free-living. Crucially, age distribution data are not required, only estimates of total population size sampled multiple times per generation. Sampling at even intervals through a generation is neither necessary nor desirable, since that makes it difficult to estimate synchrony separately from median age. We show that it is possible to quantify synchrony from as little as four generations of data, i.e., four IDCs for malaria parasites. Estimating synchrony from even shorter time series may be feasible depending on the error distribution of the count data. Ideal methods for estimating population size would produce consistent coefficients of variation whether individuals are abundant or rare (e.g., negative binomial error distribution). For the case of malaria infections, our new approach allows for sampling schedules that are entirely feasible in experimental animal models (Deharo et al., 1994, Deharo et al., 1996; O’Donnell et al., 2011; Prior et al., 2018; O’Donnell et al., 2022) and even controlled human infection trials (reviewed in Duncan and Draper, 2012). For plants and animals, a wealth of data exist that can be used to parameterize matrix models (i.e., age-specific survival and fecundity; Salguero-Gomez et al., 2015, Salguero-Gomez et al., 2016), though the goal is typically to estimate population expansion rates from age distribution data (Hernández et al., 2023) rather than the reverse.

The methods in the present study rely on time series of abundance, and there are important challenges to obtaining accurate abundance data. In particular, the question of whether synchrony can be resolved from the data depends on the level of sampling noise, which for malaria parasites requires replicate samples from the same infection at the same point in time. Quantifying sampling noise would not require time series, though it would be helpful to obtain replicate samples from infections with different iRBC abundances, so that the relationship between the mean and variance of iRBC counts can be established. These data exist for PCR-derived counts in P. chabaudi (Mideo et al., 2008; Huijben et al., 2010; Miller et al., 2010), but comparable data from microscopy-derived iRBC counts appear to be lacking in both P. chabaudi and P. falciparum. Estimates of that sampling error are needed to make specific recommendations about experimental design, including sampling frequency and duration. All else equal, per capita replication rates will approach one (the population replacement level) as sampling resolution increases, because increases in iRBC abundance will be subdivided into smaller and smaller intervals. That could make it more difficult to distinguish synchrony-driven differences in per capita replication rates and may represent another reason to aim for clustering sampling around median bursting times, when changes in iRBC abundance—and therefore in R_e values—will be the greatest. Further, if error distributions were known, it would be possible to implement a maximum likelihood approach rather than minimizing the sum absolute error as we have done here.

The other major challenge is that some developmental ages are more difficult to sample. For example, the mature stages of some Plasmodium spp. sequester where they cannot be sampled, at least in certain hosts. When sequestration occurs, it would lead to underestimates of total parasite abundance (e.g., P. falciparum, White et al., 1992; P. vivax, Carvalho et al., 2010; P. berghei in wild type but not rag1^−/− mice, Khoury et al., 2014). If the level of synchrony is changing through time, the bias introduced by sequestration may likewise fluctuate, making it difficult to compare synchrony estimates from different phases of infection or across different infections. Sequestration is not universal, so some host-parasite combinations would not suffer from this bias (including P. berghei in rag1^−/− mice (Khoury et al., 2014), any P. falciparum infections in vitro). Sequestration is also a problem for studies attempting to quantify synchrony in vivo from age distribution data (Deharo et al., 1994, Deharo et al., 1996; Touré-Ndouo et al., 2009; O’Donnell et al., 2011; Ciuffreda et al., 2020), along with the added challenge of the bias introduced by population dynamics (reviewed in Greischar et al., 2019). Accounting for the bias introduced by sequestration will be crucial for making comparisons across parasite species that also vary in their propensity to sequester as well as across environments where the potential for sequestration varies (e.g., in vivo and in vitro conditions).

Empirical advances offer some promise for sidestepping the problem of sequestration to estimate total iRBC abundance from proxy measures. Methods exist to quantify the total number of parasite genomes from volatiles in breath samples (Berna et al., 2015), plasma biomarkers (Dondorp et al., 2005), and—for experimental animal models—from transgenic parasite strains expressing luciferase or other markers (Franke-Fayard et al., 2005). Caution is warranted in using these methods, since they quantify parasite genomes, rather than iRBC abundance. Unfortunately, those two quantities are only equivalent for parasite populations that are synchronized and have not yet begun DNA replication (e.g., Cheesman et al., 2003), so studies attempting to quantify total parasite biomass must make assumptions regarding the synchrony of the parasite population being examined. The present study suggests that if biomarker abundance were quantified at multiple time points per IDC, it may be possible to fit rather than assume the underlying level of synchrony by modifying the approach introduced here. As with microscopy-estimated counts, quantifying the sampling error distribution associated with these newer methods will be crucial to determining feasibility.

Perhaps even more challenging is the fact that time series of iRBC abundance are not typically gathered when the goal is estimating synchrony; rather, stage percentages are overwhelmingly used to compare synchrony through time and across different infections (reviewed in Greischar et al., 2019). The focus on stage percentages is understandable, given that researchers often require parasite populations enriched for a particular developmental age class (for example, when testing stage-specificity of antimalarial drugs, e.g., ter Kuile et al., 1993). However, knowing the stage percentages at a particular point in time does not enable the projection of future stage percentages, especially when parasite numbers are changing (Greischar et al., 2019). Since it utilizes replication rates, our new approach can project how age distributions will change through time for a given initial level of synchrony. Projected age distributions could prove useful, for example, in predicting when parasite populations will be most vulnerable to drug treatment or immune clearance.

Existing methods for quantifying synchrony are useful when abundance is steady through time and the duration of morphologically-distinguishable developmental stages is known in advance (reviewed in Greischar et al., 2019). Those conditions are likely to hold for only a small minority of cases, but the method validated here can estimate synchrony when numbers are changing—and in fact works best when populations are undergoing dramatic changes in abundance—and requires no prior knowledge of stage durations. The key requirements of this approach are time series of iRBC abundance, sampled multiple times per developmental cycle, and knowledge about the duration of development. While the focus with stage percentage data has left a dearth of requisite data, they are obtainable (O’Donnell et al., 2022). Further, given the kind and resolution of data required, the approach is likely to have utility in other organisms, both pathogenic and free-living. Past studies quantifying changes in synchrony in free living organisms have relied on time series of abundance (e.g., insect pests, Yamanaka et al., 2012; Nelson et al., 2013; Bjørnstad et al., 2016), but we have shown that a useful signal of synchrony also emerges from per capita replication rates. Whether models fit abundance or per capita growth rates, accounting for population dynamics is necessary to identify ecological processes that maintain versus erode synchrony and to determine whether synchrony is heritable and capable of evolving. These unresolved questions carry enormous practical significance—the answers can inform the schedule of drug treatment (or for free-living pests, chemical controls) and enable predictions about evolutionary responses to intervention efforts. For malaria parasites, better understanding of variation in synchrony—which requires robust methods to estimate it—could improve understanding of pathogenesis and epidemiology and inform efforts to intervene at the individual and population level.

References

• 1

  AllenR. J. W.KirkK. (2010). Plasmodium falciparum culture: the benefits of shaking. Mol. Biochem. Parasitol.169, 63–65. doi: 10.1016/j.molbiopara.2009.09.005

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  BernaA. Z.McCarthyJ. S.WangR. X.SalibaK. J.BravoF. G.CassellsJ.et al. (2015). Analysis of breath specimens for biomarkers of plasmodium falciparum infection. J. Infect. Dis.212, 1120–1128. doi: 10.1093/infdis/jiv176

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BjørnstadO. N.NelsonW. A.TobinP. C. (2016). Developmental synchrony in multivoltine insects: generation separation versus smearing. Popul. Ecol.58, 479–491. doi: 10.1007/s10144-016-0564-z

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  CarvalhoB. O.LopesS. C.NogueiraP. A.OrlandiP. P.BargieriD. Y.BlancoY. C.et al. (2010). On the cytoadhesion of Plasmodium vivax-infected erythrocytes. J. Infect. Dis.202, 638–647. doi: 10.1086/654815

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  CheesmanS. J.de RoodeJ. C.ReadA. F.CarterR. (2003). Real-time quantitative PCR for analysis of genetically mixed infections of malaria parasites: technique validation and applications. Mol. Biochem. Parasitol.131, 83–91. doi: 10.1016/S0166-6851(03)00195-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  CiuffredaL.ZoikuF. K.QuashieN. B.Ranford-CartwrightL. C. (2020). Estimation of parasite age and synchrony status in Plasmodium falciparum infections. Sci. Rep.10, 1–10. doi: 10.1038/s41598-020-67817-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  CostaG.LoizonS.GuenotM.MocanI.HalaryF.de Saint-BasileG.et al. (2011). Control of Plasmodium falciparum erythrocytic cycle: γδ T cells target the red blood cell-invasive merozoites. Blood118, 6952–6962. doi: 10.1182/blood-2011-08-376111

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  DeanL. (2005). “Blood groups and red cell antigens,” in Blood and the cells it contains (Bethesda (MD): National Center for Biotechnology Information (US)), 1–8http://www.ncbi.nlm.nih.gov/books/NBK2263/

  • Google Scholar
  • View reference in article

• 9

  DeharoE.CoquelinF.ChabaudA. G.LandauI. (1996). The erythrocytic schizogony of two synchronized strains of Plasmodium berghei, NK65 and ANKA, in normocytes and reticulocytes. Parasitol. Res.82, 178–182. doi: 10.1007/s004360050091

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  DeharoE.GautretP.GinsburgH.ChabaudA. G.LandauI. (1994). Synchronization of Plasmodium yoelii Nigeriensis and P. y. killicki infection in the mouse by means of Percoll-glucose gradient stage fractionation: determination of the duration of the schizogonic cycle. Parasitol. Res.80, 159–164. doi: 10.1007/BF00933785

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  DondorpA. M.DesakornV.PongtavornpinyoW.SahassanandaD.SilamutK.ChotivanichK.et al. (2005). Estimation of the total parasite biomass in acute falciparum malaria from plasma PfHRP2. PloS Med.2, 0788–0797. doi: 10.1371/journal.pmed.0020204

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  DuncanC. J. A.DraperS. J. (2012). Controlled human blood stage malaria infection: current status and potential applications. Am. J. Trop. Med. Hygiene86, 561–565. doi: 10.4269/ajtmh.2012.11-0504

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  FärnertA.SnounouG.RoothI.BjörkmanA. (1997). Daily dynamics of Plasmodium falciparum subpopulations in asymptomatic children in a holoendemic area. Am. J. Trop. Med. Hygiene56, 538–547. doi: 10.4269/ajtmh.1997.56.538

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  Franke-FayardB.JanseC. J.Cunha-RodriguesM.RamesarJ.BüscherP.QueI.et al. (2005). Murine malaria parasite sequestration: CD36 is the major receptor, but cerebral pathology is unlinked to sequestration. PNAS102, 11468–11473. doi: 10.1073/pnas.0503386102

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  GnangnonB.DuraisinghM. T.BuckeeC. O. (2021). Deconstructing the parasite multiplication rate of Plasmodium falciparum. Trends Parasitol. 37, 922–932. doi: 10.1016/j.pt.2021.05.001

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  GreischarM. A.ChildsL. M. (2023). Extraordinary parasite multiplication rates in human malaria infections. Trends Parasitol.39, 626–637. doi: 10.1016/j.pt.2023.05.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  GreischarM. A.MideoN.ReadA. F.BjørnstadO. N. (2016). Predicting optimal transmission investment in malaria parasites. Evolution70, 1542–1558. doi: 10.1111/evo.12969

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  GreischarM. A.ReadA. F.BjørnstadO. N. (2014). Synchrony in malaria infections: how intensifying within-host competition can be adaptive. Am. Nat.183, E36–E48. doi: 10.1086/674357

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  GreischarM. A.ReeceS. E.SavillN. J.MideoN. (2019). The challenge of quantifying synchrony in malaria parasites. Trends Parasitol.35, 1–15. doi: 10.1016/j.pt.2019.03.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  HastingsA. (1997). Population biology: concepts and models (New York: Springer-Verlag).

  • Google Scholar
  • View reference in article

• 21

  HawkingF. (1970). The clock of the malaria parasite. Sci. Am.222, 123–131. doi: 10.1038/scientificamerican0670-123

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  HawkingF.WormsM. J.GammageK. (1968). 24- and 48-hour cycles of malaria parasites in the blood; their purpose, production and control. Trans. R. Soc. Trop. Med. Hygiene62, 731–760. doi: 10.1016/0035-9203(68)90001-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  HernándezC. M.EllnerS. P.AdlerP. B.HookerG.SnyderR. E. (2023). An exact version of Life Tabl Response Experiment analysis, and the R package exactLTRE. Methods Ecol. Evol.2023, 1–13. doi: 10.1111/2041-210X.14065

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  HuijbenS.NelsonW. A.WargoA. R.SimD. G.DrewD. R.Read (2010). Chemotherapy, within-host ecology and the fitness of drug-resistant malaria parasites. Evolution64, 2952–2968. doi: 10.1111/j.1558-5646.2010.01068.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  KhouryD. S.CromerD.AkterJ.SebinaI.ElliottT.ThomasB. S.et al. (2017). Host-mediated impairment of parasite maturation during blood-stage Plasmodium infection. PNAS114, 2–7. doi: 10.1073/pnas.1618939114

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  KhouryD. S.CromerD.BestS. E.JamesK. R.KimP. S.EngwerdaC. R.et al. (2014). Effect of mature blood-stage Plasmodium parasite sequestration on pathogen biomass in mathematical and in vivo models of malaria. Infect. Imm.82, 212–220. doi: 10.1128/IAI.00705-13

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  KitchenS. F. (1949). “Symptomatology: general considerations,” in Malariology–A comprehensive survey of all aspects of this group of diseases from a global standpoint, vol. Vol. 240 . Ed. BoydM. F. (W.B. Saunders Company, Philadelphia), 966–994.

  • Google Scholar
  • View reference in article

• 28

  LambrosC.VanderbergJ. P. (1979). Synchronization of Plasmodium falciparum erythrocytic stages in culture. J. Parasitol.65, 418–420. doi: 10.2307/3280287

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  LeslieP. (1945). On the use of matrices in certain population mathematics. Biometrika33, 183–212. doi: 10.1093/biomet/33.3.183

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  McKenzieF. E.SmithD. L.O’MearaW. P.RileyE. M. (2008). Strain theory of malaria: the first 50 years. Adv. Parasitol.66, 1–46. doi: 10.1016/S0065-308X(08)00201-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  MideoN.BarclayV. C.ChanB. H. K.SavillN. J.ReadA. F.DayT. (2008). Understanding and predicting strain-specific patterns of pathogenesis in the rodent malaria Plasmodium chabaudi. Am. Nat.172, 214–238. doi: 10.1086/591684

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  MideoN.ReeceS. E.SmithA. L.MetcalfC. J. E. (2013). The Cinderella Syndrome: Why do malaria-infected cells burst at midnight? Trends Parasitol.29, 10–16. doi: 10.1016/j.pt.2012.10.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  MillerL. H.GoodM. F.MilonG. (1994). Malaria pathogenesis. Science264, 1878–1883. doi: 10.1126/science.8009217

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  MillerM. R.RåbergL.ReadA. F.SavillN. J. (2010). Quantitative analysis of immune response and erythropoiesis during rodent malarial infection. PloS Comput. Biol.6, e1000946. doi: 10.1371/journal.pcbi.1000946

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  NeaguI. A.OlejarzJ.FreemanM.RosenbloomD. I.NowakM. A.HillA. L. (2018). Life cycle synchronization is a viral drug resistance mechanism. PloS Comput. Biol.14, 1–26. doi: 10.1371/journal.pcbi.1005947

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  NelsonW. A.BjørnstadO. N.YamanakaT. (2013). Recurrent insect outbreaks caused by temperature-driven changes in system stability. Science341, 796–799. doi: 10.1126/science.1238477

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  NerlichA. G.SchrautB.DittrichS.JelinekT.ZinkA. R. (2008). Plasmodium falciparum in ancient Egypt. Emerging Infect. Dis.14, 1317–1318. doi: 10.3201/eid1408.080235

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  O’DonnellA. J.GreischarM. A.ReeceS. E. (2022). Mistimed malaria parasites re-synchronize with host feeding-fasting rhythms by shortening the duration of intra-erythrocytic development. Parasite Immunol.44, e12898. doi: 10.1111/pim.12898

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  O’DonnellA. J.SchneiderP.McWattersH. G.ReeceS. E. (2011). Fitness costs of disrupting circadian rhythms in malaria parasites. Proc. R. Soc. B: Biol. Sci.278, 2429–2436. doi: 10.1098/rspb.2010.2457

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  OwolabiA. T. Y.ReeceS. E.SchneiderP. (2021). Daily rhythms of both host and parasite affect antimalarial drug efficacy. Evolution Medicine Public Health9, 208–219doi: 10.1093/emph/eoab013

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  PaulR. E. L.ArieyF.RobertV. (2003). The evolutionary ecology of Plasmodium. Ecol. Lett.6, 866–880. doi: 10.1046/j.1461-0248.2003.00509.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  PriorK. F.Rijo-FerreiraF.AssisP. A.HirakoI. C.WeaverD. R.GazzinelliR. T.et al. (2020). Periodic parasites and daily host rhythms. Cell Host Microbe27, 176–187. doi: 10.1016/j.chom.2020.01.005

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  PriorK. F.van der VeenD. R.O’DonnellA. J.CumnockK.SchneiderD.PainA.et al. (2018). Timing of host feeding drives rhythms in parasite replication. PloS Pathog.14, 1–21. doi: 10.1371/journal.ppat.1006900

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  R Core Team (2020). R: A language and environment for statistical computing (Vienna, Austria: R Foundation for Statistical Computing).

  • Google Scholar
  • View reference in article

• 45

  ReillyH. B.WangH.SteuterJ. A.MarxA. M.FerdigM. T. (2007). Quantitative dissection of clone-specific growth rates in cultured malaria parasites. Int. J. Parasitol.37, 1599–1607. doi: 10.1016/j.ijpara.2007.05.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  Reilly AyalaH. B.WackerM. A.SiwoG.FerdigM. T. (2010). Quantitative trait loci mapping reveals candidate pathways regulating cell cycle duration in Plasmodium falciparum. BMC Genomics11, 577–590. doi: 10.1186/1471-2164-11-577

  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  Salguero-GomezR.JonesO. R.ArcherC. R.BeinC.de BuhrH.FarackC.et al. (2016). COMADRE: A global data base of animal demography. J. Anim. Ecol.85, 371–384. doi: 10.1111/1365-2656.12482

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  Salguero-GomezR.JonesO. R.ArcherC. R.BuckleyY. M.Che-CastaldoJ.CaswellH.et al. (2015). The compadre Plant Matrix Database: An open online repository for plant demography. J. Ecol.103, 202–218. doi: 10.1111/1365-2745.12334

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  SchneiderP.RundS. S.SmithN. L.PriorK. F.O’DonnellA. J.ReeceS. E. (2018). Adaptive periodicity in the infectivity of malaria gametocytes to mosquitoes. Proc. R Soc. B285, 20181876. doi: 10.1098/rspb.2018.1876

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  SimpsonJ. A.AaronsL.CollinsW. E.JefferyG. M.WhiteN. J. (2002). Population dynamics of untreated Plasmodium falciparum malaria within the adult human host during the expansion phase of the infection. Parasitology124, 247–263. doi: 10.1017/S0031182001001202

  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  ter KuileF.WhiteN. J.HollowayP.PasvolG.KrishnaS. (1993). Plasmodium falciparum: in vitro studies of the pharmacodynamic properties of drugs used for the treatment of severe malaria. Exp. Parasitol.76, 85–95. doi: 10.1006/expr.1993.1010

  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  Touré-NdouoF. S.Zang-EdouE. S.BisvigouU.Mezui-Me-NdongJ. (2009). Relationship between in vivo synchronicity of Plasmodium falciparum and allelic diversity. Parasitol. Int.58, 390–393. doi: 10.1016/j.parint.2009.07.011

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  WackerM. A.TurnbullL. B.WalkerL. A.MountM. C.FerdigM. T. (2012). Quantification of multiple infections of Plasmodium falciparum in vitro. Malar J.11, 180. doi: 10.1186/1475-2875-11-180

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  WhiteN. J.ChapmanD.WattG. (1992). The effects of multiplication and synchronicity on the vascular distribution of parasites in falciparum malaria. Trans. R. Soc. Trop. Med. Hygiene86, 590–597. doi: 10.1016/0035-9203(92)90141-X

  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  YamanakaT.NelsonW. A.UchimuraK.BjørnstadO. N. (2012). Generation separation in simple structured life cycles: models and 48 years of field data on a tea tortrix moth. Am. Nat.179, 95–109. doi: 10.1086/663201

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  YayonA.Vande WaaJ. A.YayonM.GearyT. G.JensenJ. B. (1983). Stage-dependent effects of chloroquine on Plasmodium falciparum in vitro. J. Protozoology30, 642–647. doi: 10.1111/j.1550-7408.1983.tb05336.x

  • CrossRef
  • Google Scholar
  • View reference in article