3 Algorithm and software implementation
Here we give an outline of the algorithm for general inverse bifurcation problems. The main ingredients are applications of the projection operator F(p), as well as the adjoint, F'*(p). The computation of the former is denoted by the routine APPLYF (see Figure 4). Each time APPLYF is called, corrector steps (using, for instance, Newton's method) have to be carried out on the previously computed xinit to find the initial solution for the current value of ps. Once the corrected solution xcorr is computed, the iterative procedure LOCMINDIST (see Figure 5) is called to compute the nearest point on the bifurcation manifold of interest. This procedure is based on a series of one-parameter continuations and gradient calculations according to (6). The inputs include the following: v the initial search direction in parameter space; εproj the relative tolerance on the iterative proceure. In general, several search directions (denoted by V = {v1,⋯, vmax}) have to be used to approximate the globally closest point. However, for the examples shown in the paper the initial search space V is only 1 dimensional. Once F(p) and the corresponding solution are obtained, F'*(p) is then computed. The derivative information, together with constraints (plow, pupp, c(p)) and Hessian approximation (HessA) are then used to compute an SQP step and update the Hessian. Figure 6 describes the inverse bifurcation algorithm.
Figure 4  Algorithm APPLYF.
Figure 5  Algorithm LOCMlNDlST.
Figure 6  Algorithm INVERSE BIFURCATION. Our Inverse Bifurcation Toolbox is an implementation of the above described algorithm. It combines the capability of several packages: the MATCONT package [20] for performing the underlying one-parameter continuations, the MATLAB Optimization Toolbox [21] for performing gradient-based constrained optimization, the MathSBML package [22] for reading in biological models in the SBML format and finally the Mathematica Symbolic Toolbox for MATLAB [23] 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 fx, fp, fxp and fxx. 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).