Search tips
Search criteria 


Logo of plosbiolPLoS BiologySubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)View this Article
PLoS Biol. 2010 September; 8(9): e1000488.
Published online 2010 September 21. doi:  10.1371/journal.pbio.1000488
PMCID: PMC2943438

Stochastic E2F Activation and Reconciliation of Phenomenological Cell-Cycle Models

Andre Levchenko, Academic Editor


The transition of the mammalian cell from quiescence to proliferation is a highly variable process. Over the last four decades, two lines of apparently contradictory, phenomenological models have been proposed to account for such temporal variability. These include various forms of the transition probability (TP) model and the growth control (GC) model, which lack mechanistic details. The GC model was further proposed as an alternative explanation for the concept of the restriction point, which we recently demonstrated as being controlled by a bistable Rb-E2F switch. Here, through a combination of modeling and experiments, we show that these different lines of models in essence reflect different aspects of stochastic dynamics in cell cycle entry. In particular, we show that the variable activation of E2F can be described by stochastic activation of the bistable Rb-E2F switch, which in turn may account for the temporal variability in cell cycle entry. Moreover, we show that temporal dynamics of E2F activation can be recast into the frameworks of both the TP model and the GC model via parameter mapping. This mapping suggests that the two lines of phenomenological models can be reconciled through the stochastic dynamics of the Rb-E2F switch. It also suggests a potential utility of the TP or GC models in defining concise, quantitative phenotypes of cell physiology. This may have implications in classifying cell types or states.

Author Summary

Mammalian cells enter the division cycle in response to appropriate growth signals. For each cell, the decision to do so is critically dependent on the interplay between environmental cues and the internal state of the cell and is influenced by random fluctuations in cellular processes. Indeed, experimental evidence indicates that cell cycle entry is highly variable from cell to cell, even within a clonal population. To account for such variability, a number of phenomenological models have been previously proposed. These models primarily fall into two types depending on their fundamental assumptions on the origin of the variability. “Transition probability” models presume that variability in cell cycle entry originates from the fact that entry in each individual cell is random but also governed by a fixed probability. In contrast, “growth-controlled” models assume that the growth rates across a population are variable and result in cells that are out of phase developmentally. While both kinds of models provide a good fit to experimental data, their lack of mechanistic details limits their predictive power and has led to unresolved debate between their practitioners. In this study, we developed a mechanistically based stochastic model of the temporal dynamics of activation of the E2F transcription factor, which is used here as a marker of the transition of cells from quiescence to active cell cycling. Using this model, we show that “transition probability” and “growth-controlled” models can be reconciled by incorporation of a small number of basic cellular parameters related to protein synthesis and turnover, protein modification, stochasticity, and the like. Essentially our work shows that each kind of phenomenological model holds true for describing a particular aspect of the cell cycle transition. We suggest that incorporation of basic cellular parameters in this manner into phenomenological models may constitute a broadly applicable approach to defining concise, quantitative phenotypes of cell physiology.


Cell-to-cell variability in the timing of cell-fate commitment is widely observed in biological settings [1][4]. In particular, the variable timing of transition from the quiescent to the proliferative state is a well-documented phenomenon [5][8]. In a population of proliferating cells, such variability is reflected in the partitioning of the population into subpopulations at various phases of the cell cycle. This phenomenon is observed even in a population of isogenic cells that have been synchronized by serum starvation. Upon growth stimulation, cells reenter the cell cycle from quiescence and undergo the G1/S transition, but not all cells in the population proceed at the same rate. This rate also differs among different cell types [9],[10] and can be modulated by external conditions [11].

To account for the variable transition timing in cell cycle progression, two major types of models have been proposed: the transition probability (TP) models [11][15] and the growth-controlled (GC) models [16][18]. The TP models attributed temporal variability to random state transitions through different phases of the cell cycle. One of the earliest TP models was proposed to account for the inter-mitotic variability by assuming a single random transition from a non-proliferative A-state to a proliferative B-phase [15]. It was subsequently extended to account for the timing variability in cell cycle reentry starting from quiescent (G0) cells [11],[14]. In this case, the exponential drop in the fraction of G0 cells over time was suggested to indicate a probabilistic nature of the transition. The original model and its subsequent variants have provided excellent fits to various types of experimental data [11][15]. However, a major criticism of the TP models is that the transition probability from the A-state was assumed to be time-invariant, despite uneven cell division at mitosis and obvious cell growth or metabolism through the cell cycle [19]. As an alternative, the GC models proposed that the observed temporal variability arises from growth rate heterogeneity within a cell population, rather than random state transitions. Remarkably, this line of models has been able to provide equally good fits to various experimental data [16],[20]. Integrating these two lines of thinking, hybrid models proposed cell-size control and random transitions as regulatory elements for progression to cell division [21],[22]. However, understanding of the underlying mechanisms for cell-size control and random transitions was limited at the time. Consequently, although they provided excellent fits to experimental data, these models remain descriptive to date.

There has been an active debate between these two lines of thinking since their initial propositions. While never fully resolved, the debate gradually faded as the focus in the field of cell cycle studies moved to identifying the dynamical basis for various cell cycle regulations, including the restriction point (R-point) [23], which we have shown to be controlled by a bistable Rb-E2F switch [7]. We also showed that activation of this switch is correlated with the cell's reentry from quiescence into the cell cycle. Interestingly, cell cycle reentry was explored by both the TP and GC models, which were originally developed to describe actively growing cells. For example, the TP models ascribe quiescence and proliferation to low and high transition probabilities, respectively [11],[14]. In addition, the GC models have recently been proposed as an alternative explanation for the “R-point” [18].

The temporal variability described by the GC and TP models is based on the distribution of inter-mitotic times and may differ from temporal variability in E2F activation from quiescence. However, we suggest that the stochastic Rb-E2F model embodies the concepts assumed in the TP models and the GC models. Our model predictions and experiments suggest that stochastic activation of E2F can account for temporal variability in cell cycle entry, and the degree of such variability is determined by environmental cues and the regulatory network parameters. These results suggest that the TP and GC models are not mutually exclusive but rather reflect different aspects of the same temporal dynamics in cell cycle entry, as has been speculated [21],[24]. In addition, we show that stochastic activation of the Rb-E2F bistable switch under various environmental conditions can readily be mapped into both TP and GC models with a small number of parameters (Figure 1). We propose that these parameters can potentially serve as concise, quantitative phenotypes of cell states.

Figure 1
Temporal variability in cell cycle reentry.


Our recent work has shown that traverse of the R-point is regulated by the Rb-E2F bistable switch [7]. This bistability results from interlocked positive feedback loops embedded in the Myc-Rb-E2F network (Figure S1, see Text S1 for further details). Given the bistable switching property of the Myc-Rb-E2F network, we hypothesized that this network, when subjected to noise, might demonstrate variable timing in E2F activation, which in turn might account for the temporal variability observed in cell cycle entry. This hypothesis is based on the strong correlation we previously observed between E2F activation and DNA synthesis [7]. To test this hypothesis, we developed a stochastic model for the Myc-Rb-E2F network using the chemical Langevin formulation [25],[26] as detailed in Materials and Methods. This formulation allows for implementation of intrinsic and extrinsic noise while retaining the deterministic framework. In this stochastic model, the intrinsic noise arises from the stochastic nature of the biochemical interactions among small numbers of signaling molecules. The extrinsic noise results from heterogeneity in cell size and shape, cell division, or cell cycle stage [27][32].

Modulation of E2F Activation by Serum Stimulation

The fluctuations in the bistable switch result in significant discrepancies between stochastic and deterministic simulations [33][36]. Given a set of initial conditions and parameters in the Myc-Rb-E2F network, the simulated time-courses from a deterministic model are fixed (black line in Figure 2A), but those from a stochastic model show drastically variable trajectories (gray lines in Figure 2A). For example, the stochastic Rb-E2F model can generate two modes of E2F at Time = 50 h when stimulated with weak input as shown in Figure 2A. We define a switching threshold (horizontal red line in Figure 2A) to distinguish the low E2F mode, which corresponds to a non-activated subpopulation of cells, from the high E2F mode, which represents an activated population. This threshold can be used to calculate the percentage of activated cells over time. The minimum time required for E2F to reach the switching threshold is defined as the switching time (vertical red line in Figure 2A, for the deterministic simulation). Similarly, for strong input, the deterministic time-course simulations are fixed and stochastic time-course simulations again show variable trajectories (unpublished data). The distribution of E2F activity in stochastic simulations, however, exhibits a single mode (high E2F level) at strong inputs, rather than two modes as with weak inputs.

Figure 2
Simulated temporal dynamics in E2F activation using the stochastic Rb-E2F model.

Based on our simulations and definitions in Figure 2A, we obtained G0 exit curves for weak and strong input conditions as shown in Figure 2B. These G0 exit curves are analogous to the α-curve in the TP model, which represents the frequency distribution of inter-mitotic times [15]. Both a G0 exit curve and an α-curve can be fitted by an exponential curve with two parameters (black dotted curve in Figure 2B, see Text S2 for further details): transition rate (KT) and time delay (TDP). This is because both exhibit an initial time delay followed by an exponential drop [11],[14],[15]. The transition rate of the G0 exit curve is inversely proportional to the temporal variability of the cell population. For example, a population of cells with more-synchronous E2F activation would have a higher transition rate than that of a population with less-synchronous E2F activation. If cells were completely synchronized, the G0 exit curve would have an infinite transition rate.

Our simulated E2F activation dynamics predict serum-dependence of both transition rate and time delay. For a weak input (KT = 0.029±0.0014 h−1 and TDP = 18.0±1.2 h, blue line in Figure 2B), most cells were expected to remain inactivated and the percentage of G0 cells would decrease slowly. This is because the impact of noise acting on the Rb-E2F bistable switch was only significant enough to activate E2F in some cells, but not in other cells. This would lead to a bimodal distribution of E2F activity (Figure S2A), which is consistent with previous experimental observations in mouse fibroblasts [13],[37],[38]. In contrast, the impact of noise was negligible in the case of strong input and all cells were predicted to be activated at a high transition rate (KT = 0.16±0.0076 h−1 and TDP = 7.7±0.27 h, red curve in Figure 2B).The response of the Rb-E2F bistable switch to noise would cause an increase in KT with increasing input strength (Figure 2C) as the population moves from a bimodal to a monomodal distribution at the high mode (Figure S2A). At sufficiently high input strength, further increase in input strength may have a negligible effect on KT (Figure 2C). In contrast, TDP may decrease as the population transitions from a bimodal to a monomodal distribution, and TDP may bottom out at sufficiently high input strength (Figure 2D). The dependence of KT and TDP on input strength can be recapitulated with a minimal bistable model (Figure S2B–D), suggesting that such dependence may be an intrinsic property of bistable systems.

To validate our model predictions, we measured E2F activity in the E2F-d2GFP cell line, which is derived from a rat embryonic fibroblast REF52 cell line and carries a destabilized green fluorescent protein reporter (d2GFP) under the E2F1 promoter. We have shown that this reporter system can be used to monitor E2F activity in response to serum stimulation previously [7]. Prior to serum stimulation, the E2F-d2GFP cells were synchronized at quiescence by serum-starvation (0.02% bovine growth serum, BGS) with basal E2F-GFP expression (Figure 3A). Upon weak serum stimulation (0.3% BGS), only a subpopulation of the cells switched to the high E2F mode over time. At earlier time points (0~15th h), the difference in E2F level between the non-activated and activated cells was small. The difference between the two modes became increasingly clear, resulting in distinctive bimodality starting at 18th h. In contrast, upon strong serum stimulation (5% BGS), E2F activation was more synchronous. The whole cell population gradually switched to the high mode with greater temporal synchrony without demonstrating detectable bimodality at any tested time point (Figure 3A). It is possible that noise may partition the cell population into two subsets (active and inactive towards proliferation) temporarily even at high serum stimulation. However, simulations suggest that accumulation of E2F in the activated cells at earlier time points may not be significant enough to result in any detectable difference between the two subsets (Figure S2A).

Figure 3
Experimental validation of the predicted temporal dynamics.

Based on the distribution of E2F in Figure 3A, we calculated the percentage of non-activated cells and obtained a G0 exit curve for each serum condition (Figure 3B). Consistent with predictions in Figure 2B, we observed an increase in KT and decrease in TDP for increasing serum concentration (KT = 0.031±0.0036 h−1 and TDP = 5.1±1.1 h at 0.3% serum, and 0.16±0.011 h−1 and 1.1±0.27 h at 5% serum), reminiscent of modulation of the α-curve by serum [11],[15]. Consistent with model predictions in Figure 2C–D, we observed increase in KT and decrease TDP for increasing serum concentrations (Figure 3C and 3D). An independent experiment under the same conditions on a different day exhibited similar dependence of KT and TDP on serum (Figure S3).

Modulation of Stochastic E2F Activation by Strength of CycE-Mediated Feedback

The temporal dynamics of biological systems often depend strongly on network parameters [39].Consequently, the transition rate of cell cycle entry may be modulated by nodal perturbations. This is exemplified in a recent study on the yeast cell cycle [40], which demonstrated that a positive feedback by G1 cyclins is responsible for temporal coherence in gene expression and proper division timing of yeast cells. Loss of this feedback control in the cell cycle machinery was shown to promote incoherent gene expression and abnormal delays of yeast budding. Interestingly, a similar feedback module through a G1 cyclin (CycE) can be found in the Myc-Rb-E2F network also, suggesting its potential role in the control of temporal dynamics.

To investigate modulation of the transition rate by nodal perturbations in the Myc-Rb-E2F network, we introduced in silico perturbations of one particular node: the CycE/Cdk2 complex, which forms the CycE-mediated positive feedback loop. Our bifurcation analyses predict that weakening of the CycE-mediated positive feedback loop will desensitize the Rb-E2F bistable switch to serum stimulation, requiring a higher critical serum concentration (Figure 4A) for E2F activation. Similarly, we predict desensitization to serum when CycD is down-regulated or when Rb is up-regulated (unpublished data). Such desensitization is expected to modulate the temporal dynamics of E2F activation. When the positive-feedback strength by CycE is weakened, our simulations in Figure 4B (corresponding simulated distributions in Figure S4) predicted increase in the time delay and decrease in the transition rate. For strong feedback strength, KT was estimated to be 0.17±0.0090 h−1. This value was reduced to 0.15±0.006 h−1 and 0.12±0.0054 h−1 for intermediate and weak feedback strength, respectively. In contrast, TDP for strong feedback input ( = 7.4±0.25 h) was predicted to increase to 8.5±0.53 h for intermediate feedback strength, and to extend further to 10.8±0.15 h for weak feedback strength. Similar dependence of KT and TDP on the feedback strength was predicted for all serum concentrations (Figure 4C and 4D).

Figure 4
Nodal perturbation at CycE alters temporal dynamics of E2F activation.

To test these predictions experimentally, we perturbed the Myc-Rb-E2F network by applying varying concentrations of a cyclin-dependent kinase inhibitor (CVT-313), which has a much higher affinity towards Cdk2 than to other Cdks (Figure S5) [41],[42]. In the context of the current study, which focuses on the cellular dynamics leading to E2F activation, the impact of the Cdk2 inhibitor is primarily the inhibition of the CycE/cdk2 complex. We note that the inhibitor would also affect other components of cell cycle regulation, (e.g., the CycA/cdk2 complex), which were not considered in the model due to their activity mainly downstream of the cell cycle entry point. When the CycE node was perturbed experimentally, we observed inhibitor dose-dependent changes in E2F activity, as measured by GFP fluorescence in the E2F-d2GFP cells. As shown in Figure 5A, increasing dose of the inhibitor drug reduced the fraction of cells in the high E2F mode at 24 h. For example, without the Cdk2 inhibitor, less than 1% serum was required for E2F activation in half of the cell population. With 2 µM Cdk2 inhibitor, 2% serum was required to achieve a similar fraction of E2F activation. Such desensitization to serum stimulation was seen for all inhibitor concentrations tested (Figure 5A).

Figure 5
Experimental desensitization of the Rb-E2F switch to serum by a Cdk2 inhibitor drug.

Next, we tested modulation of temporal dynamics by the Cdk2 inhibitor. At 2% serum, we applied the Cdk2 inhibitor (CVT-313, 2 µM) to monitor its effect on E2F activation over time. Our results in Figure 5B show that the transition rate of the cell population decreased (from KT = 0.078±0.0073 to 0.058±0.0070 h−1) and time delay increased (from TDP = 9.1±0.70 to TDP = 12.0±0.86 h) with addition of the Cdk2 inhibitor. Such a decrease in KT with the inhibitor drug is consistent with our model predictions in Figure 4 and was observed for all serum concentrations tested, as shown in Figure 5C (distributions of E2F in Figure S6). As predicted, time delay generally decreased with increasing serum concentrations and it increased in the presence of the Cdk2 inhibitor (Figure 5D). It is noteworthy that the estimated time delay has a large error at low serum concentrations, leading to a non-monotonic dependence of TDP on serum concentrations. This is most likely due to the small number of E2F-activated cells at low serum at earlier time points. This makes estimation of parameters using least squares challenging, giving rise to large errors. We conducted another experiment on a different day under the same experimental conditions and observed similar trends in KT and TDP, as shown in Figure S7.

Mapping Simulated Stochastic E2F Activation into TP and GC Model

Throughout this study, we have analyzed the temporal dynamics of E2F activation by extracting a set of parameters defining the TP model (transition rate and time delay). This parameter extraction establishes a connection with the mechanistic Rb-E2F model. Similarly, the GC model parameters (mean growth rate An external file that holds a picture, illustration, etc.
Object name is pbio.1000488.e003.jpg and its variance An external file that holds a picture, illustration, etc.
Object name is pbio.1000488.e004.jpg, see Text S2 for further details) can be extracted from the stochastic dynamics of E2F activation, and a connection between the GC model and the mechanistic Rb-E2F model can also be established. The GC parameters were estimated from both simulation results in Figure 2 and experimental data in Figure 3, as shown in Table 1. These results show increasing mean growth rate and decreasing variance (normalized to the mean) with increasing input strength.

Table 1
Extraction of GC and TP model parameters from simulated and experimentally measured G0 exit curves.

In addition, we predicted the dependence of the strength of the CycE-mediated positive feedback on the GC model parameters, as shown in Figure 6. Consistent with Table 1, our results predicted increasing growth rate (Figure 6A) and decreasing normalized variance (Figure 6B) for increasing input strength. However, decreasing the strength of the CycE-mediated positive feedback was predicted to reduce mean growth rate without affecting its normalized variance significantly. Such parameter extraction defining the phenomenological models provides a quantitative mapping between the phenomenological models and the mechanistic Rb-E2F model. It is noteworthy that both TP and GC models fit the data with comparable levels of uncertainty in the estimated parameters, suggesting that both models may provide similarly good fits to the stochastic dynamics of E2F activation.

Figure 6
Mapping the stochastic dynamics of E2F activation with the GC model.


Focusing on E2F activation, we have shown that the temporal variability in cell cycle entry from quiescence can be quantitatively modeled by stochastic activation of a bistable Rb-E2F switch [7]. In addition, we have shown that the degree of such variability can be modulated by varying the input strength or by perturbing the network parameters.

Our model predictions are overall consistent with experimental measurements. In particular, our analysis indicates that serum and a Cdk2 inhibitor drug exert opposite influences on the temporal dynamics of E2F activation: transition rate increases and time delay decreases with increasing serum, but transition rate decreases and time delay increases with increasing Cdk2 inhibitor concentrations. We suggest that such a well-calibrated stochastic model for the Rb-E2F switch may guide further experimental analyses to gain insights into the system-level dynamics underlying cell cycle entry. For example, our model predicts that reducing the CycD/Cdk4,6 activity may have similar effects on temporal dynamics of E2F activation as the Cdk2 inhibitor, while knocking down Rb may increase transition rate (unpublished data). In addition, we can predict stochastic dynamics of E2F activation under combinatorial perturbations including growth factors, inhibitor drugs targeting the Myc-Rb-E2F network, or mutations within this network.

Throughout this study, we have focused on a single transition during cell cycle progression (quiescence to proliferation) due to its experimental and computational tractability. To further simplify analysis, we have chosen not to model cell division or growth explicitly. Instead, the variability associated with these processes is lumped into the extrinsic noise terms in our SDE model. More explicit mechanisms to account for such variability may further improve the quantitative agreement between the modeling and the experiment. For example, our simulation results suggest that the major source of noise is extrinsic noise, while variability in the initial conditions can lead to minor yet discernable change in the temporal dynamics of E2F activation. This is evident when E2F activation dynamics are compared under two conditions: varying initial conditions and varying variance of the extrinsic noise (ω) in the stochastic model (see Materials and Methods). At a fixed variance of the extrinsic noise, increasing variability in the initial conditions (Gaussian-distributed with the mean being the base initial conditions and various variance values, Var) is predicted to decrease transition rate and time delay (Figure S8A–B). Similarly, increasing ω without any variability in the initial conditions (Var = 0) is predicted to decrease transition rate and time delay (Figure S8C–D), but these changes by extrinsic noise are predicted to be significantly greater than those by initial conditions. These decreases in KT and increases in TDP reflect the loss of synchrony in E2F activation due to increasing extrinsic noise or initial condition variability. This may explain reduced time delay in actively growing cells compared to that in quiescent cells [14].

Equally important, we further show that these predicted stochastic dynamics of the Rb-E2F model can be quantitatively mapped into two lines of phenomenological models reflecting seemingly conflicting views: the TP model and the GC model. For a given set of parameters defining the stochastic model, the simulated stochastic E2F activation at the population level can be uniquely described by a set of parameters defining the TP model or the GC model (compare Figure 4C–D and Figure 6A–B). Furthermore, different sets of parameters in the stochastic model would lead to different parameters in the TP or the GC models. We propose that this mapping provides a simple conceptual framework that reconciles the different views reflected in the TP and GC models, which have been a source of unresolved debate over the last several decades. In other words, the stochastic model can be considered as a common mechanistic basis for the two seemingly different models.

During the mapping from our stochastic model to the TP or GC models, details associated with individual signaling reactions are necessarily lost in the resulting TP or GC models, pointing to their limitations in offering mechanistic insights. However, a by-product of this mapping is a potential, unappreciated utility of the TP and GC models. On one hand, these phenomenological models are simple and are able to provide quantitative description of the population-level dynamics associated with variable cell cycle entry. On the other, specific changes in the underlying reaction networks can be manifested in changes in the parameters in these simple models. As such, together with a mechanistically based model, the TP and GC models can serve as a concise platform to define quantitative phenotypes that facilitate classification of cell types or cell states.

This utility may be particularly useful for cancer diagnosis, since most cancers have defects in the Myc-Rb-E2F signaling pathway [43],[44]. Recent approaches for cancer classification involve microarray-based gene expression profiling to develop cancer signatures [45], which have been used to reveal the activation status of oncogenic signaling pathways [46]. Here we suggest that oncogenic phenotypes resulting from deregulation in these pathways may also serve as cancer signatures. Using the mapping technique defined in this work, we can develop a library of predicted phenotypes (defined as TP or GC model parameters) based on the Myc-Rb-E2F network under various nodal mutations or stimulatory inputs. This library can be correlated with the oncogenic phenotypes (defined as TP or GC model parameters) of an unknown cancer cell type. In principle, this correlation can be used to infer the activation status of the Myc-Rb-E2F network of the cancer cell type. For a small number of test conditions, this may be challenging owing to the stochastic dynamics of cell cycle entry. However, increasing the number of test conditions may enhance the diagnostic potential of this approach.

Materials and Methods

Development of a Stochastic Rb-E2F Model

The deterministic version of the Rb-E2F model, developed in our previous work [7], served as a basis for the stochastic Rb-E2F model. To capture stochastic aspects of the Rb-E2F signaling pathway, we adopt the chemical Langevin formulation [25],[26],[47] as shown in Eqn (1).

equation image

where Xi(t) represents the number of molecules of a molecular species i (i = 1, …, N) at time t, and X(t) = (X1(t), …, XN(t)) is the state of the entire system at time t. X(t) evolves over time at the rate of aj[X(t)] (j = 1, …, M), and the corresponding change in the number of individual molecules is described in vji. An external file that holds a picture, illustration, etc.
Object name is pbio.1000488.e011.jpg and An external file that holds a picture, illustration, etc.
Object name is pbio.1000488.e012.jpg are temporally uncorrelated, statistically independent Gaussian noises. This formulation retains the deterministic framework (the first term), and reaction-dependent and reaction-independent noise. The concentration units in the deterministic model were converted to molecule numbers, so that the mean molecule number for E2F would be approximately 1,000. We assumed a mean of 0 and variance of 1 for Γj (t), and a mean of 0 and appropriately determined variance for ωj (t) (see Text S1 for more details). The resulting stochastic differential equations (SDEs) were implemented and solved in Matlab.

Cell Culture Conditions

Actively growing E2F-d2GFP cells [7] were serum-starved in Dulbecco's modified Eagle's medium (DMEM) with 0.02% of bovine growth serum (BGS). After 24 h of serum starvation, they were stimulated with varying serum concentrations for cell cycle entry in the presence or absence of Cdk2 inhibitor CVT-313 (from Calbiochem: Cat #238803). Cell cycle progression was blocked at the DNA synthesis stage by hydroxyurea (HU block), which we have shown has insignificant impact on the GFP signal [7]. At various time points, these cells were collected and fixed in 1% formaldehyde for fluorescence assay.

Fluorescence Assay with Flow Cytometry

E2F-d2GFP rat embryonic fibroblasts were assayed for a destabilized green fluorescent protein reporter (d2GFP) for E2F activity. The intensity of d2GFP was measured with a flow cytometry system (BD FACSCanto II).

Western Blots

E2F-d2GFP cells were serum-starved (BGS = 0.02%) for 24 h before they were treated with varying concentration of the Cdk2 inhibitor (CVT-313, EMD # 238803) and serum. After 24 h of serum/inhibitor drug treatment, cell lysates were collected and Western blotting was conducted with primary antibodies recognizing Rb phosphorylation at Cdk4-specific serine 780 (Santa Cruz, #sc-12901-R) and at Cdk2-specific threonine 821 (Santa Cruz, #sc-16669-R). These were conjugated with anti-rabbit secondary antibodies (GE Healthcare, #NA934) for detection. As a loading control, actin was measured with actin-recognizing primary antibodies (Santa Cruz, #sc-8432) conjugated with anti-mouse secondary antibodies (GE Healthcare, #NA9310).

Supporting Information

Figure S1

A schematic of the Rb-E2F bistable switch. Here, Rb represents the entire Rb family (pRB, p107, and p130) and E2F represents all activating E2Fs (E2F1, E2F2, and E2F3a). In quiescent cells E2F is bound by Rb and its transcriptional activities are repressed. Growth stimulation removes the Rb repression by upregulating cyclin D (CycD), which, in complex with Cdk4,6, phosphorylates Rb to release E2F. In addition, growth stimulation induces a transcription factor Myc that upregulates CycD. The free form of E2F synergizes with Myc to induce its own transcription, forming feed-forward and positive feedback loops. Subsequently, E2F activates the transcription of Cyclin E (CycE), which forms a complex with Cdk2 to further remove Rb repression by phosphorylation, constituting another positive feedback loop.

(0.39 MB TIF)

Figure S2

Simulated temporal dynamics of E2F activation by the full model and a minimal model. (A) The Rb-E2F bistable switch was stimulated with weak (S = 0.5) and strong (S = 5) input strengths. E2F distributions from 1,000 simulations were sampled at various time points for both conditions. For weak input strength, bimodality was predicted to emerge at around 16th hour. At strong input strength, however, bimodality was expected to be less clear. (B) A minimal model can be used to recapitulate the temporal dynamics of the bistable Rb-E2F switch. The model describes activity of a molecule X: An external file that holds a picture, illustration, etc.
Object name is pbio.1000488.e013.jpg, where S is the input strength, ka ( = 5) is the lumped rate term for synthesis and feedback strength, and kb ( = 0.1) is a basal synthesis term. Bifurcation analysis of the minimal bistable model shows hysteresis. (C) This minimal model was converted to a stochastic model using the chemical Langevin formulation. The transition rates were calculated for cell populations stimulated at various input strengths. The transition rate increased with input strength and reached a plateau at sufficiently high input strength. (D) In the minimal bistable model, the time delay decreased with increasing input strength and reached a plateau at sufficiently high input strength.

(0.71 MB TIF)

Figure S3

Independent time-course measurements of GFP signal reporting activity. GFP signal under the same experimental conditions as Figure 3 was measured at varying time points with flow cytometry. (A) For 0.3% serum, the transition rate and time delay were estimated to be 0.022±0.0041 h−1 and 10.0±1.8 h, respectively. At high serum concentration ( = 5%), the transition rate increased to 0.11±0.0099 h−1 and time delay decreased to 7.9±0.55 h. (B) KT increased with serum concentration. (C) TDP decreased with serum concentration.

(0.28 MB TIF)

Figure S4

Predicted modulation of the temporal dynamics of E2F activation. Temporal dynamics of E2F activation were simulated at varying input strengths (weak→S = 0.5, intermediate→S = 1, and strong→S = 5) and varying CycE-mediated positive feedback strengths (strong→kP4 = 18 h−1 and weak→kP4 = 9 h−1). With strong positive feedback (PFB), bimodality was predicted for weak input while monomodality (E2F ON) was predicted for intermediate and strong stimulations. With weak positive feedback, the percentage of E2F activation was predicted to decrease for weak and intermediate input strengths. For strong input, however, the effect of the positive feedback strength was minor.

(0.54 MB TIF)

Figure S5

Specificity of the Cdk2 inhibitor. To demonstrate the effect of the Cdk2 inhibitor on Cdk2 kinase activity, we measured Rb phosphorylation at the Cdk2-specific and Cdk4-specific residues for varying inhibitor concentrations. An isogenic population of serum-starved E2F-d2GFP cells was used for Western blotting. In serum-starvation condition (serum = 0.02%), Rb phosphorylation at either residue was negligible. With serum stimulation (serum = 10%), a significant increase in Rb phosphorylation at both residues was observed. For increasing Cdk2 inhibitor concentration, Rb phosphorylation efficiency decreased at the Cdk2-specific residue, but no significant change was observed at the Cdk4-specific residue.

(0.48 MB TIF)

Figure S6

Experimentally measured E2F time courses for varying serum concentrations, in the absence or presence of the Cdk2 inhibitor drug (at 2 µM). At 0th h E2F-d2GFP cells were synchronized in quiescence by serum-starvation (24 h at 0.02% serum), stimulated with varying serum concentrations (with or without the Cdk2 inhibitor drug), and measured for GFP (reporting E2F activity) by flow cytometry at the indicated time points.

(1.62 MB TIF)

Figure S7

E2F time-courses under varying serum concentrations in the absence or presence of Cdk2 inhibitor (at 2 µM) performed on a separate day. (A) G0 exit curves in the absence and presence of the Cdk inhibitor (2% BGS). Addition of the inhibitor reduced the transition rate from 0.10±0.0081 to 0.090±0.0091 h−1 and increased the time delay from 7.7±0.55 to 9.1±0.67 h. (B) Transition rate as a function of serum concentration in the absence or presence of the Cdk2 inhibitor. (C) Time delay as a function of serum concentration in the absence or presence of the Cdk2 inhibitor. (D) E2F distribution over time in the presence or absence of the Cdk2 inhibitor.

(0.80 MB TIF)

Figure S8

Variability in the initial conditions versus in the rates of the chemical reactions. The effects of variability in the initial conditions and in the rates of the chemical reactions were evaluated on the temporal dynamics of E2F activation. With all else the same, our simulation results predicted that transition rate (A) and time delay (B) would decrease significantly as ω was increased from 25 to 50. To describe variability in the initial condition, we assumed that the initial concentrations for Rb and the Rb-E2F complex were Gaussian-distributed with the mean being their base value and varying variance levels. At a fixed variance of extrinsic noise (ω = 50), our simulation results predicted that transition rate (C) and time delay (D) would decrease slightly with increasing variance of the initial conditions. Overall, the activation dynamics of E2F is much more sensitive to changes in extrinsic variability than those in the initial condition.

(0.25 MB TIF)

Text S1

Model development.

(0.27 MB DOC)

Text S2

Calculation of metrics: TP and GC model parameters.

(0.04 MB DOC)


We would like to thank Tony Huang, Shwetadwip Chowdhury, and William Kim for preliminary analysis and discussions, and Jin Wang and Jeff Wong for helpful suggestions.


growth control
transition probability


The authors have declared that no competing interests exist.

This work is partially supported by National Institutes of Health (1P50GM081883,, a David and Lucile Packard Fellowship (, a DuPont Young Professor Award (, and the Duke Vertical Integration Program for summer research funded through Howard Hughes Medical Institute ( The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


1. Nachman I, Regev A, Ramanathan S. Dissecting timing variability in yeast meiosis. Cell. 2007;131:544–556. [PubMed]
2. Spencer S. L, Gaudet S, Albeck J. G, Burke J. M, Sorger P. K. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature. 2009;459:428–432. [PMC free article] [PubMed]
3. Bean J. M, Siggia E. D, Cross F. R. Coherence and timing of cell cycle start examined at single-cell resolution. Mol Cell. 2006;21:3–14. [PubMed]
4. Colman-Lerner A, Gordon A, Serra E, Chin T, Resnekov O, et al. Regulated cell-to-cell variation in a cell-fate decision system. Nature. 2005;437:699–706. [PubMed]
5. Brooks R. F. Cytoplasmic origin of variability in the timing of S-phase in mammalian-cells. Cell Biol Int Rep. 1979;3:707–716. [PubMed]
6. Brooks R. F. Regulation of fibroblast cell cycle by serum. Nature. 1976;260:248–250. [PubMed]
7. Yao G, Lee T. J, Mori S, Nevins J. R, You L. C. A bistable Rb-E2F switch underlies the restriction point. Nat Cell Biol. 2008;10:476–U255. [PubMed]
8. Smith J. A. The cell-cycle and related concepts in cell-proliferation. Journal of Pathology. 1982;136:149–166. [PubMed]
9. Zetterberg A, Larsson O. Kinetic-analysis of regulatory events in G1 leading to proliferation or quiescence of Swiss 3t3 cells. Proc Natl Acad Sci U S A. 1985;82:5365–5369. [PubMed]
10. Zetterberg A, Larsson O, Wiman K. G. What is the restriction Point. Curr Opin Cell Biol. 1995;7:835–842. [PubMed]
11. Shields R, Smith J. A. Cells regulate their proliferation through alterations in transition-probability. J Cell Physiol. 1977;91:345–355. [PubMed]
12. Shields R. Further evidence for a random transition in cell-cycle. Nature. 1978;273:755–758. [PubMed]
13. Shields R. Transition-probability and origin of variation in cell-cycle. Nature. 1977;267:704–707. [PubMed]
14. Brooks R. F, Bennett D. C, Smith J. A. Mammalian-cell cycles need 2 random transitions. Cell. 1980;19:493–504. [PubMed]
15. Smith J. A, Martin L. Do cells cycle. Proc Natl Acad Sci U S A. 1973;70:1263–1267. [PubMed]
16. Castor L. N. A G1 rate model accounts for cell-cycle kinetics attributed to transition-probability. Nature. 1980;287:857–859. [PubMed]
17. Cooper S. On G0 and cell-cycle controls. Bioessays. 1987;7:220–222. [PubMed]
18. Cooper S. Reappraisal of serum starvation, the restriction point, G0, and G1 phase arrest points. Faseb Journal. 2003;17:333–340. [PubMed]
19. Koch A. L. The re-incarnation, re-interpretation and re-demise of the transition probability model. J Biotechnol. 1999;71:143–156. [PubMed]
20. Koch A. L. Does the variability of the cell cycle result from one or many chance events? Nature. 1980;286:80–82. [PubMed]
21. Tyson J. J, Hannsgen K. B. Cell growth and division: a deterministic/probabilistic model of the cell cycle. J Math Biol. 1986;23:231–246. [PubMed]
22. Tyson J. J, Hannsgen K. B. Global asymptotic stability of the size distribution in probabilistic models of the cell cycle. J Math Biol. 1985;22:61–68. [PubMed]
23. Pardee A. B. A restriction point for control of normal animal cell proliferation. Proc Natl Acad Sci U S A. 1974;71:1286–1290. [PubMed]
24. Nurse P. Cell-cycle control - both deterministic and probabilistic. Nature. 1980;286:9–10. [PubMed]
25. Gillespie D. T. The chemical Langevin equation. J Chem Phys. 2000;113:297–306.
26. Tanouchi Y, Tu D, Kim J, You L. Noise reduction by diffusional dissipation in a minimal quorum sensing motif. PLoS Comput Biol. 2008;4:e1000167. doi: 10.1371/journal.pcbi.1000167. [PMC free article] [PubMed]
27. Elowitz M. B, Levine A. J, Siggia E. D, Swain P. S. Stochastic gene expression in a single cell. Science. 2002;297:1183–1186. [PubMed]
28. Rosenfeld N, Young J. W, Alon U, Swain P. S, Elowitz M. B. Gene regulation at the single-cell level. Science. 2005;307:1962–1965. [PubMed]
29. Raser J. M, O'Shea E. K. Control of stochasticity in eukaryotic gene expression. Science. 2004;304:1811–1814. [PMC free article] [PubMed]
30. Paulsson J. Summing up the noise in gene networks. Nature. 2004;427:415–418. [PubMed]
31. Raser J. M, O'Shea E. K. Noise in gene expression: origins, consequences, and control. Science. 2005;309:2010–2013. [PMC free article] [PubMed]
32. Volfson D, Marciniak J, Blake W. J, Ostroff N, Tsimring L. S, et al. Origins of extrinsic variability in eukaryotic gene expression. Nature. 2006;439:861–864. [PubMed]
33. Srivastava R, You L, Summers J, Yin J. Stochastic vs. deterministic modeling of intracellular viral kinetics. J Theor Biol. 2002;218:309–321. [PubMed]
34. Vellela M, Qian H. Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlogl model revisited. J R Soc Interface. 2009;6:925–940. [PMC free article] [PubMed]
35. Kaern M, Elston T. C, Blake W. J, Collins J. J. Stochasticity in gene expression: from theories to phenotypes. Nature Reviews Genetics. 2005;6:451–464. [PubMed]
36. Wang X, Hao N, Dohlman H. G, Elston T. C. Bistability, stochasticity, and oscillations in the mitogen-activated protein kinase cascade. Biophys J. 2006;90:1961–1978. [PubMed]
37. Brooks R. F, Richmond F. N, Riddle P. N, Richmond K. M. V. Apparent heterogeneity in the response of quiescent Swiss 3t3 cells to serum growth-factors - implications for the transition-probability model and parallels with cellular senescence and competence. J Cell Phys. 1984;121:341–350. [PubMed]
38. Brooks R. F, Riddle P. N. The 3t3 cell-cycle at low proliferation rates. J Cell Sci. 1988;90:601–612. [PubMed]
39. Hoffmann A, Levchenko A, Scott M. L, Baltimore D. The I kappa B-NF-kappa B signaling module: temporal control and selective gene activation. Science. 2002;298:1241–1245. [PubMed]
40. Skotheim J. M, Di Talia S, Siggia E. D, Cross F. R. Positive feedback of G1 cyclins ensures coherent cell cycle entry. Nature. 2008;454:291–U212. [PMC free article] [PubMed]
41. Bhattacharjee R. N, Banks G. C, Trotter K. W, Lee H. L, Archer T. K. Histone H1 phosphorylation by cdk2 selectively modulates mouse mammary tumor virus transcription through chromatin remodeling. Mol Cell Biol. 2001;21:5417–5425. [PMC free article] [PubMed]
42. Brooks E. E, Gray N. S, Joly A, Kerwar S. S, Lum R, et al. CVT-313, a specific and potent inhibitor of CDK2 that prevents neointimal proliferation. J Biol Chem. 1997;272:29207–29211. [PubMed]
43. Malumbres M, Barbacid M. To cycle or not to cycle: a critical decision in cancer. Nat Rev Cancer. 2001;1:222–231. [PubMed]
44. Nevins J. R. The Rb/E2F pathway and cancer. Hum Mol Genet. 2001;10:699–703. [PubMed]
45. Golub T. R, Slonim D. K, Tamayo P, Huard C, Gaasenbeek M, et al. Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. Science. 1999;286:531–537. [PubMed]
46. Bild A. H, Yao G, Chang J. T, Wang Q. L, Potti A, et al. Oncogenic pathway signatures in human cancers as a guide to targeted therapies. Nature. 2006;439:353–357. [PubMed]
47. Lee T. J, Tan C. M, Tu D, You L. C. Modeling cellular networks. 2007. Bioinformatics: An Engineering Case-Based Approach Chapter 6.

Articles from PLoS Biology are provided here courtesy of Public Library of Science