Cell cycle phase transitions - Semantic Scholar
ARTICLE Identification of Age-Structured Models: Cell Cycle Phase Transitions E. Sherer,1 E. Tocce,1,2 R.E. Hannemann,1,3 A.E. Rundell,3 D. Ramkrishna1 1
School of Chemical Engineering, Forney Hall of Chemical Engineering, 480 Stadium Mall Way, Purdue University, West Lafayette, Indiana 47907; telephone: 765-494-4066; fax: 765-494-0805; e-mail: [email protected] 2 Chemical and Biological Engineering Department, University of Wisconsin, Madison, Wisconsin 3 Weldon School of Biomedical Engineering, Purdue University, West Lafayette, Indiana Received 7 May 2007; revision received 19 July 2007; accepted 9 August 2007 Published online 4 September 2007 in Wiley InterScience (www.interscience.wiley.com). DOI 10.1002/bit.21633
Introduction ABSTRACT: A methodology is developed that determines age-specific transition rates between cell cycle phases during balanced growth by utilizing age-structured population balance equations. Age-distributed models are the simplest way to account for varied behavior of individual cells. However, this simplicity is offset by difficulties in making observations of age distributions, so age-distributed models are difficult to fit to experimental data. Herein, the proposed methodology is implemented to identify an age-structured model for human leukemia cells (Jurkat) based only on measurements of the total number density after the addition of bromodeoxyuridine partitions the total cell population into two subpopulations. Each of the subpopulations will temporarily undergo a period of unbalanced growth, which provides sufficient information to extract age-dependent transition rates, while the total cell population remains in balanced growth. The stipulation of initial balanced growth permits the derivation of age densities based on only agedependent transition rates. In fitting the experimental data, a flexible transition rate representation, utilizing a series of cubic spline nodes, finds a bimodal G0/G1 transition age probability distribution best fits the experimental data. This resolution may be unnecessary as convex combinations of more restricted transition rates derived from normalized Gaussian, lognormal, or skewed lognormal transition-age probability distributions corroborate the spline predictions, but require fewer parameters. The fit of data with a single log normal distribution is somewhat inferior suggesting the bimodal result as more likely. Regardless of the choice of basis functions, this methodology can identify age distributions, age-specific transition rates, and transition-age distributions during balanced growth conditions. Biotechnol. Bioeng. 2008;99: 960–974. ß 2007 Wiley Periodicals, Inc. KEYWORDS: population balance; inverse model; cell division; BrdU; cell generations
The growth behavior of an individual cell depends on a complex system of molecular interactions, feedback regulations, and synthesis reactions (Darzynkiewicz et al., 1996; Ja´nossy et al., 2004; Nova´k and Tyson, 1993, 2004). This complexity can lead to processes that are distributed throughout the life cycle (Liou et al., 1998), such as protein (Kromenaker and Srienc, 1991) or RNA (Arino and Kimmel, 1987) production and susceptibility to certain chemotherapeutic agents (Del Bino et al., 1991). Researchers are beginning to exploit the non-uniformity of these processes with the aid of predictive models (Mantzaris et al., 2002; Sherer et al., 2006). The heterogeneity of the resulting cell population is often represented in models by a single or small set of characteristic variables. Common characteristic variables include: cell mass (Mantzaris et al., 1999; Subramanian et al., 1970; Subramanian and Ramkrishna, 1971), volume (Webb and Grabosch, 1987), DNA content, protein content (Kromenaker and Srienc, 1991), or age (Sennerstam and Stro¨mberg, 1995; Tyrcha, 2001; Zamamiri et al., 2002). A population balance model describes the cellular dynamics based on growth and transition rates that are functions of the characteristic variables (Fredrickson and Tsuchiya, 1963; Fredrickson et al., 1967; Fredrickson and Mantzaris, 2002). Age, the amount of time since a cell last made a transition from one phase of the cell cycle to the next, is a simple way to account for the varied behavior of individual cells. That is, the amount of time spent by cells in any one phase of the cell cycle is distributed, so the likelihood of a cell transitioning from one phase to the next is dependent upon its age. An
Correspondence to: D. Ramkrishna
960
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008
ß 2007 Wiley Periodicals, Inc.
age-dependent transition rate can be used to describe this variation. For example, in the cell cycle, certain intracellular tasks must be completed before a cell will proceed to the next phase, so it is unlikely that a very young cell will transition. Also, checkpoints or other disturbances can shift and bias the transition age. Lognormal distributions are most commonly used to define transition-age distributions of cell cycle phases (Steel, 1977). The bias is captured even more closely if linear combinations of distributions (Shackney and Shankey, 1999) or asymmetry factors (Baisch et al., 1995) are incorporated. Reaction probability density functions provide a more mechanistic view of cellular transitions (Hjortsø, 2006) so long as the rate of each step in the reaction sequence is known. The advantages of utilizing age as the characteristic variable to describe the transition rates are diminished by the inability to easily measure the age distributions within a population of cells. While inverse methods are available to extract individual behavior information from population observations, these require quantitative measurements of the characteristic variables (Ramkrishna, 2000; Tyrcha, 2004). Synchronizing cells would allow for age approximations, but the synchronization process can lead to altered growth patterns (Chiorino et al., 2001; Cooper, 2003, 2004; Liou et al., 1998; Liu, 2005). Continuous observation of cell cultures is currently the method of choice for extracting individual age-dependent behavior (Kutalik et al., 2005), but this method is tedious as a significant number of cells must be observed to produce relevant distributions. In this paper, we develop an integrated experimental and analytical technique to overcome the difficulty of extracting the age distributions and the associated age-dependent transition rates from experimental data. Previous derivations have estimated age distributions assuming balanced growth conditions and a fixed, constant valued transition rate (Baisch and Otto, 1993; Bernard et al., 2003; Shackney, 1973; Steel, 1977). For a clear understanding of how models are identified using our methodology, it is necessary to recognize that the approach to balanced exponential growth of the overall population is preceded by two stages: (i) physiological changes in individual cells and (ii) a transient period of population adjustment to the time-invariant distribution characteristic of the balanced growth phase. The second stage represents a situation when stage (i) has been completed, but the population distributions in the different cell cycle phases are evolving to the constant proportion in balanced growth. Stage (i) can be achieved by maintaining a constant culture environment and is characterized by timeindependent transition rates. Balanced growth is achieved upon completion of stage (ii) so the age densities are also time independent (Dyson et al., 2000; Fredrickson et al., 1967). This feature of balanced growth implies that the system has reached a state such that the population densities of the phases will maintain constant ratio to one another (Dyson et al., 2002). Herein, we show that assuming balanced growth enables the calculation of age distributions
based only on knowledge of the transition rates. Thus, the only requirement for attaining the age distributions associated with the cell cycle dynamics are the transition rates. A problem arises while trying to determine the vectoral unknowns describing the age distributions associated with each phase of the cell cycle from scalar population density data. We demonstrate how an experimental strategy that labels a subpopulation can overcome this difficulty. In utilizing this method, cells are subjected to pulse labeling with bromodeoxyuridine (BrdU). Only cells in the S-phase of the cell cycle at the time of BrdU application will be labeled with BrdU (Dien and Srienc, 1991; Dolbeare et al., 1983; Sena et al., 1999). This creates two independent subpopulations, one labeled and the other unlabeled, both of which continue around the cell cycle. Two important features of this technique are: (i) that each subpopulation begins in unbalanced growth as their age distributions are changing with time and (ii) the BrdU does not affect the growth behavior of individual cells, so the total population remains in balanced growth (Baisch and Otto, 1993; Ormerod and Kubbies, 1992). This transient behavior is a feature of how the total population densities in each of the cell cycle phases are evolving to the eventual constant proportion in fully balanced growth. Thus, the measurement of labeled and unlabelled population densities provides the vectoral transient data that has the potential to identify vectoral transition rates without explicitly measuring age distributions. While BrdU is used herein to determine the age-structured cell cycle phase transition rates, similar experimental approaches can be used to determine cell division rates. The labeled mitosis method using 3H-thymidine (Steel, 1977) or the transient cell generation distributions from pkh dyes or carboxyfluorescein diacetate succinimidly ester (Lyons, 1999; Nordon et al., 1999) also produce transient data that can be used to determine the age-structured cell division rate. In presenting the methodology for such identification, the model and mathematical foundation are first developed and followed by an in vitro illustration for its use. The Analytical Foundation Section presents the cell growth model describing BrdU labeling with a proof of extractions of transition rates from subpopulations in unbalanced growth while the total populated maintains balanced growth. The Experimental Materials and Methods Section outlines the experimental materials and methods used to collect in vitro data. The age-dependent transition rates that result from the experimental BrdU results are discussed in the Results and Discussion Section, and the fitness is compared with normalized Gaussian, lognormal, and skew lognormal age density derived transition rates. The Conclusions Section summarizes the major conclusions and findings.
Analytical Foundation This section introduces an age-structured cell cycle model (Model Structure and Derivation of Initial Conditions
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
961
Sections), proves that the model predictions are sufficient to identify age-structured transition rates using BrdU labeling (Extraction of Age-Dependent Parameters From Partitioned Population Data Section), and shows the simulation methodology for comparison with the BrdU pulse-labeling experiments (Implementation Section).
The age density is subject to the initial condition nðt; 0Þ ¼ n0 ðtÞ and a boundary condition that consists of cells that have just transitioned out of a phase to become new cells in the next phase. Z 1 nð0; tÞ ¼ Din ðtÞnðt; tÞ dt 0
Model Structure
where
The division cycle of a cell is discretized into three phases in our model: (1) a combined G0 and G1 phase, (2) an S phase, and (3) a combined G2 and M phase (see Fig. 1). Newly divided cells may either be in the quiescent G0 phase or preparing to undergo DNA synthesis in the G1 phase. These two states were indistinguishable experimentally and are combined into a single phase. The S phase is a period when the DNA and many proteins are replicated before cell division. Cells then enter the second gap phase, G2, before physically dividing by undergoing mitosis, M. As the G2 and M states were also indistinguishable experimentally, they are lumped into a single phase. The cells were maintained in a nutrient rich environment throughout the in vitro experiments so a low level of quiescence is anticipated. The number density and transition rate of cells in the ith phase can be written as ni(t, t) and Gi(t) respectively where t is the age, or time elapsed since entering the current phase, and t is the time. It is assumed that the transition rate is stable, or time independent, in balanced growth conditions but the number density will increase with time as cells divide. Neglecting cell death, the population balance equation can be written as @nðt; tÞ @nðt; tÞ þ ¼ Dout ðtÞnðt; tÞ @t @t
(1)
2
3 n1 ðt; tÞ nðt; tÞ ¼ 4 n2 ðt; tÞ 5 n3 ðt; tÞ and G1 ðtÞ Dout ðtÞ ¼ 4 0 0
0 0 Din ðtÞ ¼ 4 G1 ðtÞ 0 0 G2 ðtÞ
3 0 0 G2 ðtÞ 0 5 0 G3 ðtÞ
3 2G3 ðtÞ 0 5 0
Under balanced growth conditions, the age distribution is time invariant so that the ith phase number density can be separated into the product of the population density, Ni(t), and the age-density, fi(t) Z 1 fi ðtÞ ¼ 1 (2) ni ðt; tÞ ¼ Ni ðtÞfi ðtÞ; where 0
Substitution of Equation (2) into the population balance equation (Eq. (1)) and integration over all ages gives an ordinary differential equation for the population increase in each of the phases (see Appendix) dNðtÞ ¼ DNðtÞ dt
(3)
where Dij ¼
where
2
2
Z 0
1
h
i out Din ðtÞ D ðtÞ fj ðtÞ dt ij ij
Only two quantities are required for simulating this system: transition rates and initial number densities. As shown in the Derivation of Initial Conditions Section, the balanced growth age distributions and associated number densities can be calculated using only the transition rates under balanced growth conditions. This implies that knowledge of the transition rates is sufficient for complete model characterization.
Derivation of Initial Conditions Under balanced growth, the matrix of rate constants, D, is time invariant and its positive real eigenvalue, m, is the specific growth rate for both the total population and each cell cycle phase. The system of ordinary differential equations can be written in terms of m and the corresponding eigenvector, h
Figure 1. Phases of cell cycle model. Cells transition out of the i th phase according to an age-dependent transition rate (Gi) with cell division occurring via the G2/M to G0/G1 transition.
962
mh ¼ Dh P where the normalized eigenvector hi = j hj is the fraction of cells in the ith cell cycle phase. To find the age-distribution within each phase, the partial differential equation system is first reduced to ordinary by
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
substitution of the separated number density (Eq. (2)) into the population balance equation (Eq. (1)). As the total number of cells in each phase increases according to the same overall specific growth rate in balanced growth, ordinary differential equations for the change in the age distribution with age are obtained (see Appendix)
" # 3 X dfi ðtÞ out ¼ ðDij ðtÞfj ðtÞÞ þ mfi ðtÞ dt j¼1
(4)
Using the normalization condition stated in Equation (2), the age distribution in each phase can be written as Z
Z exp
1 00 00 ðDout ðt Þ þ mÞ dt dt0 ii 0 0 (5)
Z t out 0 0 exp ðDii ðt Þ þ mÞ dt
fi ðtÞ ¼
1
0
t0
So, during balanced growth, the age-distributed number density is given as X h ni ðt; tÞ ¼ P i fi ðtÞ Nj ðtÞ j hj j
(6)
This expression can be used as an initial condition before a perturbation is introduced through the addition of BrdU. While the Model Structure Section describes an agestructured cell cycle model subject to the initial agedistributions derived in the Derivation of Initial Conditions, the Extraction of Age-Dependent Parameters From Partitioned Population Data Section addresses the question as to whether BrdU labeling experimental data is sufficient to identify the transition rates or age-distributions. Extraction of Age-Dependent Parameters From Partitioned Population Data In a balanced growth condition there is no unique set of transition rates and initial conditions that describes the total
Figure 2. Decomposing the Jurkat cells into cell cycle phases with BrdU labeling the S phase cells 1 h after labeling. ModFit software is used to determine the percentage of cells in each phase. a: A bivariate BrdU versus DNA content density plot is used to distinguish cells that incorporated BrdU (cells above the drawn line) from those that did not (cells below the line). b: A DNA histogram of the total population is used to determine the percentage of cells in each phase for the overall cell population. c: DNA histogram of labeled (BrdUþ) cells. d: DNA histogram of unlabeled (BrdU) cells.
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
963
population data. Unbalanced growth data with a known initial age distribution and minimal fluctuations in transition rates are needed for their identification, and BrdU pulse-labeling of a balanced growth population satisfies all of these requirements. BrdU labeling (Fig. 2a) partitions the balanced growth population (Fig. 2b) into two subpopulations (Fig. 2c and d). Each of these subpopulations will undergo unbalanced growth while the total population continues in balanced growth with unaltered transition rates. BrdU is an analog to the DNA base thymidine, a uridine derivative, such that cells that are actively synthesizing DNA, the S phase cells, will incorporate greater amounts of BrdU than cells in the G0/G1 or G2/M phases. By pulsing, or removing excess BrdU after a short time period, only those cells that were in the S phase at the time of application, t ¼ 0, will incorporate the BrdU. (We found that direct removal of all excess BrdU from the media proved difficult experimentally. However, labeling an aliquot of the original population, removing excess BrdU, and then recombining with the remainder of the original cell population does result in pulse labeling. The BrdU unexposed cells likely absorb residual BrdU but in small enough amounts to appear as BrdU unlabeled.) This labeling process may be represented by a pre-multiplying projection matrix, L, which concurrently produces a second projection matrix, I L, so that there are two subpopulation distributions among the different phases described by Ln0(t) and (I L)n0(t). The BrdU technique produces the following 2
0
3
6 7 Ln0 ðtÞ ¼ 4 an0;2 ðtÞ 5; 0
2
6 ðI LÞn0 ðtÞ ¼ 4
an0;1 ðtÞ þ ð1 aÞn0;1 ðtÞ
3 (7)
7 ð1 aÞn0;2 ðtÞ 5 an0;3 ðtÞ þ ð1 aÞn0;3 ðtÞ
where the variations in age-density with time lead to a transition matrix, D, that now also becomes a function of time. The population at the time of labeling, N0, is in balanced growth, so it is contained purely in the eigenspace, Em, which corresponds to the specific growth rate eigenvalue, m. Another way of stating this is to say that the projection, Pm, of all three-dimensional vectors onto Em maps N0 into N0 (i.e., PmN0 ¼ N0). This property clearly cannot be shared by LN0 and (I L)N0 simultaneously as PmL 6¼ L. Thus, we conclude that the transients produced by either of the initial conditions in Equation (8) will be truly three-dimensional and, in principle, sufficient for the identification of the three transition rates contained in the matrix Dout(t). The pulselabeling of cells with BrdU is required for this identification to take place.
Implementation After calculation of the initial number densities, transient number densities are computed using a variation of the method of characteristics (Ramkrishna, 2000): the method of successive generations (Liou et al., 1997). The number density for the ith phase, ni(t, t), is discretized into a series of ðkÞ k number densities, ni , with k representing the number of transitions undergone or generations. The BrdU labeling further separates each initial number density into labeled and unlabeled portions so the initial generation (k ¼ 0) number densities are given by Equation (7). For example, all cells beginning the S phase (i ¼ 2, k ¼ 0) are labeled and will move to the G2/M (i ¼ 3, k ¼ 1) then G0/G1 (i ¼ 1, k ¼ 2) phases before returning to the S phase (i ¼ 2, k ¼ 3). Following this pattern, the total phase number density of labeled and unlabeled cells will be Labeled ðBrdU þ Þ Subpopulation
where a is the fraction of the total population exposed to BrdU. (In the BrdU experiments the aliquot of BrdU exposed cells is added to a population not exposed to BrdU—we found this helped achieve true pulse-labeling. The unlabeled boundary condition in Equation (7) reflects this mixture.) It is now of interest to follow the population distributions obtained by solving Equation (1) for the two initial conditions stated above. However, since only the total, rather than age-distributed, population density is measurable, the equations of interest are only those in the total population density (Eq. (3)), for the two initial conditions as shown below
n1 ðt; tÞ ¼ a n2 ðt; tÞ ¼ a n3 ðt; tÞ ¼ a
k¼0 1 P
ð3kþ2Þ
n1
ðt; tÞ
ð3kÞ
n2 ðt; tÞ ð3kþ1Þ
n3
ðt; tÞ
Unlabeled ðBrdU Þ Subpopulation n1 ðt; tÞ ¼ n2 ðt; tÞ ¼
1 P k¼0 1 P k¼0 1 P k¼0
964
k¼0 1 P
k¼0
n3 ðt; tÞ ¼ dNðtÞ ¼ DðtÞNðtÞ; dt LN0 Nð0Þ ¼ ðI LÞN0
1 P
ð3kÞ
ð3kþ1Þ
ðn1 ðt; tÞþn1 ð3kÞ
ð3kþ1Þ
ðð1 aÞn2 ðt; tÞþn2 ð3kÞ
ð3kþ2Þ
ðt; tÞÞ
ð3kþ2Þ
ðt; tÞÞ
ðt; tÞþð1 aÞn1 ðt; tÞþn2
ð3kþ1Þ
ðn3 ðt; tÞþð1 aÞn3
ð3kþ2Þ
ðt; tÞþn3
ðt; tÞÞ (9)
labeled population unlabeled population
(8) The number densities for each of the generations can be found by backtracking along the characteristics until either
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
the initial condition, for the k ¼ 0 generation, or boundary condition, for the k > 0 generations, is reached Initial Conditions : ( Rt ð0Þ ð0Þ n ðt; tÞ ¼ ni ðt t; 0Þ exp tt Gi ðt 0 Þ dt 0 t t iðkÞ ni ðt; tÞ ¼ 0; k ¼ 1; 2; . . . Boundary Conditions : ( ð0Þ ni ðt; tÞ ¼ 0 t>t Rt ðkÞ ðkÞ ni ðt; tÞ ¼ ni ð0; t tÞexp 0 Gi ðt0 Þ dt 0 ; k ¼ 1; 2; . . . (10)
The initial cell distribution is known, and the boundary conditions are calculated successively as the boundary cell density depends on the final solution in the previous phase nðkÞ ð0; tÞ ¼
Z
1
Din ðt; tÞnðk1Þ ðt; tÞ dt
(11)
0
The final number densities are calculated and the appropriate combinations of these final number densities (Eq. (9)) describe the dynamics of the subpopulations.
Experimental Materials and Methods Materials and Methods Jurkat cells were grown in 75 cm2 flasks (Corning 430725) at 378C and 5% CO2 in 20 mL of a 90/10 mixture of Iscove’s Modified Dulbecco’s Medium (Sigma I3390) and fetal bovine serum (Sigma F2442) supplemented with 1 mL Lglutamine (Sigma G2150)/10 mL medium. Exponential growth was maintained for a week prior to the addition of BrdU by passing the cells every 2 days so concentrations remained between 105 and 106 cells/mL. One day before the addition of BrdU, cells were passed to create two flasks at a concentration of approximately 2 105 cells/mL in 20 mL of media. Also 1 day before pulse labeling, 20 mL of complete media was placed in the incubator to pre-warm and pre-gas in preparation for re-suspension after pulse labeling. The BD Biosciences BrdU Flow Kit (559619) was used for the BrdU procedures. A 10 mg/mL stock BrdU in 1 DPBS solution was diluted to 1 mM by adding 31 mL of the stock to 1 mL of 1 DPBS. For pulse-labeling, a 4 mL aliquot is exposed to 40 mL of the diluted BrdU solution (10 mM final BrdU concentration) in a 50 mL centrifuge tube. Forty-five minutes after inoculation, the remaining 16 mL of cells were transferred to the centrifuge tube (so a, the fraction of unlabeled cells, ¼0.2) which was then spun at 1,000 rpm for 5 min before resuspending the cells in 20 mL of pre-warmed and pre-gassed complete media and transferring to the original flask. Samples of 1 mL were collected every 2–4 h for approximately 40 h beginning with the newly resuspended cells. Samples were washed, after centrifuging for 5 min at 1,000 rpm and discarding the supernatant, in 1 mL of
Pharmingen Stain Buffer (PSB); resuspended in 100 mL of BD Cytofix/Cytoperm buffer for 15 min at room temperature; washed in a 1 mL of 1 Perm/Wash solution; and given a final wash in 1 mL of PSB. The cells were preserved by re-suspending in 1 mL of PBS and storing in a 48C freezer. Cytometry analysis was performed within a week of sample collection. In preparation for cytometry analysis, the samples were centrifuged, resuspended in 100 mL Cytofix/Cytoperm buffer for 5 min, washed with 1 mL of 1 Perm/Wash Buffer, and resuspended in 100 mL of a 30% 1 mg/mL DNase and 70% 1 DPBS solution all at room temperature. The cells are incubated for one hour at 378C. The cells were labeled for BrdU uptake after they are washed in 1 mL of Perm/Wash Buffer and resuspended in a solution containing 1 mL of FITC and 50 mL of 1 Perm/Wash buffer for 20 min. The DNA was then labeled by washing cells in 1 mL of 1
Perm/Wash buffer, and resuspending in 1 mL of propidium iodide (PI) solution (30 mL of 0.1% (v/v) Triton X-100 (Sigma 93443) in PBS, 6 mg DNAse-free Rnase A (Sigma R5503) and 600 mL of a 1 mg/mL PI stock solution (Sigma P4864)). Samples were immediately taken to a BD-Elite flow cytometer where 10,000 cells were analyzed per sample. An argon laser was used for excitation at 488 nm with the FITC emission collected at 525 nm and the PI emission at 675 nm. There was very little emission overlap so no compensation was required. The labeled and unlabeled cells were gated using FITC emission and DNA histograms of both populations were measured with the PI. This pulse labeling BrdU experiment was performed in duplicate and the two cultures, derived the same mother culture, gave extremely consistent results at all time points.
Sample Analysis The percentage of cells in each phase for labeled and unlabeled cells was extracted from the cytometry data by differentiating phases based on the DNA content. The bivariate BrdU incorporation against DNA content density plot were then used to gate the BrdU labeled and unlabeled cells (Fig. 2a). The phase percentages were determined by plotting the DNA histogram of each of these subpopulations (see Fig. 2c and d) and using the ModFit DNA labeling cell cycle analysis software package to define the phase boundaries.
Results and Discussion Extraction of Phase Transition Rates From Jurkat Cells Using BrdU Methods Each transition rate function was fit to experimental BrdU total number density data in each labeled and unlabeled cell cycle phase using the Matlab multivariate mini-
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
965
mization function ‘‘fminsearch.’’ This function employs a Nelder–Mead Simplex method to minimize the chosen least-squares objective function. The least-squares error, E, is the summation over all data points collected—Nsamples collected (13 samples) and Nphases phases in the independent subpopulations (three unlabeled þ three labeled phases)— of the squared difference between the experimental cell cycle phase percentages, xexp i;j , and the model predicted percentages, xsim , that result from the proposed transition rates i;j E¼
NX phases NX samples i¼1
exp sim 2 ðxi;j xi;j Þ
(12)
j¼1
For each independent subpopulation, the experimental phase percentage in the ith phase and jth time sample (where the jth sample is collected at time tj) is the ratio of the number of cells in the phaseP to the total number of cells in the exp exp subpopulation: xexp ¼ N = i;j i;j k Nk;j . Similarly, the simulation phase percentages are X sim xsim N sim ðtj Þ i;j ¼ Ni ðtj Þ= k k where Nisim ðtj Þ
¼
Z
1
ni ðt; tj Þ dt:
0
In most cases, hypotheses and models are written regarding the transition age probability distribution, hi(t), rather than the transition rate, Gi(t), itself. The likelihood of changing cell cycle phases within a certain age window is represented by that phase’s transition-age probability distribution. While this distribution may be observable experimentally, a kinetic model, such as the population balance models, use rate terms such as the rate of transition for a cell of a given age. The age-dependent transition rates can be calculated from the transition-age distributions (Eakman et al., 1966) Gi ðtÞ ¼ such that
h ðtÞ Rti 1 0 hi ðt0Þ dt0 2
hi ðtÞ ¼ Gi ðtÞ41
Zt
The flow cytometry data was analyzed as described in the Sample Analysis Section. The results for pulse-labeling show that the BrdU labeled cells are initially synchronized in the S phase (Fig. 3a). As time progresses, this plug moves to the labeled G2/M phase and then divides into the labeled G0/G1 phase (Fig. 3b) before reentering the labeled S phase (Fig. 3c). The unlabeled cells are also initially synchronized due to labeling of the S phase cells, but to a lesser degree than the unlabeled population due to the addition of the cells not exposed to BrdU. The rate of loss of the synchrony in the labeled and unlabeled subpopulations is dependent upon the shape of the cell cycle phase transition rates.
ð13Þ Determination of Unrestricted Age-Dependent Transition Rates Using BrdU Pulse-Labeling
3 hi ðt 0 Þ dt 0 5;
hi ð0Þ ¼ 0
0
In this study, we evaluated and compared the transition rates derived from normalized Gaussian, lognormal, and skewed lognormal transition rates, along with a more flexible transition rate representation utilizing series of cubic splines with movable spline nodes locations. Ages up to 50 h were evaluated in 0.25 h increments. Polynomials were used to approximate all transition rates with 71 orders giving an excellent fit over the entire age ranges for the transition rates investigated.
966
Figure 3. A series of bivariate BrdU versus DNA cell density plots illustrate the progression of the labeled and unlabeled subpopulations progression around the cell cycle. Ten thousand cells were collected for each figure with speckled or empty areas representing lower densities. The line drawn through each of the plots separates the BrdU labeled population from the unlabeled population. a: Immediately after BrdU inoculation (t ¼ 1 h). b: Nine hours after BrdU inoculation. c: Twenty-four hours after BrdU inoculation.
To determine the general shape of the cell cycle phase transition rates, a flexible representation was first implemented using a series of 15 cubic splines for each of the three cell cycle phase transition rates. The transition rates were manipulated to best fit the average of the two BrdU pulselabel experiments. (Again, these two experiments were conducted simultaneously and cells from each experiment share the same mother culture.) The spline representation was flexible in that both the value of the spline (transition rate) as well as its location (age) could be manipulated to reduce errors due to improper spline node placement. The parameters associated with the age-dependent transition rates were identified by fitting the cell cycle data; the fitted
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
minimal. While the age of cells transitioning out of the S phase or G2/M phase is relatively consistent, the transition age probability distribution for the G0/G1 phase is actually bimodal. To assess the significance the bimodal phenomenon and ensure it is not an artifact of (i) over-fitting with the splines, (ii) the pulse data itself, or (iii) the balanced growth restrictions, the effect of each of these factors was evaluated. The importance of bimodal distributions is compared with single mode distributions in the Determination of Alternative Functional Forms for the Age-Dependent Transition Rates Section and the predictions based on pulse-data experiments are compared with continuous labeling data in the Comparisons With Continuous Labeling Data Section. The implication of the assumption of balanced growth was examined by including the initial phase percentages of each subpopulation as fitted parameters. (The initial phase percentages need not sum to 1 in this case.) As seen in Figures 4 and 5a and b, neither the cell cycle phase transition age probability distribution nor the fits to the experimental data were altered significantly. Determination of Alternative Functional Forms for the Age-Dependent Transition Rates
Figure 4. Best fit of cell cycle phase transition rates to BrdU pulse-labeling data using a 15 node cubic spline representation for the transition rate out of each phase. Both the rate and age locations of each node are variable in the minimization. a: Cell cycle phase transition rates—spline nodes and the polynomial approximation of the cubic spline fit. b: Cell cycle phase transition age densities resulting from those transition rates. cell cycle phase age-dependent transition rate are shown in Figure 4a and the resulting transition age distributions in Figure 4b. The fit of the thick lines shown in Figure 5a and b to the experimental pulse label BrdU data points illustrates the accuracy of the fitted model predictions. Each of the phase transition-age distributions has its own unique shape and characteristics. Both the S and G2/M phase transition-age probability distributions display a positive bias but with relatively small standard deviations (Table II) though the S phase transition is later and slightly broader than the G2/M phase. This would imply that for the S and G2/M phases: (i) cells require at least a certain amount of time to progress through the phase, (ii) cells move through the phase at nearly identical rates, and (iii) delays due to checkpoints or heterogeneity are either unlikely or
While the flexible spline representation of the phase transition rates does an excellent job of matching the pulse-data, the high number of degrees of freedom may lead to over-fitting. That is, simpler or regularized representations may fit nearly as well; in particular, the two transition age modes seen the in the G0/G1 phase may not be crucial to explaining the data. To this end, other basis functions for the age-dependent transition rates, derived from transition age densities of the lognormal, skewed lognormal, and normalized Gaussian distributions were explored. Since Lagrangian information was not used, the only quantitative method to determine the likelihood of a transition age distributions is the residual fit to the experimental data. A comparison of the least-squares errors, provided in Table I, indicates that the flexible cubic splines functional form does fit the experimental data better than the alternative individual functional forms for determining the age-dependent transition rates. Table II lists the best fit parameters for the alternative functional forms. These differences can be seen directly when the model predictions using the spline and lognormal derived transition rates are compared with the experimental phase percentages (see Fig. 5c and d). Particularly for the BrdU labeled subpopulation, the errors in the single lognormal predictions become progressively larger as time progresses while the spline fit remains accurate. In light of the extremely tight confidence intervals on the experimental data, single distributions are not as adequate as the spline fit. When normalized Gaussian distributions are used for the transition probability with respect to age, the G0/G1 representation (Fig. 6a) becomes distorted into a physiologically unlikely arrangement. The estimated G0/G1 normal distribution is positioned such that most of the young cells
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
967
Figure 5. Comparison of BrdU pulse-labeling subpopulation cell cycle phase experimental data with best-fit flexible spline model predictions. The experimental data in each phase is represented by: *, G0/G1; , S; , G2/M. The thick lines correspond to model predictions when the cell cycle phase transition rates are fit to only the pulse data and the initial phase percentages are defined by the transition rates. The thin lines are best-fit model predictions when (a and b) both pulse and continuous data (Comparisons With Continuous Labeling Data Section) is fit and the initial phase percentages are minimization input parameters or (c and d) single lognormals are fit to pulse data (Determination of Alternative Functional Forms for the Age-Dependent Transition Rates Section). a and c: Pulse-labeling dynamics for the BrdU labeled subpopulation. b and d: Pulse-labeling dynamics for the BrdU unlabeled subpopulation.
Table I. Least-squares fit, E, when the minimization is performed for each distributions examined along with the number of parameters manipulated to attain that fit. Fit to pulse data only: restricted initial condition Flexible cubic splines (defines the transition rates) Normal Lognormal Skewed lognormal 2 lognormals/phase 3 lognormals/phase
968
Total number of E ( 103), fitted pulse parameters data
E ( 103), continuous data
90
1.9553
1.0004
6 6 9 15 24
3.2681 2.5217 2.4832 1.6891 1.6893
1.3472 1.1682 1.1908 1.3655 1.3146
just entering the G0/G1 phase transition quickly while the large variance still enables the transition of older cells: that is, the increased variance leads to a likely over-prediction of younger cell division rates, even though the residual’s error when fit to the data is relatively low. When a lognormal is used as the transition-age distribution (Fig. 6b), the G0/G1 phase distribution is again extremely broad. Skewed lognormal distributions were intuitively appealing because they can enhance the bias of distributions, but this bias is of little use in capturing the bimodal behavior. The skewed lognormal fit is similar to that of the lognormal (Fig. 6c). No single distribution per phase can capture the cubic-spline estimated bimodal behavior of the G0/G1 phase, but rather,
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
Table II.
Summary of the parameter sets attained from the best fits of each of the phase transition representations.
Distribution: restricted initial condition Cubic spline (calculated)
Single transition age distributions Normal
Lognormal
Skewed lognormal
Lognormal basis functions’ transition age distributions 2 lognormal combination
3 lognormal combination
Distribution: unrestricted initial condition 3 lognormal combination ( f1,0þ ¼ [0.00 106.33 0.23], f1,0 ¼ [54.21 18.51 18.07], f1,0 ¼ [54.94 32.79 10.40])
Phase
Mean, m (h)
Standard deviation, s (h)
Weight, w, or asymmetry
G0/G1 S G2/M
16.2627 11.2515 6.3182
9.3577 2.5025 1.7564
G0/G1 S G2/M G0/G1 S G2/M G0/G1 S G2/M
11.31 10.50 5.79 16.59 10.49 5.76 16.70 10.04 5.64
13.39 0.70 0.51 14.75 1.04 0.80 14.15 1.07 0.82
1.00 1.08 0.94
G0/G1 S G2/M G0/G1 S G2/M
½ 7:57 25:61 ½ 10:52 15:22 ½ 6:09 6:26 ½ 7:75 25:39 7:93 ½ 10:47 14:45 10:53 ½ 6:03 5:28 8:77
½ 1:00 2:21 ½ 1:00 1:21 ½ 0:99 1:13 ½ 1:15 1:62 1:29 ½ 1:23 1:08 1:27 ½ 0:91 1:00 1:39
½ 0:55 0:45 ½ 0:92 0:08 ½ 0:59 0:41 ½ 0:50 0:45 0:05 ½ 0:87 0:09 0:04 ½ 0:99 0:01 0:00
G0/G1 S G2/M
½ 25:13 8:79 5:28 ½ 11:22 10:73 5:25 ½ 6:75 6:09 5:40
½ 1:87 1:39 0:99 ½ 1:24 1:65 1:28 ½ 1:30 1:14 0:89
½ 0:36 0:59 0:05 ½ 0:05 0:95 0:00 ½ 0:01 0:25 0:74
Empty items imply that that parameter does not influence that distribution.
this results in a more disperse distribution to accommodate the data as indicated by the larger standard deviations associated with the G0/G1 phase in Table II. Regardless of these differences in the G0/G1 phase, the means and variations seen in the S and G2/M phases for each of the single mode transition age distributions, are similar to those of the flexible splines (Fig. 4 vs. Fig. 6; Table II). What effect does a broad G0/G1 phase versus a bimodal G0/G1 phase have on ability of the model to fit the BrdU pulse data? Comparing the best predictions for the BrdU labeled dynamics, the broad G0/G1 phase does produce oscillations of slightly lower amplitude than the bimodal fit, but this looks to be of marginal significance in relating to the experimental BrdU pulse data (Fig. 5c). In addition, the predictions for the BrdU unlabeled subpopulation using the broad or bimodal G0/G1 phase are nearly indistinguishable (Fig. 5d). Lognormal Basis Functions The predictions of the single mode distributions were improved upon by using a weighted combination of either two or three lognormal distributions to describe the transition age distribution in each phase where the means, standard deviation, and weights are to be determined hlognormal ðtÞ ¼ i w1 hlognormal þ w2 hlognormal i;1 i;2 X þ w3 hlognormal ; wj ¼ 1 i;3 j
This description invites the possibility of multiple subpopulations or more distributed transition-age distributions using fewer parameters. We found that the combinations of lognormals orient themselves to form transition agedistributions that are similar in shape to the flexible spline nodes representation (Fig. 4b vs. Fig. 7a and b). In the S and G2/M phases, the two lognormals in each phase have considerable overlap and the superposition forms distributions with means and standard deviations that are slightly greater those derived using the splines (see Table II). However, there were significant differences in the means of the two lognormals that comprise the G0/G1 phase which results in a transition age probability distribution that is very similar to that found with the splines. Using three lognormal distribution per phase results in predicted labeled and unlabeled phase dynamics are that are similar between when cubic-splines (Fig. 4b) and two lognormals per phase (Fig. 7a) though the number parameters used was greater than five times less than that required for the spline fit.
Comparisons With Continuous Labeling Data For validation of this method, data was obtained from a continuous BrdU label where all cells eventually become BrdU labeled during progression through the S phase. This data explicitly shows the decay in the G0/G1 and G2/M phases. All cell preparation and assay experimental protocols are the same as described earlier for the pulse labeled BrdU experiments with the exception of BrdU labeling of the entire 20 mL flask of cells. (The flask contents is transferred
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
969
Figure 6.
Best fit of cell cycle transition rates to BrdU dynamic data from single transition-age distributions for each phase. a: Gaussian cell cycle phase transition-age distributions. b: Lognormal cell cycle phase transition-age distributions. c: Skewed lognormal cell cycle phase transition-age distributions.
to a centrifuge tube and exposed to 200 mL of the diluted BrdU solution; no cells are left unexposed.) As before, these continuous BrdU labeling experiments were performed in duplicate with extremely consistent results at all time points. During continuous labeling, unlabeled cells entering the S phase become labeled with BrdU as shown in the bivariate BrdU–DNA plots time sequence of Figure 8. The unlabeled population is initially in unbalanced growth. This data is used as an independent verification for the G0/G1 and G2/M predictions based on the pulse data. In light of the continuous labeling observed in the experimental data, Equation (9), which describes the phase number densities of the labeled and unlabeled cells, must be modified to reflect the accumulation of labeled cells. Equation (14) accounts for
970
this accumulation and can be used in lieu of Equation (9) for continuous labeling Labeled ðBrdU þ Þ Subpopulation 1 P ðkÞ n1 ðt; tÞ ¼ n1 ðt; tÞ n2 ðt; tÞ ¼ n3 ðt; tÞ ¼
k¼2 1 P k¼0 1 P k¼1
ðkÞ
n2 ðt; tÞ (14)
ðkÞ
n3 ðt; tÞ
Unlabeled ðBrdU Þ Subpopulation ð0Þ ð1Þ n1 ðt; tÞ ¼ n1 ðt; tÞ þ n1 ðt; tÞ n2 ðt; tÞ ¼ 0 ð0Þ n3 ðt; tÞ ¼ n3 ðt; tÞ
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
Figure 8. A series of bivariate BrdU versus DNA cell density plots illustrate the progression of the labeled and unlabeled subpopulations progression around the cell cycle. Ten thousand cells were collected for each figure with speckled or empty regions representing lower cell densities. a: Immediately after BrdU inoculation (t ¼ 1 h). b: Ten hours after BrdU inoculation. The arrow denotes movement of BrdU unlabeled cells becoming labeling during the S phase. c: Twenty-four hours after BrdU inoculation.
Figure 7. Linear combinations of lognormal distributions used to approximate the residence time distribution of each cell cycle phase. The added lognormal distributions give additional flexibility to the G0/G1, S and G2/M distributions. a: Two lognormal distributions per phase. b: Three lognormals per phase.
of balanced growth using all available data, 33 parameters are fitted to three lognormals per phase. The identified transition rates are extremely consistent (Fig. 4b) from when the pulse data is used alone, and the fit to the pulse data is also similar (see Fig. 5a and b) with a least squares fit error to the pulse data of 1.1972 103. The fit to the continuous BrdU data is also improved (Fig. 9c and d) with a least squares fit error of 1.08163 103.
Conclusions The pulse data derived transition rates for both the flexible spline and lognormal distributions were used to predict the dynamics of the continuous BrdU. The least squares error quantifying the predicted model dynamics fit to the continuous experimental data are provided in Table I for all assumed distribution forms. As illustrated in Figure 9 the spline and lognormal distribution forms both capture the trends of the continuous labeling experimental data. Throughout the continuous BrdU labeling experiment, the ratio of labeled cell fluxes into and out of the S phase remains constant hence, it impossible to determine the age structure of the S phase transition rates using this data alone. When the pulse-data (two independent subpopulations) is supplemented with the continuous data (one subpopulation) for fitting spline transition rates without the restriction
A methodology has been developed that enables predictions of cell cycle phase age-dependent transition rates based on experimental data with BrdU labeling where the crux of the identification relies on comparison of the BrdU data with model predictions based entirely on the transition rates. For making model predictions based on transition rates, an agestructured cell cycle model is presented in which the agedependent transition rates are posed. Model predictions require knowledge of initial age-distributions, which are unidentifiable directly, so an expression that yields initial age-distributions from transition rates is derived under the assumption of a balanced growth setting. Since the transition rates and an initial condition are sufficient for model identification, and the transition rates imply the initial condition, the age-dependent transition rates are the only
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
971
Figure 9. Comparison of BrdU continuous labeling subpopulation cell cycle phase experimental data with best-fit flexible spline model predictions. The experimental data in each phase is represented by: *, G0/G1; , S; , G2/M. The thick lines correspond to model predictions when the cell cycle phase transition rates are fit to only the pulse data and the initial phase percentages are defined by the transition rates. The thin lines are model predictions when (a and b) single lognormals are fit to pulse data while in (c and d) both pulse and continuous data is fit and the initial phase percentages are minimization input parameters. a and c: Continuous labeling dynamics for the labeled population. b and d: Continuous labeling dynamics for the unlabeled population.
input for making model predictions. We then illustrate how pulse BrdU labeling, which introduces instances of unbalanced growth while maintaining individual cells and the total population in balanced growth, is sufficient for the identification of the unique set of transition rates from the unbalanced growth dynamics. The identified age transition distributions were validated through a comparison of the predicted dynamics with those observed for continuous BrdU labeling. In identifying the transition rates, a general representation of the transition rates using a series of cubic splines, where both the spline node values as well as locations were minimization parameters, resulted in transition probability distributions with both expected and unexpected features. Very few transitions occur before a threshold age and after
972
which there is a window of ages when transitions are very likely; however, for the G0/G1 phase, the best fit transition age probability distribution suggested multiple subpopulations through its bimodal form. This bimodal behavior may indicate the presence of multiple subpopulations resulting from either continuous, systematic differences between cells or individual, discrete events. With multiple distinct modes, each subpopulation is governed by an independent set of cell cycle kinetics. Or, the subpopulations could result from a random event which may or may not occur each time through the cycle such as a checkpoint control. While the existence of multiple subpopulations is able to be identified and their behavior characterized using age-dependent transition rates, a Lagrangian approach where individual cells are followed would be required to distinguish between
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
the two mechanisms. Fit of the age probability distribution to the data assuming single normalized Gaussian, lognormal, and skew-lognormal transition probability distributions was reasonable, but not as accurate due to the inability to capture the G0/G1 bimodal transition age distribution. But, describing each phase transition with a convex combination of two or three lognormal distributions improved the overall fit of the BrdU data relative to the spline representation while requiring only 15 or 24 parameters rather than the 90 used for the spline fit. We see this approach as an efficient alternative with adequate flexibility to capture the major features of data while minimizing the number of parameters to be identified. The foregoing inversion involved no regularization as the error bounds on the data were very small. However, larger error bounds would call for the use of appropriate regularization (Tikhonov, 1977).
R1
0 fi ðtÞ dt ¼ 1) will give the rates of change in the total number density Z 1 Z 1 dNi ðtÞ dfi ðtÞ þ Ni ðtÞ dt ¼ Ni ðtÞ Dout i;i ðtÞfi ðtÞ dt dt dt 0 0
The boundary and containment conditions can be used to simplify the term Z 1 XZ 1 dfi ðtÞ dt ¼ fi ðtÞj1 fi ðtÞj0 ¼ Din i;j fj ðtÞ dt dt 0 0 j to give Equation (3) dNðtÞ ¼ DNðtÞ dt where Dij ¼
The authors are grateful for the assistance of the Purdue University Cytometry Laboratories and support from the Purdue Cancer Center.
Derivation of Initial Age-Distributions (Age-Distribution During Balanced Growth) Balanced growth was earlier defined as the growth condition that satisfies two criterions: (i) intransient growth conditions and (ii) time-invariant cell-state distributions. Timeindependent transition rates satisfy the first condition, and the second condition implies that the age-distribution in the ith cell cycle phase is described by the time-independent distribution fi ðtÞ when balanced growth is achieved. So, during balanced growth, the number density in the ith phase is ni ðt; tÞ ¼ Ni ðtÞfi ðtÞ (Eq. (2)). Substitution of this balanced growth number density into the population balance equation for the ith phase (Eq. (1)) @ni ðt; tÞ @ni ðt; tÞ þ ¼ Dout i;i ðtÞni ðt; tÞ @t @t
out ½Din ij ðtÞ Dij ðtÞ fj ðtÞ dt
0
dNi ðtÞ ¼ mNi ðtÞ dt The use of this expression in Equation (15) allows the total number density to be removed, leaving a differential equation (Eq. 4) for change in the age-distribution with age 3 P i ðtÞ fi ðtÞmNi ðtÞ þ Ni ðtÞ dfdt ¼ Dout ij ðtÞNj ðtÞfj ðtÞ j¼1 " # 3 P dfi ðtÞ out ðDij ðtÞfj ðtÞÞ þ mfi ðtÞ dt ¼ j¼1
dfi ðtÞ dt
¼ ðDout ii ðtÞ þ mÞfi ðtÞ
The general solution of this type of differential equation is
Z t 0 0 ðDout ðt Þ þ mÞ dt fi ðtÞ ¼ fi ð0Þ exp ii
gives @ðNi ðtÞfi ðtÞÞ @ðNi ðtÞfi ðtÞÞ þ ¼ Dout i;i ðtÞNi ðtÞfi ðtÞ @t @t
0
Since the total product density is age-independent and the balanced growth age-distribution is time-independent, this equation simplifies to dNi ðtÞ dfi ðtÞ þ Ni ðtÞ ¼ Dout i;i ðtÞNi ðtÞ fi ðtÞ dt dt
1
The time independent age-distribution gives a timeindependent coefficient matrix D. The specific growth rate of each cell cycle phase during balanced growth will be the positive real eigenvalue of this matrix, m, so
Appendix
fi ðtÞ
Z
where the regulatory condition
0
fi ðtÞ ¼ 0 defines fi ð0Þ
Rt R1 0 0 fi ðtÞ dt ¼ 1 ¼ fi ð0Þ 0 exp 0 ðDout ii ðt Þ þ mÞ dt dt i ) fi ð0Þ ¼ R 1 h R 1 1 R1 0
0
(15)
A further implication of the time-independent agedistributions is identical specific growth rates for each cell cycle phase during balanced growth. Returning to the population balance equation, substitution of the balanced growth separated number density ni ðt; tÞ ¼ Ni ðtÞfi ðtÞ into Equation (1) and integrating over all ages (where
R1
exp
ðDout ii ðt0ÞþmÞ dt0 dt
0
So, the age-distribution in the ith cell cycle phase during balance growth is given by Equation (5) fi ðtÞ ¼
Z
Z t0 1 00 00 exp ðDout ðt Þ þ mÞ dt dt 0 ii 0
Z t 0 0 0
exp ðDout ðt Þ þ mÞ dt ii 1
0
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
973
References Arino O, Kimmel M. 1987. Asymptotic analysis of a cell cycle model based on unequal division. SIAM J Appl Math 47:128–145. Baisch H, Otto U. 1993. Intratumoral heterogeneity of S phase transition in solid tumours determined by bromodeoxyuridine labelling and flow cytometry. Cell Prolif 26:439–448. Baisch H, Otto U, Hatje U, Fack H. 1995. Heterogeneous cell kinetics in tumors analyzed with a simulation model for bromodeoxyuridine single and multiple labeling. Cytometry 21:52–61. Bernard S, Pujo-Menjouet L, Mackey MC. 2003. Analysis of cell kinetics using a cell division marker: Mathematical modeling of experimental data. Biophys J 84:3414–3424. Chiorino G, Metz JA, Tomasoni D, Ubezio P. 2001. Desynchronization rate in cell populations: Mathematical modeling and experimental data. J Theor Biol 208:185–199. Cooper S. 2003. Rethinking synchronization of mammalian cells for cell cycle analysis. Cell Mol Life Sci 60:1–8. Cooper S. 2004. Is whole-culture synchronization biology’s ‘‘perpetualmotion machine’’? Trends Biotechnol 22:266–269. Darzynkiewicz Z, Gong J, Juan G, Ardelt B, Traganos F. 1996. Cytometry of cyclin proteins. Cytometry 25:1–13. Dien BS, Srienc F. 1991. Bromodeoxyuridine labeling and flow cytometric identification of replicating Saccharomyces cerevisiae: Lengths of cell cycle phases and population variability at specific cell cycle positions. Biotechnol Prog 7:291–298. Del Bino G, Lassota P, Darzynkiewicz Z. 1991. The S-phase cytotoxicity of camptothecin. Exp Cell Res 193:27–35. Dolbeare F, Gratzner H, Pallavicini MG, Gray JW. 1983. Flow cytometric measurement of total DNA content and incorporated bromodeoxyuridine. Proc Natl Acad Sci USA 80:5573–5577. Dyson J, Villella-Bressan R, Webb G. 2000. A nonlinear age and maturity structured model of population dynamics. I. Basic Theory. J Math Anal Appl 242:93–104. Dyson J, Villella-Bressan R, Webb GF. 2002. Asynchronous exponential growth in an age structured population of proliferating and quiescent cells. Math Biosci 177–178:73–83. Eakman JM, Fredrickson AG, Tsuchiya HM. 1966. Statistics and dynamics of microbial cell populations. Chem Eng Prog Symp Ser 62:37–49. Fredrickson AG, Mantzaris NV. 2002. A new set of population balance equations for microbial and cell cultures. Chem Eng Sci 57:2265–2278. Fredrickson AG, Tsuchiya HM. 1963. Continuous propagation of microorganisms. AIChE J 9:459–468. Fredrickson AG, Ramkrishna D, Tsuchiya HM. 1967. Statistics and dynamics of procaryotic cell populations. Math Biosci 1:327–374. Hjortsø MA. 2006. Population balances in biomedical engineering: Segregation through the distribution of cell states. New York: McGraw-Hill. 182 p. Ja´ nossy J, Ubezio P, Apa´ ti A´ , Mago´ csi M, Tompa P, Friedrich P., 2004. Calpain as a multi-site regulator of cell cycle. Biochem Pharmacol 67:1513–1521. Kromenaker SJ, Srienc F. 1991. Cell-cycle-dependent protein accumulation by producer and nonproducer murine hybridoma cell lines: A population analysis. Biotechnol Bioeng 38:665–677. Kutalik Z, Razaz M, Elfwing A, Ballagi A, Baranyi J. 2005. Stochastic modelling of individual cell growth using flow chamber microscopy images. Int J Food Microbiol 105:177–190. Liou JJ, Srienc F, Fredrickson AG. 1997. Solutions of population balance models based on a successive generations approach. Chem Eng Sci 52:1529–1540.
974
Liou JJ, Fredrickson AG, Srienc F. 1998. Selective synchronization of tetrahymena pyriformis cell populations and cell growth kinetics during the cell cycle. Biotechnol Prog 14:450–456. Liu SV. 2005. Debating cell-synchronization methodologies: Further points and alternative answers. Trends Biotechnol 23:9–10. Lyons AB. 1999. Tracking cell proliferation with carboxyfluorescein diacetate succinimidyl ester. Immunol Cell Biol 77:509–515. Mantzaris NV, Liou JJ, Daoutidis P, Srienc F. 1999. Numerical solution of a mass structured cell population balance model in an environment of changing substrate concentration. J Biotechnol 71:157– 174. Mantzaris NV, Srienc F, Daoutidis P. 2002. Nonlinear productivity control using a multi-staged cell population balance model. Chem Eng Sci 57:1–14. Nordon RE, Nakamura M, Ramirez C, Odell R. 1999. Analysis of growth kinetics by division tracking. Immunol Cell Biol 77:523– 529. Nova´ k B, Tyson J. 1993. Modeling the cell division cycle: M-phase trigger, oscillations, and size control. J Theor Biol 165:101–134. Nova´ k B, Tyson J. 2004. A model for restriction point control of the mammalian cell cycle. J Theor Biol 230:563–579. Ormerod MG, Kubbies M. 1992. Cell cycle analysis of asynchronous cell populations by flow cytometry using bromodeoxyuridine label and Hoechst-propidium iodide stain. Cytometry 13:678–685. Ramkrishna D. 2000. Population balances: Theory and applications to particulate systems in engineering. San Diego: Academic Press. 355 p. Sena G, Onado C, Cappella P, Montalenti F, Ubezio P. 1999. Measuring the complexity of cell cycle arrest and killing of drugs: Kinetics of phasespecific effects induced by taxol. Cytometry 37:113–124. Sennerstam R, Stro¨ mberg , JO. 1995. Cell cycle progression: Computer simulation of uncoupled subcycles of DNA replication and cell growth. J Theor Biol 175:177–189. Shackney SE. 1973. A cytokinetic model for heterogeneous mammalian cell populations. J Theor Biol 38:305–333. Shackney SE, Shankey TV. 1999. Cell cycle models for molecular biology and molecular oncology: Exploring new dimensions. Cytometry 35:97– 116. Sherer E, Hannemann RE, Rundell A, Ramkrishna D. 2006. Analysis of resonance chemotherapy in leukemia treatment via multi-staged population balance models. J Theor Biol 240:648–661. Steel GG. 1977. Growth kinetics of tumours: Cell population kinetics in relation to the growth and treatment of cancer. Oxford: Clarendon Press. 351 p. Subramanian G, Ramkrishna D. 1971. On the solution of statistical models of cell populations. Math Biosci 10:1–23. Subramanian G, Ramkrishna D, Fredrickson AG, Tsuchiya HM. 1970. On the mass distribution model for microbial cell populations. Bull Math Biophys 32:521–537. Tikhonov AN. 1977. Solution of ill-posed problems. Washington: Winston & Sons. 272 p. Tyrcha J. 2001. Age-dependent cell cycle models. J Theor Biol 213:89– 101. Tyrcha J. 2004. Cell cycle progression. C.R. Biol 327:193–200. Webb GF, Grabosch A. 1987. Asynchronous exponential-growth in transition-probability models of the cell-cycle. SIAM J Math Anal 18:897– 908. Zamamiri AM, Zhang Y, Henson MA, Hjortsø MA. 2002. Dynamic analysis of an age-distribution model of oscillating yeast cultures. Chem Eng Sci 57:2169–2181.
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
School of Chemical Engineering, Forney Hall of Chemical Engineering, 480 Stadium Mall Way, Purdue University, West Lafayette, Indiana 47907; telephone: 765-494-4066; fax: 765-494-0805; e-mail: [email protected] 2 Chemical and Biological Engineering Department, University of Wisconsin, Madison, Wisconsin 3 Weldon School of Biomedical Engineering, Purdue University, West Lafayette, Indiana Received 7 May 2007; revision received 19 July 2007; accepted 9 August 2007 Published online 4 September 2007 in Wiley InterScience (www.interscience.wiley.com). DOI 10.1002/bit.21633
Introduction ABSTRACT: A methodology is developed that determines age-specific transition rates between cell cycle phases during balanced growth by utilizing age-structured population balance equations. Age-distributed models are the simplest way to account for varied behavior of individual cells. However, this simplicity is offset by difficulties in making observations of age distributions, so age-distributed models are difficult to fit to experimental data. Herein, the proposed methodology is implemented to identify an age-structured model for human leukemia cells (Jurkat) based only on measurements of the total number density after the addition of bromodeoxyuridine partitions the total cell population into two subpopulations. Each of the subpopulations will temporarily undergo a period of unbalanced growth, which provides sufficient information to extract age-dependent transition rates, while the total cell population remains in balanced growth. The stipulation of initial balanced growth permits the derivation of age densities based on only agedependent transition rates. In fitting the experimental data, a flexible transition rate representation, utilizing a series of cubic spline nodes, finds a bimodal G0/G1 transition age probability distribution best fits the experimental data. This resolution may be unnecessary as convex combinations of more restricted transition rates derived from normalized Gaussian, lognormal, or skewed lognormal transition-age probability distributions corroborate the spline predictions, but require fewer parameters. The fit of data with a single log normal distribution is somewhat inferior suggesting the bimodal result as more likely. Regardless of the choice of basis functions, this methodology can identify age distributions, age-specific transition rates, and transition-age distributions during balanced growth conditions. Biotechnol. Bioeng. 2008;99: 960–974. ß 2007 Wiley Periodicals, Inc. KEYWORDS: population balance; inverse model; cell division; BrdU; cell generations
The growth behavior of an individual cell depends on a complex system of molecular interactions, feedback regulations, and synthesis reactions (Darzynkiewicz et al., 1996; Ja´nossy et al., 2004; Nova´k and Tyson, 1993, 2004). This complexity can lead to processes that are distributed throughout the life cycle (Liou et al., 1998), such as protein (Kromenaker and Srienc, 1991) or RNA (Arino and Kimmel, 1987) production and susceptibility to certain chemotherapeutic agents (Del Bino et al., 1991). Researchers are beginning to exploit the non-uniformity of these processes with the aid of predictive models (Mantzaris et al., 2002; Sherer et al., 2006). The heterogeneity of the resulting cell population is often represented in models by a single or small set of characteristic variables. Common characteristic variables include: cell mass (Mantzaris et al., 1999; Subramanian et al., 1970; Subramanian and Ramkrishna, 1971), volume (Webb and Grabosch, 1987), DNA content, protein content (Kromenaker and Srienc, 1991), or age (Sennerstam and Stro¨mberg, 1995; Tyrcha, 2001; Zamamiri et al., 2002). A population balance model describes the cellular dynamics based on growth and transition rates that are functions of the characteristic variables (Fredrickson and Tsuchiya, 1963; Fredrickson et al., 1967; Fredrickson and Mantzaris, 2002). Age, the amount of time since a cell last made a transition from one phase of the cell cycle to the next, is a simple way to account for the varied behavior of individual cells. That is, the amount of time spent by cells in any one phase of the cell cycle is distributed, so the likelihood of a cell transitioning from one phase to the next is dependent upon its age. An
Correspondence to: D. Ramkrishna
960
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008
ß 2007 Wiley Periodicals, Inc.
age-dependent transition rate can be used to describe this variation. For example, in the cell cycle, certain intracellular tasks must be completed before a cell will proceed to the next phase, so it is unlikely that a very young cell will transition. Also, checkpoints or other disturbances can shift and bias the transition age. Lognormal distributions are most commonly used to define transition-age distributions of cell cycle phases (Steel, 1977). The bias is captured even more closely if linear combinations of distributions (Shackney and Shankey, 1999) or asymmetry factors (Baisch et al., 1995) are incorporated. Reaction probability density functions provide a more mechanistic view of cellular transitions (Hjortsø, 2006) so long as the rate of each step in the reaction sequence is known. The advantages of utilizing age as the characteristic variable to describe the transition rates are diminished by the inability to easily measure the age distributions within a population of cells. While inverse methods are available to extract individual behavior information from population observations, these require quantitative measurements of the characteristic variables (Ramkrishna, 2000; Tyrcha, 2004). Synchronizing cells would allow for age approximations, but the synchronization process can lead to altered growth patterns (Chiorino et al., 2001; Cooper, 2003, 2004; Liou et al., 1998; Liu, 2005). Continuous observation of cell cultures is currently the method of choice for extracting individual age-dependent behavior (Kutalik et al., 2005), but this method is tedious as a significant number of cells must be observed to produce relevant distributions. In this paper, we develop an integrated experimental and analytical technique to overcome the difficulty of extracting the age distributions and the associated age-dependent transition rates from experimental data. Previous derivations have estimated age distributions assuming balanced growth conditions and a fixed, constant valued transition rate (Baisch and Otto, 1993; Bernard et al., 2003; Shackney, 1973; Steel, 1977). For a clear understanding of how models are identified using our methodology, it is necessary to recognize that the approach to balanced exponential growth of the overall population is preceded by two stages: (i) physiological changes in individual cells and (ii) a transient period of population adjustment to the time-invariant distribution characteristic of the balanced growth phase. The second stage represents a situation when stage (i) has been completed, but the population distributions in the different cell cycle phases are evolving to the constant proportion in balanced growth. Stage (i) can be achieved by maintaining a constant culture environment and is characterized by timeindependent transition rates. Balanced growth is achieved upon completion of stage (ii) so the age densities are also time independent (Dyson et al., 2000; Fredrickson et al., 1967). This feature of balanced growth implies that the system has reached a state such that the population densities of the phases will maintain constant ratio to one another (Dyson et al., 2002). Herein, we show that assuming balanced growth enables the calculation of age distributions
based only on knowledge of the transition rates. Thus, the only requirement for attaining the age distributions associated with the cell cycle dynamics are the transition rates. A problem arises while trying to determine the vectoral unknowns describing the age distributions associated with each phase of the cell cycle from scalar population density data. We demonstrate how an experimental strategy that labels a subpopulation can overcome this difficulty. In utilizing this method, cells are subjected to pulse labeling with bromodeoxyuridine (BrdU). Only cells in the S-phase of the cell cycle at the time of BrdU application will be labeled with BrdU (Dien and Srienc, 1991; Dolbeare et al., 1983; Sena et al., 1999). This creates two independent subpopulations, one labeled and the other unlabeled, both of which continue around the cell cycle. Two important features of this technique are: (i) that each subpopulation begins in unbalanced growth as their age distributions are changing with time and (ii) the BrdU does not affect the growth behavior of individual cells, so the total population remains in balanced growth (Baisch and Otto, 1993; Ormerod and Kubbies, 1992). This transient behavior is a feature of how the total population densities in each of the cell cycle phases are evolving to the eventual constant proportion in fully balanced growth. Thus, the measurement of labeled and unlabelled population densities provides the vectoral transient data that has the potential to identify vectoral transition rates without explicitly measuring age distributions. While BrdU is used herein to determine the age-structured cell cycle phase transition rates, similar experimental approaches can be used to determine cell division rates. The labeled mitosis method using 3H-thymidine (Steel, 1977) or the transient cell generation distributions from pkh dyes or carboxyfluorescein diacetate succinimidly ester (Lyons, 1999; Nordon et al., 1999) also produce transient data that can be used to determine the age-structured cell division rate. In presenting the methodology for such identification, the model and mathematical foundation are first developed and followed by an in vitro illustration for its use. The Analytical Foundation Section presents the cell growth model describing BrdU labeling with a proof of extractions of transition rates from subpopulations in unbalanced growth while the total populated maintains balanced growth. The Experimental Materials and Methods Section outlines the experimental materials and methods used to collect in vitro data. The age-dependent transition rates that result from the experimental BrdU results are discussed in the Results and Discussion Section, and the fitness is compared with normalized Gaussian, lognormal, and skew lognormal age density derived transition rates. The Conclusions Section summarizes the major conclusions and findings.
Analytical Foundation This section introduces an age-structured cell cycle model (Model Structure and Derivation of Initial Conditions
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
961
Sections), proves that the model predictions are sufficient to identify age-structured transition rates using BrdU labeling (Extraction of Age-Dependent Parameters From Partitioned Population Data Section), and shows the simulation methodology for comparison with the BrdU pulse-labeling experiments (Implementation Section).
The age density is subject to the initial condition nðt; 0Þ ¼ n0 ðtÞ and a boundary condition that consists of cells that have just transitioned out of a phase to become new cells in the next phase. Z 1 nð0; tÞ ¼ Din ðtÞnðt; tÞ dt 0
Model Structure
where
The division cycle of a cell is discretized into three phases in our model: (1) a combined G0 and G1 phase, (2) an S phase, and (3) a combined G2 and M phase (see Fig. 1). Newly divided cells may either be in the quiescent G0 phase or preparing to undergo DNA synthesis in the G1 phase. These two states were indistinguishable experimentally and are combined into a single phase. The S phase is a period when the DNA and many proteins are replicated before cell division. Cells then enter the second gap phase, G2, before physically dividing by undergoing mitosis, M. As the G2 and M states were also indistinguishable experimentally, they are lumped into a single phase. The cells were maintained in a nutrient rich environment throughout the in vitro experiments so a low level of quiescence is anticipated. The number density and transition rate of cells in the ith phase can be written as ni(t, t) and Gi(t) respectively where t is the age, or time elapsed since entering the current phase, and t is the time. It is assumed that the transition rate is stable, or time independent, in balanced growth conditions but the number density will increase with time as cells divide. Neglecting cell death, the population balance equation can be written as @nðt; tÞ @nðt; tÞ þ ¼ Dout ðtÞnðt; tÞ @t @t
(1)
2
3 n1 ðt; tÞ nðt; tÞ ¼ 4 n2 ðt; tÞ 5 n3 ðt; tÞ and G1 ðtÞ Dout ðtÞ ¼ 4 0 0
0 0 Din ðtÞ ¼ 4 G1 ðtÞ 0 0 G2 ðtÞ
3 0 0 G2 ðtÞ 0 5 0 G3 ðtÞ
3 2G3 ðtÞ 0 5 0
Under balanced growth conditions, the age distribution is time invariant so that the ith phase number density can be separated into the product of the population density, Ni(t), and the age-density, fi(t) Z 1 fi ðtÞ ¼ 1 (2) ni ðt; tÞ ¼ Ni ðtÞfi ðtÞ; where 0
Substitution of Equation (2) into the population balance equation (Eq. (1)) and integration over all ages gives an ordinary differential equation for the population increase in each of the phases (see Appendix) dNðtÞ ¼ DNðtÞ dt
(3)
where Dij ¼
where
2
2
Z 0
1
h
i out Din ðtÞ D ðtÞ fj ðtÞ dt ij ij
Only two quantities are required for simulating this system: transition rates and initial number densities. As shown in the Derivation of Initial Conditions Section, the balanced growth age distributions and associated number densities can be calculated using only the transition rates under balanced growth conditions. This implies that knowledge of the transition rates is sufficient for complete model characterization.
Derivation of Initial Conditions Under balanced growth, the matrix of rate constants, D, is time invariant and its positive real eigenvalue, m, is the specific growth rate for both the total population and each cell cycle phase. The system of ordinary differential equations can be written in terms of m and the corresponding eigenvector, h
Figure 1. Phases of cell cycle model. Cells transition out of the i th phase according to an age-dependent transition rate (Gi) with cell division occurring via the G2/M to G0/G1 transition.
962
mh ¼ Dh P where the normalized eigenvector hi = j hj is the fraction of cells in the ith cell cycle phase. To find the age-distribution within each phase, the partial differential equation system is first reduced to ordinary by
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
substitution of the separated number density (Eq. (2)) into the population balance equation (Eq. (1)). As the total number of cells in each phase increases according to the same overall specific growth rate in balanced growth, ordinary differential equations for the change in the age distribution with age are obtained (see Appendix)
" # 3 X dfi ðtÞ out ¼ ðDij ðtÞfj ðtÞÞ þ mfi ðtÞ dt j¼1
(4)
Using the normalization condition stated in Equation (2), the age distribution in each phase can be written as Z
Z exp
1 00 00 ðDout ðt Þ þ mÞ dt dt0 ii 0 0 (5)
Z t out 0 0 exp ðDii ðt Þ þ mÞ dt
fi ðtÞ ¼
1
0
t0
So, during balanced growth, the age-distributed number density is given as X h ni ðt; tÞ ¼ P i fi ðtÞ Nj ðtÞ j hj j
(6)
This expression can be used as an initial condition before a perturbation is introduced through the addition of BrdU. While the Model Structure Section describes an agestructured cell cycle model subject to the initial agedistributions derived in the Derivation of Initial Conditions, the Extraction of Age-Dependent Parameters From Partitioned Population Data Section addresses the question as to whether BrdU labeling experimental data is sufficient to identify the transition rates or age-distributions. Extraction of Age-Dependent Parameters From Partitioned Population Data In a balanced growth condition there is no unique set of transition rates and initial conditions that describes the total
Figure 2. Decomposing the Jurkat cells into cell cycle phases with BrdU labeling the S phase cells 1 h after labeling. ModFit software is used to determine the percentage of cells in each phase. a: A bivariate BrdU versus DNA content density plot is used to distinguish cells that incorporated BrdU (cells above the drawn line) from those that did not (cells below the line). b: A DNA histogram of the total population is used to determine the percentage of cells in each phase for the overall cell population. c: DNA histogram of labeled (BrdUþ) cells. d: DNA histogram of unlabeled (BrdU) cells.
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
963
population data. Unbalanced growth data with a known initial age distribution and minimal fluctuations in transition rates are needed for their identification, and BrdU pulse-labeling of a balanced growth population satisfies all of these requirements. BrdU labeling (Fig. 2a) partitions the balanced growth population (Fig. 2b) into two subpopulations (Fig. 2c and d). Each of these subpopulations will undergo unbalanced growth while the total population continues in balanced growth with unaltered transition rates. BrdU is an analog to the DNA base thymidine, a uridine derivative, such that cells that are actively synthesizing DNA, the S phase cells, will incorporate greater amounts of BrdU than cells in the G0/G1 or G2/M phases. By pulsing, or removing excess BrdU after a short time period, only those cells that were in the S phase at the time of application, t ¼ 0, will incorporate the BrdU. (We found that direct removal of all excess BrdU from the media proved difficult experimentally. However, labeling an aliquot of the original population, removing excess BrdU, and then recombining with the remainder of the original cell population does result in pulse labeling. The BrdU unexposed cells likely absorb residual BrdU but in small enough amounts to appear as BrdU unlabeled.) This labeling process may be represented by a pre-multiplying projection matrix, L, which concurrently produces a second projection matrix, I L, so that there are two subpopulation distributions among the different phases described by Ln0(t) and (I L)n0(t). The BrdU technique produces the following 2
0
3
6 7 Ln0 ðtÞ ¼ 4 an0;2 ðtÞ 5; 0
2
6 ðI LÞn0 ðtÞ ¼ 4
an0;1 ðtÞ þ ð1 aÞn0;1 ðtÞ
3 (7)
7 ð1 aÞn0;2 ðtÞ 5 an0;3 ðtÞ þ ð1 aÞn0;3 ðtÞ
where the variations in age-density with time lead to a transition matrix, D, that now also becomes a function of time. The population at the time of labeling, N0, is in balanced growth, so it is contained purely in the eigenspace, Em, which corresponds to the specific growth rate eigenvalue, m. Another way of stating this is to say that the projection, Pm, of all three-dimensional vectors onto Em maps N0 into N0 (i.e., PmN0 ¼ N0). This property clearly cannot be shared by LN0 and (I L)N0 simultaneously as PmL 6¼ L. Thus, we conclude that the transients produced by either of the initial conditions in Equation (8) will be truly three-dimensional and, in principle, sufficient for the identification of the three transition rates contained in the matrix Dout(t). The pulselabeling of cells with BrdU is required for this identification to take place.
Implementation After calculation of the initial number densities, transient number densities are computed using a variation of the method of characteristics (Ramkrishna, 2000): the method of successive generations (Liou et al., 1997). The number density for the ith phase, ni(t, t), is discretized into a series of ðkÞ k number densities, ni , with k representing the number of transitions undergone or generations. The BrdU labeling further separates each initial number density into labeled and unlabeled portions so the initial generation (k ¼ 0) number densities are given by Equation (7). For example, all cells beginning the S phase (i ¼ 2, k ¼ 0) are labeled and will move to the G2/M (i ¼ 3, k ¼ 1) then G0/G1 (i ¼ 1, k ¼ 2) phases before returning to the S phase (i ¼ 2, k ¼ 3). Following this pattern, the total phase number density of labeled and unlabeled cells will be Labeled ðBrdU þ Þ Subpopulation
where a is the fraction of the total population exposed to BrdU. (In the BrdU experiments the aliquot of BrdU exposed cells is added to a population not exposed to BrdU—we found this helped achieve true pulse-labeling. The unlabeled boundary condition in Equation (7) reflects this mixture.) It is now of interest to follow the population distributions obtained by solving Equation (1) for the two initial conditions stated above. However, since only the total, rather than age-distributed, population density is measurable, the equations of interest are only those in the total population density (Eq. (3)), for the two initial conditions as shown below
n1 ðt; tÞ ¼ a n2 ðt; tÞ ¼ a n3 ðt; tÞ ¼ a
k¼0 1 P
ð3kþ2Þ
n1
ðt; tÞ
ð3kÞ
n2 ðt; tÞ ð3kþ1Þ
n3
ðt; tÞ
Unlabeled ðBrdU Þ Subpopulation n1 ðt; tÞ ¼ n2 ðt; tÞ ¼
1 P k¼0 1 P k¼0 1 P k¼0
964
k¼0 1 P
k¼0
n3 ðt; tÞ ¼ dNðtÞ ¼ DðtÞNðtÞ; dt LN0 Nð0Þ ¼ ðI LÞN0
1 P
ð3kÞ
ð3kþ1Þ
ðn1 ðt; tÞþn1 ð3kÞ
ð3kþ1Þ
ðð1 aÞn2 ðt; tÞþn2 ð3kÞ
ð3kþ2Þ
ðt; tÞÞ
ð3kþ2Þ
ðt; tÞÞ
ðt; tÞþð1 aÞn1 ðt; tÞþn2
ð3kþ1Þ
ðn3 ðt; tÞþð1 aÞn3
ð3kþ2Þ
ðt; tÞþn3
ðt; tÞÞ (9)
labeled population unlabeled population
(8) The number densities for each of the generations can be found by backtracking along the characteristics until either
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
the initial condition, for the k ¼ 0 generation, or boundary condition, for the k > 0 generations, is reached Initial Conditions : ( Rt ð0Þ ð0Þ n ðt; tÞ ¼ ni ðt t; 0Þ exp tt Gi ðt 0 Þ dt 0 t t iðkÞ ni ðt; tÞ ¼ 0; k ¼ 1; 2; . . . Boundary Conditions : ( ð0Þ ni ðt; tÞ ¼ 0 t>t Rt ðkÞ ðkÞ ni ðt; tÞ ¼ ni ð0; t tÞexp 0 Gi ðt0 Þ dt 0 ; k ¼ 1; 2; . . . (10)
The initial cell distribution is known, and the boundary conditions are calculated successively as the boundary cell density depends on the final solution in the previous phase nðkÞ ð0; tÞ ¼
Z
1
Din ðt; tÞnðk1Þ ðt; tÞ dt
(11)
0
The final number densities are calculated and the appropriate combinations of these final number densities (Eq. (9)) describe the dynamics of the subpopulations.
Experimental Materials and Methods Materials and Methods Jurkat cells were grown in 75 cm2 flasks (Corning 430725) at 378C and 5% CO2 in 20 mL of a 90/10 mixture of Iscove’s Modified Dulbecco’s Medium (Sigma I3390) and fetal bovine serum (Sigma F2442) supplemented with 1 mL Lglutamine (Sigma G2150)/10 mL medium. Exponential growth was maintained for a week prior to the addition of BrdU by passing the cells every 2 days so concentrations remained between 105 and 106 cells/mL. One day before the addition of BrdU, cells were passed to create two flasks at a concentration of approximately 2 105 cells/mL in 20 mL of media. Also 1 day before pulse labeling, 20 mL of complete media was placed in the incubator to pre-warm and pre-gas in preparation for re-suspension after pulse labeling. The BD Biosciences BrdU Flow Kit (559619) was used for the BrdU procedures. A 10 mg/mL stock BrdU in 1 DPBS solution was diluted to 1 mM by adding 31 mL of the stock to 1 mL of 1 DPBS. For pulse-labeling, a 4 mL aliquot is exposed to 40 mL of the diluted BrdU solution (10 mM final BrdU concentration) in a 50 mL centrifuge tube. Forty-five minutes after inoculation, the remaining 16 mL of cells were transferred to the centrifuge tube (so a, the fraction of unlabeled cells, ¼0.2) which was then spun at 1,000 rpm for 5 min before resuspending the cells in 20 mL of pre-warmed and pre-gassed complete media and transferring to the original flask. Samples of 1 mL were collected every 2–4 h for approximately 40 h beginning with the newly resuspended cells. Samples were washed, after centrifuging for 5 min at 1,000 rpm and discarding the supernatant, in 1 mL of
Pharmingen Stain Buffer (PSB); resuspended in 100 mL of BD Cytofix/Cytoperm buffer for 15 min at room temperature; washed in a 1 mL of 1 Perm/Wash solution; and given a final wash in 1 mL of PSB. The cells were preserved by re-suspending in 1 mL of PBS and storing in a 48C freezer. Cytometry analysis was performed within a week of sample collection. In preparation for cytometry analysis, the samples were centrifuged, resuspended in 100 mL Cytofix/Cytoperm buffer for 5 min, washed with 1 mL of 1 Perm/Wash Buffer, and resuspended in 100 mL of a 30% 1 mg/mL DNase and 70% 1 DPBS solution all at room temperature. The cells are incubated for one hour at 378C. The cells were labeled for BrdU uptake after they are washed in 1 mL of Perm/Wash Buffer and resuspended in a solution containing 1 mL of FITC and 50 mL of 1 Perm/Wash buffer for 20 min. The DNA was then labeled by washing cells in 1 mL of 1
Perm/Wash buffer, and resuspending in 1 mL of propidium iodide (PI) solution (30 mL of 0.1% (v/v) Triton X-100 (Sigma 93443) in PBS, 6 mg DNAse-free Rnase A (Sigma R5503) and 600 mL of a 1 mg/mL PI stock solution (Sigma P4864)). Samples were immediately taken to a BD-Elite flow cytometer where 10,000 cells were analyzed per sample. An argon laser was used for excitation at 488 nm with the FITC emission collected at 525 nm and the PI emission at 675 nm. There was very little emission overlap so no compensation was required. The labeled and unlabeled cells were gated using FITC emission and DNA histograms of both populations were measured with the PI. This pulse labeling BrdU experiment was performed in duplicate and the two cultures, derived the same mother culture, gave extremely consistent results at all time points.
Sample Analysis The percentage of cells in each phase for labeled and unlabeled cells was extracted from the cytometry data by differentiating phases based on the DNA content. The bivariate BrdU incorporation against DNA content density plot were then used to gate the BrdU labeled and unlabeled cells (Fig. 2a). The phase percentages were determined by plotting the DNA histogram of each of these subpopulations (see Fig. 2c and d) and using the ModFit DNA labeling cell cycle analysis software package to define the phase boundaries.
Results and Discussion Extraction of Phase Transition Rates From Jurkat Cells Using BrdU Methods Each transition rate function was fit to experimental BrdU total number density data in each labeled and unlabeled cell cycle phase using the Matlab multivariate mini-
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
965
mization function ‘‘fminsearch.’’ This function employs a Nelder–Mead Simplex method to minimize the chosen least-squares objective function. The least-squares error, E, is the summation over all data points collected—Nsamples collected (13 samples) and Nphases phases in the independent subpopulations (three unlabeled þ three labeled phases)— of the squared difference between the experimental cell cycle phase percentages, xexp i;j , and the model predicted percentages, xsim , that result from the proposed transition rates i;j E¼
NX phases NX samples i¼1
exp sim 2 ðxi;j xi;j Þ
(12)
j¼1
For each independent subpopulation, the experimental phase percentage in the ith phase and jth time sample (where the jth sample is collected at time tj) is the ratio of the number of cells in the phaseP to the total number of cells in the exp exp subpopulation: xexp ¼ N = i;j i;j k Nk;j . Similarly, the simulation phase percentages are X sim xsim N sim ðtj Þ i;j ¼ Ni ðtj Þ= k k where Nisim ðtj Þ
¼
Z
1
ni ðt; tj Þ dt:
0
In most cases, hypotheses and models are written regarding the transition age probability distribution, hi(t), rather than the transition rate, Gi(t), itself. The likelihood of changing cell cycle phases within a certain age window is represented by that phase’s transition-age probability distribution. While this distribution may be observable experimentally, a kinetic model, such as the population balance models, use rate terms such as the rate of transition for a cell of a given age. The age-dependent transition rates can be calculated from the transition-age distributions (Eakman et al., 1966) Gi ðtÞ ¼ such that
h ðtÞ Rti 1 0 hi ðt0Þ dt0 2
hi ðtÞ ¼ Gi ðtÞ41
Zt
The flow cytometry data was analyzed as described in the Sample Analysis Section. The results for pulse-labeling show that the BrdU labeled cells are initially synchronized in the S phase (Fig. 3a). As time progresses, this plug moves to the labeled G2/M phase and then divides into the labeled G0/G1 phase (Fig. 3b) before reentering the labeled S phase (Fig. 3c). The unlabeled cells are also initially synchronized due to labeling of the S phase cells, but to a lesser degree than the unlabeled population due to the addition of the cells not exposed to BrdU. The rate of loss of the synchrony in the labeled and unlabeled subpopulations is dependent upon the shape of the cell cycle phase transition rates.
ð13Þ Determination of Unrestricted Age-Dependent Transition Rates Using BrdU Pulse-Labeling
3 hi ðt 0 Þ dt 0 5;
hi ð0Þ ¼ 0
0
In this study, we evaluated and compared the transition rates derived from normalized Gaussian, lognormal, and skewed lognormal transition rates, along with a more flexible transition rate representation utilizing series of cubic splines with movable spline nodes locations. Ages up to 50 h were evaluated in 0.25 h increments. Polynomials were used to approximate all transition rates with 71 orders giving an excellent fit over the entire age ranges for the transition rates investigated.
966
Figure 3. A series of bivariate BrdU versus DNA cell density plots illustrate the progression of the labeled and unlabeled subpopulations progression around the cell cycle. Ten thousand cells were collected for each figure with speckled or empty areas representing lower densities. The line drawn through each of the plots separates the BrdU labeled population from the unlabeled population. a: Immediately after BrdU inoculation (t ¼ 1 h). b: Nine hours after BrdU inoculation. c: Twenty-four hours after BrdU inoculation.
To determine the general shape of the cell cycle phase transition rates, a flexible representation was first implemented using a series of 15 cubic splines for each of the three cell cycle phase transition rates. The transition rates were manipulated to best fit the average of the two BrdU pulselabel experiments. (Again, these two experiments were conducted simultaneously and cells from each experiment share the same mother culture.) The spline representation was flexible in that both the value of the spline (transition rate) as well as its location (age) could be manipulated to reduce errors due to improper spline node placement. The parameters associated with the age-dependent transition rates were identified by fitting the cell cycle data; the fitted
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
minimal. While the age of cells transitioning out of the S phase or G2/M phase is relatively consistent, the transition age probability distribution for the G0/G1 phase is actually bimodal. To assess the significance the bimodal phenomenon and ensure it is not an artifact of (i) over-fitting with the splines, (ii) the pulse data itself, or (iii) the balanced growth restrictions, the effect of each of these factors was evaluated. The importance of bimodal distributions is compared with single mode distributions in the Determination of Alternative Functional Forms for the Age-Dependent Transition Rates Section and the predictions based on pulse-data experiments are compared with continuous labeling data in the Comparisons With Continuous Labeling Data Section. The implication of the assumption of balanced growth was examined by including the initial phase percentages of each subpopulation as fitted parameters. (The initial phase percentages need not sum to 1 in this case.) As seen in Figures 4 and 5a and b, neither the cell cycle phase transition age probability distribution nor the fits to the experimental data were altered significantly. Determination of Alternative Functional Forms for the Age-Dependent Transition Rates
Figure 4. Best fit of cell cycle phase transition rates to BrdU pulse-labeling data using a 15 node cubic spline representation for the transition rate out of each phase. Both the rate and age locations of each node are variable in the minimization. a: Cell cycle phase transition rates—spline nodes and the polynomial approximation of the cubic spline fit. b: Cell cycle phase transition age densities resulting from those transition rates. cell cycle phase age-dependent transition rate are shown in Figure 4a and the resulting transition age distributions in Figure 4b. The fit of the thick lines shown in Figure 5a and b to the experimental pulse label BrdU data points illustrates the accuracy of the fitted model predictions. Each of the phase transition-age distributions has its own unique shape and characteristics. Both the S and G2/M phase transition-age probability distributions display a positive bias but with relatively small standard deviations (Table II) though the S phase transition is later and slightly broader than the G2/M phase. This would imply that for the S and G2/M phases: (i) cells require at least a certain amount of time to progress through the phase, (ii) cells move through the phase at nearly identical rates, and (iii) delays due to checkpoints or heterogeneity are either unlikely or
While the flexible spline representation of the phase transition rates does an excellent job of matching the pulse-data, the high number of degrees of freedom may lead to over-fitting. That is, simpler or regularized representations may fit nearly as well; in particular, the two transition age modes seen the in the G0/G1 phase may not be crucial to explaining the data. To this end, other basis functions for the age-dependent transition rates, derived from transition age densities of the lognormal, skewed lognormal, and normalized Gaussian distributions were explored. Since Lagrangian information was not used, the only quantitative method to determine the likelihood of a transition age distributions is the residual fit to the experimental data. A comparison of the least-squares errors, provided in Table I, indicates that the flexible cubic splines functional form does fit the experimental data better than the alternative individual functional forms for determining the age-dependent transition rates. Table II lists the best fit parameters for the alternative functional forms. These differences can be seen directly when the model predictions using the spline and lognormal derived transition rates are compared with the experimental phase percentages (see Fig. 5c and d). Particularly for the BrdU labeled subpopulation, the errors in the single lognormal predictions become progressively larger as time progresses while the spline fit remains accurate. In light of the extremely tight confidence intervals on the experimental data, single distributions are not as adequate as the spline fit. When normalized Gaussian distributions are used for the transition probability with respect to age, the G0/G1 representation (Fig. 6a) becomes distorted into a physiologically unlikely arrangement. The estimated G0/G1 normal distribution is positioned such that most of the young cells
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
967
Figure 5. Comparison of BrdU pulse-labeling subpopulation cell cycle phase experimental data with best-fit flexible spline model predictions. The experimental data in each phase is represented by: *, G0/G1; , S; , G2/M. The thick lines correspond to model predictions when the cell cycle phase transition rates are fit to only the pulse data and the initial phase percentages are defined by the transition rates. The thin lines are best-fit model predictions when (a and b) both pulse and continuous data (Comparisons With Continuous Labeling Data Section) is fit and the initial phase percentages are minimization input parameters or (c and d) single lognormals are fit to pulse data (Determination of Alternative Functional Forms for the Age-Dependent Transition Rates Section). a and c: Pulse-labeling dynamics for the BrdU labeled subpopulation. b and d: Pulse-labeling dynamics for the BrdU unlabeled subpopulation.
Table I. Least-squares fit, E, when the minimization is performed for each distributions examined along with the number of parameters manipulated to attain that fit. Fit to pulse data only: restricted initial condition Flexible cubic splines (defines the transition rates) Normal Lognormal Skewed lognormal 2 lognormals/phase 3 lognormals/phase
968
Total number of E ( 103), fitted pulse parameters data
E ( 103), continuous data
90
1.9553
1.0004
6 6 9 15 24
3.2681 2.5217 2.4832 1.6891 1.6893
1.3472 1.1682 1.1908 1.3655 1.3146
just entering the G0/G1 phase transition quickly while the large variance still enables the transition of older cells: that is, the increased variance leads to a likely over-prediction of younger cell division rates, even though the residual’s error when fit to the data is relatively low. When a lognormal is used as the transition-age distribution (Fig. 6b), the G0/G1 phase distribution is again extremely broad. Skewed lognormal distributions were intuitively appealing because they can enhance the bias of distributions, but this bias is of little use in capturing the bimodal behavior. The skewed lognormal fit is similar to that of the lognormal (Fig. 6c). No single distribution per phase can capture the cubic-spline estimated bimodal behavior of the G0/G1 phase, but rather,
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
Table II.
Summary of the parameter sets attained from the best fits of each of the phase transition representations.
Distribution: restricted initial condition Cubic spline (calculated)
Single transition age distributions Normal
Lognormal
Skewed lognormal
Lognormal basis functions’ transition age distributions 2 lognormal combination
3 lognormal combination
Distribution: unrestricted initial condition 3 lognormal combination ( f1,0þ ¼ [0.00 106.33 0.23], f1,0 ¼ [54.21 18.51 18.07], f1,0 ¼ [54.94 32.79 10.40])
Phase
Mean, m (h)
Standard deviation, s (h)
Weight, w, or asymmetry
G0/G1 S G2/M
16.2627 11.2515 6.3182
9.3577 2.5025 1.7564
G0/G1 S G2/M G0/G1 S G2/M G0/G1 S G2/M
11.31 10.50 5.79 16.59 10.49 5.76 16.70 10.04 5.64
13.39 0.70 0.51 14.75 1.04 0.80 14.15 1.07 0.82
1.00 1.08 0.94
G0/G1 S G2/M G0/G1 S G2/M
½ 7:57 25:61 ½ 10:52 15:22 ½ 6:09 6:26 ½ 7:75 25:39 7:93 ½ 10:47 14:45 10:53 ½ 6:03 5:28 8:77
½ 1:00 2:21 ½ 1:00 1:21 ½ 0:99 1:13 ½ 1:15 1:62 1:29 ½ 1:23 1:08 1:27 ½ 0:91 1:00 1:39
½ 0:55 0:45 ½ 0:92 0:08 ½ 0:59 0:41 ½ 0:50 0:45 0:05 ½ 0:87 0:09 0:04 ½ 0:99 0:01 0:00
G0/G1 S G2/M
½ 25:13 8:79 5:28 ½ 11:22 10:73 5:25 ½ 6:75 6:09 5:40
½ 1:87 1:39 0:99 ½ 1:24 1:65 1:28 ½ 1:30 1:14 0:89
½ 0:36 0:59 0:05 ½ 0:05 0:95 0:00 ½ 0:01 0:25 0:74
Empty items imply that that parameter does not influence that distribution.
this results in a more disperse distribution to accommodate the data as indicated by the larger standard deviations associated with the G0/G1 phase in Table II. Regardless of these differences in the G0/G1 phase, the means and variations seen in the S and G2/M phases for each of the single mode transition age distributions, are similar to those of the flexible splines (Fig. 4 vs. Fig. 6; Table II). What effect does a broad G0/G1 phase versus a bimodal G0/G1 phase have on ability of the model to fit the BrdU pulse data? Comparing the best predictions for the BrdU labeled dynamics, the broad G0/G1 phase does produce oscillations of slightly lower amplitude than the bimodal fit, but this looks to be of marginal significance in relating to the experimental BrdU pulse data (Fig. 5c). In addition, the predictions for the BrdU unlabeled subpopulation using the broad or bimodal G0/G1 phase are nearly indistinguishable (Fig. 5d). Lognormal Basis Functions The predictions of the single mode distributions were improved upon by using a weighted combination of either two or three lognormal distributions to describe the transition age distribution in each phase where the means, standard deviation, and weights are to be determined hlognormal ðtÞ ¼ i w1 hlognormal þ w2 hlognormal i;1 i;2 X þ w3 hlognormal ; wj ¼ 1 i;3 j
This description invites the possibility of multiple subpopulations or more distributed transition-age distributions using fewer parameters. We found that the combinations of lognormals orient themselves to form transition agedistributions that are similar in shape to the flexible spline nodes representation (Fig. 4b vs. Fig. 7a and b). In the S and G2/M phases, the two lognormals in each phase have considerable overlap and the superposition forms distributions with means and standard deviations that are slightly greater those derived using the splines (see Table II). However, there were significant differences in the means of the two lognormals that comprise the G0/G1 phase which results in a transition age probability distribution that is very similar to that found with the splines. Using three lognormal distribution per phase results in predicted labeled and unlabeled phase dynamics are that are similar between when cubic-splines (Fig. 4b) and two lognormals per phase (Fig. 7a) though the number parameters used was greater than five times less than that required for the spline fit.
Comparisons With Continuous Labeling Data For validation of this method, data was obtained from a continuous BrdU label where all cells eventually become BrdU labeled during progression through the S phase. This data explicitly shows the decay in the G0/G1 and G2/M phases. All cell preparation and assay experimental protocols are the same as described earlier for the pulse labeled BrdU experiments with the exception of BrdU labeling of the entire 20 mL flask of cells. (The flask contents is transferred
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
969
Figure 6.
Best fit of cell cycle transition rates to BrdU dynamic data from single transition-age distributions for each phase. a: Gaussian cell cycle phase transition-age distributions. b: Lognormal cell cycle phase transition-age distributions. c: Skewed lognormal cell cycle phase transition-age distributions.
to a centrifuge tube and exposed to 200 mL of the diluted BrdU solution; no cells are left unexposed.) As before, these continuous BrdU labeling experiments were performed in duplicate with extremely consistent results at all time points. During continuous labeling, unlabeled cells entering the S phase become labeled with BrdU as shown in the bivariate BrdU–DNA plots time sequence of Figure 8. The unlabeled population is initially in unbalanced growth. This data is used as an independent verification for the G0/G1 and G2/M predictions based on the pulse data. In light of the continuous labeling observed in the experimental data, Equation (9), which describes the phase number densities of the labeled and unlabeled cells, must be modified to reflect the accumulation of labeled cells. Equation (14) accounts for
970
this accumulation and can be used in lieu of Equation (9) for continuous labeling Labeled ðBrdU þ Þ Subpopulation 1 P ðkÞ n1 ðt; tÞ ¼ n1 ðt; tÞ n2 ðt; tÞ ¼ n3 ðt; tÞ ¼
k¼2 1 P k¼0 1 P k¼1
ðkÞ
n2 ðt; tÞ (14)
ðkÞ
n3 ðt; tÞ
Unlabeled ðBrdU Þ Subpopulation ð0Þ ð1Þ n1 ðt; tÞ ¼ n1 ðt; tÞ þ n1 ðt; tÞ n2 ðt; tÞ ¼ 0 ð0Þ n3 ðt; tÞ ¼ n3 ðt; tÞ
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
Figure 8. A series of bivariate BrdU versus DNA cell density plots illustrate the progression of the labeled and unlabeled subpopulations progression around the cell cycle. Ten thousand cells were collected for each figure with speckled or empty regions representing lower cell densities. a: Immediately after BrdU inoculation (t ¼ 1 h). b: Ten hours after BrdU inoculation. The arrow denotes movement of BrdU unlabeled cells becoming labeling during the S phase. c: Twenty-four hours after BrdU inoculation.
Figure 7. Linear combinations of lognormal distributions used to approximate the residence time distribution of each cell cycle phase. The added lognormal distributions give additional flexibility to the G0/G1, S and G2/M distributions. a: Two lognormal distributions per phase. b: Three lognormals per phase.
of balanced growth using all available data, 33 parameters are fitted to three lognormals per phase. The identified transition rates are extremely consistent (Fig. 4b) from when the pulse data is used alone, and the fit to the pulse data is also similar (see Fig. 5a and b) with a least squares fit error to the pulse data of 1.1972 103. The fit to the continuous BrdU data is also improved (Fig. 9c and d) with a least squares fit error of 1.08163 103.
Conclusions The pulse data derived transition rates for both the flexible spline and lognormal distributions were used to predict the dynamics of the continuous BrdU. The least squares error quantifying the predicted model dynamics fit to the continuous experimental data are provided in Table I for all assumed distribution forms. As illustrated in Figure 9 the spline and lognormal distribution forms both capture the trends of the continuous labeling experimental data. Throughout the continuous BrdU labeling experiment, the ratio of labeled cell fluxes into and out of the S phase remains constant hence, it impossible to determine the age structure of the S phase transition rates using this data alone. When the pulse-data (two independent subpopulations) is supplemented with the continuous data (one subpopulation) for fitting spline transition rates without the restriction
A methodology has been developed that enables predictions of cell cycle phase age-dependent transition rates based on experimental data with BrdU labeling where the crux of the identification relies on comparison of the BrdU data with model predictions based entirely on the transition rates. For making model predictions based on transition rates, an agestructured cell cycle model is presented in which the agedependent transition rates are posed. Model predictions require knowledge of initial age-distributions, which are unidentifiable directly, so an expression that yields initial age-distributions from transition rates is derived under the assumption of a balanced growth setting. Since the transition rates and an initial condition are sufficient for model identification, and the transition rates imply the initial condition, the age-dependent transition rates are the only
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
971
Figure 9. Comparison of BrdU continuous labeling subpopulation cell cycle phase experimental data with best-fit flexible spline model predictions. The experimental data in each phase is represented by: *, G0/G1; , S; , G2/M. The thick lines correspond to model predictions when the cell cycle phase transition rates are fit to only the pulse data and the initial phase percentages are defined by the transition rates. The thin lines are model predictions when (a and b) single lognormals are fit to pulse data while in (c and d) both pulse and continuous data is fit and the initial phase percentages are minimization input parameters. a and c: Continuous labeling dynamics for the labeled population. b and d: Continuous labeling dynamics for the unlabeled population.
input for making model predictions. We then illustrate how pulse BrdU labeling, which introduces instances of unbalanced growth while maintaining individual cells and the total population in balanced growth, is sufficient for the identification of the unique set of transition rates from the unbalanced growth dynamics. The identified age transition distributions were validated through a comparison of the predicted dynamics with those observed for continuous BrdU labeling. In identifying the transition rates, a general representation of the transition rates using a series of cubic splines, where both the spline node values as well as locations were minimization parameters, resulted in transition probability distributions with both expected and unexpected features. Very few transitions occur before a threshold age and after
972
which there is a window of ages when transitions are very likely; however, for the G0/G1 phase, the best fit transition age probability distribution suggested multiple subpopulations through its bimodal form. This bimodal behavior may indicate the presence of multiple subpopulations resulting from either continuous, systematic differences between cells or individual, discrete events. With multiple distinct modes, each subpopulation is governed by an independent set of cell cycle kinetics. Or, the subpopulations could result from a random event which may or may not occur each time through the cycle such as a checkpoint control. While the existence of multiple subpopulations is able to be identified and their behavior characterized using age-dependent transition rates, a Lagrangian approach where individual cells are followed would be required to distinguish between
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
the two mechanisms. Fit of the age probability distribution to the data assuming single normalized Gaussian, lognormal, and skew-lognormal transition probability distributions was reasonable, but not as accurate due to the inability to capture the G0/G1 bimodal transition age distribution. But, describing each phase transition with a convex combination of two or three lognormal distributions improved the overall fit of the BrdU data relative to the spline representation while requiring only 15 or 24 parameters rather than the 90 used for the spline fit. We see this approach as an efficient alternative with adequate flexibility to capture the major features of data while minimizing the number of parameters to be identified. The foregoing inversion involved no regularization as the error bounds on the data were very small. However, larger error bounds would call for the use of appropriate regularization (Tikhonov, 1977).
R1
0 fi ðtÞ dt ¼ 1) will give the rates of change in the total number density Z 1 Z 1 dNi ðtÞ dfi ðtÞ þ Ni ðtÞ dt ¼ Ni ðtÞ Dout i;i ðtÞfi ðtÞ dt dt dt 0 0
The boundary and containment conditions can be used to simplify the term Z 1 XZ 1 dfi ðtÞ dt ¼ fi ðtÞj1 fi ðtÞj0 ¼ Din i;j fj ðtÞ dt dt 0 0 j to give Equation (3) dNðtÞ ¼ DNðtÞ dt where Dij ¼
The authors are grateful for the assistance of the Purdue University Cytometry Laboratories and support from the Purdue Cancer Center.
Derivation of Initial Age-Distributions (Age-Distribution During Balanced Growth) Balanced growth was earlier defined as the growth condition that satisfies two criterions: (i) intransient growth conditions and (ii) time-invariant cell-state distributions. Timeindependent transition rates satisfy the first condition, and the second condition implies that the age-distribution in the ith cell cycle phase is described by the time-independent distribution fi ðtÞ when balanced growth is achieved. So, during balanced growth, the number density in the ith phase is ni ðt; tÞ ¼ Ni ðtÞfi ðtÞ (Eq. (2)). Substitution of this balanced growth number density into the population balance equation for the ith phase (Eq. (1)) @ni ðt; tÞ @ni ðt; tÞ þ ¼ Dout i;i ðtÞni ðt; tÞ @t @t
out ½Din ij ðtÞ Dij ðtÞ fj ðtÞ dt
0
dNi ðtÞ ¼ mNi ðtÞ dt The use of this expression in Equation (15) allows the total number density to be removed, leaving a differential equation (Eq. 4) for change in the age-distribution with age 3 P i ðtÞ fi ðtÞmNi ðtÞ þ Ni ðtÞ dfdt ¼ Dout ij ðtÞNj ðtÞfj ðtÞ j¼1 " # 3 P dfi ðtÞ out ðDij ðtÞfj ðtÞÞ þ mfi ðtÞ dt ¼ j¼1
dfi ðtÞ dt
¼ ðDout ii ðtÞ þ mÞfi ðtÞ
The general solution of this type of differential equation is
Z t 0 0 ðDout ðt Þ þ mÞ dt fi ðtÞ ¼ fi ð0Þ exp ii
gives @ðNi ðtÞfi ðtÞÞ @ðNi ðtÞfi ðtÞÞ þ ¼ Dout i;i ðtÞNi ðtÞfi ðtÞ @t @t
0
Since the total product density is age-independent and the balanced growth age-distribution is time-independent, this equation simplifies to dNi ðtÞ dfi ðtÞ þ Ni ðtÞ ¼ Dout i;i ðtÞNi ðtÞ fi ðtÞ dt dt
1
The time independent age-distribution gives a timeindependent coefficient matrix D. The specific growth rate of each cell cycle phase during balanced growth will be the positive real eigenvalue of this matrix, m, so
Appendix
fi ðtÞ
Z
where the regulatory condition
0
fi ðtÞ ¼ 0 defines fi ð0Þ
Rt R1 0 0 fi ðtÞ dt ¼ 1 ¼ fi ð0Þ 0 exp 0 ðDout ii ðt Þ þ mÞ dt dt i ) fi ð0Þ ¼ R 1 h R 1 1 R1 0
0
(15)
A further implication of the time-independent agedistributions is identical specific growth rates for each cell cycle phase during balanced growth. Returning to the population balance equation, substitution of the balanced growth separated number density ni ðt; tÞ ¼ Ni ðtÞfi ðtÞ into Equation (1) and integrating over all ages (where
R1
exp
ðDout ii ðt0ÞþmÞ dt0 dt
0
So, the age-distribution in the ith cell cycle phase during balance growth is given by Equation (5) fi ðtÞ ¼
Z
Z t0 1 00 00 exp ðDout ðt Þ þ mÞ dt dt 0 ii 0
Z t 0 0 0
exp ðDout ðt Þ þ mÞ dt ii 1
0
Sherer et al. : Cell Cycle Phase Transitions Biotechnology and Bioengineering. DOI 10.1002/bit
973
References Arino O, Kimmel M. 1987. Asymptotic analysis of a cell cycle model based on unequal division. SIAM J Appl Math 47:128–145. Baisch H, Otto U. 1993. Intratumoral heterogeneity of S phase transition in solid tumours determined by bromodeoxyuridine labelling and flow cytometry. Cell Prolif 26:439–448. Baisch H, Otto U, Hatje U, Fack H. 1995. Heterogeneous cell kinetics in tumors analyzed with a simulation model for bromodeoxyuridine single and multiple labeling. Cytometry 21:52–61. Bernard S, Pujo-Menjouet L, Mackey MC. 2003. Analysis of cell kinetics using a cell division marker: Mathematical modeling of experimental data. Biophys J 84:3414–3424. Chiorino G, Metz JA, Tomasoni D, Ubezio P. 2001. Desynchronization rate in cell populations: Mathematical modeling and experimental data. J Theor Biol 208:185–199. Cooper S. 2003. Rethinking synchronization of mammalian cells for cell cycle analysis. Cell Mol Life Sci 60:1–8. Cooper S. 2004. Is whole-culture synchronization biology’s ‘‘perpetualmotion machine’’? Trends Biotechnol 22:266–269. Darzynkiewicz Z, Gong J, Juan G, Ardelt B, Traganos F. 1996. Cytometry of cyclin proteins. Cytometry 25:1–13. Dien BS, Srienc F. 1991. Bromodeoxyuridine labeling and flow cytometric identification of replicating Saccharomyces cerevisiae: Lengths of cell cycle phases and population variability at specific cell cycle positions. Biotechnol Prog 7:291–298. Del Bino G, Lassota P, Darzynkiewicz Z. 1991. The S-phase cytotoxicity of camptothecin. Exp Cell Res 193:27–35. Dolbeare F, Gratzner H, Pallavicini MG, Gray JW. 1983. Flow cytometric measurement of total DNA content and incorporated bromodeoxyuridine. Proc Natl Acad Sci USA 80:5573–5577. Dyson J, Villella-Bressan R, Webb G. 2000. A nonlinear age and maturity structured model of population dynamics. I. Basic Theory. J Math Anal Appl 242:93–104. Dyson J, Villella-Bressan R, Webb GF. 2002. Asynchronous exponential growth in an age structured population of proliferating and quiescent cells. Math Biosci 177–178:73–83. Eakman JM, Fredrickson AG, Tsuchiya HM. 1966. Statistics and dynamics of microbial cell populations. Chem Eng Prog Symp Ser 62:37–49. Fredrickson AG, Mantzaris NV. 2002. A new set of population balance equations for microbial and cell cultures. Chem Eng Sci 57:2265–2278. Fredrickson AG, Tsuchiya HM. 1963. Continuous propagation of microorganisms. AIChE J 9:459–468. Fredrickson AG, Ramkrishna D, Tsuchiya HM. 1967. Statistics and dynamics of procaryotic cell populations. Math Biosci 1:327–374. Hjortsø MA. 2006. Population balances in biomedical engineering: Segregation through the distribution of cell states. New York: McGraw-Hill. 182 p. Ja´ nossy J, Ubezio P, Apa´ ti A´ , Mago´ csi M, Tompa P, Friedrich P., 2004. Calpain as a multi-site regulator of cell cycle. Biochem Pharmacol 67:1513–1521. Kromenaker SJ, Srienc F. 1991. Cell-cycle-dependent protein accumulation by producer and nonproducer murine hybridoma cell lines: A population analysis. Biotechnol Bioeng 38:665–677. Kutalik Z, Razaz M, Elfwing A, Ballagi A, Baranyi J. 2005. Stochastic modelling of individual cell growth using flow chamber microscopy images. Int J Food Microbiol 105:177–190. Liou JJ, Srienc F, Fredrickson AG. 1997. Solutions of population balance models based on a successive generations approach. Chem Eng Sci 52:1529–1540.
974
Liou JJ, Fredrickson AG, Srienc F. 1998. Selective synchronization of tetrahymena pyriformis cell populations and cell growth kinetics during the cell cycle. Biotechnol Prog 14:450–456. Liu SV. 2005. Debating cell-synchronization methodologies: Further points and alternative answers. Trends Biotechnol 23:9–10. Lyons AB. 1999. Tracking cell proliferation with carboxyfluorescein diacetate succinimidyl ester. Immunol Cell Biol 77:509–515. Mantzaris NV, Liou JJ, Daoutidis P, Srienc F. 1999. Numerical solution of a mass structured cell population balance model in an environment of changing substrate concentration. J Biotechnol 71:157– 174. Mantzaris NV, Srienc F, Daoutidis P. 2002. Nonlinear productivity control using a multi-staged cell population balance model. Chem Eng Sci 57:1–14. Nordon RE, Nakamura M, Ramirez C, Odell R. 1999. Analysis of growth kinetics by division tracking. Immunol Cell Biol 77:523– 529. Nova´ k B, Tyson J. 1993. Modeling the cell division cycle: M-phase trigger, oscillations, and size control. J Theor Biol 165:101–134. Nova´ k B, Tyson J. 2004. A model for restriction point control of the mammalian cell cycle. J Theor Biol 230:563–579. Ormerod MG, Kubbies M. 1992. Cell cycle analysis of asynchronous cell populations by flow cytometry using bromodeoxyuridine label and Hoechst-propidium iodide stain. Cytometry 13:678–685. Ramkrishna D. 2000. Population balances: Theory and applications to particulate systems in engineering. San Diego: Academic Press. 355 p. Sena G, Onado C, Cappella P, Montalenti F, Ubezio P. 1999. Measuring the complexity of cell cycle arrest and killing of drugs: Kinetics of phasespecific effects induced by taxol. Cytometry 37:113–124. Sennerstam R, Stro¨ mberg , JO. 1995. Cell cycle progression: Computer simulation of uncoupled subcycles of DNA replication and cell growth. J Theor Biol 175:177–189. Shackney SE. 1973. A cytokinetic model for heterogeneous mammalian cell populations. J Theor Biol 38:305–333. Shackney SE, Shankey TV. 1999. Cell cycle models for molecular biology and molecular oncology: Exploring new dimensions. Cytometry 35:97– 116. Sherer E, Hannemann RE, Rundell A, Ramkrishna D. 2006. Analysis of resonance chemotherapy in leukemia treatment via multi-staged population balance models. J Theor Biol 240:648–661. Steel GG. 1977. Growth kinetics of tumours: Cell population kinetics in relation to the growth and treatment of cancer. Oxford: Clarendon Press. 351 p. Subramanian G, Ramkrishna D. 1971. On the solution of statistical models of cell populations. Math Biosci 10:1–23. Subramanian G, Ramkrishna D, Fredrickson AG, Tsuchiya HM. 1970. On the mass distribution model for microbial cell populations. Bull Math Biophys 32:521–537. Tikhonov AN. 1977. Solution of ill-posed problems. Washington: Winston & Sons. 272 p. Tyrcha J. 2001. Age-dependent cell cycle models. J Theor Biol 213:89– 101. Tyrcha J. 2004. Cell cycle progression. C.R. Biol 327:193–200. Webb GF, Grabosch A. 1987. Asynchronous exponential-growth in transition-probability models of the cell-cycle. SIAM J Math Anal 18:897– 908. Zamamiri AM, Zhang Y, Henson MA, Hjortsø MA. 2002. Dynamic analysis of an age-distribution model of oscillating yeast cultures. Chem Eng Sci 57:2169–2181.
Biotechnology and Bioengineering, Vol. 99, No. 4, March 1, 2008 DOI 10.1002/bit
