A systems biology dynamical model of mammalian G1 cell cycle progression.

The current dogma of G(1) cell-cycle progression relies on growth factor-induced increase of cyclin D:Cdk4/6 complex activity to partially inactivate pRb by phosphorylation and to sequester p27(Kip1)-triggering activation of cyclin E:Cdk2 complexes that further inactivate pRb. pRb oscillates between an active, hypophosphorylated form associated with E2F transcription factors in early G(1) phase and an inactive, hyperphosphorylated form in late G(1), S and G(2)/M phases. However, under constant growth factor stimulation, cells show constitutively active cyclin D:Cdk4/6 throughout the cell cycle and thereby exclude cyclin D:Cdk4/6 inactivation of pRb. To address this paradox, we developed a mathematical model of G(1) progression using physiological expression and activity profiles from synchronized cells exposed to constant growth factors and included a metabolically responsive, activating modifier of cyclin E:Cdk2. Our mathematical model accurately simulates G(1) progression, recapitulates observations from targeted gene deletion studies and serves as a foundation for development of therapeutics targeting G(1) cell-cycle progression.

Introduction

Growth factor stimulation of G0 quiescent cells drives them into the growth factor-dependent early G1 phase of the cell cycle, followed by transition into the growth factor-independent late G1 phase and progression into S phase (Ho and Dowdy, 2002; Ortega ; Sherr, 2004). This irreversible commitment point, termed the restriction point (Pardee, 1974), is correlated with inactivation of the retinoblastoma tumor suppressor protein (pRb) by hyperphosphorylation. pRb, a negative regulator of G1-phase cell-cycle progression, binds E2F transcription factors and chromatin remodeling proteins (HDAC, Suv39, BRG1) to repress E2F target gene expression (Ho and Dowdy, 2002; Ortega ; Sherr, 2004). pRb is unphosphorylated and not bound to E2Fs in G0 quiescent cells. Growth factor stimulation of cells and transition into early G1 results in hypophosphorylation of pRb by cyclin:Cdk complexes that promote its assembly with E2Fs (Ezhevsky ). At the restriction point, cyclin:Cdk complexes inactivate pRb by hyperphosphorylation resulting in E2F release and target gene induction. In both cycling normal and tumor cells, pRb oscillates between an active, hypophosphorylated form in early G1 and an inactive, hyperphosphorylated form in late G1, S and G2/M phases (Fang ; Ezhevsky , 2001; Brugarolas ; Nagahara ; Stevaux and Dyson, 2002). Thus, pRb is differentially regulated by two opposing cyclin:Cdk complexes.

The current paradigm of G1 cell-cycle progression argues that growth factor stimulation of quiescent cells results in the gradual accumulation of active cyclin D:Cdk4/6 complexes that perform at least two critical functions for cell-cycle progression: (1) partial inactivation of pRb at the restriction point, resulting in release of E2F transcription factors that drive cyclin E expression, and (2) sequestration of the p27Kip1 Cdk inhibitor (CKI), resulting in activation of cyclin E:Cdk2 complexes that complete pRb inactivation. However, recent targeted gene deletions of G1-phase cyclins and Cdks in mice have raised several concerns about the validity of this paradigm (Pagano and Jackson, 2004; Sherr and Roberts, 2004). Moreover, similar to tumor cells growing in vivo, the current paradigm remains largely untested in proliferating cells exposed to continuous growth factors. Here, using experimental data derived from highly synchronized cells exposed to continuous growth factors, we developed a mathematical model of G1 cell-cycle progression that accurately simulates the observed physiological data, kinetics and transitions. The model points to the presence of an unknown activating ‘modifier' element that responds to metabolic input to activate Cdk2 and potentially serves as a functional analogue of the yeast G1-phase activator Bck2 (Newcomb ; Costanzo ; de Bruin ).

Results and discussion

To dissect regulation of G1 cell-cycle progression, we analyzed human HCT116 colon carcinoma cells (RB, p53, p16) grown under conditions of constant growth factor exposure. HCT116 cells were synchronized in early G1 phase by density arrest (contact inhibition) for 48 h in the presence of serum and released from arrest by low-density replating. Synchronized cells consistently entered S phase at 10 h and G2/M phase by 18 h (Figure 1A). This highly reproducible cell-cycle synchronization method served as the basis for all subsequent experiments.

Figure 1

Physiological expression and activity profiles of cell-cycle regulators during transition from early G1, across the restriction point, into late G1 and S phase. Density-arrested (contact-inhibited) and released HCT116 cells were assayed for (A) kinase activity and immunoblot analysis, (B) mRNA levels by semi-quantitative RT–PCR, (C) cyclin E:Cdk2 association by co-immunoprecipitation–immunoblot analysis, (D) released cells were plated into media containing glutamine (+Gln) or deprived of glutamine (−Gln) and assayed for cyclin E, Cdk2 and Cdk4 kinase activity and pRb immunoblot analysis; (E) Cdk2 and Cdk4/6 kinase activity and pRb immunoblot following selective Cdk2 inhibition by roscovitine (15 μM) at indicated time points. All studies were performed in triplicate with similar results; arrowhead denotes slower migrating, inactive hyperphosphorylated pRb (panels B, D); panel D is derived from two regions of the same immunoblot.

We examined the expression and activity profiles of cell-cycle regulatory proteins at various time points in the synchronized HCT116 cells, which yielded several important observations: (1) cyclin D1 and D2 were constitutively expressed and cyclin D:Cdk4/6 complexes were constitutively active from start to finish (0–27 h) (Figure 1B). We did not detect cyclin D3. (2) In contrast, cyclin E:Cdk2 complexes were present throughout the cell cycle, but oscillated between inactive complexes in early G1 and active complexes in late G1, correlating with the inactivating, hyperphosphorylation of pRb and p130 (Figure 1B). (3) pRb and p130 pocket proteins were present in their active, hypophosphorylated forms in early G1 phase (0–4 h) and the inactive, hyperphosphorylated forms first appeared at 6 h (Figure 1B). We did not detect p107. (4) Consistent with hyperphosphorylated pRb/p130 being functionally inactive, the switch from hypo- to hyperphosphorylation of pRb/p130 correlated with expression of bona fide E2F target genes, cyclin A2 and DHFR (Figure 1C), in late G1. Similar profiles were obtained from S-phase-synchronized arrested/released HCT116 cells (data not shown). (5) Lastly, similar to G1-induced arrest by metabolic starvation of glucose (Singh ), deprivation of glutamine amino acids results in an early G1 cell-cycle arrest typified by continued cyclin D:Cdk4/6 activity, the presence of active hypophosphorylated pRb and inactive cyclin E:Cdk2 complexes (Figure 1D).

Thus, under constant growth factor conditions, we find constitutive cyclin D:Cdk4/6 activity, whereas cyclin E:Cdk2 and pRb activity oscillates during transition from early to late G1 phase. These observations are consistent with results from cycling human T cells and cycling human primary fibroblasts (Ezhevsky ). These observations combined with the triple cyclin D and double Cdk4/6 genetic ablation studies (Kozar ; Malumbres ) suggest that cyclin D:Cdk4/6 activity is neither necessary nor sufficient for activation of cyclin E:Cdk2 complexes. Thus, the rate-limiting activator of Cdk2 complexes at the restriction point remains undefined.

Several mathematical models have been developed that provide insights into dynamical mechanisms of the cell cycle under specific conditions (Aguda and Tang, 1999; Qu ; Novak and Tyson, 2004). Novak and Tyson's (2004) model of the restriction point reproduces earlier experimental results of cell-cycle progression in cycloheximide-treated cells performed by Zetterberg and Larsson (1995). That model as well as a model developed by Qu reflect the current paradigm, which ultimately requires cyclin E:Cdk2 complexes to oscillate in phase with inactivated pRb. In both models, cyclin E:Cdk2 also acts as part of a negative feedback loop promoting the degradation of its inhibitor p27. Although these models reproduce central features of G1 progression in mammalian cells, they do not account for our experimental results of constitutively active cyclin D and/or constitutive cyclin E expression throughout G1 phase. Modification of the Tyson model to constitutively express cyclin E resulted in continued activation of cyclin E:Cdk2 complexes and constant nonphysiological inactivation of pRb. A similar effect on pRb hyperphosphorylation was observed when p27 was eliminated from this model. In contrast, p27 genetic deletions have no effect on the timing of pRb inactivation (Pagano and Jackson, 2004; Sherr and Roberts, 2004). Therefore, we surmised that the rate-limiting mechanism governing Cdk2 activation requires mechanisms distinct from those of positive feedback loops and antagonistic interactions represented in these mathematical models.

To quantitatively describe experimentally observed dynamics of restriction point progression, we developed a mathematical model of G1 progression that incorporates additional regulation of Cdk2. For ease of interpretation, the mathematical model is presented in diagrammatic cell language (DCL) (Figure 2A) and in systems biology markup language (see Supplementary information). The model diagram includes interactions between some of the relevant players in G1 cell-cycle progression, including cyclins, Cdks, CKIs, pRb, E2F transcription factors, APC/C and Emi1 (Ezhevsky ; Ho and Dowdy, 2002). Cdk4 represents both Cdk4 and Cdk6, as both proteins have been shown to have homologous functions in G1 progression and they are activated when bound to cyclin D (D1, D2 or D3). Cyclin D:Cdk4/6 complexes are inhibited by p16-INK4a family (not present in HCT116 cells) and p27 family of proteins. Cyclin E and cyclin A complexed with Cdk2 are also inhibited by p27 binding.

Figure 2

Dynamical model of the mammalian G1 cell-cycle progression. (A) The mathematical model is represented in the DCL. Cdks and pRb are boxes, with circles inside denoting possible states. As an example we show the various Cdk2 states as vectors of integers (box read in the order of top to bottom and left to right). Proteins with a single internal state are represented by ovals. Transitions between states are represented by filled circles for binding and unbinding, and by squiggles that can also denote degradation of bound components or conformational change. The diagram describes Cdk4 and Cdk2 binding to p27. Cdk2 binds cyclin E or A, Cdk1 binds cyclin A and Cdk2/1 can undergo a conformational change by the activating ‘modifier'. pRb is hypo- to hyperphosphorylated by active forms of Cdk4 (light green) and Cdk2/Cdk1 complexes, respectively (dark green). Hyperphosphorylated pRb unbinds E2F; E2F promotes the synthesis of cyclin E, A, Emi and itself. Cyclin A is degraded by APC/C when not bound to Emi. (B–E) Simulations (solid lines) of the model fit to data from six representative parameter sets that meet the deletion criteria (dots with error bars). Time zero refers to the synchronized entry of the cells into early G1. (B) Active Cdk2, total p27 and total cyclin E. (C) hyperphosphorylated E2F and pRb and (D) total Emi and total cyclin A. (E) In silico prediction of Cdk2/1 inhibition: hyperphosphorylated pRb at 16 h with respect to Cdk2/1 inhibition, normalized to control (0% inhibition) for the population of six models that best fit the data and meet the deletion criteria.

In formulating the model, we made select simplifications to reduce the number of model parameters and focus on modeling the biology relevant for ascertaining the potential rate-limiting mechanism governing Cdk2 activation. Cdk2 and Cdk1 are present in an inhibited phosphorylated (Thr14 and Tyr15) state and require dephosphorylation by Cdc25A, followed by Thr160 or Thr161 activating phosphorylation, respectively, by CAK. We did not explicitly include these kinases and phosphatases in our mathematical model because, beyond the requirement for Cdk activation, their role in regulating G1 transition is not yet well defined. As this model concentrates on G1 progression, interactions involving APC/C were also simplified (Vodermaier and Peters, 2002; Murray, 2004).

Parameter values of processes described in the model are generally unknown in vivo, are cell-type-specific and can vary over large ranges. Physiologically correct model output is fairly robust to variations of parameters within those ranges. Most parameter values are determined by fitting to quantitative data. To add additional constraints to the possible solutions in parameter space that are consistent with the quantitative data, we fixed some parameter values to physiological ranges found in the literature. Protein unbinding rates were set to 0.1 min−1, which is of the same order of magnitude as in other models of protein interaction. Turnover for cyclins is rapid and was fixed to 0.05 min−1 (half-life of 14 min) (Maity ; Carlson ; Kahl and Means, 2004). Concentrations of proteins in ‘ground' states were set to physiologically reasonable rates found in the literature. Cdk levels were 100 000 molecules/cell (Arooz ). Components such as pRb and APC/C have concentrations lower than Cdks.

Cycling cells (normal and tumor) represent an oscillating system in a constant environment that periodically resets itself, requiring doubling of cellular mass and genome before cytokinesis (Jorgensen and Tyers, 2004). Recent evidence suggests that achieving an autonomous critical cell mass is a prerequisite for cells to progress from G1 into S phase and complete cell division (Dolznig ; Jorgensen and Tyers, 2004). Also, treatment of cells with cell-cycle inhibitors does not slow down cellular growth, whereas inhibitors of cellular growth induce cell-cycle arrest, suggesting that growth directly or indirectly regulates cell-cycle progression (Fingar ; Kim ; Dutcher, 2004). In addition, metabolic restriction experiments result in an early G1 cell-cycle arrest containing active cyclin D:Cdk4/6 complexes, active hypophosphorylated pRb and inactive cyclin E:Cdk2 complexes (Figure 1D; Singh ). Moreover, the analogous G1 regulatory pathway in yeast (START) contains a metabolically regulated activator of G1 cell-cycle progression, termed Bck2 (Newcomb ; Costanzo ; de Bruin ). Taken together, these observations led us to introduce into our G1 progression model a ‘cell growth' variable as an unspecified activating ‘modifier' of Cdk2 (Figure 2A).

The activating ‘modifier' is implemented as a switch, turned on at the end of early G1 phase inducing a change in Cdk2 required for its activation, which is the simplest functional interaction. Genetic ablation studies of cyclin E or Cdk2 demonstrate that primary and tumor cells fail to arrest (Berthet ; Ortega ; Tetsu and McCormick, 2003). In contrast to previously published models of G1 progression in mammalian cells, our mathematical model also incorporates mechanisms to account for cyclin E or Cdk2 deletion. Consistent with appropriate kinetic activation of cyclin A:Cdk2/1 activation in cyclin E1/E2-deficient cells (Kozar ), cyclin A:Cdk2 or cyclin A:Cdk1 complexes are also activated by the ‘modifier', providing compensatory inactivating hyperphosphorylation of pRb.

Global parameter optimization techniques yielded an ensemble of parameter sets that enabled our model to quantitatively reproduce the data in Figure 1 (Figure 2B–E). The result of an optimization run is an ensemble of 100 parameter sets that give the lowest cost function. To account for cyclin E and Cdk2 genetic deletion experiments, we examined each parameter set of the ensemble under the cyclin E or Cdk2 deletion condition by setting the respective protein synthesis rates to zero. The network is robust to variation in most of the parameters under the deletion condition. However, owing to variations in the parameters for interactions where cyclin A is involved, there are model solutions that do not satisfy the deletion condition. Solutions were accepted for the parameter sets for which cyclin A at 16 h rose to at least 20% of the activation level of the unperturbed simulation for each deletion. Out of the 100 parameter sets, six fulfilled the requirement of required cyclin A levels for these parameter sets under the unperturbed condition, cyclin E deletion and Cdk2 deletion (see Supplementary Figure 2). In addition, all of these solutions showed no effect on pRb inactivation when p27 was eliminated from the model (data not shown). This is consistent with the activating ‘modifier' as an essential regulator of Cdk2/1 activation in our model.

To test the ability of the model to predict the outcome of an intervention, we examined the effect of Cdk2/1 inactivation. Prediction of Cdk2/1 inhibition was accomplished by reducing the association rates of Cdks to its substrates by a specified percent fraction. The binding constant of the active Cdk complexes with the substrate pRb, kbE2pRB, kbA2pRB and kbA1pRB, were simultaneously reduced by different amounts of inhibition. The model predicts that at greater than 80% inhibition of Cdk2/1, the fraction of active hypophosphorylated pRb at 16 h increases, whereas inactive hyperphosphorylated pRb decreases (Figure 2E). These results were consistent with validation experiments by acute Cdk2/1 inactivation with Roscovitine (Senderowicz, 2003; Figure 1E). Roscovitine treatment of cells abolished Cdk2/1 activity and pRb remained in its active, hypophosphorylated form. In contrast, control cells contained active Cdk2/1 and inactive, hyperphosphorylated pRb. At the low dose used here (15 μM), no alterations of Erk or cyclin D:Cdk4/6 activity were detected (Figure 1E). Cyclin E:Cdk2 kinase activity peaks during early S phase (Figure 1B). However, pRb becomes hyperphosphorylated concomitant with the initial cyclin E:Cdk2/1 activation at/near the restriction point, suggesting that only a low level of cyclin E:Cdk2 activity is required to completely inactivate a cell's complement pRb (Figures 1B and 2E). These data also effectively exclude the potential for the modifier to act on a speculated early G1-phase pRb phosphatase that would keep pRb in an active state. In that case, the prediction would be that inactivation of downstream Cdk2/1 kinase activity would result in hyperphosphorylated (inactive) pRb, as the putative early G1 phosphatase would be inactivated by the modifier. However, the experimental data in Figure 1E clearly show that Cdk2/1 inactivation maintains pRb in the active state, and thereby excludes the involvement of a speculative pRb phosphatase. Moreover, a pRb early G1 phosphatase is also without precedent in the literature. Thus, our model accurately predicted the experimental outcome for both complete Cdk2 inactivation and partial Cdk2 inhibition (Figure 3).

Figure 3

Model of G1 cell-cycle progression. Combining the physiological data (Figure 1) with results from the data-driven mathematical model (Figure 2), we propose the following: growth factor stimulation of quiescent G0 cells activates cyclin D:Cdk4/6 complexes that hypophosphorylate both p130 and pRb resulting in E2F transition from p130 to active, hypophosphorylated pRb in early G1. Continued mitogenic signaling results in increased cellular growth and accumulation of mass, reaching a critical threshold that activates the ‘modifier' leading to activation of pre-existing cyclin E:Cdk2 complexes, followed by pRb hyperphosphorylation (inactivation), release of E2F transcription factors and induction of late G1/S-phase-specific genes and transition into late G1.

Most current explanations of G1 cell-cycle progression are based primarily on data generated from serum deprivation/restimulation experiments and overexpression studies. In serum restimulation experiments, cyclin D:Cdk4/6 activity gradually increases during the first early G1 phase. However, in proliferating cells, cyclin D:Cdk4/6 activity is essentially constitutive in early and late G1 and S phases. In contrast, pRb oscillates between the active, hypophosphorylated form in early G1 and the inactive, hyperphosphorylated form in late G1 and S phases of the cell-cycle phases (see Figure 1) (Fang ; Ezhevsky , 2001; Brugarolas ; Nagahara ). The G1 cell-cycle paradigm will need to be adjusted to account for these experimental observations. Thus, either two completely independent regulatory networks controlling G1 cell-cycle progression exist—one regulating early G1 phase in cells coming from quiescence and another regulating G1-phase transitions in cycling cells—or, there is a common mechanism that regulates G1-phase progression that has been overshadowed by the ramping up of cyclin D:Cdk4/6 activity observed at the beginning of early G1 phase in serum restimulation experiments. The data-driven mathematical approach (Christopher ) presented here supports a single, regulatory mechanism where cyclin E:Cdk2 activation is independent of cyclin D:Cdk4/6. Our mathematical model is in agreement with recent reports showing genetic ablation of either all three D-type cyclins (Kozar ) or their cognate Cdk4/6 kinases (Malumbres ) has limited effect on the timing of cyclin E:Cdk2 activation. Importantly, these reports also exclude p27Kip1 sequestration by cyclin D:Cdk4/6 complexes as the rate-limiting step for cyclin E:Cdk2 activation (Pagano and Jackson, 2004; Sherr and Roberts, 2004).

Inhibition of metabolic processes results in a G1 arrest typified by continued active cyclin D:Cdk4/6 complexes, active hypophosphorylated pRb and inactive cyclin E:CDk2 complexes, suggesting that metabolism plays a key, yet undefined role in G1 cell-cycle progression. Our mathematical model is consistent with evidence linking cell growth with regulation of proliferation (Tapon ; Saucedo and Edgar, 2002). It postulates the existence of an unspecified activating ‘modifier' pathway that responds to cellular growth and results in activation of cyclin E:Cdk2 complexes at the restriction point allowing cells to progress through the cell cycle. Consistent with the presence of an activating modifier element in the mammalian G1 cell cycle, the yeast G1 cell cycle is positively regulated by both cyclin:cdk complexes (Cln3:Cdc28) and an uncharacterized activator, Bck2 (Newcomb ; Costanzo ; de Bruin ). Bck2 is regulated by metabolic stimuli and genetically functions below Cln3 (cyclin D homologue) and activates the equivalent of cyclin E:Cdk2 complexes (Cln1/2:Cdc28) at the yeast restriction point (START). Bck2 is a non-cyclin, non-Cdk protein of unknown function that potentially allows for metabolic activation of G1 cell-cycle progression. Unfortunately, similar to the lack of homology between pRb and the equivalent yeast transcription factor repressor, Whi5, a metazoan Bck2 homologue has not yet been identified. Elucidating the biochemical mechanism of Cdk2 activation and ascertaining the identity of the activating ‘modifier' proposed here is an important step in investigating the effects of therapeutics targeting G1 cell-cycle progression to advance drug development.

Materials and methods

Cell culture and cell-cycle synchronization

Human HCT116 colon carcinoma cells (ATCC) were grown in DMEM (Life Technologies), 5% FBS (Sigma) and antibiotics, at 37°C. Cells were density-arrested by plating at 5 × 105 cells/cm2 for 48 h in 5% FBS, trypsinized and replated at low density (1 × 105 cells/cm2). Cell-cycle progression was assayed by using propidium iodide FACS. Cells were treated with 15 μM roscovitine (EMD Biosciences) at replating.

Immunoblotting, immunoprecipitation and kinase assays

Immunoblots were performed using anti-cyclin E (HE12), anti-cyclin A (H432), anti-cyclin D (H295), Cdk2 (M20 and D12), anti-p27 (C19), anti-actin (I19), anti-p130 (C-20) (Santa Cruz Biotech) and anti-pRb (554136; BD Biosciences) antibodies. Immunoprecipitations were performed using anti-Cdk6 (C-21), anti-Cdk4 (C-22) and anti-Cdk2 (M20) antibodies (Santa Cruz). Immunoprecipitation–kinase reactions were performed using GST-pRb as a substrate for Cdk4/6 and histone H1 (Sigma) for Cdk2.

Reverse transcriptase–PCR

RNA was purified (Qiagen), and RT–PCR was performed (Qiagen) using the following primers:

cyclin A2 (GGCCGAAGACGAGACGGGTTGCACC),

(CAGGCCAGCTTTGTCCCGTGAC);

DHFR (ATGCCTTTCTCCTCCTGG), (CGCTAAACTGCATCGTCGC); and

β-actin (TGAACCCCAAGGCCAACCGCGAGAA), (AAGCAGCCGTGGCCATCTCTTG);

cyclin D1 (GCCTGAACCTGAGGAGCCCC), (GTCTAAGTCAGAGATGGAAGG);

cyclin E1 (TAAAGTGGCGTTTAAGTCCCC), (ATCTTCATCAGCGACGCC).

Numerical simulations and parameter optimization

Simulation and parameter estimation were performed with a custom software written in C++ (Press ; Maimon and Browning, 2001) on a cluster of 16 Linux computers, each containing two 2.0 GHz Intel® Xeon™ processors. The software converted the diagram in Figure 2 into a list of reactions between model components and created a system of coupled differential equations for concentrations of protein species. Differential equations were solved using CVODE (Cohen and Hindmarsh, 1996). For parameter estimation, we used a differential evolution algorithm (Storn, 1999). The algorithm maintains a population of parameter sets while minimizing the cost function. At each iteration, population members with intolerably high costs are discarded, and new ones are created on the basis of lower cost. At the end of an optimization run, an ensemble of parameter sets with lowest costs was stored for further analysis and simulations (Supplementary information).

Competing interest statement

TH, RAC, AD, RM, RG, SVA and IGK declare competing financial interests: details accompany the manuscript. BM, NY and SFD declare that they have no competing financial interests.

Supplementary Figure 1

Supplementary Figure 2

Supplementary Methods

Supplementary Table 1A

Supplementary Table 2

Supplementary Data 1

CellDesigner file for Supplementary Figure 1

Supplementary Data 2

model in SBML format