Bioinformatics Advance Access originally published online on September 25, 2006
Bioinformatics 2006 22(22):2806-2812; doi:10.1093/bioinformatics/btl484
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Software for dynamic analysis of tracer-based metabolomic data: estimation of metabolic fluxes and their statistical analysis
1 Departamento de Bioquimica i Biologia Molecular, University of Barcelona Barcelona 08028, Catalunya, Spain
2 CERQT-Parc Cientific de Barcelona Spain
3 Department of Pediatrics, Harbor-UCLA Medical Center, Research and Education Institute Torrance, CA 90502, USA
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 SYSTEM, METHODS AND...3 IMPLEMENTATION AND DISCUSSION4 CONCLUSIONSREFERENCES |
|---|
Motivation: Metabolic flux analysis of biochemical reaction networks using isotope tracers requires software tools that can analyze the dynamics of isotopic isomer (isotopomer) accumulation in metabolites and reveal the underlying kinetic mechanisms of metabolism regulation. Since existing tools are restricted by the isotopic steady state and remain disconnected from the underlying kinetic mechanisms, we have recently developed a novel approach for the analysis of tracer-based metabolomic data that meets these requirements. The present contribution describes the last step of this development: implementation of (i) the algorithms for the determination of the kinetic parameters and respective metabolic fluxes consistent with the experimental data and (ii) statistical analysis of both fluxes and parameters, thereby lending it a practical application.
Results: The C++ applications package for dynamic isotopomer distribution data analysis was supplemented by (i) five distinct methods for resolving a large system of differential equations; (ii) the ‘simulated annealing’ algorithm adopted to estimate the set of parameters and metabolic fluxes, which corresponds to the global minimum of the difference between the computed and measured isotopomer distributions; and (iii) the algorithms for statistical analysis of the estimated parameters and fluxes, which use the covariance matrix evaluation, as well as Monte Carlo simulations.
An example of using this tool for the analysis of 13C distribution in the metabolites of glucose degradation pathways has demonstrated the evaluation of optimal set of parameters and fluxes consistent with the experimental pattern, their range and statistical significance, and also the advantages of using dynamic rather than the usual steady-state method of analysis.
Availability: Software is available free from http://www.bq.ub.es/bioqint/selivanov.htm
Contact: martacascante{at}ub.edu
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEM, METHODS AND...3 IMPLEMENTATION AND DISCUSSION4 CONCLUSIONSREFERENCES |
|---|
The introduction of substrates labeled with 13C into a metabolic network produces a vast number of isotopic isomers (isotopomers) in intermediates, which can be measured, e.g. by gas chromatography/mass spectrometry (GC-MS) and nuclear magnetic resonance (NMR). An analysis of a metabolic flux profile from steady-state isotopomer distribution data was described in Wiechert and de Graaf (1996) and further developed in Wiechert and de Graaf, 1997, Wiechert et al. (1997, 1999) and Mollney et al. (1999). Although these tools have been well documented, the requirement of using only isotopic steady-state data remains a serious barrier for investigation: even if global metabolism occurs in a steady state, the isotopic steady state can require more time than the experimental conditions would allow. Moreover, dynamic experimental data are generally more informative. Other restriction of the available tools is that the evaluated metabolic fluxes remain disconnected from the enzyme kinetic basis. To tackle these restrictions we developed a new method for the automatic construction of differential equations simulating the time course for all the isotopomers of selected metabolites (Selivanov et al., 2004). The link of this isotopomer simulation with respective kinetic model (Selivanov et al., 2005) allowed analyzing dynamic isotopomer data for determination not only fluxes but also the kinetic parameters of the analyzed pathways.
The practical application of our method for dynamic flux analysis linked to a kinetic model We mean: the application requires a statistical analysis of the estimated metabolic fluxes and kinetic parameters. The statistical analysis developed earlier specifically for flux estimation from isotopomer distribution data (Wiechert et al., 1997) cannot be used in our case, where the fluxes are complex functions of metabolite concentrations and kinetic parameters. Here we describe our implementation of the simulated annealing algorithm, which is used to identify the set of kinetic parameters and fluxes corresponding to the global minimum of the sum of square deviations from the experimental data, weighted by the respective standard deviations (
2). Our software evaluates covariance matrix for the statistical analysis of kinetic parameters. We also implemented Monte Carlo multiple optimizations for the statistical analysis of kinetic parameters and fluxes in the cases where the second derivative matrix of the merit function 2 with respect to the parameters is ill-conditioned and the covariance matrix is undefined. Although this indicates the presence of ambiguous parameter combinations, the range of the parameters, which contains the ‘true’ values, can be defined.
The statistical methods implementation in our software completes the development of a readily available tool. The binary executable files and brief user guide are offered as a free download from http://www.bq.ub.es/bioqint/selivanov.htm. They can be used without modification for the analysis of any experiment where the mass isotopomers in any of the following metabolites are measured: glucose, g6p, lactate, glutamate and glycogen. It is hoped that interaction with potential users will help further refine its functionality and interface.
|
2 SYSTEM, METHODS AND ALGORITHMS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEM, METHODS AND...3 IMPLEMENTATION AND DISCUSSION4 CONCLUSIONSREFERENCES |
|---|
The present software implementation, customized for the analysis of rat central metabolism, simulates the reaction network shown in Figure 1. This network includes the following metabolic pathways: glycolysis, gluconeogenesis, glycogen metabolism, pentose phosphate pathway (PPP), Krebs cycle and anaplerotic reactions. The effluxes (e19–e28) and influxes (i29–i31) connect the relevant parts of glucose metabolism to the rest of the cell's metabolic network. The software first estimates global fluxes and concentrations corresponding to a given set of kinetic parameters, and then, using this global solution, simulates isotopomer production in the corresponding reactions (Selivanov et al., 2004, 2005).
|
2.1 Kinetic model
Global fluxes and concentrations (not divided into isotopomer fractions) corresponding to a given set of kinetic parameters are computed by solving a system of ordinary differential equations (ODEs), using the Bulirsch–Stoer method (Stoer and Bulirsch, 1980). Time derivatives of global concentrations are described in this system as the sum of in- and out-fluxes for a given metabolite. Connection of a particular kinetic model with subsequent isotopomer analysis does not demand any specific description of enzyme reaction rates. The rates can be described, for example, via the Michaelis formula or through the elementary rate constants of the respective kinetic mechanism. The only addition to a classical kinetic model, which connects it with subsequent isotopomer distribution analysis, is mathematical description of multiple isotope exchange (Selivanov et al., 2005). An example of such multiple isotope exchanges in the case of transketolase is shown in the inset I1 of Figure 1. It includes six ‘normal’ fluxes corresponding to the three catalytic reactions, competing for the same enzyme.
 | (1) |
 | (2) |
Although each of the reaction cycles (1) changes the metabolite concentrations and therefore has to be accounted for any respective kinetic model, the reversible pass from E to EG in (2) does not change global concentrations, thereby usually it is ignored. However it changes isotopomer content and therefore is taken into account in our model.
The scheme for the transaldolase (TA) reaction is similar to that for TK except that it consists of only one reaction cycle comprising forward and reverse fluxes as well as two ‘invisible’ half-reactions. The expressions for all isotope exchange fluxes are reproduced on our website (‘Mathematica’ notebook ‘TK-competition.nb’).
The fluxes g3 and g5, representing a combination of enzymatic reactions, as well as the citrate synthase reaction (c14), are expressed in the model according to the mechanism shown in the inset I2 of Figure 1. Those reactions shown in Figure 1 that do not produce specific isotope exchange are described by the Michaelis formula.
2.2 Simulation of isotopomer production
The overall concentrations and fluxes of the metabolites, obtained from the kinetic model described above, were employed in constructing and resolving the much more complex ODE system for all isotopomers. It was assumed that all isotopomers of the same metabolite have the same biochemical properties. In this case contribution of each isotopomer into the total reaction rate is proportional to its concentration. This assumption was used to express the fluxes for each isotopomer through the respective global fluxes, thus connecting the kinetic model with isotopomer analysis as was described in details elsewhere (Selivanov et al., 2004). Here we made some modifications of the previously described algorithms that accelerated computing. The most time-consuming calculations of TK and TA bisubstrate–biproduct reactions, previously calculated as the simulated n x m interactions of n ketose isotopomers with m aldose isotopomers, were optimized to (n + 4 x m) simulations, as is described in the user guide featured on our website.
To solve the huge system describing dynamics of isotopomer distribution, five different methods were implemented using the algorithms of Press et al. (2002) and Milde (1998, http://www.minet.uni-jena.de/www/fakultaet/iam/ode??/ode??_e.html): two for non-stiff systems (fourth-order Runge–Kutta and Bulirsch–Stoer); and three for stiff systems [Rosenbrock, semi-implicit Bulirsch–Stoer and backward differentiation formula (BDF)].
2.3 Finding the optimal set of parameters and fluxes
Fitting the simulation to an experimental distribution of isotopomers reveals the set of parameters used in rate equations and flux profile corresponding to analyzed conditions. This iterative multi-step optimization procedure exploits the simulated annealing method of global optimization (Kirkpatrick et al., 1983) and Powell's method for downhill descent in a local area (Brent, 1973). One optimization step consists of (i) the kinetic model execution to evaluate the global fluxes corresponding to the given set of parameters; (ii) computing the respective isotopomer distribution dynamics in the linked isotopomer module, which uses the obtained global fluxes and concentrations; and (iii) comparison of the calculated and experimental pattern, 2 calculation and parameter change for the next step of optimization.
The change of parameters could be either random perturbations as simulated annealing optimization suggests or subsequent change of the chosen parameter in the direction, which decreases 2. The step with random perturbation is followed by a number of steps, which perform coordinate descent according to the Powell's method, changing all the parameters one by one until the local minimum in all coordinates is reached.
Then the algorithm performs random perturbation and descent to the new local minimum again. Finally it finds the set of parameters and fluxes corresponding to the global minimum of the 2 merit function.
The choice of optimization algorithm was defined by the necessity of underdetermined systems analysis, which restricts the number of applicable algorithms (Dennis and Schnabel, 1983; Gill et al., 1981). In particular, methods that use the second derivative matrix of 2 with respect to parameters (Hessian matrix), such as the Levenberg–Marquardt method (Marquardt, 1963) cannot be applied.
2.4 Symmetrical reorganization of labels
Fumarate molecules are symmetrical; their reversible hydration to malate does not distinguish between symmetrical positions; i.e. hydroxyl could be bound to the second or third atoms with equal probabilities. If the malate carbons are numbered in such a way that the hydroxyl is in position two, then this numbering (and hence, the labeling pattern) is reversed in 50% of the sample, compared with that of the non-symmetrical precursors (e.g. 2-oxoglutarate). In this way, the transformation of symmetrical substances reorganizes isotopomers such that the concentrations of mutually symmetrical pairs become equal. These symmetrical pairs are as follows: 0001(1)1–1000(8), 0010(2)–0100(4), 0011(3)–1100(12), 0101(5)–1010(10), 0111(7)–1110(14) and 1011(11)–1101(13). The algorithm summarizes the concentrations pairwise and redistributes the sum equally between the components. The remaining four of the total sixteen isotopomers are symmetrical—[0000(0), 1001(9), 0110(6) and 1111(15)]—i.e. equal to their symmetrical matches. Their further transformation will lead to the same labeling of oxaloacetate, whatever position the hydroxyl may take. Thus, they do not need any special treatment in simulations. In this way, the concentration of any labeling pattern in the products derived from fumarate equals that of its symmetrical pattern. Although the transformation from fumarate equals the concentrations in mutually symmetrical pairs, other pathways of oxaloacetate production (e.g. by pyruvate carboxylase) do not. However, if reversible transformation from oxaloacetate to fumarate is much faster than non-symmetrical pathways of oxaloacetate production, one can regard the entire pool of oxaloacetate to be symmetrically reorganized, an option implemented in the current version of our software. In principle, both symmetrical and non-symmetrical pathways of oxaloacetate production can be specifically taken into account, although detailed study of their relative contribution is outside the scope of the present contribution.
2.5 Statistical analysis
The fact that the Hessian matrix is ill-conditioned indicates that the experimental data could not distinguish between the different values of some kinetic parameters (Press et al., 2002), since different parameter sets match the former equally well. Although the Levenberg–Marquardt method fails because the Hessian matrix is ill-conditioned, Powell's method of coordinate descent identifies a set of parameters corresponding to the 2 minimum, although other sets could potentially fit the data equally well. In order to define the set of parameters unambiguously certain parameters should either be (i) eliminated or (ii) re-set as constants. Although the singular value decomposition, also exploited in our software, can define the maximum detectable number of parameters the model should contain, both of these simplifications are not appropriate. First, since all the used parameters have definitive physical meaning (Vmax, Km, rate constants of elementary steps) and simplification would involve substituting these known characteristics by their combinations, which cannot be clearly interpreted. Second, to assign fixed values to some of the parameters, it is necessary to know the respective values.
However, despite the formal multiplicity of the solutions, the interval in the space of parameters, in which different points yield the same 2-value, might remain small. In this particular case, there is no practical difference between identifying a reliable set of parameters and finding a reliable interval within the space of parameters, within which various points are consistent with the experimental pattern. The simulated annealing algorithm identifies the global minimum for the merit function 2 and the respective set of parameters and fluxes. If the system is underdetermined, a repetition of the optimization procedure could reveal a different set of fluxes corresponding to the same minimum of merit function. By multiple repetitions of the optimization procedure, it then defines the intervals of parameters and fluxes in which different sets correspond to the same 2-value.
Thereby, the software checks the goodness of fit based on the 2 criterion as it is described in Press et al. (2002). The probability that a model with 
degrees of freedom is correct and 2 could exceed the value given above by chance, is given as an incomplete gamma function, which is normally used for estimating goodness of fit (Press et al., 2002):
 |
If the Q-value is acceptable, the true values of fluxes are expected to be included within the boundaries of optimal sets found by the multiple repetition of optimization procedure, the usual assumption for modeling of experimental data. Inside the interval found, the probability of values to be included in an optimal set by the optimization procedure is distributed close to normal, as is explored in the next section. This probability is expected to coincide with the probability of finding the true values of fluxes and since it is close to normal, it can be presented as the mean and standard error, thus summarizing all the values obtained from independent Monte Carlo simulations.
|
3 IMPLEMENTATION AND DISCUSSION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEM, METHODS AND...3 IMPLEMENTATION AND DISCUSSION4 CONCLUSIONSREFERENCES |
|---|
3.1 Analyzed experimental data
In the present example, we used the GC-MS data pattern obtained from the liver cells of fasting rats. The cells were incubated for 2 h in media containing 50% [1,2-13C2]D-glucose (the experimental conditions are described in the legend to Table 1). The isotopomer pattern was obtained separately from eight animals, as shown in Table 1.
|
Goodness of fit. Table 1 also summarizes 200 different optimizations started from initial values of the parameters chosen randomly. Each optimization ended with the
2-value at its minimum (
250), which summarized the square deviations from the 31 mass isotopomer fractions measured in each of eight animals. The total number of parameters, which was subjected to optimization, was 43. Thus, the total number of degrees of freedom was (31 x 8)–43 = 205. For our case, the Q-value was 0.264. Although there is no strict threshold for model acceptance, it is not uncommon to deem Q > 0.001 as acceptable (Press et al., 2002). Thus, the above Q-value indicates a good fit, which is visually confirmed by Table 1 and shows that the computed fractions are close to the ones obtained experimentally. In all cases, the computed values were within the interval obtained experimentally for different animals.
3.2 Intervals for optimal values of the parameters
After finding the optimum 2 merit function and estimation of the goodness of fit, the expected next step is estimating the standard deviations of the parameters by calculating the second derivative matrix of the 2 merit function with respect to the parameters, its inversion and finding the covariance matrix. Although in principle this method is implemented by our software, in the example considered here the singular value decomposition, when applied to the Hessian matrix, had shown that it was ill-conditioned, close to singular. This indicated that ambiguous combinations of parameters, which could give approximately the same fit, could exist. Therefore the Monte Carlo global optimizations, resulted in the same 2-value, were repeated numerous times, each time recording that set of parameters which corresponds to the minimum value of merit function. Table 2 shows the mean fluxes and the standard errors of mean, summarizing 200 optimizations.
|
3.3 The type of distribution for the parameters and fluxes
The optimizations yielded a spectrum of values for each flux; this spectrum can be presented as a probability distribution of values falling within fixed intervals. This probability distribution could be fit by a Gaussian curve, thus indicating that the distribution of the probability of parameter being included in an optimal set was close to normal. Figure 2 shows such probability distributions for some fluxes. Although the formal Kolmogorov–Smirnov test shows that not all of the computed fluxes have normal distributions, the examples shown in Figure 2 nevertheless present the most common situation: all the distributions have one maximum, such that the mean value corresponds to the maximum density of probability and the formal standard error characterizes its uncertainty. Thus, instead of finding the exact point corresponding to the minimum of
2, we can find, in the case of the ill-conditioned Hessian matrix, the most probable value and its standard error. The following example shows that the standard error of mean in the sets corresponding to a particular experimental pattern could be small enough to distinguish even between metabolic fluxes in cells from individual animals under the same experimental conditions. This indicates that although formal criteria demonstrate the possibility of ambiguous combinations of parameters, for the method implemented the difference from the case of unambiguous fit is practically insignificant.
|
3.4 Analysis of the data for individual animals
In principle, the data collected for each animal allow us to define peculiarities in individual flux distributions, if the defined flux and parameter areas are different. The fit of the mass isotopomer data obtained individually from two animals chosen at random is shown in Table 1. The respective fluxes are shown in Table 2. It should be noted that despite the apparent similarity in the data from different animals, the
2-test clearly distinguishes them: simulation of the data for the second animal using the best-fit parameters for the first one increases the 2-value 6-fold. Therefore, the resulting flux intervals are different; formal application of Student's t-test showed that the means for the selected animals differed significantly not only from each other, but also from the best-fit for the entire group, defined as the minimum sum of deviations from all experimental patterns. Evidently, the values obtained in such a way for a group of animals do not necessary include all individual variations within the group. Determining the parameter and flux profiles for each individual, and then statistically analyzing the resulting individual values, would be more representative of the metabolism variation in a group of experimental animals. Thus, even when the experimental data are insufficient to define the model parameters and fluxes unambiguously, the model could be sensitive even to individual variations within the group of experimental animals.
3.5 Dynamic characteristics of isotopomer accumulation
It is commonly assumed that isotopic and metabolic steady-state conditions are met when cells undergo incubation with isotope labeled precursor for 2 h, as in the present experiment. This raises a question regarding the necessity of using dynamic analysis when steady-state analysis has already been fully elaborated, including its statistical aspects (Wiechert et al., 1997). In fact, the time course to the isotopic steady state could differ for various isotopomers of the same substance. Moreover, the isotopomer analysis could often be made only for extracellular metabolites, where due to dilution the relative rate of isotopomer accumulation could be extremely slow. When an extracellular metabolite exchanges with the intracellular one, this exchange complicates the dilution. As it may not be possible to verify whether the steady state has been reached, dynamic simulation of isotopomer accumulation may provide an answer to this question. Figure 3 reproduces a time course for the intracellular overall hexose concentration and the extracellular glutamate containing four labeled atoms, each calculated for the set of parameters, which fits the measured isotopomer pattern. Although the total intracellular concentrations reach a steady state within the first 2 min of incubation, isotopomer fractions in glutamate do not approach steady state in 2 h of incubation. Therefore, the steady-state approach was not appropriate for analyzing the results from this type of experiment.
|
Thus, even for rapidly metabolizing liver cells, 2 h is not sufficient to reach the isotopic steady state for the measured metabolite. On the other hand, it is difficult to maintain the same global metabolic fluxes during longer incubation periods. Therefore, the tool presented in this study would be the only means of analyzing such experiments.
|
4 CONCLUSIONS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEM, METHODS AND...3 IMPLEMENTATION AND DISCUSSION4 CONCLUSIONSREFERENCES |
|---|
We present the development of a new method for stable isotope tracer data analysis. The software algorithms implemented included (i) various methods to automatically construct and numerically resolve the vast system of differential equations for isotopomers; (ii) a fitting algorithm based on ‘simulated annealing’ global minimization of the deviations from experimental data; and (iii) statistical analysis of estimated parameters and fluxes, which makes this new tool applicable in practice. Several classical routines for statistical analysis are incorporated in our software. Singular value decomposition for the Hessian matrix tests its singularity (or ill-conditioning). If possible, it defines the covariance matrix and standard deviations of the parameters. In the case of singularity or ill-conditioning, we implemented Monte Carlo multiple optimizations, thereby specifying the interval wherein various combinations of parameters yielded the same experimental data fit. The present example shows that despite singularity, parameter and flux intervals can remain sufficiently small to reveal individual variability between experimental animals.
One of the main advantages of this new software tool is the possibility of conducting dynamic analysis. The accumulation of certain isotopomers could be so slow that it eliminates any possibility of achieving the isotopic steady state. Therefore, the steady-state approach is not applicable for numerous experimental situations.
The software presented on our website can be used for similar types of analysis involving isotopomer distribution data from different experiments.
|
Acknowledgments |
|---|
This work was supported by the grants: European Community Foundation FP6-2005- 037939 Biobridge; Fundacion la Caixa (ONO3-70-0) Ministerio de Educacion y Ciencia, of Spanish Government (SAF2005-01627 and AGL2004-07579-C04-03/ALI, UNBA05-33-001); European Union FEDER funds NIH DK56090-A1 (to W.N.P.L); and Generalitat de Catalunya 2005 PEIR 0051/69 and (ABM/acs/PIV2002-32) (to V.A.S). The GC/MS Facility is supported by PHS grants P01-CA42710 to the UCLA Clinical Nutrition Research Unit, Stable Isotope Core and M01-RR00425 to the General Clinical Research Center. The Bioinformatic grant program of the Foundation BBVA.
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
Associate Editor: John Quackenbush
1Here 0 states for non-labeled carbon atom (12C), 1 states for labeled carbon atom, so an isotopomer is present as a binary number (its decimal equivalent is given in parentheses), which in the same time is its reference number in the array of concentrations. 
Received on July 25, 2006; revised on August 31, 2006; accepted on September 14, 2006
|
REFERENCES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEM, METHODS AND...3 IMPLEMENTATION AND DISCUSSION4 CONCLUSIONSREFERENCES |
|---|
Brent, R.P. Algorithms for Minimization without Derivatives, (1973) , NJ Prentice-Hall, Englewood Cliffs.
Dennis, J.E. and Schnabel, R.B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations, (1983) , NJ Prentice-Hall, Englewood Cliffs.
Gill, P.E., Murray, W., Wright, M.H. Practical Optimization, (1981) , NY Academic Press.
Kirkpatrick, S., et al. (1983) Optimization by simulated annealing. Science, 220, 671–680
Marin, S., et al. (2004) Dynamic profiling of the glucose metabolic network in fasted rat hepatocytes using [1,2-13C2]glucose. Biochem. J, . 381, 287–294[CrossRef][Web of Science][Medline].
Marquardt, D.W. (1963) An algorithm for least squares estimation of non-linear parameters. J. Soc. Indust. Appl. Math, . 11, 431–441[CrossRef].
Mollney, M., et al. (1999) Bidirectional reaction steps in metabolic networks: IV. Optimal design of isotopomer labeling experiments. Biotechnol. Bioeng, . 66, 86–103[CrossRef][Web of Science][Medline].
Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T. Numerical Recipes in C: The Art of Scientific Computing, (2002) , NY Cambridge University Press.
Selivanov, V.A., et al. (2004) An optimized algorithm for flux estimation from isotopomer distribution in glucose metabolites. Bioinformatics, 20, 3387–3397
Selivanov, V.A., et al. (2005) Rapid simulation and analysis of isotopomer distributions using constraints based on enzyme mechanisms: an example from HT29 cancer cells. Bioinformatics, 21, 3558–3564
Stoer, J. and Bulirsch, R. Introduction to Numerical Analysis, . (1980) , Springer-Verlag, NY.
Wiechert, W. and de Graaf, A.A. (1996) In vivo stationary flux analysis by 13C labeling experiments. Adv. Biochem. Eng. Biotechnol, . 54, 109–154[Medline].
Wiechert, W. and de Graaf, A.A. (1997) Bidirectional reaction steps in metabolic networks: I. Modeling and simulation of carbon isotope labeling experiments. Biotechnol. Bioeng, . 55, 101–117[CrossRef].
Wiechert, W., et al. (1997) Bidirectional reaction steps in metabolic networks: II. Flux estimation and statistical analysis. Biotechnol. Bioeng, . 55, 118–135[CrossRef].
Wiechert, W., et al. (1999) Bidirectional reaction steps in metabolic networks: III. Explicit solution and analysis of isotopomer labeling systems. Biotechnol. Bioeng, . 66, 69–85[CrossRef][Web of Science][Medline].
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||