Transcriptional Regulation Is a Major Controller of Cell Cycle Transition Dynamics

DNA replication, mitosis and mitotic exit are critical transitions of the cell cycle which normally occur only once per cycle. A universal control mechanism was proposed for the regulation of mitotic entry in which Cdk helps its own activation through two positive feedback loops. Recent discoveries in various organisms showed the importance of positive feedbacks in other transitions as well. Here we investigate if a universal control system with transcriptional regulation(s) and post-translational positive feedback(s) can be proposed for the regulation of all cell cycle transitions. Through computational modeling, we analyze the transition dynamics in all possible combinations of transcriptional and post-translational regulations. We find that some combinations lead to ‘sloppy’ transitions, while others give very precise control. The periodic transcriptional regulation through the activator or the inhibitor leads to radically different dynamics. Experimental evidence shows that in cell cycle transitions of organisms investigated for cell cycle dependent periodic transcription, only the inhibitor OR the activator is under cyclic control and never both of them. Based on these observations, we propose two transcriptional control modes of cell cycle regulation that either STOP or let the cycle GO in case of a transcriptional failure. We discuss the biological relevance of such differences.


H. sapiens
Cdk2/CycE,A Wee1hu (Parker et al, 1992;Watanabe et al, 2005;Watanabe et al, 1995) hCdc25a (Boutros et al, 2006) A (Donzelli et al, 2003;Lukas et al, 2004;Mailand et al, 2000) A (Hoffmann et al, 1994) Cdk2/CycE (Lundberg and Weinberg, 1998) PP1 (Durfee et al, 1993) I (Wade Harper et al, 1993) I (Geng et al, 1996) Table S1. Cell cycle transition regulators in various organisms, including references for all claims. Similar to Table 1 of the main text, just here we add references that support our claims on the roles of activator and/or inhibitors in checkpoint regulation (ChP) and their positive feedback (PFB) interactions with the transition regulators (TR). Checkpoint regulation (ChP) and positive feedback loop (PFB) notation: Aacting through activator, I-through inhibitor, B-through both of them. Blue, bold letters note genes that are periodically expressed during the cell cycle in system wide studies (Gauthier et al, 2010); controversially to this hCdc25c was found static in individual experiments (Donzelli et al, 2003). Note that all regulations are by phosphorylation -dephosphorylation reactions, with activators being phosphatases and inhibitors being kinases, except two reverse systems, noted by # . At the M/G1 transition of fission yeast cells Wee1 is activated and Cdc25 is inactivated after dephosphorylation by Clp1 (Wolfe et al, 2004), similarly Wee1hu is activated and Cdc25c is inhibited during the human M/G1 transition (Potapova et al, 2009) these effects are lumped together in one row of this table. (INH) sign and italic letters for the whole row means the TR is an inhibitor of the cell cycle transition, thus all effects on it are acting with reverse sign to the transition.

Model development
Here we present in detail the various models we used in the paper. Figure S1 shows a complete picture of all the main species and all their interactions (0 < i < 29) the kinetics laws of which are detailed in Table S2.

Figure S1: Detailed interaction map of the generic cell cycle transition regulator.
Solid arrows represent reactions, while dashed arrows represent activating effects. In brackets we note the reactions where transcription factors (TFs) and checkpoint proteins (ChPs) act.
In the model notation we follow, TR* and TR represent, respectively, the active and inactive forms of the transition regulator, inhibitor* and inhibitor represent, respectively, the active and inactive forms of the inhibitor and activator* and activator represent, respectively, the active and inactive forms of the activator. In detail, TR* activates activator and inhibits inhibitor*. Although TR (the inactive form of the transition protein) acts on the activator and the inhibitor in the same way as TR*, its efficiency is 100 times lower, as proposed by others after theoretical and experimental observations (Ciliberto et al, 2007;Deibler et al, 2010). Moreover, note that both the activity of TR* (and TR) on inhibitor* and on activator form positive feedbacks (PFB); activator* activates TR* and inhibitor inhibits it. The checkpoint moves inhibitor to a hyperactive form, #inhibitor that is four times stronger inhibitor of TR* than the normally active (but not checkpoint affected) inhibitor* form. This way of implementing checkpoint activation was taken from an earlier model of morphogenetic checkpoint of budding yeast cells (Ciliberto et al, 2003). Table S2: Kinetic laws associated to reactions, corresponding to notation on Figure S1. The values associated to basal parameters are reported in Table S3.
Syntheses and degradations follow the law of mass action, while all the other reactions follow the Michaelis-Menten kinetics; parameters starting with character k are the catalytic constants (dimension 1/min), while parameters starting with j are the Michaelis constants (dimensionless). Moreover, since we describe systems in terms of explicit molecule counts and not in terms of concentrations, we introduce a scaling factor α which is used to represent the size of the system -volume of the cell, as described earlier (Mura and Csikasz-Nagy, 2008). All the molecular species that are present in the expressions of Tables S2 but not in Figure S1 have constant population, representing background activities of constantly present kinases and phosphatases. In particular, E1 activates the inhibitor (and also #inhibitor), while E2 inhibits the activator*, both are needed to give a threshold for the autocatalytic TR* induced autocatalysis. Moreover, S is representing a background signal, activating activator in case the PFB A is removed and there is no activation of activator by TR*. Species ChP A and ChP I are introduced to represent checkpoints on the activator and the inhibitor, respectively. While ChP A promotes the inhibition of the activator, ChP I transforms the inhibitor into the more active #inhibitor form (see above).
When we say that TF I acts on the inhibitor we introduce the synthesis reaction R 10 and all the corresponding degradation reactions R 9 , R 11 , R 23 and R 24 . In the same way, when we say that TF A acts on the activator we refer to the synthesis reaction R 16 and all the corresponding degradation reactions R 15 and R 17 . In case there is no transcriptional regulation of the activator or inhibitor then we keep the above parameters at 0 and use a constant high (500 molecules) total activator or inhibitor level (see Table S4 for details).
Explanation of model notation.   The 24 different models represent combinations of the three regulatory effects (transcription, positive feedback and checkpoint): • we considered two models, denoted by 6 7 and 6 8 , that have PFB both on the inhibitor and the activator, have no checkpoints (i.e., ChP A and ChP I set to 0) and are such that one (6 7 ) has TF acting only on the activator and the other (6 8 ) has TF acting only on the inhibitor. Then we consider models pairs 6 7 7 , 6 8 7 and 6 7 8 , 6 8 8 that are obtained from 6 7 and 6 8 by adding checkpoint ChP A or checkpoint ChP I , respectively; • we considered two models, denoted by $ 7 and $ 8 , that have PFB only on the inhibitor, have no checkpoints (i.e., ChP A and ChP I set to 0) and are such that one ($ 7 ) has TF acting only on the activator and the other ($ 8 ) has TF acting only on the inhibitor.
Then we consider models pairs $ 7 7 , $ 8 7 and $ 7 8 , $ 8 8 and $ 7 9 , $ 8 9 that are obtained from $ 7 and $ 8 by adding, respectively, checkpoint ChP A , checkpoint ChP I or both the checkpoints -here we needed to introduce molecule S to induce some activation of activator in the absence of PFB on activator; • we considered two models, denoted by -7 and -8 , that have PFB only on the activator, have no checkpoints (i.e., ChP A and ChP I set to 0) and are such that one (-7 ) has TF acting only on the activator and the other (-8 ) has TF acting only on the inhibitor. Then we consider models pairs -7 7 , -8 7 and -7 8 , -8 8 that are obtained from -7 and -8 by adding checkpoint ChP A or checkpoint ChP I , respectively.

Model implementation
All the models presented in the previous section have been implemented in BlenX (Dematté et al, 2010) and simulated by means of the Beta Workbench (Dematte et al, 2008). The peculiar characteristics of BlenX allowed us to write a complete model from which, thanks to the modularity of the language, all the different models are obtained by commenting the fragments of code we are not interested in.
The basic metaphor of BlenX is that a biological entity (i.e., a component that is able to interact with other components to accomplish some biological functions) is represented by a computational device called box. A box has a set of interfaces and an internal program and can be graphically represented. Interfaces have associated types and they are used by a box to interact with other boxes; the internal program, instead, codifies for the set of actions that a box can perform after a specific interaction with another box in the system has happened. For example, if a box is modeling a protein, its interfaces may represent sensing and effecting domains.

Boxes as abstractions of biological entities.
The small squares on the border of the box are the interfaces; x and y are the interface subjects (omitted when not necessary); t1 and t2 are the interface types (omitted when not necessary); P is the internal process and B is the name of the box. Subjects are used by the internal program to refer to interfaces.
Through sensing domains the protein receives signals that are then propagated through the internal actions to activate/inactivate a set of effecting domains. The exchange of signals can happen between boxes whose interfaces have a certain degree of affinity, which codes the strength of their interaction. These affinities are calculated by definable expressions, which can be declared as simple real numbers if the reaction that they are accounting for is an elementary mass action law, or they can be arbitrary functions (e.g., Michaelis-Menten, Hill response,…) if the reaction is an aggregated process whose elementary mechanism of interaction between entities is not known.
We use a small simplified example to introduce briefly the language and the methodology (design pattern) we used to implement our models. Consider two boxes B1 and B2 representing, respectively, activator and TR:

Boxes of activator and TR.
For both boxes interfaces with subject rec represent sensing domains, and are used by our proteins to receive activation and inhibition signals. Interfaces with subject out represent effecting domains, and are used by our proteins to send activation and inhibition signals. Types with subscript 0 represent sensing domains of proteins in active forms, while types with subscript 1 represent sensing domains of proteins in inactive forms. Moreover, types with subscript S represent effecting domains of proteins in active forms, while types with subscript I represent effecting domains of proteins in inactive forms.
Box B1 can execute a sequence of actions, starting with rec?(), in parallel (denoted with | symbol) with the action out!(). The action out!() sends a signal through interface with subject out. The primitive rec?() waits a signal on the interaction site with subject rec that enables the change of the types of both interfaces by means of the sequence of actions ch(rec, tr0).ch(out, trS). The same holds for box B2.
Given the structure of boxes B1 and B2, when types tr1 and cdcS are affine, the two boxes can interact by synchronizing on the corresponding interfaces through the corresponding out!() and rec?() actions. Specifying the affinity between types tr1 and cdcS as an expression representing the kinetics law used to describe reaction (see Table S2), we have that the B1 and B2 interaction follows the Michaelis-Menten kinetics, causing the transformation of box B2 into B2', which represents the active protein TR*:

Boxes of activator and TR* after TR activation.
Specifying the affinity between types trS and act0 as an expression representing the kinetics law used to describe reaction (see Table 1), we have that also B1 and B2' can perform an interaction, which causes the transformation of B1 in a box that represents the inactive protein activator.
This small and simplified example shows how we can implement our protein species as boxes following a common methodology, using simple and reusable internal programs; moreover, it shows how the ability of a protein to activate or inhibit another protein is specified through the specification of the affinity relation between the types. Stochastic simulations are performed by means of the Beta Workbench stochastic simulator (Dematte et al, 2008), which implements a variant of the Gillespie algorithm (Gillespie, 1977) based on the Gibson and Bruck implementation (Gibson and Bruck, 2000).
A complete and detailed description of BlenX can be found in (Dematté et al, 2010), along with a more detailed description of the implementation methodology adopted and here briefly introduced. On the contrary, when TFs are present in low amount (bottom row) then the two models behave totally differently. 1/10 of the normal TF A level seems to be too weak transcriptional input to induce a transition ever. The same 10 times lower transcription of inhibitor cannot delay the transition and TR* turns into its active form soon after it is produced. Thus a transcriptional failure of TF A acts as a STOP signal to the transition, while a failure of TF I gives a GO signal.

Figure S2: Transcriptional STOP and GO controls with positive feedback only on the inhibitor (a) or only on the activator (b). Here we plot the active forms of the key molecules of (a) model I A (right) and I I (left) and (b) A A (right) and A I (left) at high (top) and low (bottom) transcription factor levels. Black: TR*, red: inhibitor*, green; activator*
Similarly, Figure S2a plots simulations of the basic models I A and I I . When TFs are present in high amount (top row) then the transitions look qualitatively the same, the threshold of TR* = 150 is reached approximately after the same time. Note that here the transition is less sharp compared to the one in Figure 2. Although in the figure model I I tends to reach the threshold earlier, values of S signal (30 and 40 for models I A and I I ., respectively) have been chosen to let the two models having the transition in average at approximately the same time. Also in this case when TFs are present in low amount (bottom row) then the two models behave totally differently. Also in this case a transcriptional failure of TF A acts as a STOP signal to the transition, while a failure of TF I gives a GO signal.
Finally, Figure S2b plots simulations of the basic models A A and A I . Also in these cases simulation results show the same dynamics of the ones reported in Figures 2 and S2a, both when TFs are present in high amount and when TFs are present in low amount. Note however that when the TFs are present in high amount, although the transition is sharp, we have that the maximal value of TR* is lower; indeed, since the inhibitor is always present only in its active form, the system is not able to convert all the TR into TR*.
Worth to notice is the effect of the noise on transcription on the noise in protein levels. In the right panels of Fig S2a,b the inhibitor is under transcriptional control and as a result the inhibitor* level is much noisier. Similar, although less obvious increase in activator* noise can be observed when the activator is under transcriptional regulation. On Fig S2b left panel the inhibitor is not affected neither by transcription or feedback regulation, so we ended up using a constant value for its level.

Simulation methods and details on the main figures of the paper
Here we provide details on the different analyses we performed on our models and presented in the figures of the main text. (Figure 3 and Figure S3). We considered model pairs 6 7 /6 8 , $ 7 /$ 8 and -7 /-8 ; in red are depicted results of models with TF I acting on the inhibitor, while in green results of models with TF A acting on the activator. For each of these models we considered 31 different variants where TFs are assumed to be induced at different timings in the set: −45,−40, …,−5, 0, 5, …, 40, 45, 50,…,90,95,100} With time 0 we refer to the time point where the transcription of TR is initiated (TF TR is set from 0 to 500). For times greater than 0, we used a specific feature of BlenX that allowed us to start simulations with TF I and TF A equal to 0 and to change, at the desired time, the amount of TFs from 0 to 500. For times less than 0, we run simulations with initial populations of activator and inhibitor that reflect the assumption of a TF A or TF I induction in advance. In particular, to calculate these initial populations of activator and inhibitor, we used the formula: where s is the synthesis rate, d the degradation rate and t is the transcriptional induction time.

Robustness of the different models to perturbation in timing of transcriptional induction of activator or inhibitor relative to the transcriptional induction of TR
For each of these 186 models we ran 1000 stochastic simulations and measured, for each simulation, the time at which TR* reaches population 150 (which value corresponds to the approximate inflection point on the TR* activation curve, see Figure 4 top panels). We used the stochastic simulations to calculate means and standard deviations of transition times that are presented by solid lines and background shading, respectively, on Figure 3 and Figure S3. Figure 3 in the main text with all the other considered models. (Figure 4 and Figure  S4). We considered models 6 7 , 6 8 , $ 7 , $ 8 , -7 and -8 ; in red are depicted results of models with TF on the inhibitor, while in green results of models with TF on the activator. For each of the models we set parameter kms to an initial value of 0.0001 and used a specific feature of BlenX that allows to update the value of a parameter when time reaches a certain value. In each single simulation we increased the value of kms by adding 0.0001 at times 5000, 10000, 15000, …, 95000 and decreased its value by subtracting 0.0001 at times 100000, 105000, …, 190000. Hence, in a single simulation the value of kms is increased at fixed time intervals from value 0.0001 to value 0.002 and then decreased in the same way from value 0.002 to value 0.0001 again. For each of these 6 models we ran 100 stochastic simulations and for each simulation we calculated in each time interval [(i-1)*5000,i*5000] (with i=1,…,39) the mean values of TR* inside the sub-intervals [i*5000-2000,i*5000]; we calculated hence the steady state values of TR* corresponding to the different values of the kms parameter.

Bistability in cell cycle transitions under various control models
We then calculated means and standard deviations of these steady state values of TR* over the 100 stochastic simulations and obtained the results reported in Figure 4 and Figure  S4. Parameter robustness test of the various models ( Figure 5 and Figure S5). We considered models 6 7 , 6 8 , $ 7 , $ 8 , -7 and -8 ; in red are depicted results of models with TF on the inhibitor, while in green results of models with TF on the activator. Our robustness analysis is inspired by (Barkai and Leibler, 1997). Each of the six models is considered as a reference model and from each of it we generated 1000 variants by randomly modifying some parameters. For models with TF on the activator we modify parameters kcs, kcd, and for models with TF on the inhibitor we modify parameters kws, kwd, * . Each alternation of the reference system is characterized by the total parameter variation, k, which is defined as: where L is the list of the parameters subject to variation, E is the altered parameter and E ) is the corresponding parameter in the reference model. Given a parameter E ) , the altered parameters are obtained by multiplying E ) to a value x, generated using a loguniform distribution in the interval [0.1-10].
For each of the generated 6000 models we ran 100 stochastic simulations and measured, for each simulation, the time at which TR* reaches population 150. For each combination of each model we calculated the mean (over the 100 simulations) of the times at which TR* reaches population 150 and plotted results in Figure 5 and Figure S5.   Figure 6 in the main text withal the other considered models.
ChP A and ChP I . In particular, for models with checkpoint on the activator we generated different models with parameter ! with values in this in the set {0, 0.1, 0.2, 0.3, 0.4}, while for models with checkpoint on the inhibitor we generated different models with parameters ! and ! both having value (the same) in the set {0, 0.1, 0.2, 0.3, 0.4}.
For each of the generated 90 models we ran 1000 stochastic simulations and measured, for each simulation, the time at which TR* reaches population 150 for all models. We used these simulations to calculate means and standard deviations of transition times.
Full bars in the figure indicate that with this parameter value in more than the 90% of the 1000 simulations TR* do not reach population 150; the star means that the percentage is < 90% but is not 0%. Table S75 shows for all the generated models the exact percentages of simulations (over the 1000) where TR* reaches population 150.

Oscillations
We considered models 6 7 , 6 8 , $ 7 , $ 8 , -7 and -8 . To generate oscillations each of the models is extended with a minimal negative feedback loop model shown below. In the figure osc1 and osc1* represent respectively the active and inactive forms of a new intermediary enzyme and osc2 and osc2* represent respectively the active and inactive forms of another new intermediary enzyme. In detail, TR* activates osc1 that, when active, activates osc2 which, when active, induces the degradation of TR* and TR. Such combination of positive and negative feedback loops is supposed to give a robust minimal cell cycle oscillator (Ferrell et al, 2011;Pomerening et al, 2005). Kinetic laws and related parameters are shown in Table S6. Reactions ) , , % and ( follow the law of mass action, while the other reactions follow the Michaelis-Menten kinetics. Parameters have been tuned to let enzyme osc1 to activate fast only when TR* is far above the bistability threshold and, when osc2 is active, to let TR* to stay active for an amount of time that realistically can be enough for a cell to finish mitosis. For each of the generated 6 models we ran simulations of 200 stochastic cell cycles and measured the period of TR* oscillations; A period is considered starting when TR* switches on and reaches population level 150. We used these values, for all the models, to calculate mean and coefficients of variations of the TR* period (see Table S7). Examples of oscillations for all the considered model are shown in Figure S7.
Minimal negative feedback loop model following (Goldbeter, 1991