 Methodology article
 Open Access
 Published:
Detection of biological switches using the method of Gröebner bases
BMC Bioinformatics volume 20, Article number: 615 (2019)
Abstract
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, nonlinear 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 fixedpoint 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 steadystate 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 (Mitogenactivated 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 ndimensional 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. Switchlike 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 spikewave 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 speciesreaction graph. Angeli et al. [23] presented a graphical method to detect bistability for biological positivefeedback systems. Under some mild assumptions, if the openloop 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 steadystates 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 limitcycle 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 steadystates (ss) for the state x, two of which are stable (\( {x}_{ss}^{1,S}\ \mathrm{and}\ {x}_{ss}^{3,S}\Big) \) and one is unstable \( {x}_{ss}^{2,U} \) where the superscript S and U denote stable and unstable, respectively.
We refer to such bistability as fixedpoint bistability to distinguish it from other types of bistabilities mentioned above. Each stable steadystate (or fixedpoint) has its own basin of attraction (i.e. the set of initial conditions that asymptotically converge to that steadystate). These basins are separated with a boundary defined by a separatrix. Most often, the separatrix contains a steadystate that is an unstable saddlepoint [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 fixedpoint state bistability.
A onedimensional example
Figure 1 shows a bistable system resulting from a onedimensional ordinary differential equation \( \frac{dx}{dt}={x}^3+6{x}^211x+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 twodimensional example
In [24] a smallest chemical reaction system with bistability was proposed. The model consists of the following twocomponent massaction kinetic ODE system:
The system has three steadystates, two of which are stable at (0, 0), (6, 4.5) and an unstable steadystate which is a saddlepoint at (2, 0.5). Figure 2 shows the phase plane with trajectories emanating from different initial conditions. Due to the saddlepoint, the phase plane is divided into two basins of attraction which contain the trajectories approaching the two stable steadystates.
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 steadystate output values x_{1, ss} are different at the two stable steadystates of the state x i.e. at (\( {x}_{ss}^{1,S}\ and\ {x}_{ss}^{3,S}\Big) \).
Differentiation between state bistability (Definition 1) and output switchability (Definition 2) is not made in the literature. For a onedimensional 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 twocomponent massaction ODE system:
with k_{1} = 8, k_{2} = k_{3} = 1, k_{4} = 1.5.
Since the steadystates 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:
Steadystate mass balance for y gives \( y=\left(\frac{k_2}{k_1}\right){x}^2=\frac{x^2}{8}. \) The third reaction provides the cubic depletion rate \( {k}_3 xy={k}_3\left(\frac{k_2}{k_1}\right){x}^3=\frac{x^3}{8} \) . Without this reaction the system can not have three steadystates 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}− \( \frac{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 steadystate.
 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 steadystate values as the system state switches from one stable steadystate 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 steadystates.
Example 3. The smallest bistable system whose output is not switchable is given by a two dimensional system:
\( \dot{z}=f\left(z,y\right) \)
which satisfies the following conditions at steadystate:
For example, the following ODEs meet the above conditions:
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 steadystate solutions for the state \( x=\left(\begin{array}{c}Y\\ {}Z\end{array}\right) \) = \( \left(\begin{array}{c}9\\ {}1\end{array}\right),\left(\begin{array}{c}10\\ {}2\end{array}\right),\left(\begin{array}{c}9\\ {}3\end{array}\right) \) 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 steadystates \( \left(\begin{array}{c}9\\ {}1\end{array}\right),\left(\begin{array}{c}9\\ {}3\end{array}\right) \) indicating that these are the stable steadystates; one eigenvalue is positive for \( \left(\begin{array}{c}10\\ {}2\end{array}\right) \) indicating that this is the unstable steadystate. Therefore, the system is bistable.
g_{1}(Y) = 0 has two solutions, one less than the total number of steadystates for the state. Therefore, one of the roots Y = 9 is necessarily repeated in the steadystate solutions for the state: \( x=\left(\begin{array}{c}Y\\ {}Z\end{array}\right) \) = \( \left(\begin{array}{c}9\\ {}1\end{array}\right),\left(\begin{array}{c}10\\ {}2\end{array}\right),\left(\begin{array}{c}9\\ {}3\end{array}\right) \). Since the repeated roots belong to the stable steadystate solutions, the system is not switchable in the output Y. Figure 7 shows the state trajectories x(t) approaching the stable steadystates separated by the middle unstable fixed point at \( \left(\begin{array}{c}10\\ {}2\end{array}\right) \).
Figure 7 shows the state trajectories x(t) approaching the stable steadystates separated by the middle unstable fixed point at \( \left(\begin{array}{c}10\\ {}2\end{array}\right) \). However, unlike Z and X = cZ, Y is not a switchable output since its value remains the same (equal to 9) at the stable steadystates.
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):
and the constraint X + Z = c.
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 steadystate:
The total production rate of Y is maximum at the middle unstable steadystate as shown in Fig. 7. At the stable steadystate 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 steadystate 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 steadystates; 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 steadystates 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 steadystate 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 steadystate solutions for the state are easily computed: \( x=\left({x}_{ss}^{1,S},{x}_{ss}^{2,U},{x}_{ss}^{3,S}\right)=\left(\begin{array}{c}Z\\ {}X\\ {}Y\end{array}\right)=\left(\begin{array}{c}2\\ {}4\\ {}18\end{array}\right),\left(\begin{array}{c}3\\ {}9\\ {}19\end{array}\right),\left(\begin{array}{c}4\\ {}16\\ {}18\end{array}\right). \) The system is bistable as determined by the Jacobians.
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 steadystates. 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 twostate 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 righthand side of (37) is first converted to a rational polynomial function so that (37) and (38) can be expressed as:
\( \dot{\ {x}_1} \)= \( \frac{f_1\left({x}_1,{x}_2\right)}{d_1\left({x}_1,{x}_2\right)} \) and \( \dot{\ {x}_2} \) = f_{2}(x_{1}, x_{2})
where both f_{1}(x_{1}, x_{2}) and f_{2}(x_{1}, x_{2}) are polynomials. Since the steadystates 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 steadystates:
\( \left(\begin{array}{c}{x}_{1, ss}\\ {}{x}_{2, ss}\end{array}\right) \) = \( \left(\begin{array}{c}0.1689\\ {}1.2228\end{array}\right),\left(\begin{array}{c}0.310\\ {}1.083\end{array}\right),\left(\begin{array}{c}0.964\\ {}0.435\end{array}\right) \)
The sufficiency of bistability (i.e. checking the stability status of the three steadystates) 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 transporter4 (GLUT4) 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 SignalRegulated 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 RASGTP (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 twosite 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]:
\( \left(x,y,z\right)\ \mathrm{are}\ \mathrm{the}\ \mathrm{dimensionless}\ \mathrm{concentrations}\ \frac{MEK}{M_T},\frac{MEK^{pp}}{M_T}\ \mathrm{and}\ \frac{MEK^p}{M_T} \), (respectively)
M_{T} is the total concentration of MEK. The rates are given by:
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 reexpressed as:
The Gröebner basis was obtained using MATHEMATICA:
which can be solved easily to give three solutions:
\( \left(\begin{array}{c}y\\ {}x\\ {}z\end{array}\right)=\left(\begin{array}{c}0.0033\\ {}0.9655\\ {}0.0312\end{array}\right),\left(\begin{array}{c}0.1421\\ {}0.6969\\ {}0.1610\end{array}\right),\left(\begin{array}{c}0.9347\\ {}0.0065\\ {}0.0588\end{array}\right) \). 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=\frac{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 ShCGrb2SOS complex, and it is transformed to its active conformation by exchanging GDP (Guanosine Diphosphate) for GTP. Active RasGTP starts the sequential phosphorylation of the MAPK pathway that consists of the RAfMEKERK signaling cascade. Catalytic activation of RAS by the SOS (Son of Sevenless) complex ShcGrb2SOS while RASGTP is bound to its allosteric site creates a positive loop resulting in a bistable switching response of RasGTP [41]. 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 ShcGrb2SOS complex; R_{T} is RASGTP; R_{D} is RASGDP, 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 steadystates, 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 steadystates.
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, nonlinear 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 steadystates.
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.
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, nonlinear 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 steadystate 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:
where the steadystate solutions satisfy f(x) = 0, and they are denoted as
\( {x}_{ss}^i={\left[{x}_{1, ss}^i\ {x}_{2, ss}^i\dots \dots .{x}_{n, ss}^i\ \right]}^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 \( \left({x}_{ss}^1,{x}_{ss}^2,{x}_{ss}^3\right) \) for the state x:
Necessity follows from the working definition of bistability which requires three distinct steadystate 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 zerodimensional 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 x_{1}, x_{2}, x_{3}, etc.
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 steadystates \( {x}_{ss}^i \):
A steadystate \( {x}_{ss}^i \) of the dynamical system S_{f} is stable if all the eigenvalues of Jacobian matrix J^{i} have negative real parts. The steadystate 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 steadystates \( \left({x}_{ss}^1,{x}_{ss}^2,{x}_{ss}^3\right) \) 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 steadystates. 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 steadystate solution \( {x}_{ss}^{2,U} \) . This follows from the definition of output switchability which requires that the output values are different at the two stable steadystates.
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 steadystates (\( {x}_{ss}^{1,S}\ and\ {x}_{ss}^{3,S}\Big) \).
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}\ \mathrm{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 polynomial has three distinct roots, it is always possible to construct the onedimensional bistable system \( {S}_{g_1,{x}_1}. \)
It can be seen from Eq. (62) that g_{1}(x_{1}) is the imbalance between the rate at which the output species is generated and the rate at which it is depleted:
At the roots of g_{1}(x_{1}), the two rates equilibrate. Between the roots, the sign of (g_{generation} − g_{depletion}) alternates as (+ − +) to create a bistable switching output response as shown in Fig. 15.
According to (63) for a bistable system to have a switchable output x_{1}, the rate of generation must consist of a quadratic and a constant term, and the rate of depletion must consist of a cubic and a linear term:
This suggests that a biological switch for an output species can be synthetically designed by constructing a reaction network (and its corresponding ODEs) who’s univariate Gröbner basis polynomial has the above types of depletion and production terms. Examples 1, 3 and 4 were constructed in this fashion.
Analysis procedure
We have incorporated the theoretical results into the following analysis procedure:
 1.
Given a Differential Algebraic System, compute its Gröbner basis. We have used MATLAB and MATHEMATICA for this purpose.
 2.
By solving the triangular system of eqs. (58), check if three distinct nonnegative solutions exist in some region of the state space. If it does, check the eigenvalues of the Jacobian at the three steadystates to decide if the system is bistable. If the system is bistable, proceed to the next step. If three distinct nonnegative solutions for the state do not exist, the system cannot be (fixedpoint) bistable and stop.
 3.
For the bistable system, identify the output of interest y and compute the univariate basis polynomial in y i.e. g_{1}(y) by reordering the variables, if necessary. Compute the roots of g_{1}(y). If the number of roots is three, the bistable system is switchable in the output if these roots are all distinct. In all other cases where there are two repeated roots, one of the repeated roots must belong to the unstable steadystate solution \( {x}_{ss}^{2,U} \). Otherwise the system is not switchable in output y.
 4.
In case the univariate basis polynomial g_{1}(y) does not exist, solutions y are calculated using a different univariate basis polynomial g_{1}(x_{i}) where x_{i} ≠ y. The bistable system is switchable in its output y, only if the solutions at two stable steadystates (\( {x}_{ss}^{1,S}\ and\ {x}_{ss}^{3,S}\Big) \) are distinct.
Availability of data and materials
All data generated or analysed during this study are included in this manuscript. Given the data, the Gröbner bases were calculated by calling the GroebnerBasis routine available in MATHEMATICA or the gbasis routine available in the Symbolic Math Toolbox of MATLAB. No special software is needed except to enter the data reported in the manuscript in MATLAB or MATHEMATICA.
Abbreviations
 AKT:

Protein kinase B
 CRN:

Chemical Reaction Network
 DAE:

Differential Algebraic Equations
 ERK:

Extracellular SignalRegulated Kinase
 GDP:

Guanosine Diphosphate
 GLUT4:

Glucose transporter4
 Grb2:

Growth factor receptorbound protein 2
 GTP:

Guanosine Triphosphate
 ICSA:

Instability Causing Structure Analysis
 IRS:

Insulin receptor substrate
 LP:

Limit Point
 MAPK:

Mitogen activated protein kinase
 MEK:

Mitogen activated protein kinase kinase
 NC:

Necessary condition
 ODE:

Ordinary Differential Equations
 pAKT:

Phosphorylated AKT
 pIRS:

Phosphorylated Insulin Receptor Substrate
 RAF:

Rapidly Accelerated Fibrosarcoma
 RAS:

Rat Sarcoma
 RASGTP:

Active form of RAS
 RGAP:

RAS GTPase Activating Proteins
 Shc:

Src homology 2 domain containing
 SOS:

Son of Sevenless
References
 1.
Tyson JJ, Albert R, Goldbeter A, Ruoff AP, Sible J. Biological switches and clocks. J R Soc Interface. 2008;5:S1–8.
 2.
Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol. 2004;164:353–9.
 3.
Ferrell JE Jr. The biochemical basis of an allornone cell fate switch in Xenopus oocytes. Science. 1998;280(5365):895–8.
 4.
Arkun Y, Yasemi M. Dynamics and control of the ERK signaling pathway: sensitivity, bistability, and oscillations. PLoS One. 2018;13(4):e0195513.
 5.
Giri L, Mutalik VK, Venkatesh KVA. Steady state analysis indicates that negative feedback regulation of PTP1B by Akt elicits bistability in insulinstimulated GLUT4 translocation. Theor Biol Med Model. 2004;1:1–16.
 6.
Wang G. Singularity analysis of the AKT signaling pathway reveals connections between cancer and metabolic diseases. Phys Biol. 2010;7:046015.
 7.
Eissing T, Conzelmann H, Gilles ED, Allgöwer F, Bullinger E, Scheurich P. Bistability analyses of a caspase activation model for receptorinduced apoptosis. J Biol Chem. 2004;279(35):36892–7.
 8.
Nakakuki T, Birtwistle MR, Saeki Y, Yumoto N, Ide K, Nagashima T, et al. Ligandspecific cfos expression emerges from the spatiotemporal control of ErbB network dynamics. Cell. 2010;141(5):884–96.
 9.
Yao G, Lee TJ, Mori S, Nevins JR, You L. A bistable Rb–E2F switch underlies the restriction point. Nat Cell Biol. 2008;10(4):476–82.
 10.
Verdugo A, Vinod PK, Tyson JJ, Novak B. Molecular mechanisms creating bistable switches at cell cycle transitions. Open Biol. 2013;3:120179.
 11.
Trotta L, Sepulchre R, Bullinger E. Global analysis of dynamical decisionmaking models through local computation around the hidden saddle. PLoS One. 2012;7(3):e33110. https://doi.org/10.1371/journal.pone.0033110.
 12.
Zhang T, Brazhnik P, Tyson JJ. Computational analysis of dynamical responses to the intrinsic pathway of programmed cell death. Biophys J. 2009;97:415–34.
 13.
Bala SI, Ahmad NMR. Comp Appl Math. 2018;37:266. https://doi.org/10.1007/s4031401704674.
 14.
Song H, Smolen P, AvRon E, Baxter DA, Byrne JH. Bifurcation of singularity analysis of a molecular network for the induction of longterm memory. Biophys J. 2006;9:2309–25.
 15.
MartinezCorral R, Liu J, Suel G, GarciaOjalvo J. Bistable emergence of oscillations in structured cell populations. BioRxiv. 2018. https://doi.org/10.1101/276113.
 16.
Fan D, Liu S, Wang Q. Epileptic stimulusinduced epileptic spikewave discharges in thalamocortical model with disinhibition. Sci Rep. 2016;6:1–21.
 17.
Dovzhenok A, Kuznetsov AS. Exploring neuronal bistability at the depolarization block. PLoS One. 2012;7(8):e42811. https://doi.org/10.1371/journal.pone.0042811.
 18.
Byrne JH, Heidelberger R, Waxham MN. (eds)From Molecules to Networks. An Introduction to Cellular and Molecular Neuroscience: Academic. Academic Press, Elsevier; 2014.
 19.
Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in Escherichia coli. Nature. 2000;403(6767):339–42.
 20.
Ferrell JE, Ha SH. Ultrasensitivity part II: multisite phosphorylation, stoichiometric inhibitors, and positive feedback. Trends Biochem Sci. 2014;39(11):556–69.
 21.
Thomas R. On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations. Springer Ser Synergetics. 1981;9:180–93.
 22.
Craciun G, Tang Yand Feinberg M. Understanding bistability in complex enzymedriven reaction networks. Proc Natl Acad Sci. 2006;103(23):8697–702.
 23.
Angeli D, Ferrell JE, Sontag ED. Detection of multistability, bifurcations, and hysteresis in a large class of biological positivefeedback systems. Proc Natl Acad Sci U S A. 2004;101:1822–7.
 24.
Wilhelm T. The smallest chemical reaction system with bistability. BMC Syst Biol. 2009;3:90.
 25.
Feliu E, Helmer M. Multistationarity and bistability for fewnomial chemical reaction networks. M Bull Math Biol. 2019;81:1089. https://doi.org/10.1007/s1153801800555z.
 26.
Buchberger B. An Algorithm for Finding the Bases Elements of the Residue Class Ring Modulo a Zero Dimensional Polynomial Ideal (German). PhD thesis. Austria: University of Innsbruck; 1965.
 27.
Buchberger B. An Algorithmical Criterion for the Solvability of Algebraic Systems of Equations (German). Aequationes Math. 1970;4(3):374–83.
 28.
Buchberger B, Winkler F. editorsGröebner Bases and Applications, volume 251 of London MATHEMATICAl Society Series, Proc. of the International Conference “33 Years of Gröebner Bases”: Cambridge University Press. London Mathematical Society Lecture Note Series; 1998.
 29.
MartínezForero I, Pelá EzLó Pez A, Villoslada P. Steady State Detection of Chemical Reaction Networks Using a Simplified Analytical Method. PLoS One. 2010;5(6):5.
 30.
Minimair M, Barnett MP. Solving polynomial equations for chemical problems using Gröebner bases. Mol Phys. 2004;102:2521535.
 31.
Mercedes PM, Dickenstein A, Shiu A, Conradi C. Chemical Reaction Systems with Toric Steady States. Bull Math Biol. 2012;74:1027–65.
 32.
Grimbs S, Arnolda A, Koseskac A, Kurths J, Selbiga J, Nikoloski Z. Spatiotemporal dynamics of the Calvin cycle: Multistationarity and symmetry breaking instabilities. BioSystems. 2011;103:212–23.
 33.
Calandrini GL, Paolini EE, Moiola JL. Gröebner bases for designing dynamical systems. Lat Am Appl Res. 2003;33:4 Bahía Blanca.
 34.
Zhiping L, Xu L, Bose NK. A tutorial on Gröebner bases with applications in signals and systems. IEEE Trans Circuits Syst I. 2008;55(1):445–61.
 35.
Wenz M, Wörn H. Solving the inverse kinematics problem symbolically by means of knowledgebased and linear algebrabased methods, IEEE International Conference on Emerging Technologies and Factory Automation, ETFA; 2007. p. 1346–53.
 36.
Guevara MR. Bifurcations Involving Fixed Points and Limit Cycles in Biological Systems. In: Beuter A, Glass L, Mackey MC, Titcombe MS, editors. (eds) Nonlinear Dynamics in Physiology and Medicine. Interdisciplinary Applied Mathematics, vol. 25. New York: Springer; 2003.
 37.
Abraham R, Shaw CD. Dynamicsthe geometry of behavior: periodic behavior: Aerial Press. Basic Books; 1982.
 38.
Qiao L, Nachbar RB, Kevrekidis IG, Shvartsman SY. Bistability and oscillations in the HuangFerrell model of MAPK signaling. PLoS Comput Biol. 2007;3(9):1819–26.
 39.
Cizmeci D, Arkun Y. Regulatory networks and complex interactions between the insulin and angiotensin II signaling systems: Models and implications for hypertension and diabetes. PLoS One. 2013;8(12):e83640.
 40.
Ortega F, Garcés JL, Mas F, Kholodenko BN, Cascante M. Bistability from double phosphorylation in signal transduction: kinetic and structural requirements. FEBS J. 2006;273(17):3915–26.
 41.
Das J, Ho M, Zikherman J, Govern C, Yang M, Weiss A, et al. Digital signaling and hysteresis characterize Ras activation in lymphoid cells. Cell. 2009;136(2):337–51.
 42.
Adams WW, Loustaunau P. An Introduction to Gröbner Bases. Graduate Studies in Mathematics, vol. 3: American MATHEMATICAl Society. American Mathematical Society; 2000.
 43.
Winkler F. The method of Gröebner bases. In Polynomial Algorithms in Computer Algebra, Texts and Monographs in Symbolic Computation, chapter 8. Wien: Springer; 1996.
 44.
Khalil HK. Nonlinear systems: Prentice Hall. Pearson Higher Education; 2014.
Acknowledgements
The financial support through grant 117F123 by TÜBİTAK, The Scientific and Technological Research Council of Turkey, is gratefully acknowledged.
Funding
The work was supported through grant 117F123 by TÜBİTAK, The Scientific and Technological Research Council of Turkey. TÜBİTAK did not play any role in the design of the study and collection, analysis, and interpretation of data, or in writing this manuscript.
Author information
Affiliations
Contributions
YA is solely responsible for developing the methodology, collecting data, obtaining results and preparing the manuscript. The author read and approved the final manuscript.
Corresponding author
Correspondence to Yaman Arkun.
Ethics declarations
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable
Competing interests
The author declares that he has no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Arkun, Y. Detection of biological switches using the method of Gröebner bases. BMC Bioinformatics 20, 615 (2019) doi:10.1186/s1285901931550
Received
Accepted
Published
DOI
Keywords
 Bistability
 Output switchability
 The Gröebner bases
 Univariate basis polynomial
 Steadystate solutions
 Bifurcation
 Polynomial equations
 Biomolecular reactions