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 multidimensional 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 highthroughput and genomewide 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 proteinprotein 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 steadystate models of reaction stoichiometry or topology [4,5], approximate behavior and regulation through heuristicbased approaches [6,7], or reduce the solution space by applying enzyme kinetics such as in traditional MichaelisMenten or linlog kinetics, respectively [8,9]. While this enables the detailed modelling of timedependent processes [10], the underlying quasi steadystate assumption a priori neglects the dynamics of enzyme complexes. Thus, in many cases the use of massaction kinetics in the form of elementary reactions becomes necessary, for instance when some of the enzyme kinetics do not follow MichaelisMenten [11], a steadystate 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 timedependent models based on chemophysical 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 gradientbased methods and genetic and evolutionary algorithms [19]. All of these approaches however consider single, isolated points in the parameter space; this can be very timeconsuming 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 multidimensional 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 massbalance 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 allzero 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 steadystate conditions) subject to additional equality or inequality constraints (consistency with experimental data). This can be timedemanding, 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., gradientbased 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 interiorpoint methods [2730].
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 rootfinding or numerical integration, at which point the simulation can be verified or falsified for a given set of parameters k. In contrast, the SDPbased 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
= (1,1/3,1/3,1/3,1/3) and = (1/3,1/3,1/3,1/3).
A simple check shows that 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, 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 SDPbased 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 ndimensional 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 cuboidlike 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).
A formal proof of the guaranteed inconsistency of whole regions is provided in the Appendix. The pseudocode of the complete algorithm for an efficient classification of the parameter space in feasible and infeasible regions is shown in Table 1. The procedure extends the infeasibility information from an isolated set of parameters to larger regions of parameter space by using SDP optimization. All parameter sets which lie in these extended regions need not to be considered any longer for feasibility purposes. Even more, in contrast to common gridlike discretization approaches, the algorithm allows an efficient and continuous exploration of the parameter space without undergoing the danger of missing a feasible solution by accident.
Table 1. Pseudocode for efficient classification of the parameter space.
Efficient classification of the parameter space
The results in Fig. 1 correspond to the example discussed earlier (further examples are shown in the Appendix). In the figure, the contour plot representing the feasible and infeasible regions in the 3dimensional parameter space was obtained by exhaustive evaluation of sets of parameters. Particular solutions obtained with the SDPbased algorithm that prove infeasibility of whole regions are indicated in the graph. The respective regions of feasibility and infeasibility are shown, as well as a few cuboids that illustrate how we can discard the existence of solutions on whole regions of parameter space. Each box is obtained from the solution of a single small SDP optimization problem. It is clear that the SDPbased search in the infeasible region of the parameter space provides a powerful analysis tool, as it allows the efficient exploration of whole areas instead of the time consuming consistency check at single points. The figure also shows how the size of the boxes varies depending on the relative position (Fig. 1). In other words, the allowed perturbation η increases with the distance from the feasible region. This is highlighted in Fig. 2, where the distribution of η for one slice of the feasible region (in the k_{1 } k_{3 }plane) is shown. Finally, since the systems of equations and inequalities are all linear in k, the feasible region is invariant under nonnegative scalings (i.e., it is a cone), as can be seen in the threedimensional plot in Fig. 1.
Figure 1. Graph of the feasible parameter space. Contour plot of the threedimensional parameter space of k_{1}, k_{2 }and k_{3 }from the example given in the text. The cone with the black edges marks the feasible region, infeasible areas are illustrated as cuboids with gray edges. Gray circles represent the original set of parameters where the bisection search started.
Figure 2. Maximal possible parameter perturbation. The size η of the maximal possible perturbation in the example given in the text is shown as a function of k_{1 }and k_{3 }(here, k_{2 }= 1). The size of the perturbation increases with the distance from the feasible region and can be used to further refine the search directions in parameter space.
Rigorous proofs for model discrimination
The basic algorithm for exploration of the parameter space can also be used to discriminate between different model alternatives, for instance when the system structure is uncertain [14]. This is possible since SDP allows for the exploration of the complete solution space. If this is the case, the failure of parametrization (i.e., the lack of suitable values for the model parameters) directly proves model inconsistency. Moreover, we avoid the danger of missing possible solutions by accident since we consider whole regions of parameter space instead of isolated points. As an example, consider a branching point in a metabolic pathway (Fig. 3). It is known that all conversion steps are catalyzed by enzymes (E_{i}). However, while the reaction scheme for the pathway from A to C via B is well established, it may not be known whether the conversion from B to D is catalyzed by enzyme E_{3 }or whether component A exerts some cooperative effect (Fig. 3). To discriminate between these two scenarios, two model alternatives that describe possible reaction schemes can be proposed (Table 2).
Figure 3. Uncertain model structure. It is unknown whether A exerts a positive influence on the conversion of B to D. By numerical integration of model alternative II, concentrations were obtained and used in lieu of experimental data. Overall component concentrations, a step response upon an increase of v_{1 }and corresponding steady states (circle) are indicated in the inlays next to each component.
Table 2. Comparison of model alternatives for enzymatic conversion of B to D.
The consistency of both models is analyzed through defined variations in the incoming carbon flux v_{i}: An artificial data set is generated by simulation of model II, that in turn is used for parameter analysis (Figs. 3 and 4). Hence, two steady states can be simultaneously considered for SDP optimization. As expected, parameter estimation of model II identifies the feasible region in the solution space without any problems. In contrast, no feasible set of parameters is found for model I when an experimental error of ε = 0.25 is assumed (Fig. 4). This is a direct proof that model I is incorrect and needs to be modified. Only when an intolerably high deviation of ε = 0.75 is allowed, the parametrization algorithm finds solutions that are consistent with the experimental data (Fig. 4). The basic ideas underlying this small example can be directly applied to real biological problem sets where the structure is unknown. Since the search of the parameter space can be completed in a few steps, SDP provides a powerful tool for model discrimination and can be used for model invalidation and subsequent redesign of experimental setups.
Figure 4. Model discrimination based on two distinct steady states. The figure shows variations in v_{1 }simulated with model I (gray dashed line) and model II (black solid line), respectively. A large standard deviation of the experimental data (ε = 0.75, dashed error bars) has to be allowed for model I, otherwise mass balances are violated (ε = 0.25, solid error bars). For the simulations, the constants for complex formation were set to the following values: association k_{on,i }= 1; dissociation k_{off }= 0.1; catalytic step k_{cat,i }= 1. During the parametrization step, only the association constants of the complexes were estimated in a range from 0.3 to 3.
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 steadystate 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 timeconsuming 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 timeintensive 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 [2730], 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 MichaelisMenten 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 SDPbased 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.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
LK and PAP designed the study. LK conceived the modeling background. PAP derived and supervised the SDP part of the work. US participated in the design and coordinated the project. All authors read and approved the final manuscript.
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_{j0 }± Δk_{j}. Then, we can write
where δ_{j }:= (k_{j } k_{j0})/Δ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 twoelement linear reaction, where component A is converted into component B and vice versa:
We assume unimolar concentrations of species A and B that are kept within an accuracy of 25%. The twodimensional contour plot including feasible and infeasible regions is shown in Fig. 5.
Figure 5. Contour plot for additional example 1. The feasible and infeasible regions are shown in white and black, respectively. Starting from an initial set of parameters (black cross), whole areas can be proven to be inconsistent (gray).
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:
For simplicity, the overall concentrations of mRNA, ribosome and P are all equal to 1 and have to be kept within an accuracy of ε = 25%. The overall parameter space is 4dimensional. Fig. 6 shows the contour plot in the k_{3 } k_{4 }plane.
Figure 6. Contour plot for additional example 2. The feasible and infeasible regions are shown in white and black, respectively. Starting from an initial set of parameters (black cross), whole areas can be proven to be inconsistent (gray) (parameters k_{1 }and k_{2 }are fixed, k_{3 }and k_{4 }are varied in a range between 0 and 5).
References

Kitano H: Systems biology: a brief overview.
Science 2002, 295(5560):16621624. PubMed Abstract  Publisher Full Text

Bailey JE: Mathematical modeling and analysis in biochemical engineering: past accomplishments and future opportunities.
Biotechnol Prog 1998, 14:820. PubMed Abstract  Publisher Full Text

Stelling J, Sauer U, Szallasi Z, Doyle FJ 3rd, Doyle J: Robustness of cellular functions.
Cell 2004, 118(6):675685. PubMed Abstract  Publisher Full Text

Edwards JS, Covert M, Palsson B: Metabolic modelling of microbes: the fluxbalance approach.
Environ Microbiol 2002, 4(3):133140. PubMed Abstract  Publisher Full Text

Kuepfer L, Sauer U, Blank LM: Metabolic functions of duplicate genes in Saccharomyces cerevisiae.
Genome Res 2005, 15:14211430. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Varner JD: Largescale prediction of phenotype: concept.
Biotechnol Bioeng 2000, 69(6):664678. PubMed Abstract  Publisher Full Text

Lee B, Yen J, Yang L, Liao JC: Incorporating qualitative knowledge in enzyme kinetic models using fuzzy logic.
Biotechnol Bioeng 1999, 62(6):722729. PubMed Abstract  Publisher Full Text

Biochemical Engineering Fundamentals. McGrawHill chemical engineering series. 1986.

Heijnen JJ: Approximate kinetic formats used in metabolic network modeling.
Biotechnol Bioeng 2005, 91(5):534545. PubMed Abstract  Publisher Full Text

Schoeberl B, EichlerJonsson 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):370375. PubMed Abstract  Publisher Full Text

Mishra J, Bhalla US: Simulations of inositol phosphate metabolism and its interaction with InsP(3)mediated calcium release.
Biophys J 2002, 83(3):12981316. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stelling J, Gilles ED, Doyle F Jr 3rd: Robustness properties of circadian clock architectures.
Proc Natl Acad Sci U S A 2004, 101(36):1321013215. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Smith AE, Slepchenko BM, Schaff JC, Loew LM, Macara IG: Systems Analysis of Ran transport.
Science 2002, 295(5554):488491. PubMed Abstract  Publisher Full Text

Haunschild MD, Freisleben B, Takors R, Wiechert W: Investigating the dynamic behavior of biochemical networks using model families.
Bioinformatics 2005, 21(8):161725. PubMed Abstract  Publisher Full Text

Barkai N, Leibler S: Robustness in simple biochemical networks.
Nature 1997, 387(6636):913917. PubMed Abstract  Publisher Full Text

Bentele M, Lavrik I, Ulrich M, Stosser S, Heermann DW, Kalthoff H, Krammer PH, Eils R: Mathematical modeling reveals threshold mechanism in CD95induced apoptosis.
J Cell Biol 2004, 166(6):839851. PubMed Abstract  Publisher Full Text

Voit EO, Almeida J: Decoupling dynamical systems for pathway identification from metabolic profiles.
Bioinformatics 2004, 20(11):16701681. PubMed Abstract  Publisher Full Text

Polisetty PK, Voit EO, Gatzke EP: Identification of metabolic system parameters using global optimization methods.
Theor Biol Med Model 2006., 3(4) PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: a comparison of global optimization methods.
Genome Res 2003, 13(11):24672474. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Vandenberghe L, Boyd S: Semidefinite programming.
SIAM Rev 1996, 39:4995. Publisher Full Text

Flaherty P, Jordan MI, Arkin AP: Robust Design of Biological Experiments.
Proceedings of the Neural Information Processing Symposium 2005 2005.

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 GProtein Signaling in Yeast.
Proceedings of FOSBE (Foundations of Systems Biology and Engineering) 2005.

Boyd S, Vandenberghe L: Convex Optimization. Cambridge University Press; 2004.

Parrilo PA: Semidefinite programming relaxations for semialgebraic problems.
Math Prog 2003, 96(2, Ser B):293320. Publisher Full Text

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):17731785. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lofberg J: YALMIP: A Toolbox for Modeling and Optimization in MATLAB.

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:625653.

YALMIP [http://control.ee.ethz.ch/~joloef/yalmip.php] webcite

SeDuMi [http://sedumi.mcmaster.ca] webcite