How to quantify developmental synchrony in malaria parasites

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.


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 (Toure-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;Toure-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 wellresolved, 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.

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.

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.
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.

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: Time required for 99 % bursting IDC duration : (1) 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 (ŝ P ) 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.

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).
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 Differences in synchrony are dwarfed by population dynamics, even for unrealistically well-sampled time series.We used a range of starting age distributions for the initial cohort of iRBCs (A) and simulated iRBC abundance for P. chabaudi infections of mice (B) and P. falciparum infections in vitro (C), both "sampled" every 15 minutes.The initial age distribution for each cohort is described by a symmetric beta distribution and is assumed to be extremely synchronous (a narrow distribution of parasite ages, s P = 500, orange), highly synchronous (s P = 100, pink), synchronous (s P = 10, purple), or asynchronous (a uniform distribution of parasites ages, s P = 1, blue).Those beta distributions translate into 99% of the initial inoculum bursting within 2, 4, 13 or 24 hours (respectively) for P. chabaudi, or within 4, 9, 26 or 48 hours (respectively) for P. falciparum.
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 q 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 (q 1 = q 2 = q 3 =… = q 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 1 = m 2 = m 3 = … = 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 pt to indicate the numbers of iRBCs in each age class at time t, pt+1 = Lp t : (3) 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, R e 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 p0 , 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, p0 , for a given number of stage compartments (n).The initial I 0 iRBCs are sorted into n stages according to the beta distribution specified by a given shape parameter, ŝ P (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 0 = 100 to avoid errors of precision that might occur with lower iRBC counts.

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).
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 R e values that best match the observed R e values: Step 4. We project iRBC abundance for a large number of starting age distributions p0 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 (ŝ P ) 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 ( R e ) for comparison with observed per capita replication rates (R e ).
Step 6. Identify the initial age distribution associated with the best fitting R e (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 ŝ P .Calculate synchrony according to Equation 1.

Numerical optimization and estimates of synchrony
We use the Nelder-Mead algorithm to locate the value of ŝ P (and the offset parameter) that minimizes the sum of the absolute error between the observed and predicted per capita replication rates (R e and R e , 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 (ŝ P ) 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.

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 The new approach estimates synchrony while accounting for population dynamics.Example iRBC count data were obtained (Step 1) using the initial age distribution shown (s P ) and projecting forward using the B t values shown at left (black, see Figure 2).We calculate PMR values (Step 2a) and per capita replication rates (R e , Step 2b), the trajectory to be fitted.We estimate B t values by linear interpolation of PMRs (Step 3).Using estimated B t values, we guess an initial age distribution and project population dynamics (Step 4) and calculate predicted R e values to compare with observed R e values (Step 5).We repeat steps 4 and 5 and retain the best fit initial distribution (ŝ P ) to calculate synchrony using Equation 1 (Step 6).To conserve space, y-axis labels are given at the top of each column of panels.
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., Toure-Ndouo et al., 2009;Ciuffreda et al., 2020).

3.2
The new approach recovers true differences in synchrony 3.2.1The 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 R e 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- Stage percentages vary with population dynamics, making it difficult to accurately assess and compare synchrony when estimated in this way.The percentage of rings in simulated infections is shown for time points when rings would be expected (i.e., in the first half of each IDC).Ring percentages from unevenly "sampled" infections (open circles in (A, B)) show percentages that deviate substantially from values expected for perfectly synchronous, semi-synchronous, and asynchronous infections, while even sampling (closed squares in (C, D)) can only distinguish between asynchrony and some degree of synchrony.Colors and x-axis indicate the true duration of bursting and hence initial level of synchrony as in Figure 1A.P. chabaudi infections (A, C) were simulated over 16 IDCs, while P. falciparum simulations represent 4 IDCs.Vertical line segments indicate ±10%, roughly the standard error reported for stage percentage data where sequestration is not occuring (as in rag1 −/− mice; 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).
nonidentifiability-was foreshadowed by the extremely similar per capita replication rates (R e values) 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.
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.

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, The new approach estimates the duration of 99% bursting (A, C) and synchrony (B, D) very close to the true values (i.e., the values derived from the Beta distribution used to simulate the "observed" time series), despite differences in IDC length and population dynamics.The time series to which the Leslie matrix model is fitted were also used to generate Figure 4, but we used total iRBC abundance as shown in Figure 1 rather than stage percentages.Fits to simulated P. chabaudi (P.falciparum in vitro) infections are shown in panels (A-D).Open circles denote uneven sampling intervals while closed squares refer to sampling performed at even intervals twice per IDC.Vertical segments indicate the 95% high density region of predicted values; only evenly sampled time series showed any variation in the predicted values.True initial synchrony and duration of bursting is indicated with colors (as in Figure 1) and by position on the x-axis.Note that in addition to better identifiability, fitting to unevenly sampled data provides more accurate estimates of true synchrony, accounting for modest loss of synchrony over simulated infections (Supplementary Figure S7).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.

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., gd 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 With sampling error, it is still possible to discern differences in synchrony, depending on the error distribution.The time series to be fitted are identical to those shown in Figure 1 and used in Figure 5, but with simulated sampling noise added (see Supplementary Methods for details).Violin plots showing the distribution of 20 fits to time series with simulated sampling noise, with a binomial (A, C) or negative binomial (B, D) distribution.When error is binomially distributed (A, B), synchrony differences can be obscured, with asynchronous infections often appearing spuriously synchronous.Synchrony estimates are better behaved when sampling noise follows a negative binomial distribution (B, D).P. falciparum tends toward much lower percent parasitemia, especially in vitro (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.Greischar et al. 10.3389/fmala.2024.1386266Frontiers in Malaria frontiersin.orgagainst 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 ŝ P , 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 antimalarial 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., agespecific 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 (Hernańdez 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;Toure-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 synchronywhich requires robust methods to estimate it-could improve understanding of pathogenesis and epidemiology and inform efforts to intervene at the individual and population level.
FIGURE 2 Leslie matrices project how age distributions change in response to changing population growth.Panel (A) shows a hypothetical example of changing numbers of new iRBCs produced per schizont (B t ) through time.In panel (B), progression through development is shown as a life cycle diagram (left) and, equivalently, as a Leslie matrix (right).Leslie matrices-one for each B t value-are used to project how age distributions change through time (C), as in the example calculation shown at right.The Leslie matrix is denoted L while p indicates the abundance of iRBCs in each developmental class.

TABLE 1
Key terms and definitions used in our metric of synchrony, as applied to malaria parasites.