A checkpoints capturing timing-robust Boolean model of the budding yeast cell cycle regulatory network

doi:10.1186/1752-0509-6-129

BMC Syst Biol. 2012; 6: 129.

Published online 2012 September 28. doi: 10.1186/1752-0509-6-129

PMCID: PMC3573974

A checkpoints capturing timing-robust Boolean model of the budding yeast cell cycle regulatory network

Changki Hong,¹ Minho Lee,² Dongsup Kim,² Dongsan Kim,² Kwang-Hyun Cho,² and Insik Shin¹

¹Department of Computer Science, KAIST, Daejeon, Korea

²Department of Bio and Brain Engineering, KAIST, Daejeon, Korea

Corresponding author.

Changki Hong: ckhong/at/kaist.ac.kr; Minho Lee: MinhoLee/at/kaist.ac.kr; Dongsup Kim: kds/at/kaist.ac.kr; Dongsan Kim: dongsan/at/kaist.ac.kr; Kwang-Hyun Cho: ckh/at/kaist.ac.kr; Insik Shin: insik.shin/at/cs.kaist.ac.kr

Received October 26, 2011; Accepted August 30, 2012.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Cell cycle process of budding yeast (Saccharomyces cerevisiae) consists of four phases: G1, S, G2 and M. Initiated by stimulation of the G1 phase, cell cycle returns to the G1 stationary phase through a sequence of the S, G2 and M phases. During the cell cycle, a cell verifies whether necessary conditions are satisfied at the end of each phase (i.e., checkpoint) since damages of any phase can cause severe cell cycle defect. The cell cycle can proceed to the next phase properly only if checkpoint conditions are met. Over the last decade, there have been several studies to construct Boolean models that capture checkpoint conditions. However, they mostly focused on robustness to network perturbations, and the timing robustness has not been much addressed. Only recently, some studies suggested extension of such models towards timing-robust models, but they have not considered checkpoint conditions.

Results

To construct a timing-robust Boolean model that preserves checkpoint conditions of the budding yeast cell cycle, we used a model verification technique, ‘model checking’. By utilizing automatic and exhaustive verification of model checking, we found that previous models cannot properly capture essential checkpoint conditions in the presence of timing variations. In particular, such models violate the M phase checkpoint condition so that it allows a division of a budding yeast cell into two before the completion of its full DNA replication and synthesis. In this paper, we present a timing-robust model that preserves all the essential checkpoint conditions properly against timing variations. Our simulation results show that the proposed timing-robust model is more robust even against network perturbations and can better represent the nature of cell cycle than previous models.

Conclusions

To our knowledge this is the first work that rigorously examined the timing robustness of the cell cycle process of budding yeast with respect to checkpoint conditions using Boolean models. The proposed timing-robust model is the complete state-of-the-art model that guarantees no violation in terms of checkpoints known to date.

Keywords: Timing robustness, Yeast cell cycle regulatory network, Model checking, Asynchronous Boolean networks

Background

A cell must undergo the process of duplicating all its components and separating them, more or less evenly, to two daughter cells such that each daughter has the information and dynamics necessary to repeat the process. Such cell cycle dynamics are known in more detail for the budding yeast, Saccharomyces cerevisiae, compared to other eukaryotic organism [1,2]. The cell cycle process of budding yeast consists of four phases: G1, S, G2, and M. Initiated by stimulation of the G1 stationary phase, the cell cycle sequence proceeds (i.e., G1→S→G2→M) and finally returns to the G1 stationary phase. It is important to reach the final phase after completing each phase properly since any mistakes can cause significant defect to the cell cycle process. Hence, a cell verifies whether essential conditions are satisfied at checkpoints in G1, S, G2 and M phases, respectively. Before entering S phase, the cell must be large enough and have undamaged DNA (G1 phase checkpoint). Before entering M phase, DNA synthesis should be completed (S and G2 phase checkpoint). In M phase, the chromosomes need to be properly aligned and the spindles need to be oriented towards the daughter cell (M-metaphase checkpoint), and the cell should be correctly divided into two, before the end of mitosis (M-telophase checkpoint). The cell cycle can proceed to the next phase properly only if the above checkpoint conditions are satisfied.

Among several approaches [3-9] to capture the cell cycle process of budding yeast, one promising approach is to use Boolean network modeling, which is a widely used modeling framework in systems biology [10-18]. Claiming that Boolean models are useful for representing the cell cycle regulatory networks since much of the biology seems to be reflected in the on/off characteristics of the network components, the first Boolean model for the cell cycle regulatory network of budding yeast was introduced by Li, et al.[19]. Even with such a simple representation, they found that there exists a prominent dynamic gene expression trajectory satisfying the checkpoint conditions, and then it leads back to the G1 stationary state. It was also observed that essential conditions are robustly preserved at checkpoints under network perturbations. However, since the model employed synchronous update rules to represent state transitions such that all the nodes update their states synchronously at the same time, it lacks the timing robustness analysis of essential properties.

Timing robustness is the ability of a model to maintain its function in the presence of timing perturbations. Among a few ways to introduce various timing variations, Boolean modeling often uses an asynchronous updating of models. Unlike the synchronous update rule, the asynchronous update rule allows a maximum of one variable to be updated at each time instant, and if multiple variables are enabled to change, one of them is chosen arbitrarily. In this way, variations in reaction rates can be represented depending on the order in which the nodes update their values (i.e., some of the nodes update their values immediately while other nodes take longer). Here it is important to note that this way can generate a large number of distinct state transition trajectories, possibly as many as the number of all different order combinations. And such all different trajectories reflect every possible timing variation under dynamically changing environments. However, such a large number of trajectories to explore make it difficult to perform timing robustness analysis through biological simulation, since such simulation generally involves randomness for trajectory selection.

The main goal of this study is to construct a timing-robust Boolean model that properly preserves checkpoint conditions of the budding yeast cell cycle even in the presence of timing variations. Towards this goal, we used a model verification technique, ‘model checking’ [20], to cope with difficulties of timing robust analysis in conventional simulation based approach. Model checking is a method for formal verification of finite-state systems. A model checker, a software tool of model checking, explores all possible state transitions (i.e., all possible variations in reaction time) of a given model under fully asynchronous update rules. This way, we can then check exhaustively whether essential system properties always hold in the model or not. If the model contains a wrong state transition (i.e., hazard) that leads to a violation of system properties, a model checker automatically detects the hazard and produces a counter-example that can be used to pinpoint the source of the hazard.

In order to utilize model checkers, it needs to specify essential system properties in a form that model checkers can recognize, which is temporal logic [21]. Hence, based on the extensive literature studies [1,19,22-29], this study translates the critical checkpoint conditions of the budding yeast cell cycle into temporal logic formulas. It is worthy to note that such conditions should be properly preserved in order to proceed the cell cycle in the presence of timing variation. For example, an earlier wok by Mangla, et al.[30] aimed to extend the Li, et al.’s model towards a timing-robust model. Though this earlier work was able to construct a model that preserves the G1 stationary state in a timing-robust manner, it fails to capture essential checkpoint conditions completely against variations in reaction rates. It is also turned out that the checkpoint conditions violated in the Mangla, et al.’s model cannot be satisfied in the Li, et al.’s model as well under timing variations. We observed that such inadequate sequences of state transitions can cause a significant failure of the budding yeast cell cycle by violating the M phase checkpoint conditions. In particular, the both models allow a division of a budding yeast cell into two before the completion of its full DNA replication and synthesis. They have the gene expression trajectory that enables the gene Cdc20 to be activated before the activation of the gene Clb2, and such out-of-order sequential dynamics lead a transit to the M phase without completing the G2 phase.

In this study, we present a timing-robust model that properly preserves all the essential checkpoint conditions against timing variations. The proposed model is the complete state-of-the-art model that guarantees no hazard in terms of checkpoints known to date. As a result, our model naturally eliminates the hazards contained in the previous models (i.e., Li, et al.’s and Mangla, et al.’s model). Our simulation results show that the proposed timing-robust model is more robust even against network perturbations and can better represent the nature of cell cycle than previous models. The key to the success of the cell cycle process is to completely capture the checkpoint conditions phase-by-phase even in the presence of variations in reaction rates. To our knowledge this is the first work that rigorously examined the timing robustness of the cell cycle process of budding yeast with respect to checkpoint conditions using Boolean models.

Results and Discussion

In this study, model checkers, software tools of model checking, are used to examine whether or not a specified logical property holds on every possible state of a Boolean model. The inputs to the model checker consist of a Boolean model, described as a set of variables and rules that update their values, an initial state, and a logical property to check. In this study, previously published Boolean models of the budding yeast cell cycle (i.e., the Li, et al.’s model and Mangla, et al.’s model) and the stimulated G1 state are used as input Boolean models and initial state to the model checker, respectively.

We derived the logical properties to check from essential checkpoint conditions of the budding yeast cell cycle. Based on the comprehensive literature studies [1,19,22-29], we found that the key regulators of the S, G2 and M phase checkpoints, Clb2 and Cdc20, are lethal genes. The activation of Clb2 initiates the M phase, and the activation and deactivation of Cdc20 trigger the metaphase to anaphase transition and the exit from mitosis, respectively. Since regulations of Clb2 and Cdc20 are closely related to the checkpoints, any damages of these genes can cause a fatal defect of the cell cycle process. Each checkpoint condition is translated into a group of specific sequences of state transitions that can be derived by ordering pairs of state transition among Clb2, Cdc20 and their interacting genes. For example, the M-metaphase checkpoint conditions can be rendered into two essential sequences of state transitions: Clb2 activation should precede Cdc20 activation; Mcm1 activation should precede Cdc20 activation. All essential properties derived from up-to-date checkpoint conditions are described in Additional file 1 with supporting evidences. Note that no property is derived from the G1 phase checkpoint because Boolean models based on the Li, et al.’s study do not completely include genes related to the checkpoint. Any state transitions violating such essential sequences of state transitions are called hazards. In this study, we used the NuSMV model checker [31] to construct a timing-robust model that properly preserves up-to-date essential checkpoint conditions of the budding yeast cell cycle. With such logical properties as inputs, the model checker detected two hazards in the previously published Boolean models (i.e., Li, et al.’s model and Mangla, et al.’s model) after conducting automatic and exhaustive state space search based on the fully asynchronous update rule. These hazards violate the M-metaphase (property 4, see Additional file 1) and M-telophase checkpoint (property 5-7, see Additional file 1), respectively.

Timing robustness of the budding yeast cell cycle

The first hazard in the Mangla, et al.’s model is shown in Figure Figure1.1. The hazard can lead the model to a biologically undesirable situation, in which it enters the M phase even though the DNA synthesis process is not complete. It violates the property 4 of the M-metaphase checkpoint conditions (see Additional file 1). Such inadequate state transition can occur when Clb5 gets activated in the model. After Clb5 transitions to 1, both Clb2 and Mcm1 are enabled to change from 0 to 1 (Figure (Figure1B).1B). There are two cases, depending on which of Clb2 and Mcm1 changes first. If Clb2 changes first, the cell cycle normally proceeds to the M phase through the G2 phase, guaranteeing the completion of DNA synthesis [23]. However, if Mcm1 changes first, Cdc20 is enabled to change from 0 to 1 (Figure (Figure1C).1C). As shown in Figure Figure1D,1D, the first hazard is detected when Cdc20 transitions to 1 before the activation of Clb2, which means that the division of a cell into two can begin before the completion of DNA replication. Note that the activation of Clb2 is required for the proper cell progression to the M phase and that Cdc20 becomes active after Clb2 phosphorylates APC core proteins (e.g., Cdc16, Cdc23, and Cdc27) [23]. There can be a further timing gap between the activation of Clb2 and Cdc20 in reality because, even if the phosphorylated form of the APC is bound to Cdc20, Cdc20 becomes active only after the chromosomes properly align in the metaphase stage.

Figure 1

Budding yeast model. Nodes in the graph represent molecules. Lines with arrowhead represent the activation, and lines with flat ends indicate the inhibition. Thin arrows represent a low weight, a weight of 1/3, normal arrows indicate a medium weight, (more ...)

The Mangla, et al.’s model can be revised to eliminate the first hazard by adjusting the weight of an edge. When the weight of the edge from Clb5 to Mcm1 decreases to a low level, it ensures that when Clb5 is activated, only Clb2 is able to be subsequently activated, but Mcm1 is not (Figure (Figure1E).1E). Note that either of Clb2 and Mcm1 can be activated in any order in the Mangla, et al.’s model. Chen, et al.[22] supports this revision, showing that Clb5 makes a stronger reaction to Clb2 than to Mcm1. Other works [32,33] also support this revision, revealing that Clb2 and Mcm1 make positive feedbacks; Mcm1 is activated by the low level activation of Clb2, and then the activation of Mcm1 causes Clb2 to be activated at a higher level of concentration in return. Such a biological evidence was reflected in part to the Mangla, et al.’s model such that Clb2 has one of three possible values: 0, representing a negligible concentration of Clb2; 1, representing a low concentration; and 2, representing a high concentration. However, the Mangla, et al.’s model does not capture the positive feedback completely. In the model, Clb2 is able to be activated to at high concentration (i.e., the level of 2) without forming positive feedback. We revised the Mangla, et al.’s model by eliminating the first hazard, and the revised model is able to capture the proper dynamics due to the positive feedback between Clb2 and Mcm1.

The second hazard violating the properties 5 to 7 of the M-telophase checkpoint conditions (see Additional file 1) occurs after Cdc20 transitions to 1 (Figure (Figure2A).2A). Since the activation of Swi5 is mostly dependent on the activation of Cdc20 [22], Swi5 should be activated when Cdc20 transitions to 1 regardless of the activation level of Clb2. In the Mangla, et al.’s model, however, Swi5 is able to change its state only when the activation level of Clb2 decreases to 1 from 2; if the concentration level of Clb2 is 2, Swi5 remains deactivated until the level of Clb2 becomes 1 even though Cdc20 is activated. Such inadequate state transition can cause the model to inevitably delay the exit from mitosis. The complete division of a budding yeast cell related gene, Sic1, is deferred to be activated since the activation of the transcription factor of Sic1 (i.e., Swi5) became delayed.

Figure 2

Budding yeast model. (A) A subset of the Mangla, et al.’s model that highlights the second hazard. Since Clb2 is activated to the highest level, the activation of Cdc20 cannot activate Swi5. (B) The hazard can be eliminated by replacing the weights (more ...)

According to the study by Chen, et al.[22], the activation of Mcm1 causes Swi5 to be transcribed at a low level, and then Swi5 is activated at a higher level after Cdc20 is activated. To capture this understanding, we revise the model so that Swi5 also has one of three possible values: 0, 1, or 2 like Clb2. In addition, the model is revised such that the reactions from Cdc20 and Clb2 make stronger and weaker impacts on Swi5, respectively. This revision ensures that Swi5 is activated to 1 by its transcription factor, Mcm1, and then activated to 2 after Cdc20 is activated regardless of the activation level of Clb2, since reduced inhibition from Clb2 cannot dominate the level of transcription of Swi5 anymore (Figure (Figure2B).2B). In accordance with previous observation [22], the model is additionally revised to have a weaker reaction from Swi5 to Sic1. This revision guarantees the reaction from Swi5 to have impact on the activation of Sic1 only when the concentration of Swi5 increases to the highest level of 2.

Figure Figure2C2C represents the final hazard-free model for the budding yeast cell cycle. We found that the checkpoint conditions violated in the Mangla, et al.’s model cannot be satisfied in the Li, et al.’s model as well. To our knowledge, the proposed model is the first timing-robust model that properly captures up-to-date checkpoint conditions of the budding yeast cell cycle in the presence of variations in reaction rates.

Temporal evolution of gene states for budding yeast cell cycle models

Mathematical modeling based on biochemical rate equations, provides a rigorous and reliable tool for unraveling the complexities of molecular regulatory networks. However, this approach is only suited for small and well-characterized systems with known kinetic parameters since there is a lack of detailed knowledge of quantitative reaction kinetics for most of the reactions in a cell [1,34]. Fortunately, the cell cycle regulatory system of budding yeast is most fully worked out, so its control system is revealed in exquisite details. The mathematical model of the budding yeast cell cycle has been published by Chen, et al.[22], and the model is widely used due to its acceptable accuracy in explaining a real cell [35,36].

To investigate how closely Boolean models of the budding yeast cell cycle represent nature, we compared the temporal evolutions of gene states in Boolean models to that in the mathematical model. Note that the asynchronous update rule allows a maximum of one variable to be updated at each time instant, and if multiple variables are enabled to change, one of them is chosen in an arbitrary fashion. Thus, too many different state transitions can exist under the asynchronous update rule. When applying the synchronous update rule to Boolean models for simple comparison, we observed that the temporal evolution of gene states in the proposed model maintains a similar structure to those in the other Boolean models overall. However, we also found that the extension from the Mangla, et al.’s model leads some genes in the proposed model to evolve analogously to the dynamics of the corresponding genes in the mathematical model in some period, and such genes are closely related to the hazards which are eliminated during the extension.

As shown in Figure Figure3A,3A, the mathematical model clearly shows a positive feedback between Clb2 and Mcm1: Clb2 is transcribed at a low concentration in the interval from 30 to 45, which is sufficient to activate Mcm1, and then Clb2 is transcribed at a higher level at the time of 50 via the activation of Mcm1. The temporal evolution of the proposed model also forms positive feedback as shown in Figure Figure3D,3D, but the other two models cannot result in such a relationship since both Clb2 and Mcm1 are activated simultaneously at the time of 5. In the application of the asynchronous update rule, those previously published Boolean models (Figures (Figures3B3B and 3C) can generate the hazard since Cdc20 can be activated before the transcription of Clb2 if Mcm1 is activated first under the condition that both Clb2 and Mcm1 are enabled to change.

Figure 3

Temporal evolution of variables related to the first hazard in the budding yeast cell cycle model. X-axis and y-axis represent the time and the level of activation concentration in simulation, respectively. (A) The low level activation of Clb2 through (more ...)

Moreover, the temporal evolution of Swi5 in the proposed model is closer to the dynamics of the gene in reality (Figure (Figure4D).4D). Consistent with the mathematical model, Swi5 in the proposed model is first transcribed at a low concentration (i.e., a value of 1) at the time of 7 by the activation of Mcm1. It is then activated to a high concentration (i.e., a value of 2) at the next time step after Cdc20 is activated. However, in the other Boolean models, it is shown that Swi5 is still in the inactive state even after Mcm1 becomes activated. Specifically, Swi5 is deactivated until Clb2 is degraded to the low level of activation regardless of the activation of Cdc20 in the Mangla, et al.’s model.

Figure 4

Temporal evolution of variables related to the second hazard in the budding yeast cell cycle model. X-axis and y-axis represent the time and the level of activation concentration in simulation, respectively. (A) Swi5 is activated to a low level by Mcm1, (more ...)

From this simulation, it appears that the proposed model follows similar dynamics as those in the other Boolean models, even if the proposed model is extended from the others. In addition, in the proposed model, the temporal evolution of genes related to the extension shows better consistency with the one in the mathematical model of the budding yeast cell cycle. Therefore, the proposed model which is constructed by eliminating hazards in the previous models better reflects the nature of cell cycle than the previously published Boolean models.

Relative durations of each cell cycle phase

Many efforts have been made to discover checkpoints and key transcription factors responsible for phase transitions in cell cycles [1,19,22-29]. Therefore, it is useful for cell cycle models to reflect properties on the duration of each phase in a cell cycle.

Recent studies have assumed that every edge in the yeast cell cycle regulatory network proceeds with the same speed since transcriptions normally happen on similar time scales [37,38]. Following this assumption, we compared the average length of state transitions in each phase in the models, referring to the relative duration of each phase in nature. The S phase begins when the state of the model is the same as the stationary G1 phase except that Cln3 is 1. Table Table11 shows that in the Li, et al.’s model, the average length of the S phase is 9.30 and the average length of the G2/M is 7.56. In the Mangla, et al.’s model, the average length of the G2/M becomes longer, which is 10.89. These results are inconsistent with the experimental data, in which it was observed that the relative duration of the G2/M phase is twice as long as that of the S phase [39,40]. However, in the proposed model, the G2/M phase is almost twice as long on average than the S phase. In our understanding, this is because the S phase becomes shorter with our revision of the Mangla, et al.’s model with regard to positive feedback between Clb2 and Mcm1, and the G2/M phase becomes longer with our extension of Swi5’s activation level, which produced additional state transitions to reach the end of the M phase.

Table 1

Average length of state transitions for each phase and their variation

Attractor analysis

We used the proposed model to study the attractors of the network dynamics by starting from each of the 2⁹×3²=4,608 states in the 11-node proposed model with two 3-valued nodes. We found that all of the initial states eventually flow into one of the nine stationary states, also called attractors (Table (Table2).2). The basin size of an attractor is the number of initial states which reach the attractor after a finite number of time steps. As seen in Table Table2,2, there is one big fixed point which attracts 4,323 or ≈ 94% protein states from 4,608 initial states. This is consistent with a previous study [19] which reveals that the model has one big attractor, and the dominant attractor is the biological G1 stationary state. Interestingly, the basin size of the big attractor in the proposed model is much larger than the one presented in the Li, et al.’s model (i.e., 1,764 or ≈ 86%, Table Table3).3). It is also larger than the one in the Mangla, et al.’s model (i.e., 2769 or ≈ 90%, Table Table4).4). It is obvious that the proposed model is more stable than the other models of budding yeast because the big attractor represents a cell’s stationary state, and the basin size is the largest in the proposed model. Under normal conditions, a cell will be sitting at the state that the biggest attractor represents, waiting for another round of divisions.

Table 2

Attractors and their basin sizes of the proposed model

Table 3

Attractors and their basin sizes of the Li, et al.’s model

Table 4

Attractors and their basin sizes of the Mangla, et al.’s model

Number of different state transitions and timing robustness

Table Table55 shows the number of different state transitions for each cell cycle phase, and the major difference among the three models is found in the G2 phase. To investigate the relationships between the number of different state transitions and the timing robustness of the model, we introduced random mutations into the models. To perturb only an interested cell cycle phase, we applied mutations to edges related to the corresponding phase by deleting, adding, or re-weighting them. In this way, 200 mutant networks were generated, and then checked whether or not the dynamics of networks preserve the G1 stationary state as their global attractor under the fully asynchronous update rule using the NuSMV model checker. In the phase transition from the G2 phase to the M phase, the Mangla, et al.’s model was found to be vulnerable to mutations (Table (Table6),6), and the fraction of timing-robust mutants holding the cell cycle property was much smaller than in the proposed model. It was even smaller than those from the other phases in the same model. This is interpreted to indicate that incomplete construction of the positive feedback between Clb2 and Mcm1 in the Mangla, et al.’s model causes the G2 phase to be fragile to perturbations. In addition, the extension of the budding yeast cell cycle model by elimination of such hazards increases the number of viable state transitions to correctly arrive at the M phase from the G2 phase, which allows the proposed model to have about 10% greater timing robustness. The Li, et al.’s model includes hazards which drives the model to different attractors rather than the G1 stationary state. Thus, mutants of the model are also unlikely to hold the property, as shown in Table Table66.

Table 5

Number of different state transitions for each phase

Table 6

Fraction of timing-robust mutants for each phase

We conjecture that the network robustness of a model is closely related to the number of different biological state transitions which satisfy properties of the model. This is because even if some state transitions are altered by mutations, the remaining transitions are likely to hold biological properties of the model, and our simulation results demonstrate this.

Limitations

The proposed model is the timing-robust model that properly preserves all the essential checkpoint conditions against timing variations. Although the proposed model is the complete state-of-the-art model that guarantees no hazard in terms of checkpoints known to date, it still has a gap with quantitative models. Our model is not yet directly applicable to explaining and predicting the quantitative outcome of biological experiments of the budding yeast cell cycle. Quantitative models can potentially describe molecular interactions with high precision and in quantitative terms that correspond to realistic laboratory measurements. However, Boolean models can still be used by a subset of researchers because of easy understanding of the dynamics of the budding yeast cell cycle. We expect that the proposed model can be used to help them by providing a more stable and timing-robust Boolean model.

The main obstacle in application of model checking in practice is the state space explosion problem [20]. Since model checkers should examine all the possible model states, the number of states can grow exponentially in the number of program variables. When verifying a large-scale biological system, it is intractable to explore the entire state spaces because they exceed computational limits (i.e., time and memory). This problem is known as state space explosion. The last 30 years have seen various techniques for resolving the state space explosion problem [41]. In particular, several techniques have been introduced to decrease the number of states to be explored and the memory requirements needed for storing explored states, leading to substantial reduction in state space [42]. We further expect that biological experts can propose biology-specific insight and knowledge to state space reduction techniques, such as biological abstraction that allows a set of biologically equivalent states to be considered as a single symbolic state, resulting in significant state space reduction.

Conclusions

Timing robustness analysis is one of the most important and challenging problems in systems biology [43-45]. It is critical for biological systems to maintain its essential dynamics robustly under various reaction delays in dynamically changing environments [46-49]. Pointing out that time delays are common and substantial in gene regulatory networks, Chen, et al.[50] proposed a method to design robust gene regulatory networks under biochemical timing variations and molecular noises. Lopez-Aviles, et al.[51] revealed that time delays have a significant influence on fundamental dynamics by showing that the unidirectionality of eukaryotic cell cycle transitions (i.e., G1→S→G2→M→G1) can hold only after a certain amount of time delays of cyclin degradation.

Boolean network modeling, which is now a widely used modeling framework in systems biology, requires timing robustness analysis since reaction kinetic parameters inevitably vary over a certain range. However, most researchers paid little attention to timing robust analysis so far and assumed that Boolean models are updated in a synchronous manner, neglecting timing variations [30]. Considering the importance of robustness of biological systems over timing variations, timing robustness analysis needs to be carefully considered in most Boolean network models.

A number of experimental studies on budding yeast, Saccharomyces cerevisiae, illuminated its cell cycle dynamics in greater detail than any other eukaryotic organisms [1,2]. Such rich experimental results have fostered a growing attention to modeling of the budding yeast cell cycle network as it is recognized as a good benchmark example for studying the fundamental design principles behind the well-orchestrated behavior [3-9]. In particular, Li, et al.[19] presented a yeast cell cycle network model that has brought a substantial influence on other studies on the yeast cell cycle. It provided a basis for constructing more robust models under stochastic environmental fluctuation [52-55]. Other studies extended the Li, et al.’s model by incorporating activation/deactivation delays with the auxiliary nodes [24,56].

Building upon the Li, et al.’s model, Mangla, et al.[30] was able to construct a model that preserves the G1 stationary state in a timing-robust manner, it fails to capture essential checkpoint conditions completely against variations in reaction rates. Hence, we proposed here a timing-robust Boolean model that properly preserves up-to-date checkpoint conditions of the budding yeast cell cycle. Our model provides a basis for other subsequent further studies on budding yeast cell cycle analysis with Boolean network models that would result in biologically more robust results.

Boolean modeling and analysis of complex biological networks aim to provide a system-level understanding on complex biological phenomena [10-18]. In practice, however, constructing Boolean network models from biological data typically requires a significant amount of manual efforts through repetitive modeling and checking processes; models are commonly revised iteratively until they conform with targeted essential behavior of biological systems [57-60]. Model checking could facilitate such repetitive designing processes by performing automatic checking of biologically plausible models and suggesting new testable predictions upon every model-based simulation failure [11,20,61-65]. For instance, Fisher, et al.[61] inferred new regulation of inductive and lateral signaling crosstalk of C. elegans according to the testable predictions suggested by a model checker and then confirmed the newly inferred regulation with biological experiments.

In systems biology, mathematical models are becoming too complicated to be validated by examining some essential dynamics in an ad hoc way. It becomes even more difficult to check manually whether a combination of dynamics (e.g., ordered dynamics) are met simultaneously in complex biological networks [66]. On the other hand, model checking can validate a set of essential functionalities and their combinations in an automatic manner. So, model checking can be beneficial to modeling and analysis of biological systems, particularly for developmental systems [67] and cell-fate decision systems [68] where it is crucial to maintain ordered dynamics over environmental variations. In this respect, our model checking-based approach can be useful and provide a systematic framework for robustness analysis.

Methods

NuSMV

Model checking [20] is a technique for the verification of correctness properties of finite-state systems. One benefit of this technique is that systems can be automatically verified by use of a tool, called a model checker. The system is translated into a model described in the input language of the model checker. The model checker explores every state space of the given model to check whether system properties hold in the model or not. If the model violates any of the system specifications, the model checker produces an execution trace (i.e., counter-example) showing why the specification turns out false in the model. Counter-examples usually serve as good hints when a model is revised.

We used a symbolic model checker, NuSMV [31], as in the study by Mangla, et al.[30]. NuSMV can easily switch its modes in exploring state spaces between synchronous and asynchronous manners by introducing a single control variable. NuSMV, therefore, is useful to test whether a budding yeast cell cycle model satisfies its own cell cycle properties both in synchronous and asynchronous manners. We chose a Boolean decision diagram (BDD)-based implementation for the exploration of state spaces among various implementation candidates. Traditional representation of Boolean functions generates redundant state spaces explicitly when exploring state spaces. On the other hand, the BDD data structure can keep much smaller state spaces in the exploration of state spaces than the traditional ones do. Thus, BDD is beneficial in dealing with the state space explosion problem, which is one of the main challenges of model checking techniques for the verification of large-scale systems.

A NuSMV model consists of one or more modules. Each module can declare variables and their update rules. Variables can be declared to have a range of discrete values. The rules specify how to initiate variables and update them at every time step from their current values. Update rules can be non-deterministic since they can result in different values of a variable under the same condition. NuSMV checks whether a given property holds over all different possibilities. In describing our budding yeast cell cycle model in the NuSMV input language, we assign a single variable to each node of the model and specify the update rules of every variable. By default, variables in the NuSMV model are updated in a synchronous manner. There is a global clock and all modules execute in parallel every time the global clock ticks. To apply the asynchronous update rule to the NuSMV model, we define an additional control variable. At each time step, the control variable indicates which variable to update in a non-deterministic manner such that only the update rules of the chosen variable can execute in the model. In addition, we add a FAIRNESS property to every variable to prevent the control variable from repeatedly choosing the same variables for update, which would cause other variables not to be chosen for a long time.

Programs in the language can be annotated by properties expressed in temporal logic, that is, computation tree logic (CTL) [21]. Besides the properties describing the sequences of transitions between states, some properties such as “eventually” or “never” can also be specified with special temporal operators in CTL. For example, cell cycle properties (e.g., a cell starts from the G1 phase, and eventually returns to the G1 stationary phase after a division; Cdc20 activation should never begin before Clb2 activation) can be encoded by CTL formulas.

Boolean model construction

In this paper, we presented Boolean networks of the budding yeast cell cycle. Boolean models include nodes and edges for different components and interactions of the system, respectively. Each node in the Boolean model has one of two values: 1 for ON (active) and 0 for OFF (inactive). A state S of each node i at any time instant t (denoted as S_i(t)) is determined according to a Boolean function (rule) and the states of its input nodes at the previous time instant t−1. Most Boolean functions are threshold based [18]. When the value of a node is updated, it is assigned one if the weighted sum of positive and negative inputs exceeds a pre-defined threshold for that node. For example, the Li, et al.’s model describes its Boolean rule as follows [19]:

where w_ijrepresents the weight of an incoming edge to a node i from a node j, and the threshold (denoted by θ_i) is set to zero.

A synchronous Boolean model is one of the simplest implementations for the application of such Boolean update rules to nodes. In such a model, a Boolean update rule is applied to all the nodes simultaneously at each time instant. Synchronous models are deterministic since nodes are assumed to work in the same time scale, resulting in convergence to the same state from the same initial condition after the same number of time steps. As the result of applying the Boolean function described above synchronously, the first synchronous Boolean model for the budding yeast cell cycle, the Li, et al.’s model, converged into seven attractors from 2¹¹ initial states in the 11-node network model.

Although synchronous Boolean models have been widely used due to their simple nature and ease of implementation, they lack consideration of a variety of time scales in biological systems. To deal with this drawback, asynchronous models are suggested in which a maximum of one node is chosen to be updated at each time instant. Since it is usually unknown exactly how long specific biological processes take, most asynchronous algorithms are non-deterministic in a way that a single node is randomly chosen at each time unit. The Mangla, et al.’s model is constructed by the use of the NuSMV model checker such that it can preserve its desirable attractor, the G1 stationary state in every possible combination of node updating order. The model is designed to be independent of whatever time scales biological reactions are subject to; thus the model is considered timing-robust model. To construct a timing-robust budding yeast cell cycle model dealing with the hazards of the Li, et al.’s model, the Mangla, et al.’s model allows Clb2 to have one of three possible values (i.e., 0, 1 or 2) and edges assigned to one of three magnitudes of weights (i.e., equation M2

, ±1 or ±3). The Boolean rule is also defined, extending that of the Li, et al.’s model, as follows:

In the model, the thresholds for the value of 1 (denoted by θ_i,1) and 2 (denoted by θ_i,2) for a node i are set to 0.5 and 1.5, respectively.

Basically, the model proposed in this study extended a set of available values and weights for each node and edge, respectively. The proposed model also follows the same Boolean function as the Mangla, et al.’s model follows, but it refined the Mangla, et al.’s model to hold essential checkpoint conditions as well. It finally became the timing-robust model that captures up-to-date checkpoint conditions of the budding yeast cell cycle.

Competing interests

The authors declare that they have no competing interests.

Author’s contributions

CH and IS jointly conceived the study, and wrote the manuscript. ML and DK (Dongsan Kim) interpreted simulation results and helped draft the manuscript. DK (Dongsup Kim) and KC critically reviewed the manuscript. All authors read and approved the final manuscript.

Supplementary Material

Additional file 1

Essential ordered properties derived from checkpoints. The PDF file contains a list of all the essential ordered properties derived from the up-to-date checkpoint conditions.

Click here for file^{(99K, pdf)}

Acknowledgements

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MEST) (2009-0086964).

References

Chen KC, Csikasz-Nagy A, Gyorffy B, Val J, Novak B, Tyson JJ. Kinetic Analysis of a Molecular Model of the Budding Yeast Cell Cycle. Mol Biol Cell. 2000;11:369–391. [PMC free article] [PubMed]
Charvin G, Cross FR, Siggia ED. Forced periodic expression of G1 cyclins phase-locks the budding yeast cell cycle. Proceedings of the National Academy of Sciences. 2009;106(16):6632–6637. [PubMed]
Shmulevich I, Kauffman SA, Aldana M. Eukaryotic cells are dynamically ordered or critical but not chaotic. Proceedings of the National Academy of Sciences of the United States of America. 2005;102(38):13439–13444. [PubMed]
Kauffman SA. Metabolic stabiligy and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969;22:437–467. [PubMed]
Giot L, Bader JS, Brouwer C, Chaudhuri A, Kuang B, Li Y, Hao YL, Ooi CE, Godwin B, Vitols E, Vijayadamodar G, Pochart P, Machineni H, Welsh M, Kong Y, Zerhusen B, Malcolm R, Varrone Z, Collis A, Minto M, Burgess S, McDaniel L, Stimpson E, Spriggs F, Williams J, Neurath K, Ioime N, Agee M, Voss E, et al. Furtak K. A Protein Interaction Map of Drosophila melanogaster. Science. 2003;302(5651):1727–1736. [PubMed]
Lee TI, Rinaldi NJ, Robert F, Odom DT, Bar-Joseph Z, Gerber GK, Hannett NM, Harbison CT, Thompson CM, Simon I, Zeitlinger J, Jennings EG, Murray HL, Gordon DB, Ren B, Wyrick JJ, Tagne JB, Volkert TL, Fraenkel E, Gifford DK, Young RA. Transcriptional Regulatory Networks in Saccharomyces cerevisiae. Science. 2002;298(5594):799–804. [PubMed]
Shen L, Chepelev I, Liu J, Wang W. Prediction of quantitative phenotypes based on genetic networks: a case study in yeast sporulation. BMC Syst Biol. 2010;4:128. [PMC free article] [PubMed]
Ding S, Wang W. Recipes and mechanisms of cellular reprogramming: a case study on budding yeast Saccharomyces cerevisiae. BMC Syst Biol. 2011;5:50. [PMC free article] [PubMed]
Eriksson O, Andersson T, Zhou Y, Tegner J. Decoding complex biological networks - tracing essential and modulatory parameters in complex and simplified models of the cell cycle. BMC Syst Biol. 2011;5:123. [PMC free article] [PubMed]
Helikar T, Konvalina J, Heidel J, Rogers JA. Emergent decision-making in biological signal transduction networks. Proceedings of the National Academy of Sciences. 2008;105(6):1913–1918. [PubMed]
Faure A, Naldi A, Chaouiya C, Thieffry D. Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle. Bioinformatics. 2006;22(14):e124–e131. [PubMed]
Davidich MI, Bornholdt S. Boolean Network Model Predicts Cell Cycle Sequence of Fission Yeast. PLoS ONE. 2008;3(2):e1672. [PMC free article] [PubMed]
Mendoza L, Thieffry D, Alvarez-Buylla ER. Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis. Bioinformatics. 1999;15(7):593–606. [PubMed]
Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J Theor Biol. 2003;223:1–18. [PubMed]
Klamt S, Saez-Rodriguez J, Lindquist J, Simeoni L, Gilles E. A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinf. 2006;7:56. [PMC free article] [PubMed]
Chaves M, Albert R, Sontagm ED. Robustness and fragility of Boolean models for genetic regulatory networks. J Theor Biol. 2005;235(3):431–449. [PubMed]
Saez-Rodriguez J, Simeoni L, Lindquist JA, Hemenway R, Bommhardt U, Arndt B, Haus UU, Weismantel R, Gilles ED, Klamt S, Schraven B. A Logical Model Provides Insights into T Cell Receptor Signaling. PLoS Comput Biol. 2007;3(8):e163. [PubMed]
Bornholdt S. Boolean network models of cellular regulation: prospects and limitations. J R Soc Interface. 2008;5(Suppl 1):S85–S94. [PMC free article] [PubMed]
Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell cycle network is robustly designed. Proceedings of the National Academy of Sciences of the United States of America. 2004;101:4781–4786. [PubMed]
Clarke EM, Grumberg OJ, Peled DA. Model Checking. Cambridge: The MIT Press; 1999.
Heljanko K. Model Checking the Branching Time Temporal Logic CTL. Research Report A45, Helsinki University of Technology, Digital Systems Laboratory, Espoo, Finland. 1997.
Chen KC, Calzone L, Csikasz-Nagy A, Cross FR, Novak B, Tyson JJ. Integrative analysis of cell cycle control in budding yeast. Mol Biol Cell. 2004;15:3841–3862. [PMC free article] [PubMed]
Calzone L. Temporal organization of the budding yeast cell cycle: General principles and detailed simulations. PhD thesis, Virginia Polytechnic Institute and State University, Department of Biology 2003.
Irons D. Logical analysis of the budding yeast cell cycle. J Theor Biol. 2009;257(4):543–559. [PubMed]
Simon I, Barnett J, Hannett N, Harbison CT, Rinaldi NJ, Volkert TL, Wyrick JJ, Zeitlinger J, Gifford DK, Jaakkola TS, Young RA. Serial Regulation of Transcriptional Regulators in the Yeast Cell Cycle. Cell. 2001;106(6):697–708. [PubMed]
Zhu G, Spellman PT, Volpe T, Brown PO, Botstein D, Davis TN, Futcher B. Two yeast forkhead genes regulate the cell cycle and pseudohyphal growth. Nature. 2000;406:90–94. [PubMed]
Kumar R, Reynolds DM, Shevchenko A, Shevchenko A, Goldstone SD, Dalton S. Forkhead transcription factors, Fkh1p and Fkh2p, collaborate with Mcm1p to control transcription required for M-phase. Curr Biol. 2000;10(15):896–906. [PubMed]
MacKay VL, Mai B, Waters L, Breeden LL. Early Cell Cycle Box-Mediated Transcription of CLN3 and SWI4 Contributes to the Proper Timing of the G1-to-S Transition in Budding Yeast. Mol Cell Biol. 2001;21(13):4140–4148. [PMC free article] [PubMed]
Lydall D, Ammerer G, Nasmyth K. A new role for MCM1 in yeast: cell cycle regulation of Swi5 transcription. Genes and Dev. 1991;5:2405–2419. [PubMed]
Mangla K, Dill DL, Horowitz MA. Timing Robustness in the Budding and Fission Yeast Cell Cycles. PLoS ONE. 2010;5:e8906. [PMC free article] [PubMed]
Cimatti A, Clarke E, Giunchiglia E, Giunchiglia F, Pistore M, Roveri M, Sebastiani R, Tacchella A. Proc. International Conference on Computer-Aided Verification (CAV 2002) Volume 2404 of. LNCS, Copenhagen, Denmark: Springer; 2002. NuSMV Version 2: An OpenSource Tool for Symbolic Model Checking.
Amon A, Tyers M, Futcher B, Nasmyth K. Mechanisms that help the yeast cell cycle clock tick: G2 cyclins transcriptionally activate G2 cyclins and repress G1 cyclins. Cell. 1993;74:993–1007. [PubMed]
Maher M, Cong F, Kindelberger D, Nasmyth K, Dalton S. Cell cycle-regulated transcription of the CLB2 gene is dependent on Mcm1 and a ternary complex factor. Mol Cell Biol. 1995;15:3129–3137. [PMC free article] [PubMed]
von Dassow G, Meir E, Munro EM, Odell GM. The segment polarity network is a robust develpmental module. Nature. 2000;406:188–192. [PubMed]
Sobie EA. Computational Modeling of the Cell Cycle. Science Signaling. 2011;4(192):tr11. [PubMed]
Vinod PK, Freire P, Rattani A, Ciliberto A, Uhlmann F, Novak B. Computational modelling of mitotic exit in budding yeast: the role of separase and Cdc14 endocycles. J R Soc Interface. 2011;8(61):1128–1141. [PMC free article] [PubMed]
Larson DR. What do expression dynamics tell us about the mechanism of transcription? Curr Opin Genet Dev. 2011;21(5):591–599. [PMC free article] [PubMed]
Peil K, Varv S, Looke M, Kristjuhan K, Kristjuhan A. Uniform Distribution of Elongating RNA Polymerase II Complexes in Transcribed Gene Locus. J Biol Chem. 2011;286(27):23817–23822. [PMC free article] [PubMed]
Rustici G, Mata J, Kivinen K, Lio P, Penkett CJ, Burns G, Hayles J, Brazma A, Nurse P, Bahler J. Periodic gene expression program of the fission yeast cell cycle. Nat Genet. 2004;36:809–817. [PubMed]
Orlando DA, Lin CY, Bernard A, Wang JY, Socolar JES, Iversen ES, Hartemink AJ, Haase SB. Global control of cell cycle transcription by coupled CDK and network oscillators. Nature. 2008;453:944–947. [PMC free article] [PubMed]
Grumberg O, Veith H (Eds) 25 Years of Model Checking - History, Achievements, Perspectives, Volume 5000 of Lecture Notes in Computer Science, Springer; 2008.
Pelanek R. In: Formal Methods for Industrial Critical Systems, Volume 5596 of Lecture Notes in Computer Science. Cofer D, editor. Fantechi A: Springer Berlin / Heidelberg; 2009. Fighting State Space Explosion: Review and Evaluation; pp. 37–52.
Rao CV, Wolf DM, Arkin AP. Control, exploitation and tolerance of intracellular noise. Nature. 2002;420:231–237. [PubMed]
McAdams H, Arkin A. Stochastic mechanisms in gene expression. Proceedings of the National Academy of Sciences. 1997;94(3):814–819. [PubMed]
Hasty J, McMillen D, Isaacs F, Collins JJ. Computational studies of gene regulatory networks: in numero molecular biology. Nat Rev Genet. 2001;2:268–279. [PubMed]
Kitano H. Systems Biology: A Brief Overview. Science. 2002;295(5560):1662–1664. [PubMed]
Csete ME, Doyle JC. Reverse Engineering of Biological Complexity. Science. 2002;295(5560):1664–1669. [PubMed]
Kitano H. Biological robustness. Nat Rev Genet. 2004;5:826–837. [PubMed]
Kitano H. Towards a theory of biological robustness. Mol Syst Biol. 2007;3:137. [PMC free article] [PubMed]
Chen BS, Chang YT. A systematic molecular circuit design method for gene networks under biochemical time delays and molecular noises. BMC Syst Biol. 2008;2:103. [PMC free article] [PubMed]
Lopez-Aviles S, Kapuy O, Novak B, Uhlmann F. Irreversibility of mitotic exit is the consequence of systems-level feedback. Nature. 2009;459:592–595. [PMC free article] [PubMed]
Zhang Y, Qian M, Ouyang Q, Deng M, Li F, Tang C. Stochastic model of yeast cell-cycle network. Physica D: Nonlinear Phenomena. 2006;219:35–39.
Braunewell S, Bornholdt S. Superstability of the yeast cell-cycle dynamics: Ensuring causality in the presence of biochemical stochasticity. J Theor Biol. 2007;245(4):638–643. [PubMed]
Ge H, Qian H, Qian M. Synchronized dynamics and non-equilibrium steady states in a stochastic yeast cell cycle network. Math Biosci. 2008;211:132–152. [PubMed]
Okabe Y, Sasai M. Stable Stochastic Dynamics in Yeast Cell Cycle. Biophys J. 2007;93(10):3451–3459. [PubMed]
Faure A, Chaouiya C, Ciliberto A, Thieffry D. Logical modelling and analysis of the budding yeast cell cycle. BMC Bioinf. 2007;8(Suppl 8):P1.
Akutsu T, Miyano S, Kuhara S. Proceedings of the Pacific Symposium on Biocomputing, Big Island, Hawaii: World Scientific Press; 1999. Identification Of Genetic Networks From A Small Number Of Gene Expression Patterns Under The Boolean Network Model; pp. 17–28.
Ideker TE, Thorsson V, Karp RM. Proceedings of the Pacific Symposium on Biocomputing, Oahu, Hawaii: World Scientific Press; 2000. Discovery of Regulatory Interactions Through Perturbation: Inference and Experimental Design; pp. 305–316.
Maki Y, Tominaga D, Okamoto M, Eguchi Y. , Watanabe S. Proceedings of the Pacific Symposium on Biocomputing, Big Island, Hawaii: World Scientific Press; 2001. Development Of A System For The Inference Of Large Scale Genetic Networks; pp. 446–458.
LÃhdesmÃki H, Shmulevich I, Yli-Harja O. On learning gene regulatory networks under the Boolean network model. Machine Learning. 2003. pp. 147–167.
Schaub M, Henzinger T, Fisher J. Qualitative networks: a symbolic approach to analyze biological signaling networks. BMC Syst Biol. 2007;1:4. [PMC free article] [PubMed]
Fisher J, Henzinger T, Mateescu M, Piterman N. Formal Methods in Systems Biology, Volume 5054 of Lecture Notes in Computer Science. Edited by Fisher J: Springer Berlin / Heidelberg; 2008. Bounded Asynchrony: Concurrency for Modeling Cell-Cell Interactions; pp. 17–32.
Garg A, Di Cara A, Xenarios I, Mendoza L, De Micheli G. Synchronous versus asynchronous modeling of gene regulatory networks. Bioinformatics. 2008;24(17):1917–1925. [PubMed]
Fages F, Soliman S, Chabrier-rivier N. Modelling and querying interaction networks in the biochemical abstract machine biocham. J Biol Phys Chem. 2004;4:64–73.
Koh CH, Nagasaki M, Saito A, Li C, Wong L, Miyano S. MIRACH: efficient model checker for quantitative biological pathway models. Bioinformatics. 2011;27(5):734–735. [PubMed]
Li C, Nagasaki M, Ueno K, Miyano S. Simulation-based model checking approach to cell fate specification during Caenorhabditis elegans vulval development by hybrid functional Petri net with extension. BMC Syst Biol. 2009;3:42. [PMC free article] [PubMed]
Nakajima A, Isshiki T, Kaneko K, Ishihara S. Robustness under Functional Constraint: The Genetic Network for Temporal Expression in Drosophila Neurogenesis. PLoS Comput Biol. 2010;6(4):e1000760. [PMC free article] [PubMed]
Lee T, Yao G, Nevins J, You L. Sensing and Integration of Erk and PI3K Signals by Myc. PLoS Comput Biol. 2008;4(2):e1000013. [PMC free article] [PubMed]

Articles from BMC Systems Biology are provided here courtesy of
BioMed Central