Inverse bifurcation analysis: application to simple gene systems
© Lu et al; licensee BioMed Central Ltd. 2006
Received: 21 June 2006
Accepted: 21 July 2006
Published: 21 July 2006
Bifurcation analysis has proven to be a powerful method for understanding the qualitative behavior of gene regulatory networks. In addition to the more traditional forward problem of determining the mapping from parameter space to the space of model behavior, the inverse problem of determining model parameters to result in certain desired properties of the bifurcation diagram provides an attractive methodology for addressing important biological problems. These include understanding how the robustness of qualitative behavior arises from system design as well as providing a way to engineer biological networks with qualitative properties.
We demonstrate that certain inverse bifurcation problems of biological interest may be cast as optimization problems involving minimal distances of reference parameter sets to bifurcation manifolds. This formulation allows for an iterative solution procedure based on performing a sequence of eigen-system computations and one-parameter continuations of solutions, the latter being a standard capability in existing numerical bifurcation software. As applications of the proposed method, we show that the problem of maximizing regions of a given qualitative behavior as well as the reverse engineering of bistable gene switches can be modelled and efficiently solved.
The use of mathematical models provides tools for the analysis of complex molecular interactions aiming at an understanding of processes occurring in living cells. For many problems in cellular control, stochastic effects and time-delays can be ignored and systems of first-order ordinary differential equations (ODEs) can adequately model the underlying processes. Denoting by x and p the biochemical concentrations and parameters, respectively, the instantaneous change in x is described by the vector field f:
In the study of such systems, an important goal is to understand how the observed physiological behavior arises out of gene network topology and parameters p. Some of these questions may be studied via examining the bifurcation manifolds Σ of the ODE system, which partition the parameter space into regions of different qualitative behavior (see e.g.,  for a general overview to bifurcation theory). From ODE models and measured parameters, the forward problem of computing the bifurcation diagram has contributed significantly towards elucidating the complex mechanisms underlying cellular processes. For instance, mathematical and symbolic bifurcation analysis has led to an understanding of the possible dynamical behaviors that may arise out of simple gene systems (for a monograph, see , examples of more recent papers dealing with natural, designed, and model systems are [3–7]). For cell cycle models, bifurcation diagrams have given biologists a systems-level perspective of the roles played by the various constituent modules, as well as providing the ability to predict the behavior of mutant cells [8, 9]. The desire to locate regions in parameter space exhibiting interesting dynamics has led to the development of computational tools for the automatic discovery of bifurcation sets .
In contrast, inverse problems have only recently attracted attention in biology as a way to unravel the workings of cellular mechanisms. In inverse problems one looks for causes for observed or desired effects . Mathematically, such problems are typically ill-posed, in particular unstable; special mathematical techniques, called "regularization methods", have to be used to cope with this ill-posedness. Many variational and iterative regularization techniques have been developed over the years and applied to a variety of problems in science, engineering and finance (see  and some references quoted there, and ).
In the current context of cell biology, one would like to address problems such as: which parameter configurations lead to an observed qualitative behaviour of the system ("identification")? How can one introduce a certain qualitative behaviour into the system via parameter variations ("design")? We summarize such problems under the name of inverse bifurcation problems, where the task is to map the space of bifurcation diagrams back to the space of parameters. In particular, we consider inverse bifurcation problems of two types: identification and design. For the former, one would like to infer parameter values from observed bifurcation diagrams and hence the issue of uniqueness is typically of concern. For the latter, one is interested in parameter values that produce the desired outcome, hence uniqueness is not an issue. The notion of "inverse bifurcation" was first introduced into biology in ; more recently, another inverse bifurcation problem from mathematical ecology was studied in  by integral equation methods. These papers are concerned with the existence and uniqueness of nonlinear terms in equations realizing prescribed bifurcation diagrams and use analytical methods. In contrast, given the size and complexity of gene networks, we take a computational approach in addressing problems of inverse bifurcation. In this paper, we consider problems of moderate size for which the ill-posedness is not yet a crucial issue. Since ill-posedness increases with dimension, a major issue for larger problems will be to use appropriate regularization techniques.
For the design type of inverse bifurcation problems, there exists previous work in the engineering literature: The distance to bifurcation manifolds has been introduced to quantify the "parametric robustness" of system designs in ; in the context of design of chemical processes, optimization problems with constraints involving distance to bifurcation manifolds have been treated in . For a recent review see . For biological applications as we have them in mind, this issue of parametric robustness is also important. In addition, other geometric properties of the bifurcation diagram are of interest. These include the size of the parameter region resulting in bistability of solutions and the parametric distance between regions of different qualitative behavior. We will develop methodolgies by which inverse bifurcation problems involving the optimization of such quantities can be solved.
In this paper, we show the applicability of our inverse bifurcation algorithms to low dimensional gene systems. We formulate the inverse problems as constrained optimization problems, whose objective function and constraints involve geometric properties of bifurcation diagrams. We demonstrate that these problems can be solved efficiently by applying gradient-based nonlinear optimization algorithms in combination with one-parameter continuation methods to locate bifurcation points. The latter is a standard capability provided by existing bifurcation analysis software (see  for references to state-of-the-art numerical implementations).
2 Inverse bifurcation analysis
The ODE model (1) defines a mathematical relationship between the parameters p and the time course of biochemical concentrations, x(t). Of the set of all parameters, some describe the biochemical mechanisms constituting the machinery of the regulatory system, while others are parameters to whose variation the gene system should respond. In this paper, we refer to the former as the system parameters and the latter as the control inputs.
In situations where the qualitative behavior of the regulatory system changes in response to shifts in the control input, the corresponding bifurcation diagram is of importance. In particular, the proper working of the regulatory system may depend on the locations, shapes and sizes of the various regions of qualitative behavior. For instance, in the cell cycle model  the correct spatial relationship (in parameter space) of the transitions points is crucial for ensuring the survival of cells. In fact, the model predicts that "dynamically challenged" mutants with shifted bifurcation diagrams may suffer irrecoverably. The resulting effects include difficulties exiting mitosis and decreases in cell mass after each cell division. For this particular model, the inverse bifurcation problem is to map geometric as well as topological relationships in the bifurcation diagrams to conditions on the parameters. Another application of inverse bifurcation in the cell cycle model is for the so-called "Pinocchio effect" involving check points. This effect occur when the necessary conditions for safe progression to the next cell phase are not satisfied hence the regulatory system should delay the state transition as far as possible. Here, the inverse bifurcation problem is to find out how the system may be constructed so that, when triggered by the presence or absence of certain chemicals, the range of bistability is maximized.
In the case of system parameters, they are either regulated by other control mechanisms in the cell, or the behavior of the system is constructed to be insensitive to variation in these parameters . Within an operating region in the control input space, the minimal distance in the system parameter space to bifurcation manifolds provides a quantitative measure of the robustness of the system to environmental perturbations. Here, the inverse bifurcation problem would involve finding the parameter combinations so that the system bifurcation diagram is insensitive to the unregulated parameters. Biological applications of this class of problems include homeostasis and rhythmic pacemakers. Again, the problem of interest is to map some shape property of the bifurcation diagram to conditions on the parameter space.
For inverse bifurcation problems of biological interest such as those described above, the question arises as to how to formulate them mathematically so that the solution can be obtained in a computationally tractable and stable manner. Typically, in biological applications the parameter space is of high dimension. Furthermore, in carrying out inverse analysis, the (forward) bifurcation analysis usually has to be applied a large number of times. In most cases, the former condition precludes formulations based on the use of multi-parameter continuation techniques . Instead, formulations based on continuation along rays in parameter space are preferred. Below, we formulate distance to bifurcation manifolds (first introduced in ) in terms of a forward operator and an l2 functional. Subsequently, we describe the sensitivity of the minimal distance with respect to parameters by adjoint methods.
F(p) ≡ (F(p) i , F(p) s )
The notation is illustrated in Figure 1.
Using F(p), the inverse bifurcation problems examined in this paper are mathematically formulated as:
Subject to: plow ≤ p ≤ pupp
0 ≤ c(F(p) i ), (2)
where ||·|| denotes the l2-norm and c : P i → ℝ k represents k-dimensional nonlinear constraints. In the rest of this section, we discuss adjoint sensitivity analysis of F(p), which is important for applying gradient-based optimization methods for solving (2) as well as for computing F(p) iteratively. Denote the linearization of the above map at p as F'(p). The adjoint operator F'*(p) is defined to satisfy the following: for all δp, δ ∈ ℝ m ,
where the notation <·, ·> denotes the l2-inner product. Suppose <·, l> is a linearized functional of interest on P. Given a parameter perturbation δ p, the forward functional sensitivity is then given by <F'(p)δp, l>. The same sensitivity can be obtained via the adjoint solution, ψ ≡ F'*(p)l. From the definition of the adjoint operator (3), it follows that:
<F'(p)δp, l> = <δp, ψ>. (4)
That is, for all δp ∈ P the functional sensitivities can be computed via a single linear solution for ψ rather than repeated application of the linearized forward operator, F'(p)δp.
For the case where the functional of interest is the l2-distance J(p) = ||F(p) i - p i ||, the adjoint solution for the linearization J'(p)(·) is given in terms of the vector normal to the manifold Σ. Denoting N i and N s as the components of the normal vector in P i and P s respectively, it can be shown that (up to sign),
To fix the sign of ψ, the component N i is chosen to be p i - F(p) i . Figure 2 illustrates the components of the vector N normal to the manifold Σ.
Thus, obtaining expressions for the adjoint vector for various bifurcations of interest reduces to the problem of deriving the associated normal vectors, N s . Under certain transversality conditions, the normal vectors for several codimension-one and higher bifurcations have been derived [16, 17]. Here we consider the generic codimension-one bifurcations, namely saddle-node and Hopf bifurcations. Let the left and right critical eigenvectors of f x at the given bifurcation point be denoted by w and v respectively. That is, w and v solve the following eigen-systems for the critical eigenvalue ωcrit and its conjugate crit,
f x v = ωcrit v
The expressions for normal vectors are given as:
Finally, we mention that the minimal distance functional can be used to model many other problems of interest. For instance, it can serve as an estimate of the separation between a reference manifold Σref and a given region of qualitative behavior, or the size of the region of a qualitative behavior via the maximum radius of the inscribed sphere.
3 Algorithm and software implementation
Our Inverse Bifurcation Toolbox is an implementation of the above described algorithm. It combines the capability of several packages: the MATCONT package  for performing the underlying one-parameter continuations, the MATLAB Optimization Toolbox  for performing gradient-based constrained optimization, the MathSBML package  for reading in biological models in the SBML format and finally the Mathematica Symbolic Toolbox for MATLAB  for communications. The combination of MATLAB and Mathematica has the advantage of allowing building on existing freely-available software.
After the SBML model is read in via MathSBML, the vector fields are differentiated symbolically to obtain expressions for f x , f p , f xp and f xx . Subsequently, they are numerically evaluated when called by MATCONT. Reparametrization of parameters is done automatically to allow continuation in arbitrary directions. Presently, the toolbox is able to handle bifurcations of saddle-node and Hopf type. Augmentation to include bifurcations of limit cycles is currently underway.
Given the limitations of the current implementation, our software is suitable for handling problems involving tens of variables and parameters but possibly not suited for problem dimensionality on the order of hundreds or higher. Further algorithm development is needed to significantly upscale the method. In addition, in high dimensional applications instabilities may appear and hence appropriate regularization techniques need to be developed (e.g., stopping rules in the case of iterative methods).
4.1 Finding distance to oscillation in a 4 gene model
To understand design principles underlying biological systems, simple oscillatory and switch-like systems have been constructed experimentally. Several systems based on E. coli [24, 25] have been successfully demonstrated. Recently, components of the Lac and Ntr systems have been used to construct a genetic clock that upon change in connectivity, exhibits switch behavior . After ignoring the dynamics of the read-out gene, the system of three equations for the concentrations of mRNA (x1, x3, x5) and proteins (x2, x4, x6) is given as,
where β i are rate constants. The rates of transcription are represented by the following tri-phasic functions,
where g jk are kinetic-order parameters, typically less than 4.
Given a stable, steady solution for a gene system with a set of nominal parameters, a relevant engineering question is how to construct a genetic oscillator. For the current system, the analytical solutions for the bifurcation manifolds are available via the Routh-Hurwitz stability criterion (, supplementary data). However, such analysis cannot be easily carried out for more complex models. Instead, computational methods have to be used.
0.1 ≤ β1, β3 ≤ 25
0.1 ≤ β2, β4 ≤ 4.
4 gene system
4.2 Maximizing regime of oscillations in repressilator system
Motivated by biological applications such as increasing the robustness of rhythmic pacemakers, we consider maximizing the region of oscillations with respect to a given operating point. We utilize a case of the generalized repressilator  as our test system.
After non-dimensionalizing, an ODE system for the protein and RNA concentrations x i and y i is obtained, with the following dimensionless parameters: α i , β i are the ratios of degradation rates, δ i represents the leakiness of the gene transcription, h i the Hill-coefficient reflecting the degree of cooperativity of gene transcription,
In this example, we take n = 3 and examine the symmetric case, i.e., α i = α, β = β, etc. Analysis shows that for ranges of parameters α and β, the system can exhibit oscillations or be in a stable, steady-state . Thus, we take p i = (α, β) as the input parameters and seek values of p s = (δ, h) within the constraint set,
(10-4, 0) ≤ (δ, h) ≤ (10-1, 2)
4.3 Engineering bistable switch in G1/S transition
We consider the reverse engineering of check points in mammalian cell phase transition. In particular, we use a 9 state, 40 parameter differential equation model . In the publication describing the model, parameter scanning is done to find out parameter values that achieve certain behaviors, such as changing the threshold for the onset of synthesis (S) phase from the G1 state and conditions ensuring irreversibility of the gene switch. Using our inverse bifurcation analysis, specifying a desired behavior on the transition points gives rise to conditions on the parameters in a more systematic and efficient manner.
4.3.1 Example 1: engineer irreversibility
First, we mention a number of main species and their interactions. The regulatory system consists of a core double inhibitor-activator module as well as several positive-feedback loops. The transcription factor E2F1 activates itself as well as pRB, a tumor suppressor. In turn, pRB inhibits itself as well as E2F1. The following equations describe this core module:
Another important collection of species is AP-1, which is used to denote a family of transcription factors that mediate the mitogenic signal F m . The family AP-1 is also activated by E2F1, with feedback strength k25. The dynamics for the concentration of the family [AP-1] is modelled as:
In the first example, we consider the construction of an irreversible gene switch via parameter changes. It has been demonstrated  that for some value of k25 in (8), the G1/S transition becomes irreversible. That is, if F m is increased beyond a certain point, a subsequent decrease to zero will not lead to a G0 state with a low level of E2F1.
Example 1: irreversibility
Example 2: irreversibility, fixed G1/S
4.3.2 Example 2: engineer irreversibility, fixed G1/S transition
In this example, we consider an extension of the previous problem. Here, we would like to engineer an irreversible switch with the additional constraint that the x-abscissa of the G1/S transition point SN1 should be fixed at the value of the nominal system. In addition to allowing 0.1 ≤ k25 ≤ 1.5 as before, we also modify the stability of E2F1 and AP-1. That is, we relax the values of φE2F1 and φAP-1 in (7) and (8) from from their fixed nominal values 0.1 and 0.01 respectively, to allow searching within the ranges,
0.08 ≤ φE2F1 ≤ 0.12
0.008 ≤ φAP-1 ≤ 0.012.
We propose an inverse problems methodology to address a class of problems arising in the study and design of gene systems. The methodology is based on the (inverse) map from the space of bifurcation diagrams to the space of biochemical parameters. Within this general inverse bifurcation framework, we show that many questions of biological interest may be formulated as optimization problems involving minimal distances to bifurcation manifolds. Adjoint sensitivity analysis is carried out to allow for efficient solution using gradient-based optimization methods.
Several general questions remain for the methodology: what biological problems can be formulated in terms of geometric properties of the bifurcation diagrams? Given desired geometric properties for bifurcation diagrams, how best can the problem be formulated mathematically to allow for tractable solution? Finally, we remark that the development of appropriate numerical methods, novel computational strategies and regularization techniques are necessary to allow for the upscaling of the current computational tool to high dimensional biological applications.
We gratefully acknowledge Stefan Müller for providing the repressilator model, Lukas Endler for helpful discussions on the cell cycle model, as well as Andreas Hofinger and Philipp Küigler for comments on the paper. This work is supported by the WWTF, project number MA05.
- Kuznetsov YA: Elements of Applied Bifurcation Theory. 2004, New York, USA: Springer-VerlagView ArticleGoogle Scholar
- Thomas R, D'Ari R: Biological feedback. 1990, Florida, USA: CRC-PressGoogle Scholar
- Hasty J, Dolnik M, Rottschafer V, Collins JJ: Synthetic gene network for entraining and amplifying cellular oscillations. Phys Rev Lett. 2002, 88 (14): 148101-PubMedView ArticleGoogle Scholar
- Tyson JJ, Chen KC, Novak B: Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol. 2003, 15 (2): 221-231.PubMedView ArticleGoogle Scholar
- Angeli D, Ferrell JEJ, Sontag ED: Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems. Proc Natl Acad Sci USA. 2004, 101 (7): 1822-1827.PubMedPubMed CentralView ArticleGoogle Scholar
- Widder S, Schicho J, Schuster P: Dynamic patterns of gene regulation I: simple two gene systems. J Theor Biol.Google Scholar
- Müller S, Hofbauer J, Endler L, Flamm C, Widder S, Schuster P: A generalized model of the repressilator. J Math Biol.Google Scholar
- Tyson JJ, Csikasz-Nagy A, Novak B: The dynamics of cell cycle regulation. Bioessays. 2002, 24 (12): 1095-1109.PubMedView ArticleGoogle Scholar
- Battogtokh D, Tyson JJ: Bifurcation analysis of a model of the budding yeast cell cycle. Chaos. 2004, 14 (3): 653-661.PubMedView ArticleGoogle Scholar
- Chickarmane V, Paladugu SR, Bergmann F, Sauro HM: Bifurcation discovery tool. Bioinformatics. 2005, 21 (18): 3688-3690.PubMedView ArticleGoogle Scholar
- Engl HW, Hanke M, Neubauer A: Regularization of inverse problems, Volume 375 of Mathematics and its Applications. 1996, Dordrecht:Kluwer Academic Publishers GroupView ArticleGoogle Scholar
- Engl HW, Kügler P: Nonlinear inverse problems: theoretical aspects and some industrial applications. Multidisciplinary Methods for Analysis, Optimization and Control of Complex Systems. Edited by: V Capasso JP. 2005, 3-48. New York: Springer VerlagView ArticleGoogle Scholar
- Iwasaki K, Kamimura Y: An inverse bifurcation problem and an integral equation of the Abel type. Inverse Problems. 1997, 13: 1015-1031.View ArticleGoogle Scholar
- Iwasaki K, Kamimura Y: Inverse bifurcation problem, singular Wiener-Hopf equations, and mathematical models in ecology. J Math Biol. 2001, 43: 101-143.PubMedView ArticleGoogle Scholar
- Dobson I: Computing a closest bifurcation instability in multidimensional parameter space. J Nonlinear Sci. 1993, 3 (3): 307-327.View ArticleGoogle Scholar
- Mönnigmann M, Marquardt W: Normal vectors on manifolds of critical points for parametric robustness of equilibrium solutions of ODE systems. J Nonlinear Sci. 2002, 12 (2): 85-112.View ArticleGoogle Scholar
- Dobson I: Distance to bifurcation in multidimensionalparameter space: margin sensitivity and closest bifurcations. Bifurcation control, theory and application. Edited by: Chen G, Hill DJ, Yu X. 2003, 49-66. New York: Springer VerlagGoogle Scholar
- Morohashi M, Winn AE, Borisuk MT, Bolouri H, Doyle J, Kitano H: Robustness as a measure of plausibility in models of biochemical networks. J Theor Biol. 2002, 216: 19-30.PubMedView ArticleGoogle Scholar
- Hendenson ME: Multiple parameter continuation: computing implicitly defined k-anifolds. J Bifur Chaos. 2002, 12: 451-476.View ArticleGoogle Scholar
- Dhooge A, Govaerts W, Kuznetsov YA: MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Software. 2003, 29 (2): 141-164.View ArticleGoogle Scholar
- MathWorks Ltd. http://www.mathworks.com/
- Shapiro BE, Hucka M, Finney A, Doyle J: MathSBML: a package for manipulating SBML-based biological models. Bioinformatics. 2004, 20: 2829-2831.PubMedPubMed CentralView ArticleGoogle Scholar
- Mathematica Symbolic Toolbox for MATLAB – Version 2.0. http://library.wolfram.com/infocenter/MathSource/5344/
- Elowitz MB, Leibler S: A synthetic oscillatory network of transcriptional regulators. Nature. 2000, 403 (6767): 335-338.PubMedView ArticleGoogle Scholar
- Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339-342.PubMedView ArticleGoogle Scholar
- Atkinson MR, Savageau MA, Myers JT, Ninfa AJ: Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli. Cell. 2003, 113 (5): 597-607.PubMedView ArticleGoogle Scholar
- Gill PE, Murray W, Wright MH: Practical Optimization. 1981, London: Academic PressGoogle Scholar
- Swat M, Kel A, Herzel H: Bifurcation analysis of the regulatory modules of the mammalian G1/S transition. Bioinformatics. 2004, 20 (10): 1506-1511.PubMedView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.