A globally attractive cycle driven by sequential bifurcations containing ghost effects in a 3-node yeast cell cycle model
Fangting Li, Mingyang Hu, Bo Zhao, Hao Yan, Bin Wu, Qi Ouyang
11 A globally attractive cycle driven by sequential bifurcationscontaining ghost effects in a 3-node yeast cell cycle model
Fangting Li , , ∗ , Mingyang Hu , , Bo Zhao , Hao Yan , Bin Wu , Qi Ouyang , , ∗ ∗ E-mail: li [email protected], [email protected]
Abstract
Yeast cells produce daughter cells through a DNA replication and mitosis cycle associated with check-points and governed by the cell cycle regulatory network. To ensure genome stability and genetic informa-tion inheritance, this regulatory network must be dynamically robust against various fluctuations. Herewe construct a simplified cell cycle model for a budding yeast to investigate the underlying mechanismthat ensures robustness in this process containing sequential tasks (DNA replication and mitosis). Wefirst establish a three-variable model and select a parameter set that qualitatively describes the yeastcell cycle process. Then, through nonlinear dynamic analysis, we demonstrate that the yeast cell cycleprocess is an excitable system driven by a sequence of saddle-node bifurcations with ghost effects. Wefurther show that the yeast cell cycle trajectory is globally attractive with modularity in both state andparameter space, while the convergent manifold provides a suitable control state for cell cycle checkpoints.These results not only highlight a regulatory mechanism for executing successive cell cycle processes, butalso provide a possible strategy for the synthetic network design of sequential-task processes.
Author Summary
A challenge in biology is to understand how complex regulatory networks in a cell execute differentbiological functions, such as the DNA replication event and mitosis event in eukaryotic cell cycle processes,orderly and reliably. We call the process as the sequential-task process. In this work, we establisha simplified toy model for a budding yeast cell cycle, and develop the analysis methods to describe theevolving process. Our results show that the cell cycle process is an excitable system with global robustnesstowards different initial states, where the phase transitions of the cell cycle are controlled by feedback inthe network and driven by bifurcations with a slow passage that we refer to as a ghost effect. These ghosteffect bifurcations are sufficient for providing both the longer duration of each event and the modularityin state/parameter space. The modularity and decoupling of the kinetic parameters in different eventsensures the stability of the system against parameter changes; thus small parameter changes in any givenevent should not influence the general character of the cell cycle process. We describe this as the wavetransition “domino” model of the cell cycle. This wave transition model not only provides an efficientand easily controlled strategy for the cell cycle, but also suggests a guide for building synthetic networksof sequential-task processes.
Introduction
It has been suggested that dynamic robustness and modularity are the important characteristics ofa cellular regulatory network for fulfilling complex biological processes. Recent progress in quantitativemodeling and experiments of some fundamental cellular processes, including the chemotaxis in E. coli [1,2]and the cell cycle process in yeast and oocyte [3–8], have provided a better systematic and quantitativeview of these biological processes. These studies have revealed the structural and dynamic propertiesof regulatory networks, such as the basic motifs and modules of network structure [9–11], the feedback a r X i v : . [ q - b i o . M N ] S e p loops that dynamically control the genetic switches and transitions [12–15], and the dynamic stabilityand robustness of such networks [1, 16].In the cell cycle process, the cell successively executes DNA replication and mitosis events usinga checkpoint mechanism. The dynamic regulatory mechanism in cell cycle processes was revealed bythe pioneering work of Tyson, Novak, and Ferrell [17–19]. The cell cycle process is considered to bea series of irreversible transitions from one state to another [20, 21], which is regulated by system-levelfeedback [21]. For example, the cell cycle process in oocyte maturation is built upon a hysteretic switchwith two saddle-node bifurcations [8]. The skeleton models of the cell cycle in a budding yeast [22]and mammalian cells [23] have also been established to investigate the temporal order and sequentialactivation in the cell cycle. Furthermore, a DNA replication checkpoint in a mammalian cell cycle hasbeen proposed to lengthen the cell cycle period and lead to a better separation of the S and M phases [24].To reveal the essential dynamic properties of a sequential-task cell cycle process, some fundamentalquestions must be studied. For instance, how can the cell cycle network ensure the cell reliably executes asuccessive order of events? How does the network ensure that the process is robust against environmentalfluctuations? What is the function of the cell cycle checkpoint from a dynamic view? Here, we addressthese questions using a simplified cell cycle network model in a budding yeast Saccharomyces cerevisiae .We first develop a simplified yeast cell cycle model and select a set of parameters that can qualitativelydescribe the cell cycle process. We then apply perturbation methods to analyze the dynamic properties ofthe local manifold along the evolving cell cycle trajectory. We demonstrate the following: that the yeastcell cycle process is an excitable system driven by sequential saddle-node bifurcations; the yeast cell cycletrajectory is globally attractive with modularity in both state space and parameter space; and that theconvergent manifold of the trajectory corresponds to the cell cycle checkpoints. With these advantages,the yeast cell cycle process becomes effective and reliable throughout its execution.
Model
The budding yeast cell cycle process is a well-established system for investigating and exploring the controlmechanisms of regulatory networks [3–5,25,26]. Recent works integrating modeling and experiment havehighlighted the function of positive and negative feedback in the Start point [6,15,27], spindle separationprocess [28], and mitotic exit [29, 30]. To investigate the essential dynamic features of budding yeastcell cycle, especially the dynamic robustness and modularity, we construct and analyze a coarse-grainingmodel that is mainly based on the abstract architectures and functions of cell cycle and ignored themolecular details. These simplification approaches have been applied widely to study the design principlesof regulatory networks [16, 31, 32]
The simplified yeast cell cycle network and model
The cell cycle process in a budding yeast produces new daughter cells through two major events: DNAreplication in the S phase and mitosis in the M phase. We denote this kind of processes that executedifferent biological functions and events in order and reliably as the sequential-task processes. The dif-ferent cyclins and transcription factors (TF) control different successive events, and the DNA replicationcheckpoint and spindle assembly and separation checkpoint ensure the stability and hereditary validity ofthe genetic information [20]. The checkpoints ensure the completion of earlier events before the beginningof later events to maintain the orderly progression of the cell cycle [33].Based on the cell cycle regulatory network (Figure 3-34 in [20]), our previous network with an intra-S and spindle checkpoint [5], and the full cell cycle network in the supplemental information (SI), weseparated the yeast cell cycle regulatory network into the G1/S, early M and late M modules, where thepositive and negative feedback loops play an essential role in governing the yeast cell cycle process. Wethen obtained a coarse-graining simplified 3-module network, shown in Figure 1. In the figure, x representsthe activities of key regulators such as the cyclins Cln2 and Clb5, and SBF and MBF in the G1 and Sphases; y represents the activities of key regulators such as the cyclin Clb2 and the transcriptional factorMcm1/SFF in the early M phase; and z represent the activities of key inhibitors such as Cdh1, Cdc20,and Sic1 in the late M/G1 phase. For simplification, we ignore the G2 phase.Then a simplified three-variable ordinary differential equation (ODE) model can be established. The 3-node cell cycle model is not only reduced on the topological structure of budding yeast cell cycle network,but also constructed to describe the key dynamic features of yeast cell cycle, especially the geneticswitches, irreversible transitions, and the sequential-task process from DNA replication to mitosis. The3-node model is written as follows: d x d t = x j + x − k x − yx, (1a)d y d t = k a x + y j + y − k y − zy, (1b)d z d t = k a y + k s z j + z − k z − k i xz. (1c)In this 3-node model, we assume that the simplest forms to represent the interactions in the cell cyclenetwork in order to qualitatively describe the genetic switches in the START point, G2/M transitionand the late M phase, and we introduce the checkpoint mechanism by changing kinetic parameters.The positive feedback in each module is represented by a 2-order Hill function, while the inhibition andrepression component among different modules is assumed to have its simplest form, i.e., − yx , − zy , and − xz . The parameters j i and k i are the Hill coefficients and degradation rates of x , y , and z respectively; k a x is the activation rate from x to y ; and k a y is the activation rate from y to z . Furthermore, the − k i xz term in Eq. 1c represents the repression from x to z , − zy in Eq. 1b represents the repression from z to y , and − yx in Eq. 1a represents the repression from y to x , these are negative feedback. More detailconcerning network reduction, model assumptions, formulations, and the dimensionless form deductioncan be found in SI.We focus on the dynamic robustness and modularity of the cell cycle process, particularly from theexcited G1 state to the final steady G1 state. Before we quantitatively analyze the model, we qualitativelyillustrate the dynamics of x , y , and z and the corresponding yeast cell cycle events and process. At thebeginning of the cell cycle process in the G1 phase state, inhibitor z is at a high level with low levels of x and y , and the cell is in the resting G1 state. Once the system enters the cell cycle process, x is activatedand dominant, which represses the inhibitor z to a lower level. In the early M phase, y is activated anddominant, and the activated y represses x in turn. In the late M and G1 phases, a z wave is triggeredby y and increases to high levels. Sequentially, the x wave is dominant in the S phase, the y wave isdominant in the early M phase, and the z wave is dominant in the late M and G1 phases. We denotethis set of sequential steps as the wave transition “domino” model of the yeast cell cycle process.The excited G1 state of the yeast cell cycle process can be settled as follows. If we set k a = k a = 0and ignore the repression terms − xy , − yz , and − k i xz in Eq. 1, the equations for x , y , and z can bedecoupled. In this situation, the steady state of x in equation dxdt = x / ( j + x ) − k x ≡ G ( x, j , k ) = 0has three fixed points, i.e., x (1) = 0 , x (2) = − √ − j k k , and x (3) = √ − j k k , where x (1) and x (3) are stable nodes and x (2) is an unstable saddle that serves as the threshold to the excited G1 state andseparates the attractive ranges of x (1) and x (3) . These dynamic features are only determined by j and k .Similarly, if we set H ( y, j , k ) ≡ y j + y − k y = 0 and I ( z, k s , j , k ) ≡ k s z j + z − k z = 0, we have y (1) =0 , y (2) = − √ − j k k , y (3) = √ − j k k , z (1) = 0 , z (2) = k s − √ k s − j k k , and z (3) = k s + √ k s − j k k . Theexcited G1 state is defined as ( x ≥ x (2) , y = y (1) , z = z (3) ). Constraint conditions and suitable kinetic parameters
In the yeast cell cycle process, the cell should finish the DNA replication before mitosis, and there shouldbe enough time to execute both the DNA replication and mitosis events. The question then arises as towhat criteria are available to qualitatively describe the yeast cell cycle process. We use an event orderand the long duration of the x and y waves as constraints to search for suitable kinetic parameter setsin Eq. 1. The first constraint is that the x wave is followed by the y wave, while the z wave decreases inthe S phase and increases in the late M and G1 phases. The second constraint is that the duration times(calculated as the full duration at half maximum) of both the x and y waves should be sufficiently long(i.e., ≥
30 arbitrary units (au)) in order to ensure the successful execution of any DNA replication andmitosis events.Using the 10 random parameter sets we produced by Latin hypercube sampling [34] and the excitedG1 state as initial state, we obtained 239 groups of parameter sets in a simulation iteration that satisfiedthe above conditions. These parameter sets were clustered as shown in Figure 2A. We found that theabove constraint conditions only required small activation rate values for k a (from the x wave to the y wave) and k a (from the y wave to the z wave). Furthermore, k s and k —which are related to z —mustbe correlated, while j , j , and j are relatively large in value. Different simulation iterations showedsimilar results.We then set j m = 0 . k m = 0 . m = 1 , ,
3) to search for the satisfied parameters k a , k a , k i ,and k s . The results are plotted in Figure 2B, with the x wave duration T x as a function of k a . Theenvelope curve corresponds to T x ∝ ( k a − k ca ) − / , where k ca will be discussed in the following section.Furthermore, a suitable value of k a is 0 .
001 for the constraint T x > au . For the same reason, weselect k a = 0 . j m = 0 . k m = 0 . m = 1 , , k a = k a = 0 . k i = 5 .
0, and k s = 1 .
0. We denotethis set as the perfect yeast cell cycle parameter set.
Duration of G1/S and early M phases is controlled by strength of activationthrough the ghost effect
What is the relationship between the duration of the x or y wave and k a or k a ? To illustrate thisrelationship, we analyze how the x wave triggers the activation of the y wave. During the late S phaseand early M phase, x is almost fully activated and has repressed z to zero, so we can ignore z in Eq. 1to get dxdt = x j + x − k x − xy = G ( x, j , k ) − xy, (2a) dydt = k a x + y j + y − k y = k a x + H ( y, j , k ) . (2b)We showed in the previous section that G ( x, j , k ) and H ( y, j , k ) determine x ( i ) and y ( i ) ( i = 1 , , k a x in Eq. 2b will trigger the activation of y , and the repression term − xy in Eq. 2awill repress x . In the early M phase, k a x triggers the activation of y through a saddle-node bifurcation,where two smaller fixed points of y ( y (1) and y (2) ) collide and vanish, leaving only y (3) . We denote k ca asthe critical bifurcation point of k a , and its value can be estimated as follows: H ( y, j , k ) ≡ y j + y − k y is defined as part of the right hand side terms of Eq. 2b, while H ( y ∗ ) is the minimum of H ( y, j , k ). Wecan ignore the inhibition term − yx in Eq. 2a in the late S phase because x is approaching its maximum x (3) and y is not fully activated. When k a x (3) ≥ | H ( y ∗ ) | or k a ≥ k ca ≡ | H ( y ∗ ) | /x (3) , the saddle-nodebifurcation will occur. Note that x (3) is only determined by j and k , while H ( y ∗ ) is determined by j and k . Thus, k ca is determined by j , k , j , and k .Furthermore, if 0 < ( k a − k ca ) (cid:28)
1, then just after the bifurcation we have x (cid:39) x (3) , y (cid:39) y ∗ (cid:39) y (1) ,and z (cid:39)
0, so d y d t = k a x + H ( y, j , k ) (cid:39) ( k a − k ca ) x (3) + c ( y − y (1) ) (cid:28) . It is just after the two smaller fixed points collide that the saddle-node remnant, or ghost effect, witha slow passage will be observed [35]. We estimate the duration of the x wave as a function of k a bycalculating y from zero to + ∞ , i.e., T x = (cid:90) ∞ dy ( k a − k ca ) x (3) + c ( y − y (1) ) = π (cid:112) c ( k a − k ca ) x (3) . Considering the inhibition term − xy in Eq. 2a, if y is near y (2) , then x can reach x (3) (cid:48) = √ − j k (cid:48) k ,where k (cid:48) = k + − √ − j k k . When k a ≥ k ca = | H ( y ∗ ) | /x (3) (cid:48) , y can be triggered and activated, while k ca is mainly determined by parameters j , k , j , and k .The envelope curve of Figure 2B is based on this analysis, and is consistent with our simulations.Using the selected parameter set, we calculated the value of k ca based on the above estimation, andfound k ca = 5 . × − , (i.e., very close to 0 . y is activated and then evolves to y (3) , and the repression term − xy in Eq. 2amakes x (3) unstable and represses x to zero. A similar mechanism occurs in the late M phase, wherethe activation of z is triggered by k a y , and the long duration of the y wave is caused by the saddle-node remnant, i.e., the ghost effect. If we ignore the inhibition term − zy from z to y in Eq. 1b, then k ca = | I ( z ∗ ) | /y (3) and I ( z ∗ ) is the minimum of the I ( z ) curve determined by k s , j , and k , while y (3) isdetermined by j and k . The duration of the y wave is determined by T y ∝ ( k a − k ca ) − / , and k ca isdetermined by k s , j , k , j , and k . Again, using the selected parameter set, we calculated the value of k ca and found k ca = 5 . × − , a value very close to 0 . Methods
In this section, we develop methods to illustrate local stability along an evolving dynamic trajectory,then apply those methods to the yeast cell cycle process. Using the perfect parameter set, the perfectcell cycle trajectory evolving from the excited G1 state to the stable G1 state is shown in Figure 3A andB. The figure also shows x ( s ), y ( s ), and z ( s ) in the upper panels of Figure 4B, where s is the trajectorylength from the initial state. Note that V ( t ) and V ( s ) are the velocity of the trajectory as a function oftime t and of s , respectively. The normal plane of the trajectory is defined as the plane that is alwaysperpendicular to V . Local spherical surface perturbation
Points on a small spherical surface of radius r around the exited G1 state were taken as perturbed initialstates. In Figure 4 and Figure 5, we set r = 0 .
01. This sphere of perturbed states evolves together alongthe perfect trajectory with varying deviations, leading to a change in the shape of the sphere. At eachtime t or trajectory length s , the deviation of each point on the perturbed sphere from the perfect stateat the center was calculated as r i ( t ) or r i ( s ), where i is the index of the perturbed trajectory. The normalradius r n of the perturbed sphere is the average of all r i for each value of i on the cross section of thesphere that goes through the sphere center and is normal to V . The forward and backward tangentialradius r t and r t are defined as the farthest distance on the perturbed sphere from the center, in thedirection parallel and antiparallel to the velocity V , respectively. As an example, the normal radius r n ( t )and r n ( s ) of the perfect yeast cell cycle trajectory are plotted in the middle panels in Figure 4. Local velocity analysis
In the above method, the manifolds are formed by the trajectories evolving from the given small sphericalsurface at t = 0 or s = 0. In order to fairly examine the system’s dynamic properties at different t and s , circular perturbations on the normal plane of the yeast cell cycle trajectory were added at each stepof the trajectory. The velocity vector of a perturbed trajectory can be decomposed with respect to the V axis into the normal component V ⊥ and the tangential/parallel component V (cid:107) . The average velocitieson the perturbed circle were then calculated. The magnitude of the average velocities V indicates howfast the manifold evolves, while the average of the normal components V ⊥ and tangential components V (cid:107) denote how the manifolds converge and disperse on average. The V ⊥ ( t ) and V ⊥ ( s ) of the perfect cellcycle trajectory are plotted in the middle panels of Figure 4 as an example. Eigenvalues of the normal Jacobian matrix
The stability of any fixed point can be examined through the eigenvalues and eigenvectors of the Jacobianmatrix of the ODEs. For a 3-dimensional dynamic system, dx i /dt = f i ( x j ), i = 1 , ,
3, the Jacobianmatrix is defined as ∂f ∂x ∂f ∂x ∂f ∂x ∂f ∂x ∂f ∂x ∂f ∂x ∂f ∂x ∂f ∂x ∂f ∂x If any real part of the eigenvalue of the Jacobian matrix is positive, the fixed point will be unstable inthe direction of the correlative eigenvector. In other words, a fixed point remains stable only when all ofits eigenvalues are negative.Along the evolving trajectory, all state points should evolve in the velocity direction except for the finalfixed point. To analyze the convergence properties of each manifold along the yeast cell cycle trajectory,the Jacobian matrix should be projected onto the normal plane of the trajectory. We adopt two unitorthogonal vectors n and n , together with unit velocity vector V / | V | , to compose a set of orthogonalunit vectors. Then the projecting matrix P can be defined as ( n , n , J ⊥ is equal to P T JP , in the form j j j j
00 0 0 The four non-zero elements consist of a 2 × λ and λ andeigenvector dictate the convergence properties in the normal plane of the evolving trajectory.As an example, the real parts of λ and λ along the perfect cell cycle trajectory are shown in the lowerpart of Figure 4. The real parts are almost always negative except for the first part of the cell cycle andnear the two vertex in Figure 3B. At the two vertex, where x and y reach their respective maximums andbegin to decrease, there are positive peaks followed by large negative eigenvalues. Because the perfectcell cycle trajectory changes its evolving direction quickly near these vertices, we assume the positiveeigenvalue peaks are caused by a numerical error in the velocity directions at the vertices. Results
The yeast cell cycle trajectory driven by saddle-node bifurcations with ghosteffects is globally attractive
Using the perfect parameter set, we depict the perfect yeast cell cycle trajectory in Figure 3A and B, andits local dynamic analysis results in Figure 4. In Figure 4, we plot the normal radius r n to describe thedeviation from the perfect trajectory in the normal direction of the velocity of the initial state, the localaverage velocity V to depict the velocities of the circular perturbations along the evolving trajectory,and the real parts of eigenvalues λ and λ to describe the convergence properties in the normal planeof the trajectory. These results reveal some interesting dynamic properties along the perfect cell cycletrajectory.In the first part of the trajectory (G1 and S phases), Figure 4 shows that r n reaches a peak at t = 4 . s = 5 .
32, suggesting that the perfect trajectory diverges at the beginning of the S phase. When thesystem evolves at t = 20 .
09 and s = 8 .
59 to the first vertex ( x max = 4 . y = 0 . z = 0 .
00) inFigure 3B, x is fully activated and has repressed z to zero, and x also begins to trigger the activationof y by saddle-node bifurcation. The duration of the x wave is given by T x ∝ ( k a − k ca ) − / . Because( k a − k ca ) (cid:28) x wave. Theresults in Figure 4 show that, near the first vertex, the normal radius r n decreases to near zero, the realparts of eigenvalues λ and λ sharply decrease to negative values, and the average local velocity V ≈ x wave to execute a DNA replication event, whilethe converging manifold and slowly evolving trajectory near the first vertex provides a suitable state forthe DNA replication checkpoint. On the other hand, the vertex with a convergent manifold works as ahidden intrinsic checkpoint mechanism for the cell cycle process.A similar dynamic property is found in the later part of the perfect cell cycle trajectory. When t = 39 .
16 and s = 10 .
78 in the early M phase, the trajectory diverges. Near the second vertex (0, y max = 4 .
35, 0), y is fully activated and has repressed x to zero, so y is triggering the activation of z through saddle-node bifurcation. This saddle-node bifurcation and ghost effect lead to a converging,attractive and slowly evolving manifold near the second vertex. The second vertex corresponds to thespindle assembly and separation event. Modularity of state and parameter space in the perfect cell cycle model
The modularity of both state and parameter space can be found in the perfect cell cycle model. In Figure3A and B, not only do the x , y , and z waves have little overlap time, the three vertices in the phasemap have sharp pointed ends and the kinetic parameters in the x , y , and z waves are also decoupled intoindependent groups.Near the first vertex, k ca = | H ( y ∗ ) | /x (3) , where H ( y ∗ ) is the minimum of H ( y, j , k ) and x (3) = √ − j k k . So the first bifurcation near the first vertex ( x max , 0, 0) is mainly controlled by the y wave’skinetic parameters j and k through H ( y ∗ ), and the x wave’s kinetic parameters j and k through x (3) . Note that these values are almost independent of the z wave’s kinetic parameters because the eventoccurs mainly in the x − y plane with z = 0. Similarly near the second vertex, k ca = | I ( z ∗ ) | /y (3) , where I ( z ∗ ) is the minimum of I ( z, j , k , k s ) and y (3) = √ − j k k . The second bifurcation near the secondvertex (0, y max , 0) is mainly in the y − z plane with x = 0, and is controlled by the z wave’s kineticparameters j and k , and by k s through I ( z ∗ ) and the y wave’s kinetic parameters j and k through y (3) . Again, these values are almost independent of the x wave’s kinetic parameters.When 0 < ( k a − k ca ) (cid:28) < ( k a − k ca ) (cid:28)
1, the relative saddle-node bifurcations occur insufficient order. The duration of x wave is mainly controlled by the ( k a − k ca ) and the maximum of x wave x (3) shown in 2B. Similar is the y wave, in which the duration of y wave is mainly controlled by the( k a − k ca ) and the maximum of y wave y (3) . In this case, the system not only provides a long durationfor the execution of the DNA replication and mitosis events, but also separates the x , y , and z waves intorelatively independent events and decouples the kinetic parameters of different waves and phases. Theseresults provide a suitable strategy for the regulation of sequential events, especially because the changeand modification parameters of one event influence other events very little. The perfect cell cycle model and its trajectory
The robustness and manifold analysis along the cell cycle trajectory provide more dynamic informationof the cell cycle process. Novak and coworkers discussed a kinetic differential equation model in which thebudding yeast cell cycle trajectory could be separated into excitation and relaxation periods accordingto the eigenvalues of the Jacobian of the trajectory [36]. The divergence and convergence manifold hasalso been observed in a Boolean network model, where the manifold converges at the S phase statecorresponding to a DNA replication checkpoint [5]. In the perfect cell cycle model we find some newresults, including the manifold of key regulators diverges and converges in different cell cycle phases, andthe dynamic robustness and modularity of cell cycle process.Let us interpret the results from the biological perspective. Suppose that the states (activities of thekey regulators) of yeast cells can be dynamically observed, if a group of yeast cells start from differentexcited G1 states with fluctuations of biochemical parameters that tend to vary from cell to cell, how dothe yeast cells evolve during the whole cell cycle process? During the early S phase, the states of theyeast cells should be observed to separate and diverge significantly (diverge manifold) and their states arechanging quickly (large average velocity), and in this phase the yeast cells activate kinases and expressgenes for the execution of DNA synthesis and replication task. When the cells enter the late S phase thatthe DNA replication task is working, the yeast cell states begin to change slowly and gradually convergeto the same point; this corresponds to the system evolves into the first vertex ( x max , , , y max , , , z max ) waiting for another cell cycle signal. During the whole cell cycle process, fromthe DNA replication event to the mitosis, the dynamic robustness and modularity of the system provide asuitable strategy to execute the sequential events reliably. So in cell cycle process, the manifolds divergeand converge, wave after wave, cycle in cycle out. This is the dynamical picture of cell cycle.The durations of the x and y waves in the perfect cell cycle model are controlled by, and sensitive to, k a and k a , respectively. In the SI we investigate a simplified cell cycle model with inhibitor I , in which I represses y and x represses I by one-step or multi-step phosphorylation, noted as x (cid:97) I (cid:97) y . In theSI, we show that the inhibitor with phosphorylation can largely reduce the duration sensitivity of the x wave T x to changes in the relative kinetic parameters from the x to inhibitor I , and the x wave duration T x is not sensitive to the maximum of x wave x (3) . These differences can be also used to distinguish thedirect activation and indirect activation (with inhibitor) between the successive waves. The “ideal” yeast cell cycle trajectory and the cell cycle checkpoints
In a real yeast cell cycle process, yeast cells produce new daughter cells through DNA replication inthe S phase and undergo mitosis in the M phase, with each event associated with a relative checkpointmechanism. For example, when the cell is in a DNA replication event, or the DNA is damaged in theS phase, the relative pathway is activated and turns on the DNA replication checkpoint to keep the cellin the S phase; only when the damaged DNA is repaired or the yeast cell finishes the DNA replicationevent is the DNA replication checkpoint turned off and the the yeast cell allowed to enter the G2 and Mphases.To simulate the above cell cycle checkpoint function, we artificially changed k a and k a to representthe on or off states of the cell cycle checkpoints. When the DNA replication checkpoint is turned on, k a = 0; when it is turned off, k a > k ca . Similar to the spindle checkpoint, k a = 0 represents the onstate while k a > k ca is the off state. We denote this as the real or “ideal” yeast cell cycle trajectory.The ideal trajectory evolves through the 3 vertex states ( x (3) , , , y (3) , , , z (3) ) in order,and the checkpoints not only control the duration of the x and y waves, but also separate the x , y ,and z waves. Furthermore, the checkpoint mechanisms make the kinetic parameters in different wavesindependent. This kind of extrinsic checkpoint mechanism in real yeast cell cycle process provides a morestable duration for each phase.The ideal yeast cell cycle trajectory is globally attractive and stable. The following discussion providesfurther evidence for this conclusion. If we were to ignore the activations and repressions among x , y , and z in Eq. 1 as formulated in the Model section of this work, there would be 3 fixed points, one on eachof the x , y , and z axis. Thus, x (1) = 0, x (3) , y (1) = 0 , y (3) , and z (1) = 0 , z (3) are stable nodes, while x (2) , y (2) , and z (2) are unstable saddles. Consider the repression terms − xy , − yz , and − k i xz in Eq. 1,where only one of x , y , and z can be in its maximum state, and (0 , ,
0) is unstable. In this case, only( x (3) , , , y (3) , , , z (3) ) are the possible final stable attractors.When discussing the evolution of the cell cycle engine, Andrew Murray suggested to make “toy sys-tems” that are much simpler and mimic the key features of cell cycle process for the deeper understandingof cell cycle [37]. From the intrinsic checkpoint mechanism in the perfect cell cycle to the outside check-point mechanism in real yeast cell cycle process, our results suggest a possible evolution course of cellcycle network and sequential-task regulatory networks. Imperfect yeast cell cycle trajectory
In this section we investigate a counter example, the imperfect yeast cell cycle process, where we set k a = 0 .
04 and k a = 0 .
04, and k a > k ca and k a > k ca , all of which are relatively far from the bifurcationpoints. The imperfect trajectory (Figure 3C and D) and its local dynamic properties (Figure 5) are quitedifferent from the above perfect yeast cell cycle trajectories. In Figure 3C and Figure 3D, the maxima of x and y just reach half of x (3) and y (3) , there is a relatively large overlap time for the x , y , and z waves,and there is no sharp pointed end near the 1st and 2nd vertices in phase space. The duration times ofthe x and y waves are only 1/6 that of the duration times of the perfect cell cycle trajectory. All of theseresults show that the x , y , and z waves are coupled together. In Figure 5, one observes that the expansionis less than the perfect trajectory ( r n ), even though the local manifold also expands and converges alongthe trajectory. Near the vertices in the phase map, the manifold cannot converge to a narrow statespace with a certain evolving velocity ( V ), which then cannot provide a long enough duration for DNAreplication and spindle assembly and separation events. Furthermore, there is no suitable point for thecheckpoint mechanism. Discussion and Conclusion
The yeast cell cycle process executes DNA replication and mitosis in subsequent order, and shouldbe robust against any intrinsic fluctuations and environmental change. In this paper, we constructedand investigated a simplified yeast cell cycle model as a toy model to capture the dynamic essence ofsequential-task biological processes. We have shown that the cell cycle process is an excitable system withglobal robustness in the state space, and that the phase transitions are driven by dynamic bifurcationscontaining ghost effects. These bifurcations and associated ghost effect sufficiently provide the longerdurations required for each event or wave, as well as the modularity of the state and parameter spaces.This cell cycle model is robust and modular in both state and parameter space. The robustness in statespace ensures that the cell cycle process is stable against fluctuations in protein activity. The modularityand decoupling of the kinetic parameters in the different waves avoids allowing small parameter changes inany one event to influence the general cell cycle process. Our wave transition model provides an efficientand easily controlled strategy for executing the cell cycle processes.0Furthermore, our toy cell cycle model suggests a possible synthetic network design for robustly ex-ecuting other sequential-task processes. In the synthetic network, each event can be controlled by therelevant key regulators, and the duration of each event is regulated by the activation strength betweensuccessive events (similar to k a x and k a y in our toy model). These synthetic networks could be used tocheck our prediction that the long duration and modularity of the sequential-task process is caused andcontrolled by the ghost effect.Most of the insights we obtained from the toy model are independent of the specific formulation ofthe model and the number of dimensions, and the sufficient condition is that the dynamic cell cyclemodels should be governed by sequential saddle-node bifurcations containing ghost effects. We believethat our results and analysis method are applicable to other eukaryotic cell cycle processes and othersequential-task biological processes. Life is modular in parts and more than the sum of its parts, andthis robustness and modularity are the design principles for biological regulatory networks [11, 38–40].Indeed, if a network behaves along a globally attractive trajectory as well as being modular in stateand parameter space, it can execute successive irreversible events more robustly and effectively and thusobtain an evolutionary advantage. Supporting Information
Text S1
This file contains details needed to understand the main body of this work, which are arrangedas follows:I The regulatory network in budding yeast cell-cycle process and a 3-node modelII A simplified cell-cycle model with inhibitor
Figure S1
The regulatory network of key regulators in the budding yeast cell-cycle process, whichcan be separated into the G1/S, early M, and late M modules. It is contained within Section I of the SI.
Figure S2
The network of the one-step phosphorylation inhibitor that is inserted between x and y in the 3-node yeast cell-cycle model. It is contained within Section II of the SI. Figure S3
An inhibitor with 4-step phosphorylation is added between x and y in the yeast cell-cyclemodel, which decreases the sensitivity of wave duration to changes in the triggering rate b . It is containedwithin Section II of the SI. Acknowledgments
The authors are grateful to Chao Tang, Louis Tao, Tiejun Li, Xiaomeng Zhang, Tianqi Liu, and YuhangHou for helpful discussions. The work is supported by NSFC grants nos. 11174011, 11021463, and91130005 (F.Li), and nos. 11074009 and 10721463 (Q.Ouyang).
References
1. Barkai N, Leibler S (1997) Robustness in simple biochemical networks. Nature 387: 913-917.2. Hansen C, Endres RG, Wingreen NS (2008) Chemotaxis in
Escherichia coli : A Molecular Modelfor Robust Precise Adaptation. PLoS Comput Biol 4(1): e1.3. Chen K, et al. (2000) Kinetic Analysis of a Molecular Model of the Budding Yeast Cell Cycle. Mol.Biol. Cell 11: 369-391.4. Chen K, Calzone L, Csikasz-Nagy A, Cross FR, Novak B et al. (2004) Integrative Analysis of CellCycle Control in Budding Yeast. Mol. Biol. Cell 15: 3841-3862.15. Li F, Long T, Lu Y, Ouyang Q, Tang, C (2004) The yeast cell cycle is designed robustly. Proc.Natl. Acad. Sci. U S A 14: 4781-4786.6. Charvin G, Oikonomou C, Siggia ED, Cross FR (2010) Origin of Irreversibility of Cell Cycle Startin Budding Yeast. PLoS Biol 8: e1000284.7. Murray AW, Kirschner MW (1989) Cyclin synthesis drives the early embryonic cell cycle. Nature339: 275-280.8. Ferrell J, Pomerening J, Kim S, Trunnell N, Xiong W, Huang C, Machleder E (2009) Simple,realistic models of complex biological processes: Positive feedback and bistability in a cell fateswitch and a cell cycle oscillator. FEBS Letters 583: 3999-4005.9. Shen-Orr S, Milo R, Mangan S, Alon U (2002) Network motifs in the transcriptional regulationnetwork of
Escherichia coli . Nature Genetics 31: 64-68.10. Alon U (2007) Network motifs: theory and experimental approaches. Nat Rev Genet 8: 450-461.11. Alon U (2007) An introduction to system biology: design principles of biological circuits. Chapmanand Hall.12. Isaacs F, Hasty J, Cantor C, Collins J (2003) Prediction and measurement of an autoregulatorygenetic module. Proc. Natl. Acad. Sci. U S A 100: 7714-7719.13. Xiong W, Ferrell J (2003) A positive-feedback-based bistable ‘memory module’ that governs a cellfate decision. Nature 426: 460-465.14. Thattai M, Lim H, Shraiman B, van Oudenaarden A (2004) Multistability in the lactose utilizationnetwork of Escherichia coli. Nature 427: 737-740.15. Skotheim JM, Di Talia S, Siggia ED, Cross F (2008) Positive feedback of G1 cyclins ensures coherentcell cycle entry. Nature 454: 291-296.16. Ma W, Trusina A, El-Samad H, Lim W, Tang C (2009) Defining Network Topologies that CanAchieve Biochemical Adaptation. Cell 138: 760-773.17. Tyson JJ, Chen K, Novak B (2001) Network dynamics, cell physiology. Nature Reviews MolecularCell Biology 2: 908-916 .18. Csikasz-Nagy A, Battogtokh D, Chen KC, Novak B, Tyson JJ (2006) Analysis of a Generic Modelof Eukaryotic Cell-Cycle Regulation. Biophys J 90: 4361-4379.19. Ferrell J, Tsai T, Yang Q (2011) Modeling the Cell Cycle: Why Do Certain Circuits Oscillate?Cell 144: 874-885.20. Morgan D (2007) The cell cycle: principles of control. United Kingdom: New Science Press. 55 p.21. Novak B, Tyson J, Gyorffy B, Csikasz-Nagy A (2007) Irreversible cell-cycle transitions are due tosystems-level feedback. Nat Cell Biol 9: 724-728.22. Gerard C, Tyson JJ, Novak B (2013) Minimal models for cell-cycle control based on competitiveinhibition and multisite phosphorylations of Cdk substrates. Biophys, J 104: 1367-1379.23. Gerard C, Goldbeter A (2011) A skeleton model for the network of cyclin-dependent kinases drivingthe mammalian cell cycle. Interface Focus 1: 24-35.224. Gerard C, Goldbeter A (2009) Temporal self-organization of the cyclin/Cdk network driving themammalian cell cycle. Proc. Natl. Acad. Sci. U S A 106: 21643-21648.25. Cross F, Siggia E (2005) Mode locking the cell cycle. Phys. Rev. E 72: 021910(6).26. Battogtokh D, Aihara K, Tyson J (2006) Synchronization of Eukaryotic Cells by Periodic Forcing.Phys. Rev. Lett. 96: 148102(4).27. Yang X, Lau KY, Sevim V, Tang C (2013) Design Principles of the Yeast G1/S Switch. PLoS Biol11: e1001673.28. Holt L, Krutchinsky A, Morgan D (2008) Positive feedback sharpens the anaphase switch. Nature454: 353-357.29. Lopez-Aviles S, Kapuy O, Novak B, Uhlmann F (2009) Irreversibility of mitotic exit is the conse-quence of systems-level feedback. Nature 459: 592-595.30. Lu Y, Cross F (2010) Periodic Cyclin-Cdk Activity Entrains an Autonomous Cdc14 Release Os-cillator. Cell 141: 268-279.31. Ma W, Lai L, Ouyang Q, Tang C (2006) Robustness and modular design of the Drosophila segmentpolarity network. Mol Syst Biol. 2: 7032. Lim WA, Lee CM, Tang C (2013) Design principles of regulatory networks: searching for themolecular algorithms of the cell. Mol Cell 49: 202-12.33. Hartwell LH, Weinert TA (1989) Checkpoints: controls that ensure the order of cell cycle events.Science 246: 629-634.34. Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2007) Numberical recipes. CambridgeUniversity Press. 409 p.35. Strogatz Steven H (2001) Nonlinear dynamics and chaos. Westview Press. 99 p.36. Lovrics A, Csikasz-Nagy A, Zsely I, Zador J, Turanyi T, Novak B (2006) Time scale and dimensionanalysis of a budding yeast cell cycle model. BMC Bioinformatics 7: 94-505.37. Murray A (2004) Recycling the Cell Cycle: Cyclins Revisited. Cell 116: 221-234.38. Hartwell L, Hopfield J, Leibler S, Murray A (1999) From molecular to modular cell biology. Nature402: C47-C52.39. Kitano H (2002) Systems Biology: A Brief Overview. Science 295: 1662-1664.40. Sneppen K (2005) Physics in Molecular Biology. Cambridge University Press. 6 p.
Figure Legends x Cln2, Clb5SBF/MBF y z Clb2Sic1, Cdh1Cdc20Mcm1/SFF
G1/S phaseearly M phaselate M phase
Activator node Inhibitor nodeActivation Repression
Figure 1.
The simplified regulatory network in a budding yeast cell-cycle, where x , y , and z representkey regulators of the G1/S, early M, and late M modules, respectively. Different modules are connectedby activation and repression interactions. For details, see the main text and Supporting Information. k a2 k a1 k s k k k k i j j j k a1 T x A B
Figure 2.
The durations of the G1/S and early M phases are controlled the activation rates k a and k a , respectively. (A) Clustering map of the parameter sets in log scale that satisfy the fundamentalconstraint conditions of the yeast cell-cycle process. The insert gives the scale of value. (B) The x waveduration T x of the satisfied parameter sets as a function of k a . The envelope curve corresponds to T x ∝ ( k a − k ca ) − / .4 AB C x y z x y z D Figure 3.
The evolving cell cycle trajectory. (A) and (B) Using the perfect parameters ( k ai = 0 . k ai = 0 . x max , 0, 0)corresponds to the DNA replication event, while the second vertex (0, y max , 0) corresponds to thespindle assembly and separation event.5 Figure 4.
Local dynamic analysis of a perfect cell cycle trajectory ( k ai = 0 . s . The upper panels plot the evolving trajectory. The middlepanels plot the normal radius r n of the local spherical perturbation surface (red curve) and the averagelocal velocity | V | (black curve), showing that the manifolds converge to a small state space when x or y are fully activated. The lower panels plot the real parts of eigenvalues λ and λ of the normal Jacobianmatrix.6 Figure 5.
Local dynamic analysis of an imperfect cell cycle trajectory ( k ai = 0 . r n (red line) and the average local velocity | V | . The lower panels plot the real parts of theeigenvalues λ and λ2