Unstable Modes and Order Parameters of Bistable Signaling Pathways at Saddle-Node Bifurcations : A Theoretical Study Based on Synergetics

Mathematical modeling has become an indispensable part of systems biology which is a discipline that has become increasingly popular in recent years. In this context, our understanding of bistable signaling pathways in terms of mathematical modeling is of particular importance because such bistable components perform crucial functions in living cells. Bistable signaling pathways can act as switches ormemory functions and can determine cell fate. In the present study, properties of mathematical models of bistable signaling pathways are examined from the perspective of synergetics, a theory of self-organization and pattern formation founded by Hermann Haken. At the heart of synergetics is the concept of so-called unstable modes or order parameters that determine the behavior of systems as a whole close to bifurcation points. How to determine these order parameters for bistable signaling pathways at saddle-node bifurcation points is shown. The procedure is outlined in general and an explicit example is worked out in detail.

While on a conceptual level there is a good understanding of bistable signaling pathways, there is still a need for improving our methodological tools to address such bistable components from a mechanistic point of view.Within the framework of ordinary differential equation (ODE) models, bistable signaling pathways have been described by single species or multiple species models.Single species models capture the pathway activity by a single variable.A benchmark model in this regard is a model for promoter activity under the impact of a positive feedback loop.The feedback loop is established by means of a transcription factor that acts as an activator [2,4,8,9].In the case of single variable models, the bistability is given in terms of a low and high state of the species concentration under consideration or the gene expression or promoter activity of interest.In what follows, the focus will be on two-species models as described by two coupled ODEs [2,3,5].The two states of the bistable element typically describe that either species  is dominant over species  or species  is dominant over species ; that is, there are two states of interest that will in what follows be denoted by (high), (low) and (low), (high).In order to function as a switch, the bistable component must feature a control parameter and the bistable component must be driven through a saddle-node bifurcation by varying the control parameter gradually.At the saddle-node bifurcation the bistable pathway becomes monostable such that only one of the two species can dominate the signaling pathway.In doing so, bistable models can a switch from (high), (low) to (low), (high) or vice versa from (low), (high) to (high), (low).While the saddle-node bifurcation has been studied using analytical and numerical solution methods [2,3,5,10], in this context little is known about the socalled unstable mode or order parameter [11] emerging at this saddle-node bifurcation.According to the theory of synergetics, a theory of pattern formation and self-organization, the unstable mode at the bifurcation point determines the order of the whole two-species system and for this reason 2 Advances in Mathematical Physics may be referred to as order parameter [11].In the following section, the general approach to determine the unstable mode or order parameter will be described and subsequently an application to the bistable model proposed by Markevich et al. (2004) [5] will be presented.

Unstable Mode and Order
Parameter of Bistable Signaling Pathways Involving Two Species where  and  are appropriately defined functions and depend on the pathway being considered.Equation ( 1) is assumed to involve a control parameter .Furthermore, ( 1) is assumed to be bistable in a particular interval of .At the boundaries of that interval a stable fixed point merges with an unstable fixed point giving rise to a saddle-node bifurcation.
The stable and unstable fixed points merge and disappear and the pathway becomes monostable.Let us consider the saddlenode bifurcation in more detail.To this end, we "follow" the stable fixed point coordinates (st) and (st) involved in the saddle-node bifurcation along variations of .The linearized version of ( Note that in (3)  1 comes with a plus sign in front of the square root function, whereas  2 exhibits the minus sign.In general, at a saddle-node bifurcation of a two-dimensional dynamical system there is one eigenvalue equal to zero while the other is negative.Using the aforementioned plus and minus sign convention of (3), at the saddle-node bifurcation, we have det() = 0, Tr() < 0,  1 = 0, and  2 < 0. From det() = 0 it follows that Moreover, substituting det() = 0 into (3) we obtain It is useful to summarize this intermediate result.Let (crit) denote a critical value at which the saddle-node bifurcation occurs.Consequently, at  = (crit) we have  1 = 0 and  2 < 0 and (5) holds.If we shift next  away from the critical value by a small amount such that the pathway is monostable,  = (crit)+, where  is small, then there are two directions in the two-dimensional state space.In one direction the dynamics evolves relatively slowly whereas in the direction orthogonal to that direction the dynamics evolves relatively fast because for  = 0 we have  1 = 0 as compared to  2 < 0.
In line with synergetics, the direction associated with the slow dynamics corresponds to the unstable mode at  = 0 and dominates the behavior of the whole system; that is, the mode can be considered as the order parameter [11].In particular, the mode dominates the signaling pathway even when  is not zero but a small parameter.This holds for perturbations that drive the system temporarily away from the fixed point (st), (st) in the bistable parameter domain and for the transient dynamics in the monostable domain that describes the transition from the disappearing fixed point (st), (st) to the alternative fixed point that is still stable in the monostable domain.
Let us determine the unstable mode.In general, that is, for arbitrary eigenvalues  1 and  2 , the eigenvectors normalized to unity are given in components by for  = 1, 2. At the saddle-node bifurcation point, we have  1 = 0 and the corresponding eigenvector represents the unstable mode or order parameter.Consequently, at a saddlenode bifurcation the order parameter of a bistable signaling pathway described by a two-species model given in terms of two coupled ODEs is defined by Advances in Mathematical Physics with , , , , and  > 0. The concentration levels are normalized such that () and () are restricted to the interval [0, 1].The model and generalizations of it have also been discussed in subsequent studies [10,12].

Nullcline Approach.
For sake of simplicity, we introduce the new parameters: The nullclines are given by Using the nullclines we would like to illustrate the characteristic features of the pathway model.After that, a more detailed analysis will be conducted.
In order to illustrate the characteristic features of the pathway model, we consider first the symmetric case given by  1 =  2 = .In this case, the model can exhibit a single asymptotically stable fixed point (Figure 1) or three fixed points (Figure 2).In the latter case, the pathway model is bistable.The fixed point on the diagonal is unstable.Note that in the symmetric case the number of fixed points only depends on a single parameter  as we will see later.In particular, at a critical value of , there is a pitchfork bifurcation connecting the monostable and bistable parameter domains.We are interested in the switch-like functioning of bistable pathway models and the order parameter that characterizes the switching from one state to another state.With respect to the model by Markevich et al. (2004) [5], we need to consider the asymmetric case.As we will see below, if  1 is much larger than  2 the model is monostable with the species  dominating the species .In contrast, if  2 is much larger than  1 then the model is monostable with the species  dominating the species .Consequently, if we start initially with a large ratio  2 / 1 then the pathway exhibits a high state for  and a low state for .If we gradually decrease the ratio, then the fixed point for (high) and (low) will converge towards the unstable fixed point of the pathway model.The two fixed points will merge and disappear in a saddle-node bifurcation.This is shown in Figures 3 and 4. The pathway will switch to the alternative state with a relatively high state for  and a relatively low state for .
Figures 1, 2, 3, and 4 illustrate the two basic mechanisms leading to the saddle-node bifurcation underlying the switchlike functioning of the model.First, in the symmetric case, the model is bistable.Second, when breaking the symmetry, in the asymmetric case, the bistable domain can be destroyed via a saddle-node bifurcation resulting in a switch from the state (low), (high) to the state (high), (low) or vice versa.Let us analyze these mechanisms in more detail.To this end, we derive a fixed point equation.

Cubic Fixed Point Equation and Bifurcation Diagram.
In order to determine the fixed point semianalytically, we put the left-hand sides of (9) equal to zero.From the two relations in (9) we can obtain the following relations: It is useful to introduce the new variables  and V like Moreover, we introduce the new parameters  and  like The parameter  measures the degree of asymmetry.In the symmetric case we have  = 1.From ( 14) we obtain the intermediate result: and finally In deriving (18) we used  > 0. The variable  as a function of V defined by ( 18) is a monotonically decaying function for any parameter .For  → 0 we have  = SQRT(1 − V 2 ).Importantly, from (17), it follows that for any parameter  the solutions of (18) are symmetric; that is, if  =  and V =  is a solution, then  =  and V =  is a solution as well.From ( 14) it follows that with Equation ( 19) is a cubic equation in the variable V. Therefore, there are maximal three solutions.If  = 1 then the equation is symmetric in  and V. Therefore, in view of the symmetry property of (18) for any parameter  and the symmetry property of (19) in the special case  = 1, we conclude that if  = 1 then the number of solutions is either 1 or 3. From (19) it follows that in the symmetric case  = 1 (which implies by  1 =  2 = ) the parameter  = 1/ 2 acts as control  9) acting as a switch for a particular value of  when varying the degree of asymmetry.The stationary solutions (st) are plotted versus the control parameter values .The solutions (st) were computed from (18), ( 19), (20), and (21).Solid line: asymptotically stable fixed points of the state (low), (high).Dotted line: unstable fixed points in the bistable parameter domain.Dashed line: asymptotically stable fixed points of the state (high), (low).Parameters:  = 0.01,  = 0.5, and  = 3.5.
For any given set of parameters , , and , the stable and unstable fixed points (st) and (st) of the signaling pathway model ( 9) can be determined numerically by finding the solutions (roots) of the cubic equation defined by (19) in combination with (18) and (20).In doing so, the bifurcation diagram can be constructed when defining a suitable control parameter .In order to focus on switch-like functioning we vary the degree of asymmetry like with  > .In particular, for  = ( − )/2 we have the symmetric case with  = 1 and  1 =  2 =  = ( + )/2.We choose the parameters  and  sufficiently large such that in this (symmetric) case the model is bistable.Figure 5 presents the bifurcation diagram thus obtained for a fixed value of .Figures 6, 7, and 8 depict the functions   and   of the implicit equation (19) for the symmetric case (Figure 7) and the critical values at which saddle-node bifurcations occur (Figures 6 and 8).19).The functions   (solid line) and   (dashed line) of ( 19) are plotted versus V at the saddle-node bifurcation point in which the states (high) and (low) emerge or disappear.Parameters:  = 0.8; all other parameters are as in Figure 5.   19).The functions   (solid line) and   (dashed line) of ( 19) are plotted versus V at the saddle-node bifurcation point in which the state (low), (high) emerges or disappears.Parameters:  = 2.2; all other parameters are as in Figure 5.  3), (4), and (22) for the state (low), (high) (solid lines) and the state (high), (low) (dashed line) shown in Figure 5 as functions of .Magnifying (not shown) the top parts of the two panels reveals that  1 converges to zero when approaching the saddle-node bifurcation points, whereas  2 assumes a finite negative value.Parameters:  =  = 1; all other parameters are as in Figure 5.
We computed the matrix coefficients for the asymptotically stable fixed points (low), (high) and (high), (low) shown in Figure 5 as functions of the control parameter .Using (3) we computed the eigenvalues.Figure 9 presents the eigenvalues as functions of .We found that both eigenvalues were real-valued in the considered parameter domain.Importantly, as expected,  1 converged to zero when  4) and ( 22) for the state (low), (high) (solid line) and the state (high), (low) (dashed line) shown in Figure 5 and as function of .Parameters:  =  = 1; all other parameters are as in Figure 5.
approaching the saddle-node bifurcation points, whereas  2 assumed a finite negative value.The convergence of one of the eigenvalues to zero can be seen even more clearly when studying the determinant of the Jacobian matrix as function of .We computed det() from ( 4).The result is shown in Figure 10.As expected, det() converged to zero at the saddlenode bifurcation points.Finally, we computed the eigenvector associated with  1 from (7) for the asymptotically stable fixed points (low), (high) and (high), (low) shown in Figure 5 and as functions of the control parameter . Figure 11 presents the eigenvector components thus obtained.The order parameters of the signaling pathway model correspond to the eigenvectors at the two saddle-node bifurcation points.These order parameters (order parameter vectors) are highlighted in Figure 11 by circles.Figure 11 demonstrates explicitly how the general methodology outlined in Section 2.1 can be used to determine order parameters of bistable signaling pathways at saddle-node bifurcation points.

Discussion
Bistable signaling pathways that function as switches and can be described by means of two-variable ODE models were considered.It was demonstrated that in general the switching behavior involves so-called order parameters that identify directions in which the dynamics evolves relatively slowly.As it is known from synergetics [11], a theory of self-organization and pattern formation, these directions play a special role because they dominate in a particular sense the dynamics of the whole system.For a model of a bistable signaling pathway proposed by Markevich et al. (2004) [5], the order parameters at the two switching points described by the pathway model were determined explicitly.
In previous studies, much effort has been devoted to determine the bifurcation diagrams of bistable multispecies  3) and ( 4) are shown for the state (low), (high) (solid line) and the state (high), (low) (dashed line) reported in Figure 5 and as functions of .At the saddle-node bifurcation points (circles), these components constitute the respective order parameters.
signaling pathways either numerically or by means of analytical methods [3, 5-7, 10, 12].The present study adds a novel perspective to these studies and to similar studies.Accordingly, in the mathematical space spanned by the biochemical species of bistable signaling pathways exhibiting switch-like functioning there exist particular directions that determine (in the sense of synergetics) the emerging order when a switch occurs.For signaling pathways described by means of two-variable models these directions can be determined in terms of order parameter vectors (or unstable modes) as outlined in Section 2.1.
The method presented in Section 2.1 in principle can be generalized and applied to signaling pathways characterized by more than two biochemical species.For such higher dimensional problems, it might be difficult to determine analytical expressions for eigenvalues and eigenvectors.If so, it might be still possible to carry out the procedure sketched in Section 2.1 by means of numerical solution methods.Future work might be devoted to address this issue.

Figure 3 :Figure 4 :
Figure 3: Illustration of the saddle-node bifurcation (part I) for the asymmetric pathway model.The stable fixed point given by (high) and (low) is about to merge with the unstable fixed point.Parameters:  = 0.01,  1 = 2.2, and  2 = 1.8.

Figure 5 :
Figure 5: Bifurcation diagram of the bistable pathway model (9) acting as a switch for a particular value of  when varying the degree of asymmetry.The stationary solutions (st) are plotted versus the control parameter values .The solutions (st) were computed from (18), (19), (20), and (21).Solid line: asymptotically stable fixed points of the state (low), (high).Dotted line: unstable fixed points in the bistable parameter domain.Dashed line: asymptotically stable fixed points of the state (high), (low).Parameters:  = 0.01,  = 0.5, and  = 3.5.

2. 2 . 4 .Figure 6 :
Figure6: Illustration of the saddle-node bifurcation using the cubic fixed point equation (19).The functions   (solid line) and   (dashed line) of (19) are plotted versus V at the saddle-node bifurcation point in which the states (high) and (low) emerge or disappear.Parameters:  = 0.8; all other parameters are as in Figure5.

Figure 9 :
Figure 9: Eigenvalues  1 (a) and  2 (b) computed from (3), (4), and (22) for the state (low), (high) (solid lines) and the state (high), (low) (dashed line) shown in Figure5as functions of .Magnifying (not shown) the top parts of the two panels reveals that  1 converges to zero when approaching the saddle-node bifurcation points, whereas  2 assumes a finite negative value.Parameters:  =  = 1; all other parameters are as in Figure5.

Figure 10 :
Figure 10: Determinant of the Jacobian matrix  as computed from (4) and (22) for the state (low), (high) (solid line) and the state (high), (low) (dashed line) shown in Figure5and as function of .Parameters:  =  = 1; all other parameters are as in Figure5.

Figure 11 :
Figure 11: Eigenvector components of the first eigenvector and order parameters at the saddle-node bifurcation points.The eigenvector components V 1, (a) and V 1, (b) as computed from (3) and (4) are shown for the state (low), (high) (solid line) and the state (high), (low) (dashed line) reported in Figure5and as functions of .At the saddle-node bifurcation points (circles), these components constitute the respective order parameters.
Let us first consider the general case of a two-species model of a bistable signaling pathway given in terms of two coupled ODEs.Accordingly, the pathway is described by the concentrations  and  of the two species under consideration measured in arbitrary units.Variations () and () over time  are assumed to satisfy 2.1.General Case.