1. Deterministic SCM

Proteins regulate one another by controlling their abundances through rates of synthesis and degradation and their activities through post-translational modification (e.g., phosphorylation and dephosphorylation), and by associating into multisubunit complexes. These three classes of reactions often proceed on different time scales. Synthesis and degradation cause rather slow changes in the total amount of a protein (time scale > 10 min). Phosphorylation and dephosphorylation cause faster changes in protein activity (time scale = 1–10 min). Rapid association and dissociation of protein complexes bring the complexes and subunits into equilibrium on a short time scale (1 min or less). In cases where this separation of time scales is known or suspected to be the case, we can partition the components of a PRN into three classes:

1. Class-1 variables track the total amounts of proteins, which evolve rather slowly in time due to synthesis and degradation.
2. Class-2 variables track the activities of proteins, which change faster due to covalent modifications.
3. Class-3 variables track protein complexes, which turn over rapidly by association and dissociation of subunits.

In principle, variables of these three classes can represent different forms of the same protein. For example, we can use a class-1 variable to represent the total amount of a protein, while using a class-2 variable to represent the fraction of the total protein that is phosphorylated, and a class-3 variable to represent the fraction of the total protein that is bound to a stoichiometric inhibitor.

In our formalism, class-1 variables (X_i) are governed by pseudo-linear differential equations for protein synthesis and degradation (1)

The ODE is linear in X_i, the number of molecules of species X_i, but the rates of synthesis and degradation are functions of variables Y_j that may belong to any of the three classes. It is often possible to use linear functions for A_i and B_i, (2) (3) where α_i0 and β_i0 are basal rates of synthesis and degradation, and α_ij ·Y_j and β_ij ·Y_j are rates regulated by transcription factors and proteolytic enzymes, respectively. (In this case, the biochemical rate parameters α_i0, β_i0, α_ij, and β_ij are all positive constants.) In other cases—especially for transcription factors that inhibit gene expression—nonlinear functions for A_i and B_i may be required.

Class-2 variables are governed by nonlinear ODEs of the form (4) where Y_j represents the activity of protein Y_j (e.g., the phosphorylated or the active form of Y_j), Y_j,_T is the total number of molecules of protein Y_j, γ_j determines the time scale of the reaction, and H(x) = 1/(1 + e^−x) is the sigmoidal function illustrated in S1 Fig. (H is a hyperbolic tangent function shifted along the y-axis. In population biology it is known as the “logistic” function. We refer to H as the “soft-Heaviside” function, because we use it to replace the step-like Heaviside function used in the piecewise-linear models of Glass, Kauffman and others.) In the soft-Heaviside function, W_j describes the net influence of all components in the network on the component Y_j: (5)

In Eq 5, ω_jk and ω_jl are weights (always positive values) that describe the influences of variables Y_k and Y_l on the variable Y_j. K represents all variables that have positive influences on the variable Y_j, and L represents all variables that have negative influences on the variable Y_j. Y_k and Y_l can be variables of any of the three classes of species. The background influence, ω_j0, which can be preceded by either the positive or negative sign, determines the value of the soft-Heaviside function when protein Y_j is receiving no inputs from the other proteins in the network. The parameter σ_j controls the steepness of the soft-Heaviside function; see S1 Fig. In principle, the value of σ_j could be absorbed into the values of the ω’s, but we prefer to treat σ_j as a separate parameter and to think of the ω’s as relative interaction strengths. That way, we can vary the steepness of the soft-Heaviside function independently of the relative interaction strengths and vice versa.

Eq 4 shows that H (σ_j·W_j) determines the steady state of the variable Y_j as a fraction of the total amount Y_j,T. If the total amount of the protein remains constant over time, Y_j,T is a constant parameter in the model. If the total amount of the protein changes in time, we can use a class-1 variable to keep track of Y_j,T while using a class-2 variable to keep track of the fraction of the protein that is in the active form.

Class-2 variables evolve to a steady state on a time scale that is proportional to γ_j⁻¹. Therefore, if γ_j is large, we can invoke the pseudo-steady state approximation for the class-2 variable: (6)

Moreover, if both γ_j and σ_j are large, then the class-2 variable, Y_j/Y_j,T, behaves like a Boolean variable, switching rapidly between 0 and 1 in response to other components in the network (S1 Fig).

Eq 4 has been used to describe interactions in PRNs many times before [12, 20–23]. We use the soft-Heaviside function because biochemical reactions (such as phosphorylation and dephosphorylation) are nonlinear in nature and often show ultrasensitive, sigmoidal responses [24–26]. Different mechanisms have been proposed to give rise to ultrasensitivity, and different types of rate equations have been used to capture this response. For transcription factors binding to gene promoter regions, Hill functions are often used to express the highly nonlinear (ultrasensitive) response of gene transcription rate to transcription factor concentration. For post-translational modifications of proteins, a commonly used mechanism is zero-order ultrasensitivity, originally proposed by Goldbeter & Koshland (GK) [25]. The GK equation describes the steady-state activity of a protein modified by two Michaelis-Menten enzymes with opposing effects (activation and inactivation). More recently, multisite phosphorylation has received considerable attention as an alternative mechanism of ultrasensitivity in PRNs [26–28]. All these mechanisms ultimately lead to sigmoidal-like responses, similar to the soft-Heaviside function (S1 Fig). Therefore, we eschew specific assumptions about the molecular mechanisms of ultrasensitivity and simply use the soft-Heaviside function as a phenomenological law that captures the ultrasensitive responses that are characteristic of PRNs.

Class-3 variables describe the rapid association and dissociation of proteins in multi-subunit complexes. For example, if proteins Y and I bind together strongly in an inactive complex, then the number of Y molecules that are free (not bound to I) is given by (7)

This equation is a good approximation if both association and dissociation rates are fast, and the equilibrium binding constant is large. The max function returns zero if the total amount of the protein, Y_T, is less than the total amount of its inhibitor, I_T, meaning that all of protein Y is sequestered in the complex. On the other hand, if Y_T > I_T, then the free form of the protein, Y, is simply the excess of Y_T over I_T.

If an inhibitor associates with two different proteins with comparable association and dissociation rates, the amounts of the free proteins can be approximated by (8) where I_T, Y_1,T and Y_2,T are the total amounts of the proteins, and I, Y₁ and Y₂ are the amounts of the free forms of the proteins. In Eq 8, we assume that the inhibitor is distributed equally between Y₁ and Y₂, according to the availability of the two proteins. If I binds much more strongly to Y₁ than to Y₂, then Eq 8 should be replaced by (8′)

In the following sections, we show how to use these three classes of variables as building blocks to construct an SCM of the molecular network controlling the cell division cycle in budding yeast. Because the full model is very complex, we approach it in a series of simpler steps.

2. The Start transition in budding yeast

Start is an event in G₁ phase of the budding yeast cell cycle when a cell commits to a new round of DNA synthesis and mitosis. A crucial step of the Start transition is the translocation of Whi5, a stoichiometric inhibitor of SBF and MBF transcription factors, from nucleus to cytoplasm [29, 30]. SBF and MBF control the expression of CLN2 and CLB5 genes, which encode “cyclin” proteins Cln2 and Clb5, respectively. Cln2 and Clb5 bind to kinase subunits (Cdc28) to form heterodimers with “cyclin-dependent kinase” (CDK) activity. CDK activity generated at Start triggers initiation of DNA synthesis and bud emergence. Because kinase subunits are in excess over cyclin partners [31], CDK activity is determined solely by the abundance of cyclin proteins. For simplicity in illustrating the SCM approach for the Start transition, we combine Cln2- and Clb5-dependent kinase activities into a single variable, called ClbS. We also treat SBF and MBF as a single variable, called SBF.

During normal cell cycle progression in budding yeast, the cell needs to grow sufficiently large to execute Start [32, 33]. The major players involved in “size control” of Start are Cln3 and Whi5. Whi5 prevents the Start transition by binding to and inhibiting SBF, and Cln3 promotes Start by phosphorylating and inactivating Whi5 [29,30]. The accumulation of Cln3 in G₁ phase seems to depend on cell growth [34], and recent evidence suggests that Whi5 concentration is diluted out by cell growth [35]. As the cell grows, Cln3-dependent kinase phosphorylates Whi5, resulting in translocation of Whi5 from nucleus to cytoplasm and the release of its inhibition on SBF. Free SBF promotes the synthesis of ClbS, which stimulates its own expression by further phosphorylating Whi5. This positive feedback loop is thought to enforce the irreversible commitment of cells to the Start transition [36]. A schematic diagram illustrating the molecular basis of the Start transition is shown in Fig 1A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. The Start transition.

(A) Schematic diagram of the Start transition in budding yeast. In early G₁, SBF is inactivated by its stoichiometric inhibitor, Whi5. As cell size increases, Cln3 accumulates and begins to phosphorylate Whi5. Phosphorylated Whi5 loses its ability to bind to SBF. As a result, SBF is free and promotes the production of ClbS (Cln2 and Clb5). ClbS exerts positive feedback on its own accumulation by further phosphorylating Whi5. The activation of SBF correlates with the onset of the Start transition. Subsequent accumulation of ClbS promotes both bud emergence and the G₁/S transition. (B) Wiring diagram of the MultiP model. The first three forms of Whi5 (Whi5, Whi5P₁, and Whi5P₂) bind to SBF and inhibit its ability to activate the synthesis of ClbS. The higher phosphorylated forms are inactive and do not bind to SBF. The model also includes mRNA species for each protein component. (C) Wiring diagram of the standard component model. The 10 distinct forms of Whi5 in the MultiP model are replaced by two forms of Whi5 (active and inactive). For panels B and C, solid lines indicate chemical reactions (synthesis and degradation, phosphorylation and dephosphorylation, association and dissociation) and dashed lines indicate activatory or inhibitory influences of components on chemical reactions.

https://doi.org/10.1371/journal.pone.0153738.g001

Before constructing an SCM of the Start transition, we first describe a multisite phosphorylation (MultiP) model that will serve as a “reference point” for judging the adequacy of the SCM.

3. A multisite phosphorylation model of the Start transition

Our MultiP model is a simplified version of a model developed by Barik et al. [27] (see Fig 1B), who used mass-action kinetics to describe all the reactions involved in the Start transition and carried out detailed stochastic simulations based on Gillespie’s algorithm [13]. Our version of the model (see the equations in S1 Text) is governed by 20 molecular species, including both proteins and mRNAs (m_n3, m_bS, m_i5, and m_hi5 are mRNAs for Cln3, ClbS, Whi5, and Hi5, respectively). The 20 molecular species participate in 50 reactions, representing synthesis, degradation, phosphorylation, dephosphorylation, association, and dissociation of the species. Mass-action rate laws are used for all reactions. Each differential equation specifies how the number of molecules of each species changes with time.

Cell size (V) is assumed to increase exponentially. The rate of synthesis of each protein is assumed to be proportional to volume V × number of mRNA molecules m encoding the protein, because we assume that the number of ribosomes per cell increases proportionally to cell size V. In this way, we ensure that the concentration (N/V) of every constitutively expressed protein remains constant as the cell grows. The only exception is Cln3, whose synthesis rate—we assume—is proportional to the square of cell size (V²), as in Barik’s model [27]. This assumption, whose consequence is that the concentration of Cln3 increases exponentially as the cell grows, introduces a “size dependence” on the Start transition, which is needed to account for many properties of cell cycle progression in budding yeast. Although this assumption is sufficient to account for the observed size-dependence of yeast cell division, it is certainly not necessary. Other explanations are possible. Indeed, Schmoller et al. [35] have recently shown that Cln3 synthesis rate increases in direct proportion to V, so that its total concentration is nearly constant during G₁ phase; and size dependence of the yeast cell cycle depends on diluting out Whi5, an inhibitor of Start, as the cell grows. This alternative will be explored in future versions of the model.

In Barik’s model [27], multisite phosphorylation of Whi5 is the source of nonlinearity necessary for the ultrasensitive Start transition. Whi5 is phosphorylated in vivo on ~10 CDK phosphorylation sites [37]. In Barik’s model Whi5 has seven phosphorylated states: Whi5, Whi5P₁, Whi5P₂, …, Whi5P₆. In the model, the sites are phosphorylated sequentially and distributively [26]. The first three forms bind rapidly and strongly to SBF; the higher phosphorylated states (Whi5P₃, …, Whi5P₆) are inactive and unable to bind to SBF. Free SBF binds to and activates the ClbS gene (G_i + SBF ↔ G_a). Cln3 and ClbS phosphorylate Whi5 (both free and in complex with SBF), while Whi5P_i species are dephosphorylated by an unspecified phosphatase, called Hi5 (“H” for phosphatase, “i5” for Whi5).

In Barik’s model [27], Whi5:SBF complexes are called Cmp, CmpP₁, and CmpP₂. Whi5P₂ in the complex (CmpP₂) is assumed not to get further phosphorylated. For Whi5P₂ in the complex to be phosphorylated and thereby inactivated, Barik’s model supposes that Whi5P₂ must first dissociate from CmpP₂. We think this requirement is unnecessarily restrictive, so we allow the doubly-phosphorylated Whi5 in the complex (CmpP₂) to be further phosphorylated, and the triply-phosphorylated Whi5 is assumed to immediately release SBF. We also assume, in the MultiP model, that the dissociation rate of SBF:Whi5P_i complexes is negligible.

4. An SCM for the Start transition

Fig 1C is the SCM representation of the Start transition. We assign the components in this diagram to the three classes of variables proposed above.

Cln3 and ClbS are described by class-1 variables: (9) (10)

In these equations, V = cell size, which we assume increases exponentially in time, V(t) = V₀e^μt, where μ is the specific growth rate of the cell. For each protein, its rate of synthesis is the product of (rate of translation, k_s,…) × (number of mRNAs, m_…) × (cell size, V). We explicitly account for the number of mRNAs in our synthesis rates (the factors labeled m_n3 and m_bS), as the mRNA number will be important later for representing noise terms. For simplicity, we assume (in the deterministic SCM) that mRNAs are always at their steady state levels. As in the MultiP model, we assume that the synthesis rate of Cln3 is proportional to the square of cell size (V²) so that the concentration of Cln3 (number/volume) will increase exponentially as the cell grows. The activity of the gene encoding ClbS, G_a, depends on active SBF binding to the promoter region of the gene.

Since multisite phosphorylation of Whi5 gives rise to its ultrasensitive response to total CDK activity, we use a class-2 variable for the active form of Whi5 and a class-1 variable for the total amount of Whi5: (11) (12) (13) (14) where H(x) is the soft-Heaviside function defined earlier (see S1 Fig). We assume that the phosphorylation and dephosphorylation of Whi5 are independent of its binding state to SBF. In this way, we are able to use two differential equations (Eqs 11 and 13) to represent 10 distinct states of Whi5 and Whi5:SBF complexes. Here, Whi5_A (Eq 11) represents the total amount of active Whi5 both free and bound to SBF (it includes Whi5, Whi5P₁, Whi5P₂, Cmp, CmpP₁, and CmpP₂ in the MultiP model); whereas (Whi5_T –Whi5_A) represents the other four inactive forms of Whi5 (Whi5P₃ to Whi5P₆). Hi5 (the phosphatase that re-activates Whi5) is described by a class-1 variable in Eq 14.

In Eq 12, the interaction coefficients ω_p,i5, ω′_p,i5 and ω_dp,i5 are all positive numbers. The signs—positive or negative—in front of each term determine whether the interaction is activating or inhibiting (dephosphorylation or phosphorylation, in this case).

Assuming that binding between active forms of Whi5 and SBF is rapid and strong, we describe free SBF as a class-3 variable (15)

This equation indicates that free SBF is equal to the excess of the total SBF over the total active Whi5, where SBF_T (Eq 16) is represented by a class-1 variable.

(16)

Because the original model of Barik et al. [27] did not have an mRNA species for SBF, we have not included mRNA for SBF in the SCM.

Using the SCM approach, we reduce the complexity of the Start model to seven differential equations plus one algebraic equation (from 20 equations in S1 Text for the MultiP model). The parameter values used in the SCM (Table 1) are inherited, for the most part, from the parameter values in Barik et al. [27]. The parameter values assigned by Barik et al. were estimated from experimental data wherever possible. For example, protein degradation rates were calculated from protein half-life measurements in the literature, and synthesis rate constants were assigned to agree with the average number of protein molecules observed for an asynchronous population of yeast cells growing on glucose medium. For phosphorylation and dephosphorylation reaction, there are no experimentally measured rate constants. Barik et al. assigned values to these rate constants so that their model compared well with the observations of Di Talia et al. [38], and we have done the same in assigning values to γ, σ and the ω’s in Eqs 11 and 12. (S1 Table specifies four additional parameter values needed for the MultiP model.)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameter values for the standard component model of the Start transition.

https://doi.org/10.1371/journal.pone.0153738.t001

In both models, initial conditions are set as follows: Cln3 = ClbS = SBF = 0, and initial conditions for all other variables are set at their steady state levels (see Table 2 for the SCM and S2 Table for the MultiP model). In the MultiP model, at the beginning of the simulation, all Whi5 is in the unphosphorylated form and all SBF is complexed with Whi5. In the SCM, at the beginning of the simulation, all Whi5 is in the active form and all SBF is complexed with active Whi5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Initial conditions for simulations of the standard component model of the Start transition in Figs 3 and 5.

https://doi.org/10.1371/journal.pone.0153738.t002

5. Stochastic version of SCM

Deterministic models are usually sufficient to predict the average behavior of a population of cells. However, at the level of individual cells, molecular regulatory networks operate under noisy conditions. A major source of noise inside cells are fluctuations of the numbers of molecules of biochemical species undergoing random events of synthesis and degradation. These inevitable fluctuations are referred to as “molecular noise”. Under the influence of such noise, the number of molecules, N, of a biochemical species fluctuates around the value predicted by the deterministic model. For a simple synthesis-degradation process, the variance of these fluctuations is equal to the mean value, ⟨N⟩, of the number of molecules [39]. Hence, the coefficient of variation (CV = standard deviation/mean) of the number of molecules is expected to be , becoming proportionally larger as the mean number of molecules becomes smaller. For yeast cells, which carry only a few copies of some macromolecules (e.g., mRNA species), the fluctuations are potentially large.

To adapt an SCM for stochastic simulations, we add Langevin-type fluctuations to our equations for class-1 variables: (17)

In this equation, molecular noise is attributed to random fluctuations at two levels: protein and mRNA. Fluctuations at the protein level are described by the chemical Langevin approximation [40]. A_i and B_i represent the protein synthesis and degradation rates as described in Eqs 2 and 3, respectively. ζ₁ (t) and ζ₂ (t) are independent random variables, each chosen from a normal distribution N(0, 1) with mean = 0 and standard deviation = 1. Δt is the step size of the numerical integration. The last term in Eq 17 represents noise at the protein level inherited from fluctuations in mRNA molecules, and we next explain the origin of this term.

Since the number of mRNA molecules in a cell is typically much less than the numbers of protein molecules encoded by the mRNA, stochastic effects arising from mRNA fluctuations could contribute significantly to fluctuations in protein abundance. Instead of explicitly including mRNA dynamics in our SCM, as done in the Barik et al. model [27], we take an alternative approach. Pedraza and Paulsson [39] have shown that the square of the CV of protein numbers caused by mRNA fluctuations (at steady state) follows the equation (18) Where ⟨m⟩ is the average number of mRNA molecules at steady state, and τ_m and τ_p are half-lives of mRNAs and proteins, respectively. The last term of Eq 17 was derived to approximate the CV² predicted by Eq 18 under steady state conditions for the protein and mRNA levels (see S2 Text for the derivation). In Eq 17, k_dm is the rate constant for mRNA degradation; so the term B_i /(B_i + k_dm) is equal to τ_m /(τ_m + τ_p) in Eq 18.

There is still a problem with the mRNA noise term in Eq 17. When ⟨m_i⟩ is small (e.g., the case of a low level of the transcription factor SBF in Eq 10), the Langevin approximation breaks down and the noise term becomes very large. To avoid this, we replace ⟨m_i⟩ with ⟨m_i⟩ + ⟨m_min,i⟩ in our simulation, where ⟨m_min,i⟩ is assumed to be a minimum number of mRNA molecules always present in the cell. Eq 17 then becomes (19)

We treat ⟨m_min,i⟩ as a parameter, and its value may vary from one mRNA species to another.

We do not incorporate protein and mRNA noise terms into variables of classes 2 and 3. We assume that protein phosphorylation and dephosphorylation reactions and the association and dissociation of protein complexes occur fast enough to bring any fluctuations quickly back to the average dynamics.

To create a stochastic SCM for the budding yeast Start transition, we modify the equations of class-1 variables by adding Langevin-type noise terms, as in Eq 19. For example, the stochastic equation for Cln3 is: (20) where A_n3 and B_n3 represent the synthesis and degradation rates of Cln3 protein, and ⟨m_n3⟩ is the average number of CLN3 mRNA molecules calculated at steady state. ⟨m_min,n3⟩ is a parameter representing the minimum number of CLN3 mRNA molecules present in the cell. The full list of stochastic equations is given in S3 Text.

Barik et al. [27] used Gillespie's SSA to simulate all chemical reactions in their stochastic model. They assumed that hyper-phosphorylation of Whi5 corresponds to its nuclear export (as an indication of the Start transition). Their model accurately predicted the timing of the transition and its dependence on cell size at birth, when compared with experimental observations of Di Talia et al. [38].

In the next section, we compare deterministic and stochastic simulations of our SCM to deterministic and stochastic simulations of the MultiP model in S1 Text, and we show that the SCM, a model considerably less complex than MultiP, accounts for most of the dynamic properties of the system.

6. Comparisons of the Start model simulations

To validate our approach, we compare qualitative and quantitative aspects of the SCM to the MultiP model. Fig 2 shows the one-parameter bifurcation curves for both models, in which we plot the steady-state number of ClbS molecules as a function of cell volume (V) as the bifurcation parameter. Both models exhibit bistability within a range of V from ~6 fL to ~30 fL. Hence, both models agree on cell size at the Start transition (~30 fL) for wild-type cells. Newborn cells in G₁ phase of the cell cycle are captured by the stable steady state with ClbS level very low. They must grow to a size of ~30 fL, before the lower stable steady state disappears at a saddle-node bifurcation point. For V > 30 fL, the number of ClbS molecules rises rapidly to the high steady-state level, corresponding to the Start transition. In S2 Fig, we use two-parameter bifurcation diagrams to show how the SCM and the MultiP model behave when the synthesis rates of Cln3 and ClbS are varied at various (fixed) cell size.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. One-parameter bifurcation diagram.

The steady-state number of ClbS molecules is plotted as a function of (fixed) cell volume. Solid line: stable steady states; dashed line: unstable steady states; blue lines: multisite phosphorylation (MultiP) model; red lines: standard component model (SCM). Both models exhibit a region of bistability between V ≈ 6 fL and V ≈ 30 fL. The right bifurcation point (at V ≈ 30 fL) corresponds to the threshold size for the Start transition.

https://doi.org/10.1371/journal.pone.0153738.g002

Next we compare time-series dynamics of the two models in their deterministic formulations (i.e., simulations of the nonlinear ODEs). In Fig 3, we show simulations of a cell that starts with V = 10 fL at t = 0. The two models show very similar time courses for the control proteins. As in the multisite phosphorylation paper [27], we mark the Start transition as the time when the concentration of active SBF increases above 15 nM and the G₁/S transition as the time when ClbS concentration increases above 37.5 nM. Concentration in nM is calculated from molecule number N and volume V in fL by the formula where N_A is Avogadro’s number, 6.02 × 10²³.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Deterministic trajectories simulated by the MultiP model and the SCM of the Start transition.

Both models are simulated for 300 min, starting with V = 10 fL at t = 0. Initial conditions for the SCM are specified in Table 2 and for the MultiP model in S2 Table. In these panels we plot concentration (in nM), which is calculated from molecule number and volume (in fL) by the formula “concentration” = “number”/(0.6 × “volume in fL”). In the MultiP model, Whi5_A = Whi5 + Whi5P₁ + Whi5P₂ + Cmp + CmpP₁ + CmpP₂. In both models, Whi5_i = Whi5_T –Whi5_A. The changes in Whi5_A and Whi5_i over the first 100 min are quite different in the two models because their descriptions of the Whi5 activation process are quite different. Nonetheless, they predict similar timing for the cell cycle transitions. We presume that the Start and G₁/S transitions occur when [SBF] = 15 nM and [ClbS] = 37.5 nM, respectively. (These values are 50% of the maximum concentrations from the original MultiP model [27], not 50% of the final concentrations shown in this figure). The MultiP model (left panels) executes the Start and G₁/S transitions at t = 142 min and t = 152 min, respectively. The SCM (right panels) executes the Start and G₁/S transitions at t = 145 min and t = 153 min, respectively.

https://doi.org/10.1371/journal.pone.0153738.g003

Since it is well known that the smaller is cell size at birth the longer is the time to initiate bud formation and to enter S phase, we measure the durations from birth (t = 0) to the Start transition (T₁) and to the G₁/S transition (T_G1) for various values of initial volume, V₀. In Fig 4, we plot T₁ and T_G1 as functions of initial cell size. (The left-most bars of the figure correspond to the time-series simulations in Fig 3). The results show that the SCM is quantitatively comparable to the more complex MultiP model. The figure also shows that T₁ and T_G1 decrease as birth size increases, consistent with experimental observations [38]. The gap between Start and the G₁/S transition, T₂ = T_G1 –T₁, is small in both models (S3 Fig) and nearly independent of birth size [38].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Deterministic simulations of the relation between initial cell size (V₀) and cell age at the Start transition (T₁) (upper panel) and cell age at the G₁/S transition (T_G1) (lower panel) for the Start models.

Red bars: MultiP model; green bars: SCM. The left-most bars of the figure correspond to the time-series simulations in Fig 3.

https://doi.org/10.1371/journal.pone.0153738.g004

Next we study stochastic properties of the models. The MultiP model is simulated by SSA [13] while the SCM is simulated by the CLE approximation described above. S4 Fig compares the steady-state distributions of Cln3, SBF and ClbS at fixed volume for the MultiP model and the SCM. Fig 5 shows envelopes of sample trajectories from stochastic simulations of both models, starting with V = 10 fL at t = 0. Distributions of T₁ and T_G1 can be computed for the stochastic models as functions of V₀. For each value of V₀, we do 100 independent simulations of each model and compute the average values and standard deviations of T₁ and T_G1. The average values show similar patterns to the deterministic simulations in Fig 4 (see S5 Fig). The coefficient of variation (CV) is shown in Fig 6 and the distributions of T₁ and T_G1 are shown in S6 Fig. The results show that the noise intensities of both models agree well at small V₀ but not at large V₀ (red and green bars in Fig 6). Our SCM without mRNA fluctuations (i.e., discarding the mRNA-inherited noise term from the equations) shows less intensity of noise (blue bars in Fig 6), indicating that mRNA fluctuations are important. By expanding the SCM to include mRNA synthesis and degradation explicitly and simulating mRNA fluctuations using CLEs (see equations in S4 Text), we obtained much better agreement with the fully stochastic MultiP model (Fig 7), but at a cost of considerably more complexity.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Stochastic trajectories generated by the MultiP model and the SCM of the Start transition.

Means (lines) and standard deviations (shadows) calculated from 100 independent trajectories are shown for each model, starting with V = 10 fL at t = 0. Initial conditions for the SCM are specified in Table 2 and for the MultiP model in S2 Table. As in Fig 3 we plot “concentration in nM” = “number of molecules”/(0.6 × “volume in fL”). The MultiP model is simulated by Gillespie’s SSA and the SCM is simulated by the chemical Langevin approach described in the text.

https://doi.org/10.1371/journal.pone.0153738.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Stochastic simulations of T₁ and T_G1 for the Start models.

As in Fig 4, we compute the times from birth (t = 0, V = V₀) to the Start transition (T₁, when [SBF] = 15 nM for the first time; upper panel) and to the G₁/S transition (T_G1, when [ClbS] = 37.5 nM for the first time; lower panel). For each model we compute 100 stochastic trajectories and calculate the mean and standard deviation of the time to the event. The mean times agree well with the deterministic simulations in Fig 4. Here we plot the coefficient of variation (CV = standard deviation/mean) of the times, in order to judge how well stochastic CLE simulations of the SCM (green bars) compare with SSA simulations of the MultiP model (red bars). Clearly the SCM under-estimates the variability of the transitions at large birth size. Removing mRNA noise from the SCM (blue bars) makes matters worse, as expected.

https://doi.org/10.1371/journal.pone.0153738.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Stochastic simulations of T₁ and T_G1 for the Start SCM with explicit account of fluctuating mRNA species.

As in Fig 6, except now we have added synthesis and degradation of mRNA species explicitly to the SCM. The remaining discrepancies are attributable in part to differences between the models and in part to simulating mRNA noise by CLE rather than SSA.

https://doi.org/10.1371/journal.pone.0153738.g007

With this simple model of the Start transition in budding yeast, we have shown that the SCM approach can capture essential dynamical features of a complex control network in quantitative detail. Encouraged by these results, we built an SCM of the full network of molecular controls of the budding yeast cell cycle, as described in the next section.

7. An SCM for the budding yeast cell cycle control system

Our cell cycle model (Fig 8) is based on the interaction network proposed in [41] with some modifications in the Start transition (by adding the inhibition of SBF by Whi5) and mitotic exit (by adding the stimulatory effects of Polo kinase and separase Esp1 on Cdc14 release), according to the models proposed in [36] and [42], respectively. In the model of Chen et al. [41] (Chen-2004, hereafter), each component is represented in terms of a concentration that has been scaled to a dimensionless number, called its normalized concentration, [] _n. Cell volume, V_n, is also dimensionless. In order to compare SCM simulations directly to Chen-2004, we switch from numbers of molecules to normalized concentrations, using the “characteristic concentrations” listed in Table 3. To be more specific, [S]_n = [S]/c_S = N_S/(0.6 V∙c_S), where [S] = concentration of species S in nM, c_S = characteristic concentration of S in nM, N_S = number of S molecules in volume V (in fL), and V = V_n c_vol.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Variables, initial values and characteristic concentrations for the standard component model of the full cell cycle control system.

https://doi.org/10.1371/journal.pone.0153738.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Wiring diagram of the full cell cycle control network in budding yeast.

The network consists of three major modules: Start, S/G₂/M, and Exit. Red and blue icons: active forms of components; orange icons: inactive forms. Solid lines: chemical reactions (synthesis and degradation, phosphorylation and dephosphorylation, association and dissociation); dashed lines: activatory or inhibitory influences of components on chemical reactions. T-shaped reaction arrows with black circles on the reactants side of the arrow indicate reversible association of two proteins to form a complex. T-shaped arrows without black circles represent irreversible reactions. Not all reactions are shown on this figure; see the equations in Table 4 for complete details.

https://doi.org/10.1371/journal.pone.0153738.g008

In budding yeast, as in all eukaryotes, cyclin/Cdk complexes are the main regulators of cycle progression. The kinase subunit, Cdc28, being far more abundant than all cyclin proteins, is always available to bind to any cyclins present in the cell [31]. Therefore, we use concentrations of cyclins to indicate the activities of cyclin/Cdc28 complexes (neglecting Cdc28 itself in the model). As in previous models, components that have similar functions are treated as single variables; for instance, “SBF” stands for SBF and MBF, “Cln2” for Cln1 and Cln2, “Clb2” for Clb1 and Clb2, “Clb5” for Clb5 and Clb6, and “CKI” for Sic1 and Cdc6.

For the cell cycle control system, class-1 variables include proteins whose concentrations change slowly over time due to the activities of transcription factors and proteolytic enzymes. This class includes, for example, the total concentrations of cyclin proteins (Cln2, Cln3, Clb2, and Clb5) and the total concentration of the stoichiometric inhibitor CKI. The synthesis of Cln2 and Clb5 is regulated by the transcription factor SBF, the synthesis of Clb2 is regulated by the transcription factor Mcm1, and the synthesis of CKI is controlled by the transcription factor Swi5. The degradation of Clb2 and Clb5 is regulated by the proteolytic factors Cdc20 and Cdh1 in association with an E3 ubiquitin ligase called the Anaphase Promoting Complex (APC).

The properties of many cell cycle-control proteins can be regulated by phosphorylation and dephosphorylation. For instance, phosphorylated Whi5 is a less potent inhibitor of SBF, and phosphorylated CKI is more susceptible to proteolysis. We use class-2 variables to represent the fractions of these proteins that are in “active” forms. The regulation of protein activity by phosphorylation and dephosphorylation is believed to be essential for the bistable switches that govern cell cycle transitions [6, 27, 36, 41, 43–45]. Class-2 variables capture the nonlinear behavior of these reactions in terms of the “soft-Heaviside” functions built into the SCM formulation.

Class-3 variables represent the free forms of proteins that form tight complexes with stoichiometric binding partners. For example, both Clb5 and Clb2 are inhibited by binding to CKI, as is Cdc14 inhibited by binding to Net1. In these cases, the “max” function employed by SCMs is a reasonable approximation because stoichiometric inhibitors usually bind strongly and rapidly to their partners. This strong-binding assumption is part of the original ODE model as well [41].

In Fig 8 we have divided the regulatory network into three modules (Start, S/G₂/M, and Exit), which provides a useful framework for constructing the SCM. In Table 3 we classify the regulatory proteins into the three classes of SCM variables, and in Table 4 we translate Fig 8 into the differential and algebraic equations that define the dynamics of the network according to the SCM approach.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Equations for the standard component model of the full cell cycle system.

https://doi.org/10.1371/journal.pone.0153738.t004

Formulation of the Exit module.

For a budding yeast cell to exit from mitosis and return to G₁ it must transiently activate a phosphatase, Cdc14, which dephosphorylates many proteins that had been phosphorylated by cyclin-dependent kinases in S/G₂/M. In particular, Cdc14 activates Cdh1 and CKI. These two factors stabilize G₁ phase of the cell cycle by repressing the activity of Clb-dependent kinases. The activation of Cdc14, as described by the Exit module in Fig 8C, is accomplished by two pathways: FEAR (Cdc fourteen early anaphase release) and MEN (mitotic exit network) (for a review, see [70]). In the FEAR pathway, active Cdc20 cleaves Pds1 (securin) [66] and releases Esp1 (separase) [67]. Free Esp1 is responsible for chromatid separation [67] and inactivation of PPX (PP2A^Cdc55), a phosphatase that dephosphorylates Net1 [71]. Subsequent phosphorylation of Net1 causes Cdc14 to be released from Net1:Cdc14 complexes in the nucleolus, and free Cdc14 is the phosphatase that drives exit from mitosis [71–75]. A schematic diagram of the FEAR pathway is: Cdc20 → Esp1—| PPX → Net1—| Cdc14.

We use a class-1 variable to describe the total concentration of Pds1 (Eq 55) and a class-3 variable to represent free Esp1, i.e., Esp1 that is not bound to Pds1 (Eq 56). The active forms of PPX and Net1 are described by class-2 variables (Eqs 57 and 59, respectively). The active (unphosphorylated) form of Net1 is a stoichiometric inhibitor of Cdc14, so we use a class-3 variable to represent free Cdc14 (Eq 61). In Eq 61 we introduce a stoichiometric factor, ρ_14,net1, in order to simulate two mutants, net1-ts and TAB6-1, that show reduced association between Net1 and Cdc14 [74–76]. For wild-type cells we set ρ_14,net1 = 1 and for net1-ts and TAB6-1 cells we set ρ_14,net1 < 1.

The FEAR pathway is not sufficient to return cells to G₁ phase because Cdc14 activates Cdh1 which degrades Clb2, allowing Net1 to regain its ability to sequester Cdc14 in the nucleolus. To get robust phosphorylation of Net1 and full release of Cdc14 from the nucleolus, the FEAR pathway must be supplemented by the MEN pathway, whose role is to activate Cdc15 and Tem1 (see Fig 8C). Active Cdc15 then supports continued phosphorylation of Net1 even as Clb2-dependent kinase activity is dropping. Premature activation of Cdc15 is prevented by Clb2-dependent phosphorylation (inactivation) of Cdc15 in metaphase. Hence, Cdc15 activity is low as cells enter anaphase and rises abruptly as Clb2 is degraded by Cdc20_A:APC_P and Cdc14 is activated by the FEAR pathway (Eqs 67 and 68) [77].

Meanwhile, Polo kinase, which was activated in prometaphase when Clb2-kinase activity was high (Eqs 62–64) [78], is able to activate Tem1 after PPX is inactivated by the FEAR pathway (Eqs 65 and 66) [71, 79]. Together, Tem1 and Cdc15 form an active complex (Tem1_A:Cdc15_A) (Eq 69) that phosphorylates Dbf2/Mob1 (not modeled explicitly), which then phosphorylates and inactivates Net1, resulting in full Cdc14 release. In this way, the transient release of Cdc14 by the FEAR pathway initiates a positive feedback loop between the activation of Tem1_A:Cdc15_A (MEN) and further release of Cdc14 [42].

Full release of Cdc14 re-sets the cell to G₁ by two complementary actions. First of all, Cdc14 dephosphorylates and activates Cdh1, and active Cdh1 (in association with APC) fully degrades Clb2 [68, 80, 81] and Polo kinase [82]. In addition, Cdc14 stabilizes CKI and activates its transcription factor, Swi5 [62, 83]. Abundant CKI and active Cdh1 are characteristic molecular signatures of G₁ phase in budding yeast.

We assume that cell division occurs when the concentration of Clb2 drops below a threshold value, K_EZ. Because budding yeast cells divide asymmetrically, we set the sizes of mother and daughter cells after division equal to 58% and 42% of the size of the dividing cell in glucose medium, which are close to the proportions observed by Di Talia et al. [36]. In galactose medium, for which cells grow more slowly and divide more asymmetrically [84], we take these proportions to be 61% and 39%.

BUD, ORI, and SPN, the dimensionless variables used to mark bud emergence, onset of DNA synthesis, and spindle assembly, are reset to 0 as the cell exits mitosis. We reset [BUD]_n and [SPN]_n to 0 at cell division, when [Clb2]_n = K_EZ. We reset [ORI]_n to 0 when [Clb2]_n+[Clb5]_n drops below K_EZ2, because origins of replication are relicensed only if all Clb-dependent kinase activity is extinguished in G₁ phase.

Full lists of initial conditions and kinetic constants used to simulate wild-type cells are given in Tables 3 and 5, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Parameter values for wild-type cells.

https://doi.org/10.1371/journal.pone.0153738.t005

The full SCM (Table 4) tracks 33 protein species (not counting V, BUD, ORI and SPN) and requires ~100 parameter settings (Table 5). By comparison, Chen-2004 tracks the same number of protein species and requires ~120 parameter settings.

8. Deterministic simulations of wild-type and mutant cells

To simulate cell cycle progression in wild-type and mutant yeast cells, we solve the ODEs in Table 4 with the parameter values in Table 5. To simulate a broad selection of mutant yeast strains (S3 Table), we made appropriate changes to the “wild-type” parameter values, as outlined in S4 Table.

Our choice of wild-type parameter values (Table 5) was guided initially by the rate constant assignments in Chen et al. [41] and then adjusted manually to give a good fit to the observed phenotypes of the 133 mutant strains in S3 Table. We were able to account for the observed phenotypes of 125 of the 133 mutant strains (94%) in the data set we used to constrain the model. We tried automatic exploration of parameter space by a genetic algorithm (“Differential Evolution”) but could not find a set of parameter values that improves on the 94% success rate achieved manually [85].

The simulation of wild-type cells (Fig 9) shows oscillating patterns of cell cycle control variables that are in agreement with observations and expectations. The Start transition (Fig 9A) is initiated by the inactivation of Whi5 as Cln3 accumulates, followed by the positive feedback-driven expression of Cln2. Fig 9B shows an alternating pattern of G₁ phase (during which G₁ stabilizers CKI and Cdh1 are abundant) and S/G₂/M phase (during which Clb5 and Clb2 are abundant). During mitotic exit, the release of Cdc14 resets the system back to G₁ phase. Quantitative comparisons to experimental data will be discussed in the next section on stochastic simulations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Simulations of wild-type cells.

(A) Start: As the cell grows (increasing V_n), Cln3 accumulates and phosphorylates Whi5. At a critical cell size, SBF is abruptly released from the inactive SBF:Whi5 complex and initiates a positive feedback loop between the accumulation of Cln2 and the phosphorylation of Whi5. (B) G₁/S/G₂/M: SBF also promotes the synthesis of Clb5. Once Clb5 titrates out CKI, then Clb5 and Cln2 together inactivate Cdh1, resulting in the accumulation of Clb2. Exit: Clb2 triggers many mitotic events, eventually leading to the release of Cdc14 during mitotic exit. When Clb2 drops below a normalized concentration of 0.4, the cell divides asymmetrically between daughter and mother cells. The daughter cell receives 42% of the cell size at division, and the mother cell (not shown here) receives the remaining 58%. (C and D) The stochastic model shows the typical fluctuations of protein concentrations around the average dynamics predicted by the corresponding deterministic model. For easier comparison to the deterministic simulation (A and B), we converted the numbers of molecules reported by the stochastic simulation to normalized concentrations. Start and division events are indicated by up-pointing and down-pointing black triangles, respectively.

https://doi.org/10.1371/journal.pone.0153738.g009

To verify our model, we computed the expected phenotypes of 133 mutant strains (S3 Table) and compared our results to the original model of Chen et al. [41] and to available experimental observations. (See our web site at http://mpf.biol.vt.edu for the observed phenotypes of mutant strains and references to the original literature.) Simulations that exhibit periodic cell division correspond to viable mutants, and simulations that arrest at some point in the cell cycle correspond to inviable mutants, with phenotypes assigned according to the rules in S5 Table. In our simulations, 125 mutants show phenotypic characteristics in agreement with the original model and experiments. (Notice that our SCM achieves the same accuracy as Chen-2004.) The eight mutant simulations that are not in agreement with observed phenotypes are listed in S6 Table and discussed in detail in S5 Text.

9. A stochastic SCM of the budding yeast cell cycle

In this section, we construct a stochastic version of our SCM, in order to address three questions. (1) Is our model of the cell cycle control system robust in the face of inevitable molecular fluctuations within single cells? (2) Is a stochastic version of our model consistent with quantitative measurements of cell cycle variability among wild-type cells? (3) Do stochastic models behave differently than deterministic models under some circumstances? With regard to the third question, we have in mind situations where mutations may make the regulatory mechanism less robust and more sensitive to molecular noise. For example, in the double mutant strain CLB2-dbΔ clb5Δ (CLB2 destruction box deletion and CLB5 gene deletion) some cells are able to complete the cell cycle and divide whereas other cells become arrested in telophase and eventually die, as originally observed by Cross [86] and confirmed by Ball et al. [87]. This sort of partial viability is clearly incompatible with a deterministic model of the cell division cycle but may be accommodated by a stochastic model.

With regard to question (2), single-cell experimental techniques (e.g., flow cytometry and live-cell imaging of fluorescent proteins) have revealed much information about cell-to-cell variation among genetically identical cells. In budding yeast, fluorescence microscopy has been used to study proteins that regulate cell cycle progression and to determine the onset of cell cycle events such as Start, bud emergence, and cell division in individual cells [38, 88, 89]. In addition, fluorescence in situ hybridization (FISH) has been used recently to quantify mRNA levels in individual yeast cells [90], determining that a yeast cell carries roughly 5–15 mRNA molecules of each gene measured, including a sample of cell cycle control genes [91]. Since the number of mRNA molecules directly determines the production rate of its protein product, noise due to low, fluctuating numbers of mRNA molecules will be significantly amplified as fluctuations in protein abundance [39]. Consistent with this expectation, stochastic simulations by Kar et al. [92] suggested that mRNA noise is a major source of fluctuations in the budding yeast cell cycle control system. Therefore, to understand stochasticity in the system, inclusion of mRNA noise seems to be necessary. The stochastic version of SCM described above was designed specifically to deal with molecular noise derived from both protein and mRNA fluctuations.

10. Stochastic simulations of wild-type and mutant cells

Stochastic simulations of wild-type cells (e.g., Fig 9C and 9D) show high variability in progression through the cell cycle. Following the lead of Cross and his group [38], we divide the cell cycle into three periods: T₁, T₂, and T_b. T₁ is the duration from cell birth to Start, which we identify with SBF reaching 50% of its maximum value. T₂ is the period from Start to bud emergence, which we identify with [BUD]_n = 1. Finally, T_b is the duration of the “budded phase,” from bud emergence to the next cell division.

From stochastic simulations of wild-type cells, we computed means and variances of T₁, T₂, T_b, cell cycle duration, and cell size at birth, which we compare to experimental observations [38] in Fig 10. (The full distributions of T₁, cell cycle duration and V_birth, for both mother and daughter cells, are shown in S7 Fig.) Overall, our results show a reasonable agreement with the experimental data, although there are some quantitative discrepancies. Our simulated T₁ is longer and our simulated T₂ is shorter than observed in both mother and daughter cells, which may be attributed to a difference between our theoretical criterion for Start (SBF reaching 50% of its maximum activity) and the experimental criterion for Start (Whi5 disappearance from the nucleus). T₁ + T₂ = T_G1 = duration of the “unbudded phase” is longer than observed in mother cells, but quite accurate for daughter cells. The predicted and observed CVs of T_G1 are all ~50%. In our simulations the average delay from bud emergence (when [BUD]_n = 1) to the initiation of DNA synthesis (when [ORI]_n = 1) is ~12 min, which is inconsistent with the observation that the two events occur nearly simultaneously in wild-type cells.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Statistical properties of cell cycle progression.

Experimental observations for wild-type cells [38] (blue bars) are compared to simulations from four different cases of our model: the full stochastic SCM (red bars), the stochastic SCM without the mRNA-inherited noise term (green bars), the deterministic SCM with “extrinsic” noise only (magenta bars), and the completely deterministic SCM (black bars). T_c is cell cycle duration (min) = T_G1 + T_b = T₁ + T₂ + T_b. T₁ is the duration from cell birth to Start (min), which we identify with SBF reaching 50% of its maximum value. T₂ is the period from Start to bud emergence (min), which we identify with [BUD]_n = 1. T_G1 = T₁ + T₂ = duration of the “unbudded phase” (min). T_b is the duration of the “budded phase,” from bud emergence to the next cell division (min). V_birth is cell size at birth (fL). Asterisks indicate unreported data.

https://doi.org/10.1371/journal.pone.0153738.g010

The duration of G₁ is negatively correlated with birth size in experiment and in simulation.

Next, we study the correlation between size at birth and duration of G₁ phase (T_G1). If cell cycle progression during G₁ is perfectly controlled by cell size, then a plot of μ T_G1 against ln(V_birth) should give a straight line with slope = –1 [38], where μ is the specific growth rate and V_birth is the birth size of cells. In the experiment of Di Talia et al. [38], small daughter cells show the expected negative correlation (slope = –0.7), indicating a strong size-control mechanism operating in small budding yeast cells. Larger daughter cells show less negative correlation (slope = –0.3). Mother cells seem to operate with little or no size control (slope = –0.1). In good accord with experimental observations, our simulations (Fig 11) show that the small and large daughter cells exhibit correlations with slopes –0.67 and –0.30, respectively. The break between strong size control and weak size control occurs at a birth size of , where = average size of mother cell at birth. The mother cells show less correlation (slope = –0.24).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. The joint distributions of cell size at birth and G₁ phase duration predicted by our stochastic model of the full cell cycle control system.

(Upper panel) Daughter cells. (Lower panel) Mother cells. T_G1 is the duration from cell birth to bud emergence, and μ is the specific growth rate of the cell culture. In the plots, cell size is normalized by the average size of mother cells at birth and is plotted on a log scale. Small dots in the background represent data collected from simulations of individual cells. Large dots represent average μ·T_G1 of the small dots binned in 2 fL intervals. For daughter cells, the binned data can be divided into two groups (small and large cells; break point at ln(V_birth) = −0.37) and fitted by two straight lines with slopes of –0.67 and –0.30 (respectively). Mother cells can be fitted by one straight line with slope –0.24. The experimental slope values, as reported by Di Talia et al. [38], are shown in parentheses. The binned data for mother cells is slightly better fit by two straight lines with a break point at ln(V_birth) = 0.

https://doi.org/10.1371/journal.pone.0153738.g011

It is beyond the scope of this paper to provide a statistically rigorous analysis of the negative correlation between G₁ duration and size at birth predicted by the SCM. In a separate study, we have fitted our simulations to the full experimental distributions (kindly provided by Stefano Di Talia). We subdivided the empirical joint distribution into 17 rectangles (“bins”) so that no bin contained more than about 20 cells or less than 3, and we subdivided the simulated joint distribution in an identical fashion. Then we computed the Hellinger distance between the two distributions (empirical and simulated) and minimized the distance by a quasi-Newton algorithm for stochastic optimization [95]. By using the full distribution (rather than the two-slope representation) we were able to improve the fit of the model to the data, but there remained a noticeable discrepancy in that the model predicts longer tails of G₁ durations than are observed experimentally for both mother and daughter cells. Hence, although the model is qualitatively in accord with the observed negative correlation between T_G1 and birth size, it exhibits long tails in the T_G1 distributions that appear to be statistically significant deviations from observed G₁ durations. These results will be reported in full in a later publication.

In S8 Fig we plot the correlations of cell size at birth with both T_G1 and T₁ and confirm the experimental observation [38] that T₂ = T_G1 –T₁ is relatively constant. This observation indicates that the size control mechanism in budding yeast operates at the Start transition, i.e., at the phosphorylation and inactivation of Whi5. Whi5 inactivation frees up SBF to drive Cln2 and Clb5 production. The initial phosphorylation of Whi5 is the job of Cln3 and Bck2, and since their rates of accumulation are proportional to cell size, V, small cells take a longer time to inactivate Whi5, and large cells take a shorter time.

Simulated numbers of protein molecules for an asynchronous population are in agreement with observations.

Next, we study the abundances of protein species in an asynchronous population, as might be measured experimentally. To this end, for each simulated cell we select a random time point between birth and division and record the number of each protein species present at that time. To obtain a sample of time points representative of an asynchronous population of cells, we must take into account that the number of newborn cells is twice the number of dividing cells. Hence, we select a random number θ, 0 ≤ θ ≤ 1, from an exponential probability density function (71) and then calculate t_birth + (t_division−t_birth)·θ as a random time point in each complete cell division cycle. By this prescription, the probability of choosing a cell at birth is twice the probability of choosing a cell at division. (Eq 71 is not exactly correct for asymmetrically dividing cells, like budding yeast, but the errors introduced by this approximation are minor.)

This set of simulated protein abundances should be comparable to numerical values collected from an asynchronous population of yeast cells. In Table 6 we show that our computed average protein abundances are in good agreement with experimental values [31, 96], except for Clb5 and Clb2, whose predicted abundances are two-fold larger and smaller, respectively, than observed.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Average protein abundances.

https://doi.org/10.1371/journal.pone.0153738.t006

Stochastic simulations of the mutant strain CLB2-dbΔ clb5Δ explain its partial viability when growing on raffinose medium.

Next we study viability of the mutant strain CLB2-dbΔ clb5Δ. In this strain, Clb2’s destruction box has been deleted from the CLB2 gene, so Clb2 protein cannot be degraded by Cdc20 and only partially so by Cdh1 [86]. This mutant strain is inviable on glucose medium (mass-doubling time = mdt = 100 min) but can be partially rescued on raffinose medium, which supports slower growth (mdt = 150 min) [86, 87]. Obviously, “partial viability” (some cells complete the cell cycle and some do not) cannot be explained by a deterministic model in which all cells behave the same. A stochastic version of Chen's 2004 model [87] predicts that some mutant cells, growing in raffinose, divide properly while others fail to exit from mitosis. However, the model in [87] applied Gillespie's SSA directly to Chen's 2004 model without unpacking the complex rate laws into elementary steps, and it did not take into account transcription-translation coupling. Hence, the results of the model in [87], though suggestive, are not very reliable. It is challenging to test whether our SCM approach provides a better account of the properties of this mutant strain.

Since the mutant cells cannot survive on glucose medium and are only partially viable on raffinose, the mutant strain is kept viable in the laboratory by introducing the GAL-SIC1 gene into the CLB2-dbΔ clb5Δ cells and growing the strain in galactose, where Sic1 is overexpressed. In Fig 12, top panel, we show predicted growth curves for CLB2-dbΔ clb5Δ GAL-SIC1 cells grown in media containing galactose (GAL-SIC1 is “on”) or raffinose (GAL-SIC1 is “off”). These cells proliferate in galactose with a population number-doubling time (ndt) ≈ 150 min, which is equal to the mass-doubling time (mdt) of individual cells. On the other hand, for the triple mutant cells growing in raffinose, ndt > 150 min, indicating that some cells are inviable. Cell viability is better illustrated in Fig 12, bottom panel, where we plot cumulative distribution functions for triple-mutant cells in galactose and in raffinose. Greater than 90% of the cells growing in galactose (mdt = 150 min and Sic1 overexpressed) complete the cell cycle within 250 min of birth, indicating that they are mostly viable. On the other hand, for the same cells growing in raffinose (mdt = 150 min but Sic1 not overexpressed), ~ 25% of the cells are undivided even 300 min after birth, indicating that ~25% of the cells growing in raffinose never complete the cell division cycle. This result agrees well with the experimental observation [87] that ~15–40% of these mutant cells never divide in raffinose growth medium. Furthermore, our simulations of mutant cells grown in glucose medium (mdt = 100 min and Sic1 not overexpressed) show less than 10% viability, which is consistent with the fact that these cells are inviable in glucose medium [86]. In Table 7, we compare the statistical properties of this mutant strain grown in raffinose from simulations and from experimental measurements [87].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 7. Statistical properties of simulated CLB2-dbΔ clb5Δ cells in raffinose.

Mean cycle time (in minutes) with standard deviation in parentheses.

https://doi.org/10.1371/journal.pone.0153738.t007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Growth curves and cycle time distributions for CLB2-dbΔ clb5Δ GAL-SIC1 cells.

(Upper panel) Logarithm of the total number of cells is plotted against time. The increase in cell number for CLB2-dbΔ clb5Δ GAL-SIC1 cells growing in galactose (red crosses; Sic1 overexpressed) indicates exponential growth (black line) with the number doubling time = 150 min. The number doubling time of the same cells growing in raffinose (blue crosses; normal level of Sic1) is greater than 150 min, indicating that some cells do not complete the cell cycle when GAL-SIC1 is not expressed. (In our simulations, the mutant cells are assumed to have mass doubling time = 150 min in both galactose and raffinose media.) (Lower panel) Cumulative distributions of cycle times for CLB2-dbΔ clb5Δ GAL-SIC1 cells growing in galactose medium (red line) or in raffinose medium (blue line). The ordinate, P(t), is the probability that a simulated cell has a cycle time greater than t.

https://doi.org/10.1371/journal.pone.0153738.g012

As in the original model of Chen et al. [41], we assume that the synthesis of Clb2 is dependent on cell size. By this assumption, cell size at the onset of mitotic exit is critical for the cell's fate. Cells with larger size at the onset of mitotic exit have more Clb2 and are less likely to exit from mitosis. (Cells exit only if Clb2 is degraded below a threshold concentration.) In wild-type cells, the presence of both Cdc20 and Cdh1 ensures that Clb2 is fully degraded during mitotic exit. However, Clb2 protein encoded by the CLB2-dbΔ gene lacks the sequence targeted by Cdc20 and can be only partially degraded by Cdh1. To exit from mitosis, CLB2-dbΔ cells depend critically on the activation of Cdh1 and Sic1, but both Cdh1 and Sic1 are also inhibited by Clb2-dependent kinase. When grown on a rich medium such as glucose, the mutant strain has a high growth rate and Clb2 accumulates to such a high level that it keeps Cdh1 and Sic1 inactive, resulting in telophase arrest. On a slower growth medium such as raffinose or galactose, a cell may exit from mitosis if Cdh1 and Sic1 can suppress Clb2. This explains why cells show high viability when Sic1 is overexpressed in galactose. In raffinose, however, without the help of Sic1 overexpression, noise plays a major role in determining whether or not a cell can exit from mitosis and divide. Cells with higher activities of Sic1 and/or Cdh1 or lower levels of Clb2 may be able to divide, but cells at the other extreme will be arrested in anaphase.

In S9 Fig we show that CLB2-dbΔ GAL-SIC1 (CLB5 in place) mutant cells growing in raffinose (Sic1 not produced) have comparable viability to CLB2-dbΔ clb5Δ GAL-SIC1 cells growing in raffinose, because Clb5 protein is mostly degraded by APC/Cdc20 by the time the CLB2-dbΔ CLB5 cells are exiting (or not) from mitosis. Although CLB2-dbΔ CLB5 GAL-SIC1 cells are inviable in glucose medium (rapidly growing cells that are not producing Sic1), there is anecdotal evidence that these cells are partially viable when cultured on poor growth medium, like raffinose [86]. So, our stochastic model is consistent with a conclusion that CLB2-dbΔ CLB5 cells are on the cusp of viability/inviability when growing on poor carbon sources.