Stability Analysis for the Cellular Signaling Systems Composed of Two Phosphorylation-Dephosphorylation Cyclic Reactions

The regulatory mechanisms in cellular signaling systems have been studied intensively from the viewpoint that the malfunction of the regulation is thought to be one of the substantial causes of cancer formation. On the other hand, it is rather difficult to develop the theoretical framework for investigation of the regulatory mechanisms due to their complexity and nonlinearity. In this study, more general approach is proposed for elucidation of characteristics of the stability in cellular signaling systems by construction of mathematical models for a class of cellular signaling systems and stability analysis of the models over variation of the network architectures and the parameter values. The model system is formulated as regulatory network in which every node represents a phosphorylation-dephosphorylation cyclic reaction for respective constituent enzyme. The analysis is performed for all variations of the regulatory networks comprised of two nodes with multiple feedback regulation loops. It is revealed from the analysis that the regulatory networks become mono-stable, bi-stable, tri-stable, or oscillatory and that the negative mutual feedback or positive mutual feedback is favorable for multi-stability, which is augmented by a negatively regulated node with a positive auto-regulation. Furthermore, the multi-stability or the oscillation is more likely to emerge in the case of low value of the Michaelis constant than in the case of high value, implying that the condition of higher saturation levels induces stronger nonlinearity in the networks. The analysis for the parameter regions yielding the multi-stability and the oscillation clarified that the stronger regulation shifts the systems toward multi-stability. How to cite this paper: Sueyoshi, C. and Naka, T. (2017) Stability Analysis for the Cellular Signaling Systems Composed of Two Phosphoryltion-Dephosphorylation Cyclic Re-actions. Computational Molecular Bioscience, 7, 33-45. https://doi.org/10.4236/cmb.2017.73003 Received: June 30, 2017 Accepted: August 13, 2017 Published: August 16, 2017 Copyright © 2017 by authors and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/


Introduction
Cellular signaling systems have been studied extensively from the recent viewpoint that their disorder is thought to be one of the causes for cancer formation since the systems are known to regulate biochemical reactions operating in cells for various functions such as cell differentiation, cell proliferation, and homeostasis.Figure 1(a) displays MAPK (Mitogen-activated Protein Kinase) cascade which is one of the typical cellular signaling systems, and has been studied intensively, elucidating the various regulatory characteristics for activities of living cells such as the switch-like responses, bi-stability, oscillations, robustness, and so on [1]- [6].The cellular signaling systems are comprised of enzymatic reactions such as phosphorylation-dephosphorylation cyclic reactions which are primary components of MAPK cascade as seen in Figure 1(a).Cyclic reaction node i; i P : active form of enzyme at node i; j P : activating enzyme for i U ; k P : inacti- vating enzyme for i P .The rate constant for the activation and the inactivation are de- noted by i a and i d , respectively.Computational Molecular Bioscience seems to be primary reaction system for other cellular signaling systems such as Rac1, PAK, and RhoA signaling networks [7].
Many studies have employed simulation analysis with the given values for the parameters in the signaling systems because of the difficulty in developing the analytical method for the systems analysis due to their nonlinearity.In this study, more general approach is proposed for characterization of the stability in cellular signaling systems by construction of mathematical models for a class of cellular signaling systems and simulation analysis of the models over variation of the network architectures and the values of parameters.The model system is formulated as regulatory network in which every node represents an activation-inactivation cyclic reaction for respective constituent enzyme of the network and the regulatory interactions between the activated enzyme and the reaction are depicted by arcs between nodes.The Michaelis-Menten mechanism is assumed for the reaction paths in each cyclic reaction and the emergence of the stable point in steady states of the network is analyzed.
It is biologically significant to analyze the characteristics of the stability in cellular signaling systems since stable points are convergent states of the relaxation process in dynamic changes due to random noises, and seem to correspond to the distinct biochemical states such as normal states or malfunctional states.
In this study the stability analysis is performed for all variations of the regulatory networks comprised of two nodes with multiple feedback regulation loops.

Formulation of Regulatory Networks for Cellular Signaling Systems
The regulatory network is formulated as follows.Every node represents an activation-inactivation cyclic reaction for respective constituent enzyme of the network.
The activated enzyme in a node acts on another node for the positive regulation which increases the active enzyme or for the negative regulation which increases the inactive enzyme through the reverse path.The regulatory interactions between the activated enzyme and the reaction are depicted by arcs between nodes.
The MAPK cascade is comprised of several cyclic reactions and has a feedback regulation from MAPK-PP as shown in Figure 1(a).Figure 1(d) shows a cyclic reaction in node i which is regulated positively by node j, and negatively by node k.We assume the Michaelis-Menten mechanism as the reaction mechanism in the cyclic reactions and the reaction rate equations are formulated as Equations ( 1)-(5): , , ( ) It should be noted that the Michaelis-Menten equation is adopted as the reaction mechanism, and therefore, the enzyme-substrate complex does not appear in the reaction rate equations.i T is the total concentration of enzyme i P in node i and i R is the relative concentration of i P .j L and k L are the norma- lized Michaelis constants of enzymes j P and k P , respectively.In addition to j L and k L , the normalized rate Equation (4) has a parameter i K , which is the ratio of maximum inactivation velocity to maximum activation velocity.The value of j R is set to be unity in the case of no positive regulation for node i, whereas the value of k R is set to be unity in the case of no negative regulation.
We analyze the characteristics of stability for all variations of two-node regulatory networks with multiple feedback regulation loops under the restriction that each node has at most one positive regulation and one negative regulation as designated in Figure 2.

Characteristics of Multi-Stability in the Regulatory Networks
In Figure 3 the effects of Michaelis constant on the emergence ratio of multistable state are demonstrated for each regulatory network examined.It follows from this figure that the networks are divided into two groups; that is, the networks with high emergence ratio (F, G, H, I, J) and those with zero or slight emergence ratio (A, B, C, D, E).The former group is further categorized into those with very high emergence ratio (I, J) and those with moderate emergence ratio (F, G, H).The latter group also is further assorted into the networks with zero emergence ratio (A, B, C) and slight emergence ratio (D, E).It is observed that the emergence ratio decreases when the Michaelis constant gets higher, and then it converges at certain levels for the networks with very high emergence ratio.On the other hand, multi-stable state arises when the Michaelis constant is low, and then the ratios become zero for the networks with the moderate emergence ratio.The highest emergence ratio at Michaelis constant of 5 2 − is 100%, 69.4%, 25.6%, 19.0%, and 16.5% for the networks J, I, H, G, and F, respectively.
Regarding the networks E and D, the highest ratio of 1.7% is attained at Michaelis constant of 5 2 − and 1 2 − , respectively.In Figure 4 the distributions of multi-stable equilibrium points in the parameter space of 1 K and 2 K are illustrated for each regulatory network and various Michaelis constants examined.The heat maps in Figure 4 mean the sensitivities which quantify the effect of the change of parameter values on the set of stable equilibrium points.The precise definition of the sensitivity is provided in the method section.The blue area corresponds to the zero or slight sensitivities, implying that the set of stable equilibrium points hardly change in the area.In contrast, the red area indicates high sensitivity where the set changes drastically.
It should be noted that, as the networks G, H, and J are symmetric for each node, the graphs for these networks are symmetric to the diagonal line.It is revealed that bi-stable states appear in the area with low values of i K for the network G which is composed of mutual negative feedbacks.In contrast, bi-stable states are likely to appear in the area with high values of i K for the , R R -plane in subgraphs in the right four graphs.Computational Molecular Bioscience networks H and J which contain mutual positive feedbacks.Since the parameter i K is the maximum velocity of inactivation divided by that of activation, higher and lower values of i K enforce negative regulations and positive regulations, respectively.Therefore, bi-stable state is likely to appear in the area where the regulations are enforced.
The bi-stable area for lower value of the Michaelis constant include the area for higher value of the Michaelis constant for the network F, G, H, I, and J, that is consistent with the tendency that lower value of the Michaelis constant is favorable to emergence of multi-stable state shown in Figure 3.Likewise, tri-stable states in the area with low values of i K of network J gets smaller with lower value of the Michaelis constant, and then vanishes at L = 2 0 .
It follows from Figure 4 that the parameter space is divided into the robust several subspaces according to the change in i K and bifurcation occurs on the high sensitive purple area across the low sensitive light blue area.The set of stable equilibrium points in the light blue areas contains the combinations of nearly zero or nearly unity values of 1 R or 2 R .The graphs in the bottom row explain these combinations for each graph in the upper rows.The values under 0.1 and the values over 0.9 are regarded as zero and unity, respectively, in these graphs.
It is noted that the networks with mutual positive feedbacks yield the bi-stable state in which the both of two nodes are fully activated or fully inactivated, while the networks with mutual negative feedbacks yield the bi-stable state in which either of the nodes is fully activated and the other node is fully inactivated in the parameter space with high or moderate values of i K .On the other hand, both of the two nodes are fully activated in the areas with low values of i K for all the networks examined.This feature is predictable since the low value of i K en- forces the activation, and it seems capable to yield the tri-stable equilibrium points in network J.
The white area in the graph at for the network H indicates that the set of stable equilibrium points is empty since the eigenvalues of this point are −2 and 0. Likewise, the white areas in the graphs for the network D depict the empty set of stable equilibrium points because of the oscillations.

Characteristics of Oscillations in the Regulatory Networks
Figure 5 displays the limit cycles for the oscillations indicated by green cross marks in the graphs for network D in Figure 4.It is observed that the oscillation area for lower value of the Michaelis constant includes the area for higher value of the Michaelis constant as in the case of multi-stability.Furthermore, oscillation arises in the case that 2 K is less than unity.It seems difficult to find out the simple rules for the locations and the sizes of limit cycles with respect to the value of Michaelis constant L; however, the size of the limit cycle for each L value is getting bigger with higher 2 K , and then getting smaller, suggesting that the optimal value of 2 K seems to exist to maximize the size of the limit cycle for each L value.It is also seen that the almost of all limit cycles appear in the area with small values of ( ) , R R , and go through the neighbor points of the origin ( ) It is found that all of the limit cycles are counterclockwise.The limit cycles arise in the case of 1 2 K K > , implying that the inactivation reaction is more dominant than the activation reaction at node 1, while the activation reaction is more dominant than the inactivation reaction at node 2. As for the network D, node 1 has a positive auto-regulation and a negative regulation from node 2, while node 2 has a positive regulation from node 1.Therefore, if the system is operated around ( ) ( ) , 0.5, 0.5 R R = , the value of 1 R would decrease under the condition of 1 2

K K >
; that is, the state would move left in ( ) When the value of 1 R becomes small enough, then the value of 2 R is getting lower since the activating power of 1 R to 2 R is getting weak; that is, the state would move down along with left edge of the ( ) becomes small enough, and then the value of 1 R would be getting higher since the inactivating power of 2 R to 1 R is getting weak; that is, the state would move right along with bottom edge of the ( ) , R R -plane.When the value of 1 R becomes high enough, the value of 2 R is getting higher since the activating power of 1 R to 2 R is retrieved; that is, the state would move up, and return to the starting point of the state.Subsequently, the next cycle would start again.
This might be the mechanism for yielding the limit cycles.

Discussion
The characteristics of the emergence of stable equilibrium points are analyzed quantitatively by the emergence ratio and the sensitivity of equilibria which are newly proposed and defined in this study.These quantities seem to successfully detect the effects of the architectures and the values of parameters on the characteristics.
Comparison of the emergence ratios of bi-stable state for the regulatory networks as shown in Figure 2 suggests that the networks with high ratios have either the mutual positive or negative feedbacks.The auto-regulation seems to further affect the emergence of multi-stability.The effects of the mutual regulations and the auto regulations could be summarized as in Table 1; that is, positive mutual regulations or negative mutual regulations are likely to raise the emergence ratio of multi-stable state, while positive and negative mutual regulations tend to make the emergence ratio lower.With respect to the auto regulations, containment of node with both of positive auto regulation and negative regulation from Table 1.Effects of network structure on emergence of multi-stable state.

Boost Decline
The partial network structures in the upper row make the emergence ratio of multi-stable state higher, while those in the bottom row work reversely.Computational Molecular Bioscience the other node could boost multi-stable systems, while a network containing node with both of negative auto regulation and positive regulation from the other node curtail the multi-stable systems.
As a common feature for the regulatory networks examined, smaller normalized Michaelis constant L yields higher emergence ratio of multi-stability or oscillation, implying that the condition of the higher saturation levels induces stronger nonlinearity.In the parameter spaces multi-stability arises in the region where i K is small in the case that the node is regulated positively and i K is large in the case that the node is regulated negatively, suggesting that the stronger regulation shifts the systems toward multi-stability.Multi-stable state is observed to emerge very robustly on the regulatory network J which undergoes mutual negative feedbacks and positive auto regulation on both of two nodes.It has been suggested that this robust reaction system works as the sequential Bayesian filter to reduce the noise in cells [16] [17].
It is difficult to visualize and analyze the number of the stable equilibrium points and the values of the points, since high dimensional representation is required due to the existence of plural stable equilibrium points in two-dimensional space such as ( ) , R R -plane at each point in the parameter space such as ( ) , K K -plane.The sensitivity of equilibria is introduced to resolve this difficulty in this study.At first, the sensitivity is utilized to divide the parameter space to the robust area and high sensitive area, then the number of the stable equilibrium points and the values of the points are analyzed in the robust areas.
It is revealed that the variations of these points are scarce and sorted in the robust areas as shown in the graphs in the bottom row in Figure 4.
Furthermore, the emergence of oscillations are analyzed to elucidate that the emergence of the oscillations are rare from the viewpoints of the regulatory structures and the values of parameters comparing to the high emergence of multi-stable states, such as 100% for the network J at the small values of L.
It is suggested that the method proposed in this study utilizing the emergence ratio and the sensitivity of equilibria are useful for this kind of analysis.On the other hand, it is unclear if these results are available to the higher order networks such as three-node or four-node regulatory networks due to the nonlinearity of the systems.Therefore, the similar analysis for the higher order networks are undergoing in our laboratory, facing the difficulty for the high dimension to visualize.It is also required to invent the new way to visualize and analyze the high dimensional dynamics.

Range of Parameter Space for the Analysis
The analysis is performed with variation of the parameter values in appropriate range to cover the assumed values for reactions of MAPK cascade described in other studies [18] [19] [20] [21] [22].From physiological point of view, the parameter set serves as surrogates of varying cellular activities across different cell i j i j i j i j i j i j i j i j

S i j D e e D e e D e e D e e
where , i j e is a set of the stable equilibrium points at 1 2 i K = and 2 2 j K = in the parameter space.Namely, the sensitivity is the mean value of distances between a set of the stable equilibrium points and a set of the stable equilibrium points at Neumann neighborhood.The system with low sensitivity is robust to the change of parameter values, while the high sensitivity might imply that the stability characteristics of the system could be regulated by the parameters.

Figure 1 .
Figure 1.Regulatory network for MAPK cascade.(a): The MAPK cascade which is composed of several cyclic reactions with a feedback regulation from MAPK-PP; (b): Simplified version of the MAPK cascade. 1 P , 2 P , 3 P , and 4 P correspond to active Ras, Raf, Mek, and Erk, respectively, while 1 U , 2 U , 3 U , and 4 U are inactive forms of those en- zymes; (c): Regulatory network representing the MAPK cascade.Each node represents the cyclic reaction in (b); (d): Cyclic reaction in node i. i U : inactive form of enzyme at

Figure 1 (
b) is a simplified representation of the MAPK cascade of Figure 1(a).

Figure 1 (
c) illustrates the regulatory network of the MAPK cascade, which is analyzed in this study.

Figure 2 .
Figure 2. Regulatory networks with feedback regulations.Red arrows and blue arrows depict the positive and negative regulations, respectively.The number in each node corresponds to the parameter 1 K or 2 K in the graphs of parameter space as shown in Figure 4 and Figure 5.

Figure 3 .
Figure 3. Emergence ratio of multi-stability for the regulatory networks.The abscissa depicts each regulatory network and the ordinate represents logarithm values of base 2 for the normalized Michaelis constant L. The applicate axis represents the emergence ratio.

Figure 4 . 2 − and 2 2
Figure 4. Distributions of multi-stable equilibrium points and the sensitivities for each regulatory network in the parameter space.Each column corresponds respectively to the individual regulatory networks of D, E, F, G, H, I, and J for which multi-stability emerges.L denotes the normalized Michaelis constant.The abscissa and the ordinate in each small graph depict the logarithm values of base 2 for 1 K and 2 K , respectively.The subscripts of K correspond to the numbers in the nodes in Figure 2. Orange circle marks and blue triangle marks locate the bi-stable points and tri-stable points, respectively.Green cross marks in the graphs for network D depict the oscillation points for the values between 5 2 − and 2 2 − .The sensitivities of a set of equilibrium points to the change of parameter values of 1 K and 2 K are displayed as the heat maps in which light blue shows low sensitivity and pur- ple indicates high sensitivity.The white areas in the graph for network H and

Figure 5 .
Figure 5.The characteristics of the limit cycles in Network D. The abscissa depicts the parameter values of 2 K and the ordinate denotes the values of the normalized Michaelis constant L. The parameter values of 1 K is 2 0 except the three graphs in second row and two graphs in forth row, where the values of 1 K is 2 1 as indicated in the respective rows.The abscissa in each small graph depicts 1 R which is the activation level of node 1, while the ordinate shows 2 R which is the activation level of node 2. The subscripts of K and R correspond to the numbers in the nodes in Figure 2. The light blue arrows demonstrate the dynamic flows and the series of red arrows indicate the courses of the limit cycles.