Detection of biological switches using the method of Gröebner bases

Background Bistability and ability to switch between two stable states is the hallmark of cellular responses. Cellular signaling pathways often contain bistable switches that regulate the transmission of the extracellular information to the nucleus where important biological functions are executed. Results In this work we show how the method of Gröebner bases can be used to detect bistability and output switchability. The method of Gröebner bases can be seen as a multivariate, non-linear generalization of the Gaussian elimination for linear systems which conveniently seperates the variables and drastically simplifies the simultaneous solution of polynomial equations. A necessary condition for fixed-point state bistability is for the Gröbner basis to have three distinct solutions for the state. A sufficient condition is provided by the eigenvalues of the local Jacobians. We also introduce the concept of output switchability which is defined as the ability of an output of a bistable system to switch between two different stable steady-state values. It is shown that bistability does not necessarily guarantee switchability of every state variable of the system. We further show that, for a bistable system, the necessary conditions for output switchability can be derived using the Gröebner basis. The theoretical results are incorporated into an analysis procedure and applied to several systems including the AKT (Protein kinase B), RAS (Rat Sarcoma) and MAPK (Mitogen-activated protein kinase) signal transduction pathways. Results demonstrate that the Gröebner bases can be conveniently used to analyze biological switches by simultaneously detecting bistability and output switchability. Conclusion The Gröebner bases provides a novel methodology to analyze bistability. Results clarify the distinction between bistability and output switchability which is lacking in the literature. We have shown that theoretically, it is possible to have an output subspace of an n-dimensional bistable system where certain variables cannot switch. It is possible to construct such systems as we have done with two reaction networks.


Background
Bistable dynamical systems are frequently encountered in cellular processes. Information processing within cells is carried out by a complex network of switches and oscillators [1]. A bistable system is a system with two attractors. The system can switch between two distinct stable states without resting in intermediate states. Switch-like bistable responses have been observed in many applications including signal transduction [2][3][4][5][6], cell cycle control [7][8][9][10][11][12][13], learning and memory [14], growing bacterial biofilms [15], epileptic spike-wave discharges [16], neurons [17] and synaptic transmission [18]. Bistable switches have been designed synthetically as well. A genetic toggle switch in Escherichia coli has been constructed in [19]. The bistable switch forms an addressable cellular memory unit and has implications for biotechnology, biocomputing and gene therapy.
Considering its biological importance, significant research has been devoted to explaining the physical origin of bistability, to develop necessary conditions for its existence and to construct algorithms for its detection. In particular positive feedback and ultrasensitivity have been proposed as two necessary conditions for the physical appearance of bistability [20]. It is also well known that adding negative feedback to positive feedback can turn bistability into oscillations [21]. The theory of chemical reaction network (CRN) [22] proposes conditions for bistability by making use of the properties of a species-reaction graph. Angeli et al. [23] presented a graphical method to detect bistability for biological positive-feedback systems. Under some mild assumptions, if the open-loop response (when the positive feedback loop is opened) is monotone and has a sigmoidal shape, the system is guaranteed to be bistable for some values of the feedback gains. Finally, Wilhelm [24] proposed a smallest chemical reaction system with bistability. Two reactions constitute a positive feedback loop; a third reaction filters out small stimuli, and a fourth reaction prevents explosions. Analysis is based on the method of Instability Causing Structure Analysis (ICSA) which is based on the Jacobian of the reaction network. Recently a new method was proposed to study multistationarity and bistability of chemical reaction networks with few chemical complexes. The method uses polynomial systems with few distinct monomials and Gale duality [25].
The method of Gröebner bases was introduced by Buchberger in [26,27] as a powerful computational tool to address fundamental questions in commutative algebra (polynomial ideal theory, algebraic geometry). Since its original inception, the method of Gröebner bases was applied to simplify the algorithmic solution of many difficult problems expressed in terms of multivariate polynomials. These include [28]: solving polynomial equations, coding theory, integer programming, partial differential equations, symbolic summation, graph theory and statistics. Today exploring its applicability to many diverse fields such as computational biology [29], chemical kinetics [30][31][32] and systems theory [33][34][35] is an active area of research.
The objective of this work is to explore the method of Gröebner bases to analyze the bistability and output switchability of biological signaling systems. The utility of the method is demonstrated using several examples from the cellular signaling systems literature. Next we give working definitions and examples of bistability and output switchability.
Bistability Consider the dynamical system S f expressed as DAE (Differential Algebraic Equations): where Eq. (1) represents the dynamic mass balances with the nonnegative concentrations of species x and y. Eq.
(2) is a set of algebraic constraints due to the species conservation laws. Substituting y from Eq. (2) into Eq.
(1), one gets: The steady-states of S f are the solutions of 0 = f(x).
A bistable system is a system with two attractors. In general, a plethora of interesting possibilities with different attractors and separatrices can result in different types of bistabilities which coexist in the parameter space of interest. Some of the popular bistabilities are between two stable fixed points [10,11,14], between a stable fixed point and a limit-cycle oscillator [15][16][17]; or between two stable periodic orbits [36,37]. In this paper, we adopt the following definition for bistability: Definition 1: The dynamical system S f is state bistable if it has three nonnegative distinct real steady-states (ss) for the state x, two of which are stable ( x 1;S ss and x 3;S ss Þ and one is unstable x 2;U ss where the superscript S and U denote stable and unstable, respectively. We refer to such bistability as fixed-point bistability to distinguish it from other types of bistabilities mentioned above. Each stable steady-state (or fixed-point) has its own basin of attraction (i.e. the set of initial conditions that asymptotically converge to that steady-state). These basins are separated with a boundary defined by a separatrix. Most often, the separatrix contains a steady-state that is an unstable saddle-point [11]. Upon perturbations in the medium reflected by the parameter changes in the model, other types of attractors like limit cycles can be born from these fixed points via Hopf bifurcation [36]. For example, it was shown in [38] that bistability is a necessary condition for the emergence of oscillations in the MAPK cascade signaling.
Throughout the paper bistability will mean fixed-point state bistability.
A one-dimensional example Figure 1 shows a bistable system resulting from a onedimensional ordinary differential equation dx dt ¼ −x 3 þ 6 x 2 −11x þ 6: There are three positive fixed points at x = 1, 2 and 3 . The unstable fixed point at x= 2 separates the basins of attraction of the stable fixed points. The trajectories starting from the initial conditions to the left of 2 approach the stable fixed point x = 1, and the trajectories starting from the initial conditions to the right of 2 approach the stable fixed point at x = 3.

A two-dimensional example
In [24] a smallest chemical reaction system with bistability was proposed. The model consists of the following two-component mass-action kinetic ODE system: The system has three steady-states, two of which are stable at (0, 0), (6, 4.5) and an unstable steady-state which is a saddle-point at (2, 0.5). Figure 2 shows the phase plane with trajectories emanating from different initial conditions. Due to the saddle-point, the phase plane is divided into two basins of attraction which contain the trajectories approaching the two stable steady-states.
Bistability exits in models with higher dimensions n > 2 as well. For example, in [11] existence of a saddle point and two stable fixed points are highlighted with an apoptosis model that consists of 8 states. The authors present a local analysis to identify the saddle point that helps to understand the global properties of biological switches.

Output switchability
Many biological events are binary with certain variables switching on and off between active and inactive states to perform important biological functions. For a bistable dynamical system S f , one is often interested in the switching response of the concentration of some species. Therefore, we include an output variable y in the model and the analysis. In general output y is taken to be any of the system states x i i = 1 : m. Without loss of generality, the concentration of the first species, or the first state variable x 1 , is defined as the output. The new dynamical system with output x 1 is denoted by S f ;x 1 and expressed as: Next, we introduce the following definition: Definition 2: A bistable dynamical system S f ;x 1 is called output switchable if the steady-state output values x 1, ss are different at the two stable steady-states of the state x i.e. at (x 1;S ss and x 3;S ss Þ. Differentiation between state bistability (Definition 1) and output switchability (Definition 2) is not made in the literature. For a one-dimensional system with a single output, bistability implies output switchability. However, for higher dimensional systems, bistability does not necessarily guarantee switchability for every output variable or state variable x i i = 1 : m. Theoretically, it is possible to have an output subspace containing certain variables that do not switch.

Results
The method of Gröebner bases (see Methods section) is applied to several systems to detect biological switches.
a. Bistable systems with switchable outputs.
Example 1. In [24] a smallest bistable system is given that consists of the following four reactions: The system is described by a two-component massaction ODE system: x y with k 1 = 8, k 2 = k 3 = 1, k 4 = 1.5.
Since the steady-states are determined by the solutions f 1 (x, y) = 0 and f 2 (x, y) = 0, the Gröbner basis is computed for these two polynomials using the reduced Gröebner basis program gbasis available in the Symbolic Math Toolbox of MATLAB: Solving this triangular system and checking the eigenvalues of the Jacobians confirms that the system is bistable with three steady states: Considering x as the output, the system is output switchable since g 1 (x) satisfies the necessary and sufficient conditions (Eqs. 59-61) for switchable outputs given in the Methods section. Plot of g 1 (x) with its three distinct roots is given in Fig. 3.
Note that the inverse problem of given a "desirable" univariate basis polynomial such as (8), reconstruction of a corresponding reaction network is possible although this network is not unique in general. The cubic depletion term −x 3 suggests a bilinear term (−xy) where y is proportional to x 2 so that (−xy) = − x 3 . This can be realized by the following set of reactions: The third reaction provides the cubic depletion rate −k 3 Without this reaction the system can not have three steady-states and bistability is not possible. The first two reactions also provide the quadratic production term x 2 for the species x. Finally, a first order reaction: gives the linear depletion rate −1.5x. Summing up all the terms yields g 1 (x) = x 2 − x 3 8 −1:5x which has the same solutions as the targeted univariate basis (8). Note that in [24] the above smallest bistable system is constructed in a similar way but without introducing the Gröebner basis.
In order to check if the output y is switchable, we change the lexicographic monomial order of the unknown variables (x, y) and recompute the Gröebner basis where the univariate polynomial is now a function of y: Since g 1 (y) satisfies the conditions for switchable outputs (Eqs. (59-61)), the system is switchable in the output y as well.
Example 2. Consider the Edelstein reaction scheme analyzed in [29] and given in Fig. 4.
The system is described by the following DAEs: Substituting the parameter values (see Fig. 4) and eliminating x 3 , one gets: The Gröbner basis is calculated as The system is bistable with three steady states for the state x: In the Methods section we derive that for a bistable system to have a switchable output, the rate of generation of the output must consist of a quadratic and a constant term, and the rate of depletion must consist of a cubic and a linear term. Considering x 1 as the output in the current example, the system is output switchable since g 1 (x 1 ) has quadratic and constant generation, and linear plus cubic depletion terms (see (17)), and it satisfies the necessary and sufficient conditions (Eqs. 59-61). The dynamic output responses of the original system (x 1, f ) and the univariate basis polynomial dynamical system ( x 1;g 1 ) are compared in Fig. 5 for one initial condition. Trajectories converge to the same stable steady-state.
b. Bistable systems with unswitchable outputs.
It is difficult to find physical examples in the literature for bistable systems with unswitchable output(s). There are several reasons for this seemingly lacking data. First, only the switching variables are analyzed to show bistability; therefore, even if there are some outputs (or states) that do not switch, they are not reported. Second, it is possible that some outputs lose their switchability under abnormal conditions only (e.g. disease states due to mutations etc.) that create the right conditions for the emergence of unswitchable output(s). As a result, no distinction is made in the literature between bistability and output switchability. However, as we have presented above, the two concepts are not the same. In addition, at a practical level, for a complex system with many states, it is plausible for some output(s) of a bistable sytem to keep the same steady-state values as the system state switches from one stable steady-state to another. These special outputs can be acting as chaperons that change their values in transient only in order to help the other outputs to switch, and when they complete their tasks they return to their steady-states. Example 3. The smallest bistable system whose output is not switchable is given by a two dimensional system: which satisfies the following conditions at steady-state: Þ¼0 has three distinct nonnegative solutions for z and y ¼ h z ð Þ has two repeated solutions both of which belong to the stable subspace: For example, the following ODEs meet the above conditions: Y In this example we construct a reaction network that satisfies (19) and (20). Consider the reaction network shown in Fig. 6.
The conservation equations with mass action kinetics are given by two ODEs and one algebraic equation: Using the values for the rate constants given in Table 1 and eliminating X via (23) gives: The Gröebner basis G of (f 1 , f 2 ) is calculated using the GroebnerBasis program under Polynomial Algebra of MATHEMATICA and the triangular system of basis polyniomials is given by (compare with Eq. 58): Solving this triangular system of equations yields three distinct steady-state solutions for the state x as shown in Fig. 7.
Therefore, the necessary condition NC for bistability stated in the Methods section (see Eq. 58) is satisfied. Next, one proceeds with the calculation of the Jacobians to establish sufficiency. Both eigenvalues of the Jacobian are negative at steady-states 9 1 ; 9 3 indicating that these are the stable steady-states; one eigenvalue is positive for 10 2 indicating that this is the unstable steady-state. Therefore, the system is bistable. g 1 (Y) = 0 has two solutions, one less than the total number of steady-states for the state. Therefore, one of the roots Y = 9 is necessarily repeated in the steady-state solutions for the state: Since the repeated roots belong to the stable steady-state solutions, the system is not switchable in the output Y.  is not a switchable output since its value remains the same (equal to 9) at the stable steady-states. The reason for Y not to be a switchable output can be physically explained as follows.
The mass balance for Y is determined by the following set of reactions (see Fig. 6): The first two reactions produce Y and the third reaction depletes Y. Conservation of Y is given by the following ODE: The first term k 1 (c − Z)Z is the rate of production of Y by the first reaction, and the second term k 2 (c − Z) is the rate of production of Y by the second reaction. The last term k 3 Y is the rate of consumption of Y which is equal to the sum of the two production rates at steady-state: The total production rate of Y is maximum at the middle unstable steady-state as shown in Fig. 7. At the stable steady-state to the left of the maximum, the first production rate is greater than the second production rate k 1 (c − Z)Z > k 2 (c − Z) and the total production rate is 9. At the stable steady-state to the right of the maximum, the reverse is true i.e. k 1 (c − Z)Z < k 2 (c − Z) but the total production rate k 3 Y remains the same. Since k 3 = 1 (see Table 1), k 3 Y = Y = 9 at the stable steady-states; thus, it cannot switch. The above result shows that, if an output species (Y) is produced by two reactions and the sum of the reactants is constant X + Z = c; then, for some values of the rate constants, the total production rate can remain the same at the stable steady-states leading to unswitchability of Y while both X and Z can switch.
Example 4. Consider the reaction network given in Fig. 8.
The conservation equations with mass action kinetics are given by: The parameter values are given in Table 2.
At steady-state Yf 3 (X, Y, Z) = 0 and one solution is Y = 0, but for this value of Y, there is only one real positive solution for Z at 6.56; therefore, the system cannot be bistable. Thus, we consider the other solutions that satisfy f 3 (X, Y, Z) = 0. If Y is designated as the output variable, the univariate polynomial in Y does not exist. Choosing Z as the output variable and using the lexicographic order (X, Y, Z), the Gröbner basis is computed for the polynomials f 1 (Y, Z), f 2 (X, Z), f 3 (X, Y, Z) and the triangular system of basis polyniomials is given by The steady-state solutions for the state are easily com- The univariate basis polynomial (34) satisfies the output switchability conditions; thus, the system is switchable in output Z. Due to (35) it is switchable in X as well. But it is not switchable in the output Y, since the solution Y = 18 determined from (36) is repeated in the stable steady-states. This is also shown in the bifurcation diagrams in Fig. 9.
c. Cellular signaling pathways: AKT, RAS and MAPK signal transduction systems.
Example 5. AKT signaling pathway AKT signaling pathway plays a key role in the most significant metabolic action of insulin, which is the glucose uptake. Insulin resistance can develop through impairments in the signaling events involved in the activation of AKT. We use the following minimal dimensionless two-state model which we have derived earlier [39] from the original model presented in [6]: where the states are the dimensionless concentrations x 1 = pAKT (active AKT) and x 2 = pIRS (insulin receptor substrate). The input is the amount of insulin represented by λ. The parameter values are taken from [39].
Since the Gröebner basis is defined for polynomials, the right-hand side of (37) is first converted to a rational polynomial function so that (37) and (38) can be expressed as: where both f 1 (x 1 , x 2 ) and f 2 (x 1 , x 2 ) are polynomials. Since the steady-states are determined by the solutions f 1 (x 1 , x 2 ) = 0 and f 2 (x 1 , x 2 ) = 0, the Gröbner basis is computed for these two polynomials using MATHEMATICA, and for the insulin level λ =0.4 it is given by: Since g 1 (x 1 ) given by (39) satisfies all the conditions for three distinct roots (Eqs. 59-61), the output x 1 = pAKT is switchable, if the system is bistable. The necessary condition NC for bistability is satisfied since the triangular system (39)- (40) gives three steady-states: The sufficiency of bistability (i.e. checking the stability status of the three steady-states) is established by bifurcation analysis of eqs. (37)-(38) using XPPAUTO and is given in Fig. 10. The stable and unstable branches show that AKT is bistable for the range of λ between LP1 = 0.38 and λ = LP2 = 0.65. This confirms that the system is indeed bistable for λ = 0.4 and output x 1 = pAKT is switchable.
In order for insulin to perform its function, AKT has to switch between its inactive and active states. Activated AKT (pAKT) enables the translocation of glucose transporter-4 (GLUT-4) from cytosol to the plasma membrane, thus glucose is taken into the cell.
The bistable behavior of pAKT is shown in Fig. 11. pAKT resides on either its active or inactive stable state depending on the initial condition.
The MAPK cascade is an integral part of the ERK (Extracellular Signal-Regulated Kinase) signaling pathway which plays a key role in cell cycle control. In the first stage of the cascade, RAF (Rapidly Accelerated Fibrosarcoma) gets activated by RAS-GTP (Rat Sarcoma Guanosine Triphosphate), and it triggers the second stage where MEK (Mitogen activated protein kinase kinase) gets double phosphorylated [2,38]. This is followed by the activation ERK in the last stage. Here we will focus on the second stage which is shown in Fig. 12.
The two-site MAPK phosphorylation and dephosphorylation cycle with a distributive kinetic mechanism for the kinase and phosphatase possesses the necessary properties to exhibit bistable response [2,40]. MEK and MEK p compete for the same kinase (RAF) for phosphorylation; MEK pp and MEK p compete for the same phosphatase (MEK P'ase) for dephosphorylation. Through this competition, MEK inhibits the production of MEK pp , and MEK pp inhibits the production of MEK. This double inhibition results in a positive feedback loop which leads to bistability under the right set of operating conditions or parameter values. Next we detect and confirm this bistability by using the method of Gröebner bases.
The model is taken from [40]: with the parameters: where c and p are the concentrations of RAF kinase and the MEK phosphatase, respectively. Pertinent data is listed in Table 3. First (41)- (43) are re-expressed as: The Gröebner basis was obtained using MATHEMATICA: A . Evaluation of the eigenvalues of the Jacobians shows that the system is bistable. Since the univariate basis polynomial g 1 (y) meets the conditions for three distinct roots, the system is switchable for the output y ¼ MEK pp M T as well. In fact, all the states are switchable. Example 7. RAS signaling RAS, which is a small GTP (Guanosine Triphosphate) binding protein, serves as an important molecular switch in signaling pathways. For example, in ERK signaling pathway, RAS interacts with the ShC-Grb2-SOS complex, and it is transformed to its active conformation by exchanging GDP (Guanosine Diphosphate) for GTP. Active Ras-GTP starts the sequential phosphorylation of the MAPK pathway that consists of the RAf-MEK-ERK Fig. 11 Bistability of the AKT signaling pathway. AKT system has two stable fixed points (red circles) and an unstable saddle-point (green circle). Trajectories are separated into two basins of attraction of the stable steady-states. λ = 0.4 The model is taken from [41] and it is converted to a dimensionless form. It consists of the following equations: The variables are defined as follows: S is the Shc-Grb2-SOS complex; R T is RAS-GTP; R D is RAS-GDP, SR D and SR T are the complexes formed by the reactions; [.] denotes the concentration. The total concentration of S molecules is α, and the total concentration of RAS molecules is β. The values for the parameters are given in Table 4.
The model has three steady-states, two of which are stable representing the active and inactive states of RAS, and a saddle point with a positive eigenvalue. Bistability    is illustrated in Fig. 13. Trajectories first converge to the unstable manifold, and then they are attracted to either of the two steady-states. MATLAB computes the univariate basis polynomial as a cubic polynomial: with three distinct roots, thus R T is a switchable output.

Discussion
We have developed a new method to detect and analyze biological switches by simultaneously treating bistability and output switchability using the Gröebner bases. As demonstrated by several examples, the proposed methodology offers the following: First the method provides computational advantages due to its nice properties. Specifically, the method of Gröebner bases to solve polynomial systems can be seen as a multivariate, non-linear generalization of the Gaussian elimination for linear systems. Multistationarity is easily checked by solving a triangular set of equations which facilitates the root finding. Bistability is confirmed by local stability analysis using the Jacobian which is also straightforward. It provides a theoretical framework and a systematic methodology that analyzes both bistability and output switchability simultaneously. Output switchability conditions follow immediately from the univariate Gröebner polynomial basis and are easy to check. We show by Examples 3 and 4 that some bistable systems can have outputs that do not switch their steady-states. The univariate Gröebner basis polynomial provides useful biological insight which can help in the design of biological switches. A bistable dynamical system with output x, is output switchable, if its univariate basis polynomial g 1 (x) = − x 3 + bx 2 − cx + d with b > 0, c > 0, d ≥ 0 has three distinct nonnegative roots. This result provides some biological insight. The terms −x 3 − cx represent the rate of depletion of species x, and the terms bx 2 + d represent the rate of production of species x. This suggests that a biological switch for an output species x can be designed by constructing a reaction network (with its corresponding ODEs) whose univariate Gröbner basis polynomial in x has the above types of depletion and production terms. In fact Examples 1, 3 and 4 were constructed in this fashion.

Conclusions
We have presented a new method to detect biological switches by analyzing their bistability and output switchability properties. The methodology is based on the Gröebner bases. Conditions are established to make the connections between the Gröebner bases, bistability and output switchability. Various examples are given to elucidate the theoretical results. We show that the method can analyze bistability and output switchability while providing useful insight into the underlying mechanisms. The method is easy to apply since significant software such as MAPLE, MATLAB and MATHEMATICA exists to perform the Gröebner bases computation. It goes without saying that high dimensionality can pose computational problems as in other methods. As a remedy, techniques such as lumping, network complexity reduction can be used to reduce the number of ODEs before the Gröebner basis calculation is carried out. In general the method can be applied to other types of polynomial differential equations derived from data instead of first principles. Such potential models include the nonlinear polynomial regression models.

Methods
Bistability and output switchability analysis is based on the Gröebner bases. Fig. 13 Bistability of a three-dimensional RAS model. The system has two stable fixed points (blue stars) and an unstable saddle-point (magenta circle). Trajectories reach the stable fixed points after following the unstable manifold of the saddle point. S is the Shc-Grb2-SOS complex. R T is RASGTP. SR T is the SOS-RASGTP complex

The Gröebner bases
Any finite set of multivariate polynomials F can be transformed by an algorithm (see Buchberger's algorithm [26]) into another set of basis polynomials, called the Gröebner basis G. Many problems that are difficult to handle by the original set of polynomials can be easily solved by using the method of Gröebner bases due to its "nice" properties. Readily available computer software such as MAPLE, MATHEMATICA and MATLAB are equipped with the computational machinery of the Gröebner bases. The most basic definitions and properties of the Gröebner bases are presented in the Additional file 1. In this paper we explore the Gröebner bases within the context of bistability analysis. To this end, we state some of the useful properties of the Gröebner bases for solving polynomial equations.

Solution of polynomial equations by the Gröebner bases
The method of Gröebner bases to solve polynomial systems can be seen as a multivariate, non-linear generalization of the Gaussian elimination for linear systems [42].
The ideal I = < f 1 , f 2 , . . , f n > is the set of all possible linear combinations of f i ′ s where the coefficients are polynomials p i (Additional file 1). Since F = (f 1 , f 2 , , .., f n ) and its Gröebner basis G generate the same ideal, they have the same solutions [42,43]. The advantageous property of the Gröebner basis G is that it yields a triangular system which conveniently seperates the variables and drastically simplifies the calculation. This triangular system is like the reduced row echelon form obtained by pivoting in Gaussian elimination in the case of linear systems.
Consider the steady-state solutions of the dynamical system (Eqs. [4][5] which satisfy following set of polynomial equations: The Gröebner basis G for these two polynomials with respect to the lexicographic ordering is given in a triangular form: First the univariate basis polynomial g 1 (x) is easily solved for its roots x: (0, 6,2). Next these x values are substituted into the bivariate basis g 2 (x, y) to determine its corresponding roots y: (0,4.5,0.5). Thus, the solutions (x, y) of the original set of polynomials F are obtained: (0, 0), (6,4.5) and (2,0.5).

Detection of Bistability
Consider the dynamical system S f given by: x where the steady-state solutions satisfy f(x) = 0, and they are denoted as x i ss ¼ ½x i 1;ss x i 2;ss ……:x i n;ss T i ¼ 1 : m with m the number of solutions. Let G = [g 1 (x 1 ), g(x)] be the Gröebner basis for the ideal I = < f 1 , f 2 , …, f n >, where g 1 (x 1 ) is the univariate basis polynomial, and g(x) is the vector of remaining polynomials arranged in triangular form: A necessary condition for bistability (NC) The dynamical system S f is bistable only if the following set of equations have three distinct real nonnegative solutions ðx 1 ss ; x 2 ss ; x 3 ss Þ for the state x: Necessity follows from the working definition of bistability which requires three distinct steady-state solutions for the state vector x, and the fact that the steadystate solutions of S f and the solutions of the Gröebner basis polynomials (58) are the same. Moreover, for a zero-dimensional ideal, the triangular structure in the Gröebner basis always exists [42]. For Gröbner bases, unlike other triangular systems, it is guaranteed that each partial solution can be extended to a full solution. This means that every solution x 1 of the first polynomial can be extended to a solution (x 1 , x 2 ) of the polynomials in x 1 and x 2 , and each of these solutions can be further extended to a solution (x 1 , x 2 , x 3 ) of the polynomials in It is important to note that bistability cannot be detected by checking the number of solutions of the univariate basis polynomial g 1 (x 1 ) alone but the whole basis must be considered. This follows from the fact that the number of solutions of g 1 (x 1 ) = 0 can be less than the number of solutions for the state x.
Systems that fail to meet the necessary condition cannot be bistable; thus, they are easily eliminated from further consideration. However, satisfying the necessary condition does not guarantee bistability. Further tests should be applied to confirm it. The most common approach is to compute the eigenvalues of the Jacobian of f(x). The Jacobian matrix is obtained by linearizing the dynamical system S f at its steady-states x i ss : ss of the dynamical system S f is stable if all the eigenvalues of Jacobian matrix J i have negative real parts. The steady-state is unstable if at least one of the eigenvalues has a positive real part [44]. Bistability can be easily ascertained by checking the stability status of each of the three distinct steady-states ðx 1 ss ; x 2 ss ; x 3 ss Þ of S f .

Detection of switchable outputs
A bistable dynamical system with output x 1 ; S f ;x 1 ; is output switchable, if its univariate basis polynomial satisfies the following conditions: These conditions guarantee that g 1 (x 1 ) has three distinct nonnegative roots so that, the output x 1 can take different values when the state x changes between its stable steady-states. Existence of three distinct nonnegative roots can be shown as follows. According to the Descartes' rule of signs, the maximum number of negative real roots of a polynomial f(x) is equal to the number of changes in sign of the coefficients of the terms of f(−x). When (60) is true, there are no sign changes in the coefficients of g 1 (−x 1 ) in (59); thus, g 1 (x 1 ) can have maximum three nonnegative roots. The inequality (61) is the cubic discriminant condition which guarantees that there are three distinct real roots. The general graph of g 1 (x 1 ) satisfying conditions (59)-(61) is shown in Fig. 14.
If the univariate basis polynomial g 1 (x 1 ) does not have three distinct roots; then, the bistable system S f ;x i is switchable in its output x i , only if one of the repeated solutions belongs to the unstable steady-state solution x 2;U ss . This follows from the definition of output switchability which requires that the output values are different at the two stable steady-states.
It is possible that for an output variable of interest y, the univariate basis polynomial g 1 (y) may not exist. If this happens, by changing the lexicographic ordering, the triangular system of basis polyniomials (58) is calculated using a different univariate basis polynomial g 1 (x i ) where x i ≠ y. In such cases, the bistable system is switchable in its output y, only if the solutions y obtained from the triangular system of basis polynomials are different at the two stable steady-states (x 1;S ss and x 3;S ss Þ.
The univariate basis polynomial dynamics For a bistable system with a switchable output x 1 , we define the univariate basis polynomial dynamics describing the output as: S g 1 ;x 1 is bistable as seen from the signs of the local Jacobians depicted in Fig. 14. Note that the output trajectories of S g 1 ;x 1 and the original system S f ;x 1 are different in transient but both converge to the same stable steadystate, since they have the same output solutions due to the property of Gröbner basis. When the univariate basis