Efficient classification of complete parameter regions based on semidefinite programming
- Lars Kuepfer^{1},
- Uwe Sauer^{1}Email author and
- Pablo A Parrilo^{2}
DOI: 10.1186/1471-2105-8-12
© Kuepfer et al; licensee BioMed Central Ltd. 2007
Received: 22 September 2006
Accepted: 15 January 2007
Published: 15 January 2007
Abstract
Background
Current approaches to parameter estimation are often inappropriate or inconvenient for the modelling of complex biological systems. For systems described by nonlinear equations, the conventional approach is to first numerically integrate the model, and then, in a second a posteriori step, check for consistency with experimental constraints. Hence, only single parameter sets can be considered at a time. Consequently, it is impossible to conclude that the "best" solution was identified or that no good solution exists, because parameter spaces typically cannot be explored in a reasonable amount of time.
Results
We introduce a novel approach based on semidefinite programming to directly identify consistent steady state concentrations for systems consisting of mass action kinetics, i.e., polynomial equations and inequality constraints. The duality properties of semidefinite programming allow to rigorously certify infeasibility for whole regions of parameter space, thus enabling the simultaneous multi-dimensional analysis of entire parameter sets.
Conclusion
Our algorithm reduces the computational effort of parameter estimation by several orders of magnitude, as illustrated through conceptual sample problems. Of particular relevance for systems biology, the approach can discriminate between structurally different candidate models by proving inconsistency with the available data.
Background
Systems biology is a framework integrating diverse disciplines such as molecular biology and genetics with mathematical modelling and simulation, where computational models assume an increasingly important role in recent years. Simulations are an essential element for quantitative understanding because the highly nonlinear behavior of complex biological systems can sometimes be counterintuitive [1, 2]. System identification, which comprises both parameter estimation and structural network analysis in biological systems, can be extremely difficult since the governing principles and topological interactions are often not known [3]. The flood of data from novel high-throughput and genome-wide analyses further adds to the problem, and their diversity poses additional challenges for consistency testing and integration. Exact values of kinetic parameters (e.g., association or dissociation constants for protein-protein interaction) are very difficult to determine experimentally, because they are a function of the heterogeneous spatial and developmental cellular conditions.
Most approaches to mathematical modelling of biological systems either avoid parametrization by focusing on static steady-state models of reaction stoichiometry or topology [4, 5], approximate behavior and regulation through heuristic-based approaches [6, 7], or reduce the solution space by applying enzyme kinetics such as in traditional Michaelis-Menten or linlog kinetics, respectively [8, 9]. While this enables the detailed modelling of time-dependent processes [10], the underlying quasi steady-state assumption a priori neglects the dynamics of enzyme complexes. Thus, in many cases the use of mass-action kinetics in the form of elementary reactions becomes necessary, for instance when some of the enzyme kinetics do not follow Michaelis-Menten [11], a steady-state of the system cannot be assumed [12], transport functions have to be described [13] or the network structure is uncertain [14, 15].
The systems of equations that arise in time-dependent models based on chemo-physical kinetics of any kind are usually hard to parameterize, since their components are tightly coupled and there is limited information about the time evolution of the concentrations of all the species. Several approaches for efficient system identification have been developed, ranging from reduction of system dimensionality [16] to decoupling of the differential equations [17]. The actual parametrization algorithm, here referred to as the wrapping algorithm, is of particular importance [18] and was for example highlighted in a comparison between gradient-based methods and genetic and evolutionary algorithms [19]. All of these approaches however consider single, isolated points in the parameter space; this can be very time-consuming due to the iterative nature of the parametrization process [19] and furthermore cannot guarantee that the algorithm finds possible solutions.
We present here a conceptually novel approach based on techniques from numerical convex optimization, in particular semidefinite programming (SDP) [20], that can efficiently partition the parameter space into feasible and infeasible regions. Basically, our approach allows to analyze sets of mass action kinetics, i.e., ordinary differential equations (ODEs) at steady state, under consideration of additional algebraic equality and inequality constraints (e.g., mass balances or formation rates). It is thus possible to simultaneously consider all the available information during the parameter estimation itself, rather than checking consistency a posteriori as in most current methods.
Previous applications used SDP either for rational experimental design based on covariance analysis [21] or for the derivation of barrier certificates by construction of surrogate models [22]. In contrast, the algorithm presented here needs no problem reformulation but is based directly on the original parameter estimation problem. By exploiting convexity in the search for feasible steady state concentrations, our methods enable a multi-dimensional perturbation analysis for sets of parameters. Our algorithm therefore provides rigorous proofs for the classification of the parameter space into feasible and infeasible regions, thus reducing the number of points needed for consideration by several orders of magnitude.
Methods
Steady state analysis
The natural starting point in the analysis of dynamical systems in biology is the determination of a steady state equilibrium of the model, since it represents the reference point for any kind of perturbation introduced to the system. It is therefore of utmost importance that the species' concentrations at an equilibrium point, which is henceforth synonymously referred to as steady state, are realistic and validated against as much data as possible.
In the present study we consider models in the form of mass action kinetics, hence all reaction rates are direct functions of the concentrations y and the kinetic constants k. At an equilibrium point, the corresponding system of n ODEs reduces to a set of polynomial equations
f_{ i }(k, y) = 0, i = 1, ..., n. (1)
Besides the equations above, the state concentrations also have to satisfy various kinds of constraints based on experimental data, e.g., formation rates or overall mass balances. Since these will in practice inevitably include certain errors, inequality constraints also arise naturally. Consider for instance the general mass balance equations, here described by a matrix M, which when multiplied with the array of state variables y satisfy
M·y = b, M ∈ ℝ^{m × n}, (2)
where b ∈ ℝ^{ m }is the total amount of each species. Hence, Eqn. (2) represents overall mass balances while Eqn. (1) comprises the mass action kinetics, i.e., the formation and consumption rates of each single species. Since components can exist both unbound or bound in a complex, respectively, there are in general fewer mass balance equations than differential equations (n ≥ m). An experimental error ε can be taken into account by relaxing the exact mass balance equations to inequalities where
(1 - ε)·b ≤ M·y ≤ (1 + ε)·b. (3)
Thus, the parameter estimation problem can be rewritten as a nonlinear optimization
that depends on the steady state concentrations y and the set of kinetic parameters k. Here, f_{0} (k, y) is an arbitrary objective function, e.g., the overall model error, the equality constraints f_{ i }(k, y) represent the steady state condition and g (k, y) is a set of inequality constraints, e.g., the mass balances.
Semidefinite programming
In the following, we will introduce a reformulation of the generalized parameter estimation problem (4) based on semidefinite programming (SDP). SDP is a specific kind of convex optimization problem [20, 23, 24], with very appealing numerical properties. An SDP problem corresponds to the optimization of a linear function subject to a matrix inequality. The only prerequisite for the applicability of the SDP is a quadratic (or polynomial) representation of the original set of equalities and inequalities (4), which allows a reformulation/relaxation in terms of symmetric matrices. For clarity of exposition, we focus on the case where the f_{ i }are quadratic and the g_{ i }are linear, although our methods extend to the fully polynomial case; see [25] for details.
For the SDP relaxation we define new variables x, X in terms of the original state variables y by:
Based on these definitions, the original problem (4), including both steady state equations and data consistency inequalities, can be rewritten as:
Here, the symmetric matrix Q_{ obj }defines an objective function (e.g., the identity matrix), that selects a specific solution out of the possibly many equilibria. The symmetric matrices Q_{ i }correspond to the set of ODEs at steady state (Eqn. (1)) and the matrix L to linear inequality constraints derived, for instance, from approximate mass-balance equations (Eqn. (3)). Note, that in general the matrices Q_{ i }are a function of the set of (generally unknown) rate parameters k_{ j }, while L, which represents the mass balances, is not. The basic convex relaxation described in the Appendix then yields the following SDP relaxation of (6):
where Tr is the trace operator, which adds up the diagonal elements, and the inequality in the last line indicates that the matrix X must be positive semidefinite, i.e., all its eigenvalues should be greater than or equal to zero. The vector e_{1} is an all-zero vector, except for its first entry, which equals one to enforce X_{11} = 1. The set of feasible solutions, i.e., the set of matrices X that satisfy the constraints, is always a convex set. Recall that the matrix X as defined in Eqn. (5) is by construction a rank one matrix. This rank condition, however, is not guaranteed for the X obtained from the optimization, and thus this property has to be checked independently after a solution to the SDP is computed. This is necessary because the relaxed problem formulation (7) is less strict then the original form (6), hence the set of feasible solutions becomes larger. Note, that introduction of additional nonlinear constraints helps to reduce the set of " false positive" solutions, an aspect discussed in more detail in an additional section below. Finally, in the particular case of Q_{ obj }= 0, the problem reduces to whether or not the inequality can be satisfied for some matrix X. In this case, the SDP is referred to as a feasibility problem, where we are interested in proving the mere existence of solutions rather than finding any particular one.
Dual SDP problems
The convexity of SDP has made it possible to develop sophisticated and reliable analytical and numerical methods to solve them [20]. A very important feature of SDP problems, from both the theoretical and applied viewpoints, is the associated duality theory (see also Appendix). For every SDP of the form (7) (usually called the primal problem), there is another associated SDP, called the dual problem, which can be derived via Lagrangian duality:
where the constraint represents a matrix inequality with r ≥ 0, S ≥ 0, S_{ ii }= 0. The dual variables γ, λ, r, S are Lagrange multipliers associated to the different constraints in the primal problem. In the case of feasibility problems (i.e., Q_{ obj }= 0), the dual problem can be used to certify the nonexistence of solutions of the primal. This property will be crucial in our developments.
Results
Parameter estimation for nonlinear mass action kinetics at steady state
Identifying a satisfactory equilibrium point of a set of ODEs requires the solution of a system of algebraic equations (those that define the steady-state conditions) subject to additional equality or inequality constraints (consistency with experimental data). This can be time-demanding, because the system must often be simulated from a given set of initial conditions until it settles to an equilibrium [11, 26]. This is particularly troublesome when the computed solution violates the experimental constraints and must hence be discarded. Traditional heuristic tools for this optimization problem require, for instance, an iterative procedure where the candidate sets of parameters are generated by some kind of wrapping parametrization algorithm (e.g., gradient-based or evolutionary), and a subsequent consistency check by integrating the system equations for these parameter values [19]. This trial and error method can be very time consuming because it only allows to check consistency in a subsequent step. Hence, an algorithm that searches for steady state solutions but is guaranteed to be consistent with all experimental data would be extremely valuable.
SDP and nonlinear systems of equations
As opposed to these "indirect" techniques, our method is a conceptually novel direct approach based on a convex relaxation of the generalized parameter estimation problem at steady state (4). Our techniques apply whenever the model presentation is in polynomial form, since in this case the resulting system can be relaxed into a semidefinite programming problem (6), which in turn can be solved using efficient interior-point methods [27–30].
We illustrate the application of our approach with the following nonlinear system given by
Here, two components A and B form a complex, A • B, that transforms into a modified form, A • B*, and finally decays into its building blocks. This system yields a set of four ODEs that at steady state reduces to a polynomial system:
If it is furthermore known from experimental data that in equimolar concentrations of A and B one third of them is unbound while the rest is associated in one of both complexes, we can write mass balances as
which have to be kept within a certain experimental accuracy ε = 20%.
The example above is thus a system of nonlinear polynomial relations based on mass action kinetics that comprises both equalities (Eqns. (10)) and inequalities (Eqns. (11)). As explained earlier, we define the state variables
y^{ T }= ([A], [B], [A • B], [A • B*]), (12)
which are then used to produce a new matrix variable X according to Eqn. (5). To write the SDP problem, the equality and inequality constraints, Eqns. (10) and Eqns. (11), respectively, have to be rewritten into the form of Eqn. (6). For example, the steady state equation for species B becomes
while the rows of L corresponding to the mass balance inequalities for component A are
The feasible set of the relaxed SDP problem (7) is always convex and includes all the equilibria of the original problem. Thus, SDP optimization enables the possibility of directly solving the underlying nonlinear system of algebraic equations, or of proving its inconsistency. Furthermore, this property makes possible a direct consistency check based on feasibility or infeasibility of the SDP optimization. Typically, a parameterization approach would require identification of a stable steady state, either by root-finding or numerical integration, at which point the simulation can be verified or falsified for a given set of parameters k. In contrast, the SDP-based approach bypasses this time consuming step: if the problem (7) is infeasible the parameters k are inconsistent with the experimental data, and if a feasible solution of rank one can be found, the given parameters are valid.
For our running example, consider a set of parameters, k_{1} = 3, k_{2} = 1 and k_{3} = 1, respectively, for which the SDP is feasible. Let the solution be
which is rank one. Using the largest eigenvalue of X to scale the corresponding eigenvector yields
A simple check shows that ${y}_{sol}^{T}$ directly fulfills the steady state condition and all additional constraints.
Additional nonlinear constraints
The decision variable X in Eqn. (5) is by construction a rank one matrix. This property of X, however, is not necessarily guaranteed when solving the SDP relaxation (7) (including this constraint would destroy convexity). Hence, it has to be verified independently once the optimization is completed. From our numerical experience, a very limited number of solutions fails the rank one condition such that the corresponding results have to be discarded. A simple but effective combination of constraints, however, allows to reduce the number of these "false positive" solutions: since most of the inequality constraints are linear (e.g., mass balance equations), they can be used to tighten the description of the feasible set by using redundant constraints. By exhaustive multiplication of m pairs of linear inequality constraints, i.e. each upper bound is multiplied with every possible lower bound, $\left(\begin{array}{c}m\\ 2\end{array}\right)$ new quadratic inequalities are generated. This multiplication of constraints corresponds to the term L·X·L^{ T }appearing in (7). Note, that even the product of a lower and an upper bound of the same constraint can be helpful, because the thus generated matrix has nonzero diagonal (i.e., quadratic) elements. These additional nonlinear algebraic constraints improve the performance of the algorithm significantly, because the set of feasible solutions that do not meet the rank one condition becomes much smaller.
Analyzing whole regions of parameter space
The approach presented in the previous paragraphs can directly handle nonlinear algebraic equations in polynomial form. However, it considers only single sets of parameters, that have to be provided by an external parameter estimation algorithm. Since in general the parameter space can be quite large, it would be very helpful to be able to discard a priori large regions of the space, where we know for sure that no consistent rates can be found. As we will see, we can achieve this goal in a very efficient way by considering the dual optimization problem (8) presented earlier.
As in general convex optimization, infeasibility of the primal problem can be proven by the existence of a dual feasible solution, and conversely, dual infeasibility follows from a primal solution. Thus, the existence of dual variables satisfying the dual constraints in (8) directly proves primal infeasibility. This criterion can be used to yield an efficient procedure for the analysis of the parameter space. Therein, large regions can be detected where the given model is inconsistent with the experimental data. The feasible region will usually represent a small fraction in the overall parameter space. Therefore an efficient search of the parameter space will focus rather on negative (inconsistent) than positive (consistent) regions.
The use of SDP-based parametrizations allows to obtain exactly this information. Once dual feasibility has been proven for a given set of parameters, this point in the solution space is associated to specific values of the dual variables. This information can then be used to extend the feasibility check to larger regions in the parameter space. Recall that only the matrices Q_{ i }, i.e., the mass action kinetics, are actually dependent on the rate parameters k. We thus used a bisection approach in which the current parameter set is an n-dimensional box. Since this is a polytope, it follows that if the conditions are valid on the corners then they are also valid on the interior, and thus can be neglected for further analysis. In particular, we determine a factor η by which each parameter can be perturbed such that dual feasibility holds. We remark that the choice of cuboid-like regions is only a matter of convenience, the results can be easily extended in a direct fashion to any other polyhedral partition scheme (for instance, using simplices).
Pseudo-code for efficient classification of the parameter space.
let K ⊂ ℝ^{ n }be the parameter space |
---|
begin: set K* = K and K_{ inconsistent } = ∅ |
while (K* ≠ ∅) |
consider loop i |
generate parameter subset (k_{ i }± Δk) ∈ K* |
perform dual SDP optimization |
if solution (γ, λ_{ i }, P, r, S) exists |
// k_{ i }is dual feasible/primal infeasible |
while (dual feasibility holds) |
increase Δk_{ i } |
end |
K* = K*\(k_{ i }± Δk_{ i }) |
K_{ inconsistent }= K_{ inconsistent }∪ (k_{ i }± Δk_{ i }) |
end |
end |
Efficient classification of the parameter space
Rigorous proofs for model discrimination
Comparison of model alternatives for enzymatic conversion of B to D.
Model alternative I: |
---|
$\begin{array}{lllll}[B]+[{E}_{3}]\hfill & \underset{{k}_{off,3}}{\overset{{k}_{on,3}}{\rightleftarrows}}\hfill & [B\xb7{E}_{3}]\hfill & \stackrel{{k}_{cat,3}}{\to}\hfill & [D]+[{E}_{3}]\hfill \end{array}$ |
Model alternative II: |
$\begin{array}{lll}[B]+[{E}_{3}]\hfill & \underset{{k}_{off,3}}{\overset{{k}_{on,3}}{\rightleftarrows}}\hfill & [B\xb7{E}_{3}]\hfill \end{array}$ |
$\begin{array}{lllll}[B\xb7{E}_{3}]+[A]\hfill & \underset{{k}_{off,3b}}{\overset{{k}_{on,3b}}{\rightleftarrows}}\hfill & [B\xb7{E}_{3}\xb7A]\hfill & \stackrel{{k}_{cat,3b}}{\to}\hfill & [D]+[B\xb7{E}_{3}]\hfill \end{array}$ |
Discussion
We introduced a conceptually novel approach for parameter estimation and parameter space classification based on a direct search for solutions that are simultaneously consistent with the steady-state concentrations and the inequalities arising from experimental data. The method is based on semidefinite optimization, whose convexity properties allow for the direct solution (of relaxations) of nonlinear optimization problems in polynomial form. It is thus possible to establish the possible infeasibility of steady state concentrations for a given set of parameters under direct consideration of experimental constraints. Our approach avoids the time-consuming numerical identification of a stable steady state, where feasibility can only be validated in retrospective [11, 26]. Besides the direct consistency analysis for single sets of parameters, the duality properties of SDP optimization problems allow the direct proof of feasibility (or infeasibility) of whole regions in the parameter space instead of the mere consideration of isolated spots. This significantly reduces the total number of possibilities that have to be evaluated. Moreover, our approach is based on a simple convex relaxation of the original parametrization problem and can hence easily be applied without time consuming problem reformulation.
The possibility of discarding whole regions of the parameter space from further exploration is extremely valuable for model discrimination, where determination of model consistency or inconsistency with experimental data is the desired goal [14]. This approach can be very time-intensive due to the trial and error method that includes several integrations per set of parameters. Our algorithm thus provides an valuable tool for model discrimination.
Despite the immediate practical applications of our method, further work is necessary to fully exploit the total possibilities of the algorithm presented. Since the problem size increases with the number of variables of the system, the corresponding optimization problems will inevitably result in large matrices. Current versions of SDP solvers [27–30], however, slow down considerably when larger problem sets with more than about 30 state variables were analyzed. Stiffness of the underlying system of equations is another difficulty, which remains unsolved to date. Hence, problems which large differences in the kinetic parameters, e.g., if Michaelis-Menten kinetics are transformed to mass action kinetics, can probably not be solved without sophisticated preprocessing. An adequate way of scaling therefore needs to be developed. Additionally, the efficiency of the wrapping parametrization algorithm itself could also be improved, since currently there is no control level algorithm that supervises the search direction in the parameter space. Valuable information, however, is available, since the shape of the boxes indicates the maximal possible perturbation η, and this could be used to determine the direction of the next step (Fig. 2). A promising approach would thus be a method that simultaneously uses primal and dual information.
Conclusion
In conclusion, we believe the results shown here are an important first step towards the integration of SDP as a tool to solve and analyze polynomial systems in chemical and biochemical engineering. Since our SDP-based algorithm allows to increase the efficiency of the pivotal steps in parameter estimation, it has great potential for the identification of nonlinear systems that prevail in biology.
Appendix
Convex relaxation
Consider a set of quadratic equations and linear inequalities for the vector x as defined in (5).
x^{ T }·Q_{ i }·x = 0, L·x ≥ 0. (16)
Defining X = x·x^{ T }, we find that these equations can be rewritten as affine expressions in the matrix X. A semidefinite relaxation of Eqns. (16) is then obtained by replacing the nonconvex constraint X = x·x^{ T }with the weaker (but convex) alternative:
X_{11} = 1, X ≥ 0,
yielding the system
Tr(Q_{ i }·X) = 0, L·X·e_{1} ≥ 0, (17)
with X_{11} = 1, X ≽ 0.
In general, the original (Eqns. (16)) and the new (Eqns. (17)) formulations are equivalent only if the matrix X has rank one. However, since such rank constraint is nonconvex, it cannot be directly included in the SDP optimization. This relaxation causes the set of solutions to becomes larger, and as a consequence the rank condition must be checked independently after the optimization is completed. Notice also that the formulation can be strengthened by adding the redundant constraints L·X·L^{ T }≥ 0.
Weak Duality in SDP problems
A typical SDP problem in primal form is
min Tr(C·X)
s.t Tr(A_{ i }·X) = b_{ i }, i = 1, ..., m (18)
X ≽ 0.
The associated dual problem can be stated as
where b = (b_{1},...,b_{ m }), and the vector z = (z_{1},...,z_{ m }) contains the dual decision variables.
The key relationship between the primal and the dual problem is the fact that feasible solutions of one can be used to bound the values of the other problem. Indeed, let X and z be any two feasible solutions of the primal and dual problems respectively. Then we have the following inequality:
Tr(C·X) ≥ b^{ T }z. (20)
This property is known as weak duality. Thus, we can use any feasible X to compute an upper bound for the optimum of b^{ T }z, and we can also use any feasible z to compute a lower bound for the optimum of Tr(C·X).
Proof for regions of inconsistency
If a feasible dual solution can be found, Eqn. (8) will be satisfied. As noticed earlier, the matrices Q_{ i }explicitly depend on the unknown rates k in an affine way. Thus, we want to be able to disprove the consistency of not just one particular fixed value of the rate constants, but to simultaneously discard whole sets of rates at the same time. In other words, we want to find solutions (γ, λ_{ i }, r, S) of the dual form (8) that work for all rate constants k on a given set.
Notice that the dependence of the matrices Q_{ i }on the rate constants k is affine. Assume the nominal value and deviation of rate constant j is k_{j 0}± Δk_{ j }. Then, we can write
where δ_{ j }:= (k_{ j }- k_{j 0})/Δk_{ j }, and |δ_{ j }| ≤ 1. To guarantee that the dual form (8) holds for all allowable values of the rate constants, we can use the following result:
Lemma 1 Consider the linear matrix inequality given by:
If there exist η, W _{ i } , such that
for i = 1,..., n, then Eqn. (21) holds for all δ such that |δ_{ i }| ≤ η.
The lemma follows easily from the identity:
Since the expressions in (22) are affine in the unknowns η, W_{ i }, we can hence directly use this lemma to obtain guaranteed regions of inconsistency.
Further case studies
Additional example 1:
Consider a simple two-element linear reaction, where component A is converted into component B and vice versa:
Additional example 2:
Consider another simple nonlinear system, where mRNA binds to a ribosome to form a protein P, which decays at a certain rate:
Declarations
Authors’ Affiliations
References
- Kitano H: Systems biology: a brief overview. Science 2002, 295(5560):1662–1624. 10.1126/science.1069492View ArticlePubMedGoogle Scholar
- Bailey JE: Mathematical modeling and analysis in biochemical engineering: past accomplishments and future opportunities. Biotechnol Prog 1998, 14: 8–20. 10.1021/bp9701269View ArticlePubMedGoogle Scholar
- Stelling J, Sauer U, Szallasi Z, Doyle FJ 3rd, Doyle J: Robustness of cellular functions. Cell 2004, 118(6):675–685. 10.1016/j.cell.2004.09.008View ArticlePubMedGoogle Scholar
- Edwards JS, Covert M, Palsson B: Metabolic modelling of microbes: the flux-balance approach. Environ Microbiol 2002, 4(3):133–140. 10.1046/j.1462-2920.2002.00282.xView ArticlePubMedGoogle Scholar
- Kuepfer L, Sauer U, Blank LM: Metabolic functions of duplicate genes in Saccharomyces cerevisiae . Genome Res 2005, 15: 1421–1430. 10.1101/gr.3992505PubMed CentralView ArticlePubMedGoogle Scholar
- Varner JD: Large-scale prediction of phenotype: concept. Biotechnol Bioeng 2000, 69(6):664–678. 10.1002/1097-0290(20000920)69:6<664::AID-BIT11>3.0.CO;2-HView ArticlePubMedGoogle Scholar
- Lee B, Yen J, Yang L, Liao JC: Incorporating qualitative knowledge in enzyme kinetic models using fuzzy logic. Biotechnol Bioeng 1999, 62(6):722–729. 10.1002/(SICI)1097-0290(19990320)62:6<722::AID-BIT11>3.0.CO;2-UView ArticlePubMedGoogle Scholar
- Bailey J, Ollis D: Biochemical Engineering Fundamentals. McGraw-Hill chemical engineering series. 1986.Google Scholar
- Heijnen JJ: Approximate kinetic formats used in metabolic network modeling. Biotechnol Bioeng 2005, 91(5):534–545. 10.1002/bit.20558View ArticlePubMedGoogle Scholar
- Schoeberl B, Eichler-Jonsson C, Gilles ED, Muller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol 2002, 20(4):370–375. 10.1038/nbt0402-370View ArticlePubMedGoogle Scholar
- Mishra J, Bhalla US: Simulations of inositol phosphate metabolism and its interaction with InsP(3)-mediated calcium release. Biophys J 2002, 83(3):1298–1316.PubMed CentralView ArticlePubMedGoogle Scholar
- Stelling J, Gilles ED, Doyle F Jr 3rd: Robustness properties of circadian clock architectures. Proc Natl Acad Sci U S A 2004, 101(36):13210–13215. 10.1073/pnas.0401463101PubMed CentralView ArticlePubMedGoogle Scholar
- Smith AE, Slepchenko BM, Schaff JC, Loew LM, Macara IG: Systems Analysis of Ran transport. Science 2002, 295(5554):488–491. 10.1126/science.1064732View ArticlePubMedGoogle Scholar
- Haunschild MD, Freisleben B, Takors R, Wiechert W: Investigating the dynamic behavior of biochemical networks using model families. Bioinformatics 2005, 21(8):1617–25. 10.1093/bioinformatics/bti225View ArticlePubMedGoogle Scholar
- Barkai N, Leibler S: Robustness in simple biochemical networks. Nature 1997, 387(6636):913–917. 10.1038/43199View ArticlePubMedGoogle Scholar
- Bentele M, Lavrik I, Ulrich M, Stosser S, Heermann DW, Kalthoff H, Krammer PH, Eils R: Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis. J Cell Biol 2004, 166(6):839–851. 10.1083/jcb.200404158PubMed CentralView ArticlePubMedGoogle Scholar
- Voit EO, Almeida J: Decoupling dynamical systems for pathway identification from metabolic profiles. Bioinformatics 2004, 20(11):1670–1681. 10.1093/bioinformatics/bth140View ArticlePubMedGoogle Scholar
- Polisetty PK, Voit EO, Gatzke EP: Identification of metabolic system parameters using global optimization methods. Theor Biol Med Model 2006., 3(4):Google Scholar
- Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Res 2003, 13(11):2467–2474. 10.1101/gr.1262503PubMed CentralView ArticlePubMedGoogle Scholar
- Vandenberghe L, Boyd S: Semidefinite programming. SIAM Rev 1996, 39: 49–95. 10.1137/1038003View ArticleGoogle Scholar
- Flaherty P, Jordan MI, Arkin AP: Robust Design of Biological Experiments. Proceedings of the Neural Information Processing Symposium 2005 2005.Google Scholar
- Tau JF, Fazel M, Liu X, Otitoju T, Papachristodoulou A, Prajna S, Doyle J: Application of Robust Model Validation Using SOSTOOLS to the Study of G-Protein Signaling in Yeast. Proceedings of FOSBE (Foundations of Systems Biology and Engineering) 2005.Google Scholar
- Todd MJ: Semidefinite Optimization. Acta Numerica 2001, 10: 515–560.View ArticleGoogle Scholar
- Boyd S, Vandenberghe L: Convex Optimization. Cambridge University Press; 2004.View ArticleGoogle Scholar
- Parrilo PA: Semidefinite programming relaxations for semialgebraic problems. Math Prog 2003, 96(2, Ser B):293–320. 10.1007/s10107-003-0387-5View ArticleGoogle Scholar
- Kremling A, Fischer S, Gadkar K, Doyle FJ, Sauter T, Bullinger E, Allgower E, Gilles ED: A benchmark for methods in reverse engineering and model discrimination: problem formulation and solutions. Genome Res 2004, 14(9):1773–1785. 10.1101/gr.1226004PubMed CentralView ArticlePubMedGoogle Scholar
- Lofberg J: YALMIP: A Toolbox for Modeling and Optimization in MATLAB. Proceedings of the CACSD Conference 2004.Google Scholar
- Sturm JF: Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software (Special issue on Interior Point Methods) 1999, 11/12: 625–653.View ArticleGoogle Scholar
- YALMIP[http://control.ee.ethz.ch/~joloef/yalmip.php]
- SeDuMi[http://sedumi.mcmaster.ca]
Copyright
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.