## 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 *X*_{i}, *i* = 1:*N*, and by the parameters *p*_{k}, *k* = 1:*m*,

*f*

_{j}(

*X*

_{i},

*p*

_{k}) are, in general, nonlinear functions of all the

*X*

_{i},

*p*

_{k}. The eigenvalues of the Jacobian (∂

*f*

_{i}/∂

*X*

_{j}) determine the stability of the system evaluated at a steady state. If λ

_{i}= λ

^{R}

_{i}+

*i** λ

^{I}

_{i}are the complex eigenvalues of the Jacobian for a particular steady state, then instability will result if any of

*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, *p*_{i}, 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.

**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 ∂

*f*

_{i}/∂

*p*

_{j}, 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:

*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.

^{1}In 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**

**Fig. 1**

**Fig. 2**

**Fig. 2**

**Table 1**

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 | 6 |

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) | 5 |

Heinrich et al., 1977 | Positive feedback oscillator | 4 |

Seno et al., 1978 | Modified edelstein relaxation oscill | 5 |

Kholodenko, 2000 | MAPK feedback oscillator | 8 |

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 | 6 |

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) | 5 |

Heinrich et al., 1977 | Positive feedback oscillator | 4 |

Seno et al., 1978 | Modified edelstein relaxation oscill | 5 |

Kholodenko, 2000 | MAPK feedback oscillator | 8 |

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

^{a}Approximate 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.

## Comments