Identification of metabolic system parameters using global optimization methods
 Pradeep K Polisetty^{1},
 Eberhard O Voit^{2} and
 Edward P Gatzke^{1}Email author
DOI: 10.1186/1742468234
© Polisetty et al; licensee BioMed Central Ltd. 2006
Received: 26 November 2005
Accepted: 27 January 2006
Published: 27 January 2006
Abstract
Background
The problem of estimating the parameters of dynamic models of complex biological systems from time series data is becoming increasingly important.
Methods and results
Particular consideration is given to metabolic systems that are formulated as Generalized Mass Action (GMA) models. The estimation problem is posed as a global optimization task, for which novel techniques can be applied to determine the best set of parameter values given the measured responses of the biological system. The challenge is that this task is nonconvex. Nonetheless, deterministic optimization techniques can be used to find a global solution that best reconciles the model parameters and measurements. Specifically, the paper employs branchandbound principles to identify the best set of model parameters from observed time course data and illustrates this method with an existing model of the fermentation pathway in Saccharomyces cerevisiae. This is a relatively simple yet representative system with five dependent states and a total of 19 unknown parameters of which the values are to be determined.
Conclusion
The efficacy of the branchandreduce algorithm is illustrated by the S. cerevisiae example. The method described in this paper is likely to be widely applicable in the dynamic modeling of metabolic networks.
Background
The past few years have witnessed an enormous increase in the availability and quality of highthroughput data characterizing the status of cells at the genomic, proteomic, metabolic and physiological levels. In most cases, these data were interpreted as simple snapshots or in a comparative setting with the goal of differentiating between normal and perturbed or diseased cells. It is now becoming feasible to use the same methods to record the status of cells over time. The resulting time series data contain enormous amounts of information on the dynamics of functioning cells. Several groups of scientists around the world have begun to develop methods for inferring from these profiles the underlying functional networks at the genomic or metabolic level. In principle, this task is a straightforward matter of defining a suitable model and estimating its structure, but numerous conceptual and computational difficulties have made the implementation of this inverse problem challenging. The difficulties fall into several categories.
First, it is necessary to select a mathematical modeling framework that is rich enough to capture the observed dynamics with sufficient accuracy but that is also structured in a fashion that allows interpretation of results beyond pure parameter estimation. For instance, if one selected a highorder polynomial and estimated parameters capturing the observed time profiles, then the resulting coefficients would not have much meaning and could hardly be translated into biological insight. Several groups [1–10] have therefore focused on Ssystem models within the modeling framework of Biochemical Systems Theory (BST), which itself has a rich history of successful analysis and application. Ssystems, and the alternative variant of Generalized Mass Action (GMA) systems, have the advantage for inverse tasks that their parameter values map essentially onetoone onto structural and regulatory features of biological networks. Thus, structure inference is reduced to the simpler task of parameter estimation. BST has been the subject of several hundred articles, reviews, books, chapters, and presentations [11–21], which permits us to review only a few features that are of particular interest here.
The second difficulty arising in the inference problem is the preparation of data and the preprocessing of the task itself. Clearly, noise in the data complicates the estimation and often leads to local minima in the search space, as well as to unwanted redundancies in inference. Furthermore, the fact that essentially any model of a dynamic biological system involves differential equations necessitates efficient integrators, because these methods may consume in excess of 95% of the time required to estimate parameters in systems of differential equations [9]. In order to reduce this computational cost, several groups have devised methods addressing some of these problems. One efficient strategy is the estimation of slopes of the profiles, which permits the replacement of differentials with the estimated slopes at many time points and consequently the conversion of systems of differential equations into larger systems of more easily computed algebraic equations [3, 9, 18, 22, 23]. A prerequisite of this method is obviously the reliable estimation of slopes, for which various smoothing methods have been proposed, including neural network smoothing [3, 9], filtering [6] and collocation methods [8]. It was also shown that decoupling of systems of n differential equations may be achieved by treating n1 of the data sets as inputs in the remaining equation [2]. In complementary approaches, the search process was simplified by making use of auxiliary information about the biological system, which was translated into constraints on the parameters that had to be estimated [3], and by priming the search process with reasonable initial guesses, which were obtained directly from the topology of the systems [24] or through various means of linearization [25].
The third difficulty of the inverse problem is the estimation of parameters itself. Because of the inherent nonlinearities, the task is hampered by trapping of the search algorithm in local minima, lack of convergence, or a convergence speed that makes inference of larger systems infeasible. In the past, this step has probably received the least attention, and the various groups addressing the inverse problem have resorted to standard methods of nonlinear regression, genetic algorithms or simulated annealing. In this article, we address this subtask of the inference problem by replacing the regression step with a global optimization algorithm. Specifically, we formulate the nonconvex estimation problem as a global optimization task that uses branchandbound principles to identify the best set of model parameters given observed time profiles. The greatest advantage of this method is that it guarantees that the optimum obtained is global within predefined bounds on the parameter search space. Note that global optimization does not guarantee that the resulting solution is unique; the method guarantees that no other points exist with a better objective function than the global solution. Multiple degenerate solution points could exist with identical unique objective function values. As an example, we estimate the parameters of a model describing the fermentation pathway in Saccharomyces cerevisiae, as described in [26]. This system has five dependent states and a total of 19 unknown parameters. It is manageable in size, yet representative of the nonlinearities typically encountered in metabolic modeling and has therefore been used for a variety of analyses in the past [18, 27–29].
Model formulation
Metabolic pathway analysis is concerned with the modeling, manipulation and optimization of biochemical systems. While valuable insights may be gained from the formulation of these systems as stoichiometric networks [30], whose functioning may be constrained as described in Flux Balance Analysis [31], it will ultimately be necessary for many purposes to formulate the processes as dynamical systems that account for detailed kinetic features, such as the regulation and modulation of enzymecatalyzed steps and transport processes. The default approach for this purpose may seem to be a model representation in the tradition of Michaelis and Menten. However, it was recognized early on that this representation is not particularly well suited for the analysis of large networks [32–34], and this led to the development of alternative methods, among which BST, Metabolic Control Analysis [35] and the "loglinear" approach [36] have received the most attention. In particular, it was shown that the Ssystem variant within BST has favorable features for the optimization of nonlinear metabolic systems [15, 37]. As an alternative to the Ssystem formulation, which may be criticized for its manner of flux aggregation, the GMA variant overcomes this issue, even though it loses some of the advantages of the Ssystem form such as its linearity at steady state [12] and its sometimes slightly higher accuracy [38, 39]. GMA representations have the advantage that they are closer to biochemical intuition and that production and degradation terms do not vanish if one of the contributing fluxes disappears, as is the case with Ssystems. GMA systems are also interesting in that they include both stoichiometric systems and Ssystems as special cases so that they allow a seamless transition from linear to fully kinetic models. Our task in this article is thus the estimation of GMA parameters from time series.
In the GMA formulation within BST, the change in each dependent pool (state variable) is described as a difference between the sums of all fluxes entering the pool and all fluxes leaving the pool. Each flux is individually linearized in logarithmic coordinates, which in Cartesian coordinates corresponds to a product of powerlaw functions that contains those and only those variables that directly affect the flux, raised to an exponent called its kinetic order. The product also contains a rate constant that determines the magnitude of the flux or speed of the process. The mathematical formulation of any GMA model is thus
where γ_{i 1},...,γ_{ ik }are rate constants corresponding to k reactions of production/consumption, and ζ_{ ijk }are kinetic orders for species i in reaction k involving species j. In cases where species j does not have any influence on a given powerlaw term, ζ_{ ijk }= 0. The number of reactions in one differential equation, k, may be different for each species. Reaction terms for consumption of one species may appear as a production term for another species. The system consists of n differential equations, representing the timedependent variables, but also contains m timeindependent variables that affect the system but are not affected by the system and are typically constant from one experiment to the next. The powerlaw terms in Eq. (1) are the result of a straightforward Taylor approximation, which is applicable to an essentially unlimited variety of underlying processes and may include different types of interactions, activation, inhibition and processes associated with dilution and growth.
It is interesting to note that every parameter of a GMA model has its unique role and interpretation. This situation is significantly different from using unstructured fitting models such as higherorder polynomials or splines. In a generic polynomial representation, every coefficient is likely to change if additional data points are used for fitting or if points are removed. Thus, except for the fact that the higherorder coefficients are associated with higher derivatives, which otherwise are not very meaningful, not much can be said about their biological role in the modeled process. In the GMA model, by contrast, every parameter has a unique meaning in the subject area of the model. Each kinetic order quantifies solely the effect that a particular variable has on a given process. For example, the first kinetic order in a later example is ζ_{121} = 0.2344 (see Figure 7). Thus, it uniquely describes the effect of metabolite X_{2} on the first production process of X_{1}. The effect is inhibitory, which is indicated by the negative sign, and is only moderately strong, which is reflected in the small magnitude of the parameter. In this fashion, there is a onetoone relationship between kinetic orders and structural features of the model. The interpretability of the parameters can also been seen from a different point of view: in principle, each kinetic order can be obtained directly from local information on the system. Namely, if it possible to vary X _{ j }while keeping all other variables constant, and to measure the consequent changes in the production of X_{ j }, then the slope of the production process as a function of X_{ j }, in log space, is exactly the kinetic order in question. This type of interpretation is usually not possible in global fitting models such as highorder polynomials. The constant multipliers are rate constants that, as in elemental chemical kinetics, quantify the turnover rate in each process and are always nonnegative. Their magnitudes depend on the scales (time, concentration, etc.) of the modeled system.
The Ssystem representation within BST is formally a special case of a GMA system with at most one positive and one negative term. The correspondence between the biological network and the mathematical representation is slightly different in an Ssystem because all fluxes entering a pool are collectively linearized in logarithmic coordinates and the same is done with all fluxes leaving the pool. Thus, each Ssystem equation has at most one influx and one efflux term and thus reads
where α_{ i }, β_{ i }are the rate constants and g_{ ij }, h_{ ij }are kinetic orders. These Ssystem parameters are directly tied to the corresponding GMA system through constraints [18].
Both forms are nonlinear and rich enough to capture any dynamic behavior that can be represented by any set of ordinary differential equations [38]. If set up as alternative descriptions of the same biological system, the two are equivalent at one operating point of choice and typically differ if the system deviates from this point, though the differences are often small in realistic situations [18].
Working with metabolic models consists of three phases: Model design, model analysis and model application. The present work focuses on the first phase of model design. This step is usually executed by assembling a topological map of the phenomenon of interest based on biological knowledge. Kinetic orders and rate constants are then estimated from measured or published kinetic information. In the case of Ssystems, it would also be feasible to use steadystate data from multiple experiments with different values of independent variables. Because the steadystate equations of Ssystems are linear (in logarithmic coordinates), such data would allow use of a simple matrix inversion leading to optimal parameters, a pseudoinverse method or a linear programming approach [37]. In the alternative GMA approach, logarithmic transformation does not completely transform the estimation problem into a linear formulation, thus requiring nonlinear methods.
In addition to this traditional bottomup approach, a topdown approach is becoming increasingly feasible. This complementary approach is based on time profiles that are used for the determination of system parameters through some type of estimation and their interpretation in terms of structural and regulatory information. We describe in the following how this estimation is facilitated with a branchandbound method, which has not previously been used for this purpose.
Optimization formulation
The goal of model identification is to estimate the "best" set of parameter values, which minimizes the error between the process data and the model response. This parameter estimation problem can be formulated as a nonconvex nonlinear optimization problem and can therefore be solved using global optimization techniques. The parameters to be estimated in GMA systems are rate constants γ and kinetic orders ζ as shown in Equation 1 leading to the following formulation:
The parameter r denotes the rnorm considered for minimization. It usually takes the values 1 (min sum absolute), 2 (min sum square), or ∞ (min max error). P is the number of data points sampled at times t, h_{ i }are the nonlinear rate expressions from Eq. (1) that define the production and consumption rates for species i given vectors of model parameters γ and ζ, e_{ i }(t) are the errors associated with each constraint equation for species i at time t, and n is the number of dependent variables. In the formulation described above, the objective function is linear because the maximum absolute error is minimized. The nonconvexity arises from the equality constraints, which are nonlinear. It is useful to split these nonlinear equality constraints into two inequality constraints, at least one of which will be nonconvex.
Following strategies proposed in the literature [6, 8, 9], it is possible to smooth the raw time profiles, which subsequently allows the computation of slopes at many data points and thus the replacement of differentials on the lefthand side of Equation 1 with estimated slopes. Thus, assuming that the rate of change of each species is obtainable at each desired time point, (t), and the values of the concentrations X_{ i }(t) at time i are known, the optimization task in Eq. (3) can be formulated as a general nonconvex nonlinear programming problem in the form shown in Eq. (4):
where the vector x ∈ R^{ N }is a vector of N unknowns including the error terms and unknown parameters and f(x), g_{ k }(x) : R^{ N }→ R^{1}. Here, the index k represents constraint k of m total constraints.
This formulation is rather general in that the functions f(x) and g_{ k }(x) may be nonlinear and nonconvex. In particular, the formulation allows the estimation of globally optimal GMA systems, given time profiles. Deterministic methods for global optimizations of this type depend on the generation of convex function relaxations of the nonconvex nonlinear functions. Numerous methods have been proposed for constructing such relaxations. For this work, we use a reformulation method for factorable expressions [40]. This method converts the original factorable nonconvex nonlinear problem into an equivalent form through the introduction of new variables z_{ ik }for every product of powerlaw terms at measurement sample time t in the system [15]:
Note that the species concentration X_{ i }(t) in the preceding equation is assumed to be known from the observed data. If it is not, it may be obtained through interpolation from a preparatory smoothing of the observed data [3, 9]. If no information at all is available on the variable, the GMA representation should probably be reduced in complexity [18]. As an example for such a reduction, suppose that X_{1} is converted into X_{2} and X_{2} is converted into X_{3}. If X_{2} is not observable, then one would probably formulate the system without X_{2}, and make X_{3} a function of X_{1}. Since the mathematics underlying the GMA representation is directly based on Taylor's theorem, X_{3} then simply becomes a powerlaw term containing X_{1}. Such omissions of variables have been discussed in the literature [31, 41].
For simplicity of discussion, we assume that a complete data set is available. The problem specified in Eq. (3) is then of the form:
where s_{ ik }for reaction k of production or consumption of species i is either +1 or 1. The variables in this reformulated problem can be collected into the vectors z, e, γ, and ζ
To obtain efficacious solutions, the problem must be formulated such that it contains only linear and "simple" nonlinear constraint functions for which one can construct relaxations that use the convex envelopes already known for simple algebraic functions. The nonlinear equality constraints in the original GMA formulation become simple (linear) sums when they are expressed in the new variables, z. Furthermore, taking the logarithm of each definition for z from Eq. (5) leads to a new set of linear equations and simple nonlinear equations. The connection between these two sets is a set of simple logarithmic constraints. To streamline notation, the additional new variables are defined as: w_{ ik }= ln(z_{ ik }) and Γ_{ ik }= ln(γ_{ ik }). Thus, we obtain a sum of logarithmic functions for each term z_{ ik }(t), namely:
For the purpose of convexification, one can omit the logarithm of the rate constant because it simply shifts the optimal solution by a constant. This results in the following formulation:
As s_{ ij }, X_{ i }(t), and (t) are known at time values t, the only unknowns are w Γ, z, ζ and e. Values for γ can be determined easily when the solution value for Γ is found. All constraints other than w_{ ik }(t) = ln(z_{ ik }(t)) are linear.
As mentioned earlier, the parameter r in Eq (3) denotes the rnorm considered for minimization. In the case with r = 1, the sum absolute value of the error, which is linear, is desired for minimization. In the case of r = 2, the objective function is nonlinear, but at least convex. This nonlinearity increases the complexity of the task during the generation of problem relaxations, requiring additional variables and constraints in the linear relaxation. It is well known that when minimizing the absolute value of some variable, constraints of the form f(x) = e can be written as two inequality constraints, f(x) ≤ e, and f(x) ≤ e. The formulation then may be written as follows:
Note that the variables Γ_{ ik }, ζ_{ ijk }and e_{ i }(t) can be consolidated into the vector y. The variables in this vector y only appear in linear constraints, while w and z are related through a simple nonlinear expression. Thus, in comparison with the formulation of Equation 4, the variables x include w, z and y. The objective function f(x) and many of the constraint functions g_{ k }(x) are linear. The objective function (1) represents a sum of the absolute value of errors on the balance equations. The balance equations in (2) and (3) relate the rate of change, (t) of the species at time t, to the individual reaction rates that consume or produce that species, z_{ ik }(t), as well as the absolute error for that equation at that point in time. Constraints (4) and (5) result from the transformation of the powerlaw expressions for the rate equations.
As bounds are known for the concave nonconvex function, the secant can be used as a linear underestimation function. Multiple outer approximation linearizations can be used as linear overestimation functions. The intersection of these linear constraints is a relaxation of the original nonlinear function.
where all linear inequality constraints are represented by A_{1}[w^{ T }z^{ T }y^{ T }]^{ T }≤ b_{1}, and A_{2}[w^{ T }z^{ T }y^{ T }]^{ T }= b_{2} defines the new linear constraints obtained from reformulation, while w = η(z) provides the relationship between w and z. Bounds on w are determined from the bounds on z using interval methods. Note that η consists of simple nonlinear (logarithmic) terms and that the formulation in Equation 10 is the same as that in Equation 4 with the vector of unknowns x made up of w, z and y. Additionally, f(x) and many of the g_{ k }(x) are linear functions in Eq. (4). Note that any equality constraint g_{ k }(x) = 0 can be written equivalently with two inequalities as 0 ≤ g_{ k }(x) ≤ 0 or g_{ k }(x) ≤ 0, g_{ k }(x) ≤ 0.
Convex relaxations for this task are constructed with DAEPACK [42, 43], an automated code generation tool. The advantage of using the DAEPACK tool for this purpose is that it can be applied directly to legacy models coded in standard FORTRAN. The convex relaxations can be denoted as:
A linearization strategy [42, 44] is then used to generate a Linear Programming (LP) relaxation of the convex NLP created using DAEPACK. The resulting LP is of the form:
where A_{3} [w^{ T }z^{ T }y^{ T }]^{ T }≤ b_{3} expresses the new linear constraints resulting from the linearization process as illustrated in Figure 1. This linearization technique is ideal because it yields a linear program for which robust solvers exist (e.g. ILOG CPLEX 8.0 [45] and the IBM OSL library [46]). Note that A_{3}, b_{3}, w^{ l }and w^{ u }are updated as z^{ l }and z^{ u }change in the spatial branchandbound algorithm.
Solution methodology
1. If the relaxed problem associated with the partition is infeasible, adding additional constraints will not make the problem feasible. The partition itself is infeasible and hence can be removed from further consideration.
2. If the objective function value of the relaxed problem associated with the current partition is greater or equal to the best solution found so far, then the partition can be removed from further consideration.
Any feasible solution to the original problem may serve as an upper bound for the global solution. The algorithm terminates when the lower bounds for all partitions either exceed or are sufficiently close to the best upper bound. At this point, a global optimum has been determined within the originally preset bounds on the parameter search space. This global optimum is the best value of the objective function. It is noted that multiple points in the parameter space may lead to equivalent values of the objective function.
For the optimization problem shown in Eq. (3), a branchandreduce method [48] was implemented. This is an extension of the traditional branchandbound method with bound tightening techniques for accelerating the convergence of the algorithm. Within this branchandreduce algorithm, infeasible or suboptimal parts of the feasible region are eliminated by using range reduction techniques such as optimalitybased and feasibilitybased range reduction tests [48–50] or interval analysis techniques [51]. These techniques render tighter variable bounds for a given partition in the search tree, thereby leading to more rapid convergence.
Didactic example
Time series data of the states and slopes for the didactic example
Data Points  X _{1}  X _{2}  X _{3} 




1  5.00e1  5.00e1  1.00  2.66  1.21e  3.81 
2  2.91e1  5.86e1  6.498e1  1.52  5.59e1  3.17 
3  1.96e1  6.22e1  3.73e1  3.82e1  2.08e1  2.34 
4  2.20e1  6.42e1  1.90e1  9.11e1  2.69e1  1.28 
5  3.65e1  6.87e1  1.16e1  1.69  6.39e1  2.59e1 
6  4.88e1  7.63e1  1.15e1  6.24e1  8.42e1  1.08e1 
7  5.04e1  8.46e1  1.24e1  1.66e1  7.85e1  5.54e2 
8  4.75e1  9.18e1  1.25e1  3.31e1  6.49e1  3.05e2 
9  4.45e1  9.76e1  1.197e1  2.52e1  5.25e1  6.29e2 
10  4.26e1  1.02e1  1.14e1  1.43e1  4.35e1  5.59e2 
11  4.15e1  1.06e1  1.09e1  8.07e2  3.73e1  3.81e2 
12  4.08e1  1.099e1  1.06e1  5.99e2  3.26e1  2.51e2 
The parameters were estimated by means of the branchandreduce global optimization algorithm for various scenarios that differed in both the initial guesses and the bounds for the parameter values. The initial guesses were selected in different ways. In the first series of experiments they were randomly chosen (with uniform distribution) within a predefined range between lower bounds and upper bounds on the parameters. Computational results are presented here with the initial guesses selected as the lower bounds on the parameter search space. Initial parameter guesses can be based on collective experience with GMA and Ssystems as described in [18]. Any reasonable initial solution may serve as the initial upper bound on the solution. For our illustration, the upper and lower bounds were selected at 10%, 100%, 200% and 500% around the true value. For instance, in the latter case, they were set at 500% and +500% of the true parameter values, which we knew for this didactic example. The initial bounds on the parameters can thus be computed as:
[k_{ true } 500% × k_{ true }, k_{ true }+ 500% × k_{ true }] (13)
where k_{ true }is the true parameter value. This technique leads to search regions in parameter space centered on the nominal values. However, in the branchandreduce algorithm, various range reduction techniques are implemented prior to the global search for the solution. For instance, one may derive tighter variable bounds on a given partition that are not centered around the nominal parameter values. In a realistic situation the true values are of course unknown, but considerable experience has been amassed throughout the years suggesting natural default values. We also included a 0% interval as a check that the true solution was recouped.
Estimated parameters using local and global optimization algorithms for the simple GMA system
Actual Parameters  Estimated parameters with varying bounds  

0 %  10 %  100 %  200 %  500 %  
Local  Global  Local  Global  Local  Global  Local  Global  Local  Global  
0.8  0.8  0.8  0.8  0.8  0.799  0.799  0.8  0.8  0.8  0.8 
1.0  1.0  1.0  1.0  1.0  1.0  1.0  3.0  1.0  1.0  1.0 
1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 
3.0  3.0  3.0  3.0  3.0  2.99  2.99  1.09  3.0  3.0  3.0 
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.45  0.5  0.5  0.5 
0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.3  0.1  0.1  0.1 
2.0  2.0  2.0  2.0  2.0  2.0  2.0  1.019  2.0  2.0  1.999 
0.75  0.75  0.75  0.75  0.75  0.75  0.75  0.40  0.75  0.75  0.749 
0.2  0.2  0.2  0.2  0.2  0.199  0.199  0.2  0.2  0.2  0.2 
1.5  1.5  1.5  1.5  1.5  1.499  1.499  1.5  1.5  1.5  1.5 
0.5  0.5  0.5  0.5  0.5  0.499  0.499  0.5  0.499  0.5  0.499 
5.0  5.0  5.0  5.0  5.0  5.0  5.0  6.074  5.0  5.0  4.999 
0.5  0.5  0.5  0.5  0.5  0.499  0.499  0.611  0.499  0.5  0.5 
Total time taken for global solution, number of partitions created during the BranchandReduce algorithm, nonconvex and convex problems solved, and objective function values for both local and global solution methods for the small GMA system.
Time for Global Sol (sec)  No of Partitions  Nonconvex Problems  Convex Problems  Obj fun Value (Global)  Obj fun Value (Local)  

0 %  0.252  1  1  1  0.0  0.0 
10 %  0.662  1  1  1  0.0  0.0 
100 %  1.453  1  1  1  0.0  0.0 
200 %  12.005  11  13  11  0.0  0.917 
300 %  14.412  13  15  13  0.0  0.884 
400 %  10.372  9  11  9  0.0  0.867 
500 %  20.042  15  17  15  0.0  0.846 
The numerical results demonstrate that both the local and global solvers give the same solution when the parameter space is small and the bounds are tight. When the bounds are increased beyond 100%, the local solver may fail to find the true parameters, while the global solver still succeeds. For instance, one notes that the local search with 200% bounds yielded a different solution, which however was inferior, as indicated by the value of the objective function.
Case study
According to Galazzo and Bailey [52] and Curto et al. [26], the observed concentrations (mM) of the dependent variables at steady state are: X_{1}(G_{ In }) – Internal Glucose = 0.0346, X_{2}(G6P) – Glucose6phosphate = 1.011, X_{3}(FDP) – Fructose1,6diphosphate = 9.1876, X_{4}(PEP) – Phosphoenolpyruvate = 0.0095, and X_{5} – Adenosine triphosphate (ATP) = 1.1278. The values of the independent variables (mM min^{1}) are: X_{6} – Glucose uptake = 19.7, X_{7} – Hexokinase = 68.5, X_{8} – Phosphofructokinase = 31.7, X_{9} – Glyceraldehyde3phosphate dehydrogenase = 49.9, X_{10} – Pyruvate kinase = 3,440, X_{11} – Polysaccharide production (glycogen + trehalose) = 14.31, X_{12} – Glycerol production = 203, X_{13} – ATPase = 25.1, and X_{14} – NAD^{+}/NADH ratio = 0.042.
Time series data of the states and slopes for the anaerobic fermentation case study (In every cell, the first value represents the state and the second value represents the slope)
Data Points  X_{6} External Glucose  X _{1} (G_{In})  X _{2}
(G6P)  X _{3} (FDP)  X _{4} (PEP)  X _{5} (ATP) 

1  19.7  3.46e2  1.011  9.188  9.53e3  1.1278 
1.78e15  1.10e13  1.97e9  3.93e9  3.93e9  
2  19.9  3.46e2  1.011  9.188  9.53e3  1.1278 
1.62e1  1.10e13  19.7e9  3.93e9  3.93e9  
3  19.9  3.498e2  1.017  9.196  9.51e3  1.1175 
7.399e5  1.021e2  1.38e1  9.79e6  4.68e2  
4  19.9  3.498e2  1.017  9.21  9.52e3  1.1167 
3.69e5  2.04e3  1.32e1  2.04e4  1.76e2  
5  19.9  3.497e2  1.018  9.22  9.54e3  1.1193 
1.177e4  8.28e3  1.06e1  2.16e4  2.91e2  
6  19.5  3.495e2  1.019  9.232  9.56e3  1.1221 
3.24e1  9.67e3  8.59e2  1.821e4  2.568e2  
7  19.5  3.41e2  1.0075  9.22  9.63e3  1.1449 
5.68e5  1.36e2  2.03e1  1.595e4  1.131e1  
8  19.5  3.41e2  1.008  9.201  9.616e3  1.148 
1.34e5  2.35e3  2.03e1  2.85e4  1.79e2  
9  19.5  3.412  1.0072  9.182  9.58e3  1.145 
1.47e4  1.02e2  1.63e1  3.26e4  4.32e2  
10  19.7  3.41e2  1.006  9.168  9.55e3  1.1406 
1.62e1  1.38e2  1.32e1  2.77e4  3.91e2 
The parameters were estimated with the MINOS nonconvex NLP solver for local searches and the branchandreduce algorithm for the global searches. As in the didactic example, we explored various scenarios that again differed with respect to the initial guesses for the parameters as well as the bounds on the parameter space to be searched. The lower and upper bounds on the parameters were selected at 100%, 200%, 300 %, and 500% around the true value. The initial guesses for the parameters were selected to be the lower bounds of the parameters search space.
It is known that the performance of any MINOS NLP solver is significantly affected by the specified tolerances. The convergence tolerances such as row, optimality and feasibility tolerances were set to 10^{5}. The maximum number of iterations was set to 5000, and the maximum number of major and minor iterations between successive linearizations of nonlinear constraints was set to a value of 60. During the branchandreduce algorithm, only the parameters (γ and p) were considered for branching. The ratios of the difference in the current bounds and the difference in the original bounds for all variables (parameters) were computed. The particular variable with the worst (largest) ratio was then selected for branching. The algorithm was implemented using bestfirst search strategy and was able to guarantee the global solution within an ε = 0.0001 tolerance.
Estimated parameters for anaerobic fermentation pathway using BranchandReduce algorithm
Actual Parameters (θ)  Estimated parameters with varying bounds  

100 %  300 %  500 %  
Local  Global  Local  Global  Local  Global  
0.8122  0.8112  0.8122  0.8112  0.8122  0.8112  0.8122 
0.2344  0.2393  0.2344  0.2396  0.2344  0.2398  0.2344 
2.8632  2.8557  2.8632  2.8521  2.8632  2.8485  2.8632 
0.7464  0.7460  0.7464  0.7456  0.7464  0.7452  0.7464 
0.0243  0.0259  0.0243  0.0256  0.0243  0.0253  0.0243 
0.5232  0.5221  0.5232  0.5214  0.5232  0.5206  0.5232 
0.7318  0.7364  0.7318  0.7373  0.7318  0.7382  0.7318 
0.3941  0.3941  0.3941  0.3949  0.3941  0.3954  0.3941 
0.0009  0.0018  0.0009  0.0036  0.0009  0.0054  0.0009 
8.6107  0.0000  8.6107  0.0000  8.6107  0.0000  8.6107 
0.011  0.0109  0.011  0.0109  0.011  0.0109  0.011 
0.6159  0.6169  0.6159  0.6177  0.6159  0.6185  0.6159 
0.1308  0.1302  0.1308  0.1305  0.1308  0.1308  0.1308 
0.04725  0.0337  0.034  0.0174  0.0467  0.0089  0.0464 
0.05  0.1  0.1  0.2  0.0517  0.3000  0.0528 
0.533  0.4852  0.486  0.3915  0.531  0.2977  0.5303 
0.0822  0.0638  0.063  0.026  0.0816  0.011  0.0811 
1.0  0.9975  1.0000  0.9937  1.0000  0.9899  1.0000 
1.0  0.9979  1.0000  1.0013  1.0000  1.0047  1.0000 
Result of estimation analysis with the anaerobic fermentation pathway model. Total time required for global solution, number of partitions created during the BranchandReduce search, number of nonconvex and convex problems solved, and the value of the objective function for both local and global solution methods.
Time for Global Sol (sec)  No of Partitions  Nonconvex Problems  Convex Problems  Obj fun Value (Global)  Obj fun Value (Local)  

100 %  38.13  27  15  27  0.0  1.36 × 10^{4} 
200 %  19.28  9  6  9  0.0  1.34 × 10^{4} 
300 %  56.98  27  15  27  1.61 × 10^{7}  1.31 × 10^{4} 
400 %  58.39  25  14  25  0.0  1.31 × 10^{4} 
500 %  57.15  25  14  25  3.51 × 10^{8}  1.25 × 10^{4} 
Conclusion
The first step of any modeling effort is the definition of symbolic equations and the numerical definition of their parameter values. In the past, the latter has typically been accomplished by reformulating literature information on local processes, such as enzymecatalyzed reaction steps, within the biological system of interest. Highthroughput data are beginning to change this situation dramatically. Instead of working from the bottom up, it is becoming possible to infer the structure and regulation of biological systems from measured time profiles. An important component of this inference is the efficient estimation of parameters, which in the case of GMA and Ssystems within BST can directly be interpreted as biological features, thereby potentially yielding novel insights.
We have shown here that the parameter estimation task may be posed as a nonconvex nonlinear optimization problem. Consequently, deterministic global optimization techniques can be applied, and among them the branchandreduce algorithm appears to be a very suitable choice, because it is guaranteed to find the global solution. This is in stark contrast to local solvers, such as nonlinear regression algorithms, which may not be able to converge to the global solution when the parameter search space is large or the error surface is ragged. As a proof of concept, we demonstrated the features of the global search method with a didactic example and a simple yet realistic biological case study, which has been used by others for the exploration of new methods [18, 26–29]. In order not to cloud the demonstration of the functioning of the branchandreduce search algorithm, both examples were selected free of experimental error. While this may seem somewhat artificial, the omission of noise was intentional (a) in order to see how the method performs under ideal conditions and (b) because a preparatory smoothing step can significantly reduce if not eliminate noise. Nonetheless, it is now necessary to consider data that are corrupted by noise, and to compare the efficacy of the branchandreduce method with alternate methods such as nonlinear regression and genetic algorithms.
For the limited range of illustrative examples shown here, the computational results demonstrate that the branchandreduce algorithm is fast and reliable. Of course, it is difficult to judge in general to what degree this algorithm, combined with estimation of species concentration rate of change, outperforms other methods, but it seems clear that direct parameter estimation in systems of differential equations can be extremely time consuming. As a case in point, Kikuchi et al. [4] used a genetic algorithm to identify the parameters of a fivevariable Ssystem from noisefree data and showed that every cycle of the genetic algorithm took about 10 hours on a cluster of over 1,000 CPUs and that seven cycles were needed to complete the identification. Voit and Almeida [9] demonstrated that this computation time may be reduced drastically if the estimation of slopes, as we used here, is employed to replace the differential equations with sets of nonlinear algebraic equations. Still they found the subsequent regression to be slow and not necessarily reliable, because of convergence to local minima or no convergence at all. Irrespective of computational speed, the branchandreduce algorithm employed here was shown to be successful in rather quickly finding global solutions within an ε tolerance, which is a great advance, because local search methods cannot guarantee convergence to the true solution.
Declarations
Acknowledgements
The Authors would like to acknowledge financial support from the South Carolina Collaborative Grants Program (E.P. Gatzke, PI), from an NSF Quantitative Systems Biotechnology (BES0120288; E.O. Voit, PI) grant, and from a National Heart, Lung and Blood Institute Proteomics Initiative (Contract N01HV28181; D. Knapp, PI). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the sponsoring institutions.
Authors’ Affiliations
References
 Seatzu C: A Fitting Based Method for Parameter Estimation in SSystems. Dynamic Systems Applications. 2000, 9 (1): 7798.Google Scholar
 Maki Y, Tominaga D, Okamoto M, Watanabe S, Eguchi Y: Development of a system for the inference of large scale genetic networks. Pac Symp Biocomput. 2001, 6: 446458.Google Scholar
 Almeida J, Voit EO: Neural networkbased parameter estimation in complex biochemical systems. Genome Informatics. 2003, 14: 114123.PubMedGoogle Scholar
 Kikuchi S, Tominaga D, Arita M, Takahashi K, Tomita M: Dynamic modeling of genetic networks using genetic algorithm and Ssystem. Bioinformatics. 2003, 19: 643650. 10.1093/bioinformatics/btg027.View ArticlePubMedGoogle Scholar
 Lall R, Rutes A, Santos H, Almeida J, Voit EO: A New Approach to Parameter Estimation using Ssystems: Modeling the Glycolytic Pathway of Lactococcus lactis. Georgia TechUGA Conference on Bioinformatics, Atlanta, GA. 2003 November 13–16
 Borges CCH, Voit EO, Almeida J: Signal extraction for numerical decoupling of Ssystems. International Conference on Molecular Systems Biology (ICMSB '04), Tahoe, CA. 2004 August 21–25
 Voit EO, Marino S, Lall R: Challenges for the identification of metabolic pathways from time series data. In Silico Biology. 2004, 5: 110.Google Scholar
 Wang FS, Tsai KY: A global/local optimization approach for dynamic system modeling of biological networks. Bioinformatic. 2005, 21: 11801188.View ArticleGoogle Scholar
 Voit EO, Almeida J: Decoupling dynamical systems for pathway identification from metabolic profiles. Bioinformatics. 2004, 20 (11): 16701681. 10.1093/bioinformatics/bth140.View ArticlePubMedGoogle Scholar
 Naval P, Gonzalez O, Mendoza E, Sison L: Heuristic Parameter Estimation Methods for SSystem Models of Biochemical Networks. Second Humanoid, Nanotechnology, Information Technology, Communications and Control, Environment and Management (HNICEM), Manila, Philippines. 2005Google Scholar
 Savageau MA: Biochemical Systems Analysis, I. Some Mathematical Properties of the Rate Law for the Component Enzymatic Reactions. J Theor Biol. 1969, 25: 365369.View ArticlePubMedGoogle Scholar
 Savageau MA: Biochemical Systems Analysis, II. The steadystate solutions for an npool system using a powerlaw approximation. J Theor Biol. 1969, 25: 370379.View ArticlePubMedGoogle Scholar
 Savageau MA: Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology. AddisonWesley, Massachusetts. 1976Google Scholar
 Savageau MA, Voit EO: Recasting Nonlinear Differential Equations as SSystems: A Canonical nonlinear Form. Math Biosci. 1987, 87: 83115. 10.1016/00255564(87)900356.View ArticleGoogle Scholar
 Torres NV, Voit EO: Pathway Analysis and Optimization in Metabolic Engineering. 2002, Cambridge University Press, Cambridge, UKView ArticleGoogle Scholar
 Voit EO: Canonical Nonlinear Modeling: SSystem Approach to Understanding Complexity. 1991, Van Nostrand Reinhold, New YorkGoogle Scholar
 Voit EO: Canonical modeling: A review of concepts with emphasis on environmental health. Environmental Health Perspectives, Mathematical Modeling in Environmental Health Studies. 2000, 108 (Suppl 5): 895909.View ArticleGoogle Scholar
 Voit EO: Computational Analysis of Biochemical Systems. 2000, Cambridge University Press, New YorkGoogle Scholar
 Voit EO: Metabolic modeling : a tool of drug discovery in the postgenomic era. Drug Discov Today. 2002, 7: 621628. 10.1016/S13596446(02)022808.View ArticlePubMedGoogle Scholar
 Voit EO: The Dawn of a New Era of Metabolic System Analysis. Drug Discovery Today BioSilico. 2004, 2 (5): 182189. 10.1016/S17418364(04)024199.View ArticleGoogle Scholar
 Voit EO, Torres NV: Canonical modeling of complex pathways in biotechnology. Recent Res Devel in Biotech and Bioeng. 1998, 1: 321341.Google Scholar
 Voit EO, Almeida J: Dynamic Profiling and Canonical Modeling: Powerful Partners in Metabolic Pathway Identification. Metabolite Profiling : Its Role in Biomaker Discovery and Gene Function Analysis. 2003, Kluwer Academic Publishing, Dordrecht, The NetherlandsGoogle Scholar
 Voit EO, Savageau MA: Powerlaw approach to modeling biological systems, 3. Methods of Analysis. J Ferment Technol. 1982, 60 (3): 233241.Google Scholar
 Torralba AS, Yu K, Shen P, Oefner PJ, Ross J: Experimental test of a method for determining casual connectivities of species in reactions. Proc Natl Acad Sci. 2003, 100: 14941498. 10.1073/pnas.262790699.PubMed CentralView ArticlePubMedGoogle Scholar
 Veflingstad SR, Almeida J, Voit EO: Priming nonlinear searches for pathway identification. BMC Theoretical Biology and Medical Modelling. 2004, 1: 810.1186/1742468218.View ArticleGoogle Scholar
 Curto R, Sorribas A, Cascante M: Comparative characterization of the fermentation pathway of Saccharomyces cerevisiae using biochemical systems theory and metabolic control analysis. Model definition and nomenclature. Math Biosc. 1995, 130: 2550. 10.1016/00255564(94)00092E.View ArticleGoogle Scholar
 Cascante M, Curto R, Sorribas A: Comparative characterization of the fermentation pathway of Saccharomyces cerevisiae using biochemical systems theory and metabolic control analysis. Steadystate analysis. Math Biosc. 1995, 130: 5169. 10.1016/00255564(94)00093F.View ArticleGoogle Scholar
 Sorribas A, Curto R, Cascante M: Comparative characterization of the fermentation pathway of Saccharomyces cerevisiae using biochemical systems theory and metabolic control analysis. Model validation and dynamic behavior. Math Biosc. 1995, 130: 7184. 10.1016/00255564(94)00094G.View ArticleGoogle Scholar
 Torres NV, Voit EO, GlexAlcon C, Rodriguez F: An indirect optimization for biochemical systems: Description of method and application to the maximization of the rate of ethanol, glyceerol, and carbohydrate production in Saccharomyces cerevisiae. Biotech Bioeng. 1997, 55 (5): 758772. 10.1002/(SICI)10970290(19970905)55:5<758::AIDBIT6>3.0.CO;2A.View ArticleGoogle Scholar
 Stephanopoulos G, Aristidou A, Nielsen J: Metabolic Engineering. Principles and Methadologies. 1998, Avademic Press, San Diego, CAGoogle Scholar
 Reed JL, Palsson BO: Thirteen years of building constraintbased in silico models of Escherichia coli. J Bacteriol. 2003, 185 (9): 26929. 10.1128/JB.185.9.26922699.2003.PubMed CentralView ArticlePubMedGoogle Scholar
 Heinrich R, Rapoport TA: A linear steadystate treatment of enzymatic chains: General properties, control and effector strength. Eur J Biochem. 1974, 42: 8995. 10.1111/j.14321033.1974.tb03318.x.View ArticlePubMedGoogle Scholar
 Hill CM, Waight RD, Bardsley WG: Does any enzyme follow the MichaelisMenten equation?. Molec Cell Biochem. 1977, 15: 173178. 10.1007/BF01734107.View ArticlePubMedGoogle Scholar
 Savageau M: Enzyme kinetics in vitro and in vivo: MichaelisMenten revisited. 1995, IAI Press Inc, Greenwich, Connecticut, 4: 93146.Google Scholar
 Fell DA: Understanding the Control of Metabolism. 1977, Portland Press, LondonGoogle Scholar
 Hatzimanikatis V, Bailey JE: MCA has more to say. J Theor Biol. 1996, 182: 233242. 10.1006/jtbi.1996.0160.View ArticlePubMedGoogle Scholar
 Voit EO: Optimization in integrated biochemical systems. Biotechn Bioengin. 1992, 40: 572582. 10.1002/bit.260400504.View ArticleGoogle Scholar
 Voit EO, Savageau M: Accuracy of alternative representations for integrated biochemical systems. Biochemistry. 1987, 26: 68696880. 10.1021/bi00395a042.View ArticlePubMedGoogle Scholar
 Sorribas A, Savageau M: Strategies for representing metabolic pathways within biochemical systems theory: Reversible Pathways. Math Biosci. 1989, 94: 239269. 10.1016/00255564(89)900667.View ArticlePubMedGoogle Scholar
 McCormick GP: Computability of Global Solutions to Factorable Nonconvex Programs: Part I – Convex Underestimating Problems. Mathematical Programming. 1976, 10: 147175. 10.1007/BF01580665.View ArticleGoogle Scholar
 Curto R, Voit EO, Sorribas A, Cascante M: Mathematical models of purine metabolism in man. Math Biosc. 1998, 151: 149. 10.1016/S00255564(98)100019.View ArticleGoogle Scholar
 Gatzke EP, Tolsma JE, Barton PI: Construction of Convex Function Relaxations Using Automated Code Generation Techniques. Optimization and Engineering. 2002, 3 (3): 305326. 10.1023/A:1021095211251.View ArticleGoogle Scholar
 Tolsma JE, Barton PI: DAEPACK: An Open Modeling Environment for Legacy Code. Ind Eng Chem Res. 2000, 39 (6): 18261839. 10.1021/ie990734o.View ArticleGoogle Scholar
 Tawarmalani M, Sahinidis NV: Global Optimization of Mixed Integer Nonlinear Programs: A Theoretical and Computational Study. 2000, Technical report, University of IllinoisGoogle Scholar
 ILOG: ILOG CPLEX 8.1: User's Manual. 2002, Mountain View, CAGoogle Scholar
 IBMOSL: IBM Optimization Solutions and Library Linear Programming Solutions. Technical report I B M. 1997Google Scholar
 Falk JE, Soland RM: An Algorithm for Separable Nonconvex Programming Problems. Management Science. 1969, 15 (9): 550569.View ArticleGoogle Scholar
 Ryoo HS, Sahinidis NV: Global Optimization of Nonconvex NLPS and MINLPs with Application to Process Design. Comput Chem Eng. 1995, 19 (5): 551566. 10.1016/00981354(94)000978.View ArticleGoogle Scholar
 Adjiman CS, Androulakis IP, Floudas CA: Global Optimization of MixedInteger Nonlinear Problems. AIChE J. 2000, 46 (9): 17691797. 10.1002/aic.690460908.View ArticleGoogle Scholar
 Smith EMB: On the Optimal Design of Continuous Processes. PhD thesis. 1996, Imperial College, LondonGoogle Scholar
 Moore RE: Methods and Applications of Interval Analysis. SIAM Philadelphia. 1979Google Scholar
 Galazzo JL, Bailey JE: Fermentation pathway kinetics and metabolic flux control in suspended and immobilized Saccharomyces cerevisiae. Enzyme Microb Technol. 1990, 12: 162172. 10.1016/01410229(90)90033M.View ArticleGoogle Scholar
 Voit EO: Computational Analysis of Biochemical Systems. 2000, Cambridge University PressGoogle Scholar
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.