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.
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 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 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 concentrations and of the two species under consideration measured in arbitrary units. Variations and over time are assumed to satisfywhere 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 saddle-node bifurcation in more detail. To this end, we “follow” the stable fixed point coordinates and involved in the saddle-node bifurcation along variations of . The linearized version of (1) at that stable fixed point readswith and . The variables ε and δ denote small deviations from the fixed point (st), (st). The matrix coefficients correspond to the coefficients of the Jacobian matrix at and . The eigenvalues and of the matrix are given bywhere the trace and the determinant of are given byNote that in (3) comes with a plus sign in front of the square root function, whereas 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 , , , and . From it follows thatMoreover, substituting into (3) we obtainIt is useful to summarize this intermediate result. Let denote a critical value at which the saddle-node bifurcation occurs. Consequently, at we have and and (5) holds. If we shift next away from the critical value by a small amount such that the pathway is monostable, , 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 we have as compared to . In line with synergetics, the direction associated with the slow dynamics corresponds to the unstable mode at 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 and , the eigenvectors normalized to unity are given in components byfor . At the saddle-node bifurcation point, we have 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
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:with , , , , and . The concentration levels are normalized such that and are restricted to the interval . 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:The nullclines are given bywith and . 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 . 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 is much larger than the model is monostable with the species dominating the species . In contrast, if is much larger than then the model is monostable with the species dominating the species . Consequently, if we start initially with a large ratio 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 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 (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.
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:It is useful to introduce the new variables and likeMoreover, we introduce the new parameters and likeThe parameter measures the degree of asymmetry. In the symmetric case we have . From (14) we obtain the intermediate result:and finallyIn deriving (18) we used . The variable as a function of defined by (18) is a monotonically decaying function for any parameter . For we have . Importantly, from (17), it follows that for any parameter the solutions of (18) are symmetric; that is, if and is a solution, then and is a solution as well. From (14) it follows thatwithEquation (19) is a cubic equation in the variable . Therefore, there are maximal three solutions. If then the equation is symmetric in and . Therefore, in view of the symmetry property of (18) for any parameter and the symmetry property of (19) in the special case , we conclude that if then the number of solutions is either 1 or 3. From (19) it follows that in the symmetric case (which implies by ) the parameter acts as control parameter and induces a pitchfork bifurcation. For example, for () and numerical evaluation of (19) shows (see below) that for (19) exhibits three admissible solutions (i.e., solutions in the interval between 0 and 1).
For any given set of parameters , and , the stable and unstable fixed points and 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 likewith . In particular, for we have the symmetric case with and . 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).
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 Figure 5. From (1), (2), and (9) it follows that the coefficients of the Jacobi matrix readWe 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, converged to zero when approaching the saddle-node bifurcation points, whereas 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 from (4). The result is shown in Figure 10. As expected, converged to zero at the saddle-node bifurcation points. Finally, we computed the eigenvector associated with 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.
(a)
(b)
(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 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.
Competing Interests
The author declares that there is no conflict of interests regarding the publication of this paper.