Robustness of cell cycle control and flexible orders of signaling events

The highly robust control of cell cycles in eukaryotes enables cells to undergo strictly ordered G1/S/G2/M phases and respond adaptively to regulatory signals; however the nature of the robustness remains obscure. Specifically, it is unclear whether events of signaling should be strictly ordered and whether some events are more robust than others. To quantitatively address the two questions, we have developed a novel cell cycle model upon experimental observations. It contains positive and negative E2F proteins and two Cdk inhibitors, and is parameterized, for the first time, to generate not only oscillating protein concentrations but also periodic signaling events. Events and their orders reconstructed under varied conditions indicate that proteolysis of cyclins and Cdk complexes by APC and Skp2 occurs highly robustly in a strict order, but many other events are either dispensable or can occur in flexible orders. These results suggest that strictly ordered proteolytic events are essential for irreversible cell cycle progression and the robustness of cell cycles copes with flexible orders of signaling events, and unveil a new and important dimension to the robustness of cell cycle control in particular and to biological signaling in general.


Supplementary Methods (1) Non-dimensionalization of equations
The concentration of a protein [V] is determined by Transformation includes phosphorylation and dephosphorylation (for ubiquitination, activation, or repression), association and dissociation. We use the following substitutions to non-dimensionalize and scale all equations by the level of APC By assuming both the synthesis and decay rates of APC 1.0, non-dimensionalization and scaling facilitate parameter estimation. The s, k, kk, p, u, and d parameters indicate rates of synthesis, association, dissociation, phosphorylation and dephosphorylation, ubiquitination, and decay. The a and r parameters in Hill functions indicate activation and repression. (

2) Equations with biological basis
The equations are formulated upon several assumptions, which are made based on experimental findings and widely adopted in previous studies. First, because cyclins bind to Cdks with high binding affinities and their dissociation half-lives are in the order of hours (Hochegger et al 2008), we do not handle dissociation and physial decay of cyclin/Cdk complexes. Instead, cyclin/Cdk complexes are destroyed by timely ubiquitination. Other proteins undergo a physical decay, and may also be ubiquinated if clear evidence exists. Second, because Cdks are presence in excess during the cell cycle (Reis et al 2004;Tyson et al 2008;Morgan 2007), we let the whole amount of Cdk1 and Cdk2 at the maximal level (Cdk1tot=Cdk2tot=1.0), and CycE/Cdk2, CycA/Cdk1, and CycB/Cdk1 binding depends on available free Cdk2 and Cdk1. Third, because competition between substrates, instead of the cooperation among multiple phosphorylation sites, makes main contribution to the ultrasensitive mutual protein phosphorylation and dephosphorylation (Kim and Ferrell 2007), we do not handle multiple phosphorylation states of Rb (Yang et al 2006;Haberichter et al 2007;Barik et al 2010), but simply let Rb have two states -active and inactive, and the transition between the two states is controlled by the competition between its modifiers. Fourth, similarly, we just consider CDKB's and Stg's active and inactive states (corresponds to all sites phosphorylated or dephosphorylated), and the transition between the two states is controlled by the competition between their modifiers (Wee, Stg, and active CDKB).

The CycE/Cdk2/Dap system
GF (Growth factor) and E2F1 stimulate cyclin E expression independently (Yao et al 2008). Skp2, one of SCF (the name of the key ubiquitin-protein ligase including Skp1, cullin, and the F-box)'s main F-box proteins, is used to ubiquitinate G1/S targets (CycE and CDKE) (Moberg et al 2001). Decay of p21 and p27 is also through SCF-Skp2 (Amati and Vlach 1999;Guardavaccaro and Pagano 2006). In Drosophila, Skp2 interacts physically with Dap and targets Dap for ubiquitination (Dui et al 2013). Dap binds to CDKE to repress CDKE's activity (de Nooij et al 1996). Since in Drosophila eye dap is expressed widely ahead of the MF (morphogenetic furrow) (Baker 2007), dap is given a constant expression level. F-box proteins have basal expression level (Moberg et al 2001;see Gerard and Goldbeter 2009). Clear evidence indicates that at least in mammalian cells APC (the anaphase-promoting complex, the other main ubiquitin-protein ligase) degrades Skp2 in G1, which represents a principal mechanism by which APC maintains the G1 state (Bashir et al 2004;Wei et al 2004;Guardavaccaro and Pagano 2006).
Another recently identified master regulator of protein destruction in G1/S phases is CRL4Cdt2 (Haven and Walter 2011). In Drosophila, the relative roles of CRL4Cdt2 and Skp2 in the destruction of Dap remain poorly resolved, the destruction of E2Fs by CRL4Cdt2 may not occur in other species (Haven and Walter 2011), and how CRL4Cdt2 itself is regulated is unclear, except that it is quickly down-regulated after S phase (Zielke et al 2011;Haven and Walter 2011). In mammals, multiple pathways restrain E2F activities, and the function of CRL4Cdt2 depends absolutely on the DNA-bound PCNA that is not explicitly described (Haven and Walter 2009). Due to these, we allowed interactions to be concise without including CRL4Cdt2 at this stage.
In mammals, unlike E2F1, the negative E2F protein E2F4 lacks a nuclear localization domain, is free to exit the nucleus, and needs p27/p130 (but not Rb) to co-locate to E2F-responsive promoters of target genes (Rayman et al 2002). p27 can bind to Cdk complexes by its N-terminal domain and to E2F4 by its C-terminal domain (Pippa et al 2012). At mid-late G1 when the expression of these genes is needed, the p27 (associated with p130/E2F4) would recruit Cdk complexes and these Cdk complexes (at late G1 cyclin D-Cdk4/6, and subsequently cyclin E-Cdk2) phosphorylate p130 and disrupt the p27/p130/E2F4 complexes (Malumbres and Barbacid, 2005;Pippa et al 2012). In p130 double mutants, transcription of the B-myb, cyclin A, cdc2 (cdk1), and E2F1 genes is significantly derepressed in comparison with wild-type controls, but the cyclin E gene is not (Rayman et al 2002). Thus, p27 not only maintains cyclin E/cdk2 complexes inactive during early G1 to prevent premature entry into S phase, but also involves in the repression of genes necessary for DNA replication (Macaluso et al 2006;Plesca et al 2007;Pippa et al 2012). Upon these findings in mammalian cells (instead of the interaction between E2F2 and DP in Drosophila), in this model (E2F2 represents the negative E2F protein), we use Dap to represent p27 and p130, let DapE2F2 repress CycA, CycB, and E2F1, and let CDKE dissociate DapE2F2. Thus, CDKE activates E2F1 by breaking RbE2F1 and represses E2F2 by breaking DapE2F2.
It was suggested that after released from its association with Rb, E2F1 is phosphorylated by CDKE and CDKB (Reis and Edgar 2004) and CDKA (Xu et al 1994;Kitagawa et al 1995) for degradation. In Drosophila, it was later found that E2F1 is periodically destructed by CRL4Cdt2, and, by checking putative Cdk phosphorylation sites in E2F1, it was suggested that this destruction does not require CDKs (Shibutani et al 2008 E2F1 transcriptionally activates CycA, which binds to Cdk1 to form CDKA. CycA and CDKA are ubiquitinated by APCFzr and APCFzy (Wolthuis et al 2008). Since rux is expressed widely in Drosophila eye disc (Escudero 2007), Rux has a constant expression level. Rux reversibly binds to and inhibits CDKA (Foley et al 1999;Avedisov et al 2000). CDKE down regulates Rux in vivo (Thomas et al 1997).

The CycB/Cdk1/Stg/Wee system
As in previous models CycB has a constant expression level (Tsai et al 2008;Gerard and Goldbeter 2009), but its expression is also repressed by p130/E2F4 (Plesca et al 2007). CycA and CycB bind competitively to Cdk1. CDKB is inactivated by Wee, activated by the active Stg, and ubiquitinated by APCFzy and APCFzr (Morgan 2007). The loss of kinase activity of CDKB at the end of mitosis depends on the destruction of the cyclin subunits (Felix et al 1990). Wee has a constant expression level, and newly expressed Wee is active. The activities of Wee are high during most of the cell cycle, but then decrease abruptly during mitosis after phosphorylation first by CDKA and then by CDKB (Watanabe et al 2005, Fung et al 2007. stg has a base level of expression (Edgar et al 1994), is also activated by E2F1 (Lehman et al 1999;Reis and Edgar 2004;Yao et al 2011), and newly expressed Stg is inactive. Stg is initially activated by CDKA (Fung et al 2007) and later by active CDKB.
Plx is needed for activation of APC by CDKB (Tsai et al 2008). CDKB turns Plx into the active form upon available inactive Plx, which is assumed at a high level (Tsai et al 2008 (

3) Parameters with biological basis
Supplementary  (Trunnell et al 2011). We choose a moderate value for all nonlinear processes. 2. We scale variables by the syhtnesis and decay levels of APC. 3. Cdk2, Cdk1, and Plx have the same total amount. 4. Since the expression of E2F1, Stg, and CycB is also regulated by E2F1 and E2F2, these contant syhthesis rates should be relatively small to allow the regulation to play significant roles. 5. These proteins are assumed to be expressed at the high levels, and these settings allow these parameters to be removed. 6. These settings allow some parameters to be removed. Because E2F1's level is dynamically degraded by multiple proteins, its decay rate should be small. In contrast, for E2F2, which is not targeted for degradation by other proteins, a large decay rate is necessary. 7-8. Binding rates should be larger than unbinding rates, and these settings allow unbinding rates to be removed. 9. These settings allow some parameters to be removed, to let the nonlinear processes controlled mainly by the thresholds in the Hill functions. 10. Thresholds in negative Hill functions should be small. 11. Skp2 should ubiquitenates G1 cyclins and Cdks timely. 12. These parameters are tuned based on the levels of CycA, CDKA, CycB, CDKBi, and CDKBa. 13. The CDKE and CDKB mediated degradation of E2F1 should be strong to allow cell cycle phase compensation to occur. Parameters beginning with a and r are coefficients in Hill functions controlling non-linearity of protein interactions. 14. We assume that Skp2  We deliberately adopted a set of very simple initial conditions (Supplementary Table 4). Unrealistic though, they indicate that the model does not demand specific initial conditions. These initial conditions do not allow the program XPPAUT to identify a fixed point for bifurcation analysis. We identified another set of initial conditions that are biologically more reasonable and allow XPPAUT to identify a fixed point of the model (Supplementary Table 4).
To reliably analyze the robustness of the model against changes of parameters and of events, we carefully explored constraints among parameters and identified two more sets of parameters that make the model generate both oscillating protein concentrations and all signaling events (Supplementary Table 3). Notice that, while many parameter values make the model generate oscillating protein concentrations, very few make the model generate all signaling events. This indicates that the three sets of parameters are biologically more reasonably constrained. We find that to generate all signaling events, significant changes of some half-maximal activating and inhibiting coefficients in Hill functions are not allowed, this explains why considerable half-maximal activating and inhibiting coefficients in the second and third set of parameters do not differ much. Supplementary

Robustness and bifurcation analysis
Bifurcation analysis has been widely used to analyze the qualitative dynamics of the cell cycle control system (Tyson et al 2002). Previous studies have revealed hysteresis and bistability of the system caused by cell mass, CycA, CycB, and cycle number (Csikasz-Nagy et al 2006;Gerard and Goldbeter 2009;Pomerening et al 2003;Calzone et al 2007). In this work we focused on the influence of half-maximal activating and inhibiting coefficients (in Hill functions, which control the timing of events) in the cell cycle control system, especially E2F1-related events, because E2F1 is the core of regulatory feedbacks, and E2F2-related events, because this negative E2F protein is firstly quantitatively examined in this model.
We used the program oscill8 and the model's default parameters to perform bifurcation analysis. First, we found that many half-maximal activating and inhibiting coefficients do not cause the model to generate hysteresis and bistability indicated by the CycE level, that fixed points occur when these coefficients are near their default values, and, especially, that Fzr-related coefficients show no influence on CycE (Supplementary Figure 6). Second, the responses of CycE to changes of a CDKEE2F1 and a CDKBaE2F1 that control the degradation of E2F1 by CDKE and CDKB, to changes of a E2F1E2F1 that controls the activation of E2F1 by E2F1, and to changes of d E2F1 that controls the decay of E2F1, show hysteresis (Supplementary Figure  8). Third, DapE2F2 negatively regulates the production of E2F1, CycA, and CycB. Only changes of r DapE2F2E2F1 make the model generate complex CycE responses (Supplementary Figure 7). We also examined the responses of E2F1, CycA, and CycB to changes of r DapE2F2E2F1 , r DapE2F2CycA , and r DapE2F2CycB , and found that r DapE2F2E2F1 makes not only CycE but also CycA generate hysteresis-like responses. Equation (12) indicates that the production of CycA is sensitive to both E2F1 and E2F2. Together, these bifurcation analyses suggest that E2F1 and E2F2 decisively influence the dynamic properties of the cell cycle control system.
Next, we analyzed the robustness of the model against changed timing of signaling events. By running the model with the second and third set of parameters we repeated the analysis shown in Table 3 and generated Supplementary Table 3, in which ↗ and ↘ indicate 66.6% increase and decrease of the control parameter, 0 indicates that oscillating protein concentrations were not generated (cell cycle fails), 1 indicates that oscillating protein concentrations were generated but some events were absent or present persistently, and 2 indicates that oscillating protein concentrations were generated and all signaling events occurred periodically. Thus, the result of 0/0 indicates that the model is highly sensitive to the changed timing of the event, and the result of 2/2 indicates that the model is robust against the changed timing of the event.
We classified all events into 9 groups (Supplementary Table 5 in addition that the means of the CDKB group is also large. Although Levene's test indicates that not all populations have a common variance, the Student-Newman-Keuls test and Duncan's multiple range test equally indicate that Fzy, Fzr, and CDKBa are in the same group whose means are significantly larger than others. These statistics therefore support that the proteolysis of Cyclin and Cdk complexes by Fzr and Fzy are more robust compared with other events. Supplementary