Numerical modelling of label-structured cell population growth using CFSE distribution data
© Luzyanina et al; licensee BioMed Central Ltd. 2007
Received: 10 April 2007
Accepted: 24 July 2007
Published: 24 July 2007
The flow cytometry analysis of CFSE-labelled cells is currently one of the most informative experimental techniques for studying cell proliferation in immunology. The quantitative interpretation and understanding of such heterogenous cell population data requires the development of distributed parameter mathematical models and computational techniques for data assimilation.
Methods and Results
The mathematical modelling of label-structured cell population dynamics leads to a hyperbolic partial differential equation in one space variable. The model contains fundamental parameters of cell turnover and label dilution that need to be estimated from the flow cytometry data on the kinetics of the CFSE label distribution. To this end a maximum likelihood approach is used. The Lax-Wendroff method is used to solve the corresponding initial-boundary value problem for the model equation. By fitting two original experimental data sets with the model we show its biological consistency and potential for quantitative characterization of the cell division and death rates, treated as continuous functions of the CFSE expression level.
Once the initial distribution of the proliferating cell population with respect to the CFSE intensity is given, the distributed parameter modelling allows one to work directly with the histograms of the CFSE fluorescence without the need to specify the marker ranges. The label-structured model and the elaborated computational approach establish a quantitative basis for more informative interpretation of the flow cytometry CFSE systems.
Understanding the dynamics of cell proliferation, differentiation and death is one of the central problems in immunology . A cell population is an ensemble of individual cells, all of which contribute in a different way to the overall observed behavior . A quantitative characterization of this heterogeneity is provided by flow cytometry. Flow cytometry is a technique based on the use of fluorescence activated cell sorter (FACS) for a quantitative single cell analysis of the suspensions of cells, which are labelled with fluorescent substance(s). Once the labelled cells are run through the cell sorter machine, the computer collects data on the fluorescence intensity for each cell . The FACS is capable of analyzing up to a dozen parameters per cell at rates up to 105 cells per second. Therefore it represents a versatile tool with an enormous potential to describe the complex nature of cell populations .
A thorough interpretation and comprehensive understanding of CFSE-labelled lymphocytes population data requires both the development of quantitatively consistent mathematical models, e.g. based on distributed parameter systems such as hyperbolic partial differential equations, and efficient computational techniques for the solution and identification of these models. The heterogeneity of the dividing cell populations can be described by a wide range of characteristics, e.g. the number of divisions made, the position in the cell cycle, the mass, the label expression, the doubling time, the death rate. The mathematical modelling approaches for the analysis of cell growth from CFSE assay data developed so far consider the cell populations as a mixture of cells which differ only in the mean level of the CFSE expression per generation [7–11]. The cells within each generation (compartment) are assumed to possess the same constant level of CFSE fluorescence which is reduced by a factor of 2 after one division. Most of the models ignore the heterogeneity of cell populations with respect to the division and death rates, except for the naive versus dividing cells. The effect of cell heterogeneity with respect to the division times in the context of CFSE data analysis is explored in . An extended comparative analysis of the existing compartmental models for CFSE-labelled cell growth has recently been presented in . These models, formulated using ordinary or delay differential equations, consider the dynamics of the consecutive generations of dividing cells but not the single cell identity. Hence they can be referred to as unstructured and non-corpuscular, following the definitions in .
Distributed population balance models, which use partial differential equations (PDEs), are regarded as the most general way of describing heterogenous cell systems. Such models are considerably more difficult to analyze mathematically and numerically than their unstructured counterparts. The most extensively studied distributed parameter models for population dynamics are the age-structured models [14–16]. The only example of application of the age-maturity structured model for the CFSE data analysis is presented in . The cell population is considered to be continuously structured with respect to the cell age, but the maturity variable (the CFSE fluorescence) is discrete, i.e., k distinct cell generations are considered, each characterized by some average CFSE fluorescence per cell, M/2 k , with M the initial fluorescence. The division and death rates are assumed to be independent of the maturity and they are estimated by fitting experimental data with the model visually. In general, for cell growth problems the age-structured population models are considered to be of limited practical value due to the fact that the cell age is difficult to measure experimentally .
A class of distributed parameter models for cell populations growth, which allows direct reference to the experimentally measurable properties of cells, is represented by so-called size- or mass-structured cell populations models [4, 18–20]. The terms ”size” and ”mass” refer to any cell property which satisfies a conservation law, e.g. volume, protein content, fluorescence label, etc. A rigorous mathematical analysis of such models was presented in . The mass-structured population balance models are considered to provide a consistent way to estimate the fundamental physiological functions from flow cytometry data in the area of biotechnology [4, 13].
In this study we formulate a one-dimensional first order hyperbolic PDE model for the dynamics of cell populations structured according to the CFSE fluorescence level. This structure variable defines the division age of the cell. We let the fluorescence intensity of the initial cell population and, therefore, of the consecutive generations to range continuously in some interval, thus relaxing a restricting assumption of an equal expression of CFSE by cells which have undergone the same number of divisions.
The proposed CFSE label-structured model potentially has the following advantages with respect to existing compartmental models: (i) it allows one to estimate the turnover parameters directly from the distributions of CFSE-labelled cells followed over time by flow cytometry; (ii) it does not require an ad hoc assumption on the relationship between the label expression level and the number of divisions cells undergone. Notice that this is an important aspect for a long-term follow up of the CFSE-labelled populations as the correspondence between the CFSE intensity range and the division generation can be heavily biased by the overall loss of the label over time and by the initial heterogeneity of the labelled cell population; (iii) it allows to estimate the kinetic parameters of cell proliferation and death as functions of the marker expression level (and hence of the number of cell divisions).
Modelling with hyperbolic PDEs, being used in the context of data-driven parameter identification, presents a significant computational challenge due to the hyperbolic nature of the equations and due to the large size of the discretized problem. To our knowledge, no publicly available software package exists which deals with optimization of hyperbolic PDE models. We estimate the distributed parameters of the proposed model following the maximum likelihood approach and using the direct search Nelder-Mead simplex method applied to a finite dimensional approximation of the original infinite dimensional optimization problem. The initial-boundary value problem is solved with a Matlab program by Shampine , which implements the well established second order Richtmyer's two-step variant of the Lax-Wendroff method. Because this program is fully vectorized, it allows very fast execution, which is otherwise difficult to achieve in Matlab. This is especially important when solving a PDE in an optimization loop. Using two original CFSE data sets, we demonstrate the biological consistency of the proposed label-structured model and compare its predictions with the predictions of the ODE (ordinary differential equation) compartmental model published recently .
The outline of this paper is as follows. In the next section we formulate the label-structured cell populations model. In section ”CFSE data” we describe two original sets of data on in vitro growth of human CFSE-labelled T-lymphocytes and the preprocessing of the corresponding CFSE histograms used in this study. The major aspects and the numerical treatment of the distributed parameter identification problem are presented in sections ”Parameter estimation” and ”Numerical procedure”. Results of the application of the proposed model to the analysis of the turnover parameters of proliferating cells from the CFSE intensity histograms for the two data sets are presented in section ”Applications to CFSE assay”. Here we also compare the performance of the proposed PDE model and the compartmental ODE model. Finally, we discuss the major advantages and the bottlenecks of the proposed approach.
Label-structured cell populations model
In this section we introduce the mathematical model for the dynamics of lymphocyte populations in the CFSE proliferation assay. We consider a population of cells which are structured according to a single variable x that characterizes the CFSE expression level in terms of units of intensity, UI. Therefore the amount of CFSE label is treated as a continuous variable. The state of the population at time t is described by the distribution (density) function n(t, x)(cell/UI), so that the number of cells with the CFSE intensity between x1 and x2 is given by
At the beginning of the follow-up experiment, the lymphocyte population is stained with CFSE giving rise to the initial (starting) distribution of cells with respect to the CFSE fluorescence. The following phenomenological features of the label-structured lymphocyte proliferation have to be taken into account by the model for the dynamics of the distribution of labelled cells ([5–7, 23]):
• During cell division CFSE is partitioned equally between daughter cells;
• Each CFSE division peak represents a cohort of cells that entered their first division at approximately the same time;
As the cells proliferate, the initially bell-shaped distribution of the CFSE fluorescence in the population becomes multimodal, moving over time to lower values of x. The histograms of the CFSE intensity provide profiles for cell divisions;
• As the dividing cell population approaches the autofluorescence level of unlabelled cells, the division peaks start to compress, thus limiting the number of divisions that can be followed. Usually cells are stained to an intensity of about 103 times brighter than their autofluorescence, so that up to 10 divisions can be permitted while maintaining both the parental and the final generation intensities all on scale.
The label-structured cell population behavior can be expressed using a modification of the model proposed originally by Bell & Anderson for size-dependent cell population growth when reproduction occurs by fission into two equal parts . We assume that the physiological parameters of cells (division and death rates) strongly correlate with the label expression level.
Let the initial CFSE distribution of cells at time t0 be given by the density functionn(t0, x) =: n 0 (x), x ∈ [xmin, xmax].
The first equation consists of the following terms:
v(x)∂n(t, x)/∂x, the advection term, describes the natural decay of the CFSE fluorescence intensity of the labelled cells with the rate v(x), UI/hour;
(α(x) + β(x))n(t, x) describes the local disappearance of cells with the CFSE intensity x due to the division associated CFSE dilution and the death with α(x) ≥ 0 and β(x) ≥ 0 being the proliferation and death rates, respectively, both having the same unit 1/hour;
2γα(γx)n(t, γx) represents the birth of two cells due to division of the mother cell with the label intensity γx. The first factor accounts for the doubling of numbers, and the second for the difference by a factor γ in the size of the CFSE intervals to which daughter and mother cells belong. Indeed, those cells which originate from division of cells with CFSE in the range (γx, γ(x + dx)) enter into the range (x, x + dx).
Under the assumption of equal partition of the label between the two daughter cells and no death during the division one expects that γ = 2. This would ensure conservation of CFSE label, similar to the conservation of volume-size [19, 20]. However, we allow the label partitioning parameter γ to take values smaller than 2 so that x <γx ≤ 2x, in order to check the consistency of the assumptions with experimental data.
The above consideration applies to cells with levels of CFSE below the maximal initial staining xmax divided by γ. The population dynamics of the cells with xmax/γ <x ≤ xmax is governed by the second equation of model (2) without the source term. The division, death and transition rates, α(x), β(x) and v(x), of the structured population are assumed to be functions of (i.e., correlate with) the CFSE intensity. The precise dependence on x is not known a priori and will be estimated from the flow cytometry data.
The initial data for model (2) are given by (1) specifying the distribution of cells at time t0. The lack of cells with CFSE intensity above the given maximal value xmax for all t > t0 is taken into account by the boundary conditionn(t, xmax) = 0, t > t0.
where ν(z) = v(10 z )/log(10)10 z . The structured population balance model (4) is used for the description of the evolution of CFSE histograms and to estimate the division, death and transfer rates of labelled cell populations from CFSE proliferation assays.
CFSE intensity histograms of proliferating cell population
To investigate the appropriateness of the label-structured cell population model (4) and the developed parameter estimation procedure, two original data sets characterizing the evolution of CFSE distribution of proliferating cell cultures were used. The data sets were obtained from in vitro proliferation assay with human peripheral blood mononuclear cells (PBMC) as follows. The cells were labelled with CFSE at day 0. To induce the proliferation of T cells, two different activation stimuli were used:
• the mitogen stimulator phytohemagglutinin (PHA), which activates the T lymphocytes unspecifically, i.e., independent of a signal transduced by the T cell receptor (data set 1, considers the total CD4 and CD8 T cells);
• the antibodies against CD3 and CD28 receptors on T cells which provide signals similar to those transduced by the T cell receptor (data set 2, considers the CD4 T cells).
At regular times after the onset of cell proliferation the cells were harvested, stained with antibodies to CD4 or CD8 and analyzed by flow cytometry for CFSE expression level on individual cells. The total cell number in the proliferation culture was also quantified. The combination of CFSE labelling and flow cytometry allows one to generate the time series of histograms of CFSE distribution .
In a standard approach, the CFSE fluorescence histograms are used to evaluate the fractions of T cells that have completed certain number of divisions [6, 7]. This type of 'mean fluorescence intensity' data can be obtained either manually or by using various deconvolution techniques implemented in programs, such as ModFit (Verity Software), CellQuest (Becton Dickinson), CFSE Modeler (ScienceSpeak). The corresponding computer-based procedures require setting of the spacing between generations, i.e., marking the CFSE fluorescence intensities that separate consecutive generations of dividing cells. Note that when the starting population of cells exhibits a broad range of CFSE fluorescence, the division peaks can be not easily identifiable, making conventional division tracking analysis problematic [3, 23, 25]. The number of divisions which can be followed is limited by the autofluorescence of unlabelled cells. For the data we consider, the resolution of the division peaks is not possible after about 7 division cycles. We present and make use of the division number lumped CFSE distribution data, i.e., 'mean fluorescence intensity', in the last section for comparison of the parameter estimation results for the PDE and ODE based models of cell proliferation.
Preprocessing of CFSE intensity histograms for parameter estimation
A direct use of such fluctuating histogram data for numerical parameter estimation might lead to the following major difficulties: (i) the possibility of overfitting, when the measurement noise rather than the true dynamics is approximated; (ii) the emergence of discontinuities in the computed model solution due to a discontinuous initial cell distribution function, as suggested by the flow cytometry histogram. Overall, for the parameter estimation we need to infer the underlying cell distribution densities n(t i , z) from which the histograms of CFSE counts were sampled. The functional approximation allows one to make predictions about the CFSE-labelled cell density for the z coordinate where cells have not been observed. Because the density distribution is supposed to be a continuous function, the corresponding estimation problem involves some regularization procedure.
Here q defines the ”global” level of smoothing and m i stands for the number of measurements with ni, j> a i in the histogram being smoothed. The performance of the continuous smoothing procedure is presented in Fig. 3 for two choices of the parameter q. Note that a moderate level of smoothing (q = 0.03) preserves important features of the data (the division associated peaks), while q = 0.05 leads to oversmoothing (information loss) as manifested by the disappearance of the division cohort structure presented in the histogram. In our study we used q = 0.03.
The histograms obtained by flow cytometry cover the whole range of the CFSE fluorescence x from 1 to 104. In particular, the starting population of undivided cells can spread up to the upper end of 104units. We did not consider the tiny fraction of cells which differ substantially in their CFSE intensity from the bulk population of homogeneously stained cells. These CFSE bright cells might represent a measurement noise rather than genuine cells as they remain in the same area of the histogram at later observation times. Therefore, for parameter estimation we assumed that there is some maximum CFSE intensity zmax, which depends on the initial staining of cells. This upper level of fluorescence was prescribed specifically for data sets 1 and 2.
The population balance model (4), describing the distribution of cells n(t, z) structured according to the log10-transformed CFSE intensity, depends on the unknown rate functions of cell division α(z), death β(z) and the label loss ν(z). The identification of these functions from the observed CFSE histograms, using some measure of closeness of the model solution to the observations, represents an inverse problem. This problem is characterized by a finite set of observations ni, jand an infinite-dimensional space of the functions to be estimated. Following a general approach to the numerical solution of the parameter estimation problem for distributed parameter systems [28–33], we need to parameterize the elements of the function space in order to represent them by a finite set of parameters and to select the cost functional.
Here φ j are cubic polynomials, such that φ j (z j ) = 1, φ j (z k ) = 0 for j ≠ k, and hence α L (z k ) = a k , β L (z k ) = b k , k = 1, ..., L. Elements of the vectors and are the unknowns to be estimated.
In terms of the CFSE fluorescence level x, cf. model (2), the first case assumes that the rate of label decay is directly proportional to the amount of label expressed on the cell: v(x) = cx log 10, while the second one implies that the CFSE loss does not depend on its level on the cells: v(x) ≡ c, x ∈ [xmin, xmax].
Using the above parametrization, the original infinite dimensional problem of identifying the rate functions reduces to a finite dimensional one over a vector of parameters,p := [a, b, c, γ] ∈ ℝ2L+2.
The implementation details of the rate functions approximation are presented in the section ”Applications to CFSE assay” below.
To estimate the vector of best-fit parameters p*, we follow a maximum likelihood approach and seek for the parameter values which maximize the probability of observing the experimental data ni, jprovided that the true values are specified by the model solution n(t, z; p*). The choice of the probability function should take into account the statistical nature of the observation errors. Because the statistical characterization of the CFSE fluorescence histograms for growing populations of cells is a poorly analyzed issue, we follow the principle stated in : ”...in the absence of any other information the Central Limit Theorem tells us that the most reasonable choice for the distribution of a random variable is Gaussian.” Therefore, we assume that (i) the observational errors, i.e., the residuals defined as a difference between observed and model-predicted values, are normally distributed; (ii) the errors in observations at successive times are independent; (iii) the errors in cell counts for consecutive label bins are independent ((ii) – (iii) imply that the errors in the components of the state vector are independent); (iv) the variance of observation errors (σ2) is the same for all the state variables, observation times and label expression level.
provided that σ2 is assigned the value = Φ(p*)/n d , where p* is the vector which gives a minimum to Φ(p) and is the total number of scalar measurements. Relevant details of the computational treatment of the parameter estimation problem for the PDE model (4) are presented in the next section.
The parameter estimation problem for hyperbolic PDEs is non-trivial due to the hyperbolic nature of the equations (possible discontinuity of solutions) and due to the large size of the discretized problem. Moreover, model (4) is not a standard differential equation due to the solution term n(t, z + log10 γ) with the transformed argument z + log10 γ To our knowledge, no publicly available software package exists which deals with optimization (parameter estimation in particular) of models described by hyperbolic PDEs. For parabolic PDEs, which, after a suitable space discretization, can be treated as large systems of ODEs, available optimization tools (software, numerical methods) for large-scale problems can be used.
Solutions of a hyperbolic PDE can be discontinuous at the characteristic curve. Due to the solution term n(t, z + log10 γ) in model (4), the discontinuity of solutions at a point (t, z0) on the characteristic curve propagates to the points (t, z j ), z j = z0 - j log10 γ, j = 1, 2, .... A discretization of the initial-boundary value problem (4) should take into account the hyperbolicity of the equations and it should be robust and efficient since it is used in an optimization loop during model parameter identification. Moreover, available optimization tools for large-scale problems are based on some variants of Newton's method, which involves the computation of derivatives of the objective function with respect to the parameters to be estimated. These derivatives may not exist for discontinuous solutions. Note also that the optimization technique based on variants of Newton's method is efficient only if a good initial guess for the estimated parameters is available. For our problem, a derivative free minimization method which is robust with respect to the initial guess is preferable. Below we outline the numerical methods used and computational details of the problem under study.
The initial-boundary value problem
To solve the initial-boundary value problem (IBVP) for model (4), we use the Matlab program hpde by L. Shampine developed for systems of first order hyperbolic PDEs in one space variable . This program implements the well established second order Richtmyer's two-step variant of the Lax-Wendroff method (LxW) . This method is dispersive and therefore the software contains the possibility to apply after each time step a nonlinear filter  to reduce the total variation of the numerical solution. When the solution is smooth, filtering has little effect, but the filter is helpful in dealing with the oscillations which are characteristic of the LxW scheme when the solution is discontinuous or has large gradients. The choice of this method was also influenced by its ability to be fully vectorized, which allows to speed up computations in Matlab significantly. This is especially important when solving a PDE in an optimization loop. To compute the solution term with the transformed argument z + log10 γ, we modified the code hpde so that this term is interpolated, through its closest neighbors, preserving the second order accuracy of the LxW scheme.
The time step is recomputed at each iteration of the optimization procedure since it depends on the estimated function ν(z).
It is well known that solutions of a hyperbolic PDE are discontinuous if the compatibility condition for the initial and boundary conditions is not fulfilled. In our case the compatibility condition reads asn(0, zmax) = n0(zmax) = 0.
Our experience with the solution of the IBVP for model (4), using the code hpde, has shown that oscillations in the computed solution, occurring due to the discontinuity of the exact solution, do not propagate significantly with respect to z. Hence, the accuracy of the computed solution is only influenced locally, see Fig. 4. With the mesh refinement, the amplitude of the oscillations grows, while the interval of the propagation of the oscillations decreases, cf. Fig. 4 (right, top). The filtering procedure of the hpde smoothes the oscillations, see Fig. 4 (right, bottom).
If the exact solution of model (4) is smooth, the order of accuracy of the computed solution on the interval [zmin, zmax] is uniform and corresponds to the order of the LxW scheme. This is the case for data set 1, for which the initial function is compatible with the boundary condition, n0(zmax) = 0 for t0 = 72 hours. For N = 1000 the accuracy of the best-fit solution is about 10-3 - 10-2 and slowly decreases with time. For data set 2 the compatibility condition (13) is not fulfilled as n0(zmax) ≠ 0 for t0 = 0. In this case the solution is discontinuous at points = - j log10 γ, j = 0, 1, ..., 10, see Fig. 4, and the above level of accuracy can only be achieved outside some small intervals around the discontinuity points.
Since model (4) is linear with respect to n(t, z), we scaled it by the factor 10-5 to avoid the possible accuracy loss when dealing simultaneously with very large and small numbers in computations. To speed up the computations, the parameter estimation problem was treated in two stages. First we used a coarser mesh Z with N = 500 to solve the IBVP. Then the obtained best-fit parameter values were taken as a starting point to minimize the objective function using a finer mesh with N = 1000 to solve the IBVP.
Parameterization of the estimated functions
According to the proposed parameterization (7) of the functions α(z) and β(z), the parameters to be estimated are elements of the vectors and . Each pair (a k , b k ) approximate the corresponding rate function at some value z k ∈ [zmin, zmax] so that α L (z k ) = a k and β L (z k ) = b k , k = 1, ..., L. Values z k should be chosen such that all the consecutive divisions of cells could be captured properly. Hence the minimal value of L has to be larger than the maximal number of divisions cells have undergone. On the other hand, L should not be very large to treat the minimization problem efficiently. Values of α L (z) and β L (z) for z ≠ z k were evaluated with the code interp1 by ashape-preserving piecewise cubic interpolation. In the following we omit the subscript L for simplicity.
For the initial parameterization we used L = 8. After the best-fit solution was found, the parameterization of α(z) and β(z) was updated as follows. For α(z), we added new points, thus introducing additional parameters to be estimated. The increase of L was restricted by the requirement that adding new parameters should allow one a better fit of the data, i.e., lead to a significant improvement in the computed minimum of the objective function. For data set 1, all estimated b k were close to some constant value. Therefore, we assumed that β(z) can be treated as a constant function. This simplifying assumption leads to a minor change in the values of the objective function (1%). For data set 2, all b k corresponding to z k < 2.5 were zeros and we fixed them to be zero.
To solve the minimization problem, we use the Matlab code fminsearch implementing the Nelder-Mead simplex method. This method is a classical direct search algorithm that is widely used in case when the gradient of the objective function with respect to the estimated parameters cannot be evaluated. In our case the gradient, if it exists (i.e., if the solution of model (4) is continuous), can be computed numerically, but the computational cost is too large for the parameter estimation problem. As this method can trap in local minima for nonconvex objective functions, a number of runs with different initial guesses are necessary.
Applications to CFSE assay
In this section we investigate the appropriateness of the proposed label-structured PDE model (4), using the two original data sets introduced in section ”CFSE data”. The performance of this model with respect to the data sets is further compared with that of the compartmental ODE model developed recently in .
Mitogen-induced T cell proliferation
CD3/CD28 antibody induced T cell proliferation
Comparison to compartmental ODE model
The total number of live lymphocytes, N i , and the distribution of lymphocytes with respect to the number of divisions they have undergone, , at the indicated times t i
Time days t i
Total number of live cells N i
Numbers of cells w.r.t. the number of divisions (j) they undergone
data set 1
1.4 × 105
2.5 × 105
4.4 × 105
5.0 × 105
5.7 × 105
data set 2
The term 2αj - 1Nj - 1(t) for j ≥ 1 represents the cell birth (influx from the previous compartment because of division), whereas the term (α j + β j )N j (t) represents cell loss (outflux from the compartment) due to division and death. The model (15) allows an analytical solution which was used in the parameter estimation procedure.
which corresponds, under a set of assumptions similar to those presented above for the PDE model, to the maximum likelihood approach, see  for details. Here and N j (t i ; p) specify the data and the model solution for the cell generation j at time t i , p is the vector of estimated parameters α j and β j .
The best-fit parameter values α j and β j for data sets 1 and 2 are presented in Figs. 6 and 8. Note that these values are plotted at the middle of the intervals of the log10-transformed CFSE intensity which correspond to cells divided j times. Quantitatively, the rate functions α(z), β(z) and parameters α j , β j provide a different characterization of the cell kinetics. The qualitative bell-shaped behavior of the division rate with respect to the structure variable is similar for the PDE and ODE models (CFSE level and division number, respectively).
Many immunological phenomena result from cell proliferation. To quantify the cell proliferation, the technology based upon flow cytometry in conjunction with fluorescent dye (such as CFSE) that stain cell membrane or cytoplasm is extensively used in experimental and clinical research. It provides large amounts of data on the evolution of the histograms of fluorescence intensity of the cell population growing in response to a perturbing agent. The challenge is not only to collect the data, but also to analyze them in a way that enhances our understanding of the kinetics of the cellular responses. To this end, mathematical models are needed that quantitatively describe and interpret the data, in particular allowing one to estimate the rates of cell division and death.
Recently a number of different mathematical approaches have been proposed for the analysis of CFSE data [7–10, 12, 17, 25]. The corresponding models were instructive in appreciating the complexity of the parameter estimation problem from CFSE distribution data. These models take into account the heterogeneity of the growing cell populations with respect to the division number and are based upon systems of ordinary, delay or (age-structured) partial differential equations. None of these models considered the label intensity as a structure variable thus the CFSE histograms must be transformed into simplified descriptions of the generation structure of the population. This can be a vaguely defined procedure if the initial staining is not homogeneous. Such models, although easier to solve, cannot describe cell growth accurately enough due to the lack of structure information included.
In this study we developed a computational approach which allows a direct reference to the CFSE distributions. The label-structured cell population dynamics is described by a first order hyperbolic PDE model similar to those proposed by Bell and Anderson  for heterogenous cell populations structured by volume or size. The proposed model characterizes cell populations using the rate functions of cell division, death, label decay and the label dilution factor. We showed that this model provides a consistent mathematical tool for the analysis of CFSE-structured lymphocyte populations.
We presented a numerical approach for the parameter estimation implemented in the widely used package Matlab . The major elements of this approach are: (i) the smoothing of the histograms of CFSE data, which generates a continuous functional approximation of the distribution density of the cell population with a reduced level of noise; (ii) the software for the solution of the initial-boundary value problem for the proposed PDE model using the second-order Lax-Wendroff scheme; (iii) the parameterization of the rate functions in order to reduce the variational problem of CFSE data assimilation to a finite-dimensional parameter optimization task; (iv) the maximum likelihood approach to the parameter estimation.
Two original data sets from in vitro CFSE proliferation assays with human T lymphocytes were used to evaluate the performance of the proposed approach. It was shown that the model quantitatively describes the kinetics of the cell populations both at the global level and with respect to the fine structure of the CFSE distribution. The estimated rate functions provide a deeper insight into the turnover kinetics of the growing T cell populations. By computing the mean values of the rate functions for consecutive CFSE ranges corresponding to the sequential generations, one can characterize the division number dependent T cell turnover rates. In particular, the estimated rate functions α(z) imply the following division number dependent proliferation rates = [0.0023, 0.014, 0.020, 0.023, 0.010, 0.0032, 0.00002, 0] for data set 1 and = [0.0033, 0.0088, 0.033, 0.016, 0.0082, 0.0011] for data set 2, where j stands for the division number. The best-fit estimate of the reduction factor for CFSE per cell after division is smaller than two, ranging from 1.9 (data set 1) to 1.7 (data set 2). A possible interpretation might be that the CFSE molecules bonded to proteins upon release from cells dying in the process of division can be taken up actively or adhere to the live cells. Interestingly, the a quantitative CFSE data analyses published recently also indicates that the factor difference in median fluorescence intensity of adjacent CFSE peaks is typically not exactly 2  and there might be a few percent difference among siblings in the CFSE fluorescence inherited from the mother cell .
A number of issues require further systematic analysis: (i) the statistical error model underlying the fluctuations in the CFSE histograms; (ii) the level of noise smoothing used in the generation of the continuous distributions from the histogram data; (iii) the convergence of the finite-dimensional approximation of the rate functions estimation problem; (iv) the analysis of the confidence bounds for the estimated rate functions; (v) the application of Tikhonov regularization for the function identification (inverse) problem.
Overall, our study suggests that the label-structured modelling of cell population balance could become a component of the CFSE flow cytometry analysis software. The model's modifications can be used as building blocks for integrative mathematical description of complex in vivo labelling experiments in infected subjects such as those presented recently in [40, 41] for investigation of T cell activation and homeostasis.
The authors would like to acknowledge with thanks financial support of this work provided by the Russian Foundation of Basic Research (RFBR) and the Ministry of Flanders (MF) Programme ”Bilateral scientific cooperation” (RFBR-MF-05-01-02853-MF-a), the RFBR Grant RFBR-05-01-00732 and the U.S. Civilian Research and Development Foundation (Award RUX1-2710-MO-06).
- Doherty PC, Christensen JP: Accessing complexity: The dynamics of virus-specific T cell responses. Annu Rev Immunol. 2000, 18: 561-592. 10.1146/annurev.immunol.18.1.561.View ArticlePubMedGoogle Scholar
- Abu-Absi NR, Zamamiri A, Kacmar J, Balogh SJ, Srienc F: Automated flow cytometry for acquisition of time-dependent population data. Cytom Part A. 2003, 51: 87-96. 10.1002/cyto.a.10016.View ArticleGoogle Scholar
- Liu D, Yu J, Chen H, Reichman R, Wu H, Jin X: Statistical determination of threshold for cellular division in the CFSE-labeling assay. J Immunol Methods. 2006, 312: 126-136. 10.1016/j.jim.2006.03.010.View ArticlePubMedGoogle Scholar
- Srienc F: Cytometric data as the basis for rigorous models of cell population dynamics. J Biotechnol. 1999, 71: 233-238. 10.1016/S0168-1656(99)00026-7.View ArticleGoogle Scholar
- Lyons AB, Doherty KV: Flow cytometric analysis of cell division by dye dilution. Current Protocols in Cytometry. 2004, 9.11.1-9.11.10. Suppl 27Google Scholar
- Lyons AB: Divided we stand: tracking cell proliferation with carboxyfluorescein diacetate ester. Immunol Cell Biol. 1999, 77: 509-515. 10.1046/j.1440-1711.1999.00864.x.View ArticlePubMedGoogle Scholar
- Gett AV, Hodgkin PD: A cellular calculus for signal integration by T cells. Nat Immunol. 2000, 1: 239-244. 10.1038/79782.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: 455-476. 10.1016/j.jtbi.2004.04.011.View ArticlePubMedGoogle Scholar
- De Boer RJ, Perelson AS: Estimating division and death rates from CFSE data. J Comput Appl Math. 2005, 184: 140-164. 10.1016/j.cam.2004.08.020.View ArticleGoogle Scholar
- Ganusov VV, Pilyugin SS, De Boer RG, Murali-Krishna K, Ahmed R, Antia R: Quantifying cell turnover using CFSE data. J Immunol Methods. 2005, 298: 183-200. 10.1016/j.jim.2005.01.011.View ArticlePubMedGoogle Scholar
- De Boer RJ, Ganusov VV, Milutinovic D, Hodgkin PD, Perelson AS: Estimating lymphocyte division and death rates from CFSE data. Bul Math Biology. 2006, 68 (5): 1011-1031. 10.1007/s11538-006-9094-8.View ArticleGoogle Scholar
- Luzyanina T, Mrusek S, Edwards JT, Roose D, Ehl S, Bocharov G: Computational analysis of CFSE proliferation assay. J Math Biol. 2007, 54: 57-89. 10.1007/s00285-006-0046-6.View ArticlePubMedGoogle Scholar
- Mantzaris NV, Liou JJ, Daoutidis P, Srienc F: Numerical solution of a mass structured cell population balance model in an environment of changing substrate concentration. J Biotechnol. 1999, 71: 157-174. 10.1016/S0168-1656(99)00020-6.View ArticleGoogle Scholar
- Gurtin ME, MacCamy RC: Some simple models for nonlinear age-dependent population dynamics. Math Biosci. 1979, 43 (3/4): 199-211. 10.1016/0025-5564(79)90049-X.View ArticleGoogle Scholar
- Murray JD: Mathematical Biology. 1993, SpringerView ArticleGoogle Scholar
- Webb GF: Theory of Nonlinear Age-Dependent Population Dynamics. 1985, New York: Marcel DekkerGoogle Scholar
- Bernard S, Pujo-Menjouet L, Mackey MC: Analysis of cell kinetics using a cell division marker: mathematical modelling of experimental data. Biophys J. 2003, 84: 3414-3424.PubMed CentralView ArticlePubMedGoogle Scholar
- Sinko JW, Streifer W: A model for populations reproducing by fission. Ecology. 1971, 52: 330-335. 10.2307/1934592.View ArticleGoogle Scholar
- Bell GI, Anderson EC: Cell growth and division I. A mathematical model with application to cell volume distributions in mammalian suspension cultures. Biophys J. 1967, 7: 330-351.Google Scholar
- Diekmann O, Heijmans HJM, Thieme HR: On the stability of the cell size distribution. J Math Biol. 1984, 19: 227-248. 10.1007/BF00277748.View ArticleGoogle Scholar
- Metz JAJ, Diekmann O: The dynamics of physiologically structured populations. Lecture Notes in Biomathematics. 1986, Berlin: Springer, 68:Google Scholar
- Shampine L: Solving hyperbolic PDEs in MATLAB. Appl Numer Analysis and Comput Math. 2005, 2 (3): 346-358. 10.1002/anac.200510025.View ArticleGoogle Scholar
- Givan AL, Fisher JL, Waugh M, Ernstoff MS, Wallace PK: A flow cytometric method to estimate the precursor frequencies of cells proliferating in response to specific antigens. J Immunol Methods. 1999, 230: 99-112. 10.1016/S0022-1759(99)00136-2.View ArticlePubMedGoogle Scholar
- Matera G, Lupi M, Ubezio P: Heterogeneous cell response to topotecan in a CFSE-based proliferation test. Cytometry A. 2004, 62 (2): 118-128. 10.1002/cyto.a.20097.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 R Soc B. 2006, 273: 1165-1171. 10.1098/rspb.2005.3432.PubMed CentralView ArticlePubMedGoogle Scholar
- Volodine T, Vanderstraeten D, Roose D: Smoothing of meshes and point clouds using weighted geometry-aware bases. Lecture Notes in Computer Science. Edited by: Kim MS, Shimada K. 2006, Springer, 4077: 687-693.Google Scholar
- Hansen PC: Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. 1998, Philadelphia, PA, USAView ArticleGoogle Scholar
- Bertuzzi A, Gandolfi A, Vitelli R: Influence of cell loss in the analysis of proliferating populations by flow cytometry. Math Biosc. 1986, 82: 257-260. 10.1016/0025-5564(86)90005-2.View ArticleGoogle Scholar
- Banks HT, Fitzpatrick BG: Statistical methods for model comparizon in parameter estimation problems for distributed systems. J Math Biol. 1990, 28: 501-527. 10.1007/BF00164161.View ArticleGoogle Scholar
- Banks HT, Fitzpatrick BG: Estimation of growth rate distributions in size structured population models. Q Appl Math. 1991, XLIX (2): 215-235.Google Scholar
- Banks HT, Reich S, Rosen IG: An approximation theory for the identification of nonlinear distributed parameter systems. SIAM J Control Optim. 1990, 28: 552-569. 10.1137/0328033.View ArticleGoogle Scholar
- Banks HT: Some remarks on estimation techniques for size-structured population models. Lecture Notes in Biomathematics. 1994, Brlin: Springer, 100:Google Scholar
- Schittkowski K: Numerical Data Fitting in Dynamical Systems. 2002, Dordrecht: Kluwer Academic PublishersView ArticleGoogle Scholar
- Gershenfeld N: The Nature of Mathematical Modelling. 2002, Cambridge University PressGoogle Scholar
- Baker CTH, Bocharov GA, Paul CAH, Rihan FA: Computational modelling with functional differential equations: identification, selection and sensitivity. Appl Numer Math. 2005, 53: 107-129. 10.1016/j.apnum.2004.08.014.View ArticleGoogle Scholar
- Richtmyer RD, Morton KW: Difference Methods for Initial-Value Problems. 1967, New York: Wiley InterscienceGoogle Scholar
- Engquist B, Lötstedt P, Sjögreen B: Numerical filters for efficient shock computations. Math Comput. 1989, 52: 509-537. 10.2307/2008479.View ArticleGoogle Scholar
- Yates A, Chan C, Strid J, Moon S, Callard R, George AJ, Stark J: Reconstruction of cell population dynamics using CFSE. BMC Bioinformatics. 2007, 8: 196-10.1186/1471-2105-8-196.PubMed CentralView ArticlePubMedGoogle Scholar
- Grossman Z, Meier-Schellersheim M, Paul WE, Picker LJ: Pathogenesis of HIV infection; what the virus spares is as improtant as what it destroys. Nature Medicine. 2006, 12: 289-295. 10.1038/nm1380.View ArticlePubMedGoogle Scholar
- Hataye J, Moon JJ, Khoruts A, Reilly C, Jenkins MK: Naive and memory CD4+ T cell survival controlled by clonal abundance. Science. 2006, 312: 114-116. 10.1126/science.1124228.View ArticlePubMedGoogle Scholar
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.