

ORIGINAL ARTICLE 

Year : 2020  Volume
: 10
 Issue : 2  Page : 94104 

Extended robust boolean network of budding yeast cell cycle
Sajad Shafiekhani^{1}, Mojtaba Shafiekhani^{2}, Sara Rahbar^{3}, Amir Homayoun Jafari^{3}
^{1} Department of Biomedical Engineering, School of Medicine; Research Center for Biomedical Technologies and Robotics; Students' Scientific Research Center, Tehran University of Medical Sciences, Tehran, Iran ^{2} Department of Biomedical Engineering, Amirkabir University of Technology, Tehran, Iran ^{3} Department of Biomedical Engineering, School of Medicine; Research Center for Biomedical Technologies and Robotics, Tehran University of Medical Sciences, Tehran, Iran
Date of Submission  14Aug2019 
Date of Decision  22Oct2019 
Date of Acceptance  20Oct2019 
Date of Web Publication  25Apr2020 
Correspondence Address: Amir Homayoun Jafari Department of Medical Physics and Biomedical Engineering, School of Medicine, Tehran University of Medical Sciences, Tehran Iran
Source of Support: None, Conflict of Interest: None
DOI: 10.4103/jmss.JMSS_40_19
Background: How to explore the dynamics of transition probabilities between phases of budding yeast cell cycle (BYCC) network based on the dynamics of protein activities that control this network? How to identify the robust structure of protein interactions of BYCC Boolean network (BN)? Budding yeast allows scientists to put experiments into effect in order to discover the intracellular cell cycle regulating structures which are well simulated by mathematical modeling. Methods: We extended an available deterministic BN of proteins responsible for the cell cycle to a Markov chain model containing apoptosis besides G1, S, G2, M, and stationary G1. Using genetic algorithm (GA), we estimated the kinetic parameters of the extended BN model so that the subsequent transition probabilities derived using Markov chain model of cell states as normal cell cycle becomes the maximum while the structure of chemical interactions of extended BN of cell cycle becomes more stable. Results: Using kinetic parameters optimized by GA, the probability of the subsequent transitions between cell cycle phases is maximized. The relative basin size of stationary G1 increased from 86% to 96.48% while the number of attractors decreased from 7 in the original model to 5 in the extended one. Hence, an increase in the robustness of the system has been achieved. Conclusion: The structure of interacting proteins in cell cycle network affects its robustness and probabilities of transitions between different cell cycle phases. Markov chain and BN are good approaches to study the stability and dynamics of the cell cycle network.
Keywords: Boolean network, budding yeast cell cycle, genetic algorithm, Markov chain model
How to cite this article: Shafiekhani S, Shafiekhani M, Rahbar S, Jafari AH. Extended robust boolean network of budding yeast cell cycle. J Med Signals Sens 2020;10:94104 
Introduction   
One of the main goals in biology is to discover the fundamental design principles that govern the structures and functions of cells.^{[1],[2]} To have a better understanding of complex biological behaviors, mathematicians in collaboration with biologists have designed computer algorithms and mathematical relations that imitate biological phenomena.^{[3]} Computational models are based on computer algorithms that mimic a natural process with every level of complexity such as models of thymocyte development, biochemical processes, and cell fate specification during Caenorhabditis elegans development.^{[4],[5],[6]} The quantitative relationship between variables (gene activity level and molecular concentration) have been designed to indicate cell signaling pathways in a biologically and physically realistic manner and generate innovative and useful hypotheses.^{[1],[7]} In situ ations whereby the quantitative relationship is unknown, computational models are effective alternatives because they can qualitatively and without experimental data or with missing data describe the biological process and predict quantitative responses.^{[4],[8],[9],[10],[11]} Moreover, since they can be nondeterministic or stochastic, they produce outputs with a range of uncertainties which are natural in biological processes.^{[12],[13]} The goal of computational modeling is to comprehend the general properties of complex networks by quantitative or qualitative terms in order to address the structure of cellular networks and modeling subtle dynamics from molecular levels that contribute to biological functions to intracellular levels that show the average behavior of biological systems.^{[14],[15],[16]}
Biological network modeling has largely been developed by new computational methods that comprise the identification and mathematical definition of restrictions due to cellular regulations as well as biochemical and biophysical laws.^{[17],[18],[19]} Systematic characterization of biochemical reactions and prediction of novel reactions and pathways largely developed with computational methods for attaining automatic reconstruction of regulatory networks.^{[20],[21],[22],[23]} Lee et al. revealed that eukaryotic cellular function in a regulatory network is highly related to a network of transcriptional regulations and they are assessed in the eukaryote saccharomyces cerevisiae, where most of the transcriptional factors are encoded.^{[24]}
Yeast is an organism widely used as a model system and allows the scientist to put the essential experiments into effect.^{[25],[26]} Budding yeast (Saccharomyces cerevisiae) is a simple model system with a single cell eukaryote which grows and then divides into two daughter cells.^{[27],[28]} It is one of the best understood model systems with varieties of tools which are accessible or under development. The cell cycle process consists of four phases: G1, S, G2, and M. There are several checkpoints that guarantee every step in cell cycle has been fully achieved.^{[29]} In the cell cycle process of budding yeast, nearly 800 genes interact with each other. However, the main genes that are responsible for the regulation and dynamic control are very little.^{[30]}
Modeling approaches have demonstrated benefits in unraveling the dynamics of cell proliferation in fission yeast, fruit fly embryos, frog eggs, and budding yeast.^{[31],[32]} Comprehensive literature studies about quantitative modeling of budding yeast cell cycle (BYCC) network have been conducted.^{[33]} For answering many biological questions, such as stem cell control and differentiation, cell commitment (e.g., to apoptosis) and cell cycle progression, the sole prediction of consecutive patterns of states of the control circuit of a cell without detailed information about the sojourn time in each state, or exact time period of the control circuit dynamics could progress our knowledge efficiently.^{[34]} Modeling the path that a cell takes depends on the aim of the study. We did not require the exact time periods of the control circuit dynamics. Therefore, without a knowledge about the numbers of biochemical parameters related to time evolutions which are hardly obtained, a model can be built.^{[34]} Recently, a lot of attention has been paid to mathematical modeling of yeast cell cycle regulation.^{[34],[35],[36]} For instance, Chen et al. constructed a comprehensive mathematical model of the yeast cell cycle with ordinary differential equations (ODE) which consists of many parts such as: equations that describe cyclindependent kinase dynamics, inhibitors of Clbdependent kinase dynamics, Clb degradation dynamics, DNA synthesis, spindle formation, and transcription factor dynamics.^{[37]} In another study conducted by Li et al., they constructed a Boolean network (BN) of BYCC and compared it with a random network and mathematically demonstrated that the biological network of BYCC is extremely stable. Hence, robustness of system with respect to little perturbations and noise will be preserved.^{[30]} BN is widely used for studying the robustness and stability of biological systems.^{[30],[34],[38]} The main disadvantage of BN is that they cannot be used to model large networks.^{[4]} A study conducted by Mura and CsikászNagy defined a stochastic Petri net (PN) model of BYCC from deterministic ODE model.^{[36]} PN is very appropriate for the quantitative and qualitative modeling of concurrent, asynchronous and distributed systems and widely used for modeling of biological networks.^{[39],[40],[41],[42],[43],[44]} Steggles et al. introduced a method for automatic translation of BN into PN model.^{[45]} Zhang et al. created a stochastic model of BYCC and defined 2048 cell states whose transition probabilities between these states are affected by noise parameter. They found out that when the noise (perturbation) is larger than the amount of the order of interaction intensity, the network dynamics rapidly becomes noise dominating and unstable.^{[33]}
Randomness in gene regulation and expression is a natural subject in biological systems.^{[46],[47],[48],[49],[50]} It includes cell cycle since stochastic models that can capture this randomness matches some characteristics of the cell cycle that cannot be achieved with the deterministic models.^{[51],[52]} ODE models of the BYCC has many kinetic parameters that can hardly be approximated. These models are deterministic and cannot capture randomness in gene interactions. Similarly, BN is deterministic. Despite the fact that previous stochastic models of BYCC can model randomness in gene regulation, however, we will consider cell cycle phases as states of Markov model (MM) and we will construct a new computational model. Hence, we can compute the transition probabilities between cell cycle phases. In addition, we will consider apoptosis state in cell cycle according to biological evidences.
Estimating the weight of interactions between genes or proteins as kinetic parameters is an interesting aspect for researchers.^{[53],[54]} Moles et al. estimated the parameters of a biochemical pathway using optimization methods.^{[55]} Estimation of kinetic parameters can be done with global optimization methods. These methods can be classified into deterministic and stochastic strategies. Stochastic methods have a weak assurance of convergence to the global optima. There are many stochastic methods such as: Adaptive random search, clustering methods, multistart methods (that identify the vicinity of local optima) and evolutionary computation (biologically inspired methods which uses the notion of reproduction, mutation, and principle of survival in finding the best solution by iteratively constructing new generation). Evolutionary computation methods can be classified into three: evolutionary programming, evolutionary strategies, and genetic algorithm (GA).^{[56]} In this study, we aim to optimize the kinetic parameters of BYCC by GA that fully will be described in the next sections.
In this section, the importance of modeling in biology is pointed out and some computational models of the yeast cell cycle have been reviewed. In this study, our aim was to construct an extended BN based on an existing BN of BYCC to achieve more robustness and also to construct an MM for analyzing the state space of cell cycle network which results from extended BN. By extended BN, the dynamics of protein activities that control cell cycle network can be achieved while the MM is constructed to find the dynamics of transition probabilities of the cell between its phases (states of MM).
Structure
Our model is constructed according to BN of a protein interaction model of Li et al., which is depicted in [Figure 1]. This model contains four classes of members: (1) transcription factors (SBF, MBF, complex Mcm1/SFF, Swi5); (2) the inhibitors, degrades of the cyclin/Cdc28 complexes (Cdc20, Cdc14, Cdh1, Sic1); (3) cyclins (Clb1, 2, 5, and 6, as well as Cln1, 2, and 3 which bind to the kinase Cdc28); and (4) checkpoints (the cell size, DNA replication and injury, and spindle assembly).^{[30]} SBF and MBF are important cell cycle transcription factors. It has been proven that either SBF or MBF by cyclin Clb5 is enough for budding yeast cells to initiate DNA replication and duplication of the spindle pole body. The activation of SBF and MBF initiates the synthesis of Clb5,6, and Cln1,2 cyclins, and stimulates the transcription of G1/S genes (Cln1,2, Cb5,6). Afterwards, the following actions by other factors occur respectively: DNA synthesis, nucleus movement to the bud neck, constitution of an intranuclear spindle and a gathering of repeated chromosomes to the metaphase plate. Monitoring successful completion of these events is performed at the metaphase checkpoint. Anaphase, telophase, and cell division will go after metaphase respectively. In the S phase, cyclin Clb5 begins DNA replication, then G2/M genes such as Clb2 are transcribed by means of activation of transcription factor complex Mcm1/SFF. Finally, with inhibition and degradation of Clb2 by Cdc20, Cdh1 and the Sic1 cell will exit from mitosis into two cells.  Figure 1: The cell cycle network of budding yeast. Nodes represent proteins or protein complexes, arrow represent activation effect and dotted arrow represent inhibition effect
Click here to view 
BNs are deterministic computational models. They were first represented by Kauffman in 1970. BNs estimate the dynamics of biological networks by assigning active state (“on” or logic 1) or inactive state (“off” or logic 0) to each molecule (e.g., gene or protein) while intermediate expression levels are neglected. State of each molecule is active if the sum of its activation becomes larger than its inhibition and is inactive if the sum of its inhibition becomes larger than its activation.
In [Figure 1], arrows and dotted arrows represent positive regulation (activation effect) and negative regulation (deactivation, repression, inhibition, or degradation) respectively. In [Figure 1], 11 key regulators that are responsible for the regulation and control of this complex network are considered.^{[30]} The dynamics of cell cycle network is modeled by assigning active state (“on” or logic 1) or inactive state (“off” or logic 0) to each of 11 nodes. Therefore, each of the 11 nodes namely MBF, SBF, Cln3, Cdh1, Swi5, Mcm1/SFF, Cln1, Sic1, Clb5, Clb1, and Cdc20 has two states, S_P_{i} = 1 (state of i^{th} protein is 1 or protein is active) and S_P_{i} = 0 (state of i^{th} protein is zero or protein is inactive). Since the model has 11 binary nodes, it has S_C={0,1}^{11}_{i} state spaces. Cell state is a row vector with 11 elements and 2048 (2^{11}) different initial states for executing the rule 1. The cell state can be determined as: S_C= (Cln3, MBF, SBF, Cln1, Cdh1, Swi5, Cdc20, Clb5, Sic1, Clb1, Mcm1) so that each element of vector can be zero or one.
In BN, the time evolution of protein states is determined via the following rule (30):
In BN the value of W_{ij} is equal to 1 or1 for arrows and dotted arrows in [Figure 1]. In the extended BN, W_{ij} represents the weight of interaction of protein j to protein i. It represents stoichiometric coefficients of a biochemical reaction between protein (or gene) j and protein (or gene) i in BYCC. In [Figure 1], arrows represent activation of protein i by protein j and the value of W_{ij} is positive while the dotted arrows represents inhibition of protein j by protein i, and the value wij is negative. For each of the initial states of extended BN, we have implemented the rule 1. For all arrows, W_{ij} is a positive integer ϵ (1, 10), and in all dotted arrows, wij is a negative integer ϵ (−10, −1), until the cell states become stable and reach the steady state. By further execution of rule 1, the cell state does not change. Li et al. have shown that special biological phases of cells in the BYCC (G1 phase, S phase, G2 phase and M phase) correspond to the special states of proteins.^{[30]} We used the MM to investigate the transitions of cell phases due to extended BN during the dynamics of the BYCC.
We defined six states: (1) stationary state equal to stationary G1 state (SG1), (2) G1 state (Gap 1 phase), (3) S state (synthesis phase), (4) G2 state (Gap 2 phase), (5) M state (Mitosis phase) and (6) “other” state respectively as first to sixth states of MM. Stationary state corresponds to stationary G1 while “other” state represents the state of cell at some point between both of the five predefined states like when the cell evolves from G2 to M phase. Before reaching M phase, the state of cell is not the same as G2 and M state, it corresponds to “other” state. This MM has six states. An addition of apoptosis state to MM makes a total of 7 states. BN model of yeast cell cycle was introduced by Li et al. It consists of 34 arrows between genes or proteins that describe their interaction. Li et al. considered same intensity for all interactions. Our aim to optimize the intensity or weight of these interactions Wij≠0 for i, j =1, 2,....,11 (to find extended BN) by GA led to maximum dedicated probability of a special biological pathway that can be computed with MM (with seven states) and to achieve more stable extended BN. The GA is based on an interpretation of Darwinian evolution. In this approach, a population of individuals was created. Each individual is indicated by a “chromosome” that encodes a possible answer of the optimization problem. Each chromosome contains some genes. The number of genes is equal to the number of unknown parameters that must be estimated. We defined a fitness function, as a measure of the correctness of the answer. By applying the data encoded on the chromosome, we calculated the value of fitness function. The fitness of the individual represents the likelihood that the individual genes advance to the next generation. The value of each gene determines the value of the corresponding unknown parameter. The evolution of solution can be done similar to real organisms, with techniques such as “crossover” and “mutation”.^{[56]} Therefore, the weight of interactions between genes or proteins of BYCC extended BN can be estimated by GA.
We considered the population of chromosomes with each chromosome having 34 genes. Each gene represents one of the kinetic parameters (the value of W_{ij}) of yeast cell cycle extended BN model of Li et al. in [Figure 1]. Atfirst, we considered a vector (with length 34) of random integer numbers in a range ( [1, 10] for activation effects and [−10, −1] for inhibition effects) of chromosomes. Then, by calculating the value of fitness function computed with MM and extended BN corresponding to each chromosome and applying GA techniques such as crossover and mutation of best chromosomes, we found the best solution. The fitness function is a summation of the subsequent probability of transition between cell cycle phases (to maximize occurrence of the G1>S>G2>M>G1>SG1 cycle) and the relative basin size of SG1 as the major attractor of extended BN. Therefore the fitness function of GA is:
P(G1 > S > G2 > M > G1 > SG1) + normalized basin size of SG 1. For creating an MM with six states (and then add the seventh state as apoptosis), the following actions must be done. By executing rule 1 for each of the 2048 initializations until cell state reaches steady state. Afterwards, checking that the cell states are equal to which of the six predefined states (SG1, G1, S, G2, M, and “other” state).
Hence, the data sequence {x^{(n)}} containing subsequent cell states are taken from a deterministic extended BN model with weights proposed by GA results.
Therefore, this data sequence has the form:
{x^{(i)}} = Cell state for the i^{th} initialization (until steady state), i = 1, 2,3..., 2048 In the data sequence, finding the transition frequency rates F_{ij} is done by counting the number of transitions from state i to j in one step. For the sequence {x^{(n)}}, one can create the onestep transition matrix F and the estimation for P_{ij} as follows:
Therefore, each element pij of the matrix p represents the probability of transition from state i to j in one step. Actually, this probability is as a result of using the concept of frequency rate.
So far, we have assumed six cell states according to deterministic extended BN. When the cell is arrested in each of the cell cycle checkpoints, it may undergo apoptosis. Therefore, the seventh state of MM as apoptosis state must be added and MM is improved to seven states. According to [Figure 2] and [Figure 3], after G1 phase, the cell may evolve to S or to “other” state and then to S, or may be arrested in G1 by checkpoint that checks DNA damage. When arrested, the cell may undergo apoptosis or recover and evolve to next state, S. After S phase, the cell may evolve to G2, or to “other” state and then to G2, or may be arrested in S by checkpoint that checks correct replication of DNA. When arrested, the cell may undergo apoptosis or recover and continue to the next state, G2. After G2 phase, the cell may evolve to M, or to “other” state and then to M, or may be arrested in G2 by DNA damage checkpoint. When arrested, the cell may undergo apoptosis or recover and evolve to the next state, M. Finally, after M phase, the cell may evolve to G1, or to “other” state and then to G1, or may be arrested in M because of improper spindle formation. When arrested, the cell may undergo apoptosis or recover and evolve to the next state, G1. According to the result of deterministic extended BN model of yeast cell cycle, the stationary state is only accessible after G1 state. Therefore, after G1 state, the next state may be stationary G1. The i^{th} state is named absorbing state of Markov chain if P_{ij} in transition probability matrix, and actually, once the system hits i^{th} state, it loses ability to escape. Therefore, stationary G1 state and apoptosis state are absorbing states of Markov chain because the stationary state is an attractor of this system. After attaining this state, cell state will not change and after apoptosis, the cell dies. At the absorbing state, a new routine of cell cycle program will start. Therefore, these elements of transition probability matrix must be zero. The transition probability matrix P is given as follows:
P12 = P13 = P14 = P15 = P16 = P17 = P72 = P73 = P74 = P75 = P76 = 0
 Figure 2: Cell cycle checkpoints, in certain places of cell cycle may cell to be arrested and then undergo apoptosis or evolve to next phase
Click here to view 
 Figure 3: Seven states of Markov model and possible interactions between states. Apoptosis state and SG1 state are absorbing states of Markov chain
Click here to view 
Some elements of this transition probability matrix must be zero because they represent an unusual phenomenon in biology. For example, the probability of transition from S to G1 is zero. This probability matrix P which results from the definition of six predefined states do not have apoptosis state. Therefore, we built a new stochastic chain and considered the apoptosis state according to transition probability matrix p. The next state of the cell during time steps was determined according to the structure of extended BN that deterministically regulate biochemical interactions between genes or proteins, can be affected by randomness.^{[33],[46],[47],[48],[49]} Randomness in protein interactions by transition probability matrix is as a result of deterministic extended BN and roulette wheel selection technique
By assigning the numbers one to six to states: SG1, G1, S, G2, M and “other” respectively, we have state space A = {1, 2, 3, 4, 5, 6}. If the present state in stochastic Markov chain corresponds to i∈{1, 2, 3, 4, 5, 6}, creating a vector V_{i} of nonzero elements of the i^{th} row of matrix P is as follow: j= 1, 2,...,6 i = 1, 2,...,6n_{j} = number of elements of vector V_{i}.
Then, generating a random number between minimum and maximum of the numbers of the vector V_{i}, which is pE7 (that is (n_{E}+1 an element of the vector V_{i}) and then normalizing the numbers of the vector V_{i} to their sum, thereby generating the numbers between zero and one that are elements of the vector V_{i} that has (n_{E}+1 elements.
Then creating a roulette wheel with (n_{E}+1 sections in which the length of each section is equal to probabilities of vector V_{i} and generating a random number with the uniform distribution between zero and one. As shown in [Table 1], by examining that this number corresponds to which section of the roulette wheel, the state corresponding to that section wins and next state of the cell is determined.
For example, if the present state in stochastic chain is S (corresponds to state 3). For creating probability of transition from S state to apoptosis state, we first generated a random number between minimum and maximum of the numbers such that p 37 normalizes the numbers to their summation. Thereby generating the numbers between zero and one that are elements of vector:.
That represent the probability of transition from S to G2, “other” and finally to apoptosis state in one step, then creating a roulette wheel with four sections (equal to the number of elements of the vector followed by generating a random number with a uniform distribution between zero and one and examining that this section of a roulette wheel and state corresponding to that section wins and the next state of the cell is determined.
Generally, if the current state of the cell is G1, S, G2 or M, exactly similar to predefined procedure, we respectively compute these probabilities and implement roulette wheel selection technique. Therefore, repeating this procedure results into a new data sequence {Y^{(n)}}. This data sequence or stochastic chain is created based on the notion that the next state of the cell is determined stochastically according to the current state of the cell while the data sequence {X^{(n)}} is as a result of the deterministic rules of extended BN. For the data sequence {Y^{(n)}}, one can create the onestep transition probability matrix P_{2 } using rule 2. Transition probability matrix P_{2} is different from transition probability matrix P because it regards apoptosis as the 7^{th} state and describes probabilities of transition between all states to apoptosis. Furthermore, this matrix is as a result of randomness in gene regulation while the transition matrix P is due to deterministic extended BN. Transition matrix P_{2} can be used for constructing the fitness function of GA. The network's evolution has Markov property because statistically, cell state at the next time step only depends on cell state at the present time step. The stochastic process presumed to be timehomogeneous. Since all states are not accessible from each other, Markov chain is nonreducible. The time steps are logical steps that describe causality rather than actual time steps. Under the Markov assumption, transition probability of the Markov chain {Y^{(n)}} is computed. Therefore, transition matrix P_{2} according to Markov assumption is given as:
Our aim was to extend an existing BN model of Li et al. by optimizing its kinetic parameters according to one criterion computed from MM introduced in this paper to produce a gold standard behavior. The purpose of gold standard behavior is that the probability of special biological pathway is maximized, and based on attraction of stationary G1 state (S_G1), the biggest attractor of BYCC network is increased.
Let be a Markov chain with state space S ∈{1, 2, 3, 4, 5, 6, 7} and transition probability matrix p. According to MM, the probability of transition from state i to j after n steps provided that this transition has a certain way to win is computed via rule 4:
Therefore, we find the probability of special pathway in the BYCC via rule 4. The known pathway is the subsequent transition from G1 to S, S to G2, G2 to M and finally G1 to stationary G1. These probabilities from transition matrix p 2 are obtained as:
Finding the probabilities of the first arrival to each of these seven states from other states can be interesting for biologists to compute using MM rules. According to MM, the first passage time from the state i to j is the number of time steps, Tij required to attain state j for the first time, given that the chain was initially in state i. The first passage time is a random variable with probability distribution defined with rule 6. It can be computed recursively using rule 7 as follows:
The kinetic parameters of extended BN of the cell cycle are optimized by GA as illustrated in [Table 2]. By applying these kinetic parameters to the network, the relative basin size of stationary G1 state increased from 86% to 96.48% while the number of attractors decreased from 7 to 5. The probability of special pathway in the BYCC via rule 3, which is a known pathway and is the subsequent transition from G1 to S, S to G2, G2 to M, M to G1, and finally G1 to stationary G1 is maximized. This represents the increment of robustness of the system. Despite the randomness and stochastic gene regulations modeled by MM which enables random cell state selection, the convergence to stationary G1 state as seen in [Table 3] shows five fixed points of cell cycle network and their basin of attractions. Execution of deterministic extended BN begins from each of the 2048 initial states using the optimized parameters of [Table 2] results to one of the five attractors. This extended BN is extremely robust because in 1976 cases from 2048 different initializations of extended BN, cell state reaches stationary G1 state (relative basin size is 96.48%).  Table 2: The kinetics parameters of budding yeast cell cycle network optimized by genetic algorithm for reaching to robust structure with maximum probability of special pathway in the cell cycle
Click here to view 
 Table 3: The fixed points of the cell cycle network and their basin sizes
Click here to view 
[Figure 4], shows the histogram of the incidence of each of the cell states in the extended BN of the cell cycle with optimized kinetic parameters for all 2048 initializations until cell state reaches steady state. The X and Y axes represent the decimal equivalent of the protein states in the network and the number of incidence of each cell state respectively. The decimal equivalent number is attained by the composition of the 11bit vector number that each of its elements corresponds to an active state or inactive state of the proteins in the network. The two maximums of [Figure 4] depict G1 and SG1 states respectively and are subsequent states in cell cycle process. [Figure 4] shows that in 1976 from 2048 initializations of cell cycle processes, the subsequent transition of G1 to SG1 and trap of process in SG1 is seen.  Figure 4: Histogram of occurrence each of cell states in extended BN of budding yeast cell cycle with optimized kinetic parameters for all 2048 initialization until cell state reaches to steady state. X axis and Y axis represent the states and their histograms, respectively
Click here to view 
As illustrated in [Figure 2], the checkpoints which control decency and integrity of the cell cycle progression may interrupt the cell cycle in certain phases and arrest the cycle.^{[57]} Hence, the cell may undergo apoptosis or may recover and evolve to next phase of cell cycle if the problem is fixed.
The probability of the first arrival to each of the seven states after exiting from other states by rule 5 can be computed. [Figure 5] depicts the probabilities of the first arrival to apoptosis state after exit from each of the G1, S, G2 and M states while [Figure 6] depicts the probabilities of the first arrival to SG1 state after exit from each of the G1, S, G2 and M states. These probabilities decrease with time. [Figure 5] and [Figure 6] also depicts the dynamics of the probability of the first arrival to apoptosis (to SG1) state after exiting other states. By these probabilities, we can predict the behavior of BYCC network. As future work, you can apply some interventions to the structure of BYCC network and find its relation with these probabilities or identifying the best intervention in the structure of abnormal (unregulated) mammalian cell cycle network that cause increment in the probability of first arrival to apoptosis, therefore, the proliferation of cells will be prevented.  Figure 5: X axis represents the time steps to reach apoptosis state after exit from each of the four states: G1, S, G2 and M. Y axis represents the probability of first arival to apoptosis state after exit from each of the G1, S, G2, M respectively
Click here to view 
 Figure 6: X axis represents the time steps to reach stationary G1 state fter exit from each of the four states: G1, S, G2 and M. Y axis represents the probability of first arival to stationary G1 state after exit from each of the G1, S, G2, M respectively
Click here to view 
Discussion   
In most gene regulatory networks, kinetic parameters are unknown and can hardly be achieved. Obtaining kinetic parameters requires designing expensive experiments. We have designed a procedure to optimize the kinetic parameters of gene regulatory network of BYCC. For this aim and for modeling randomness in gene regulation, we have constructed an MM whose states correspond to four phases of cell cycle (G1, S, G2 and M phases), stationary G1 state and apoptosis state (in each of the checkpoints of cell cycle, the cell may be arrested and undergo apoptosis). In addition to MM, probabilistic BN (PBN) is one of the best methods for modeling intrinsic noise in molecular interactions. PBN is subset of several BNs with various probabilities and predict the next state of each gene as well as the inherent stochastic nature of molecular interactions. In each time step, the state of the target gene is predicted by a specific vector of functions, such that each element of this vector is specified by a specific BN. These vectors for different genes are numerous. Each gene according to particular BNs with specific probabilities predicts its future state. This paper has described a method for modeling molecular interactions based on single BN. By MM, we can model intrinsic noise (randomness) in gene regulation that in deterministic models are not considered. Our stochastic model is constructed based on transition rates derived from deterministic extended BN model. The dynamics of the network is as a result of randomness in gene regulation which causes a cell in each of the cell cycle phases, according to probabilities derived from deterministic extended BN model select its next state, while in deterministic BN model of yeast cell cycle or in ODE model of BYCC, the time evolution of cell states according to deterministic rules as a result of randomness in gene regulations are neglected.
Li et al. introduced a BN model of yeast cell cycle network that genes or proteins interact with each other by deterministic rules and the dynamics of the network during the time steps is achieved.^{[30]} Hence, the cell evolves from one cell cycle phase to the next. They considered all biochemical kinetic parameters of this network to have the same value. In other words, they considered that all genes or proteins interact with each other with same weights. We have optimized the biochemical kinetic parameters of deterministic extended BN model with GA that causes the probability of special biological pathway to be maximized and basin size of the stationary G1 state to rise. Therefore, the occurrence of the G1>S>G2>M>G1>SG1 cycle is maximized. In addition, the number of attractors of deterministic extended BN decreased from 7 to 5 while basin size of stationary G1 state increased from 86% to 96.48%. This results in increasing robustness of extended BN with optimized kinetic parameters. The present study, similar to the studies,^{[30],[33]} was designed to analyze the stability and dynamics of the BYCC network and like these study, which lacks experimental data, according to expert knowledge of genegene or proteinprotein interactions of BYCC constituents simulates the cell cycle behavior and demonstrates the stability of the system.
Zhang et al. in an interesting study constructed a stochastic MM of BYCC.^{[33]} Their stochastic model having 2048 states and in each time step, cell state with certain probabilities evolved to the next state. The BYCC network of their study is similar to our model and consists of 11 proteins (with two active/inactive states), resulting in a total of 2048 states. Their model consist of two noise parameters α and β to describe the level of noise in cell state transitions. If these parameters are infinite (α=β=∞), their probabilistic model will behave similar to the deterministic BN model. They considered randomness in gene regulation by neglecting kinetic parameters and assumed that all biochemical kinetic parameters have the same value except negative feedbacks of genes on themselves that have another weight. In our MM, we considered seven states. In our model, the time evolution of protein states considering apoptosis state was studied.
In an interesting study, Kraikivsk et al. describe the cell cycle process from DNA replication to mitosis (cell division) by a system of differentialalgebraic equations.^{[31]} Their model is constructed by the hypothetical regulatory mechanism of cell cycle control in budding yeast that has good agreement with experimental data. The model is calibrated by data set of observed phenotypes of 257 mutant yeast strains. The calibrated model is used to predict the phenotypes of 30 novel combinations of mutant alleles (the mutations are simulated by parameter changes). Therefore, the aim of their study was to predict the state of genes or proteins in different phenotypes of BYCC by a deterministic model, while the aim of our study is to develop an extended BN model of BYCC for the dynamic and stability analysis of cell cycle constituents (genes or proteins). In our study, we optimized the weight of proteinprotein interactions by GA to create an extended BN of the BYCC that is robustly designed. In general, the purpose of our study was to evaluate the stability of the cell cycle system and to estimate the rate of proteinprotein or genegene interactions to obtain a more robust and stable network that has a reliable function in the context of noise. Noise refers to the uncertainty in the transition from one state to another in the cell cycle process that is simulated by Markov chain.
In another study, Barik et al. constructed a computational model to describe how a cell control cell cycle progression and prevents daughter cell birth with abnormal genetic features. The simulation of the model is performed by deterministic (to capture the events of bud emergence and DNA synthesis) and stochastic methods (to simulate chemical reactions of cell cycle network). Deterministic simulation of the model is carried out by parameter estimation toolkit that generates the ODE of model variables and stochastic simulation is performed using Gillespies' stochastic simulation algorithm. This model is used to simulate variable phenotypes of 20 mutant yeast strains and to investigate the role of feedbacks regulations in the reduction of noise in cell cycle progression. In addition to the difference in modeling strategies of our study with respect to this study, the constituent proteins of the cell cycle network of this study are also different from our study. In present study, we used MM to analyze the state space created by extended BN and all steps of model construction are described in detail. The importance of MM for achieving the dynamics of the transition probability of states in the extended BN are referred to. In addition, the importance of using optimization methods to increase the stability of BN has been mentioned. Although the proposed method was applied to BYCC, it is not confined to modeling and analysis of this network. It can be utilized for modeling and analysis of state space of any signaling pathway or biological network. We believe that the major contribution of our survey, in addition to the offer of the use of MM and GA in biological network modeling, is that a deeper understanding and comprehensive analysis of the biological network can be achieved and besides dynamical behavior of BN, dynamics of transition probabilities between states of BN can be achieved. Finally, by optimizing the kinetic parameters of BNs, extended BNs with more robustness can be achieved.
Acknowledgment
My profound gratitude to my devoted parents for their love and support all through my life. Thank you both for giving me the strength to reach my dreams. My special appreciation also goes to my siblings who have always supported me and lastly to the cancer modeling group at Tehran University of Medical Sciences.
Financial support and sponsorship
None.
Conflicts of interest
There are no conflicts of interest.
Sajad Shafiekhani is a Ph.D student in Physics and Biomedical Engineering department of Tehran University of Medical Sciences(TUMS). His interests are mathematical and computational models in cancer modeling. He works on tumorimmune system by fuzzy logic, stochastic petri net, agent based models, Markov model and boolean networks.
Email: sshafikhani@razi.tums.ac.ir
Mojtaba Shafiekhani is a master student at Amirkabir university of Tehran (Ploytechnic) and graduated with a Bachelor's degree in Electrical Engineering from"the University of"Khaje Nasir Toosi in 2018. For His masterthesis, he is working under supervision of Abbas Nasiraei in the imagereconstruction of magnetic resonance imaging. He built"skills"in"Matlab, python, C and Deep learning.
Email: mojtaba9375@gmail.com
Sara Rahbar is a Ph.D. candidate in Physics and Biomedical Engineering department of Tehran University of Medical Sciences (TUMS). Her project is focused on optimizing combinational immunotherapy protocol by using Agentbased modeling of tumor and the immune system. Her interests are also the application of structural and interactionbased modeling techniques like Game Theory and Petri Nets in biology. She loves spending time with her 6yearold daughter playing with Lego bricks.
Email: sarahbar2002@yahoo.com
Amir Homayoun Jafari is currently an associate professor of Biomedical Engineering in Tehran University of Medical Sciences. He has received his PhD degree in Bioelectric Engineering from AmirKabir University of Technology. His works in the field of modeling are basically focused on cancer modeling, immunodeficiency and neuromuscular system. He is also working on nonlinear analysis methods mainly with the aim of applications in BCI and neuroscience. His interests are signal processing as well chaos and game theory used in biological applications. He enjoys strategy games especially chess.
Email: h_jafari@tums.ac.ir
References   
1.  Kiyoumarsioskouei A, Shamloo A, Azimi S, Abeddoust M, Saidi M. A computational model for estimation of mechanical parameters in chemotactic endothelial cells. Sci Iran Trans B Mech Eng 2016;23:260. 
2.  Waltemath D, Karr JR, Bergmann FT, Chelliah V, Hucka M, Krantz M, et al. Toward community standards and software for wholecell modeling. IEEE Trans Biomed Eng 2016;63:200714. 
3.  Yaghoubi S, Pishevar A, Saidi M, Shirani E. Modeling selfassembly of the surfactants into biological bilayer membranes with special chemical structures using dissipative particle dynamics method. Sci Iran Trans B Mech Eng 2016;23:942. 
4.  Fisher J, Henzinger TA. Executable cell biology. Nat Biotechnol 2007;25:123949. 
5.  Tanay A, Regev A. Scaling singlecell genomics from phenomenology to mechanism. Nature 2017;541:3318. 
6.  Bargaje R, Trachana K, Shelton MN, McGinnis CS, Zhou JX, Chadick C, et al. Cell population structure prior to bifurcation predicts efficiency of directed differentiation in human induced pluripotent cells. Proc Natl Acad Sci U S A 2017;114:22716. 
7.  Davidson EH, Rast JP, Oliveri P, Ransick A, Calestani C, Yuh CH, et al. A genomic regulatory network for development. Science 2002;295:166978. 
8.  Huang L, Jiang Y, Chen Y. Predicting drug combination index and simulating the networkregulation dynamics by mathematical modeling of drugtargeted EGFRERK signaling pathway. Sci Rep 2017;7:40752. 
9.  Wittwehr C, Aladjov H, Ankley G, Byrne HJ, de Knecht J, Heinzle E, et al. How adverse outcome pathways can aid the development and use of computational prediction models for regulatory toxicology. Toxicol Sci 2017;155:32636. 
10.  Shariatpanahi SP, Shariatpanahi SP, Madjidzadeh K, Hassan M, AbediValugerdi M. Mathematical modeling of tumorinduced immunosuppression by myeloidderived suppressor cells: Implications for therapeutic targeting strategies. J Theor Biol 2018;442:10. 
11.  Andreozzi S, Chakrabarti A, Soh KC, Burgard A, Yang TH, Van Dien S, et al. Identification of metabolic engineering targets for the enhancement of 1,4butanediol production in recombinant E. coli using largescale kinetic models. Metab Eng 2016;35:14859. 
12.  Sipahi R, Zupanc GK. Stochastic cellular automata model of neurosphere growth: Roles of proliferative potential, contact inhibition, cell death, and phagocytosis. J Theor Biol 2018;445:15165. 
13.  Lin YT, Hufton PG, Lee EJ, Potoyan DA. A stochastic and dynamical view of pluripotency in mouse embryonic stem cells. PLoS Comput Biol 2018;14:e1006000. 
14.  Zhu X, Yuan R, Hood L, Ao P. Endogenous molecularcellular hierarchical modeling of prostate carcinogenesis uncovers robust structure. Prog Biophys Mol Biol 2015;117:3042. 
15.  Bordbar A, Monk JM, King ZA, Palsson BO. Constraintbased models predict metabolic and associated cellular functions. Nat Rev Genet 2014;15:10720. 
16.  StewartOrnstein J, Chen S, Bhatnagar R, Weissman JS, ElSamad H. Modelguided optogenetic study of PKA signaling in budding yeast. Mol Biol Cell 2017;28:2217. 
17.  Neuert G, Munsky B, Tan RZ, Teytelman L, Khammash M, van Oudenaarden A. Systematic identification of signalactivated stochastic gene regulation. Science 2013;339:5847. 
18.  Cheng TM, Heeger S, Chaleil RA, Matthews N, Stewart A, Wright J, et al. A simple biophysical model emulates budding yeast chromosome condensation. Elife 2015;4:e05565. 
19.  Serohijos AW, Shakhnovich EI. Merging molecular mechanism and evolution: Theory and computation at the interface of biophysics and evolutionary population genetics. Curr Opin Struct Biol 2014;26:8491. 
20.  Briat C, Gupta A, Khammash M. Antithetic integral feedback ensures robust perfect adaptation in noisy biomolecular networks. Cell Syst 2016;2:1526. 
21.  Zickenrott S, Angarica VE, Upadhyaya BB, del Sol A. Prediction of diseasegenedrug relationships following a differential network analysis. Cell Death Dis 2016;7:e2040. 
22.  Teodosiu C, David D, Teodosiu R. Numerical prediction of thermoaeraulic behavior for a cavity with internal linear heat source. Sci Iran Trans B Mech Eng 2016;23:618. 
23.  O'Brien EJ, Monk JM, Palsson BO. Using genomescale models to predict biological capabilities. Cell 2015;161:97187. 
24.  Lee TI, Rinaldi NJ, Robert F, Odom DT, BarJoseph Z, Gerber GK, et al. Transcriptional regulatory networks in Saccharomyces cerevisiae. Science 2002;298:799804. 
25.  Spellman PT, Sherlock G, Zhang MQ, Iyer VR, Anders K, Eisen MB, et al. Comprehensive identification of cell cycleregulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol Biol Cell 1998;9:327397. 
26.  Nurse P. The Josef Steiner lecture: CDKs and cellcycle control in fission yeast: Relevance to other eukaryotes and cancer. Int J Cancer 1997;71:7078. 
27.  Nasmyth K. Viewpoint: Putting the cell cycle in order. Science 1996;274:16435. 
28.  Mendenhall MD, Hodge AE. Regulation of Cdc28 cyclindependent protein kinase activity during the cell cycle of the yeast Saccharomyces cerevisiae. Microbiol Mol Biol Rev 1998;62:1191243. 
29.  Hartwell LH, Weinert TA. Checkpoints: Controls that ensure the order of cell cycle events. Science 1989;246:62934. 
30.  Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cellcycle network is robustly designed. Proc Natl Acad Sci U S A 2004;101:47816. 
31.  Kraikivski P, Chen KC, Laomettachit T, Murali TM, Tyson JJ. From START to FINISH: Computational analysis of cell cycle control in budding yeast. NPJ Syst Biol Appl 2015;1:15016. 
32.  Ahmadian M, Tyson J, Cao Y, editors. A Stochastic Model of Size Control in the Budding Yeast Cell Cycle. Proceedings of the 2018 ACM International Conference on Bioinformatics, Computational Biology, and Health Informatics, ACM; 2018 
33.  Zhang Y, Qian M, Ouyang Q, Deng M, Li F, Tang C. Stochastic model of yeast cellcycle network. Phys D Nonlin Phenom 2006;219:359. 
34.  Davidich MI, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PLoS One 2008;3:e1672. 
35.  Barik D, Baumann WT, Paul MR, Novak B, Tyson JJ. A model of yeast cellcycle regulation based on multisite phosphorylation. Mol Syst Biol 2010;6:405. 
36.  Mura I, CsikászNagy A. Stochastic petri net extension of a yeast cell cycle model. J Theor Biol 2008;254:85060. 
37.  Chen KC, CsikaszNagy 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:36991. 
38.  Chaves M, Sontag ED, Albert R. Methods of robustness analysis for boolean models of gene control networks. Syst Biol (Stevenage) 2006;153:15467. 
39.  Liu F, Heiner M, Gilbert D. Coloured petri nets for multilevel, multiscale and multidimensional modelling of biological systems. Brief Bioinform 2019;20:87786. 
40.  Liu F, Heiner M, Yang M. Fuzzy stochastic petri nets for modeling biological systems with uncertain kinetic parameters. PLoS One 2016;11:e0149674. 
41.  Werthmann B, Marwan W. Developmental switching in Physarum polycephalum: Petri net analysis of single cell trajectories of gene expression indicates responsiveness and genetic plasticity of the Waddington quasipotential landscape. J Phys D Appl Phys 2017;50:464003. 
42.  Herajy M, Liu F, Heiner M. Efficient modelling of yeast cell cycles based on multisite phosphorylation using coloured hybrid petri nets with markingdependent arc weights. Nonlin Anal Hybrid Syst 2018;27:191212. 
43.  Laomettachit T, Chen KC, Baumann WT, Tyson JJ. A Model of yeast cellcycle regulation based on a standard component modeling strategy for protein regulatory networks. PLoS One 2016;11:e0153738. 
44.  Castro C, Flores DL, CervantesVásquez D, VargasViveros E, GutiérrezLópez E, MuñozMuñoz F. An agentbased model of the fission yeast cell cycle. Curr Genet 2019;65:193200. 
45.  Steggles LJ, Banks R, Wipat A, editors. Modelling and Analysing Genetic Networks: From Boolean Networks to Petri Nets. International Conference on Computational Methods in Systems Biology, Springer; 2006. 
46.  Kaufmann BB, van Oudenaarden A. Stochastic gene expression: From single molecules to the proteome. Curr Opin Genet Dev 2007;17:10712. 
47.  Elowitz MB, Levine AJ, Siggia ED, Swain PS. Stochastic gene expression in a single cell. Science 2002;297:11836. 
48.  McAdams HH, Arkin A. Stochastic mechanisms in gene expression. Proc Natl Acad Sci U S A 1997;94:8149. 
49.  Paulsson J. Models of stochastic gene expression. Phys Life Rev 2005;2:15775. 
50.  Soltani M, VargasGarcia CA, Antunes D, Singh A. Intercellular variability in protein levels from stochastic expression and noisy cell cycle processes. PLoS Comput Biol 2016;12:e1004972. 
51.  IyerBiswas S, Wright CS, Henry JT, Lo K, Burov S, Lin Y, et al. Scaling laws governing stochastic growth and division of single bacterial cells. Proc Natl Acad Sci U S A 2014;111:159127. 
52.  Ghusinga KR, VargasGarcia CA, Singh A. A mechanistic stochastic framework for regulating bacterial cell division. Sci Rep 2016;6:30229. 
53.  Liepe J, Kirk P, Filippi S, Toni T, Barnes CP, Stumpf MP. A framework for parameter estimation and model selection from experimental data in systems biology using approximate bayesian computation. Nat Protoc 2014;9:43956. 
54.  Gábor A, Banga JR. Robust and efficient parameter estimation in dynamic models of biological systems. BMC Syst Biol 2015;9:74. 
55.  Moles CG, Mendes P, Banga JR. Parameter estimation in biochemical pathways: A comparison of global optimization methods. Genome Res 2003;13:246774. 
56.  Shaw O, Steggles J, Wipat A. Automatic parameterisation of stochastic petri net models of biological networks. Electron Notes Theor Comput Sci 2006;151:11129. 
57.  Scialdone A, Natarajan KN, Saraiva LR, Proserpio V, Teichmann SA, Stegle O, et al. Computational assignment of cellcycle stage from singlecell transcriptome data. Methods 2015;85:5461. 
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6]
[Table 1], [Table 2], [Table 3]
