Extended Robust Boolean Network of Budding Yeast Cell Cycle.

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. Copyright:

Introduction

One of the main goals in biology is to discover the fundamental design principles that govern the structures and functions of cells.[12] 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.[456] 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.[17] In situations 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.[4891011] Moreover, since they can be nondeterministic or stochastic, they produce outputs with a range of uncertainties which are natural in biological processes.[1213] 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.[141516]

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.[171819] Systematic characterization of biochemical reactions and prediction of novel reactions and pathways largely developed with computational methods for attaining automatic reconstruction of regulatory networks.[20212223] 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.[2526] Budding yeast (Saccharomyces cerevisiae) is a simple model system with a single cell eukaryote which grows and then divides into two daughter cells.[2728] 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.[3132] 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.[343536] 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 cyclin-dependent kinase dynamics, inhibitors of Clb-dependent 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.[303438] The main disadvantage of BN is that they cannot be used to model large networks.[4] A study conducted by Mura and Csikász-Nagy 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.[394041424344] 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.[4647484950] 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.[5152] 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.[5354] 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, multi-start 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 intra-nuclear 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

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 = 1 (state of ith protein is 1 or protein is active) and S_P = 0 (state of ith protein is zero or protein is inactive). Since the model has 11 binary nodes, it has S _ C={0,1}11 state spaces. Cell state is a row vector with 11 elements and 2048 (211) 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):

s

In BN the value of W is equal to 1 or-1 for arrows and dotted arrows in Figure 1. In the extended BN, W 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 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 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) of yeast cell cycle extended BN model of Li et al. in Figure 1. At first, 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 SG1. 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(} 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(} = Cell state for the ith initialization (until steady state), i =1, 2,3..., 2048 In the data sequence, finding the transition frequency rates F is done by counting the number of transitions from state i to j in one step. For the sequence {x(}, one can create the one-step transition matrix F and the estimation for P 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 Figures 2 and 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 ith state is named absorbing state of Markov chain if P in transition probability matrix, and actually, once the system hits ith 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:

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

Figure 3

Seven states of Markov model and possible interactions between states. Apoptosis state and SG1 state are absorbing states of Markov chain

P12 = P13 = P14 = P15 = P16 = P17 = P72 = P73 = P74 = P75 = P76 = 0

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.[3346474849] 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 of nonzero elements of the ith row of matrix P is as follow: j= 1, 2,...,6 i = 1, 2,...,6 n number of elements of vector V.

Then, generating a random number between minimum and maximum of the numbers of the vector V, which is pi7 (that is (n + 1) an element of the vector V) and then normalizing the numbers of the vector V to their sum, thereby generating the numbers between zero and one that are elements of the vector V that has (n +1) elements.

Then creating a roulette wheel with (n + 1) sections in which the length of each section is equal to probabilities of vector V 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.

Table 1

Roulette wheel selection technique

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 {p33, p34, p36, p37} 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(}. 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(} is as a result of the deterministic rules of extended BN. For the data sequence {Y(}, one can create the one-step transition probability matrix P using rule 2. Transition probability matrix P is different from transition probability matrix P because it regards apoptosis as the 7th 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 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 time-homogeneous. Since all states are not accessible from each other, Markov chain is non-reducible. 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(} is computed. Therefore, transition matrix P 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 p2 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:

Results

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

W1=2	W6=3	W11=10	W16=4	W21=10	W26=7	W31=2
W2=2	W7=9	W12=5	W17=8	W22=6	W27=6	W32=1
W3=3	W8=6	W13=6	W18=1	W23=10	W28=1	W33=5
W4=6	W9=10	W14=10	W19=10	W24=7	W29=6	W34=9
W5=10	W10=10	W15=2	W20=4	W25=8	W30=8

Table 3

The fixed points of the cell cycle network and their basin sizes

Number of attractors	Basin size	Cln3	MBF	SBF	Cln1, 2	Cdh1	Swi5	Cdc20	Clb5, 6	Sic1	Clb1, 2	Mcm1
1	1976	0	0	0	0	1	0	0	0	1	0	0
2	60	0	0	1	1	0	0	0	0	0	0	0
3	7	0	0	0	0	0	0	0	0	1	0	0
4	4	0	0	0	0	0	0	0	0	0	0	0
5	1	0	0	0	0	1	0	0	0	0	0	0

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 11-bit 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

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. Figures 5 and 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

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

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,[3033] 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 gene-gene or protein-protein 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 differential-algebraic 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 protein-protein 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 protein-protein or gene-gene 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.

Financial support and sponsorship

None.

Conflicts of interest

There are no conflicts of interest.

BIOGRAPHIES

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: s-shafikhani@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 interaction-based modeling techniques like Game Theory and Petri Nets in biology. She loves spending time with her 6-year-old 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