Reconstruction of cell population dynamics using CFSE
- Andrew Yates†^{1}Email author,
- Cliburn Chan†^{2}Email author,
- Jessica Strid^{3},
- Simon Moon^{4, 5},
- Robin Callard^{6},
- Andrew JT George^{7} and
- Jaroslav Stark^{4, 5}
DOI: 10.1186/1471-2105-8-196
© Yates et al; licensee BioMed Central Ltd. 2007
Received: 25 September 2006
Accepted: 12 June 2007
Published: 12 June 2007
Abstract
Background
Quantifying cell division and death is central to many studies in the biological sciences. The fluorescent dye CFSE allows the tracking of cell division in vitro and in vivo and provides a rich source of information with which to test models of cell kinetics. Cell division and death have a stochastic component at the single-cell level, and the probabilities of these occurring in any given time interval may also undergo systematic variation at a population level. This gives rise to heterogeneity in proliferating cell populations. Branching processes provide a natural means of describing this behaviour.
Results
We present a likelihood-based method for estimating the parameters of branching process models of cell kinetics using CFSE-labeling experiments, and demonstrate its validity using synthetic and experimental datasets. Performing inference and model comparison with real CFSE data presents some statistical problems and we suggest methods of dealing with them.
Conclusion
The approach we describe here can be used to recover the (potentially variable) division and death rates of any cell population for which division tracking information is available.
Background
Quantifying the dynamics of cell populations involves measuring rates of division and death. On a practical level, knowledge of these rates can be important for the clinical assessment of diseases characterised by dysregulated cell populations such as neoplasias. Perhaps more fundamentally, quantifying cell dynamics is important for testing hypotheses regarding the population biology of cells.
Studies of cell proliferation have benefited in recent years from the development of a method to measure the number of divisions single cells have undergone using CFSE (Carboxy Fluoroscein Succinimidyl Ester), a fluorescent and cell-membrane impermeable dye. CFSE is now used widely in immunology to study lymphocyte dynamics [1] but also in oncology [2], stem cell research [3, 4] and to study the kinetics of bacterial division [5]. Briefly, the procedure is as follows. A population of cells is stained with CFSE, and the dye contained in each cell is shared approximately equally among daughter cells upon division. The fluorescence intensities of the population of CFSE-labeled cells can then be measured at a later time using flow cytometry. Cohorts of cells that have undergone the same number of divisions are usually observed to have approximately log-normally distributed intensities, with median decreasing roughly two-fold with each division. Analysis of CFSE profiles allows the estimation of the proportions of cells in culture that are in each generation. These proportions can indicate the extent of division in a population, but CFSE information can also be used to simultaneously quantify division and death if the total numbers of live cells in each generation are known at two or more timepoints. In in vitro experiments, these can be estimated by adding known numbers of fluorescent beads to the culture, sampling from it, counting both cells and beads in the sample using flow cytometry and scaling the generation proportions appropriately.
The information CFSE provides regarding this generational structure augments methods of pulse-labelling with markers such as BrDU (5-bromo-2'-deoxyuridine) or tritiated thymidine, which have traditionally been used to quantify proliferation. These compounds are taken up during DNA synthesis and allow the measurement of the proportion of the population undergoing mitosis during the labelling period. This technique has been used in conjuction with mathematical models to quantify the turnover of populations that are essentially homogeneous (see, for example, [6]). Models have been used to quantify turnover from CFSE data in similar situations [7–11]. In these studies, all cells are considered to be identical, and death or entry into division are represented as Poisson processes. ODEs are usually used, providing the expected numbers of cells in each division. While these models are useful as a starting point, in their simplest form they allow for arbitrarily short inter-division times. This is a biologically unrealistic artifact which can lead to difficulties in the interpretation of estimates of average division and death rates [12]. Other CFSE modeling studies have overcome this by turning to the classic Smith-Martin model of the cell cycle [13]. In this model cells are assumed to spend exponentially-distributed times in a quiescent A-phase before progressing deterministically through an 'actively dividing' B-phase (roughly corresponding to DNA synthesis and mitosis) of finite duration. However, if different susceptibilities to death are allowed in the two phases, as might reasonably be expected given the metabolic differences between quiescence and mitosis, it has been shown that CFSE data alone is not sufficient to identify all parameters of the general Smith-Martin model [9, 10], and additional information (such as the proportion of cells in each generation that are in the A- and B- phases) is required.
As a further complication, it has increasingly been recognised that rates of division and death are usually not homogeneous, and that it is essential to consider this if CFSE is to be used as a practical tool for studying cell dynamics in any depth. Rates of division and death typically vary systematically at a population level. This variation might occur with the number of divisions a cell has undergone; with time, for example as the availability of nutrients, inter-cellular signalling molecules or pro- or anti-apoptotic factors changes over the course of an experiment; or both. Some of these issues were tackled in a series of elegant studies by Gett and Hodgkin [14], Deenick et al. [15], and the subsequent extension of their analysis by de Boer and colleagues [12, 16]. They quantified the kinetics of in vitro stimulation of CFSE-labeled T cells, using a hybrid model in which entry into the first division is stochastic and subsequent divisions are deterministic. They discuss the estimation of the distribution of entry times into the first division, and showed a significant improvement in fit using a division-dependent death rate. Towards a more general approach, Leon et al. [17] proposed a framework for modeling asynchronous division with CFSE data and used this to determine the parameters of probability distributions of inter-division times, allowing for heterogeneity in cell kinetics with respect to division history. However, their analytic approach and the lack of treatment of the sources of discrepancy between model and data make the fitting and comparison of models difficult, and so limits its practical usefulness.
The method we present here has at least two advantages over existing approaches. Firstly, in many cases even time-series of CFSE data may be insufficient to identify the parameters of more detailed models of cell division, and in some cases (as in the general Smith-Martin model discussed above) unique identification of all parameters with CFSE alone is not possible. In contrast, branching processes make minimal assumptions regarding the cell cycle – essentially, the finite timestep imposes a lower bound on the time required to complete a division – and in general all of their parameters are identifiable. In particular this allows useful dynamical information to be recovered even from limited CFSE datasets, such as a single timepoint. Secondly, the inference procedure we propose provides a statistically sound basis for model fitting. Many studies (implicitly) ascribe the discrepancies between the model and the counts of cells in each generation recovered from CFSE profiles as measurement error terms of constant variance. In this paper we challenge this assumption and use a standard stochastic description of cell population dynamics, along with a more realistic treatment of the sources of discrepancy between model and data, to provide the appropriate weighting to each observation when fitting models. Specifically, when estimating parameters of stochastic models from data it is important to assess the relative contributions of fluctuations arising from the intrinsically probabilistic nature of cell dynamics and measurement error or other forms of experimental noise. In this paper we describe two frameworks for parameter estimation; one when fluctuations are the most important form of discrepancy between model and data, and the other when other forms of measurement error dominate. In the latter case, the procedure we describe in this paper can be applied to any model used to describe CFSE data that provides the expected cell counts in each generation.
Using a likelihood-based estimation method requires calculating the probability (likelihood) of a set of observations arising given a model. The generating-function approach we describe allows us in principle to write an exact likelihood given a specification of a branching process model, initial cell numbers, and experimentally observed cell counts at one or more timepoints. However, this method becomes impractical when used with more than a few cells or one or two cell divisions, and is essentially impossible to apply to experimental situations which involve typically tens of thousands of cells. We propose a solution to this problem with the use of a Quasi-Likelihood estimation method. This requires only the first two moments of the probability distribution of the total numbers of cells in each generation – that is, their expectation values and their variance-covariance matrix. We will show that this key simplification allows the model parameters to be inferred from CFSE information.
Results
In Section 1 we describe the theory underlying the parameter estimation and in Section 2 we validate it using synthetic datasets. In Section 3 we describe how to deal with statistical issues that may arise with the application of the method to experimental data, and illustrate this with an analysis of data from an in vitro T cell proliferation experiment.
1. Cell kinetics as a branching process
Calculating the probability distribution of cell counts
To apply a maximum likelihood method to estimate parameters of a stochastic model of cell division and death from CFSE data, we need to characterise the probability distribution of cell counts predicted by the model. In this section we outline this calculation for a general branching process model in discrete time, or a Galton-Watson process [23].
In these models, during each timestep a cell can do one of the following: divide, with probability γ; survive without dividing, with probability δ; or die, with probability 1 - γ - δ (Figure 1). A particular model of the kinetics of a cell population specifies these probabilities, which in the simplest case might be assumed to be constant. In general they may depend on either the number of divisions the cell has undergone (which we refer to as the generation number), explicitly on time, or both. The key assumptions are that all cells act independently, their offspring generate their own branching processes according to the same rules, and that cells retain no memory of events in previous timesteps other than the total number of divisions they have undergone.
The parameters of biological interest are usually γ and α (the probabilities of division and death). However, in the formalism we use here it proves simpler to work with the quantities γ and δ (the probability of survival without division). The probability of death α can then be calculated from 1 - γ - δ. A particular branching process model of cell division is specified by a choice of timestep, a starting condition – the number of cells in each generation at a given time, usually all in generation 0 – and a set of parameters that determine the probabilities γ_{ i }(t) and δ_{ i }(t) for each generation at each subsequent timestep.
Let the state of the cell population at timestep t be the vector ${Z}_{t}=\left({Z}_{t}^{0},{Z}_{t}^{1},\mathrm{...}{Z}_{t}^{n}\right)$, where the components ${Z}_{t}^{i}$ are random variables that represent the number of live cells that have divided i times. The maximum division number n is chosen to be at the limit of detectability on a CFSE profile, or the maximum division number of interest. Given a model and a dataset consisting of the cell counts in each generation at two or more timepoints, we wish to estimate the model parameters. To do this we use the data and the joint probability distribution of Z_{ t }at each timepoint to construct a likelihood. Maximising this with respect to the model parameters and the timestep provides us with best-fit estimates.
We use a probability-generating function (pgf) approach, described in detail in Methods, which allows us to calculate the moments of the distribution of cell numbers in each generation at one timestep given knowledge of their numbers in the previous timestep. Derivatives of the pgf are used to construct a transition matrix M which maps a measured set of cell counts Z_{ t }to their expected values E (Z_{t+1}) at the following timestep. For stationary (time-independent) parameters, we show in the Methods section that given any set of initial cell counts ${Z}_{0}=\left({Z}_{0}^{0},{Z}_{0}^{1},\mathrm{...}{Z}_{0}^{n}\right)$
E (Z_{ t }|Z_{0}) = Z_{ 0 }M^{ t },
and the entries in M are the probabilities of a cell in generation i dividing (γ_{ i }) or surviving without dividing (δ_{ i }), and γ_{ i }+ δ_{ i }≤ 1. Typically an experiment begins with a population of undivided cells and so Z_{0} = (N_{0}, 0, ..., 0).
This stochastic approach also provides the covariance matrix of cell counts in each generation at time t, V_{ t }, in terms of Z_{0}, the E (Z_{ t }) and M (see Methods). The framework is easily extended to calculate the quantities E (Z_{ t }) and V_{ t }when the parameters governing cell kinetics are also functions of time. In the analyses we present below, we used Mathematica [24] to generate E (Z_{ t }) and V_{ t }given initial cell counts Z_{0} and a set of parameters that specify a branching process model – i.e., how the probabilities γ and δ vary with division and/or time.
This approach can also be applied to a qualitatively different class of models, Markovian branching processes in continuous time. In these models cells have exponentially distributed lifetimes, at the end of which they either divide or die. We describe this in Appendix 1. Indeed the method we discuss in the following section applies to any stochastic model which provides the quantities E (Z_{ t }) and V_{ t }given a set of initial cell counts Z_{ 0 }.
Parameter estimation using quasi-likelihood
In principle a likelihood can be computed exactly for any branching process and a dataset. While this is feasible for small cell populations or one or two divisions, with the cell numbers encountered in most experimental situations this becomes intractable for combinatorical reasons (see Appendix 2 for a discussion). As a solution, we take a Quasi Likelihood (QL) approach which requires only the first two moments of the cell counts [25]. QL yields consistent parameter estimates, (that is, the estimates converge to their true values for large sample sizes or large numbers of cells) with minimal confidence intervals [26]. Given the large numbers of cells typically observed in experiments, one might intuitively expect that by the central limit theorem the distribution of cell counts might be well specified by their means and covariances alone.
We find convergence with this algorithm is robust to the choice of initial guess. To speed convergence, particularly with complex models, we select an initial condition by randomly generating a large sample of candidate parameter vectors and choose the one that maximises the likelihood as defined in the following section.
and use these in the expressions above, with the parameter vector β now including the independent unknowns p^{(1)}, ..., p^{(m - 1)}.
The covariance matrix of the parameter estimates cov (β*) is asymptotically the inverse of the information matrix i (β). Since U is (asymptotically) the derivative of a log likelihood, i^{-1} (β) is an estimate of the curvature of the log likelihood surface in parameter space. This provides confidence intervals directly if we assume no error in the cell counts Y_{ t }– that is, if all uncertainty in our parameter estimates comes from the underlying stochasticity of cell behaviour expressed by the model. These confidence intervals are typically rather small given the large numbers of cells usually observed in proliferation assays.
We also note that when the observations are generated by a true branching process the weighting to datapoints provided by the covariance structure is not required for generating point estimates of parameters, since the fitting procedure is essentially a minimisation of a sum of squared residuals, each of which is non-negative and is strictly zero (along with the score function) at the QL estimate of the parameters. The covariance structure is important, however, for the correct estimation of confidence intervals on branching process parameters using the information matrix, and for model discrimination using likelihood ratio tests (see below).
A Mathematica notebook which implements the calculation of the mean and covariances of the cell counts, the generation of the initial parameter estimate and the QL estimation procedure is available on request from the authors (AY and CC).
Model comparison
Typically there may be several candidate branching process models that might describe the biology and we want to assess the relative support for each. Again, assuming no measurement error in the observed cell counts Y_{ t }, the usual procedure for comparing two nested models A and B, A with n additional parameters is to use the residual deviance [25], defined as twice the difference between the maximum achievable log likelihood given the data and the log likelihood at the QL estimate of the parameters -
D (Y; μ) = 2 L (Y; Y) - 2 L (Y; μ),
where L (Y; μ) is the logarithm of the likelihood of a model with expected cell counts μ generating the observations Y. The quantity D_{ A }– D_{ B }for models A and B is asymptotically χ^{2}-distributed with n degrees of freedom. This is the standard likelihood ratio test.
This is essentially a multivariate normal approximation to the true log likelihood.
To compare non-nested models, the simplest approach is to compare the absolute values of likelihoods (see, for example, [20]) or to use the Akaike Information Criterion. This is necessary when comparing the fits with different timesteps, of which there are usually a restricted set of discrete choices; these are dictated by the maximum division number observed at each timepoint, and the intervals between these timepoints. It can also be used to compare members of a family of models with the same number of parameters – for example, when division or death probabilities are assumed to change at a given, but unknown, division number.
2. Validation of the method
Testing the validity of the QL estimator
Validation of the method in the presence of measurement noise
We make three simple observations here. First, the uncertainty in parameters scales approximately linearly with the amplitude of the noise, and a given fractional uncertainty σ in cell counts translates into a comparable fractional uncertainty in parameter estimates. Second, the division probabilities strongly in fluence the shape of the CFSE profile and so in general are estimated more accurately when total counts are subject to noise than when cell counts in each generation are subject to independent error. Third, the division and death probabilities that apply to more CFSE peaks or measurements (in this example, γ_{1+} and α_{1+}, which determine the division and death probabilities for all cells in generations 1 and above) can be estimated more accurately than those constrained by fewer measurements (here, γ_{0} and α_{0} for undivided cells). This effect is again more pronounced when the proportions of cells in each generation are known more accurately than the total numbers.
Relation of parameters to more complex models
The parameters in the branching process (BP) and Smith-Martin (SM) models can be related with some approximations. In this instance of the SM model the probability of a cell dying during a finite interval τ, the branching process parameter α, is independent of the cell being in the A or B phase and so we predict that the QL estimate α should be given by
α = 1 – e^{-μτ}.
The QL procedure identifies the homogeneous model correctly and the estimated death probability α agrees closely with the predicted value for all division rates. The QL estimate of division probability γ agrees well with the predicted value (6) when the SM division rate λ is low, but the two diverge as λ increases. The discrete time process does not specify the true (continuous) distribution of interdivision times, but instead 'coarse-grains' this distribution by allowing division at any time within each timestep. For constant probabilities of division and death, this generates a geometric distribution in discrete time, such that (in the absence of cell death) the probability that a given cell observed since t = 0 divides during the interval t' = nτ and t = (n + 1) τ is P (n) = γ (1 - γ)^{ n }; while for the SM model with constant parameters the probability density for the interdivision time t, P (t), is exponential with a delay, or P (t) = 0 for 0 ≤ t ≤ Δ and P (t) = λ exp (-λ (t - Δ)) for t > Δ. These distributions converge for t = nτ when division rates are low; that is, when the timestep τ is smaller than the average time spent in the A-phase (τ << 1/λ) and when the average time spent in the A-phase is much longer than the B-phase (11/λ >> Δ).
3. Dealing with experimental CFSE data
An important issue when quantifying the dynamics of CFSE-labeled cells is assessing our confidence in the observed cell counts Y_{ t }. In this section we discuss how to deal with various sources of uncertainty in the cell counts and how these impact on model fitting and comparison. Another significant source of disagreement between model and observations, of course, is that the underlying model may not represent the biology well. With this in mind, what we discuss here applies not only to the discrete time branching models we describe here but also to any stochastic model of cell division that can be used to provide likelihood-based parameter estimates.
Uncertainties in the assignment of cells to generations from CFSE profiles
The process of assigning a division number to cells in a CFSE profile can be a significant source of error, particularly if the peaks corresponding to cells in one generation are ill-defined. The distributions of neighbouring peaks usually overlap significantly, and cells in the tails of these distributions may be mis-assigned to neighbouring generations. Further, the factor difference in median fluorescence intensity of adjacent peaks is typically not exactly 2, and this error can amount to uncertainties of as much as a whole division for cells that have divided multiple times. This is particularly noticeable in CFSE profiles which contain distinct subpopulations of cells separated by several divisions and with few cells to mark the location of intermediate generations. In many circumstances, then, the 'gating' or assignment of cells to different divisions is itself a process of inference.
We used a standard algorithm to perform this, based on the Expectation-Maximisation (EM) algorithm [29]. EM is a bounded optimisation technique for the computation of maximum likelihoods typically used in incomplete-data problems. CFSE histograms generated in experiments (i.e., the plot of event counts against the logarithm of fluorescence intensity) can usually be approximated well by normal mixtures (i.e. a superposition of Gaussian distributions) and estimating the parameters for such a normal mixture is a standard application of the EM algorithm. In practice, we find that the algorithm works well only if we provide good initial conditions for the modes (maxima) of each normal component in the mixture, as well as some constraint on the variance of each component. Initial locations for modes are found by first specifying the data range which contains 99% of the total events, then calculating the offset (alignment of entire fit) and stride (the average fold reduction in fluorescent intensity between peaks) that produce the average largest event count. This works well because the inter-peak distances for CFSE profiles tend to be similar, as we would expect if CFSE is equally distributed between daughter cells. As a result, the initial modes are regularly spaced; however, the EM algorithm is then free to adjust the modes to produce the best fit. We heuristically set a constraint such that the variance of each component is less than or equal to that of the component with the tallest peak. Counts are then estimated using the relative area under each normal component scaled by the total number of cells.
- 1.
Use the EM method to identify a maximum-likelihood set of log-normal profiles from a raw CFSE profile containing N_{0} cells. We refer to the resulting set of counts of cells in each generation as Y^{(0)}, where the sum of the elements of Y^{(0)} equals N_{0}.
- 2.
Using Y^{(0)} and a model characterised by a set of parameters β, calculate a best-fit (QL) set of parameter estimates β_{0}.
- 3.
Generate P artificial CFSE profiles, as follows. For each generation or peak k in the original profile, draw ${Y}_{k}^{(0)}$ random numbers from the log-normal probability distribution used to fit that peak. This generates a population of N_{0} cells with fluorescent intensities drawn from the predicted distributions. Use this to re-estimate the numbers of cells in each division using the EM method. Repeat this P times. This generates a set of new, artificial CFSE fluorescence profiles (Y^{(1)}, Y^{(2)}, ..., Y^{(P)}) derived from the original counts Y^{(0)}.
- 4.
For each artificial dataset Y^{(i)}calculate a parameter set estimate β_{ i }.
- 5.
We now have P samples from a probability distribution of parameter estimates representing our uncertainty in the assignment of division numbers to cells in the original CFSE profile. Calculate confidence limits on β_{0} from this distribution.
As noted above, if the procedure provides estimates of the division and quiescence probabilities γ and δ, probabilities of death α can be calculated using α = 1 - γ - δ. It is then straightforward to calculate confidence intervals on α given the distribution of estimates of γ and δ.
We also note that each estimate β_{ i }comes with its own confidence limits, stemming from the stochasticity of the branching process. We thus have at least two independent sources of uncertainty in parameters – one that stems from the uncertainty in the assignment of cells to different generations, which we estimate with the Monte Carlo procedure above; and the other from the underlying stochasticity of the branching processes – that is the range of parameter values that could reasonably (i.e. with some significant probability) have generated each of the datasets (Y^{(0)}, Y^{(1)}, ..., Y^{(P)}).
This procedure assumes high levels of confidence in the measured total cell numbers. If only a single experimental replicate is available, one may have little a priori knowledge of the uncertainty in total cell counts and its effect on parameter estimates. This may be significant in in vitro experiments, but is particularly important when tracking CFSE-labeled cells in vivo. For example, if labeled cells are transferred to an animal and recovered blood and/or lymphoid tissues at a later timepoint, there may be both loss of cells in the recovery procedure as well as uncertainties in the number transferred successfully (e.g. the initial 'take' after intravenous transfer). We suggest that in the absence of experimental replicates, one approach to this problem is to make a heuristic estimate of the error in total counts, and then apply noise at this level to the total cell counts in the Monte Carlo procedure described above. We describe this in the example that follows.
Application to an experimental dataset
To illustrate our method of estimation with branching processes, we apply it to an experimental CFSE dataset (Figure 6). We modeled the response of a polyclonal population of CD8^{+} T cells to stimulation in vitro with anti-CD3 and and anti-CD28 antibodies, in the presence of IL-2 (a growth factor). CFSE profiles from independent cultures were obtained at days 1–4. Little cell death or division was observed in the first 24 h so the 24 h timepoint was taken as the initial condition. The majority of T cells were expected to respond to this stimulus and so we modeled the system as a single population with division or death parameters varying (possibly) with time and/or generation number.
We fitted a variety of models to this data, allowing parameters to vary with time and/or division. The optimal timestep for all models (as measured by the absolute value of the likelihood) was 12 h, and assuming no divisions took place before 36 h. A reasonable fit was obtained with a four-parameter model that allowed undivided cells (generation 0) and divided cells (generations 1+) to have distinct probabilities of division and death; an extension to six parameters allowed different division and death probabilities in generations 0, 1–3 and 4+. The extended model gave a significantly better fit (χ^{2} test on the difference in log likelihoods, on 2 degrees of freedom, p < 10^{-6}). The best fit using the six-parameter model and the corresponding parameter estimates are shown in Figure 6 and Table 2.
These results suggest slow recruitment of undivided cells into division after 36 hours, with a significant probability of apoptosis in the undivided population. Cells that have divided once divide again with approximately 40% probability in each 12 h interval, with increased susceptibility to apoptosis; division slows significantly in the fourth generation. Thus the method identifies the slow first division commonly observed in T cell proliferation assays; it also suggests that cells dividing rapidly have an associated high probability of death.
We quote confidence intervals on the parameter estimates using (i) the asymptotic properties of the QL estimator; (ii) the Monte Carlo (MC) method, taking into account the uncertainty of assigning cells to CFSE peaks, and (iii) the more conservative MC method, applying an additional estimate of measurement error (5% Gaussian noise applied to total cell numbers) to each of the MC replicates. We note that the parameters governing the 4th division are not well constrained as their estimation depends on the single measurement of the cell counts in generation 5 at 96 h.
Comparing models using estimation of measurement error
where the next-to-diagonal elements ρ represent the misassignment between generations, and the diagonal elements σ represent the combination of misassignment and error in total cell counts, if any. The matrix Λ may also be expected to vary between timepoints (Λ = Λ_{ t }). We refer to the parameters that characterize these matrices collectively as η.
where now the sum is over all timepoints t and over all replicates (Monte Carlo or experimental), i.
This quantity can be used directly for model comparison, either with likelihood ratio tests or information criteria statistics such as the AIC [30], although obtaining the estimate of $\mathcal{L}$ by numerical maximisation of (7) may be difficult for complex models. This has the flavour of a mixed-effects approach [31, 32] in which the original 'fixed' effects represented by the quantities μ_{ t }(β) and V_{ t }(β) are augmented by the 'random' effects Λ. The parallel remains to be explored further, however, as in our case random effects are represented at the level of the variance in the predicted response μ rather than in the underlying parameters β as in standard mixed-effects models. The estimation of additional parameters in the variance function has previously been discussed as an extension to the quasi-likelihood method [33], but the non-integrability of the score function that we note above prevents the use of this formalism directly.
Estimation of timestep
A natural choice of timestep for single CFSE measurements is provided by the number of divisions observed during the experiment. That is, if it is clear that cells have undergone at most n divisions times over a time t, this suggests a timestep of t/n.
The timesteps are not required to be of equal length, although dividing the duration of the experiment into equal intervals provides the most intuitive interpretation of the probabilities of division and death per time-step as 'rates' for these processes. The method we use here is to generate a discrete set of candidate timesteps that are consistent with the number of significant peaks observed at each timepoint in a CFSE dataset, and simply search systematically for the combination of model and timestep that maximises the (quasi-) likelihood.
- 1.
Start with a dataset that contains zero cells in generation n + 1.
- 2.
Generate QL parameter estimates using this dataset.
- 3.
Calculate the expected numbers of cells in the 'missing' peak using these parameter values and construct a new dataset with this number of cells in generation n + 1.
- 4.
Repeat steps 2 and 3 until parameter estimates converge.
Discussion
In this paper we have presented a method for fitting and comparing classical branching process models of cell division and death to data from CFSE labeling experiments. All parameters of this class of models can be estimated from CFSE data alone. To do this, we take a Quasi-Likelihood approach, overcoming the problem of non-computability of the exact likelihood. Further, by modeling explicitly the different sources of uncertainty in present in CFSE data, the methods we describe here can improve on existing approaches to estimating parameters, which use least squares fitting with the assumption that cell counts in each generation are subject to errors of constant variance.
Many other approaches to modeling CFSE data characterise the continuous distribution of interdivision times explicitly (exponential for simple ODE models, delayed exponential for Smith-Martin models, lognormal for the first division in the model used by Gett and Hodgkin). In contrast, the discrete-time branching process models we discuss here deal only with the average probabilities of division or death during a finite time interval. While this might be seen as a limitation, in many cases the true distribution of interdivision times may be unknown and the discrete-time approach may provide more robust predictions than with other models. The discrete timestep also provides a lower bound on the time taken for a cell to divide without modeling transit through the cell cycle explicitly, and the parameters of these models are identifiable given sufficient CFSE data. We suggest that the branching process approach is particularly suitable for analysis of data in which prior information regarding division kinetics is limited, as well as providing a method of dealing simultaneously with stochastic fluctuations and measurement errors.
Differences in the expected cell counts predicted by any model and the data, assuming the model is a faithful representation of the cell dynamics, stem from (1) the contributions of stochastic fluctuations (if any) from the model and (2) other forms of experimental noise. In the limit where the contribution of (2) overwhelms that from (1), we suggest the Monte Carlo method described here can be used to estimate confidence intervals on model parameters, and that any covariance structure predicted by the stochastic model can reasonably be neglected and only the expectation values need be used. In this case, proper model comparison using the likelihood, which explicitly contains the weightings (variances) associated with each CFSE peak, becomes very dependent on reasonable estimates being obtained for these weights. These are best estimated simply and empirically with replicate datasets. On the other hand, if cell numbers are small and and measurement errors are smaller than or comparable to fluctuations, or when only a single replicate is available, the covariance structure is important as a basis for inference.
By using knowledge of initial cell numbers and total cell counts at subsequent timepoints, models applied to CFSE data allow the estimation of death rates (averaged over all phases of the cell cycle) without the requirement of an assay for dead cells. This is particularly useful as dead cells do not persist in culture for long, and are particularly difficult to identify in vivo, so direct estimates of their numbers are error-prone. However, a limitation of all current methods of estimating death rates from CFSE alone is that cells may remain CFSE-positive for a short time after death and be counted as live. To improve the reliability of these estimates, for example, cells can be co-stained with Propidium Iodide to exclude those whose DNA content has been degraded.
We note that the moment-based QL estimation procedure can be applied to any stochastic model of cell dynamics which provides a covariance structure, and is not restricted to branching processes. We also emphasise that the Monte Carlo procedure for quantifying errors in the counts derived from CFSE fluorescence profiles can be applied directly to parameter estimation with deterministic models. Whatever description of the dynamics is used, treating the different sources of uncertainty in cell population data in the way we describe here allows us to more carefully test and discriminate between models of cell dynamics.
Methods
Detailed derivation of the moments of the distribution of cell counts
Here we show the calculation of moments of the probability distribution of the cell counts Z_{ t }for a stationary branching process, one in which the probabilities δ_{ i }and γ_{ i }are independent of time. We use a probability-generating function (pgf) approach.
where p_{ i }is the probability that a cell will provide i offspring in the next generation and s is a dummy variable. That is, a cell divides with probability p_{2} = γ to produce two cells; survives without dividing (that is, provides one 'offspring') with probability p_{1} = δ; and dies with probability p_{0} = 1 - γ - δ. The pgf enumerates all the possible outcomes after one timestep, and this is contained in the property f (1) = 1, or ∑p_{ i }= 1.
with a similar expression to eqn. (10) for the variance, and where f^{(t)}(s) is f iterated t times (that is, f^{(t)}(s) = f (f^{(t - 1)}(s)) = f (f (f^{(t - 2)}(s))), etc.)
This models the total number of cells in the population. To keep track of the numbers of cells in each division we need to extend this procedure to a multi-type branching process in which a cell's 'type' or 'generation' is the number of divisions it has undergone, with undivided cells in generation 0. To calculate the probability distribution of cells in each generation after t timesteps requires a pgf that accounts for the type-label now associated with each cell and the probabilities of transition between types or generations. To do this, the pgf and the dummy variable s become vector-valued quantities with number of components equal to the number of cell types – in our case, the number of divisions we wish to follow using CFSE. In addition, this allows us to specify different probabilities of division and death for cells in different generations.
For example, taking the first entry in f, f_{0}, in one timestep a cell in generation 0 produces an expected number of cells in generation 1 of ∂f_{0}/∂s_{1} = 2 γ_{0}, and an expected number of cells in generation 0 of ∂f_{0}/∂s_{0} = δ_{0} (all derivatives evaluated at s= 1). Notice that we assume that cells in generation n simply die or divide further with probability 1 - δ_{ n }. This would correspond, for example, to cells dividing beyond the range of generations of experimental interest, or to their CFSE fluorescence intensity becoming so low that they become indistiguishable from the CFSE-negative population in the culture – typically after 8 or 9 divisions.
The branching process we describe here is 'memoryless' or a discrete-time Markov process with live cells making probabilistic transitions between the n + 1 possible states or generations. The matrix M is related to the transition matrix of this Markov process, but includes not only the transition probabilities per timestep for cells in different generations, but also the expansion in population size associated with division (transition from generation i to i + 1). In other words, it maps the cell counts at one timestep to their expected values at the following timestep. Note that we do not include dead cells as a state here – an advantage of our approach is we do not require assays for dead cells, and so do not include them as an observable in our models.
or in more compact (matrix multiplication) notation
E (Z_{1}/Z_{0}) = Z_{0} M.
After t timesteps, the expectation values and higher moments of the cell counts in each generation can be calculated from the pgf f^{(t)}(s) (eqn. (12)) using the recursive definition [23]
f^{(t)}(s) = f(f^{(t - 1)}(s))
As a consistency check, each component of the pgf at time t must satisfy the property ${f}_{i}^{(t)}$ (s= 1) = 1. Since f(1) = 1 from the definition (12), it follows from (14) that f^{(t)}(1) = 1 as required.
By simple extension, expected cell counts at later timepoints can be calculated with repeated matrix multiplication using M -
E (Z_{ t }/Z_{0}) = Z_{ 0 }M^{ t }.
For a general initial state Z_{0} = (c_{0}, c_{1}, ..., c_{ n }), the assumption of independence of the branching processes initiated by each cell gives V_{1} = ∑_{ i }c_{ i }v_{ i }. Again, a typical CFSE experiment might start with N cells in generation 0, yielding V_{1} = N v_{0}.
Appendices
1. The Continuous-Time Analogue
The continuous-time analogue of a Galton-Watson process is the Markov age dependent process. This is characterised by cells having life spans that are exponentially distributed random variables with parameter λ [21]. This is conceptually a quite different model of cell behaviour to that described above. It can be compared to a limit of the Smith-Martin model in which death can only occur during the B phase (during which cells are actively dividing) and the duration of this phase approaches zero. Thus, for example, it may be a reasonable model for slow homeostatic division in which the average time spent in the cell cycle is negligible compared to the time spent in quiescence.
for k = η. Using the integrating factor e^{ λt }this system of equations can easily be solved by back substitution yielding
F_{ η }= 1 + e^{-λαt}(s_{ η }- 1)
If we start with one undivided cell at time zero the expectation of the normalized cell count E (N_{ k }) for generation k is obtained by differentiating F_{0} with respect to s_{ k }and subsequently setting all s_{ k }= 1. In this simple case the derivatives of F_{0} do not include terms in s_{ k }and this last step can be omitted.
where θ = λ γ t. As before, given an initial state Z_{0} = (c_{0}, c_{1}, ..., c_{ n }) the covariance matrix at time t will be V(t) = ∑_{ i }c_{ i }v_{ i }(t). Further extension of this model can be achieved through altering constraints on λ and γ. In the case of the probability of division γ this parameter can either become a function of generation k or a function of t. In the latter case the system of equations may become inhomogeneous with respect to time and therefore its solution may prove difficult. A much more general approach to continuous-time models is the use of Bellman-Harris processes where the distribution of lifetimes is not restricted to the exponential. However, many such processes are non-Markovian and so also become significantly harder to analyse.
2. Computing an exact likelihood
Parameter estimates with synthetic data.
Model 1 | ||||
---|---|---|---|---|
Proportion of simulations within | ||||
Par | True | Mean (SD) | 95% CI | 99% CI |
γ _{0} | 0.2 | 0.200 (0.003) | 0.947 | 0.989 |
δ _{0} | 0.7 | 0.700 (0.002) | 0.951 | 0.990 |
γ _{1} | 0.7 | 0.700 (0.006) | 0.949 | 0.990 |
δ _{1} | 0.25 | 0.250 (0.007) | 0.950 | 0.991 |
Model 2 | ||||
Proportion of simulations within | ||||
Par | True | Mean (SD) | 95% CI | 99% CI |
γ _{0} | 0.2 | 0.200 (0.002) | 0.950 | 0.989 |
δ _{0} | 0.7 | 0.700 (0.003) | 0.952 | 0.990 |
γ _{1} | 0.7 | 0.700 (0.005) | 0.951 | 0.989 |
δ _{1} | 0.25 | 0.250 (0.004) | 0.950 | 0.990 |
Model 3 | ||||
Proportion of simulations within | ||||
Par | True | Mean (SD) | 95% CI | 99% CI |
f _{ A } | 0.1 | 0.101 (0.013) | 0.953 | 0.987 |
${\gamma}_{0}^{A}$ | 0.15 | 0.150 (0.041) | 0.943 | 0.984 |
${\delta}_{0}^{A}$ | 0.70 | 0.706 (0.055) | 0.951 | 0.981 |
${\gamma}_{1}^{A}$ | 0.70 | 0.699 (0.020) | 0.951 | 0.990 |
${\delta}_{1}^{A}$ | 0.20 | 0.200 (0.035) | 0.942 | 0.981 |
γ ^{ B } | 0.40 | 0.400 (0.003) | 0.955 | 0.992 |
δ ^{ B } | 0.35 | 0.350 (0.003) | 0.946 | 0.983 |
Parameter estimates for the best fit description of the T cell proliferation data.
95% confidence intervals | ||||
---|---|---|---|---|
Parameter | QL estimate | From QL alone | Monte Carlo with EM | Monte Carlo with EM + 5% noise |
γ _{0} | 0.221 | (0.211, 0.230) | (0.203, 0.339) | (0.181, 0.345) |
α _{0} | 0.232 | (0.222, 0.242) | (0.175, 0.253) | (0.163, 0.283) |
γ _{1–3} | 0.419 | (0.409, 0.430) | (0.365, 0.443) | (0.356, 0.444) |
α _{1–3} | 0.427 | (0.412, 0.442) | (0.379, 0.582) | (0.348, 0.595) |
γ _{4} | 0.086 | (0.077, 0.096) | (0.067, 0.679) | (0.063, 0.478) |
α _{4} | 0.340 | (0.244, 0.437) | (0.027, 0.691) | (0.094, 0.731) |
Notes
Declarations
Acknowledgements
AY was supported by NIH grant RO1 AI 49334 to Rustom Antia, and portions of this study were undertaken while AY was supported by a Wellcome Trust Research Training Fellowship in Mathematical Biology, and by CoMPLEX, University College London. CC was supported by the UK BBSRC. AJTG was a BBSRC Research Development Fellow. SM was funded by an MRC studentship.
Authors’ Affiliations
References
- Lyons AB: Analysing cell division in vivo and in vitro using flow cytometric measurement of CFSE dye dilution. J Immunol Methods. 2000, 243 (1–2): 147-54. 10.1016/S0022-1759(00)00231-3.View ArticlePubMedGoogle Scholar
- Holyoake T, Jiang X, Eaves C, Eaves A: Isolation of a highly quiescent subpopulation of primitive leukemic cells in chronic myeloid leukemia. Blood. 1999, 94 (6): 2056-64.PubMedGoogle Scholar
- Groszer M, Erickson R, Scripture-Adams DD, Lesche R, Trumpp A, Zack JA, Kornblum HI, Liu X, Wu H: Negative regulation of neural stem/progenitor cell proliferation by the Pten tumor suppressor gene in vivo. Science. 2001, 294 (5549): 2186-9. 10.1126/science.1065518.View ArticlePubMedGoogle Scholar
- Prudhomme WA, Duggar KH, Lauffenburger DA: Cell population dynamics model for deconvolution of murine embryonic stem cell self-renewal and differentiation responses to cytokines and extracellular matrix. Biotechnol Bioeng. 2004, 88 (3): 264-72. 10.1002/bit.20244.View ArticlePubMedGoogle Scholar
- Ueckert JE, Nebe von Caron G, Bos AP, ter Steeg PF: Flow cytometric analysis of Lactobacillus plantarum to monitor lag times, cell division and injury. Lett Appl Microbiol. 1997, 25 (4): 295-9. 10.1046/j.1472-765X.1997.00225.x.View ArticlePubMedGoogle Scholar
- Bonhoeffer S, Mohri H, Ho D, Perelson AS: Quantification of cell turnover kinetics using 5-bromo-2'-deoxyuridine. J Immunol. 2000, 164 (10): 5049-54.View ArticlePubMedGoogle Scholar
- Veiga-Fernandes H, Walter U, Bourgeois C, McLean A, Rocha B: Response of naive and memory CD8+ T cells to antigen stimulation in vivo. Nat Immunol. 2000, 1: 47-53. 10.1038/76907.View ArticlePubMedGoogle Scholar
- Bernard S, Pujo-Menjouet L, Mackey MC: Analysis of cell kinetics using a cell division marker: mathematical modeling of experimental data. Biophys J. 2003, 84 (5): 3414-24.PubMed CentralView ArticlePubMedGoogle Scholar
- Pilyugin SS, Ganusov VV, Murali-Krishna K, Ahmed R, Antia R: The rescaling method for quantifying the turnover of cell populations. J Theor Biol. 2003, 225 (2): 275-83. 10.1016/S0022-5193(03)00245-5.View ArticlePubMedGoogle Scholar
- Ganusov VV, Pilyugin SS, de Boer RJ, Murali-Krishna K, Ahmed R, Antia R: Quantifying cell turnover using CFSE data. J Immunol Methods. 2005, 298 (1–2): 183-200. 10.1016/j.jim.2005.01.011.View ArticlePubMedGoogle Scholar
- Asquith B, Debacq C, Florins A, Gillet N, Sanchez-Alcaraz T, Mosley A, Willems L: Quantifying lymphocyte kinetics in vivo using carboxyfluorescein diacetate succinimidyl ester (CFSE). Proc Biol Sci. 2006, 273 (1590): 1165-71. 10.1098/rspb.2005.3432.PubMed CentralView ArticlePubMedGoogle Scholar
- de Boer R, Perelson AS: Estimating division and death rates from CFSE data. J Comp Appl Math. 2005, 184: 140-164. 10.1016/j.cam.2004.08.020.View ArticleGoogle Scholar
- Smith JA, Martin L: Do cells cycle?. Proc Natl Acad Sci USA. 1973, 70 (4): 1263-7. 10.1073/pnas.70.4.1263.PubMed CentralView ArticlePubMedGoogle Scholar
- Gett AV, Hodgkin PD: A cellular calculus for signal integration by T cells. Nat Immunol. 2000, 1 (3): 239-44. 10.1038/79782.View ArticlePubMedGoogle Scholar
- Deenick EK, Gett AV, Hodgkin PD: Stochastic model of T cell proliferation: A calculus revealing IL-2 regulation of precursor frequencies, cell cycle Time, and survival. J Immunol. 2003, 170 (10): 4963-72.View ArticlePubMedGoogle Scholar
- De Boer RJ, Ganusov VV, Milutinovic D, Hodgkin PD, Perelson AS: Estimating lymphocyte division and death rates from CFSE data. Bull Math Biol. 2006, 68 (5): 1011-31. 10.1007/s11538-006-9094-8.View ArticlePubMedGoogle Scholar
- Leon K, Faro J, Carneiro J: A general mathematical framework to model generation structure in a population of asynchronously dividing cells. J Theor Biol. 2004, 229 (4): 455-76. 10.1016/j.jtbi.2004.04.011.View ArticlePubMedGoogle Scholar
- Jagers P: Branching Processes with Biological Applications. 1975, London: WileyGoogle Scholar
- Yakovlev AY, Yanev NM: Transient processes in cell proliferation kinetics. 1989, Springer-VerlagView ArticleGoogle Scholar
- Hardy K, Spanos S, Becker D, Iannelli P, Winston RM, Stark J: From cell death to embryo arrest: mathematical models of human preimplantation embryo development. Proc Natl Acad Sci USA. 2001, 98 (4): 1655-60. 10.1073/pnas.98.4.1655.PubMed CentralView ArticlePubMedGoogle Scholar
- Kimmel M, Axelrod D: Branching Processes in Biology, of Interdisciplinary Applied Mathematics. 2002, Springer, 19:View ArticleGoogle Scholar
- Hyrien O, Mayer-Proschel M, Noble M, Yakovlev A: A stochastic model to analyze clonal data on multi-type cell populations. Biometrics. 2005, 61: 199-207. 10.1111/j.0006-341X.2005.031210.x.View ArticlePubMedGoogle Scholar
- Harris T: The Theory of Branching Processes. 1963, Springer-VerlagView ArticleGoogle Scholar
- Mathematica 5.2, Wolfram Research, Inc. 2005Google Scholar
- McCullagh P, Nelder J: Generalized Linear Models, of Monographs on Statistics and Applied Probability. 1989, Chapman and Hall/CRC, 37:View ArticleGoogle Scholar
- Stuart A, Ord J: Classical Inference and Relationship (Kendall's Advanced Theory of Statistics). 1991, Oxford University Press, 2:Google Scholar
- Li B: A deviance function for the quasi-likelihood method. Biometrika. 1993, 80 (4): 741-753. 10.1093/biomet/80.4.741.View ArticleGoogle Scholar
- Smyth G: Pearson's goodness of fit statistic as a score test statistic. Science and Statistics: A Festschrift for Terry Speed, of IMS Lecture Notes – Monograph Series. Edited by: Goldstein DR. 2003, Institute of Mathematical Statistics, Beachwood, Ohio, 40: 115-126.View ArticleGoogle Scholar
- Dempster AP, Laird NM, Rubin DB: Maximum likelihood from incomplete data via the EM algorithm. Journal Royal Stat Soc B. 1977, 39: 1-38.Google Scholar
- Burnham KP, Anderson DR: Model Selection and Multimodel Inference. 2002, Springer, 2Google Scholar
- Lindstrom MJ, Bates D: Nonlinear Mixed Effects Models for Repeated Measures Data. Biometrics. 1990, 46: 673-687. 10.2307/2532087.View ArticlePubMedGoogle Scholar
- Pinheiro JC, Bates DM: Approximations to the Log-Likelihood Function in the Nonlinear Mixed-Effects Model. Journal of Computational and Graphical Statistics. 1995, 4: 12-35. 10.2307/1390625.Google Scholar
- Nelder J, Pregibon D: An Extended Quasi-Likelihood Function. Biometrika. 1987, 74 (2): 221-232. 10.1093/biomet/74.2.221.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.