Estimating Lymphocyte Division and Death Rates from CFSE Data

The division tracking dye, carboxyfluorescin diacetate succinimidyl ester (CFSE) is currently the most informative labeling technique for characterizing the division history of cells in the immune system. Gett and Hodgkin [Nat. Immunol. 1:239–244, 2000] have pioneered the quantitative analysis of CFSE data. We confirm and extend their data analysis approach using simple mathematical models. We employ the extended Gett and Hodgkin [Nat. Immunol. 1:239–244, 2000] method to estimate the time to first division, the fraction of cells recruited into division, the cell cycle time, and the average death rate from CFSE data on T cells stimulated under different concentrations of IL-2. The same data is also fitted with a simple mathematical model that we derived by reformulating the numerical model of Deenick et al. [J. Immunol. 170:4963–4972, 2003]. By a non-linear fitting procedure we estimate parameter values and confidence intervals to identify the parameters that are influenced by the IL-2 concentration. We obtain a significantly better fit to the data when we assume that the T cell death rate depends on the number of divisions cells have completed. We provide an outlook on future work that involves extending the Deenick et al. [J. Immunol. 170:4963–4972, 2003] model into the classical smith-martin model, and into a model with arbitrary probability distributions for death and division through subsequent divisions.     

Measuring lymphocyte proliferation, survival and differentiation using CFSE time-series data

Cellular proliferation is an essential feature of the adaptive immune response. The introduction of the division tracking dye carboxyfluorescein diacetate succinimidyl ester (CFSE) has made it possible to monitor the number of cell divisions during proliferation and to examine the relationship between proliferation and differentiation. Although qualitative examination of CFSE data may be useful, substantially more information about division and death rates can be extracted from quantitative CFSE time-series experiments. Quantitative methods can reveal in detail how lymphocyte proliferation and survival are regulated and altered by signals such as those received from co-stimulatory molecules, drugs and genetic polymorphisms. In this protocol, we present a detailed method for examining time-series data using graphical and computer-based procedures available to all experimenters.     

Computational analysis of CFSE proliferation assay

CFSE based tracking of the lymphocyte proliferation using flow cytometry is a powerful experimental technique in immunology allowing for the tracing of labelled cell populations over time in terms of the number of divisions cells undergone. Interpretation and understanding of such population data can be greatly improved through the use of mathematical modelling. We apply a heterogenous linear compartmental model, described by a system of ordinary differential equations similar to those proposed by Kendall. This model allows division number-dependent rates of cell proliferation and death and describes the rate of changes in the numbers of cells having undergone j divisions. The experimental data set that we specifically analyze specifies the following characteristics of the kinetics of PHA-induced human T lymphocyte proliferation assay in vitroL (1) the total number of live cells, (2) the total number of dead but not disintegrated cells and (3) the number of cells divided j times. Following the maximum likelihood approach for data fitting, we estimate the model parameters which, in particular, present the CTL birth- and death rate “functions”. It is the first study of CFSE labelling data which convincingly shows that the lymphocyte proliferation and death both in vitro and in vivo are division number dependent. For the first time, the confidence in the estimated parameter values is analyzed by comparing three major methods: the technique based on the variance–covariance matrix, the profile-likelihood-based approach and the bootstrap technique. We compare results and performance of these methods with respect to their robustness and computational cost. We show that for evaluating mathematical models of differing complexity the information-theoretic approach, based upon indicators measuring the information loss for a particular model (Kullback–Leibler information), provides a consistent basis. We specifically discuss methodological and computational difficulties in parameter identification with CFSE data, e.g. the loss of confidence in the parameter estimates starting around the sixth division. Overall, our study suggests that the heterogeneity inherent in cell kinetics should be explicitly incorporated into the structure of mathematical models.      

Carboxyfluorescein diacetate succinimidyl ester and the virgin lymphocyte: a marriage made in heaven

Carboxyfluorescein diacetate succinimidyl ester (CFSE) labelling of na?ve lymphocyte populations provides unique insights into the immune response. The clonal nature of immune responses, necessitating clonal expansion to achieve a sufficiently large number of Ag-reactive effector cells, combined with the dependence of lymphocyte differentiation on cell division, underlie the usefulness of CFSE in understanding the factors that regulate responses both in vitro and in vivo. We have combined CFSE labelling with Ag receptor transgenic models, using seven channel flow cytometry to track the correlation between cell division and a number of other parameters, such as surface expression of activation markers, cytokine receptors and homing receptors, cytokine production, cytotoxic activity and indicators of apoptosis. Our data have allowed us to classify and understand immune responses in novel ways, suggesting many further avenues of enquiry and indicating previously unrecognized relationships between cell division and eventual cell fate.     

Interpreting CFSE Obtained Division Histories of B Cells in Vitro with Smith–Martin and Cyton Type Models

The fluorescent dye carboxyfluorescin diacetate succinimidyl ester (CFSE) classifies proliferating cell populations into groups according to the number of divisions each cell has undergone (i.e., its division class). The pulse labeling of cells with radioactive thymidine provides a means to determine the distribution of times of entry into the first cell division. We derive in analytic form the number of cells in each division class as a function of time using the cyton approach that utilizes independent stochastic distributions for the time to divide and the time to die. We confirm that our analytic form for the number of cells in each division class is consistent with the numerical solution of a set of delay differential equations representing the generalized Smith–Martin model with cell death rates depending on the division class. Choosing the distribution of time to the first division to fit thymidine labeling data for B cells stimulated in vitro with lipopolysaccharide (LPS) and either with or without interleukin-4 (IL-4), we fit CFSE data to determine the dependence of B cell kinetic parameters on the presence of IL-4. We find when IL-4 is present, a greater proportion of cells are recruited into division with a longer average time to first division. The most profound effect of the presence of IL-4 was decreased death rates for smaller division classes, which supports a role of IL-4 in the protection of B cells from apoptosis.     

Modeling T Cell Proliferation and Death in Vitro Based on Labeling Data: Generalizations of the Smith–Martin Cell Cycle Model

The fluorescent dye carboxyfluorescein diacetate succinimidyl ester (CFSE) classifies proliferating cell populations into groups according to the number of divisions each cell has undergone (i.e., its division class). The pulse labeling of cells with radioactive thymidine provides a means to determine the distribution of times of entry into the first cell division. We derive in analytic form the number of cells in each division class as a function of time based on the distribution of times to the first division. Choosing the distribution of time to the first division to fit thymidine labeling data for T cells stimulated in vitro under different concentrations of IL-2, we fit CFSE data to determine the dependence of T cell kinetic parameters on the concentration of IL-2. As the concentration of IL-2 increases, the average cell cycle time is shortened, the death rate of cells is decreased, and a higher fraction of cells is recruited into division. We also find that if the average cell cycle time increases with division class then the qualify of our fit to the data improves.     

On analyzing circadian rhythms data using nonlinear mixed models with harmonic terms

Wang, Ke, and Brown (2003, Biometrics59, 804-812) developed a smoothing-based approach for modeling circadian rhythms with random effects. Their approach is flexible in that fixed and random covariates can affect both the amplitude and phase shift of a nonparametrically smoothed periodic function. In motivating their approach, Wang et al. stated that a simple sinusoidal function is too restrictive. In addition, they stated that "although adding harmonics can improve the fit, it is difficult to decide how many harmonics to include in the model, and the results are difficult to interpret." We disagree with the notion that harmonic models cannot be a useful tool in modeling longitudinal circadian rhythm data. In this note, we show how nonlinear mixed models with harmonic terms allow for a simple and flexible alternative to Wang et al.'s approach. We show how to choose the number of harmonics using penalized likelihood to flexibly model circadian rhythms and to estimate the effect of covariates on the rhythms. We fit harmonic models to the cortisol circadian rhythm data presented by Wang et al. to illustrate our approach. Furthermore, we evaluate the properties of our procedure with a small simulation study. The proposed parametric approach provides an alternative to Wang et al.'s semiparametric approach and has the added advantage of being easy to implement in most statistical software packages.     

Quantitative analysis of lymphocyte differentiation and proliferation in vitro using carboxyfluorescein diacetate succinimidyl ester

Mature T and B lymphocytes respond to receptor-delivered signals received during and following activation. These signals regulate the rates of cell death, growth, differentiation and migration that ultimately establish the behaviour patterns collectively referred to as immune regulation. We have been pursuing the philosophy that in vitro systems of lymphocyte stimulation, when analysed quantitatively, help reveal the logical attributes of lymphocyte behaviour. The development of carboxyfluorescein diacetate succinimidyl ester (CFSE) to track division has enabled the variable of division number to be incorporated into these quantitative analyses. Our studies with CFSE have established a fundamental link between differentiation and division number. Isotype switching, expression of T cell cytokines, surface receptor alterations and changes to intracellular signalling components all display independent patterns of change with division number. The stochastic aspects of these changes and the ability of external signals to independently regulate them argue for a probabilistic modelling framework for describing and understanding immune regulation.     

Evaluation of two models for human topoisomerase I interaction with dsDNA and camptothecin derivatives

Human topoisomerase I (Top1) relaxes supercoiled DNA during cell division. Camptothecin stabilizes Top1/dsDNA covalent complexes which ultimately results in cell death, and this makes Top1 an anti-cancer target. There are two current models for how camptothecin and derivatives bind to Top1/dsDNA covalent complexes (Staker, et al., 2002, Proc Natl Acad Sci USA 99: 15387-15392; and Laco, et al., 2004, Bioorg Med Chem 12: 5225-5235). The interaction energies between bound camptothecin, and derivatives, and Top1/dsDNA in the two models were calculated. The published structure-activity-relationships for camptothecin and derivatives correlated with the interaction energies for camptothecin and derivatives in the Laco et al. model, however, this was not the case for several camptothecin derivatives in the Stacker et al. model. By defining the binding orientation of camptothecin and derivatives in the Top1/dsDNA active-site these results allow for the rational design of potentially more efficacious camptothecin derivatives.     

Topology and robustness in the Drosophila segment polarity network

A complex hierarchy of genetic interactions converts a single-celled Drosophila melanogaster egg into a multicellular embryo with 14 segments. Previously, von Dassow et al. reported that a mathematical model of the genetic interactions that defined the polarity of segments (the segment polarity network) was robust (von Dassow et al. 2000). As quantitative information about the system was unavailable, parameters were sampled randomly. A surprisingly large fraction of these parameter sets allowed the model to maintain and elaborate on the segment polarity pattern. This robustness is due to the positive feedback of gene products on their own expression, which induces individual cells in a model segment to adopt different stable expression states (bistability) corresponding to different cell types in the segment polarity pattern. A positive feedback loop will only yield multiple stable states when the parameters that describe it satisfy a particular inequality. By testing which random parameter sets satisfy these inequalities, I show that bistability is necessary to form the segment polarity pattern and serves as a strong predictor of which parameter sets will succeed in forming the pattern. Although the original model was robust to parameter variation, it could not reproduce the observed effects of cell division on the pattern of gene expression. I present a modified version that incorporates recent experimental evidence and does successfully mimic the consequences of cell division. The behavior of this modified model can also be understood in terms of bistability in positive feedback of gene expression. I discuss how this topological property of networks provides robust pattern formation and how large changes in parameters can change the specific pattern produced by a network.     

Within-host population dynamics and the evolution of microparasites in a heterogeneous host population

Abstract Why do parasites harm their hosts? The general understanding is that if the transmission rate and virulence of a parasite are linked, then the parasite must harm its host to maximize its transmission. The exact nature of such trade‐offs remains largely unclear, but for vertebrate hosts it probably involves interactions between a microparasite and the host immune system. Previous results have suggested that in a homogeneous host population in the absence of super‐ or coinfection, within‐host dynamics lead to selection of the parasite with an intermediate growth rate that is just being controlled by the immune system before it kills the host (Antia et al. 1994). In this paper, we examine how this result changes when heterogeneity is introduced to the host population. We incorporate the simplest form of heterogeneity–random heterogeneity in the parameters describing the size of the initial parasite inoculum, the immune response of the host, and the lethal density at which the parasite kills the host. We find that the general conclusion of the previous model holds: parasites evolve some intermediate growth rate. However, in contrast with the generally accepted view, we find that virulence (measured by the case mortality or the rate of parasite‐induced host mortality) increases with heterogeneity. Finally, we link the within‐host and between‐host dynamics of parasites. We show how the parameters for epidemiological spread of the disease can be estimated from the within‐host dynamics, and in doing so examine the way in which trade‐offs between these epidemiological parameters arise as a consequence of the interaction of the parasite and the immune response of the host.     

Estimating the basic reproductive number from viral sequence data

Epidemiological processes leave a fingerprint in the pattern of genetic structure of virus populations. Here, we provide a new method to infer epidemiological parameters directly from viral sequence data. The method is based on phylogenetic analysis using a birth-death model (BDM) rather than the commonly used coalescent as the model for the epidemiological transmission of the pathogen. Using the BDM has the advantage that transmission and death rates are estimated independently and therefore enables for the first time the estimation of the basic reproductive number of the pathogen using only sequence data, without further assumptions like the average duration of infection. We apply the method to genetic data of the HIV-1 epidemic in Switzerland.     

Probability models and the applicability of statistical procedures in the identification of chromosomal fragile sites

Summary. Böhm et al. (1995, Human Genetics 95 , 249–256) introduced a statistical model (named FSM–fragile site model) specifically designed for the identification of fragile sites from chromosomal breakage data. In response to claims to the contrary (Hou et al., 1999, Human Genetics 104 , 350–355; Hou et al., 2001, Biometrics 57 , 435–440), we show how the FSM model is correctly modified for application under the assumption that the probability of random breakage is proportional to chromosomal band length and how the purportedly alternative procedures proposed by Hou, Chang, and Tai (1999, 2001) are variations of the correctly modified FSM algorithm. With the exception of the test statistic employed, the procedure described by Hou et al. (1999) is shown to be functionally identical to the correctly modified FSM and the application of an incorrectly modified FSM is shown to invalidate all of the comparisons of FSM to the alternatives proposed by Hou et al. (1999, 2001). Last, we discuss the statistical implications of the methodological variations proposed by Hou et al. (2001) and emphasize the logical and statistical necessity for fragile site identifications to be based on data from single individuals.     

Performance of a kinetic model for intracellular ice formation based on the extent of supercooling.

Cryomicroscopy was used to study the incidence of intracellular ice formation (IIF) in protoplasts isolated from rye (Secale cereale) leaves during subfreezing isothermal periods and in in vitro mature bovine oocytes during cooling at constant rates. IIF in protoplasts occurred at random times during isothermal periods, and the kinetics of IIF were faster as isothermal temperature decreased. Mean IIF times decreased from approximately 1700 s at -4.0 degrees C to less than 1 s at -18.5 degrees C. Total incidence of IIF after 200 s increased from 4% at -4.0 degrees C to near 100% at -15.5 degrees C. IIF behavior in protoplasts was qualitatively similar to that for Drosophila melanogaster embryos over the same temperature ranges (Myers et al., Cryobiology 26, 472-484, 1989), but the kinetics of IIF were about five times faster in protoplasts. IIF observations in linear cooling of bovine oocytes indicated a median IIF temperature of -11 degrees C at 16 degrees C/min and total incidences of 97%, 50%, and 19% at 16, 8, and 4 degrees C/min, respectively. A stochastic model of IIF was developed which preserved certain features of an earlier model (Pitt et al. Cryobiology 28, 72-86, 1991), namely Weibull behavior in IIF temperatures during rapid linear cooling, but with a departure from the concept of a supercooling tolerance. Instead, the new model uses the osmotic state of the cell, represented by the extent of supercooling, as the independent variable governing the kinetics of IIF. Two kinetic parameters are needed for the model: a scale factor tau 0 dictating the sensitivity to supercooling, and an exponent rho dictating the strength of time dependency. The model was fit to the data presented in this study as well as those from Myers et al. and Pitt et al. for D. melanogaster embryos with and without cryoprotectant, and from Toner et al. (Cryobiology 28, 55-71, 1991) for mouse oocytes. In protoplasts, D. melanogaster embryos, and mouse oocytes, the parameters were estimated from IIF times in the early stages of isothermal periods, while the osmotic state of the cell was relatively constant. In bovine oocytes, the parameters were estimated from linear cooling data. Without further calibration, the model was used to predict total IIF incidence under different cooling regimes. For protoplasts, D. melanogaster embryos, and bovine oocytes, the model's predictions were quite accurate compared to the actual data. In mouse oocytes, adjustment of the hydraulic permeability coefficient (Lp) at 0 degree C was required to yield realistic behavior.(ABSTRACT TRUNCATED AT 400 WORDS)     

Estimating average cellular turnover from 5-bromo-2'-deoxyuridine (BrdU) measurements

Cellular turnover rates in the immune system can be determined by labelling dividing cells with 5-bromo-2'-deoxyuridine (BrdU) or deuterated glucose ((2)H-glucose). To estimate the turnover rate from such measurements one has to fit a particular mathematical model to the data. The biological assumptions underlying various models developed for this purpose are controversial. Here, we fit a series of different models to BrdU data on CD4(+) T cells from SIV(-) and SIV(+) rhesus macaques. We first show that the parameter estimates obtained using these models depend strongly on the details of the model. To resolve this lack of generality we introduce a new parameter for each model, the 'average turnover rate', defined as the cellular death rate averaged over all subpopulations in the model. We show that very different models yield similar estimates of the average turnover rate, i.e. ca. 1% day(-1) in uninfected monkeys and ca. 2% day(-1) in SIV-infected monkeys. Thus, we show that one can use BrdU data from a possibly heterogeneous population of cells to estimate the average turnover rate of that population in a robust manner.     

A kinetic model of telomere shortening in infants and adults

We have previously demonstrated that telomeres shorten more rapidly in peripheral mononuclear cells (PBMC) of infants than in adults (Zeichner et al., Blood 93 (1999) 2824). Here we describe a mathematical model that allows quantification of telomere dynamics both in infants and in adults. In this model the dependence of the telomere dynamics on age is accounted by assuming proportionality between the body growth, as approximated by the Gompertz equation, and the increase in the number of PBMCs. The model also assumes the existence of two subpopulations of PBMC with significantly different rates of division. This assumption is based on the results from a previous analysis of in vitro data for telomere dynamics in presence of telomerase inhibitors and our recent data obtained by measurements of BrdU incorporation in T lymphocytes in humans (Kovacs et al., J. Exp. Med. 194 (2001) 1731). The average telomere length of PBMC was calculated as the average length of these two subpopulations. The model fitted our experimental data well and allowed to derive a characteristic time of conversion of the rapidly proliferating cells to slowly proliferating cells on the order of 20 days. The half-life of the slowly proliferating cells was estimated to be about 6 months, which is in good agreement with data obtained by independent methodologies. Comparison of the one-population and two-subpopulations models demonstrated that one population model cannot explain the observed parameters of the terminal restriction fragment (TRF) dynamics while two-subpopulations model does. These results suggest that the rapid telomere shortening in infants is largely determined by the faster PBMC turnover compared to adults. This may have major implications for elucidation of the HIV pathogenesis in infants. One can speculate that the more rapid course of the HIV disease in infants is due to the existence of rapidly dividing cells, which are susceptible to HIV infection. In addition, these results could have implications for understanding of mechanisms of aging.     

Calculation of Steel's cell loss factor and of the parameters of cellular kinetics for a population with an exponential growth of the cell number and cell death at G0 phase with a probability equal to 1

The method of calculation of three cell kinetics parameters (the Steel's cell loss factor phi, the proliferative pool Pc, and the mean number m of the proliferating cells after mitotic division of one cell) was shown to be the same for the exponential growth state of cell number with cell death at the G0-phase, and for the exponential growth state with cell death occurring immediately after mitosis. The value of the mean number delta of non-proliferating cells that appeared after mitotic division of one cell is different for these two models of the exponential growth state with the equal values of the other three parameters (phi, Pc, and m). A method is proposed for calculating the parameter delta on the data of the percentage of labeled cells obtained in the experiments with continuous cultivation of cells in the nutrient medium containing 3H-thymidine. The kinetics of cell line HL-60 (the experimental data of Foa et al., 1982) can be described at the first approximation, by a model of the exponential growth state with the cell death at the G0-phase, with Pc = 0.80, phi = 0.24, m = 1.61, delta = 0.39, and the life time of the non-proliferating cells tQ = 24 hours.     

Re-evaluation of phytohormone-independent division of tobacco protoplast-derived cells

We have used a [³H] thymidine incorporation assay and microscopic observation in order to reassess recently published data dealing with the response of tobacco protoplasts to phytohormones, lipochitooligosaccharides and peptides ( Harling et al . 1997 ; Hayashi et al . 1992 ; Miklashevichs et al . 1996 ; Miklashevichs et al . 1997 ; Röhrig et al . 1995 ; Röhrig et al . 1996 ; van de Sande et al . 1996 ; Walden et al . 1994 ). These proliferation assays reveal that, in contrast to published data, isolated cells of the investigated mutant plant lines axi159 ( Hayashi et al . 1992 ; Walden et al . 1994 ), axi4/1 ( Harling et al . 1997 ) and cyi1 ( Miklashevichs et al . 1997 ), which were generated by activation T-DNA tagging, were unable to grow in the absence of auxin or cytokinin. Furthermore, lipochitooligosaccharides which play a key role in the induction of nodules on roots of legumes were unable to promote auxin- or cytokinin-independent cell division in tobacco protoplasts as claimed by Röhrig et al . (1995 , 1996 ). The finding of van de Sande et al . (1996 ) that ENOD40 confers tolerance of high auxin concentration to wild-type tobacco protoplasts was also reinvestigated. The results of our investigations show that we were unable to reproduce the proliferation data presented in this study, which were obtained by counting tobacco protoplast-derived cells undergoing division. In total, none of the published data on phytohormone-independent division of tobacco cells could be reproduced.     

Inferring speciation and extinction rates under different sampling schemes

The birth-death process is widely used in phylogenetics to model speciation and extinction. Recent studies have shown that the inferred rates are sensitive to assumptions about the sampling probability of lineages. Here, we examine the effect of the method used to sample lineages. Whereas previous studies have assumed random sampling (RS), we consider two extreme cases of biased sampling: "diversified sampling" (DS), where tips are selected to maximize diversity and "cluster sampling (CS)," where sample diversity is minimized. DS appears to be standard practice, for example, in analyses of higher taxa, whereas CS may occur under special circumstances, for example, in studies of geographically defined floras or faunas. Using both simulations and analyses of empirical data, we show that inferred rates may be heavily biased if the sampling strategy is not modeled correctly. In particular, when a diversified sample is treated as if it were a random or complete sample, the extinction rate is severely underestimated, often close to 0. Such dramatic errors may lead to serious consequences, for example, if estimated rates are used in assessing the vulnerability of threatened species to extinction. Using Bayesian model testing across 18 empirical data sets, we show that DS is commonly a better fit to the data than complete, random, or cluster sampling (CS). Inappropriate modeling of the sampling method may at least partly explain anomalous results that have previously been attributed to variation over time in birth and death rates.     

Image analysis shows that variations in actin crossover spacings are random, not compensatory.

A recent paper by Bremer et al. (1991. J. Cell Biol. 115:689-703) has argued that the random angular disorder model for actin is wrong, and that the variations in crossover spacing observed in electron micrographs of F-actin filaments can be best explained by a compensatory disorder caused by the lateral slipping of the twin (or two-start) strands which comprise the actin filament. We have analyzed the images of F-actin presented in Bremer et al. and show that their data argues against compensatory disorder and in favor of random disorder, independent of the cause of the disorder. We also revise our estimate of the angular component and show that the magnitude of this disorder is about 5-6 degrees per subunit, which is less than the 10-12 degrees that we originally proposed.