Abstract

Using a simplified budding yeast cell cycle model perturbed by intrinsic noise, we systematically explore these issues from an energy landscape point of view by constructing an energy landscape for the considered system based on large deviation theory. Analysis shows that the cell cycle trajectory is sharply confined by the ambient energy barrier, and the landscape along this trajectory exhibits a generally flat shape. We explain the evolution of the system on this flat path by incorporating its non-gradient nature. Furthermore, we illustrate how this global landscape changes in response to external signals, observing a nice transformation of the landscapes as the excitable system approaches a limit cycle system when nutrients are sufficient, as well as the formation of additional energy wells when the DNA replication checkpoint is activated. By taking into account the finite volume effect, we find additional pits along the flat cycle path in the landscape associated with the checkpoint mechanism of the cell cycle. The difference between the landscapes induced by intrinsic and extrinsic noise is also discussed. In our opinion, this meticulous structure of the energy landscape for our simplified model is of general interest to other cell cycle dynamics, and the proposed methods can be applied to study similar biological systems.

Author Summary

By constructing a global energy landscape for a simplified yeast cell-cycle regulatory network, we provide a systematic study of this issue. Our results demonstrate that the cell cycle trajectory is sharply confined as a canal bounded by ambient energy barriers, with the landscape adaptively reshaping itself in response to external signals, such as the nutrients improving and the activation of DNA replication checkpoint in our work. After performing quantitative analysis based on the landscape, We found that along the cell cycle trajectory, the typical

Interestingly, this is also basically in accordance With the force strength of the dynamics. Additionally, in places Where the driVing force strength is comparable to the noise level, some additional pits form that are associated With the checkpoint mechanisms. Overall, our energy landscape study shows that the yeast cell cycle is a robust, adaptive and multistage dynamical process.

Introduction

However, it is still difficult to quantify the robustness and adaptivity of cellular networks, even for a small cellular network perturbed by intrinsic random fluctuations, due to the massive cross regulations and nonlinear nature of such biological systems. As the size of the network grows, determining how to characterize the global stochastic dynamics of the system becomes a tough problem. “Waddington’s epigenetic landscape,” which utilizes potential energy to pictorially illustrate the dynamics and evolution of cellular networks, has been widely and repeatedly used for several decades [5, 7]. Some beautiful efforts and frameworks aiming to quantify this landscape have been made [8—13], but investigation into typical biological models still remains to be done. Furthermore, the energy landscape usually reshapes itself due to a variety of changes such as environmental signals [14], cell-cell interactions [15] and the growth rate dependence of protein concentrations [16]. Determining how to explicitly quantify this transformation for specific systems is also a major task.

Recent work has revealed the dynamic regulatory mechanisms of the cell cycle, and the cell cycle process is now considered a series of irreversible transitions from one state to another [19—21]. The cell-cycle regulatory network must also be robust and adaptive to external stresses and signal changes. To quantitatively characterize this robustness, and provide a global description of the cell cycle regulatory system, some fundamental questions must be studied. For example, how does the energy landscape reflect the robustness and successive phases of the cell cycle? How does the landscape adaptively change in response to external signals? Is there any information that the energy landscape cannot provide? If so, does any other supplemental description exist?

Our results demonstrate that the energy landscape of the cell cycle is globally attractive, and we show how the cell cycle regulatory network reduces fluctuations from its upstream process and enables long durations in the transition regime. We also describe how the landscape changes in response to external signals when nutrients become sufficient and the DNA replication checkpoint is activated. We also discuss the dynamic information provided by the non-gradient nature that the pure energy landscape cannot explain, and provide other approaches to take into account this non-gradient effect. In addition, we compare the difference between landscapes induced by intrinsic and extrinsic noise and discuss the finite volume effect. Overall, our energy landscape study shows that the budding yeast cell cycle is a robust, adaptive and multistage dynamical process.

Methods

Models

That is, once the yeast cell passes the Start checkpoint, it will proceed through the whole cell cycle process spontaneously. Based on the key regulatory network [17] and our previous study on budding yeast [22] , the cell cycle regulatory network can be simplified and separated into G1 /S, early M and late M modules, as shown in Fig. 1A. We ignore the G2 phase for simplification. Each module has a positive feedback, and different modules are connected with activation and repression interactions. The deterministic equations describing this three-module yeast cell cycle network are where x represents the concentrations of key regulators such as cyclins Cln1, 2, Clb5, 6 and transcriptional factors SBF and MBF in the excited G1 and S phases; y represents the concentrations of key regulators such as cyclins Clb1, 2 and transcriptional factor Mcm1/SFF in the early M phase; and 2 represents the concentrations of key inhibitors such as th1, Cdc20 and Sic1 in the late M/G1 phase.

Thus, the first term in each equation, the second order Hill functions, represent the positive feedback in each module [23, 24]. The second term represents the degradation rate of each regulator, while the third term represents the repression or inhibition interaction between different modules. The parameter a0 in Equation (1a) characterizes the environmental nutrition condition [25, 26], and kal x and kaz y are the trigger signals from x to y and y to 2 respectively. Starting from the excited G1 state, the system will finally evolve to a stable fixed point, the G1 state; if a0 is large enough, the system will enter the cell cycle process continually. With a proper parameter set, the model in Equation (1a) ensures a successive event order from DNA replication in the S phase to mitosis in the M phase, as well as a long duration for both events in the cell cycle process. For this work, we will simply denote the set of equations in Equation (1a) as dx/dt = b(x), where x = (x, y, z).

Here P1 (x = (O, O, zmax)) represents the G1 state—which is the globally stable state of our model where a0 = 0.001—and P2 is a saddle point used to represent the excited G1 state from which the yeast cell passes the Start checkpoint and enters the cell cycle process. The trajectory in Fig. 1C can be separated into three parts. The first part is from the excited G1 (P2) to the S phase (P3), where x = (xmax, O, O); the second part is from P3 to the early M state before the metaphase/anaphase transition (P4), where x = (O, ymax, O); and the third part evolves from P4 to the stable G1 state (P1). Compared with the model used in [9], our model does not rely on a quasi-steady-state assumption related to cell mass. More details about the network and model can be found in $1 Text and $1 Fig. Starting from the deterministic descriptions above, we now address the stochastic setup of the system. The noise can be classified into the intrinsic and extrinsic types [1]. Here we model the intrinsic noise through a Gillespie jump process [27], in which the strength of the noise is determined by the reaction network structure. Denote the state of the system X = (X, Y, Z) where each component represents the number of molecules for the corresponding specie. We then translate each term in Equation (1a) into a chemical reaction. Taking Equation (1a) as an example, we have four associated reactions for the four terms. The state change vector v for each reaction channel has the form v1 2 v4 = [1, O, O] and v2 2 v3 = [—1, O, O], which corresponds to the plus or minus sign in the equation. Once a reaction fires, the state of the system X would be updated to X+v. The reaction propensity function is determined by each term and the volume size (or system size) V, where e E V1 characterizes the magnitude of intrinsic fluctuations [28]. In Equation (1a), the four propensity functions are

Equation (1b) and (1c) can be treated similarly. There are 12 reactions in total. The above setup is suitable for the intrinsic noise. For the extrinsic noise Whose magnitude is independent of the considered system, we simply take the stochastic model as )2: = b(x) —|— \/Ev'v, Where w is the standard temporal Gaussian White noise. More details can be referred to 81 Text.

Algorithm

To study the robustness and adaptivity of our cell cycle model, we construct the Waddington-type energy landscape based on the concept quasi-potential from large deviation theory [29]. For any stable steady state x0 of Equation (1a), the local quasi-potential S(x; x0) With respect to x0 is defined as

(p is any connecting path and L is called the Lagrangian [29—31]. The concrete form of L is determined by the setup of the stochastic process defined through intrinsic or extrinsic fluctuations. In case that the driVing noise is of White noise type, L can be obtained from path integral formulation [32]. S(x; x0) tells us the difficulty of transition from state x0 to x under the noise perturbation. The local quasi-potentials starting from different stable steady states can be suitably integrated together to form a global quasi-potential S(x), Which is exactly our proposal to rationalize the Waddington landscape for any non-gradient system, i.e. the dynamic system Whose driVing force can not be simply represented by the gradient of a potential function [29, 31]. To better understand the quasi-potential, let us consider a special case. We suppose the dynamics is simply a gradient system With a single-well potential driven by small noise, i.e. For the Gillespie jump processes, there is no explicit form of 8(x). However, the invariant distribution P(x) satisfies the chemical master equation Where the Hamiltonian has the form for the standard Gillespie jump process. From classical mechanics, the S(x) obtained from WKB ansatz is exactly the quasi-potential defined through Equation (2). So even in the non-gradient case, we can still define S(x) as a generalization of the potential. That is Why it is called quasi-potential. S(x) is also a Lyapunov function of the original deterministic system [34] Equation (1a). The quasi-potential inherits key properties as the real potential guarantees in a gradient system. Besides logarithmically equivalent to the invariant distribution S(x) ~ —6 ln P(x), the mean eXit time 1‘ that the system escapes from an attractive basin has the asymptotic form 1‘ oc exp(VAS), where AS is the energy barrier height between the boundary of the basin and the stable state. The deeper the quasi-potential well is, the harder the system leaves a stable state. So the quasi-potential energy landscape describes the robustness of a system. This is similar for the extrinsic noise. For more properties about the quasi-potential, one may refer to [29—31] and 81 Text. In later text we will call S(x) the potential energy for simplicity.

Since we are interested in both the energy landscape and the transition path, we choose to compute the energy landscape through the gMAM method [31, 37]. The idea is to directly minimize the action functional Equation (2) through Maupertuis principle for the space of curves (See SI for details). In our work, the constructed energy landscape S(x) is a function of three variables. For convenience of visualization and analysis, we plot S(x) in two dimensional planes with a certain dimension f1xed. Therefore the global landscape is cut into different slices from various directions.

Results

Energy landscape of the yeast cell cycle network

Using the described method, we constructed the energy landscape S(x) for a budding yeast cell cycle network. We will state our findings from the energy landscape along the evolving path, i.e., from the excited G1 state (P2) to the final steady G1 state (P1).

Fig. 2A illustrates the slice of the energy landscape on the x-z plane where y = 0. We can see that the G1 state is the global minimum on the energy landscape, and there exists an energy barrier between P1 and P2 to prevent small noise activation. Outside this potential well, the energy function S(x) along the first part of the trajectory (from P2 to P3) is relatively flat, while the energy cost to deviate from the cycling path is high. Intuitively, we will call the cell cycle trajectory a “canal” to illustrate its flatness along the path.

The S phase canal is quite narrow when it evolves near the vertex P3. As the system evolves through P3, x gradually triggers the activation of the early M phase variable y, at which point the activated y begins to repress x. This corresponds to the S/M transition of the yeast cell cycle and we denote it as the early M phase canal. The landscape of the S phase and S/ M transition is illustrated in Fig. 2D where z = O. In the bottom right corner of Fig. 2D, the energy barrier between the canals of the S and early M phases greatly decreases the probability that the system passes the S/M transition without crossing P3, hence ensuring the robustness of the S/M transition.

In Fig. 2B, the system is shown to be temporarily restricted to a small area on the x-y plane, isolated by barriers separating it from the G1 state and the early M canal. As 2 gradually decreases to 0.05, Fig. 2C shows that this restricting area is still small but slowly shifts to a place with a larger x value. In addition, the early M canal looks more apparent. When 2 finally falls to zero (Fig. 2D), the energy barrier between the S canal and the early M canal disappears. Now the cell can execute DNA replication events with a long duration across P3, after which it successively evolves to the M phase.

When the system passes through P4, y triggers the activation of the late M variable 2, and the activated 2 begins to repress y. This is the transition from the early M phase to the late M phase, corresponding to the metaphase/ana-phase transition in yeast cell cycle process. The landscape on the x = 0 plane looks very similar

2D) and is shown in S6 Fig. Finally, in the third part of the cell cycle, the activated 2 represses y to zero and the system evolves to a G1 stable state (P1) and waits for another cell cycle division signal.

Non-gradient force and pseudo energy landscape

Since the main canal along the cell cycle trajectory in Equation (1a) after P2 is flat, the above energy landscape itself cannot tell us the moving direction of the system on the canal, Which impels us to investigate its non-gradient nature. From large deviation theory, we know that the most probable cycling path under Gillespie’s stochastic jump dynamics satisfies where S(x) is the energy landscape under discussion, 1; E V8, and H is the Hamiltonian of the considered system. In general, F(x) is not parallel to V8, and so we have additional non-gradi-ent effects for the transition paths. This point is also emphasized in [8, 9]. However, if the gradient force becomes zero, we have F(x) 2 b(x), which exactly corresponds to the dynamical driving force in the deterministic Equation (1a) (See SI Text for details). This is the case for the flat landscape along the canal in our results.

Along the cell cycle trajectory, we observe that the force strength F(x) is large in both S phase and early M phase canals, but extremely small near P3 and P4. Furthermore, we find that this tendency basically corresponds to the restriction width of the canal in the transverse direction (Fig. 3B). This means that, in the S phase and early M phase canals where the driving force is large, the cell cycle process will progress quite quickly, and hence does not require a strong restriction. As the process passes through P3 and P4, however, the system evolves more slowly (i.e. with more duration), and the driving force decreases and the canal width narrows in order to restrict fluctuations in the transverse direction. Thus, the dynamic properties around P3 work as an analogous DNA replication checkpoint. Similarly, the same characteristics around point P4 act as an analogous M phase checkpoint. Therefore, we suggest that this fast-slow dynamic and corresponding wide-narrow geometry of the landscape act as the dynamical mechanism keeping the cell cycle robust, precise and efficient.

The main idea behind this approach is to temporarily remove the globally stable state P1 from our model and only focus on the downhill flat canal after the saddle point P2. Thus the pseudo energy landscape is only a local landscape and no longer reflects the global stationary probability distribution (See SI Text for details). In Fig. 3C and 3D we illustrate the pseudo energy landscape constructed using this method. The result is a bit like combining the original landscape and the non-gradient effect together. From Fig. 3C, we can see that the effect of the force strength is replaced by the steepness of the canal in the tangential direction, while the landscape in the transverse direction remains the same. We emphasize that this point is essential to explain the directionality of the dynamic path along the flat canal, where we have observed a phenomenological fact that the gradient of the global potential S(x) becomes zero and the driving force b(x) gives the moving direction. In Fig. 3D we also show that the pseudo landscape on the z = 0.3 plane. Compared with Fig. 2B, we can see that the appearance of the S and early M canals is more apparent in the pseudo landscape. Furthermore, there exists a barrier between the S and early M canals (in the bottom right corner of Fig. 3D), which further guarantees the complete disappearance of the specie 2 before the cell enters the M phase.

Energy landscape adaptively responds to external signals

In the previous results, we obtained and analyzed the energy landscape of the yeast cell cycle model in Equation (1a) using a0 = 0.001, kal = 0.001, kaz = 0.001, and so on. However, When the external signal or stress changes, how does the energy landscape adaptively make such a transformation? First, let us discuss the energy landscape’s response to DNA damage in the S phase of the yeast cell cycle. When there is DNA damage in the S phase, the DNA replication checkpoint is activated [38]. In the cell cycle process, the checkpoint mechanism ensures the completion of early events before the beginning of later events so as to maintain the progression order of the cell cycle [18]. We can decrease the parameter kal to 0.0001 to simulate this effect (See 81 Text for more details), with the resulting new energy landscape on the z = 0 plane shown in Fig. 4A and 4B. The results show that the flat canal around P3 now turns into a small pit, which will keep the system in the P3 state until the DNA damage is repaired. Similarly, kaz = 0.0001 can be used to simulate the M phase checkpoint.

The resulting energy landscape is shown in Fig. 4C, Where the global three-dimension energy landscape is projected onto the y-z plane by choosing the potential minimum With respect to x for fixed y and z. The results demonstrate that the system shifts from an excitable system to a stable limit cycle system through bifurcation (a0 = 0.0025) When nutrition is sufficient, and the stable G1 state disappears. As a comparison, we show a similar projected energy landscape With a0 = 0.001 in Fig. 4D, Where the cell cycle process is an excitable system.

Finite volume effect

In other words, we only studied the fluctuation effects for perturbations by noises that are small compared to the reaction rates I) (x) in Equation (1a). This is a reasonable assumption for most biological systems. However, if the system nearly stagnates at some transition areas where its reaction rates are so small that they are of the same magnitude as the noise, the previous picture does not hold. In this case, some special phenomena will appear that do not fall into our traditional analysis.

From Fig. 5 we can see that the landscape is generally unchanged except around two corners. The original flat landscape with long-lasting slow dynamical behavior near the checkpoints now turns into pits. In other words, these nearly degenerate points play the role of small potential wells along the canal. This is intuitively true because the stationary probability for each point on this unidirectional path is approximately determined by the speed with which the system passes across. Consequently the system transitions to having a multipeak probability distribution, where the additional peaks do not correspond to the stable points of ODEs in the usual case. These additional peaks are also found in the cell cycle process of mammalian cells [13]. This interesting finite volume effect is unique when our driving force strength is comparable to the noise strength in some places. The additional pits can be treated as a “finger-print” for a multistage biological processes, in Which the pits act as another kind of analogous checkpoint. We speculate that a similar phenomena also holds for the classical limit

cycle dynamics.

Intrinsic versus Extrinsic noise

However, we also investigated the energy landscape of our cell cycle model perturbed by extrinsic noise using a similar approach, and the result is shown in fig. Compared with the landscape obtained in M and fi, the general shape of the landscape is almost the same, but the width of the canal does not change as significantly as the one perturbed by intrinsic noise. This result coincides with the lower insensitive fluctuation strength in the extrinsic noise case (82 Fig). This point, which indicates our cell cycle model is more tolerant with respect to intrinsic noise, can be used to help distinguish between intrinsic randomness and any environmental perturbations of the system, especially for a multistage biological process that periodically changes its reaction rates in time.

Discussion

The energy landscape of the budding yeast cell cycle is mainly comprised of two parts on a global scale: a deep pit that holds a cell in its G1 state when the environment is not suitable for division, and one unidirectional flat canal that performs a robust and accurate cell cycle progress once the system is excited (as summarized in Fig. 6A). When nutrients are sufficient, the original excitable system evolves into a stable limit cycle. Correspondingly, the deep pit lifts up, and the cell cycle can proceed efficiently without waiting in the G1 state (Fig. 6B). Once a cell meets any accident during its division, the corresponding cell cycle checkpoint is activated. In this case, the flat canal is dug to form an additional pit, and the system is held there until the accident is resolved (Fig. 6C). Furthermore, the super-slow dynamic stage on the canal corresponds to small pits if we take into account the finite volume effect. Those pits reduce the fluctuations from the cycle’s upstream processes and provides a longer stay duration when the system passes by (Fig. 6D).

The unidirectional canal and the guardrail on both sides guarantee that each event occurs only once and in the right order. Although the global energy landscape defined through the stationary probability distribution does not contain the non-gradient effect, the local pseudo energy landscape we proposed clearly visualizes this unidirectionality brought by the non-gradient force. Now the strength of the non-gradient force along the canal is characterized by the steepness of the pseudo energy in the tangential direction. For a system perturbed by intrinsic noise, the fluctuation restriction ability in the transverse direction of the canal is highly related to the driving force strength at that point. However, for a system perturbed by extrinsic noise, this relationship is much weaker.

In our model, only the G1 phase cyclins are controlled by the cell mass, and the mitosis event in the M phase is triggered by the completion of the DNA replication event. This “domino” mechanism is also found in the cell cycle regulation of higher eukaryotic organisms [43]. With our model and new methods, we clearly identified the DNA replication and M phase checkpoints in the constructed energy landscape. This landscape can reshape itself in response to environmental nutrients (similar and consistent results in mammalian cell cycle [13]) and checkpoint signals adaptively. Furthermore, we proposed the concept local pseudo energy landscape to characterize the irreversibility of the dynamic path along the flat canal in the landscape. These points, to the authors’ knowledge, have not been revealed in the previous studies.

When the problem is considered in higher dimensions, the computational cost increases exponentially. Even if this computational cost issue is resolved, however, determining how to save and eXhibit such high-dimensional information remains a tough question. There is still a need to develop a systematic reduction method to analyze high dimensional problems. Even with such limitations, we believe our meticulous study of the energy landscape of the simplified budding yeast cell cycle model is of general interest to those studying other complicated cell cycle dynamics. Most of the insights we gained studying this simple model are independent of the number of dimensions and the specific formulation of the model, and therefore will be valuable to other systems and studies.

Supporting Information

This file contains details that needed to understand the main body. It is arranged as follows: I. The three-node Budding Yeast Cell Cycle Model, 11. Stochastic Model, 111. Large Deviation Theory and the Hamiltonian, IV. Overview of the Construction of the Landscape, V. Introduction of the gMAM, VI. Construction of the Quasi-potential energy landscape, VII. Force Strength and Landscape Canal Width, VIII. Analysis of the Non-gradient Force, IX. Signals Regulation in the Cell Cycle model, X. Extrinsic and Intrinsic Noise, XI. Parameters Sensitivity Analysis.

(A) The regulatory network of key regulators in budding yeast cell-cycle process. It can be separated into Gl/ S, early M and late M modules, Where the nodes represent cyclins, transcriptional factors and inhibitors, and the green and red lines represent activation (transcription) and the inhibition, respectively. (B) The essential network of yeast cell-cycle network, Where X, Y and Z represent key regulators of the G1 /S, early M and late M modules respectively, different modules are connected by activation and inhibition interactions. More details can be found in the main text [44].

(A) The quasi-potential energy landscape under extrinsic noise perturbation, the x-y plane Where z = O. (B) The ODE driving flux strength (black dashed line) in the ODE path and the fluctuation strength in the vertical direction of the ODE path perturbed by intrinsic noise (red solid line) and extrinsic noise (blue solid line). The letters ‘a’, ‘b’ and ‘c’ mark three points on the ODE trajectory With

(A) The energy landscape in x-y plane With z = O. (B) The pseudo energy landscape on x-y plane With z = O. The red dashed boxes mark the energy barriers that contradict With the deter-S4 Fig. Parameters sensitivity analysis of the system. Effects of parameter perturbations on global evolution trajectory (AK) and the depth of energy well corresponding to the stable G1 state (L). In each subfigure, we either increase (blue solid line) or decrease (red dashed line) the value of a certain parameter by 20%. Subfigures from A to K, and the numbers on x-aXis in sub-figure (L) correspond to parameters: j1,j2,j3, k1, k2, k3, ki, k5, kal, kaz, 610, respectively. (L) AS denotes the depth of energy well corresponding to the stable G1 state with wild-type parameters, and 6S denotes the changes of depth after each parameter is increased or decreased.

In the second part of the cell cycle, the system evolves through the early M flat canal to the vertex P4 over a sufficient duration for the mitosis event. When the system passes through P4, y triggers the activation of the late M variable 2, and the activated 2 begins to repress y. This is the transition from the early M phase to the late M phase. The “early M” and “late M” in bold refer to early M phase and late M phase respectively. (EPS) Acknowledgments The authors are grateful to Weinan E and Qi Ouyang for helpful discussions.

Author Contributions

Performed the experiments: CL XL. Analyzed the data: CL XL FL TL. Contributed reagents/materials/ analysis tools: CL XL FL TL. Wrote the paper: CL XL FL TL.