## Abstract

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@kgi.edu

## 1 INTRODUCTION

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

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,

$\frac{\hbox{ d }{X}_{j}}{\hbox{ d }t}={f}_{j}\left({X}_{i},{p}_{k}\right),$
where fj (Xi, pk) are, in general, nonlinear functions of all the Xi, pk. The eigenvalues of the Jacobian (∂fi/∂Xj) determine the stability of the system evaluated at a steady state. If λi = λRi + i * λIi are the complex eigenvalues of the Jacobian for a particular steady state, then instability will result if any of
$${\lambda }_{i}^{R} > 0$$
. 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 ε, Equation (1), where the numerator is the product of the eigenvalues λi and ∏λMin refers to the product of all the eigenvalues except the minimum eigenvalue1.

(1)
$\varepsilon =\frac{\Pi {\lambda }_{i}}{\left(1-0.99\times exp\left(-\left|\Pi {\lambda }_{\hbox{ Min }}\right|\right)\right)}.$
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 ∂fi/∂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)
$\varepsilon =\frac{\Pi {{\lambda }^{R}}_{i}}{\Pi \left(1-0.99\times exp\left(-\left|{\lambda }_{i}^{I}\right|\right)\right)}.$
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

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.

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 ΠλMin (product of the the eigenvalues leaving out the minimum) are punished, and this prevents all the eigenvalues from moving towards the imaginary axis.

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.

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.

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.

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.

Table 1

Test models

Model Type Approximate no. of iterationsa
Tyson et al., 2003 (a) Irreversible bistable switch 45
Tyson et al., 2003 (b) Reversible bistable switch 36
Edelstein, 1970 Autocatalytic bistable switch 41
Hervagault and Canu, 1987 Bistable switch 54
Angeli et al., 2004 Bistable switch cdc2/Wee1
Tyson et al., 2003 (c) Negative feedback oscillator 20
Tyson et al., 2003 (d) Activator-inhibitor oscillator 15
Tyson et al., 2003 (e) Substrate-depletion oscillator 21
Nicolis and Prigogine, 1977 Autocatalytic oscillator (brusselator)
Heinrich et al., 1977 Positive feedback oscillator
Seno et al., 1978 Modified edelstein relaxation oscill
Kholodenko, 2000 MAPK feedback oscillator
Model Type Approximate no. of iterationsa
Tyson et al., 2003 (a) Irreversible bistable switch 45
Tyson et al., 2003 (b) Reversible bistable switch 36
Edelstein, 1970 Autocatalytic bistable switch 41
Hervagault and Canu, 1987 Bistable switch 54
Angeli et al., 2004 Bistable switch cdc2/Wee1
Tyson et al., 2003 (c) Negative feedback oscillator 20
Tyson et al., 2003 (d) Activator-inhibitor oscillator 15
Tyson et al., 2003 (e) Substrate-depletion oscillator 21
Nicolis and Prigogine, 1977 Autocatalytic oscillator (brusselator)
Heinrich et al., 1977 Positive feedback oscillator
Seno et al., 1978 Modified edelstein relaxation oscill
Kholodenko, 2000 MAPK feedback oscillator

(Available at www.sys-bio.org).

aApproximate number of iterations taken to converge to the solution with random initial conditions using the default settings in the tool.

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.

## REFERENCES

Angeli, D., et al.
2004
Detection of multistability, bifurcations and hysteresis in a large class of biological positive-feedback systems.

101
1822
–1827
Doedel, E.J., et al.
1991
Numerical analysis and control of bifurcation problems. Part I.
Int. J. Bifurcat. Chaos

1
493
–520
Edelstein, B.B.
1970
Biochemical model with multiple steady states and hysteresis.
J. Theor. Biol.

29
57
–62
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
Goldbeter, A., et al.
2001
From simple to complex oscillatory behavior in metabolic and genetic control networks.
Chaos

11
247
–260
Goldberg, D.E.
Genetic Algorithms in Search, Optimization and Machine Learning

1989
Heinrich, R., et al.
1977
Metabolic regulation and mathematical models.
Prog. Biophys. Mol. Biol.

32
1
–82
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
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
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
Seno, M., et al.
1978
Instability and oscillatory behavior of membrane-chemical reaction systems.
J. Theor. Biol.

72
577
–588
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