Skip Navigation


Bioinformatics Advance Access originally published online on August 4, 2005
Bioinformatics 2005 21(18):3688-3690; doi:10.1093/bioinformatics/bti603
This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (Print PDF) Freely available
Right arrow All Versions of this Article:
21/18/3688    most recent
bti603v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Right arrow Search for citing articles in:
ISI Web of Science (8)
Right arrowRequest Permissions
Google Scholar
Right arrow Articles by Chickarmane, V.
Right arrow Articles by Sauro, H. M.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Chickarmane, V.
Right arrow Articles by Sauro, H. M.
Social Bookmarking
Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?

© The Author 2005. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions{at}oupjournals.org

Bifurcation discovery tool

Vijay Chickarmane 1,*, Sri R. Paladugu 1, Frank Bergmann 1 and Herbert M. Sauro 1,2

1Keck Graduate Institute 535 Watson Dr, Claremont, CA 91711, USA
2Control and Dynamical Systems, California Institute of Technology Pasadena, CA 91125

*To whom correspondence should be addressed.


   Abstract
TOP
 Abstract
1 INTRODUCTION
 2 APPROACH
 3 EXAMPLES
 REFERENCES
 

Motivation: Biochemical networks often yield interesting behavior such as switching, oscillation and chaotic dynamics. This article describes a tool that is capable of searching for bifurcation points in arbitrary ODE-based reaction networks by directing the user to regions in the parameter space, where such interesting dynamical behavior can be observed.

Results: We have implemented a genetic algorithm that searches for Hopf bifurcations, turning points and bistable switches. The software is implemented as a Systems Biology Workbench (SBW) enabled module and accepts the standard SBML model format. The interface permits a user to choose the parameters to be searched, admissible parameter ranges, and the nature of the bifurcation to be sought. The tool will return the parameter values for the model for which the particular behavior is observed.

Availability: The software, tutorial manual and test models are available for download at the following website: http:/www.sys-bio.org/ under the bifurcation link. The software is an open source and licensed under BSD.

Contact: vijay_chickarmane{at}kgi.edu


   1 INTRODUCTION
TOP
 Abstract
 1 INTRODUCTION
2 APPROACH
 3 EXAMPLES
 REFERENCES
 
Biological systems in the form of chemical reaction networks are known to display non-trivial qualitative behaviors. Examples include switching, oscillation (Fall et al., 2004) and chaotic dynamics (Goldbeter et al., 2001). One of the issues in model building is understanding the scope of behaviors that the model can elicit. Often the parameter space in these models is large and searching for potentially interesting dynamics in this space is difficult. The onset of qualitative changes in a model occurs at points in the parameter space called bifurcation points. These points, which mark the transition between one behavior and another, are characterized by particular changes in the eigenvalues that describe the local dynamics of the model. A possible approach to locating interesting qualitative properties in a model is to use an optimization method to search for bifurcation points. This paper describes an algorithm that is capable of searching for interesting qualitative behavior in biochemical reaction models.


   2 APPROACH
TOP
 Abstract
 1 INTRODUCTION
 2 APPROACH
3 EXAMPLES
 REFERENCES
 
The essence of the technique is to minimize an objective function that is a test function for a particular bifurcation type (Strogatz, 1994; Kuznetsov, 1995). In the current implementation we have provided tests for turning points, Hopf bifurcation points and bistable switches; these are arguably the most common bifurcations in biochemical networks. Given the size of the parameter search space we chose to implement the minimization algorithm using a standard Genetic Algorithm (GA) (Goldberg, 1989). The search begins by creating a population of networks having the desired topology but with randomly generated parameter values. Therefore, each network represents a point in the parameter space. After that the fitness value for each of these networks are evaluated and the networks are ranked accordingly. The algorithm uses tournament selection to select the mating pairs of networks and then employs arithmetic crossover and mutation operators to create a new population. The mutation operator involves randomly selecting a parameter and multiplying its current value by a factor between 0.5 and 1.5 (Herrara et al., 1998). Finally, elitism is used to guarantee that the best network is carried into the next generation. The search is continued until the desired fitness value is reached or the total number of generations reaches the maximum.

2.1 Objective function
Consider a dynamical system described by N nonlinear, coupled, ordinary differential equations, whose state is described by N variables Xi, i = 1:N, and by the parameters pk, k = 1:m,

where fj (Xi, pk) are, in general, nonlinear functions of all the Xi, pk. The eigenvalues of the Jacobian ({partial}fi/{partial}Xj) determine the stability of the system evaluated at a steady state. If {lambda}i = {lambda}Ri + i * {lambda}Ii are the complex eigenvalues of the Jacobian for a particular steady state, then instability will result if any of . Hence the interesting behavior occurs when the real part of the eigenvalues cross the imaginary axis from the negative side. The algorithm is restricted to co-dimension one bifurcations, and hence if the cross over occurs for one eigenvalue we get either a turning point or a static bifurcation point (transcritical or pitchfork). If, however, a complex pair of eigenvalues become purely imaginary, then we get a Hopf bifurcation, and the system will exhibit limit cycle behavior. In the following tests, all parameters in the system are examined for possible bifurcation behaviour.

2.1.1 Detecting switch-like behavior
The condition for a turning point bifurcation with respect to one of the parameters, pi, is that one of the eigenvalues should be zero. A simple function to minimize is therefore given by {varepsilon}, Equation (1), where the numerator is the product of the eigenvalues {lambda}i and {lambda}Min refers to the product of all the eigenvalues except the minimum eigenvalue1.

(1)
The Jacobian is also singular for pitchfork and transcritical bifurcations. A genuine turning point, however, has the additional property that the rank of the matrix, formed by replacing one of the columns of the Jacobian with the column vector {partial}fi/{partial}pj, is N. When the algorithm signals a location in parameter space, where the system is thought to be at a turning point, we test the above condition, and only if it is successful do we proceed to test whether the system exhibits bistability. This last step involves a simple continuation of the steady states of one of the species, which is computed for each parameter, to identify the parameter(s) for which one obtains switch-like behavior. Although we obtain bistability for some of these networks, for the majority we obtain ultrasensitive (memory-less switches) curves.

2.1.2 Detecting Hopf bifurcations
For the Hopf bifurcation, we attempt to minimize the following function:

(2)
The numerator of the objective function simply states that at a Hopf bifurcation the real part of one of a pair of complex eigenvalues goes to zero. This condition is, however, unable to differentiate between systems that have a complex conjugate pair of eigenvalues approaching the imaginary axis and a real eigenvalue moving towards the imaginary axis. The latter is not a Hopf bifurcation but is either a turning point or a static bifurcation point. We therefore differentiate between these two alternatives by inserting a denominator in the objective function Equation (2), which rewards a network with eigenvalues that have imaginary parts and whose real parts go to zero, and punishes those networks that only have a real eigenvalue going to zero. Once the tool has located a putative bifurcation point, it is necessary to test for which parameter the eigenvalues cross the imaginary axis. For these sets of parameters we declare them as the bifurcation parameters. The tool can either pass the parameter set to a model editor for simulation or the model can be passed to a bifurcation package such as AUTO (Doedel et al., 1991).

A number of other tools are being developed for discovering bifurcation points. For example Oscill8 is a GUI wrapper around AUTO (http://oscill8.sourceforge.net/) and can be used for bifurcation analysis. Gepasi (Mendes, 1993) can also be used to carry out limited bifurcation discovery by using its built-in optimization capability.


   3 EXAMPLES
TOP
 Abstract
 1 INTRODUCTION
 2 APPROACH
 3 EXAMPLES
REFERENCES
 
We have tested the discovery tool on a variety of example models that are listed in Table 1. Figure 1 shows the progress of the fitness in a typical optimization run for locating a Hopf bifurcation. In Figure 2 we plot the results obtained from two typical search runs. The largest model (unpublished) that we have tested on has 25 state variables and 37 parameters. For the models that were tested we find that the major factor that affects the speed of the search process is the inherent stiffness of the model. Currently, we are developing an advanced simulator using CVODE to deal with stiff systems.


View this table:
[in this window]
[in a new window]
  		Table 1 Test models

 


View larger version (13K):
[in this window]
[in a new window]
  		Fig. 1 Fitness score versus iteration for activator–inhibitor oscillator [Tyson et al., 2003 (d)]. The subplot shows the fitness graph magnified from iteration 8 to 15, reaching a level where a Hopf bifurcation is detected.

 


View larger version (32K):
[in this window]
[in a new window]
  		Fig. 2 Locating bifurcation points in a relaxation oscillator model (Seno et al., 1978). In the top left subplot, the time-series for one of the variables of the model shows oscillations after the parameters have been optimized. The next subplot shows the bifurcation diagram for the optimized parameters, exhibiting two supercritical Hopf bifurcations. The bottom two plots are for the bistable cdc2/wee1 model (Angeli et al., 2004), which shows the time-series of the variables exhibiting switch like behavior, depending on their initial conditions. The bifurcation diagram shows hysteresis.

 


   Acknowledgments
 
V.C. would like to acknowledge Professor John Tyson for fruitful discussions and encouragement. The authors thank Anastasia Deckard for a critical review of the manuscript. H.M.S., V.C. and F.B. are grateful to generous support from the DARPA/IPTO BioCOMP program, contract number MIPR 03-M296-01. SRP is supported by a grant from the DOE GTL Program.

Conflict of Interest: none declared.


   Footnotes
 
1In the above equation, the denominator was added to prevent the other eigenvalues from becoming very small, which usually occurs if the objective function is the determinant of the Jacobian. Hence networks that have small {Pi}{lambda}Min (product of the the eigenvalues leaving out the minimum) are punished, and this prevents all the eigenvalues from moving towards the imaginary axis. Back

Received on April 8, 2005; revised on July 1, 2005; accepted on July 27, 2005

   REFERENCES
TOP
 Abstract
 1 INTRODUCTION
 2 APPROACH
 3 EXAMPLES
 REFERENCES
 

    Angeli, D., et al. (2004) Detection of multistability, bifurcations and hysteresis in a large class of biological positive-feedback systems. Proc. Natl Acad. Sci. USA, 101, 1822–1827[Abstract/Free Full Text].

    Doedel, E.J., et al. (1991) Numerical analysis and control of bifurcation problems. Part I. Int. J. Bifurcat. Chaos, 1, 493–520[CrossRef].

    Edelstein, B.B. (1970) Biochemical model with multiple steady states and hysteresis. J. Theor. Biol., 29, 57–62[CrossRef][ISI][Medline].

    Fall, C.P., et al. Computational Cell Biology, (2004) , Springer-Verlag New York.

    Mendes, P. (1993) GEPASI: A software package for modelling the dynamics, steady states and control of biochemical and other systems. Comput. Appl. Biosci., 9, 563–571[Abstract/Free Full Text].

    Goldbeter, A., et al. (2001) From simple to complex oscillatory behavior in metabolic and genetic control networks. Chaos, 11, 247–260[CrossRef][ISI][Medline].

    Goldberg, D.E. Genetic Algorithms in Search, Optimization and Machine Learning, (1989) , Addison-Wesley.

    Heinrich, R., et al. (1977) Metabolic regulation and mathematical models. Prog. Biophys. Mol. Biol., 32, 1–82[ISI][Medline].

    Herrara, F., et al. (1998) Tackling real-coded GA's:operator and tools for behavioural analysis. Art. Intel. Rev., 12, 265–319.

    Hervagault, J.F. and Canu, S. (1987) Bistability and irreversible transitions in a simple substrate cycle. J. Theor. Biol., 127, 439–449[CrossRef][ISI][Medline].

    Kholodenko, B.N. (2000) Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur. J. Biochem., 267, 1583–1588[ISI][Medline].

    Kuznetsov, Y.A. Elements of Applied Bifurcation Theory, (1995) , Springer-Verlag New York.

    Nicolis, G. and Prigogine, I. Self-Organization in Non-Equilibrium Systems, (1977) , Wiley-Interscience New York.

    Sauro, H., et al. (2003) Next generation simulation tools: the Systems Biology Workbench and BioSPICE integration. OMICS, 7, 355–372[CrossRef][Medline].

    Seno, M., et al. (1978) Instability and oscillatory behavior of membrane-chemical reaction systems. J. Theor. Biol., 72, 577–588[CrossRef][ISI][Medline].

    Strogatz, S.H. Nonlinear Dynamics and Chaos, (1994) Westview Press.

    Tyson, J.J., et al. (2003) Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr. Opin. Cell. Biol., 15, 221–231[CrossRef][ISI][Medline].


Add to CiteULike CiteULike   Add to Connotea Connotea   Add to Del.icio.us Del.icio.us    What's this?


This article has been cited by other articles:


Home page
 BioinformaticsHome page
G. Rodrigo, J. Carrera, and A. Jaramillo
Genetdes: automatic design of transcriptional networks
Bioinformatics, July 15, 2007; 23(14): 1857 - 1858.
[Abstract] [Full Text] [PDF]

Home page
 BioinformaticsHome page
R. R. Vallabhajosyula, V. Chickarmane, and H. M. Sauro
Conservation analysis of large biochemical networks
Bioinformatics, February 1, 2006; 22(3): 346 - 353.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (Print PDF) Freely available
Right arrow All Versions of this Article:
21/18/3688    most recent
bti603v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Right arrow Search for citing articles in:
ISI Web of Science (8)
Right arrowRequest Permissions
Google Scholar
Right arrow Articles by Chickarmane, V.
Right arrow Articles by Sauro, H. M.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Chickarmane, V.
Right arrow Articles by Sauro, H. M.
Social Bookmarking
Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?