Stochastic robustness and relative stability of multiple pathways in biological networks
Multiple dynamic pathways always exist in biological networks, but their robustness against internal fluctuations and relative stability have not been well recognized and carefully analyzed yet. Here we try to address these issues through an illustrative example, namely the Siah-1/beta-catenin/p14/19 ARF loop of protein p53 dynamics. Its deterministic Boolean network model predicts that two parallel pathways with comparable magnitudes of attractive basins should exist after the protein p53 is activated when a cell becomes harmfully disturbed. Once the low but non-neglectable intrinsic fluctuations are incorporated into the model, we show that a phase transition phenomenon is emerged: in one parameter region the probability weights of the normal pathway, reported in experimental literature, are comparable with the other pathway which is seemingly abnormal with the unknown functions; whereas, in some other parameter regions, the probability weight of the abnormal pathway can even dominate and become globally attractive. The theory of exponentially perturbed Markov chains is applied and further generalized in order to quantitatively explain such a phase transition phenomenon, in which the nonequilibrium “activation energy barriers” along each transiting trajectory between the parallel pathways and the number of “optimal transition paths” play a central part. Our theory can also determine how the transition time and the number of optimal transition paths between the parallel pathways depend on each interaction’s strength, and help to identify those possibly more crucial interactions in the biological network.
Biological entities have to function robustly against internal stochastic fluctuations for the sake of survival and evolution (2); (1). At the system level, robustness is a complex and emergent property, which could not be simply inferred through only examining the characteristics of each isolated component.
Hence, mathematical models have always been used to investigate the robustness of biological networks. Local stability can be predicted by the deterministic model: each stable state of the biological networks corresponds to an attractor in the state space, which means trajectories of the system originating from certain initial conditions nearby move towards it as time increases. However, it is common for a deterministic dynamical system to have multiple attractors, each of which has its own basin of attraction, defined as the set of initial states approaching that attractor asymptotically. Therefore, how to quantify the relative stability between different attractors becomes an important issue.
Internal stochastic fluctuations are unavoidable in biological networks. After stochasticity is incorporated, the relative stability can be quantitatively addressed according to the comparison of their probabilities, rather than attributing to the magnitudes of their own basins of attraction (3).
There are mainly two types of attractors: fixed points and limit cycles, which correspond to steady states and oscillatory pathways respectively in biological networks. The relative stability of multiple steady states have been carefully studied in many different stochastic models recently, and a nonequilibrium landscape theory has been proposed (4). However, the relative stability and stochastic robustness of parallel pathways have not been comparably recognized and analyzed, which indeed has already caused much interests from biologists (5).
In the present letter, we take the Siah-1/beta-catenin/p14/19 ARF loop of protein p53 dynamics as an illustrative example (Fig. 1)(6). The p53 protein after activated can positively regulate the transcription of the ubiquitin ligase Siah-1, which turns to degrade beta-catenin protein. On the other hand, the increasing of the beta-catenin level can activate the p14/19 ARF gene, which in turn negatively control MDM-2, resulting in further activation of p53 protein. We build both deterministic and stochastic Boolean network models for such a typical biological network, in which two parallel pathways emerge after the cell is damaged. We shall show that in different parameter regions, the relative stability of the parallel pathways behaves very differently, analogous to phase transition in statistical mechanics. These phenomenon can be well explained by a generalization of the theories of exponentially perturbed Markov chains.
Here, we apply the approach of Boolean network model, which has already been used to explore general and global properties of large gene expression networks (8); (7); (9); (13); (14); (15); (10); (11); (12). The Boolean description of cellular states is due to the on-off characteristics of the components in the biological networks. The advantage of such a qualitative model is that the predicted results are often not very sensitive to unknown kinetic parameters (16).
Recently, the stability of several protein p53 pathways has been investigated through stochastic Boolean network models (17), in which the unique pathway predicted in the deterministic model is globally attractive and has a dominant cycle flux in the stochastic model, which is the natural generalization of the deterministic limit cycle. Here we apply the same Boolean network model to the biological diagram in Fig. 1(a), where the sequence of the nodes is (p53, Siah-1, -catenin, p14/19 ARF, Mdm-2). In such a Boolean network model, each involved key factor (e.g. proteins, DNAs or RNAs) in the biological network of interest only has two states, or , representing its active and inactive states respectively. The states in the next time step are determined by the present states via the following rules (13); (10); (11); (17); (15):
Denote a matrix , in which (activation) for an arrow from protein to protein in Fig. 1 (a), and (inhibition) for a horizontal bar instead of arrowhead from protein to protein in Fig. 1 (a). Other elements are set as zero. Given , the deterministic model can be described as follows: if , then
where the function , and is the input to the -th node.
Another parameter indicates whether the cell is under normal environment () or suffered from certain external perturbation (signal) (), e.g. DNA damage. According to the diagram in Fig. 1 (a), inside a normal cell, the level of p53 protein is quite low, which in turn results in a low expression of Siah-1; the inhibition of Siah-1 then activates the -catenin as well as p19/14ARF, leading to the suppression of Mdm-2. Hence, the state should be a fixed point of the deterministic dynamics with , i.e. if . Once the cells encounter certain external damage signal as indicated by , the protein p53 is activated first, i.e. if .
|Time||p53||Siah-1||-catenin||p14/19 ARF||Mdm-2||Node no.|
|Time||p53||Siah-1||-catenin||p14/19 ARF||Mdm-2||Node no.|
There are always two parallel attractors in the deterministic models, no matter or (Fig. 1 (b) and (c) with node numbers 0-31 corresponding to the 32 states here by means of replacing “0” with “-1” in their respective binary representation). When , the two attractors are both limit cycles (see tables 1 and 2 for cycle 1 and 2 respectively). The magnitude of the two cycles’ attractive basins is comparable to each other (12 and 20 nodes respectively), preventing us from determining their relative stability at this stage. When , Cycle 1 in the deterministic model is replaced by the fixed point (Node number 6).
However, in fact only Cycle 1 matches the reported Siah-1/beta-catenin/p14/19 ARF loop of the p53 pathways in the biological literature (6). Therefore, we must appeal to the stochastic approach if we would like to understand the relative stability of the parallel attractors and the transition time between them.
if , and or ,
if , and or , and
if . Here .
To better illustrate the behavior of the stochastic Boolean network model in the limit of vanishing fluctuation, we set in advance that and ( and are nonnegative) when the stochastic model approaches the deterministic counterpart with , as , and are all very large. The cycle fluxes, which is the average times of occurrences of parallel pathways can be calculated accurately according to the expression in (18). As long as , the net flux of Cycle 1 is about regardless of the value of given a fixed and sufficiently large value of (Fig. 2 (a)). Hence the probability of Cycle 1 is about that is two times more than that of the other cycle. But when , the net flux of Cycle 1 nearly vanishes, while Cycle 2 becomes globally attractive with probability , i.e. the forward transition time from Cycle to Cycle is exponentially small compared to the backward one. It is analogous to the phase transition in statistical mechanics, and the line of serves as the critical surface, which separates phase I (where Cycle 2 is global attractive) and phase II (where there is a positive but finite ratio between the two cycle fluxes) (Fig. 2 (a)).
Similar behavior can be observed in the stochastic model approaching the deterministic counterpart with , as and are sufficiently large. Setting and , we can show that as long as , the Cycle 2 is globally attractive; while as , the probability of the fixed point is about (Fig. 2 (b)), given . Analogous to phase transition, the two phases are separated by the critical surface .
The phase transition in such a Boolean network model can be rigorously and quantitatively explained by the mathematical theory of exponentially perturbed Markov chains (19); (20). Within this type of stochastic models, we can properly define the concepts of nonequilibrium “activation energy barrier” along each transiting trajectory from one attractor, either fixed point or limit cycle, to the other one. Then the theorem in (19) tells that the mean transiting time from the -th attractor to the -th one, satisfies
where is the minimum activation energy barriers along all the transiting trajectories. See (20) for the expression of .
Fig. 2(c), (d) present the minimum activation energy barrier difference with and respectively. Comparing these figures with Fig. 2(a), (b), we can easily conclude the direct relationship between minimum activation energy barrier difference and the limiting cycle fluxes: means the dominance of Cycle 2 (Phase I), implying tends to positive infinity with an exponentially increasing rate as ; while indicates a nontrivial probability ratio between the two cycles (Phase II). In this example, we do not consider the case , but it is definitely possible with a different matrix (See below).
Corresponding to the case , the minimum activation energy barriers and , even independent of . Hence phase I in which the Cycle 2 is globally attractive, occurs if and only if . On the other hand, corresponding to the case with , the minimum energy barriers and , therefore phase I occurs in the region of . The theory clearly explains the origin of the two phases in such stochastic Boolean network models.
The theorem above cannot give the quantitative ratio of cycle fluxes in phase II, as the minimum activation energy barriers forward and backward are equal. However, according to the theory in (19), we have already known that as the noise is sufficiently low, the transiting events between the two attractors almost surely concentrate on the optimal transition paths, namely the transiting trajectories with the minimum activation energy barrier.
In such specific stochastic Boolean network models, we can prove something more: the probability of taking each optimal transition path is exactly the same, hence the ratio of mean transiting time satisfies (20)
where denotes the number of states in the -th attractors of the deterministic models and represents the corresponding number of optimal transition paths.
The numbers of optimal transition paths and the corresponding net fluxes of the two cycles in Phase II are listed in Table 3, which perfectly agree with the numerical values obtained from Fig. 2 (a)(b).
Our theory can also determine how the transition time and the number of optimal transition paths between parallel pathways depend on each interaction’s strength, as well as can identify those possibly more crucial interactions in the biological networks. Because there are multiple optimal transition paths between the parallel pathways, small changes of most single nonzero components in the interaction matrix can not significantly change the behaviors of the network. For example, if we consider the probability of the normal steady state and normal pathway of cells before or after being perturbed, which is crucial for the maintaining of the cellular functions, none of the interactions in dramatically affect the qualitative behaviors of the network in the global attractive region (Phase I) for both cases of and as well as the non-global attractive region (Phase II) for the case , as long as the strength of the interaction only varies within a small fraction (Fig. 3 (a))(20). However, in the non-global attractive region (Phase II) when , a slightly varying of single nonzero component in can result in a qualitative change of the relative stability: from Phase II transiting to a new phase (Phase III) in which the normal steady state is globally attractive, which is shown in Fig. 3 (b). It indicates that decreasing the interaction strengths from p53 to Siah-1, Mdm-2 to p53 and -catenin to p19/14ARF or increasing the interaction strength from Siah-1 to -catenin can possibly make the normal cell more robustly against stochastic perturbations. It results from the fact that some interactions only affect the transition from one attractor to another but not on the opposite direction. Fig. 3 (c), (d) show the value of minimum activation energy barriers along both directions with respect to several nonzero elements of under the conditions of Fig. 3 (a) and (b) respectively. Both of the forward and backward minimum activation energy barriers increase with each interaction strength but possibly having different slopes. Also, even the topology of the networks can be significantly changed due to the modification of certain interaction strengths, such as the and (Fig. 1 (b)). See (20) for more details.
As a conclusion, parallel pathways often exist in biological systems, and their relative stability can be of great importance for implementing the correct cellular functions. These properties can not be understood and analyzed through simply observing the topology of the network or even studying deterministic models. We demonstrate here through stochastic Boolean network models of Siah-1/beta-catenin/p14/19 ARF loop that, there can possibly be various phases with qualitatively different behaviors of relative stabilities in such a nonequilibrium system. It is closely related to the nonequilibrium activation energy barriers along the transiting trajectories between these pathways, as well as the number of optimal transition paths with the minimum activation energy barrier. It can also help to identify which interactions in the network can dramatically affect the behaviors of the network, either making it better or worse. Our theory indicates that stochastic models, combined with the idea of nonequilibrium statistical mechanics, can help to address the stochastic robustness and relative stability in biological networks, at least semi-quantitatively. More implication of our theory remains to be further elucidated, by performing suitable experiments as well as by modeling of the other important biological systems.
We thank Das Biswajit, Xiao Jin and Chen Jia for carefully reading the manuscript. This research is financially supported by NSFC (Nos. 10901040, 21373021) and the Foundation for the Author of National Excellent Doctoral Dissertation of China (No. 201119), to which we would like to express our genuine gratitude.
- Stelling J, Sauer U, Szallasi Z, Doyle F J 3rd and Doyle J. Cell, 2004, 118(6): 675-685
- Li H, Helling R, Tang C and Wingreen N. Science, 1996, 273: 666-669
- Zeeman E C. Nonlinearity, 1988, 1: 115-155
- Zhu X-M, Yin L, Hood L and Ao P. Funct. Integr. Genomics, 2004, 4: 188-195; Zhou J X, Aliyu M D S, Aurell E and Huang S. J. R. Soc. Interface, 2012, 9: 3539-3553; Lu M Y, Onuchic J and Ben-Jacob E. Phys. Rev. Lett., 2014, 113: 078102; Walczak A M, Onuchic J N and Wolynes P G. Proc. Natl. Acad. Sci., 2005, 102: 18926-18931; Dykman M I, Mori E, Ross J and Hunt P M. J. Chem. Phys., 1994, 100: 5735; Assaf M, Roberts E and Luthey-Schulten Z. Phys. Rev. Lett., 2011, 106: 248102; Lv C, Li X, Li F and Li T. Plos One, 2014, 9(2): e88167; Ge H and Qian H. Phys. Rev. Lett. 2009, 103, 148103; Ge H, Qian H and Xie X S. Phys. Rev. Lett. 2015, 114, 078101; Lei X, Tian W, Zhu HY, Chen TQ and Ao P. Scientific Reports 2015, 5, 13597
- Huang R L, Wallqvist A and Covell D G. Mol. Cancer Ther., 2006, 5(9): 2417-2427
- Harris S L and Levine A J. Oncogene, 2005, 24: 2899-2908
- Kauffman S A. The Origins of Order: Self-Organization and Selection in Evolution. New York: Oxford University Press, 1993
- Kauffman S A. Sci. Am., 1991, 265: 78-84
- Akutsu T, Miyano S and Kuhara S. Identification of genetic networks from a small number of gene expression patterns under the Boolean network model. In: Altman R B, et al. eds. Proceedings of the Pacific Symposium on Biocomputing 99. Singapore, 1999, pp 17-28
- Hopfield J J. Proc. Natl. Acad. Sci., 1982, 79: 2554-2558
- Hopfield J J. Proc. Natl. Acad. Sci., 1984, 81: 3088-3092
- Amit D J. Modeling Brain Function: The World Of Attractor Neural Networks. Cambridge University Press, 1989
- Li F, Long T, Lu Y, Ouyang Q, Tang C. Proc. Natl. Acad. Sci., 2004, 101(14): 4781-4786
- Zhang Y, Qian M P, Ouyang Q, Deng M, Li F, Tang C. Physica D, 2006, 219: 35-39
- Ge H, Qian H and Qian M. Math. Biosci., 2008, 211: 132-152
- Von Dassow G, Meir E, Munro E M and Odell G M. Nature, 2000, 406: 189-192
- Ge H. and Qian M. J. Comput. Biol. 2009, 16: 119-132
- Jiang D-Q, Qian M and Qian M-P. Mathematical Theory of Nonequilibrium Steady States: On the Frontier of Probability and Dynamical Systems (Lecture Notes in Mathematics). Berlin: Springer, 2004
- Chen D, Feng J and Qian M. Science in China, 1994, 39: 7-28
- Supplementary Material, including Mathematical details.