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

Authors: Till D. Frank
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

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 or memory 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.


---

## Body

## 1. Introduction

Bistable signaling pathways on the cellular level play a key role [1–7]. They can act as switches in general and determine cell fate and function as memory storage in particular [1–7]. 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 A is dominant over species B or species B is dominant over species A; that is, there are two states of interest that will in what follows be denoted by A(high), B(low) and A(low), B(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 A(high), B(low) to A(low), B(high) or vice versa from A(low), B(high) to A(high), B(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 so-called 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 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.


## 2. Unstable Mode and Order Parameter of Bistable Signaling Pathways Involving Two Species

### 2.1. General Case

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 concentrationsx and y of the two species under consideration measured in arbitrary units. Variations x(t) and y(t) over time t are assumed to satisfy(1)ddtx=fx,y,ddty=gx,y,where f and g 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 saddle-node bifurcation in more detail. To this end, we “follow” the stable fixed point coordinates x(st) and y(st) involved in the saddle-node bifurcation along variations of α. The linearized version of (1) at that stable fixed point reads(2)ddtδ=a11δ+a12ε,ddtε=a21δ+a22εwith x=x(st)+δ and y=y(st)+ε. The variables ε and δ denote small deviations from the fixed point x(st), y(st). The matrix coefficients ajk correspond to the coefficients of the Jacobian matrix J at x(st) and y(st). The eigenvalues λ1 and λ2 of the matrix J are given by(3)λ1,2=TrJ2±TrJ24-det⁡J,where the trace Tr(J) and the determinant det⁡(J) of J are given by(4)TrJ=a11+a22,det⁡J=a11a22-a12a21.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⁡(J)=0, Tr(J)<0, λ1=0, and λ2<0. From det⁡(J)=0 it follows that(5)a11a22=a12a21.Moreover, substituting det⁡(J)=0 into (3) we obtain(6)λ2=TrJ.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 x(st), y(st) in the bistable parameter domain and for the transient dynamics in the monostable domain that describes the transition from the disappearing fixed point x(st), y(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(7)vj,x=11+a11-λj/a122,vj,y=1-vj,x2for j=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 saddle-node 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(8)v1,x=11+a11/a122,v1,y=1-v1,x2.


### 2.2. The Signaling Pathway Model by Markevich et al. (2004)

#### 2.2.1. Definition of the Model

Markevich et al. (2004) [5] discussed the following two-species signaling pathway model:(9)ddtx=fx,y=A1-x-yK+1-x-BxK+1-y,ddty=gx,y=C1-x-yK+1-y-DxK+1-xwith A, B, C, D, and K>0. The concentration levels are normalized such that x(t) and y(t) are restricted to the interval [0,1]. The model and generalizations of it have also been discussed in subsequent studies [10, 12].


#### 2.2.2. Nullcline Approach

For sake of simplicity, we introduce the new parameters:(10)r1=AB,r2=CD,K∗=K+1.The nullclines are given by(11)x=-h1y2±h12y4-r1y2-1+K∗y+K∗,(12)y=-h2x2±h22x4-r2x2-1+K∗x+K∗with h1(y)=(r2y-K∗(1+r2))/2 and h2(x)=(r1x-K∗(1+r1))/2. 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 byr1=r2=r. 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 r as we will see later. In particular, at a critical value of r, 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 r1 is much larger than r2 the model is monostable with the species x dominating the species y. In contrast, if r2 is much larger than r1 then the model is monostable with the species y dominating the species x. Consequently, if we start initially with a large ratio r2/r1 then the pathway exhibits a high state for y and a low state for x. If we gradually decrease the ratio, then the fixed point for y(high) and x(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 x and a relatively low state for y.Figure 1
Nullclines as defined by (11) (dashed line) and (12) (solid line) for the monostable case. Parameters: K=0.01 (as in [10]) and r=0.7. All parameters are in arbitrary units.Figure 2
It is as in Figure1 but for the bistable case. Parameters: K=0.01 and r=1.5.Figure 3
Illustration of the saddle-node bifurcation (part I) for the asymmetric pathway model. The stable fixed point given byy(high) and x(low) is about to merge with the unstable fixed point. Parameters: K=0.01, r1=2.2, and r2=1.8.Figure 4
Illustration of the saddle-node bifurcation (part II). The pathway has become monostable. Magnification (not shown) of the top left corner of the phase space reveals that the nullclines do not intersect each other in this area. This is because the asymmetry of the pathway model has increased as compared to Figure3. Parameters: K=0.01, r1=2.8, and r2=1.2.Figures1, 2, 3, and 4 illustrate the two basic mechanisms leading to the saddle-node bifurcation underlying the switch-like 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 x(low), y(high) to the state x(high), y(low) or vice versa. Let us analyze these mechanisms in more detail. To this end, we derive a fixed point equation.


#### 2.2.3. 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:(13)ABK∗-yK∗-xx=11-x-y=CDK∗-xK∗-yy⟹r1K∗-yK∗-xx=r2K∗-xK∗-yy,(14)AB1-x-yx=K∗-xK∗-y=DCy1-x-y⟹r11-x-yx=1r2y1-x-y.It is useful to introduce the new variables u and v like(15)x=u2,y=v2⇕u=x,v=y.Moreover, we introduce the new parameters ξ and s like(16)ξ=BDAC=1r1r2>0,s=BCAD=r2r1>0.The parameter s measures the degree of asymmetry. In the symmetric case we have s=1. From (14) we obtain the intermediate result:(17)1-u2-v2=ξuvand finally(18)uv;ξ=-12ξv+1+ξ4-1v2.In deriving (18) we used u>0. The variable u as a function of v defined by (18) is a monotonically decaying function for any parameter ξ. For ξ→0 we have u=SQRT(1-v2). Importantly, from (17), it follows that for any parameter ξ the solutions of (18) are symmetric; that is, if u=a and v=b is a solution, then u=b and v=a is a solution as well. From (14) it follows that(19)fLv=fRuv;ξ;swith(20)fLv=K∗-v2v,fRuv;ξ;s=sK∗-u2u.Equation (19) is a cubic equation in the variable v. Therefore, there are maximal three solutions. If s=1 then the equation is symmetric in u and v. Therefore, in view of the symmetry property of (18) for any parameter ξ and the symmetry property of (19) in the special case s=1, we conclude that if s=1 then the number of solutions is either 1 or 3. From (19) it follows that in the symmetric case s=1 (which implies by r1=r2=r) the parameter ξ=1/r2 acts as control parameter and induces a pitchfork bifurcation. For example, for K=0.01 (K∗=1.01) and s=1 numerical evaluation of (19) shows (see below) that for r=2 (19) exhibits three admissible solutions (i.e., solutions in the interval between 0 and 1).For any given set of parametersξ,s, and K, the stable and unstable fixed points x(st) and y(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(21)r1=A+α,r2=B-α,α∈0,B-Awith B>A. In particular, for α=(B-A)/2 we have the symmetric case with s=1 and r1=r2=r=(A+B)/2. We choose the parameters A and B 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 K. Figures 6, 7, and 8 depict the functions fL and fR 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).Figure 5
Bifurcation diagram of the bistable pathway model (9) acting as a switch for a particular value of K when varying the degree of asymmetry. The stationary solutions x(st) are plotted versus the control parameter values α. The solutions x(st) were computed from (18), (19), (20), and (21). Solid line: asymptotically stable fixed points of the state x(low), y(high). Dotted line: unstable fixed points in the bistable parameter domain. Dashed line: asymptotically stable fixed points of the state x(high), y(low). Parameters: K=0.01, A=0.5, and B=3.5.Figure 6
Illustration of the saddle-node bifurcation using the cubic fixed point equation (19). The functions fL (solid line) and fR (dashed line) of (19) are plotted versus v at the saddle-node bifurcation point in which the states x(high) and y(low) emerge or disappear. Parameters: α=0.8; all other parameters are as in Figure 5.Figure 7
The functionsfL (solid line) and fR (dashed line) of (19) in the symmetric case. The signaling pathway model is bistable. Parameters: α=1.5 (r1=r2=r=2); all other parameters are as in Figure 5.Figure 8
Saddle-node bifurcation illustrated by means of the cubic fixed point equation (19). The functions fL (solid line) and fR (dashed line) of (19) are plotted versus v at the saddle-node bifurcation point in which the state x(low), y(high) emerges or disappears. Parameters: α=2.2; all other parameters are as in Figure 5.


#### 2.2.4. Order Parameters

Our next objective is to determine the order parameters or unstable modes at the saddle-node bifurcation points shown in the bifurcation diagram; see Figure5. From (1), (2), and (9) it follows that the coefficients of the Jacobi matrix J read(22)a11B=r11-K∗-yK∗-x2-1K∗-y,a12B=-r1K∗-x+xK∗-y2,a21D=-r2K∗-y+yK∗-x2,a22D=r21-K∗-xK∗-y2-1K∗-x.We computed the matrix coefficients for the asymptotically stable fixed points x(low), y(high) and x(high), y(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 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⁡(J) from (4). The result is shown in Figure 10. As expected, det(J) converged to zero at the saddle-node bifurcation points. Finally, we computed the eigenvector associated with λ1 from (7) for the asymptotically stable fixed points x(low), y(high) and x(high), y(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.Figure 9
Eigenvaluesλ1 (a) and λ2 (b) computed from (3), (4), and (22) for the state x(low), y(high) (solid lines) and the state x(high), y(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: B=D=1; all other parameters are as in Figure 5.
(a) (b)Figure 10
Determinant of the Jacobian matrixJ as computed from (4) and (22) for the state x(low), y(high) (solid line) and the state x(high), y(low) (dashed line) shown in Figure 5 and as function of α. Parameters: B=D=1; all other parameters are as in Figure 5.Figure 11
Eigenvector components of the first eigenvector and order parameters at the saddle-node bifurcation points. The eigenvector componentsv1,x (a) and v1,y (b) as computed from (3) and (4) are shown for the state x(low), y(high) (solid line) and the state x(high), y(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.
(a) (b)


## 2.1. General Case

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 concentrationsx and y of the two species under consideration measured in arbitrary units. Variations x(t) and y(t) over time t are assumed to satisfy(1)ddtx=fx,y,ddty=gx,y,where f and g 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 saddle-node bifurcation in more detail. To this end, we “follow” the stable fixed point coordinates x(st) and y(st) involved in the saddle-node bifurcation along variations of α. The linearized version of (1) at that stable fixed point reads(2)ddtδ=a11δ+a12ε,ddtε=a21δ+a22εwith x=x(st)+δ and y=y(st)+ε. The variables ε and δ denote small deviations from the fixed point x(st), y(st). The matrix coefficients ajk correspond to the coefficients of the Jacobian matrix J at x(st) and y(st). The eigenvalues λ1 and λ2 of the matrix J are given by(3)λ1,2=TrJ2±TrJ24-det⁡J,where the trace Tr(J) and the determinant det⁡(J) of J are given by(4)TrJ=a11+a22,det⁡J=a11a22-a12a21.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⁡(J)=0, Tr(J)<0, λ1=0, and λ2<0. From det⁡(J)=0 it follows that(5)a11a22=a12a21.Moreover, substituting det⁡(J)=0 into (3) we obtain(6)λ2=TrJ.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 x(st), y(st) in the bistable parameter domain and for the transient dynamics in the monostable domain that describes the transition from the disappearing fixed point x(st), y(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(7)vj,x=11+a11-λj/a122,vj,y=1-vj,x2for j=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 saddle-node 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(8)v1,x=11+a11/a122,v1,y=1-v1,x2.


## 2.2. The Signaling Pathway Model by Markevich et al. (2004)

### 2.2.1. Definition of the Model

Markevich et al. (2004) [5] discussed the following two-species signaling pathway model:(9)ddtx=fx,y=A1-x-yK+1-x-BxK+1-y,ddty=gx,y=C1-x-yK+1-y-DxK+1-xwith A, B, C, D, and K>0. The concentration levels are normalized such that x(t) and y(t) are restricted to the interval [0,1]. The model and generalizations of it have also been discussed in subsequent studies [10, 12].


### 2.2.2. Nullcline Approach

For sake of simplicity, we introduce the new parameters:(10)r1=AB,r2=CD,K∗=K+1.The nullclines are given by(11)x=-h1y2±h12y4-r1y2-1+K∗y+K∗,(12)y=-h2x2±h22x4-r2x2-1+K∗x+K∗with h1(y)=(r2y-K∗(1+r2))/2 and h2(x)=(r1x-K∗(1+r1))/2. 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 byr1=r2=r. 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 r as we will see later. In particular, at a critical value of r, 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 r1 is much larger than r2 the model is monostable with the species x dominating the species y. In contrast, if r2 is much larger than r1 then the model is monostable with the species y dominating the species x. Consequently, if we start initially with a large ratio r2/r1 then the pathway exhibits a high state for y and a low state for x. If we gradually decrease the ratio, then the fixed point for y(high) and x(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 x and a relatively low state for y.Figure 1
Nullclines as defined by (11) (dashed line) and (12) (solid line) for the monostable case. Parameters: K=0.01 (as in [10]) and r=0.7. All parameters are in arbitrary units.Figure 2
It is as in Figure1 but for the bistable case. Parameters: K=0.01 and r=1.5.Figure 3
Illustration of the saddle-node bifurcation (part I) for the asymmetric pathway model. The stable fixed point given byy(high) and x(low) is about to merge with the unstable fixed point. Parameters: K=0.01, r1=2.2, and r2=1.8.Figure 4
Illustration of the saddle-node bifurcation (part II). The pathway has become monostable. Magnification (not shown) of the top left corner of the phase space reveals that the nullclines do not intersect each other in this area. This is because the asymmetry of the pathway model has increased as compared to Figure3. Parameters: K=0.01, r1=2.8, and r2=1.2.Figures1, 2, 3, and 4 illustrate the two basic mechanisms leading to the saddle-node bifurcation underlying the switch-like 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 x(low), y(high) to the state x(high), y(low) or vice versa. Let us analyze these mechanisms in more detail. To this end, we derive a fixed point equation.


### 2.2.3. 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:(13)ABK∗-yK∗-xx=11-x-y=CDK∗-xK∗-yy⟹r1K∗-yK∗-xx=r2K∗-xK∗-yy,(14)AB1-x-yx=K∗-xK∗-y=DCy1-x-y⟹r11-x-yx=1r2y1-x-y.It is useful to introduce the new variables u and v like(15)x=u2,y=v2⇕u=x,v=y.Moreover, we introduce the new parameters ξ and s like(16)ξ=BDAC=1r1r2>0,s=BCAD=r2r1>0.The parameter s measures the degree of asymmetry. In the symmetric case we have s=1. From (14) we obtain the intermediate result:(17)1-u2-v2=ξuvand finally(18)uv;ξ=-12ξv+1+ξ4-1v2.In deriving (18) we used u>0. The variable u as a function of v defined by (18) is a monotonically decaying function for any parameter ξ. For ξ→0 we have u=SQRT(1-v2). Importantly, from (17), it follows that for any parameter ξ the solutions of (18) are symmetric; that is, if u=a and v=b is a solution, then u=b and v=a is a solution as well. From (14) it follows that(19)fLv=fRuv;ξ;swith(20)fLv=K∗-v2v,fRuv;ξ;s=sK∗-u2u.Equation (19) is a cubic equation in the variable v. Therefore, there are maximal three solutions. If s=1 then the equation is symmetric in u and v. Therefore, in view of the symmetry property of (18) for any parameter ξ and the symmetry property of (19) in the special case s=1, we conclude that if s=1 then the number of solutions is either 1 or 3. From (19) it follows that in the symmetric case s=1 (which implies by r1=r2=r) the parameter ξ=1/r2 acts as control parameter and induces a pitchfork bifurcation. For example, for K=0.01 (K∗=1.01) and s=1 numerical evaluation of (19) shows (see below) that for r=2 (19) exhibits three admissible solutions (i.e., solutions in the interval between 0 and 1).For any given set of parametersξ,s, and K, the stable and unstable fixed points x(st) and y(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(21)r1=A+α,r2=B-α,α∈0,B-Awith B>A. In particular, for α=(B-A)/2 we have the symmetric case with s=1 and r1=r2=r=(A+B)/2. We choose the parameters A and B 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 K. Figures 6, 7, and 8 depict the functions fL and fR 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).Figure 5
Bifurcation diagram of the bistable pathway model (9) acting as a switch for a particular value of K when varying the degree of asymmetry. The stationary solutions x(st) are plotted versus the control parameter values α. The solutions x(st) were computed from (18), (19), (20), and (21). Solid line: asymptotically stable fixed points of the state x(low), y(high). Dotted line: unstable fixed points in the bistable parameter domain. Dashed line: asymptotically stable fixed points of the state x(high), y(low). Parameters: K=0.01, A=0.5, and B=3.5.Figure 6
Illustration of the saddle-node bifurcation using the cubic fixed point equation (19). The functions fL (solid line) and fR (dashed line) of (19) are plotted versus v at the saddle-node bifurcation point in which the states x(high) and y(low) emerge or disappear. Parameters: α=0.8; all other parameters are as in Figure 5.Figure 7
The functionsfL (solid line) and fR (dashed line) of (19) in the symmetric case. The signaling pathway model is bistable. Parameters: α=1.5 (r1=r2=r=2); all other parameters are as in Figure 5.Figure 8
Saddle-node bifurcation illustrated by means of the cubic fixed point equation (19). The functions fL (solid line) and fR (dashed line) of (19) are plotted versus v at the saddle-node bifurcation point in which the state x(low), y(high) emerges or disappears. Parameters: α=2.2; all other parameters are as in Figure 5.


### 2.2.4. Order Parameters

Our next objective is to determine the order parameters or unstable modes at the saddle-node bifurcation points shown in the bifurcation diagram; see Figure5. From (1), (2), and (9) it follows that the coefficients of the Jacobi matrix J read(22)a11B=r11-K∗-yK∗-x2-1K∗-y,a12B=-r1K∗-x+xK∗-y2,a21D=-r2K∗-y+yK∗-x2,a22D=r21-K∗-xK∗-y2-1K∗-x.We computed the matrix coefficients for the asymptotically stable fixed points x(low), y(high) and x(high), y(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 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⁡(J) from (4). The result is shown in Figure 10. As expected, det(J) converged to zero at the saddle-node bifurcation points. Finally, we computed the eigenvector associated with λ1 from (7) for the asymptotically stable fixed points x(low), y(high) and x(high), y(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.Figure 9
Eigenvaluesλ1 (a) and λ2 (b) computed from (3), (4), and (22) for the state x(low), y(high) (solid lines) and the state x(high), y(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: B=D=1; all other parameters are as in Figure 5.
(a) (b)Figure 10
Determinant of the Jacobian matrixJ as computed from (4) and (22) for the state x(low), y(high) (solid line) and the state x(high), y(low) (dashed line) shown in Figure 5 and as function of α. Parameters: B=D=1; all other parameters are as in Figure 5.Figure 11
Eigenvector components of the first eigenvector and order parameters at the saddle-node bifurcation points. The eigenvector componentsv1,x (a) and v1,y (b) as computed from (3) and (4) are shown for the state x(low), y(high) (solid line) and the state x(high), y(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.
(a) (b)


## 2.2.1. Definition of the Model

Markevich et al. (2004) [5] discussed the following two-species signaling pathway model:(9)ddtx=fx,y=A1-x-yK+1-x-BxK+1-y,ddty=gx,y=C1-x-yK+1-y-DxK+1-xwith A, B, C, D, and K>0. The concentration levels are normalized such that x(t) and y(t) are restricted to the interval [0,1]. The model and generalizations of it have also been discussed in subsequent studies [10, 12].


## 2.2.2. Nullcline Approach

For sake of simplicity, we introduce the new parameters:(10)r1=AB,r2=CD,K∗=K+1.The nullclines are given by(11)x=-h1y2±h12y4-r1y2-1+K∗y+K∗,(12)y=-h2x2±h22x4-r2x2-1+K∗x+K∗with h1(y)=(r2y-K∗(1+r2))/2 and h2(x)=(r1x-K∗(1+r1))/2. 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 byr1=r2=r. 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 r as we will see later. In particular, at a critical value of r, 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 r1 is much larger than r2 the model is monostable with the species x dominating the species y. In contrast, if r2 is much larger than r1 then the model is monostable with the species y dominating the species x. Consequently, if we start initially with a large ratio r2/r1 then the pathway exhibits a high state for y and a low state for x. If we gradually decrease the ratio, then the fixed point for y(high) and x(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 x and a relatively low state for y.Figure 1
Nullclines as defined by (11) (dashed line) and (12) (solid line) for the monostable case. Parameters: K=0.01 (as in [10]) and r=0.7. All parameters are in arbitrary units.Figure 2
It is as in Figure1 but for the bistable case. Parameters: K=0.01 and r=1.5.Figure 3
Illustration of the saddle-node bifurcation (part I) for the asymmetric pathway model. The stable fixed point given byy(high) and x(low) is about to merge with the unstable fixed point. Parameters: K=0.01, r1=2.2, and r2=1.8.Figure 4
Illustration of the saddle-node bifurcation (part II). The pathway has become monostable. Magnification (not shown) of the top left corner of the phase space reveals that the nullclines do not intersect each other in this area. This is because the asymmetry of the pathway model has increased as compared to Figure3. Parameters: K=0.01, r1=2.8, and r2=1.2.Figures1, 2, 3, and 4 illustrate the two basic mechanisms leading to the saddle-node bifurcation underlying the switch-like 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 x(low), y(high) to the state x(high), y(low) or vice versa. Let us analyze these mechanisms in more detail. To this end, we derive a fixed point equation.


## 2.2.3. 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:(13)ABK∗-yK∗-xx=11-x-y=CDK∗-xK∗-yy⟹r1K∗-yK∗-xx=r2K∗-xK∗-yy,(14)AB1-x-yx=K∗-xK∗-y=DCy1-x-y⟹r11-x-yx=1r2y1-x-y.It is useful to introduce the new variables u and v like(15)x=u2,y=v2⇕u=x,v=y.Moreover, we introduce the new parameters ξ and s like(16)ξ=BDAC=1r1r2>0,s=BCAD=r2r1>0.The parameter s measures the degree of asymmetry. In the symmetric case we have s=1. From (14) we obtain the intermediate result:(17)1-u2-v2=ξuvand finally(18)uv;ξ=-12ξv+1+ξ4-1v2.In deriving (18) we used u>0. The variable u as a function of v defined by (18) is a monotonically decaying function for any parameter ξ. For ξ→0 we have u=SQRT(1-v2). Importantly, from (17), it follows that for any parameter ξ the solutions of (18) are symmetric; that is, if u=a and v=b is a solution, then u=b and v=a is a solution as well. From (14) it follows that(19)fLv=fRuv;ξ;swith(20)fLv=K∗-v2v,fRuv;ξ;s=sK∗-u2u.Equation (19) is a cubic equation in the variable v. Therefore, there are maximal three solutions. If s=1 then the equation is symmetric in u and v. Therefore, in view of the symmetry property of (18) for any parameter ξ and the symmetry property of (19) in the special case s=1, we conclude that if s=1 then the number of solutions is either 1 or 3. From (19) it follows that in the symmetric case s=1 (which implies by r1=r2=r) the parameter ξ=1/r2 acts as control parameter and induces a pitchfork bifurcation. For example, for K=0.01 (K∗=1.01) and s=1 numerical evaluation of (19) shows (see below) that for r=2 (19) exhibits three admissible solutions (i.e., solutions in the interval between 0 and 1).For any given set of parametersξ,s, and K, the stable and unstable fixed points x(st) and y(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(21)r1=A+α,r2=B-α,α∈0,B-Awith B>A. In particular, for α=(B-A)/2 we have the symmetric case with s=1 and r1=r2=r=(A+B)/2. We choose the parameters A and B 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 K. Figures 6, 7, and 8 depict the functions fL and fR 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).Figure 5
Bifurcation diagram of the bistable pathway model (9) acting as a switch for a particular value of K when varying the degree of asymmetry. The stationary solutions x(st) are plotted versus the control parameter values α. The solutions x(st) were computed from (18), (19), (20), and (21). Solid line: asymptotically stable fixed points of the state x(low), y(high). Dotted line: unstable fixed points in the bistable parameter domain. Dashed line: asymptotically stable fixed points of the state x(high), y(low). Parameters: K=0.01, A=0.5, and B=3.5.Figure 6
Illustration of the saddle-node bifurcation using the cubic fixed point equation (19). The functions fL (solid line) and fR (dashed line) of (19) are plotted versus v at the saddle-node bifurcation point in which the states x(high) and y(low) emerge or disappear. Parameters: α=0.8; all other parameters are as in Figure 5.Figure 7
The functionsfL (solid line) and fR (dashed line) of (19) in the symmetric case. The signaling pathway model is bistable. Parameters: α=1.5 (r1=r2=r=2); all other parameters are as in Figure 5.Figure 8
Saddle-node bifurcation illustrated by means of the cubic fixed point equation (19). The functions fL (solid line) and fR (dashed line) of (19) are plotted versus v at the saddle-node bifurcation point in which the state x(low), y(high) emerges or disappears. Parameters: α=2.2; all other parameters are as in Figure 5.


## 2.2.4. Order Parameters

Our next objective is to determine the order parameters or unstable modes at the saddle-node bifurcation points shown in the bifurcation diagram; see Figure5. From (1), (2), and (9) it follows that the coefficients of the Jacobi matrix J read(22)a11B=r11-K∗-yK∗-x2-1K∗-y,a12B=-r1K∗-x+xK∗-y2,a21D=-r2K∗-y+yK∗-x2,a22D=r21-K∗-xK∗-y2-1K∗-x.We computed the matrix coefficients for the asymptotically stable fixed points x(low), y(high) and x(high), y(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 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⁡(J) from (4). The result is shown in Figure 10. As expected, det(J) converged to zero at the saddle-node bifurcation points. Finally, we computed the eigenvector associated with λ1 from (7) for the asymptotically stable fixed points x(low), y(high) and x(high), y(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.Figure 9
Eigenvaluesλ1 (a) and λ2 (b) computed from (3), (4), and (22) for the state x(low), y(high) (solid lines) and the state x(high), y(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: B=D=1; all other parameters are as in Figure 5.
(a) (b)Figure 10
Determinant of the Jacobian matrixJ as computed from (4) and (22) for the state x(low), y(high) (solid line) and the state x(high), y(low) (dashed line) shown in Figure 5 and as function of α. Parameters: B=D=1; all other parameters are as in Figure 5.Figure 11
Eigenvector components of the first eigenvector and order parameters at the saddle-node bifurcation points. The eigenvector componentsv1,x (a) and v1,y (b) as computed from (3) and (4) are shown for the state x(low), y(high) (solid line) and the state x(high), y(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.
(a) (b)


## 3. 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 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 Section2.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.



---

*Source: 8938970-2016-10-31.xml*