Whi5 phosphorylation embedded in the G1/S network dynamically controls critical cell size and cell fate

In budding yeast, overcoming of a critical size to enter S phase and the mitosis/mating switch—two central cell fate events—take place in the G1 phase of the cell cycle. Here we present a mathematical model of the basic molecular mechanism controlling the G1/S transition, whose major regulatory feature is multisite phosphorylation of nuclear Whi5. Cln3–Cdk1, whose nuclear amount is proportional to cell size, and then Cln1,2–Cdk1, randomly phosphorylate both decoy and functional Whi5 sites. Full phosphorylation of functional sites releases Whi5 inhibitory activity, activating G1/S transcription. Simulation analysis shows that this mechanism ensures coherent release of Whi5 inhibitory action and accounts for many experimentally observed properties of mitotically growing or conjugating G1 cells. Cell cycle progression and transcriptional analyses of a Whi5 phosphomimetic mutant verify the model prediction that coherent transcription of the G1/S regulon and ensuing G1/S transition requires full phosphorylation of Whi5 functional sites.

T he ability of G 1 cells to select specific cell fates in response to external and internal cues is crucial for life, in both lower and higher eukaryotes 1 . In the presence of nutrients-and in the absence of the complementary mating factors-haploid G 1 cells of the budding yeast Saccharomyces cerevisiae-a widely used model for the study of the eukaryotic cell cycle-grow to the critical cell size required to enter S phase 2,3 and are committed to proliferation. The complementary mating factor induces haploid G 1 cells to undergo cell cycle arrest, differentiate and conjugate producing diploid cells 4,5 . Nutrient-limiting conditions lead to growth arrest and entry into stationary phase 6,7 .
Several features of critical cell size control have been reported. The average cell size of yeast populations increases with ploidy 3,8 . Smaller G 1 cells require a longer G 1 phase before entering S phase as compared to larger G 1 cells [9][10][11][12][13] . The large size variability of newborn daughter cells is substantially reduced when they enter S phase 13 . Cells exponentially growing in rich media supporting faster growth rate have a larger critical cell size than that of cells growing in poor media which support slower growth 11,14 . During a nutritional shift up from a poor to a rich medium, changes in growth rate and cycle progression cause the average cell size of the population to increase, slowly adapting to the size characteristic of the new medium 15 . Finally, several signalling pathways modulate critical cell size [16][17][18][19][20][21][22][23] .
The number of gene products involved in execution and control of growth and cell cycle has been estimated at around 300 for the cell cycle and around 600 for growth 24,25 . Gene products reported to modulate the critical cell size include proteins involved in the G 1 /S transition network (such as Cln3, Whi3, Ydj1, Far1, Whi5, Cln1, Cln2, Sic1and so on) and in ribosome biogenesis (such as Sfp1 and Sch9) 3,26 .
A well-accepted critical cell size control theory assumes two distinct, but interconnected, functions: sensing of actual cell size and setting of the critical cell size according to growth conditions 2,26 . Following this view, two mechanisms-possibly involving a (partially) different set of molecules-would be devoted to 'weigh' the protein content of a cell and to set the actual critical cell size for each growth condition. Cln3, the most upstream cyclin regulating the cyclin-dependent kinase Cdk1, has been proposed as the sensing component 19,21,27,28 . The most recent model in this class holds that a yeast cell may titrate Cln3 molecules against the number of genomic binding sites for SBF, a heterodimeric transcription factor that regulates transcription of the G 1 /S regulon, encoding genes required to execute the G 1 /S transition 28 . Genome-wide screenings for mutants with altered cell size have suggested that ribosome biogenesis may set the cell size value characteristic for each growth rate 26 .
Single-cell analysis suggested a sizer plus timer model, in which Cln3 is involved in a noisy sizer control (whose complementary molecular component is unspecified), which sets a period within G 1 (called T 1 ) of considerable length in small daughter cells, but much shorter in parent cells 13 . T 1 terminates with the export of Whi5 from the nucleus and the consequent transcriptional activation of the G 1 /S regulon. The following T 2 period, which ends the G 1 phase, has similar length in both parent and daughter cells and terminates with bud emergence and the onset of DNA replication 13 . The sensing/setting theory and the sizer plus timer models are reciprocally compatible and considerable effort has been dedicated to identify the counterpart of Cln3 in the sensing mechanism. A recently proposed model substantially departs from this common background to hold that the critical cell size is robustly set by the rate of linear growth during G 1 and that the chaperone Ydj1 has an important role in regulating critical cell size 29 .
In this paper we focus on the core mechanism of cell size control in steady state growth, no consideration being given yet on how changes in nutrients or growth rate affect cell size. We present a novel and conceptually simple mathematical model of the G 1 /S network, centered on the phosphorylation of Whi5, an inhibitor of G 1 /S-specific transcription. Although our model includes all known molecular players of the G 1 /S transition, only few critical parameters (that affect Whi5 phosphorylation) determine its most important functions. Model predictions have been validated by analysing cell cycle entry and transcriptional activation in mutants expressing a Whi5 protein with constitutive (pseudo)phosphorylation. The model accounts for different data sets-including pheromone effects on cell cycle-related eventsthat were not considered in model design and parameter optimization, thereby providing a new unifying, comprehensive, molecular mechanism for the core critical cell size control and for the mitosis/mating switch.

Results
Outline of the model of the G 1 /S transition. In our model ( Fig. 1a; Supplementary Notes 1-10), synthesis of the upstream cyclin Cln3 links cell growth-modelled according to ref. 30-to molecular events promoting the G 1 /S transition. Cln3 synthesis is exquisitely sensitive to alterations in growth rate and ribosomal content 31 and its average G 1 concentration is higher in fast growing cells 19,21,31 . Conflicting evidences regarding the precise pattern of Cln3 protein accumulation during the G 1 phase have been reported [32][33][34][35] . Our simulations show that the G 1 /S transition kinetics are quite insensitive to the specific pattern of Cln3 accumulation ( Supplementary Figs 1-2), while they are dependent on average Cln3 level during G 1 (cln3D and Cln3 overexpression cases). Therefore, we used the simpler hypothesis, that is, that Cln3 concentration is constant throughout the G 1 phase. In early G 1 , the chaperone Ydj1 facilitates import of Cln3 in the nucleus where it binds the catalytic kinase subunit Cdk1 (ref. 27): the nuclear amount of Cln3-Cdk1 is therefore larger in larger cells and smaller in smaller cells. Coherent transcriptional activation of the G 1 /S regulon 36 leads to synchronous expression of hundreds of gene products, orderly driving the G 1 /S transition that starts the pathway towards mitosis. It is driven by two heterodimeric transcription factors, SBF (composed by the coactivator Swi6 and the DNA binding protein Swi4) and MBF (composed by Swi6 and Mbp1 (ref. 37)) and inhibited by the transcriptional repressor Whi5 that regulates the G 1 /S transition in both mitotic 36,38,39 and mating factor-treated cells 5 through inhibition of several hundreds of SBF molecules.
Multisite phosphorylation of a regulatory protein is an effective molecular device able to produce a coherent response in the specific function controlled by the same protein, when the function involves hundreds of identical molecular players, as observed in the control of firing of DNA replication origins by phosphorylation of Sld2/Sld3 (refs 40,41). As Cln3-Cdk1 accumulates, it phosphorylates Whi5 on 12 Cdk phosphoacceptor sites that are mostly occupied in vivo in asynchronously growing cells 42 . Whi5 is released from SBF when its four specific functional phosphosites are fully phosphorylated 42,43 . The other 8 sites-whose phosphorylation does not affect binding activityact as decoys that compete with the functional sites for the available Cdk1 kinase activity. SBF-Whi5 dissociation may also derive from the full phosphorylation of four critical phosphosites 42 of Swi6, a critical target of Cln3 (ref. 44). When Whi5 is released, transcription of the SBF-driven genes of the G 1 /S regulon may start 36,42 . Full phosphorylation of the four Swi6 phosphosites of MBF drives transcription of the MBF-regulated genes of the regulon.
The order of transcription of the genes encoding the more relevant cycle regulatory proteins-which may be fine-tuned in different experimental conditions-is: CLN1, CLN2, CLB5, CLB6 and then NRM1 (ref. 36). A positive-feedback loop that involves CLN1 and CLN2 transcription ultimately commits yeast cells to S phase 45 , while phosphorylated Whi5 is released from SBF and exported from the nucleus 5,13,29,38 . MBF is repressed by a product of the same regulon, Nrm1, which terminates the transcription of the regulon as cells progress towards S phase 46 .
The model includes cytoplasm, endoplasmic reticulum and nuclear sub-cellular compartments. Ordinary Differential Equations (ODEs) based on the mass-action law, modified when required to account for nonlinear behaviour (for example, for Cln3-Ydj1 diffusion into the nucleus, Cln3-dependent Far1 degradation, Sic1 phosphorylation and so on), describe synthesis, degradation, activity and sub-cellular localization of proteins and protein complexes. Coherent regulon activation depends on multisite phosphorylation of Whi5, SBF and MBF, catalysed by Cln-Cdk1 complexes. A discrete stochastic module calculates the probability distributions of the different phosphorylation states of the DNA-bound SBF (SBF-Whi5 complex) and MBF transcriptional activators and turns on G 1 /S transcription accordingly. The model computes percentage of the activated genes of the regulon, and dedicated ODEs compute synthesis of Cln1, Cln2, Clb5, Clb6 and Nrm1. G 1 ends when 50% of Sic1-the inhibitor the Clb-Cdk1 complexes-has left the nucleus following phosphorylation by Cln1,2,3-and Clb5,6-Cdk1 complexes.
The mathematical model-constructed using a previously reported model 47 as a stepping stone-is described in Supplementary Notes 1-10. Supplementary Tables 1-7 report initial conditions and input parameters. Supplementary Fig. 3 reports simulations of various molecular players in an average newborn daughter cell. In small elutriated cells, the simulated kinetics of translocation of Whi5 from the nucleus-an emergent property of the network-shows close agreement with experimental findings (Fig. 1b) 29 .
Multisite Whi5 phosphorylation controls the G 1 /S transition. A convincing test of the usefulness of a mathematical model describing a complex biological process is given by its ability, expressed through simulation analysis, to quantitatively account for different experimental data sets and to offer new insight on how the distinctive functional features of a given biological process emerge from the interactions within its underlying molecular network.
Experimental evidence 42,43 indicates that Whi5 release requires the phosphorylation by Cdk1 of a relatively small number (4) of functional phosphorylation sites within a larger pool (12) of phosphosites. We thus asked how the properties of cell cycle regulation change by altering the ratio between functional and decoy sites, fixing the total number of sites to 12. As the number of functional sites increases from 1 to 3, the coherence of Whi5 exclusion from the nucleus and of the activation of the G 1 /S regulon genes markedly increases (Fig. 1c,d), as indicated by their Hill coefficients, whose values more than double ( Fig. 1e; Supplementary Fig. 4; Supplementary Movie 1). The length of the T 1 period (Fig. 1f) substantially increases. The 3/12, 4/12 and 5/12 configurations had very similar parameters. A further increase in the number of functional phosphorylation sites (up to the 12/12 configuration) does not substantially modify the Hill coefficient nor further lengthens T 1 (Fig. 1e,f). The fraction of activated genes for different phosphorylation configurations of Whi5 are reported in Supplementary Fig. 5, where the ratio of functional/total phosphorylation sites is kept fixed at 1/3. These outputs are unaffected by the specific pattern of Cln3 accumulation ( Supplementary Figs 1-2).
Another functional property of the G 1 /S transition, the significant reduction in cell size variability observed at the entrance into S phase 13 , appears affected by Whi5 multisite phosphorylation. Simulation of the G 1 /S transition for a cohort of newborn daughter cells with the P(0) size distribution reported in Fig. 1g was performed with a varying number of functional Whi5 sites to obtain the distribution of the critical cell size P s (Fig. 1h, case 4/12). Starting with a coefficient of variation (CV) of 15% at P(0), P s variability decreases by increasing the number of functional sites of Whi5 from 1 to 4 and remains constant thereafter (see Fig. 1i and Supplementary Fig. 6). These findings, together with previous genetic 42 and structural evidences 43 , indicate that the 4/12 configuration for Whi5 phosphorylation is likely to be the wild-type one.
Choice of cell fate determined by mating factor. G 1 arrest is a major pre-requisite for conjugation between haploid yeast cells of opposite mating types. Far1, the Cln-Cdk1 inhibitor presented in the model in Fig. 1a, plays a key role in this arrest. To describe cell mating, we added a module (Fig. 2a) that recapitulates the afactor-dependent cell cycle-related events by abruptly changing the relevant parameters when the mating factor is added, following a recently described mathematical formalism 5,48 (see Supplementary Note 11 and Supplementary Table 8 for parameter definitions).
The mating/mitosis switch is very sensitive to the level of nuclear Whi5 at the time of a-factor addition: when the fraction of removed Whi5 (d) is low (around 10-30%), there is a very high probability of cell cycle arrest, while when it is high (over 50%), the probability of cell cycle arrest drops close to zero 5 . Simulation with the standard (4/12) configuration for Whi5 phosphorylation and standard values for parameters neatly predicts the experimental behaviour (Fig. 2b), since the addition of a-factor when d ¼ 0.1 yielded cell cycle arrest, whereas pheromone addition when d ¼ 0.5 allows the cell cycle to proceed unperturbed towards mitosis.
Populations of daughter cells were exposed to virtual a-factor addition at different times after birth, their G 1 /S transition was simulated with the extended model (4/12 Whi5 configuration) and their cell fate was recorded as cell cycle arrest or entrance into S phase. The pattern of arrested cells as a function of parameter d is a good predictor of a-factor-induced arrest, with little overlap between the two subpopulations of arrested and committed cells as shown for experimental data (Fig. 2c,d). The insets report the probability of arrest, computed from the distribution of arrested cells, the shaded region indicating the 95% confidence intervals based on 10,000 bootstrapping iterations 5 . Our simulations refer to daughter cells only: almost 40% of them result committed, in good agreement with the 40% experimental value 5 , but lower than the value observed in the experimental mixed (daughter/parent) population ( Fig. 2d right, see also Fig. 2f in ref. 5). The pattern of arrest for different Whi5 phosphorylation configurations was very similar, with a mild increasing trend for the Hill coefficient and the median point when the number of functional sites rises from 1/12 to 4/12. (Fig. 2e; Supplementary Fig. 7). Simulation results further indicate that the parameter d is a better predictor of a-factor-induced arrest than the normalized volume or the time spent in G 1 , in agreement with experimental findings 5 (Fig. 2c,d,f-i).
The mating response presents hysteresis: cells pre-treated with a brief pulse of saturating a-factor arrest cell cycle when treated with a suboptimal a-factor concentration, which alone is normally ineffective 48 . Hysteresis was reproduced by our model (Fig. 2j, 4/12 configuration). Consistently with experimental results, CLN3 deletion reinforced hysteresis, in keeping with the notion that Cln3-Cdk1 is the major kinase involved in the recovery from a-factor treatment and thus a privileged target of Far1 inhibitory action. Hysteresis was also detected in simulations performed with the 1/12 configuration for Whi5 phosphorylation, but full arrest required a higher concentration of virtual a-factor (note different abscissa scale in Fig, 2j,k. See also Supplementary  Fig. 8). In summary, the extended model captures the major distinctive features of a-factor-induced cell cycle arrest.
Phosphorylations Inhibition Complex formation/dissociation diffusion, degradation Gene transcription/protein production P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P Ydj1 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11372 G 1 duration in daughters born with different size. Newborn daughter cells show large size variability. Smaller cells delay entry into S phase, suggesting a deterministic requirement for a critical cell size in the execution of the G 1 /S transition 13 . Statistical analysis suggested that an efficient sizer operates in smaller daughter cells, whereas in larger ones a timer appears to be active, size-independent molecular noise providing the largest quantitative contribution to G 1 variability 13 .
To account for these data, we simulated a population of daughter cells (Fig. 3a, red dots) with different birth size P(0) and found that it superimposes quite well with experimental findings (blue dots) in a plot of lT G1 (where l is the rate of exponential growth) versus ln[P(0)] 13 . The blue line in Fig. 3b (standard value curve) plots the lT G1 value obtained from simulations of about 20 single cells sharing the same model parameters, and having a different P(0), selected at even intervals in the same range covered by the population of newborn cells shown in Fig. 3a. The red line corresponds to the average value of 50 simulations for each P(0) value, with parameters allowed to vary according to a log-normal distribution, with a 12% coefficient of variation. The 75% quantile region (grey area), indicates that variability is larger for smaller cells and that larger cells converge to a constant lT G1 value. The large variability in lT G1 of smaller cells is not directly due to cell size, but is strongly related to the low amount of Cln3 present. In fact a cell population with the same initial cell size distribution as wild-type cells and a 6-fold higher Cln3 level -that allows more efficient phosphorylation of Whi5-shows a decrease in G 1 length and reduced variability (Fig. 3c,d).
The standard value curve was very well fitted by a hyperbolic curve (see Supplementary Notes 14-15; red line in Fig. 3e) rather than by the previously proposed broken-line plot 13,29 . The equation of the best-fitting hyperbola (Supplementary Note 15, eq.(C.5.1.1)) allows to estimate lT G1 from any P(0) value, providing an input-output relationship (resembling numerical results achievable by simulations from the model in Fig. 1a), which cannot be derived through an analytical approach.
In our model, the T 2 period, which begins when 50% of Whi5 has been exported from the nucleus, terminates when 50% of nuclear Sic1 is exported from the nucleus. We analysed the impact of CLN2 dosage-simulated by the appropriate alteration of SBF/MBF dependent synthesis-on T 2 length. Using the same approach utilized in Fig. 3b, we show that T 2 length and its variability are insensitive to CLN2 dosage (Fig. 3f). Di Talia et al. 13 experimentally decomposed the overall noise in G 1 length into a size-dependent and a size-independent part. Table 1 compares the simulated noise (see Supplementary Note 16) with the experimental one. The overall and decomposed values for G 1 noise compared quite well for wild-type cells and Cln3overexpressing cells, with the general tendency of the simulated cells to be less noisy than their real counterparts. Our model rationalizes the observed size variability and newly links it quantitatively to the variability of the G 1 phase duration.
The observation that increasing the number of SBF-binding sites in a cell delays budding and increases the critical cell size, as if a genomic binding site participates to the sensing mechanism, was considered connected to the relation between critical cell size and ploidy 28 . Simulation of our model in the presence of an increased number of SBF-binding sites (and hence of bound Whi5) reproduces the experimentally observed increase in cell size (red and black lines, Fig. 3g). Consistently with the notion that Cln3-Cdk1 activity is the major limiting factor in the G 1 /S transition, an increase in Cln3 quenched the effects of the increased availability of SBF-binding sites (compare the red and black lines in Fig. 3h), as experimentally reported 28 . Thus, multisite phosphorylation of SBF-bound Whi5 fully explains the findings that were interpreted to indicate titration of SBF-binding sites by Cln3 (ref. 28).
Linear growth rate and critical cell size. Ferrezuelo et al. 29 recently proposed a theory of critical cell size control that drastically departs from the sensing and setting theory 2,26 . They noticed that experimentally determined values of V s (the cell volume at the end of T 1 , that is, when 50% of Whi5 has been translocated out of the nucleus) are proportional to the linear growth rate in T 1 , called a. On these basis they proposed that 'the critical cell size is set at a single-cell level by linear growth rate' 29 . (see also Supplementary Fig. 9). In the following, we test whether the experimental findings reported in ref. 29 could find an alternative interpretation in our model.
Hence we considered a standard daughter cell population similar to that shown in Fig. 3a. Figure 4a reports the plot of computed aT 1 versus V i values, which have been used to construct the corresponding V s (volume at the end of timer T 1 , proportional to critical cell size) versus a plot (black dots) in Fig. 4b. Simulated data of cells growing according to an exponential growth kinetics well approximate experimental data taken from ref. 29. The linear relationship between V s and a can also be obtained by analytical approaches (Fig. 4c; see also Supplementary Notes 17-18 and Supplementary Fig. 10), starting from the best-fitting hyperbola shown in Fig. 3e. In conclusion, the observed linear relation between V s and a does not have as unique possible interpretation the one proposed in ref. 29 and hence it is not necessary to hold that it is the linear growth rate that determines the critical cell size.
We next simulated the behaviour of two mutants with deletions in CLN3 and YDJ1 (Fig. 4d,e, respectively). Inactivation of Cln3 yields larger cells. The V s versus a plot obtained by a simulated cln3D population is in good agreement with the experimental one 29 (compare Fig. 4d to experimental data in Fig. 4d of ref. 29). Experimental deletion of the YDJ1 gene results in a very disperse V s versus a plot that suggests the existence of two subpopulations with different growth rates contained within the red and blue ellipses in Fig. 4f. Gene-dosage experiments are The simulations were carried out by sampling the initial P(0) from a log-normal distribution (average value 1.6 Â 10 10 aa, CV 15%) and fixing R(0) ¼ r Â P(0). In each simulation, all the parameters were allowed to vary (log-normal distribution) with 5% CV over the average values (see Supplementary Table 10 for details). (i) Diagram of the CVs for the distributions of different phosphorylation schemes, same rules to assign parameters and initial conditions as in g,h (distributions reported in Supplementary Fig. 6).  Table 9 for the implementation details.
often accompanied by changes in growth rate and large transcriptional reprogramming [49][50][51] . Ydj1 interacts either physically or genetically with almost 500 unique genes or gene products and deletion of its encoding gene induces slow growth and reduced fitness and lifespan (www.yeastgenome.org; Supplementary Fig. 11). Accordingly, simulation of two subpopulations differing in growth rate is required to obtain a good fitting to experimental data ( Fig. 4e compared with Fig. 4f).
Parameter-sensitivity analysis. The quantifiable outputs of a multiparameter model that aims-as our model does-to describe complex biochemical networks may not contain enough information to ensure assignment of a univocal value to each parameter in the parameter space. Application of a formal Bayesian approach 52 is hindered (i) by the difficulties in achieving any a priori characterization of the probability density of parameters like time constants, kinetic coefficients and thresholds and (ii) by the observation that systems biology models share the sloppiness property 53,54 , with the term 'sloppiness' referring to 'the highly anisotropic structure of parameter space, wherein the behaviour of models is highly sensitive to variation along a few stiff directions (combinations of model parameters) and more or less insensitive to variation along a large number of 'sloppy' directions' 55 .  Table 10 for the implementation details). A CV value of 20% was chosen for the initial protein distribution, given that the experimental sample to be compared with simulation results includes unusually small cells 13 . (b) Correlation between lT G1 and P(0). Twenty three cells differing only in their initial protein content were simulated (blue line). The red line refers to the average results of 50 separate cells in each initial condition, parameters allowed to vary as described in Supplementary Table 10. The 75% quantile region is in grey. (c,d) G 1 phase length heterogeneity for wild-type cells (c) and cells overexpressing Cln3 (6 Â over wild-type, d). (e) Best-fitting hyperbola for the standard value curve of lT G1 (blue line, b) versus the initial protein content P(0). (f) Correlation between lT 2 and the Cln2 production rate. Twenty three cells differing only in their Cln2 production rate were simulated (blue line). The red line refers to the average results of 50 separate cells in each initial condition, with parameters that were allowed to vary as described in Supplementary Using sensitivity analysis, we assessed how parameter variation affects four major outputs of the model: the length of Timer T 1 and T 2 , the critical size P s and the Hill coefficient (N) of the best Hill function fitting the G 1 /S regulon activation curve. Alteration of only two parameters (k 20 and y 20 , which affect Whi5 phosphorylation) impacts on all four tested outputs (Fig. 5a,b). In e, clusters refer to cells simulated according to growth parameters providing an exponential growth rate of 0.0031 min À 1 (red cluster) or 0.0063 min À 1 (blue cluster). In f, the red cluster refers to cells of the experimental sample reported in ref. 29 that fall within the grid of Supplementary Fig. 10C. The blue cluster comprises the remaining cells. WT, wild type. Decomposition of G1 variability in terms of deterministic size control and a residual that is attributable to molecular noise under the assumption that lTG1 can be decomposed into the sum of a suitable deterministic function of the initial protein content f(P(0)) plus a stochastic variable. The chosen function for f(P(0)) is the hyperbola described in Fig. 3e. The experimental data are from ref. 13. The Timer lengths and the sharpness of the G 1 /S regulon activation curve are nontrivially sensitive to a small set of parameters, in keeping with the sloppiness hypothesis (see also Supplementary Figs 12-15). Only the length of T 1 , but not of P s and T 2 , shows some correlation with changes in the Hill coefficient N (Fig. 5c; Supplementary Fig. 16), notably so for parameters related to the Cln3, Ydj1, Far1 machinery and for some parameters related to G 1 /S regulon activation (Fig. 5c, black and red points, respectively). We then analysed the impact of the functional (4) to decoy (8) configuration of Whi5 on the coherence of the G 1 /S transcription (quantified by N, the Hill coefficient for transcriptional activation) and on the timing delay, monitored as duration of T 1 period. The connection between N and T 1 is shown in Fig. 5d. Increasing the total number of sites (from 4 to 24) with a variable number of functional sites, there is a sizable increase of both N and T 1 , till functional sites reach 4. Above this value, saturation is obtained. Increasing the total number of sites from 3 to 24 and keeping a constant ratio (1/3 as in 4/12 Whi5) between functional and total sites, the increase in total sites correlates with an increase of both T 1 length and Hill coefficient value.
The Whi5 4E mutant anticipates entrance into S phase. To be fully accepted, a new model is required both to account for large sets of existing data and to make predictions qualitatively distinct from its predecessors, which need to be validated by new, independent sets of experimental findings.
Glutamate may mimic phosphoserine or phosphothreonine 56 . A Whi5 4E protein in which the four functional phosphosites have been mutated to glutamate should dissociate from SBF nearly as well as Whi5 in which the functional phosphosites have been phosphorylated by Cdk1. Since (pseudo)phosphorylation is constitutive and unregulated, yeast cells harbouring the mutated Whi5 4E protein should exhibit a small (Whi) phenotype similar to the whi5D mutant. Our model further predicts that Whi5 4E G 1 cells  would anticipate the activation of the G 1 /S regulon, entering the S (budded) phase at a much smaller size than their wild-type counterparts (Fig. 6c,e-j). Since our model does not include events following the G 1 /S transition that down-regulate the G 1 /S regulon and lead to cell division, the fraction of budded cells remains high and transcription of genes of the G 1 /S phase steadily increases, except for CLB6, whose expression is downregulated by Nrm1. We constructed a strain expressing the mutated Whi5 4E protein. The transcript levels of WHI5 4E and WHI5 were undistinguishable (Supplementary Fig. 17). Cell volume and cell protein distributions (Fig. 6a,b) confirmed that the whi5 4E strain maintains a Whi phenotype, being much closer to the whi5D mutant than to wild type. When reinoculated in fresh glucose medium whi5 4E small G 1 cells (isolated by elutriation) showed an anticipated entry in the budded, S phase (Fig. 6d) and an early activation of G 1 /S regulon (Fig. 6k-p) relative to wild-type cells. In keeping with model predictions, both cell cycle entry and gene expression kinetics of the whi5 4E mutant are intermediate between those of the wild type and of the whi5D mutant, being more similar to the latter.

Hill coefficient
Taken together these results validate the central notion of our model: dissociation of Whi5 from SBF is promoted by phosphorylation of its functional phosphoacceptor sites.

Discussion
The cell-fate decision given by the mitosis/mating switch and the connected critical cell size control-the gate-keeper for the commitment of G 1 cells to cell division-has been for decades one of the most studied examples of complex biological functions. Genetic and biochemical analyses have identified many of the involved players without being able to propose a fully satisfactory, comprehensive molecular model. By utilizing a systems biology approach which integrates mathematical modelling and computational analysis with experimental validation, this paper proposes a precise, quantitative, predictive and unifying mechanism of the molecular events involved in this cell fate decision, accounting also for a large set of experimental data, which have been used to propose different and even conflicting conceptual models for critical cell size control 13,26,29 .
The novelty proposed by our paper is given by the multisite phosphorylation of Whi5 by Cln3-Cdk1 (first) and Cln1,2-Cdk1 (later). Our model assumes that Cln3 concentration is constant during G 1 . Our simulations show that the specific pattern of Cln3 accumulation has negligible effects on the outputs of the model ( Supplementary Figs 1 and 2), which are instead affected by Cln3 average content. Nuclear Cln3 import is favored by the Yjd1 chaperone 27,29 , the amount of nuclear Cln3 being larger in big cells than in small ones. We show that Cln3-Cdk1 is rate limiting for Whi5 phosphorylation, except than in very large cells. Cdk phosphoacceptor sites inside Whi5 include both decoy and functional sites 42,43 . Decoy sites fruitlessly engage Cln3-Cdk1 complexes, delaying phosphorylation of the four functional sites. Full phosphorylation of Whi5 (or of the Swi6 subunit of SBF) is required to release Whi5 from SBF. This mechanism originates synchronous dissociation of Whi5 from SBF at the hundreds of different promoters of the G 1 regulon present in each cell, ensuring the coherence of the G 1 /S regulon transcription, essential for an ordered cell cycle progression and fully accounting for the observed role of Whi5 as the integrator of the signalling generated by Cln3 (ref. 57).
As an experimental validation, we report that Whi5 4E -a mutant protein mimicking constitutive phosphorylation of the four functional residues-causes anticipated activation of the G 1 /S regulon and S phase entrance at a reduced size, close to that of the whi5D strain. Pseudophosphorylation of the functional residues likely induces a change in folding of Whi5, an intrinsically disordered protein 43 , that would act as a regulatory switch 58 , by reducing the binding affinity of Whi5 to SBF.
As anticipated in the Introduction, our paper does not cover the modulation of critical cell size by nutrients, but we may offer a working hypothesis, elaborated following the idea 2,26 that the mechanism that sets the critical cell size in response to nutrients adds molecular components to the basic mechanism (the dynamic interaction of Cln-Cdk1 and Whi5) by which cells gauge their size. Ribosome biosynthesis, regulated by cAMP/PKA and Tor pathways (sensors of nutrient availability 22,26,59,60 ), has been identified as the primary actor in the setting of the critical sell size by nutrients 26,61,62 . Mutations in ribosome biogenesis strongly affect Whi5 nuclear retention 63 . Whi5 has many phosphorylation sites, specific not only for Cdk1, but also for other pathways 43 , and may therefore integrate contrasting inputs, as reported in osmostress that may affect Whi5 binding activity 64 . On the other hand, the Cln3 level is very sensitive to changes in ribosome activity 31 .Our model (in which the exponential growth rate is dependent on the rate of protein synthesis) predicts that critical cell size becomes larger as a cell increases its rate of protein synthesis ( Supplementary Fig. 18). Thus, nutrients, through signalling, may affect ribosome biogenesis that, directly or indirectly, affects the Cln3/Whi5 interplay. The clarification of the proposed mechanism of critical cell size setting by nutrients may be brought on by new, quantitative, integrated, dynamic investigations performed taking into account previously collected and newly generated evidence.
A DNA fragment encoding the Whi5 4E CDK mutant was synthesized de novo by Eurofins (www.eurofins.com) and subcloned into the YIplac211 integrative plasmid 66 under the WHI5-545 pr native promoter 42 , yielding the construct YIplac211-545 pr -WHI5 4E . Single-copy genomic integration of the construct at the URA3 locus was verified by quantitative PCR.
Growth of cultures was monitored as increase in cell number using a Coulter Counter model Z2 (Coulter Electronics, Inc.). The fraction of budded cells was scored by direct microscopic observation on at least 400 cells, fixed in 3.6% formaldehyde and mildly sonicated. Cell size analysis was performed using a Coulter Z2 Particle Cell Analyzer (Beckman-Coulter).
For synchronization studies, yeast cells were grown till late exponential phase in 2 l of csm/YNB medium supplemented with 2% raffinose 38 . Small G 1 phase cells were isolated by centrifugal elutriation using a 40-ml chamber elutriator (Beckman Coulter) and released into fresh 2% glucose medium at 5 Â 10 6 cells per ml. Samples were collected at appropriate intervals for determination of budding index, cell volume, DNA and protein content by flow cytometry and for quantitative PCR with reverse transcription (qRT-PCR).
At least 2 Â 10 7 cells were collected and fixed in 70% ethanol before the cytofluorimetric analysis. For RNA staining, cells were washed once with cold PBS (3.3 mM NaH 2 PO 4 , 6.7 mM Na 2 HPO 4 , 127 mM NaCl. 0.2 mM EDTA, pH 7.2), resuspended in 1 ml of Propidium Iodide staining solution (0.046 mM propidium iodide in 0.05 M Tris-HCl, pH ¼ 7.7; 15 mM MgCl 2 ) and incubated for 30 min on ice/dark. To obtain protein distribution, cells were stained with fluorescein isothyocyanate solution (50 mg ml À 1 FITC in 0.5M NaHCO 3 ) for 30 min on ice/ dark and washed three times in PBS before the analysis. For DNA staining, cells were washed once in PBS, resuspended in 1 ml of PBS containing 1 mg ml À 1 RNAse and incubated at least 12 h at 37°C. Cells were then washed once with PBS, resuspended in 1 ml of Propidium Iodide solution and incubated for 30 min on ice/ dark. All cell suspensions were sonicated 30 s before the analysis, which was performed with a FACSCalibur (Becton Dickinson) instrument equipped with a 488 nm Ion-Argon laser. Sample flow rate during analysis did not exceed 1,000 cells per s. Typically, 50,000 cells were analysed for each sample.
For RNA extraction, B2 Â 10 8 cells were collected by filtration, washed with 5% Trichloroacetic acid solution and rapidly freezed at À 80°C. Cells were resuspended in 500 ml of LETS buffer (200 mM LiCl, 20 mM EDTA, 20 mM Tris pH ¼ 7.4, 0.5% SDS), plus 300 ml of phenol:chloroform:isoamyl-alcohol solution NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11372 ARTICLE (PCI) and lysed by glass beads. Two PCI extractions were performed on the recovered aqueous phase. RNA was subsequently isolated by ethanol precipitation in the of presence LiCl. Five hundred micrograms of RNA were treated with 9 units of DNase I (RNase-free, Qiagen) for 1 h at 37°C. Total RNA was then purified with the RNeasy kit (QIagen) and 0.5 mg were used for cDNA synthesis using the iScript cDNA Synthesis Kit (BIO-RAD) according to the manufacturer's instructions. qRT-PCR reactions (final volume 15 ml) were performed in 48-well reaction microplates suitable for the MiniOpticon detection system (BIO-RAD) using the SsoFast EvaGreen Supermix (BIO-RAD). Primers used for the assay were designed with the Beacon Designer software (PREMIER Biosoft). Data obtained were analysed with the CFX Manager software (BIO-RAD) and normalized by including TAF10 and UBC6 as reference genes in each plate. At least two biological replicates were analysed in triplicate for each gene 49 .
Model construction. The proposed mathematical model links a basic set of growth equations to the ones related to the molecular player dynamics promoting the G 1 /S transition. The growth equations are taken from ref. 30, and deal with the ribosome and protein content. The molecular model equations take into account the following classes of events: (i) Regulation of Cln3-Cdk1 activity, that is, Cln3 synthesis (in the cytoplasm), its subsequent retention in the endoplasmic reticulum and its transport in the nucleus by the Ydj1 chaperone as well as Cln3 inhibition by Far1. (ii) Transcriptional activation of the G 1 /S regulon, first by Cln3-Cdk1 kinase activity, which includes a computationally effective description of regulon activation-including Cln1,2-mediated positive feedback-and multisite phosphorylation of Whi5. (iii) Regulation of Clb5,6-Cdk1 activity-that triggers entrance into S phase-by binding of Sic1 whose nuclear level and inhibitory activity are regulated by Cln1,2-Cdk1 AND by Clb5,6-Cdk1 phosphorylation.
The link between the growth equations and the molecular model is given by the synthesis of the upstream cyclin, Cln3, since its total amount is constrained to be proportional to the overall protein content.
The model-which takes into account cytoplasm, endoplasmic reticulum and nucleus sub-cellular compartments-is a hybrid deterministic-stochastic model. The amount and sub-cellular localization of the molecular players involved in the activation of the regulon is modelled by a set of ODE. Regulon activation depends on a correct sequence of multisite phosphorylation steps-catalysed by Cln1,2,3-Cdk1 complexes. A stochastic model that calculates probability distributions of the different phosphorylation states describes the phosphorylation states of DNAbound SBF and MBF transcriptional activators, as well as of the transcriptional inhibitor Whi5. Due to the abundance of Cdk1 and to the fast dynamics of Cdk1 binding to its cyclin targets, the molecular model does not explicitly consider the binding dynamics for Cdk1 and Cln and Clb cyclins The model accounts for the percentage of the activated genes of the G 1 /S regulon (235 genes), as well as for the evolution of five of the expressed genes, Cln1, Cln2, Cln5, Clb6 and Nrm1.
From a numerical viewpoint, the model sums up a set of ordinary differential equation, whose order depends of the chosen phosphorylation configuration. For instance, in case of wild-type cells (that is, 4 functional phosphorylation sites for Swi6, 4 functional/12 total phosphorylation sites for Whi5), we have a set of 264 ordinary differential equations. To simulate these nonlinear dynamics requires the setting of about 100 model parameters/initial conditions. Essential biological evidence, supporting references, and an outline of the reconstruction and computational representation of the G 1 /S transition model are detailed in Supplementary Notes 1-10. Parameter setting and in silico experimental procedures are reported in Supplementary Note 12 and 13, respectively. The parameters are presented in Supplementary Tables 1-7.
Simulation of the G 1 /S model. Single-cell simulations are carried out in Matlab, by means of the built-in ODE Matlab solver 'ode15.s'. See Supplementary Table 9 for specific growth model parameters in single-cell simulations. From an implementation viewpoint, all the model parameters are sampled from random distributions before running of a population simulation, thus allowing to simulate different cells for different runs, all belonging to the same population. These parameters include the order of proper expression and the times of weak expression for the molecular players Cln1, Cln2, Cln5, Clb6 and Nrm1 (Supplementary Note 9).
The model provides in output ribosome and protein content as well as all the involved molecular player time courses according to which the following relevant outputs are computed: (i) length of timer T 1 , when cytoplasmic Whi5 exceeds 50% of total Whi5. (ii) length of Timer T 2 , when cytoplasmic Sic1 exceeds 50% of total Sic1. (iii) critical size P s , the protein content at the onset of the budded phase.
A further output is given by the time course of the percentage of activated genes. The Hill function best-fitting such a percentage provides a fourth relevant output: the Hill coefficient aiming at measuring the coherence of the G 1 /S regulon activation.
Code availability. The MatLab file is available at http://dx.doi.org/10.5281/ zenodo.46955. This link allows to obtain all the software, in a zip file, necessary to simulate our model plus a readme file containing the procedure to run the code. All the model parameters are defined in the d1_one_cell.m file. The MatLab version used was the R2014b.
Populations. Cell populations were obtained according to the following rules: (1) The initial protein content P(0) of each cell was drawn from a log-normal distribution of a given average value (Z P ) and s.d. (s P ). Different values of Z P allowed to build populations of smaller or larger cells, and different values of s P allowed to build populations with a larger or smaller cell protein dispersion. (2) The initial ribosome content R(0) was straightforwardly set from P(0) using the equation: R(0) ¼ r Â P(0).  Table 5). (4) All other model parameters reported in Supplementary Tables 1-4 and 6-7 were drawn from log-normal distributions of chosen average values, sharing the same CV. (5) Parameters specific for each simulation of cell population are reported in Supplementary Table 10.
Model complexity and simulation times are discussed in Supplementary Note 24, and Supplementary Fig. 19.
Sensitivity analysis. For sensitivity analysis, tested parameters were grouped in four sets: Cln3 production and nuclear import (all parameters in Supplementary  Table 3); G 1 /S regulon activation (all parameters in Supplementary Table 4 except for Swi6 tot , Swi4 tot , Mbp1 tot and Whi5 tot ); Cln/Clb function (all parameters in Supplementary Table 6); Sic1 function (all parameters in Supplementary Table 7 except for Sic1 tot ). Taking the standard value as 1, each parameter was either increased (up to 81-fold the standard value, step 3) or decreased (down to 1/81-fold the standard value, step 1/3). The impact of each of these changes in parameter value was tested on four significant outputs: N, the Hill coefficient for G 1 /S regulon activation curve, the length of the T 1 and T 2 period and the critical size P s . The sensitivity analysis was done by varying each parameter individually, associating it to the run of a single-cell deterministic simulation, with the proper order of activation and the weak activation times fixed to their average values, i.e. b Cln1 , b Cln2 , b Nrm1 , b Clb5 , b Clb6 (Supplementary Table 5). Results and further details are reported in Supplementary Figs 12-15 and Supplementary Notes 19-23.