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 [8] 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 [18]. 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 [19]. Instead, formulations based on continuation along rays in parameter space are preferred. Below, we formulate distance to bifurcation manifolds (first introduced in [15]) 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.
Consider the splitting of m-dimensional parameter space P ⊂ ℝm into input and system parameters, p = (pi, ps) ∈ Pi × ps . For an ODE system, let Σ denote a bifurcation manifold of interest, consisting of sets in parameter space P for which structural stability breaks down [1]. For a given system parameter ps, we further define Σ(ps) ≡ Σ ∩ {ps} as the intersection of Σ with the ps-plane. In Figure 2, the geometric relationship between Σ and Σ(ps) is illustrated. Let the forward operator F :P → P be a mapping in parameter space, taking a given point to its orthogonal projection on Σ(ps), assumed to be well-defined. That is,
Figure 1  Forward map in the input plane, Pi.
Figure 2  Illustration of normal vector N and components Ni, Ns. 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 : Pi → ℝ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,
<F'(p)δp, δ> = <δp, F'*(p)δ>,     (3)
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 - pi||, the adjoint solution for the linearization J'(p)(·) is given in terms of the vector normal to the manifold Σ. Denoting Ni and Ns as the components of the normal vector in Pi and Ps respectively, it can be shown that (up to sign),
To fix the sign of ψ, the component Ni is chosen to be pi - 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, Ns. 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 fx 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,
fxv = ωcrit v
.
The expressions for normal vectors are given as:
where superscript H denotes conjugate transpose and . These expressions above prescribe the components of the adjoint solution, thus enabling efficient gradient calculation via (4). Now, we briefly mention methods for computing the projection F(p). Methods of iterative and direct type for finding the (locally) closest bifurcation point have been derived [17]. In the current work we use the former approach, based on using the component of the normal vector in the input plane, Ni. Provided certain conditions on the principal curvatures of Σ(ps) are met, geometric convergence is assured. The algorithm is discussed in Section 3. Figure 3 illustrates the method in a simple example, producing a sequence of iterates (i) converging to the point F(p)i that is closest to a (non-convex) neighboring region with respect to pi.
Figure 3  Demonstration of iterative procedure for computing F(p). 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.