A feedback approach to bifurcation analysis in biochemical networks with many parameters

A feedback approach to bifurcation analysis in biochemical networks with many parameters


Feedback circuits in biochemical networks which underly cellular signaling pathways are important elements in creating complex behavior. A specific aspect thereof is how stability of equilibrium points depends on model parameters. For biochemical networks, which are modelled using many parameters, it is typically very difficult to estimate the influence of parameters on stability. Finding parameters which result in a change in stability is a key step for a meaningful bifurcation analysis. We describe a method based on well known approaches from control theory, which can locate parameters leading to a change in stability. The method considers a feedback circuit in the biochemical network and relates stability properties to the control system obtained by loop–breaking. The method is applied to a model of a MAPK cascade as an illustrative example.


and 1mail]Corresponding author (waldherr@ist.uni-stuttgart.de)

eedback circuit, stability, bifurcation analysis, MAPK cascade

1 Introduction

Feedback circuits are an important structural feature of biochemical networks (Tyson and Othmer, 1978). The presence of complex behavior such as bistability, i.e. the existence of several stable equilibria, and sustained oscillations can be attributed to the presence of feedback circuits (Cinquin and Demongeot, 2002; Kaufman and Thomas, 2003). These types of complex behavior are directly related to how feedback circuits influence stability properties of equilibria. In consequence, stability analysis of biochemical networks involving feedback is a recurring field of interest, and several theoretical results have been obtained (Dibrov et al., 1982; Thron, 1991; Angeli and Sontag, 2004).

Models for biochemical networks in cellular signaling typically contain a large number of parameters whose values are not exactly known and which can even vary due to differential gene expression (e.g. the concentration of an enzyme) or external influences (e.g. cofactors). These parameters often have a considerable influence on stability, which needs to be evaluated in order to understand the function of a network (Eißing et al., 2007; Kim et al., 2006).

A classical tool to study the influence of parameter variations on stability is bifurcation analysis. It has been applied to many cellular signalling systems, such as the lac operon (Yildirim and Mackey, 2003) and the MAPK cascade (Markevich et al., 2004; Chickarmane et al., 2007), to name but a few. When considering models with many parameters, one faces the difficulty that in classical bifurcation analysis, only one parameter at a time can be varied. Thus the effect of simultaneous variations in several parameters can not be evaluated properly.

In this paper, we present a new approach to locate bifurcations in systems with feedback loops containing many parameters which may be varied simultaneously. To this end, we make use of an appropriate frequency domain description of the system and of mathematical conditions representing necessary conditions for a bifurcation. The paper is structured as follows. In Section 2, we present the theoretical results required for our approach and suggest an optimization–based method to actually find interesting parameter values. In Section 3, we apply these results to a model of the MAPK cascade with a negative feedback circuit (Kholodenko, 2000). The relevance of our results is discussed in Section 4. The mathematical proofs of the theoretical results are not presented in this paper but will be provided elsewhere (Waldherr and Allgöwer, in preparation).

2 Theoretical background

2.1 The loop–breaking approach

We consider a nonlinear differential equation which may describe the biochemical network constituting a cellular signaling pathway,


with and , where is a connected subset of . Typically, will represent the concentrations of the signaling molecules, and collects parameters like reaction constants or enzyme concentrations. We assume that an equilibrium point exists for all parameter values and can be computed at least numerically, such that for all .

Mathematically, the system (1) is said to contain a feedback circuit if the influence graph of its Jacobian contains a nontrivial loop (Cinquin and Demongeot, 2002). We want to study the influence of such a feedback circuit on the dynamical properties of the system. Control theory provides efficient tools to study this problem. A useful approach in our setup is to consider the system (1) as a closed loop control system. It is then possible to study the corresponding open loop system, and one can resort to the rich stability theory developped for control systems.

An open loop control system corresponding to the closed loop system (1) is obtained by loop breaking, as defined in the following.

Definition 1

A loop breaking for the system (1) is a pair , where is a smooth vector field and is a smooth function, such that


The corresponding open loop system is then given by the equation


The closed loop system can again be obtained by “closing the loop”, i.e. setting . Notice that by the assumption that an equilibrium exists for the closed loop system, the open loop system also has the equilibrium when choosing . This input is denoted as .

Since our main interest is in stability properties of the equilibrium point , we can restrict the analysis to the linear approximation of the systems (1) and (3) around the equilibrium point. By using Laplace transformation, the linear approximation of the open loop system (3) can be represented by a linear parameter-dependent transfer function


where and are polynomials in the complex variable with coefficients depending on . As a technical restriction, we assume that the open loop system has no poles or zeros on the imaginary axis, i.e.  and are assumed to have no roots on the imaginary axis for any value of throughout this section.

2.2 Properties of the closed and open loop systems

Stability of an equilibrium point of the closed loop system depends on the position of the eigenvalues of the Jacobian . To characterize these eigenvalues from conditions on the open loop system, we have the following theorem.

Theorem 1

Let and . Assume that is not an eigenvalue of . Then is an eigenvalue of , if and only if .

The proof of Theorem 1 is based on a representation of as

Parameter values on the border of stability are characterised by the matrix having eigenvalues on the imaginary axis. To study the corresponding property in the frequency domain representation of the open loop system, we introduce the notation of critical frequencies and gains.

Definition 2

We say that is a critical frequency and a corresponding critical gain for the transfer function (4), if


In general, different critical frequencies and gains will be obtained for different values of .

The critical frequencies can be characterized independently of the critical gains. This result follows from (5), because the transfer function value at the critical frequency has to be a real number.

Proposition 1

is a critical frequency for , if and only if


There exists a unique corresponding critical gain for any critical frequency , which is given by


The equation (6) is a polynomial in , thus all critical frequencies can be computed numerically for fixed parameters . The set of all critical frequencies for the transfer function is given by


The concept of critical frequencies and critical gains can be understood intuitively when considering the Nyquist plot of the transfer function . A critical frequency is any value at which the Nyquist plot crosses the real axis. The corresponding critical gain is the value which scales the Nyquist plot in such a way that the crossing point at the critical frequency is mapped to in the complex plane. However, this intuitive way of scaling the Nyquist plot would require to keep all critical frequencies in constant when varying parameters, which would be a strong restriction. The next section presents an approach to overcome this restriction.

2.3 A minimal set of critical frequencies

The number of critical frequencies that exist for a given open loop system is often predefined by the position of the open loop poles and zeros in the left or right half complex plane. The following proposition guarantees the existence of a minimal number of critical frequencies.

Proposition 2

Let , where () is the number of poles of in the right (left) half complex plane and () is the number of zeros of in the right (left) half complex plane. Then, for any , has at least distinct elements, if is odd, and at least distinct elements, if is even.

Since we assumed that has no poles or zeros on the imaginary axis, the number is the same for all parameters . Thus, it can be used to characterise a set of critical frequencies as being minimal.

Definition 3

Under the assumptions of Prop. 2, the set of critical frequencies is called minimal, if it contains exactly the minimal number of elements according to Prop. 2.

If is minimal, we can label the roots of (6) in a consistent way, and write , where the can be identified with different solution branches of the polynomial equation (6).

2.4 Existence of critical parameter values

The concept of critical frequencies and gains is now applied to the problem of how stability depends on parameters. We study the problem of finding critical parameters on the border of stability, i.e. such that the eigenvalues of the Jacobian are located on the imaginary axis. Then there exist typically parameters and in a neighborhood of such that the equilibrium is stable and is unstable.

The following theorem uses the loop–breaking approach and the concept of critical frequencies to characterise the existence of critical parameters.

Theorem 2

Assume that is minimal for all . Then there exists such that is an eigenvalue of , if and only if there exist such that and , for and for any .

Thus, instead of having to look at how the eigenvalues of the closed–loop system change with parameters, we have reduced this to one number, given by , which contains all information about whether the system changes its stability properties when changing parameters. The result is global in the sense that the parameters and can be arbitrarily far apart from each other, still under the given conditions existence of critical parameters is guaranteed.

2.5 Searching for critical parameter values

The theoretical approach outlined above can be used to search for parameter values such that the equilibrium point of the closed loop system (1) changes its stability. For a biochemical system, there are often nominal parameters , giving rise to the equilibrium point . We want to find parameters such that and have different stability properties.

In view of the methods presented in this paper, given the open loop transfer function (4), one first needs to identify the critical frequency that is to be considered. This choice depends on the type of stability change one is looking for. When taking , it is possible to search for zero eigenvalues, and if , nonzero imaginary eigenvalues may be encountered, typically giving rise to a Hopf bifurcation in the closed loop system. The nominal transfer function value at the critical frequency is . To change stability properties, we will then define a value as either , if , or as otherwise. Then, any solution to the nonlinear equation


gives parameters such that and have different stability properties as indicated by the chosen critical frequency . This method has been implemented using a nonlinear constrained optimization algorithm from the Matlab Optimization Toolbox (The MathWorks Inc., 2006). It allows to efficiently compute parameter values for the desired transfer function value for medium sized systems, as the example presented in the following section illustrates.

Once a parameter is known, we can use a straight line going from to , defined as . The change in dynamical behaviour along this line can then be studied using classical bifurcation analysis with respect to , implemented usually via continuation methods (Kuznetsov, 1995). In this study, the software AUTO (Doedel et al., 2006) has been used for the bifurcation analysis along the parameter line .

3 Application to a MAPK signaling module

3.1 Model description

The method presented in Section 2 has been applied to an ODE model of a mitogen activated protein kinase (MAPK) signaling module. MAPK signaling is a recurring motif in cellular signaling pathways, and typically appears in a cascade involving three levels (Pearson et al., 2001).

For this study, we consider a mathematical model for the Ras/Raf signaling pathway similar to the one presented by Kholodenko (2000). The inhibition of the upstream molecule SOS by activated MAPK, the lowest level in the cascade, constitutes a negative feedback circuit around the cascade. Via the loop–breaking approach, the influence of this feedback connection on existence of sustained oscillations in kinase activity is analysed.

The structure of the model is illustrated in Fig. 1. The reaction rates as labeled in the figure are displayed in Table 1. The concentrations have been denoted as [MAPKKK*], [MAPKK*], [MAPKK**], [MAPK*] and [MAPK**]. The concentrations of unphosphorylated kinases can be computed by conservation laws and the three parameters , and for the total concentrations of the three kinases. The difference to the model from Kholodenko (2000) is that the phosphorylation reactions 3, 4, 7 and 8 are assumed to follow mass action rather than Michaelis-Menten kinetics. This is reasonable since the Michaelis-Menten kinetics assumes low enzyme concentration compared to the substrate, whereas the concentrations of the kinases are in a comparable range here. Nominal parameter values have been adopted from Kholodenko (2000), and are shown in Table 2 as .


(.5cm) \rput[bl](0,0) \pssetframesep=1pt \rput(1.4,3.9)\pscirclebox1 \rput(1.9,2.8)\pscirclebox2 \rput(1.5,2.3)\pscirclebox3 \rput(4.0,2.3)\pscirclebox4 \rput(3.9,1.3)\pscirclebox5 \rput(1.7,1.3)\pscirclebox6 \rput(3.5,0.8)\pscirclebox7 \rput(6.0,0.8)\pscirclebox8 \rput(5.8,-0.2)\pscirclebox9 \rput(3.7,-0.2)\pscirclebox10

Figure 1: Illustration of the MAPK cascade model.
Reaction Rate
Table 1: Reaction rates in the MAPK cascade model

Using the reaction rates from Table 1, the model can be written as a system of five ODEs with 20 parameters:


For the nominal parameters , the model has a stable equilibrium . Solutions of the model converge quickly to the steady state, as depicted in Fig. 2.

Figure 2: Convergence to steady state for parameters (dashed line) and sustained oscillations for parameters (solid line). The oscillations coexist with an unstable equilibrium (dotted line).

3.2 Parameters for a change in stability properties

This section describes the application of the method presented in Section 2 to the problem of finding destabilizing parameters for the MAPK cascade model (10).

The first step is to choose a suitable loop breaking. For the MAPK cascade, an intuitive approach is to break the loop at the feedback inhibition of reaction by MAPK**. Thus we choose , to select [MAPK**] as an output, and replace by the input in the reaction rate to obtain the dynamics of the open loop system .

It can be shown that there is a unique equilibrium of (10) for any parameters in the biologically meaningful range. The equilibrium can easily be computed numerically. A linearisation of the open loop system around this equilibrium point and a Laplace transformation gives the transfer function , whose graph is shown in Fig. 3. The set of critical frequencies is minimal with , which can be seen from Fig. 3 by the observation that the graph of encircles the origin monotonically. The only positive critical frequency is , and we will consider this frequency in the search for destabilizing parameters. The corresponding transfer function value is , corresponding to the equilibrium being stable in the closed loop system.

For the computational approach described in Section 2.5, we chose , such that the value of the transfer function would have to pass the point 1 when going from its inital value of to . The optimization method converges to the parameters , which give the desired value at a critical frequency . The parameter values in are shown in Table 2. The maximal single parameter change from to has been restricted in the numerical implementation to be not more than a factor of . Even with this restriction, parameters leading to sustained oscillations have been found. However, 11 out of the 20 parameters have been changed by more than 20 % to achieve this.

Param. Unit rel. change
2.5 2.4 nM/s
9 10.6 nM
10 9.4 nM
0.25 0.11 nM/s
8 1.6 nM
0.001 0.0026 1/(s nM)
0.001 1/(s nM)
0.75 0.32 nM/s
15 3.9 nM
0.75 3.7 nM/s
15 13.3 nM
0.001 0.0033 1/(s nM)
0.001 1/(s nM)
0.5 0.26 nM/s
15 14.9 nM
0.5 2.5 nM/s
15 15.0 nM
100 100.0 nM
300 300.8 nM
300 304.2 nM
Table 2: Reference parameters and parameters for instability in the MAPK cascade model.

The graph of is shown in Figure 3. For the new parameters , the graph now encircles the point 1. By the argument principle, we see that the linearisation of the closed–loop system around the equilibrium has some eigenvalues in the right half complex plane and is thus unstable. The sustained oscillations that appear in this case are shown in Fig. 2.


(8,2.6) \rput(2,1.3) \rput(6,1.3) \pssetlinewidth=0.006 \psframe(2.65,1.2)(3.3,1.62) \pssetlinewidth=0.001,linestyle=dashed \psline(2.65,1.2)(4.72,0.24) \psline(2.65,1.62)(4.71,2.6) \psline(3.3,1.2)(7.73,0.24) \psline(3.3,1.62)(7.73,2.6)

Figure 3: Nyquist plots of open–loop MAPK model for parameters (dashed line) and (solid line).

In conclusion, our method is able to compute parameters which render the stable equilibrium unstable and thus lead to the emergence of sustained oscillations. About half of the parameters are varied by a non–negligible amount, but all variations are within the physiological range.

3.3 Bifurcation analysis along a line

Let us now consider the line . By classical bifurcation analysis with as bifurcation parameter, we can see how the system changes from the stable to the unstable equilibrium. The resulting bifurcation diagram is shown in Figure 4. As expected, there is a Hopf bifurcation between and , at . The evolution of the limit cycle producing the sustained oscillations along the line in parameter space is obtained from the bifurcation diagram.

Figure 4: Bifurcation diagram along the line , showing stable equilibrium (solid line), unstable equilibrium (dashed line) and amplitude of oscillations (circles).

4 Conclusions

We introduced some theoretical tools to investigate the existence of parameters for which a bifurcation can occur in a dynamical system with a feedback circuit. These tools gave rise to a new computational method which allows to search for parameter values such that the stability properties of an equilibrium change in a specific way compared to the nominal parameter values. Our approach is particularly useful if there are many parameters in the system which can be varied simultaneously, and if the contribution of individual parameters to stability properties is not obvious. The ability to directly handle multiparametric variations is a clear advantage compared to using only classical bifurcation analysis.

We have shown the application of the proposed method to a model of a MAPK cascade. Using relatively small changes to most of the 20 parameters in the model leads to a change from a stable equilibrium to an unstable equilibrium with a stable limit cycle, producing sustained oscillations.


  1. thanks: [


  1. D. Angeli and E. D. Sontag. Interconnections of monotone systems with steady-state characteristics. In M. de Queiroz, M. Malisoff, and P. Wolenski, editors, Optimal control, stabilization, and nonsmooth analysis, pages 135–154. Springer-Verlag, 2004.
  2. V. Chickarmane, B. N. Kholodenko, and H. M. Sauro. Oscillatory dynamics arising from competitive inhibition and multisite phosphorylation. J. Theor. Biol., 244(1):68–76, January 2007. URL http://dx.doi.org/10.1016/j.jtbi.2006.05.013.
  3. O. Cinquin and J. Demongeot. Positive and negative feedback: striking a balance between necessary antagonists. J. Theor. Biol., 216(2):229–241, May 2002. doi: 10.1006/jtbi.2002.2544. URL http://dx.doi.org/10.1006/jtbi.2002.2544.
  4. B. F. Dibrov, A. M. Zhabotinsky, and B. N. Kholodenko. Dynamic stability of steady states and static stabilization in unbranched metabolic pathways. J. Math. Biol., 15(1):51–63, 1982.
  5. E. J. Doedel, R. C. Paffenroth, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. E. Oldeman, B. Sandstede, and X. Wang. AUTO 2000: continuation and bifurcation software for ordinary differential equations. Concordia University, Montreal, Canada, 2006.
  6. T. Eißing, S. Waldherr, F. Allgöwer, P. Scheurich, and E. Bullinger. Steady state and (bi-) stability evaluation of simple protease signalling networks. BioSystems, Epub ahead of print, 2007. URL http://dx.doi.org/10.1016/j.biosystems.2007.01.003.
  7. M. Kaufman and R. Thomas. Emergence of complex behaviour from simple circuit structures. Comptes rend. biol., 326:205–214, 2003.
  8. B. N. Kholodenko. Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur. J. Biochem., 267(6):1583–88, Mar 2000.
  9. J. Kim, D. G. Bates, I. Postlethwaite, L. Ma, and P. A. Iglesias. Robustness analysis of biochemical network models. IEE Proc. Syst. Biol., 153(3):96–104, May 2006.
  10. Y. A. Kuznetsov. Elements of Applied Bifurcation Theory. Springer-Verlag New York, 1995.
  11. N. I. Markevich, J. B. Hoek, and B. N. Kholodenko. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J. Cell Biol., 164(3):353–359, Feb 2004. doi: 10.1083/jcb.200308060. URL http://dx.doi.org/10.1083/jcb.200308060.
  12. Optimization toolbox. For use with Matlab. The MathWorks Inc., 2006.
  13. G. Pearson, F. Robinson, T. B. Gibson, B. E. Xu, M. Karandikar, K. Berman, and M. H. Cobb. Mitogen-activated protein (MAP) kinase pathways: regulation and physiological functions. Endocr. Rev., 22(2):153–183, Apr 2001.
  14. C. D. Thron. The secant condition for instability in biochemical feedback-control. 1. The role of cooperativity and saturability. Bull. Math. Biol., 53(3):383–401, 1991.
  15. J. J. Tyson and H. G. Othmer. The dynamics of feedback control circuits in biochemical pathways. Progr. Theor. Biol., 5:2–62, 1978.
  16. N. Yildirim and M. C. Mackey. Feedback regulation in the lactose operon: a mathematical modeling study and comparison with experimental data. Biophys. J., 84(5):2841–51, May 2003.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description