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
], 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
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
] 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
) and M-telophase checkpoint (property 5-7, see Additional file
Timing robustness of the budding yeast cell cycle
The first hazard in the Mangla, et al.
’s model is shown in Figure
. 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
). 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
B). 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
]. However, if Mcm1 changes first, Cdc20 is enabled to change from 0 to 1 (Figure
C). As shown in Figure
D, 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)
]. 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
E). Note that either of Clb2 and Mcm1 can be activated in any order in the Mangla, et al.
’s model. Chen, et al.
] supports this revision, showing that Clb5 makes a stronger reaction to Clb2 than to Mcm1. Other works
] 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
) occurs after Cdc20 transitions to 1 (Figure
A). Since the activation of Swi5 is mostly dependent on the activation of Cdc20
], 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.
], 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
B). In accordance with previous observation
], 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.
C 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
]. 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.
], and the model is widely used due to its acceptable accuracy in explaining a real cell
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
A, 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
D, 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
B 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
D). 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
]. 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
]. 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
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
]. 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.
Average length of state transitions for each phase and their variation
We used the proposed model to study the attractors of the network dynamics by starting from each of the 29
=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
). 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
, 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
] 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
). It is also larger than the one in the Mangla, et al.
’s model (i.e., 2769 or ≈ 90%, Table
). 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.
Attractors and their basin sizes of the proposed model
Attractors and their basin sizes of the Li, et al.’s model
Attractors and their basin sizes of the Mangla, et al.’s model
Number of different state transitions and timing robustness
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
), 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
Number of different state transitions for each phase
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.
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
]. 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
]. 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
]. 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.