Hybrid optimization for ^{13}C metabolic flux analysis using systems parametrized by compactification
 Tae Hoon Yang^{1, 2}Email author,
 Oliver Frick^{2} and
 Elmar Heinzle^{2}
DOI: 10.1186/17520509229
© Yang et al; licensee BioMed Central Ltd. 2008
Received: 08 October 2007
Accepted: 26 March 2008
Published: 26 March 2008
Abstract
Background
The importance and power of isotopebased metabolic flux analysis and its contribution to understanding the metabolic network is increasingly recognized. Its application is, however, still limited partly due to computational inefficiency. ^{13}C metabolic flux analysis aims to compute in vivo metabolic fluxes in terms of metabolite balancing extended by carbon isotopomer balances and involves a nonlinear leastsquares problem. To solve the problem more efficiently, improved numerical optimization techniques are necessary.
Results
For flux computation, we developed a gradientbased hybrid optimization algorithm. Here, independent flux variables were compactified into [0, 1)ranged variables using a single transformation rule. The compactified parameters could be discriminated between nonidentifiable and identifiable variables after model linearization. The developed hybrid algorithm was applied to the central metabolism of Bacillus subtilis with only succinate and glutamate as carbon sources. This creates difficulties caused by symmetry of succinate leading to limited introduction of ^{13}C labeling information into the system. The algorithm was found to be superior to its parent algorithms and to global optimization methods both in accuracy and speed. The hybrid optimization with tolerance adjustment quickly converged to the minimum with close to zero deviation and exactly reestimated flux variables. In the metabolic network studied, some fluxes were found to be either nonidentifiable or nonlinearly correlated. The nonidentifiable fluxes could correctly be predicted a priori using the model identification method applied, whereas the nonlinear flux correlation was revealed only by identification runs using different starting values a posteriori.
Conclusion
This fast, robust and accurate optimization method is useful for highthroughput metabolic flux analysis, a posteriori identification of possible parameter correlations, and also for Monte Carlo simulations to obtain statistical qualities for flux estimates. In this way, it contributes to future quantitative studies of central metabolic networks in the framework of systems biology.
Background
In recent years, metabolic flux analysis (MFA) has become an important tool for quantifying metabolic pathways which is essential for indepth understanding of biological systems. Among the developed tools, ^{13}C flux analysis utilizing ^{13}C labeling patterns of metabolic products that result from feeding ^{13}Clabeled substrates provides detailed information about intracellular pathway fluxes in vivo [1–3]
^{13}Cbased MFA requires carbon flux modeling through the metabolic network, which describes the mathematical relationship between unknown fluxes and the available measurement data set. It requires modeling two connected equation systems, which describe reaction stoichiometry between metabolites and between carbon isotopomers, respectively. Using the model, fluxes can be computed from the measurements by solving a nonlinear leastsquares problem (NLSP). The stoichiometric network is a linear equation system of material balances given for the reactions between metabolites, i.e., S·ν = 0 with S denoting the stoichiometric matrix and ν the fluxes. The fluxes consist of nonmeasured intracellular fluxes ν_{u} and measured effluxes ν_{m}, i.e., ν = (ν_{u}, ν_{m})^{ T }. Typically, realistic models are underdetermined, that is, the rank of S is smaller than the number of entries in ν_{u} [1]. This difference between rank(S) and the number of entries of ν_{u} equals the number of fluxes that have to be chosen as the design parameters Θ (independent flux variables) and are required when parametrizing the network such that ν = F_{flux}(Θ).
Here f(Θ) denotes the objective function to be minimized with respect to Θ and $\mathbb{F}$ (Θ) signifies the model function corresponding to the measured data set η = (η_{1}, η_{2},..., η_{ n })^{ T }consisting of measured ^{13}C labeling data (x_{m}). Such data are typically received as mass isotopomer distributions measured by mass spectrometry (MS) and/or fractional carbon labeling measured by nuclear magnetic resonance (NMR) techniques, and effluxes (ν_{m}). The measurement error ε is typically assumed to have a normal distribution such that ε ∈ N(0, Σ_{ η }), where Σ_{ η }is the covariance matrix of measurements.
For ^{13}Cbased MFA, the applied algorithms for numerical flux estimation are mainly gradientbased local optimization [6, 7] or gradientfree global optimization [8–11] such as simulated annealing (SA) or genetic algorithms (GAs). Also, a hybrid technique of globallocal optimization has been applied [12]. Such algorithms are described in detail elsewhere [13–15]. The stochastic global optimization methods can be inefficient due to the time required to obtain the socalled asymptotic convergence or reachability in high dimensional parameter spaces [11, 16–21]. Moreover, such algorithms of random nature may fail to find the global solution unless the number of samplings tends to infinity, which is practically impossible [21, 22]. In comparison, the gradientbased local optimizations have a much higher convergence speed, but the solution quality depends heavily on starting points [14, 19]. Further to this, reaching the global optimum is ensured only for convex problems, whereas it is nontrivial to determine the convexity of general nonlinear problems with nonlinear constraints [23]. Thus, one may obtain solutions that are not necessarily global [14] and that might vary depending on starting points.
In this regard, a robust method is highly desirable which guarantees a convergence with speed and accuracy for the nonlinear flux estimation problem. In the present work, we developed a method for efficient parametrization of metabolic network by compactification. Additionally, a mathematical method is suggested for model identification that solves a priori flux identifiability problems regarding ^{13}C labeling experiments in terms of model linearization. On this basis, an optimization algorithm was developed that hybridizes two gradientbased optimization tools. The developed approaches were evaluated using the central metabolism of Bacillus subtilis with succinate and glutamate feeding as the only carbon sources. We examined whether global flux solutions are obtainable using the gradientbased deterministic algorithm.
Results
In this section, we describe the mathematical and computational procedures developed and their application to the nonlinear problem of metabolic flux estimation in a realistic metabolic network of Bacillus subtilis. The model represents an inherent difficulty arising from only succinate and glutamate feeding as carbon sources. Succinate is a symmetric molecule having only two distinguishable carbon atoms. On the other hand, glutamate is also converted into succinate. Due to this, the information that is introduced by the ^{13}Clabels of the substrates is severely limited.
Parametrization of stoichiometric network
To formulate a nonlinear least squares problem (NLSP), the system of interest has to be parametrized. The stoichiometric network, which is a rankdeficient linear system (S·ν = 0), was parametrized using the following three steps.
Step (i)
The stoichiometric matrix S is transformed into its reduced row echelon form S_{RRE} using GaussJordan elimination with partial pivoting. Subsequently, each row and column of S_{RRE} is analyzed to identify the dependent and independent variables. The first nonzero element in each row of S_{RRE} that is always 1 is called 'leading 1'. If the i^{th} column of S_{RRE} contains only zeros and one leading 1, then the i^{th} element of ν is dependent, otherwise it is independent. The number of fluxes selected as independent equals the degrees of freedom in the stoichiometric network, that is, the difference between the number of variables and the rank of S. Once this step has been completed, we will see that certain intracellular fluxes and all effluxes (ν_{m}) are chosen as independent.
Step (ii)
Physiologically, metabolic fluxes are constrained such that 0 ≤ ν < ∞, if reversible reactions are considered as two independent reactions. For the effluxes, the corresponding measurements are available. Hence, the effluxes can be bounded, e.g., mean value ± $\sqrt{{\chi}_{\text{1,}\phi}^{\text{2}}}$ × standard deviation, where ${\chi}_{\text{1,}\phi}^{\text{2}}$ denotes the inverse of χ^{2}cumulative distribution function at a certain confidence level of φ, and the resulting range covers 100 × φ % of possible experimental observations.
from (2). If α > (1  ϕ_{ i })^{2} holds and, thus, dν_{ i }/dϕ_{ i }> 1, the sensitivity given by ∂x_{m}/∂ϕ = (∂x_{m}/∂ν)·(∂ν/∂ϕ) increases. Particularly, a higher sensitivity can always be obtained by setting α ≥ 1 due to finite values of fluxes, that is, 0 ≤ ν_{ i }< ∞ or 0 ≤ ϕ_{ i }< 1. Hence, setting the parameter scaling constant α ≥ 1 is more preferable for numerical optimization than α > 0. Moreover, the mapping such as (2) has proven to yield an extremely low curvature of ^{13}C labeling in the parameter space and is advantageous for model linearization [24]. The constant α can be adjusted during the optimization. This is described in detail in the Appendix.
Step (iii)
By symbolically solving the equation system consisting of stoichiometric balances and flux constraints such as (2) for the dependent fluxes (ν_{depend}), we get explicit expressions for all dependent fluxes such that ν_{depend} = F_{flux}(Θ) with Θ = (ϕ_{1}, ϕ_{2},..., ν_{m1}, ν_{m2},...)^{ T }.
Flux identifiability by model linearization
One important question to be answered prior to the numerical flux computation is the a priori parameter identifiability. To this end, the socalled Buchberger's algorithm [25] can be applied to compute the Gröbner basis [26]. However, this polynomial algebraic method can be very time consuming and, thus, is restricted to small systems, corresponding to the fact that the Gröbner basis can be extremely large for a metabolic carbon labeling system [27, 28].
In previous works [28, 29], the Jacobian matrix of the measurement model has been utilized for model identification: the matrix determinant or rank were applied to the flux identifiability analysis in carbon labeling systems. However, these approaches do not tell how to sort out identifiable flux variables in case not all fluxes are identifiable.
Here, null(J($\widehat{\Theta}$)) is the null space of J($\widehat{\Theta}$) and β is a vector containing arbitrary nonzero values. The i^{th} element of the search step $\Delta \widehat{\Theta}$ whose corresponding row of null(J($\widehat{\Theta}$)) consists of zeros can be determined uniquely from the measurement η. Thus, the unique solution of the corresponding i^{th} design parameter can be computed numerically.
Hybrid optimization algorithm
where f(Θ) denotes the objective function specified in (1) and Θ = (ϕ_{1}, ϕ_{2},...,ϕ_{ n }, ν_{m1}, ν_{m2},...)^{ T }the design parameter vector consisting of the [0, 1)fluxes ϕ_{ i }. The effluxes ν_{m} comprise substrate uptake as well as product and biomass formation with known standard deviation σ_{ m }, and ν_{depend} corresponds to all unknown intracellular fluxes.
Analytical gradient and Hessian
Typically, a gradientbased optimization problem is solved more accurately and with higher efficiency when analytical gradients are provided compared to numerical gradients, which are calculated by finitedifference approximation. By providing the analytical Jacobian matrix consisting of ∂x_{m}/∂ϕ and ∂ν_{m}/∂Θ, the gradient of the objective function formulated in (A.1) of the Appendix can be calculated. The analytical Hessian (A.2) is calculated when STRiN is used. In comparison to this, SQP updates the Hessian using the BroydenFletcherGoldfarbShanno (BFGS) formula based on linearization of the gradient [14].
Adjustment of tolerance and parameter scaling constant
Usually, the tolerances placed on parameters or on the objective function which are utilized as termination criteria for numerical optimization cannot clearly be defined in advance. They depend on the errors associated with usersupplied data as well as output sensitivities. Thus, tolerances are set rather empirically. Here, we suggest a method to adjust tolerance values and find a proper value by repeatedly restarting optimization trials. The tolerance adjustment by the restart is implemented in the developed algorithm as shown in Figure 1E. Due to this, the complete process consists of a series of suboptimization trials. At each k^{th} trials, tolerance values are adjusted.
Further to this, the parameter scaling constant α is also adjusted during the optimization. As mentioned in the previous section, α ≥ 1 is the preferable choice due to the finite values of metabolic fluxes. In practice, α is regulated to render the model's Jacobian matrix into better condition. This is obviously advantageous for the gradientbased method, which typically involves matrix inversion to compute the search direction as shown by (6). The exact procedure how to adjust tolerance and the parameter scaling constant is described in the Appendix.
Hybridization
SQP represents the stateoftheart for solving constrained nonlinear optimization problems. By introducing the Lagrangian function, SQP solves a series of quadratic programming subproblems, each involving the minimization of a quadratic approximation of the objective function subject to a linear approximation of the constraints. It shows its strength when solving problems with significant nonlinearity and from remote starting points [14, 33].
STRiN solves a nonlinear large problem by linear approximation and the method of preconditioned conjugate gradients. This algorithm is based on the trustregion algorithm, is computationally inexpensive, and provides rapid convergence in the vicinity of the solution. A downside of STRiN is that it accepts any direction of negative curvature, even when this direction gives an insignificant reduction in the objective function [14]
The key concept of the developed algorithm is hybridizing the merits of these two optimization algorithms, i.e., the loose startingpointdependency of SQP and the convergence speed of STRiN in the solution vicinity. During the optimization trails with the tolerance adjustment, we get a feasible point that may be a more suitable initial guess for the next trial. Accordingly, we initiate a few trials under "relaxed" tolerance conditions using a robust algorithm that is less dependent on the quality of the initial points. Subsequently, the local minima obtained can be tested by another algorithm whether it gives a significant improvement with a rapid convergence.
In this context, the hybridization is carried out as shown in the flow chart in Figure 1. The algorithm is forced to initiate within the feasible region by generating an arbitrary starting point Θ_{0} subject to the inequality constraint ν_{depend}(Θ_{0}) ≥ 0 (Figure 1A). This can be done by generating random numbers of design parameters within their bounds given by (8). For cases in which the random generation is expensive, e.g., when a design parameter has a very narrow range that satisfies the constraints, a deterministic method might be employed to get a starting point in the physically possible region. This is discussed in the Appendix.
As soon as a feasible set of initial values is obtained, the optimization starts using the SQP algorithm under relaxed tolerance criteria (Figure 1B). After a few successful trials, an upper limit (χ^{2}_{UL}) that schedules the initiation of STRiN is placed, e.g., half of the current objective function value, χ^{2}_{UL} = 1/2 f(Θ*_{ k }), where Θ*_{ k }is the local minimizer isolated from the k^{th} suboptimization trial. Afterwards, if the objective value f(Θ*_{ k }) resulting from the current SQP trial (Figure 1C) is smaller than χ^{2}_{UL}, the STRiN optimization is activated (Figure 1D); otherwise, the SQP is repeated. During the STRiN optimization, the progress of trials is monitored to prevent the nonreducing problem associated with the STRiN algorithm mentioned above. As a criterion of the progress check, the relative improvement of f(Θ*_{ k })  f(Θ*_{k1})/f(Θ*_{ k }) can be measured at each STRiN termination. If the improvement is insignificant, i.e., less than a userspecified value (e.g., 0.05), the STRiN receives penalty. In case a few successive trials show insignificant improvement, the algorithm sets a new upperlimit χ^{2}_{UL} and returns to the SQP optimization. For the restart, the tolerance of the previous SQP trial is reutilized. This prevents too large changes in tolerance, which can be caused by the STRiN trials with insignificant improvement.
As a termination criterion of the hybrid optimization, the changes in local optima are measured (Figure 1C), e.g., the absolute value of the slope resulting from the linear regression of y = (f(Θ*_{k4}), f(Θ*_{k3}),...f(Θ*_{ k }))^{ T }with respect to x = (1, 2,..., 5)^{ T }. If the absolute slope is less than a userspecified small value for a few further trials, the optimization can be terminated ultimately.
Application to Bacillus subtilis metabolic network
Parametrization
The metabolic network in Figure 2 has 27 intracellular fluxes (10 bidirectional and 17 unidirectional fluxes), 5 effluxes, and 10 anabolic fluxes expressed in terms of Y_{ XS }. For 15 intracellular metabolites defined in the network, 15 flux balances were set up that are linearly independent of each other (Appendix). The stoichiometric matrix S was obtained by symbolic differentiation of the balances with respect to the whole fluxes including Y_{XS}. When considering that the flux of glutamate uptake is unity due to the normalization, the network owns 17 degrees of freedom. Accordingly, 12 intracellular fluxes ν_{2r}, ν_{3r}, ν_{4r}, ν_{5r}, ν_{6r}, ν_{7r}, ν_{8r}, ν_{14}, ν_{14r}, ν_{15r}, ν_{16}, ν_{17r}, 4 effluxes, and Y_{XS} were recognized as independent variables from the reduced row echelon form S_{RRE} of S. These independent intracellular fluxes were transformed into [0, 1)fluxes as given by (2). Note that the counterpart of a backward flux of a bidirectional reaction can also be selected as an independent variable. This just depends on how to arrange the entries in the flux vector. For instance, when ν_{2r}is followed by ν_{2} in the flux vector, ν_{2} will be recognized as an independent variable instead of ν_{2r}.
Solving the equation system of the flux balances (Appendix) and the [0, 1)flux equations for the dependent fluxes gives their explicit analytical expressions with respect to design parameters, i.e., ν_{ i }= n_{flux}(Θ) with 19 design parameters of Θ = (ϕ_{2r}, ϕ_{3r}, ϕ_{4r}, ϕ_{5r}, ϕ_{6r}, ϕ_{7r}, ϕ_{8r}, ϕ_{14}, ϕ_{14r}, ϕ_{15r}, ϕ_{16}, ϕ_{17r}, succinate_{ex}, acetate_{ex}, αketoglutarate_{ex}, fumarate_{ex}, Y_{XS}, p CO_{2,1}, p CO_{2,2})^{ T }. The parameters p CO_{2,1} and p CO_{2,2} estimate the CO_{2} labeling pattern for each ^{13}C labeling experiment: CO_{2} has only one carbon and is incorporated into other metabolites by carboxylation. Thus, its labeling state can simply be determined from other ^{13}C labeling data as an additional parameter of NLSP in case it is not measured.
For flux reestimation studies, the default values of the design parameter were Θ_{default} = (0, 0.6341, 0.4337, 0.3148, 0.5162, 0.3436, 0.5884, 0.7066, 0.0643, 0.6835, 0.7125, 0.6756, 1.164, 0.021, 0.006, 0.02, 0.0859, 0.5601, 0.5601) at α = 1, which results in the flux values given in Figure 2.
Model identification
For the identifiability studies, the main question was whether Θ has one unique solution set that can be estimated from the available mass isotopomer measurements in terms of the gradientbased optimization. The effluxes (succinate_{ex}, acetate_{ex}, αketoglutarate_{ex}, fumarate_{ex}) and Y_{XS} are excluded for the identifiability analysis because they can be estimated as long as the corresponding measurements are available. To examine various flux scenarios, 200 stochastic simulations were performed for each of the 630 possible ^{13}Cexperimental designs. To this end, uniform random numbers between 0 and 0.95 were created for Θ. At each simulation the null space of J(Θ) was calculated.
According to the stochastic simulations, none of the designs had an empty null space and the entry of null(J(Θ)) corresponding to ϕ_{2r}was always nonzero disregarding the flux state and designs. Most designs, including the design considered here, yielded an identical null space of:
null(J(Θ)) = (1 0 0 0 0 0 0 0 0 0 0 0 0 0)^{ T }.
Accordingly, the design parameter ϕ_{2r}is not expected to have a unique solution, i.e., it has infinitively many solutions. Except ν_{2} and ν_{2r}, other fluxes were not functions of ϕ_{2r}. This means that the bidirectional glucose 6P isomerase cannot be determined from any ^{13}C substrates or dual substrate combinations. Only the net flux of the reaction can be calculated from the stoichiometry. Therefore, one can either regard the glucose 6P isomerase as a unidirectional reaction or set ϕ_{2} as an arbitrary constant. This renders the carbon flux network to have an empty null space of J(Θ), i.e., full rank for J(Θ), and all other design parameters are theoretically expected to have unique solutions. This theoretical expectation is examined later using the developed algorithm.
Hybrid optimization with tolerance adjustment
To evaluate the efficiency of the developed optimization algorithm, the flux values calculated in advance were reestimated numerically from the carbon mass isotopomer distributions (MDVs) of output metabolites matching the default flux values (ν_{default}). To this end, the default flux values, depicted in Figure 2, were calculated by the flux function F_{flux}(Θ_{default}) obtained by parametrization. Subsequently, the ^{13}C labeling data of output metabolites were calculated from ν_{default} and Θ_{default}. The default fluxes were reestimated by solving the constrained NLSP (8) whose inputs (η) were the effluxes, Y_{XS}, and the MDVs.
The worst algorithm among the local methods was the STRiN using analytical gradient and Hessian (STRiN ∇_{user} H_{user}). It was rapid at the beginning but improved the objective only from 10^{7} to 10^{5}. Afterwards, the optimization crashed with the carbon labeling system becoming completely singular. As previously mentioned, this seems to be due to the STRiN limitation of accepting any direction of negative curvature even if it does not give significant reduction in the objective function. Also, the global optimization methods tested, which were the GA and the SA, were notably timeinefficient compared to the local methods. For 1500 seconds of optimization time, the objective function was very slowly decreased by two orders of magnitude only. SA decreased f(Θ) only to 1 × 10^{3} after 5.1 hours and GA to 2 × 10^{2} after 6.7 hours. Thus, global optimization is expected to demand much greater time to get the equivalent result as HATA.
HATA optimization from different starting points
The flux reestimation was also performed from 200 different starting points created randomly and by simultaneously generating uniform random numbers for the [0, 1)flux ϕ_{2r}= ν_{2r}/(α + ν_{2r}). This was to examine whether the quality and speed of convergence has any startingpointdependency and whether ϕ_{2r}, which a priori expected to have infinite solutions, has any influence on optimization results.
All 200 optimization runs using the hybrid algorithm were terminated successfully. On average, HATA required 5.3 ± 1.6 min to converge to a weighted objective function value smaller than 1 × 10^{10}, which equals the nonweighted sum of squares less than 2 × 10^{18}. As shown in Figure 4A, 99.5 % of the flux reestimations were completed within 10 min and about 76 % in less than 6 min.
By the null space investigation (9) based on (7), the design parameters ϕ_{3r}and ϕ_{4r}were not expected to be correlated. In accordance with this, the fluxes could rarely be refitted in case one or both of these parameters were set arbitrarily. This indicates that the parameters ϕ_{3r}and ϕ_{4r}have a certain significant correlation that must be fulfilled at termination to yield unique solutions for other fluxes. As expected, the flux estimates of TK1 and those of TA were found to be nonlinearly correlated with each other as depicted in Figure 5B for ν_{3r}and ν_{4r}. The true solution, (ϕ_{3r}, ϕ_{4r}) = (0.6341, 0.4337) also lies on the curve. At the termination, it is obvious that the output sensitivities of ∂x_{m}/∂ϕ_{3r}and ∂x_{m}/∂ϕ_{4r}are low: the objective function reached a very small value of nearly zero (f($\widehat{\Theta}$) < 1 × 10^{10}) and we get arbitrary but correlated values for ϕ_{3r}and ϕ_{4r}. Due to this, we cannot practically get ϕ_{3r}and ϕ_{4r}estimated correctly and, accordingly, ν_{3}, ν_{3r}, ν_{4}, and ν_{4r}, while the corresponding net fluxes can always be calculated from other correct flux estimates in terms of stoichiometry.
We observed that a unique estimation of all four fluxes would be possible only if the mass isotopomers of sedoheptulose 7phosphate were additionally measured. It was observed that providing these mass isotopomers renders the output sensitivities of ∂x_{m}/∂ϕ_{3r}and ∂x_{m}/∂ϕ_{4r}much less dependent on starting points. In contrast, the output sensitivities vary strongly and even become zero when the mass isotopomers of sedoheptulose 7phosphate are not involved in the objective.
Further to this, the fluxes were not found to be influenced by the glucose 6phosphate isomerase fluxes represented by ϕ_{2r}. Disregarding the random values of ϕ_{2r}, the fluxes except ν_{3}, ν_{3r}, ν_{4}, and ν_{4r}converged to their true values. Also ϕ_{3r}and ϕ _{4r}did not show any correlation with ϕ_{2r}as depicted in Figure 5C and Figure 5D, respectively. Hence, we can conclude that the observed nonlinear correlation between ϕ_{3r}and ϕ_{4r}becomes visible solely from using the different starting points for the optimization. In addition, including ϕ_{2r}as an additional design parameter gives the same result for the fluxes, whereas ϕ_{2r}results in randomly distributed values without any obvious pattern. This shows that the flux identification carried out by the model linearization and null space investigation is appropriate for its practical use.
Discussion
The tradeoff between robustness and convergence speed is a central issue in numerical optimization [14]. The developed method fulfills both criteria and enables exact metabolic flux estimation. This was accomplished on the basis of parametrizing the metabolic flux network by compactification and hybridizing the merits of two different gradientbased algorithms. Interestingly, HATA, which hybridizes the worst local method with SQP ∇_{user} and utilizes the tolerance adjustment strategy, has proven to be superior to other approaches. This may be by virtue of providing STRiN with a more appropriate starting point advanced by SQP. Thus, because STRiN is computationally inexpensive and rapid in the vicinity of solution, it achieves a fast convergence. Furthermore, the interactive switching between the algorithms prevents STRiN from accepting an unfavorable search direction without significant objective improvement.
Further to this, we observed that the parametrization by the compactification of independent [0, ∞)fluxes into [0, 1)variables is advantageous compared to the noncompactified case. When the optimization was performed by selecting the [0, ∞)fluxes (ν_{2r}, ν_{3r}, ν_{4r}, ν_{5r}, ν_{6r}, ν_{7r}, ν_{8r}, ν_{14}, ν_{14r}, ν_{15r}, ν_{16}, ν_{17r}with upper bounds of 200) directly as design parameters, the convergence took drastically longer (27 ± 45 min) than the cases using the compactified fluxes. Moreover, optimization runs often failed to converge or terminated at suboptimal points. In comparison to this, introducing [0, 1)fluxes improved both the robustness and speed of convergence as demonstrated above. The same advantage is probably achieved when using the [0, 1)compactifications such as exchange fluxes or flux partitioning ratios applied elsewhere [9, 24, 35]. Our compactification approach allows a straight parametrization because independent fluxes can easily be recognized from the stoichiometric matrix and compactified using a single rule of (2). Moreover, it is straightforwardly differentiable when considering that compactified [0, 1)fluxes are continuous and smooth in the [0, ∞)flux space and vice versa.
The parameter scaling constant α introduced during the compactification provides a way to implement numerical flux estimation with efficiency. In particular, the parameter scaling constant α is updated during optimization to render the model's Jacobian matrix into better condition. This is advantageous for the gradientbased methods for computing the search direction because it typically involves the inverse of the Jacobian matrix as given by (6). We examined whether different α values affect the efficiency of optimization. The α value was observed to affect the convergence speed but not the accuracy of the optimization. When starting from the same initial guess applied to the optimization trial in Figure 4, the fastest convergence speed was obtained when α is fixed at 10 (about 200 s for f(Θ) < 10^{14}). The speed efficiency decreased if α is smaller or larger than 10. A trial with α adjustment initiated at α = 1 automatically updated α to 10 after a few suboptimization trials and took about 300 s to achieve f(Θ) < 10^{14}. Although it was slower than when α was fixed at 10, the proposed strategy with α adjustment recognized the optimal value. Among the 200 HATA optimization runs with α adjusting between 1 and 10 initialized with α = 1, 112 cases updated α to 10 after a few to several suboptimization trials while α remained 1 in 78 cases. Moreover, the time taken for the optimization runs using different starting points was not necessarily shorter when α was updated to 10. This indicates that the most adequate α value for the optimization may depend on starting conditions.
One interest of numerical optimization is whether the global solution can be obtained by the algorithm applied. We observed that the local minimizers from 200 random simulations agreed with each other and also with the true flux values, the global solution. Depending on experimental noise or metabolic state, the problem may contain a few to several local saddle points. In this case, one can consider a further hybridization with a global optimization approach. Global optimization algorithms are inefficient in terms of convergence speed. On the other hand, gradientbased algorithms can converge rapidly but lack a global perspective for nonconvex problems. Also in other cases, the combination of global and local search procedures has proven to offer the advantages of both methods while offsetting their disadvantages [17, 37–39]. Such hybridization with a global optimization algorithm is useful for nonconvex problems containing several local minima and can be done in many ways [17, 39, 40]. For instance, the local minima obtained from a gradientbased method can be used to update the individuals in the population to prepare the next generation when combining with the GA [17, 40, 39, 41]. Because we could prove the speed and accuracy of the developed optimization method, we consider it very useful for investigating metabolic fluxes and for developing efficient algorithms combined with a global optimization method.
In practice, the quality of the flux fitting can be judged by the χ^{2}test [42, 43] when using noisy measurement data. In our case, the upper critical χ^{2} limit is 206 at 95 % confidence level at the given degree of freedom of 174 (188 ^{13}C labeling data + 5 effluxes  19 design parameters). Hence, if the weighted objective value (1) is smaller than the critical value of 103 (this is due to the factor 1/2 in the objective function), the fit can be regarded as acceptable. This value can also be used as a termination criterion for the optimization control in Figure 1E. As shown in Figure 4A, the value was achieved at 150 s by HATA, at 170 s by SQP with ∇_{user}, and at 750 s of optimization time by SQP with ∇_{finite}. In comparison to this, GA or SA was expected to achieve the acceptable f(Θ) after 78 hours or more (by extrapolation of log_{10} f(Θ)with respect to optimization time). Hence, an efficient Monte Carlo simulation cannot be carried out by the global optimization approaches.
Further, we compared HATA with SQP with ∇_{user} by conducting 100 stochastic simulations under identical conditions. At each simulation run, normal random numbers were generated for the experimental data set (^{13}C mass isotopomer fractions and extracellular fluxes) by assuming a unit standard deviation of 1 × 10^{3} and uniform random numbers for the starting points of design parameters. Subsequently, parameter reestimation was performed in terms of HATA and SQP with ∇_{user}, respectively. Each optimization run was stopped if f(Θ) becomes smaller than the critical χ^{2} limit of 103 or does not give any further improvement. On average, f(Θ) at termination was 102 ± 26 for HATA and 167 ± 567 for SQP with ∇_{user}, and the time taken for the optimization was 353 ± 216 s for HATA and 1234 ± 1356 s for SQP with ∇ _{user}. The results support that the developed hybrid algorithm using tolerance adjustment will be more efficient and robust when computing metabolic fluxes from noisy experimental data in a realcase application. In contrast, the efficiency of the SQP algorithm seems to be affected by initial conditions (starting points and experimental data set) when considering the large variation in f(Θ) as well as in the computation time. In this regard, HATA will also achieve the most efficient optimization in the practical sense whenever repeated optimization runs are necessary. A further reduction in computation time can be obtained by utilizing a minimal set of isotopomer balances described by Antoniewicz et al. (2007) [3].
Conclusion
We developed and examined a developed hybrid algorithm for numerical ^{13}C flux estimation using a metabolic network model parametrized using a single compactification rule. We could prove its practical usefulness using the central metabolic network of Bacillus subtilis, which is a prokaryotic model organism and an industrially relevant strain. Based on the parametrization by compactification and model identification, the hybrid algorithm with tolerance adjustment allowed a fast and robust flux computation from the ^{13}C labeling data created from ^{13}Clabeled succinate and glutamate feeding. Such fast optimization was also found to be essential for the a posteriori identification of possible parameter correlations and also for the Monte Carlo approach to obtain statistical qualities for metabolic flux estimates. Therefore, this work represents an important contribution to the quantitative study of metabolic networks in the framework of systems biology.
Methods
For the flux reestimation simulations of the central metabolic network of B. subtilis, the following inputs and outputs were applied.
^{13}C Input Substrates
The choice of input tracer substrates was made by means of the Doptimality criterion [44]. Since parallel experiments using different ^{13}C input labels yield more information as theoretically shown for the respirometric flux analysis utilizing CO_{2} labeling measurement [4]. Therefore, experimental designs for two parallel cultivations were also considered here. Among the investigated experimental designs (630 possible designs from commercially available ^{13}C glutamate and ^{13}C succinate species; refer to Additional file 1!), a design combining data from a ^{13}C labeling experiment using nonlabeled glutamate and [2,3^{13}C_{2}] succinate with those from an experiment using [1,2^{13}C_{2}] glutamate and [1,4^{13}C_{2}] succinate was expected to provide the richest information at a reasonable expense. This design was selected to study flux identifiability and to examine the developed hybrid optimization.
^{13}C Output Metabolites
For simulation studies, we assumed that the mass isotopomers of proteinogenic amino acids such as alanine, valine, serine, threonine, glycine, aspartate, isoleucine, leucine, phenylalanine, tyrosine, arginine, histidine, and glutamate are measurable. In addition, the mass isotopomers of ribose 5phosphate from RNA hydrolysate and hexose from carbohydrate hydrolysate are also assumed to be measurable. This totally gives a mass isotopomer data size of 188 for the parallel experiment. The measurements of mass isotopomer distributions were assumed to have a unit standard deviation of 1 × 10^{4}.
Software implementation
All simulations were carried out on a PC equipped with a CPU of 2.4 GHz. The hybrid optimization algorithm depicted in Figure 1 was implemented in MATLAB (The Mathworks Inc., Natick, MA). Numerical optimization was carried out using the Optimization Toolbox (Version 3.0.4) and the Genetic Algorithm and Direct Search Toolbox (Version 2.1) of MATLAB. The symbolic operations required for parametrization and for computing partial derivatives were conducted using the Symbolic Math Toolbox (Version 3.1.4) of MATLAB. During simulations, random numbers were generated using the Statistics Toolbox (Version 5.2) of MATLAB. The random numbers for [0, 1)fluxes (ϕ_{ i }) were bounded such that the relative values of the corresponding dependent fluxes (ν_{depend}) lie between 0 and 20.
An example of the hybrid optimization for the problem described in the current work was coded using MATLAB (Hybridoptimizer.m) and supplemented (see Additional file 2). The method can be incorporated with any isotopomeric models [2, 3, 5, 8] that are analytically differentiable.
Appendix
Gradient and Hessian
In practice, the constrained problem (1) is formulated as the Lagrangian function, a linear combination of the objective function and the constraints [14].
Adjustment of tolerance and parameter scaling constant
The tolerance adjustment by the restart in Figure 1E requires a series of suboptimization trials. The first optimization trial starts with "sufficiently relaxed" tolerances both placed on parameters and objective function, e.g., of 100. When the k^{th} trial is acceptable (e.g., first order optimality conditions are satisfied to the specified tolerance or changes in parameters or function are smaller than the specified tolerance), the (k + 1)^{th} trial restarts with a 10fold decreased tolerance, where the k^{th} local minimizer Θ*_{ k }is used for the new starting point. The previous trial provides a more adequate starting point for the next, and it can be advantageous for local search methods. This is because the efficiency of local methods depends heavily on the quality of starting points [14, 19]. In case the current trial fails (e.g., by exceeding the maximum number of iterations allowed in the estimation process), the tolerance is increased and the next trial starts from the iterate Θ° recorded for the smallest function value up to the current trial in the feasible region (Figure 1E). The tolerance at which an optimization trial fails is registered. In case optimization trials successively fail, which is rarely the case, the step of tolerance increase is set smaller such as 5, 10, 50, or 100fold, etc., of the registered value. Correspondingly, the tolerance decrease is also made by this reduced step until the current tolerance matches the registered. This allows finding an appropriate tolerance value for NLSP.
Deterministic Method of Feasible Starting Point Generation
Here, N_{neg.flux} signifies the number of negativevalued fluxes, b_{flux} the steadystate stoichiometric flux balances set up for metabolites (see the next section!) evaluated at the current iterate Θ_{ k }, and n the number of the balances. By the first term in the objective and the first constraint, the minimization is directed such that all fluxes become nonnegative. The second term is to prevent flux values calculated that may not obey stoichiometry due to the limited floatingpoint accuracy (roundoff error). This approach is useful for preparing initial values while preserving random nature for Monte Carlo simulations. In practice, the upper bound of ϕ_{ i }can be set to a value smaller than 1 that gives sufficiently large flux values so that ϕ_{ i }has a closed interval.
Stationary State Stoichiometric Balances are set up around 15 metabolites specified in Figure 2 as follows. The values multiplied with Y_{XS} equals the strain specific precursor demand for growth (mmol precursor per gram biomass).
glucose 6P: ν_{2}  (ν_{1} + ν_{2r}+ 0.4217Y_{XS}) = 0
fructose 6P: (ν_{2r}+ ν_{4} + ν_{5} + ν_{6})  (ν_{2} + ν_{4r}+ ν_{5r}+ ν_{6r}) = 0
glyceraldehyde 3P: (ν_{3} + ν_{4r}+ ν_{5} + 2ν_{6r}+ ν_{7})  (ν_{3r}+ ν_{4} + ν_{5r}+ 2ν_{6} + ν_{7r}+0.4659Y_{XS}) = 0
3phosphoglycerate: (ν_{7r}+ ν_{8})  (ν_{7} + ν_{8r}+ 0.2209Y_{XS}) = 0
phosphoenolpyruvate: (ν_{16} + ν_{8r})  (ν_{8} + ν_{9} + 1.7854Y_{XS}) = 0
pyruvate: (ν_{9} + ν_{15})  (ν_{10} + ν_{15r}+ 0.9964Y_{XS}) = 0
pentose 5P: (ν_{1} + 2ν_{3r}+ ν_{5r})  (2ν_{3} + ν_{5} + 1.1622Y_{XS}) = 0
erythrose 4P: (ν_{4} + ν_{5r})  (ν_{4r}+ ν_{5} + 0.4659Y_{XS}) = 0
sedoheptulose 7P: (ν_{3} + ν_{4r})  (ν_{3r}+ ν_{4}) = 0
acetylCoA: ν_{10}  (ν_{11} + acetate_{ex} + 3.1585Y_{XS}) = 0
isocitrate: ν_{11}  ν_{12} = 0
αketoglutarate: (ν_{12} + ν_{17})  (ν_{13} + ν_{17r}+ αketoglutarate_{ex}) = 0
glutamate: (ν_{17r}+ glutamate_{ex})  (ν_{17} + 1.4769Y_{XS}) = 0
succinate: (ν_{13} + ν_{14r}+ succinate_{ex})  (ν_{14} + fumarate_{ex}) = 0
oxaloacetate: (ν_{14} + ν_{15r})  (ν_{11} + ν_{14r}+ ν_{15} + ν_{16} + 1.3131Y_{XS}) = 0
List of abbreviations used
 BFGS :

BroydenFletcherGoldfarbShanno
 GA :

genetic algorithm
 HATA :

hybrid algorithm with tolerance adjustment
 KEGG :

Kyoto Encyclopedia of Genes and Genomes
 MDV :

mass isotopomer distribution vector
 MFA :

metabolic flux analysis
 NLSP :

nonlinear leastsquares problem
 SA :

simulated annealing
 SQP :

sequential quadratic programming SQP∇_{ user }: SQP optimization with analytical gradient
 SQP∇ _{ finite } :

SQP optimization using the numerical gradient obtained by finite difference
 STRiN :

subspace trustregion method based on the interiorreflective Newton method
 STRiN∇ _{ user } H _{ user } :

STRiN using analytical gradient and Hessian
 TA :

transaldolase
 TK :

transketolase
Declarations
Acknowledgements
O. Frick thanks DFG (Deutsche Forschungsgemeinschaft), project He 3092/6 for the support of his Ph.D. We wish to thank Andrew Marsh at the University of Louisville's James Graham Brown Cancer Center for editing.
Authors’ Affiliations
References
 Wiechert W: 13C metabolic flux analysis. Metab Eng. 2001, 3 (3): 195206.View ArticlePubMedGoogle Scholar
 Yang TH, Wittmann C, Heinzle E: Metabolic network simulation using logical loop algorithm and Jacobian matrix. Metab Eng. 2004, 6 (4): 256267.View ArticlePubMedGoogle Scholar
 Antoniewicz MR, Kelleher JK, Stephanopoulos G: Elementary metabolite units (EMU): a novel framework for modeling isotopic distributions. Metab Eng. 2007, 9 (1): 6886.PubMed CentralView ArticlePubMedGoogle Scholar
 Yang TH, Heinzle E, Wittmann C: Theoretical aspects of 13C metabolic flux analysis with sole quantification of carbon dioxide labeling. Comput Biol Chem. 2005, 29 (2): 121133.View ArticlePubMedGoogle Scholar
 Wiechert W, Mollney M, Isermann N, Wurzel M, de Graaf AA: Bidirectional reaction steps in metabolic networks: III. Explicit solution and analysis of isotopomer labeling systems. Biotechnol Bioeng. 1999, 66 (2): 6985.View ArticlePubMedGoogle Scholar
 Wiechert W, Siefke C, de Graaf A, Marx A: Bidirectional reaction steps in metabolic networks: II. Flux estimation and statistical analysis. Biotechnol Bioeng. 1997, 55: 118135.View ArticlePubMedGoogle Scholar
 Wittmann C, Heinzle E: Genealogy profiling through strain improvement by using metabolic network analysis: metabolic flux genealogy of several generations of lysineproducing Corynebacteria. Appl Environ Microbiol. 2002, 68 (12): 58435859.PubMed CentralView ArticlePubMedGoogle Scholar
 ArauzoBravo MJ, Shimizu K: An improved method for statistical analysis of metabolic flux analysis using isotopomer mapping matrices with analytical expressions. J Biotechnol. 2003, 105 (12): 117133.View ArticlePubMedGoogle Scholar
 Dauner M, Bailey JE, Sauer U: Metabolic flux analysis with a comprehensive isotopomer model in Bacillus subtilis. Biotechnol Bioeng. 2001, 76 (2): 144156.View ArticlePubMedGoogle Scholar
 Forbes NS, Clark DS, Blanch HW: Using isotopomer path tracing to quantify metabolic fluxes in pathway models containing reversible reactions. Biotechnol Bioeng. 2001, 74 (3): 196211.View ArticlePubMedGoogle Scholar
 Schmidt K, Nielsen J, Villadsen J: Quantitative analysis of metabolic fluxes in Escherichia coli, using twodimensional NMR spectroscopy and complete isotopomer models. J Biotechnol. 1999, 71 (13): 175189.View ArticlePubMedGoogle Scholar
 Zhao J, Shimizu K: Metabolic flux analysis of Escherichia coli K12 grown on 13Clabeled acetate and glucose using GCMS and powerful flux calculation method. J Biotechnol. 2003, 101 (2): 101117.View ArticlePubMedGoogle Scholar
 Floudas CA, Pardalos PM: Recent advances in global optimization. Princeton series in computer science. 1992, x, 633 p.Princeton, N.J. , Princeton University PressGoogle Scholar
 Nocedal J, Wright SJ: Numerical optimization. Springer series in operations research. 1999, xx, 636 p.New York , SpringerGoogle Scholar
 Press WH: Numerical recipes in C: the art of scientific computing. 1992, xxvi, 994 p.Cambridge ; New York , Cambridge University Press, 2ndGoogle Scholar
 Brackin P, Colton SC: Using genetic algorithms to set target values for engineering characteristics in the house of quality. J Comput Inf Sci Eng. 2002, 2 (2): 106114. ASMEView ArticleGoogle Scholar
 Kelner V, Capitanescu F, Léonard O, Wehenkel L: An hybrid optimization technique coupling an evolutionary and a local search algorithm . J Comput Appl Math. 2007Google Scholar
 Lambert TW, Hittle DC: Optimization of autonomous village electrification systems by simulated annealing. Solar Energy. 2000, 68 (1): 121132.View ArticleGoogle Scholar
 Mendes P, Kell D: Nonlinear optimization of biochemical pathways: applications to metabolic engineering and parameter estimation. Bioinformatics (Oxford, England). 1998, 14 (10): 869883.View ArticleGoogle Scholar
 Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome research. 2003, 13 (11): 24672474.PubMed CentralView ArticlePubMedGoogle Scholar
 Xu P: A hybrid global optimization method: The multidimensional case. J Comput Appl Math. 2003, 155 (2): 423446.View ArticleGoogle Scholar
 Long CE, Polisetty PK, Gatzke EP: Nonlinear model predictive control using deterministic global optimization. Journal of Process Control. 2006, 16 (6): 635643.View ArticleGoogle Scholar
 Nash SG, Sofer A: Linear and Nonlinear Programming. 1996, New York , McGrawHillGoogle Scholar
 Wiechert W, de Graaf A: Bidirectional reaction steps in metabolic networks: I. Modeling and simulation of carbon isotope labeling experiments. Biotechnol Bioeng. 1997, 55: 102117.Google Scholar
 Buchberger B: An algorithmical criterion for the solvability of algebraic system of equation. Aequationes Mathematicae. 1988, 4: 4550.Google Scholar
 Saccomani MP: Some results on parameter identification of nonlinear systems. Cardiovascular Engineering: An International Journal. 2004, 4: 95102.View ArticleGoogle Scholar
 Wiechert W: Algebraic methods for the analysis of redundancy and identifiability in metabolic 13C labelling systems. Bioinformatics: From nucleic acids and proteins to cell metabolism. Edited by: Lessel U. 1995, 169184. Weinheim , Verlag ChemieGoogle Scholar
 van Winden WA, Heijnen JJ, Verheijen PJ, Grievink J: A priori analysis of metabolic flux identifiability from 13Clabeling data. Biotechnol Bioeng. 2001, 74 (6): 505516.View ArticlePubMedGoogle Scholar
 Isermann N, Wiechert W: Metabolic isotopomer labeling systems. Part II: structural flux identifiability analysis. Math Biosci. 2003, 183 (2): 175214.View ArticlePubMedGoogle Scholar
 Klamt S, Schuster S, Gilles ED: Calculability analysis in underdetermined metabolic networks illustrated by a model of the central metabolism in purple nonsulfur bacteria. Biotechnol Bioeng. 2002, 77 (7): 734751.View ArticlePubMedGoogle Scholar
 Boggs PT, Tolle JW: Sequential quadratic programming. Acta Numerica. 1995, 4: 151.View ArticleGoogle Scholar
 Coleman TF, Li Y: On the convergence of reflective Newton methods for largescale nonlinear minimization subject to bounds. 1992, 36 p.Ithaca, NY , Cornell Theory Center, Cornell UniversityGoogle Scholar
 Schittowski K: NLQPL: A FORTRANsubroutine solving constrained nonlinear programming problems. Annals of Operations Research. 1985, 5: 485500.View ArticleGoogle Scholar
 Sauer U, Hatzimanikatis V, Hohmann HP, Manneberg M, van Loon AP, Bailey JE: Physiology and metabolic fluxes of wildtype and riboflavinproducing Bacillus subtilis. Appl Environ Microbiol. 1996, 62 (10): 36873696.PubMed CentralPubMedGoogle Scholar
 Wittmann C, Heinzle E: Application of MALDITOF MS to lysineproducing Corynebacterium glutamicum: a novel approach for metabolic flux analysis. Eur J Biochem. 2001, 268 (8): 24412455.View ArticlePubMedGoogle Scholar
 Hill MC, Osterby O: Determining extreme parameter correlation in ground water models. Ground water. 2003, 41 (4): 420430.View ArticlePubMedGoogle Scholar
 Balasubramanian P, Bettina SJ, Pushpavanam S, Balaraman KS: Kinetic parameter estimation in hydrocracking using a combination of genetic algorithm and sequential quadratic programming. Ind Eng Chem Res. 2003, 42 (20): 47234731.View ArticleGoogle Scholar
 Klepeis JL, Pieja MJ, Floudas CA: Hybrid global optimization algorithms for protein structure prediction: alternating hybrids. Biophysical journal. 2003, 84 (2 Pt 1): 869882.PubMed CentralView ArticlePubMedGoogle Scholar
 Xu YG, Li GR, Wu ZP: A novel hybrid genetic algorithm using local optimizer based on heuristic pattern move. Applied Artificial Intelligence. 2001, 15: 601631.View ArticleGoogle Scholar
 Tahk MJ, Woo HW, Park MS: A hybrid optimization method of evolutionary and gradient search. Engineering Optimization. 2007, 39 (1): 87104.View ArticleGoogle Scholar
 Mahinthakumar GK, Mohamed S: Hybrid genetic algorithm  Local seach methods for solving groundwater source identification inverse problems. J Water Resour Plng and Mgmt. 2005, 131: 4557.View ArticleGoogle Scholar
 Antoniewicz MR, Kraynie DF, Laffend LA, GonzalezLergier J, Kelleher JK, Stephanopoulos G: Metabolic flux analysis in a nonstationary system: fedbatch fermentation of a high yielding strain of E. coli producing 1, 3propanediol. Metab Eng. 2007, 9 (3): 277292. 2007/04/03PubMed CentralView ArticlePubMedGoogle Scholar
 Massart DL: Handbook of chemometrics and qualimetrics. Data handling in science and technology; v 20. 1997, Amsterdam; New York , ElsevierGoogle Scholar
 Mollney M, Wiechert W, Kownatzki D, de Graaf AA: Bidirectional reaction steps in metabolic networks: IV. Optimal design of isotopomer labeling experiments. Biotechnol Bioeng. 1999, 66 (2): 86103.View ArticlePubMedGoogle 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.