Hybrid optimization for 13C metabolic flux analysis using systems parametrized by compactification
© Yang et al; licensee BioMed Central Ltd. 2008
Received: 08 October 2007
Accepted: 26 March 2008
Published: 26 March 2008
The importance and power of isotope-based 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. 13C metabolic flux analysis aims to compute in vivo metabolic fluxes in terms of metabolite balancing extended by carbon isotopomer balances and involves a nonlinear least-squares problem. To solve the problem more efficiently, improved numerical optimization techniques are necessary.
For flux computation, we developed a gradient-based 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 non-identifiable 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 13C 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 re-estimated flux variables. In the metabolic network studied, some fluxes were found to be either non-identifiable or nonlinearly correlated. The non-identifiable 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.
This fast, robust and accurate optimization method is useful for high-throughput 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.
In recent years, metabolic flux analysis (MFA) has become an important tool for quantifying metabolic pathways which is essential for in-depth understanding of biological systems. Among the developed tools, 13C flux analysis utilizing 13C labeling patterns of metabolic products that result from feeding 13C-labeled substrates provides detailed information about intracellular pathway fluxes in vivo [1–3]
13C-based 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 least-squares 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 non-measured 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 . 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 ν = Fflux(Θ).
Here f(Θ) denotes the objective function to be minimized with respect to Θ and (Θ) signifies the model function corresponding to the measured data set η = (η1, η2,..., η n ) T consisting of measured 13C labeling data (xm). 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 13C-based MFA, the applied algorithms for numerical flux estimation are mainly gradient-based local optimization [6, 7] or gradient-free global optimization [8–11] such as simulated annealing (SA) or genetic algorithms (GAs). Also, a hybrid technique of global-local optimization has been applied . 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 so-called 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 gradient-based 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 . Thus, one may obtain solutions that are not necessarily global  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 13C labeling experiments in terms of model linearization. On this basis, an optimization algorithm was developed that hybridizes two gradient-based 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 gradient-based deterministic algorithm.
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 13C-labels of the substrates is severely limited.
Parametrization of stoichiometric network
To formulate a non-linear least squares problem (NLSP), the system of interest has to be parametrized. The stoichiometric network, which is a rank-deficient linear system (S·ν = 0), was parametrized using the following three steps.
The stoichiometric matrix S is transformed into its reduced row echelon form SRRE using Gauss-Jordan elimination with partial pivoting. Subsequently, each row and column of SRRE is analyzed to identify the dependent and independent variables. The first non-zero element in each row of SRRE that is always 1 is called 'leading 1'. If the ith column of SRRE contains only zeros and one leading 1, then the ith 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.
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 ± × standard deviation, where 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 ∂xm/∂ϕ = (∂xm/∂ν)·(∂ν/∂ϕ) 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 13C labeling in the parameter space and is advantageous for model linearization . The constant α can be adjusted during the optimization. This is described in detail in the Appendix.
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 = Fflux(Θ) 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 so-called Buchberger's algorithm  can be applied to compute the Gröbner basis . 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()) is the null space of J() and β is a vector containing arbitrary non-zero values. The ith element of the search step whose corresponding row of null(J()) consists of zeros can be determined uniquely from the measurement η. Thus, the unique solution of the corresponding ith 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 gradient-based optimization problem is solved more accurately and with higher efficiency when analytical gradients are provided compared to numerical gradients, which are calculated by finite-difference approximation. By providing the analytical Jacobian matrix consisting of ∂xm/∂ϕ 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 Broyden-Fletcher-Goldfarb-Shanno (BFGS) formula based on linearization of the gradient .
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 user-supplied 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 1-E. Due to this, the complete process consists of a series of sub-optimization trials. At each kth 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 gradient-based 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.
SQP represents the state-of-the-art for solving constrained nonlinear optimization problems. By introducing the Lagrangian function, SQP solves a series of quadratic programming sub-problems, 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 trust-region 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 
The key concept of the developed algorithm is hybridizing the merits of these two optimization algorithms, i.e., the loose starting-point-dependency 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 1-A). 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 1-B). After a few successful trials, an upper limit (χ2UL) that schedules the initiation of STRiN is placed, e.g., half of the current objective function value, χ2UL = 1/2 f(Θ* k ), where Θ* k is the local minimizer isolated from the kth sub-optimization trial. Afterwards, if the objective value f(Θ* k ) resulting from the current SQP trial (Figure 1-C) is smaller than χ2UL, the STRiN optimization is activated (Figure 1-D); otherwise, the SQP is repeated. During the STRiN optimization, the progress of trials is monitored to prevent the non-reducing problem associated with the STRiN algorithm mentioned above. As a criterion of the progress check, the relative improvement of |f(Θ* k ) - f(Θ*k-1)|/f(Θ* k ) can be measured at each STRiN termination. If the improvement is insignificant, i.e., less than a user-specified value (e.g., 0.05), the STRiN receives penalty. In case a few successive trials show insignificant improvement, the algorithm sets a new upperlimit χ2UL 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 1-C), e.g., the absolute value of the slope resulting from the linear regression of y = (f(Θ*k-4), f(Θ*k-3),...f(Θ* k )) T with respect to x = (1, 2,..., 5) T . If the absolute slope is less than a user-specified small value for a few further trials, the optimization can be terminated ultimately.
Application to Bacillus subtilis metabolic network
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 YXS. 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 YXS were recognized as independent variables from the reduced row echelon form SRRE 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 ν2ris 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 = nflux(Θ) with 19 design parameters of Θ = (ϕ2r, ϕ3r, ϕ4r, ϕ5r, ϕ6r, ϕ7r, ϕ8r, ϕ14, ϕ14r, ϕ15r, ϕ16, ϕ17r, succinateex, acetateex, α-ketoglutarateex, fumarateex, YXS, p CO2,1, p CO2,2) T . The parameters p CO2,1 and p CO2,2 estimate the CO2 labeling pattern for each 13C labeling experiment: CO2 has only one carbon and is incorporated into other metabolites by carboxylation. Thus, its labeling state can simply be determined from other 13C labeling data as an additional parameter of NLSP in case it is not measured.
For flux re-estimation 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.
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 gradient-based optimization. The effluxes (succinateex, acetateex, α-ketoglutarateex, fumarateex) and YXS 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 13C-experimental 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 ϕ2rwas always non-zero 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 ϕ2ris 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 6-P isomerase cannot be determined from any 13C 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 6-P 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 re-estimated 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 Fflux(Θdefault) obtained by parametrization. Subsequently, the 13C labeling data of output metabolites were calculated from νdefault and Θdefault. The default fluxes were re-estimated by solving the constrained NLSP (8) whose inputs (η) were the effluxes, YXS, and the MDVs.
The worst algorithm among the local methods was the STRiN using analytical gradient and Hessian (STRiN ∇user Huser). It was rapid at the beginning but improved the objective only from 107 to 105. 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 time-inefficient 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 × 103 after 5.1 hours and GA to 2 × 102 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 re-estimation 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 starting-point-dependency 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 non-weighted sum of squares less than 2 × 10-18. As shown in Figure 4-A, 99.5 % of the flux re-estimations 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 ϕ3rand ϕ4rwere not expected to be correlated. In accordance with this, the fluxes could rarely be re-fitted in case one or both of these parameters were set arbitrarily. This indicates that the parameters ϕ3rand ϕ4rhave 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 5-B for ν3rand ν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 ∂xm/∂ϕ3rand ∂xm/∂ϕ4rare low: the objective function reached a very small value of nearly zero (f() < 1 × 10-10) and we get arbitrary but correlated values for ϕ3rand ϕ4r. Due to this, we cannot practically get ϕ3rand ϕ4restimated 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 7-phosphate were additionally measured. It was observed that providing these mass isotopomers renders the output sensitivities of ∂xm/∂ϕ3rand ∂xm/∂ϕ4rmuch less dependent on starting points. In contrast, the output sensitivities vary strongly and even become zero when the mass isotopomers of sedoheptulose 7-phosphate are not involved in the objective.
Further to this, the fluxes were not found to be influenced by the glucose 6-phosphate isomerase fluxes represented by ϕ2r. Disregarding the random values of ϕ2r, the fluxes except ν3, ν3r, ν4, and ν4rconverged to their true values. Also ϕ3rand ϕ 4rdid not show any correlation with ϕ2ras depicted in Figure 5-C and Figure 5-D, respectively. Hence, we can conclude that the observed nonlinear correlation between ϕ3rand ϕ4rbecomes visible solely from using the different starting points for the optimization. In addition, including ϕ2ras an additional design parameter gives the same result for the fluxes, whereas ϕ2rresults 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.
The tradeoff between robustness and convergence speed is a central issue in numerical optimization . 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 gradient-based 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 non-compactified case. When the optimization was performed by selecting the [0, ∞)-fluxes (ν2r, ν3r, ν4r, ν5r, ν6r, ν7r, ν8r, ν14, ν14r, ν15r, ν16, ν17rwith 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 gradient-based 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 sub-optimization 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 sub-optimization 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, gradient-based algorithms can converge rapidly but lack a global perspective for non-convex 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 non-convex problems containing several local minima and can be done in many ways [17, 39, 40]. For instance, the local minima obtained from a gradient-based 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 13C 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 1-E. As shown in Figure 4-A, 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 7-8 hours or more (by extrapolation of log10 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 (13C 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 re-estimation 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 real-case 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) .
We developed and examined a developed hybrid algorithm for numerical 13C 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 13C labeling data created from 13C-labeled 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.
For the flux re-estimation simulations of the central metabolic network of B. subtilis, the following inputs and outputs were applied.
13C Input Substrates
The choice of input tracer substrates was made by means of the D-optimality criterion . Since parallel experiments using different 13C input labels yield more information as theoretically shown for the respirometric flux analysis utilizing CO2 labeling measurement . Therefore, experimental designs for two parallel cultivations were also considered here. Among the investigated experimental designs (630 possible designs from commercially available 13C glutamate and 13C succinate species; refer to Additional file 1!), a design combining data from a 13C labeling experiment using non-labeled glutamate and [2,3-13C2] succinate with those from an experiment using [1,2-13C2] glutamate and [1,4-13C2] 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.
13C 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 5-phosphate 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.
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.
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 .
Adjustment of tolerance and parameter scaling constant
The tolerance adjustment by the restart in Figure 1-E requires a series of sub-optimization trials. The first optimization trial starts with "sufficiently relaxed" tolerances both placed on parameters and objective function, e.g., of 100. When the kth 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 10-fold decreased tolerance, where the kth 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 1-E). 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 100-fold, 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, Nneg.flux signifies the number of negative-valued fluxes, bflux the steady-state 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 non-negative. The second term is to prevent flux values calculated that may not obey stoichiometry due to the limited floating-point accuracy (round-off 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 YXS equals the strain specific precursor demand for growth (mmol precursor per gram biomass).
glucose 6-P: ν2 - (ν1 + ν2r+ 0.4217YXS) = 0
fructose 6-P: (ν2r+ ν4 + ν5 + ν6) - (ν2 + ν4r+ ν5r+ ν6r) = 0
glyceraldehyde 3-P: (ν3 + ν4r+ ν5 + 2ν6r+ ν7) - (ν3r+ ν4 + ν5r+ 2ν6 + ν7r+0.4659YXS) = 0
3-phosphoglycerate: (ν7r+ ν8) - (ν7 + ν8r+ 0.2209YXS) = 0
phosphoenolpyruvate: (ν16 + ν8r) - (ν8 + ν9 + 1.7854YXS) = 0
pyruvate: (ν9 + ν15) - (ν10 + ν15r+ 0.9964YXS) = 0
pentose 5-P: (ν1 + 2ν3r+ ν5r) - (2ν3 + ν5 + 1.1622YXS) = 0
erythrose 4-P: (ν4 + ν5r) - (ν4r+ ν5 + 0.4659YXS) = 0
sedoheptulose 7-P: (ν3 + ν4r) - (ν3r+ ν4) = 0
acetyl-CoA: ν10 - (ν11 + acetateex + 3.1585YXS) = 0
isocitrate: ν11 - ν12 = 0
α-ketoglutarate: (ν12 + ν17) - (ν13 + ν17r+ α-ketoglutarateex) = 0
glutamate: (ν17r+ glutamateex) - (ν17 + 1.4769YXS) = 0
succinate: (ν13 + ν14r+ succinateex) - (ν14 + fumarateex) = 0
oxaloacetate: (ν14 + ν15r) - (ν11 + ν14r+ ν15 + ν16 + 1.3131YXS) = 0
List of abbreviations used
- BFGS :
- GA :
- HATA :
hybrid algorithm with tolerance adjustment
- KEGG :
Kyoto Encyclopedia of Genes and Genomes
- MDV :
mass isotopomer distribution vector
- MFA :
metabolic flux analysis
- NLSP :
nonlinear least-squares problem
- SA :
- 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 trust-region method based on the interior-reflective Newton method
- STRiN∇ user H user :
STRiN using analytical gradient and Hessian
- TA :
- TK :
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.
- Wiechert W: 13C metabolic flux analysis. Metab Eng. 2001, 3 (3): 195-206.View ArticlePubMedGoogle Scholar
- Yang TH, Wittmann C, Heinzle E: Metabolic network simulation using logical loop algorithm and Jacobian matrix. Metab Eng. 2004, 6 (4): 256-267.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): 68-86.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): 121-133.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): 69-85.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: 118-135.View ArticlePubMedGoogle Scholar
- Wittmann C, Heinzle E: Genealogy profiling through strain improvement by using metabolic network analysis: metabolic flux genealogy of several generations of lysine-producing Corynebacteria. Appl Environ Microbiol. 2002, 68 (12): 5843-5859.PubMed CentralView ArticlePubMedGoogle Scholar
- Arauzo-Bravo MJ, Shimizu K: An improved method for statistical analysis of metabolic flux analysis using isotopomer mapping matrices with analytical expressions. J Biotechnol. 2003, 105 (1-2): 117-133.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): 144-156.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): 196-211.View ArticlePubMedGoogle Scholar
- Schmidt K, Nielsen J, Villadsen J: Quantitative analysis of metabolic fluxes in Escherichia coli, using two-dimensional NMR spectroscopy and complete isotopomer models. J Biotechnol. 1999, 71 (1-3): 175-189.View ArticlePubMedGoogle Scholar
- Zhao J, Shimizu K: Metabolic flux analysis of Escherichia coli K12 grown on 13C-labeled acetate and glucose using GC-MS and powerful flux calculation method. J Biotechnol. 2003, 101 (2): 101-117.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): 106-114. 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): 121-132.View ArticleGoogle Scholar
- Mendes P, Kell D: Non-linear optimization of biochemical pathways: applications to metabolic engineering and parameter estimation. Bioinformatics (Oxford, England). 1998, 14 (10): 869-883.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): 2467-2474.PubMed CentralView ArticlePubMedGoogle Scholar
- Xu P: A hybrid global optimization method: The multi-dimensional case. J Comput Appl Math. 2003, 155 (2): 423-446.View ArticleGoogle Scholar
- Long CE, Polisetty PK, Gatzke EP: Nonlinear model predictive control using deterministic global optimization. Journal of Process Control. 2006, 16 (6): 635-643.View ArticleGoogle Scholar
- Nash SG, Sofer A: Linear and Nonlinear Programming. 1996, New York , McGraw-HillGoogle 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: 102-117.Google Scholar
- Buchberger B: An algorithmical criterion for the solvability of algebraic system of equation. Aequationes Mathematicae. 1988, 4: 45-50.Google Scholar
- Saccomani MP: Some results on parameter identification of nonlinear systems. Cardiovascular Engineering: An International Journal. 2004, 4: 95-102.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, 169-184. Weinheim , Verlag ChemieGoogle Scholar
- van Winden WA, Heijnen JJ, Verheijen PJ, Grievink J: A priori analysis of metabolic flux identifiability from 13C-labeling data. Biotechnol Bioeng. 2001, 74 (6): 505-516.View ArticlePubMedGoogle Scholar
- Isermann N, Wiechert W: Metabolic isotopomer labeling systems. Part II: structural flux identifiability analysis. Math Biosci. 2003, 183 (2): 175-214.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): 734-751.View ArticlePubMedGoogle Scholar
- Boggs PT, Tolle JW: Sequential quadratic programming. Acta Numerica. 1995, 4: 1-51.View ArticleGoogle Scholar
- Coleman TF, Li Y: On the convergence of reflective Newton methods for large-scale nonlinear minimization subject to bounds. 1992, 36 p.-Ithaca, NY , Cornell Theory Center, Cornell UniversityGoogle Scholar
- Schittowski K: NLQPL: A FORTRAN-subroutine solving constrained nonlinear programming problems. Annals of Operations Research. 1985, 5: 485-500.View ArticleGoogle Scholar
- Sauer U, Hatzimanikatis V, Hohmann HP, Manneberg M, van Loon AP, Bailey JE: Physiology and metabolic fluxes of wild-type and riboflavin-producing Bacillus subtilis. Appl Environ Microbiol. 1996, 62 (10): 3687-3696.PubMed CentralPubMedGoogle Scholar
- Wittmann C, Heinzle E: Application of MALDI-TOF MS to lysine-producing Corynebacterium glutamicum: a novel approach for metabolic flux analysis. Eur J Biochem. 2001, 268 (8): 2441-2455.View ArticlePubMedGoogle Scholar
- Hill MC, Osterby O: Determining extreme parameter correlation in ground water models. Ground water. 2003, 41 (4): 420-430.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): 4723-4731.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): 869-882.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: 601-631.View ArticleGoogle Scholar
- Tahk MJ, Woo HW, Park MS: A hybrid optimization method of evolutionary and gradient search. Engineering Optimization. 2007, 39 (1): 87-104.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: 45-57.View ArticleGoogle Scholar
- Antoniewicz MR, Kraynie DF, Laffend LA, Gonzalez-Lergier J, Kelleher JK, Stephanopoulos G: Metabolic flux analysis in a nonstationary system: fed-batch fermentation of a high yielding strain of E. coli producing 1, 3-propanediol. Metab Eng. 2007, 9 (3): 277-292. 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): 86-103.View ArticlePubMedGoogle Scholar