Bioinformatics

Bioinformatics Advance Access originally published online on July 7, 2005
Bioinformatics 2005 21(17):3558-3564; doi:10.1093/bioinformatics/bti573

This Article Abstract FREE Full Text (Print PDF) All Versions of this Article: 21/17/3558    most recent bti573v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Search for citing articles in: ISI Web of Science (4) Request Permissions Google Scholar Articles by Selivanov, V. A. Articles by Cascante, M. Search for Related Content PubMed PubMed Citation Articles by Selivanov, V. A. Articles by Cascante, M. Social Bookmarking          What's this?

© The Author 2005. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions{at}oupjournals.org

Rapid simulation and analysis of isotopomer distributions using constraints based on enzyme mechanisms: an example from HT29 cancer cells

Vitaly A. Selivanov ¹, Ludmilla E. Meshalkina ², Olga N. Solovjeva ², Philip W. Kuchel ³, Antonio Ramos-Montoya ¹, German A. Kochetov ², Paul W. N. Lee ⁴ and Marta Cascante ^1,*

¹Departamento de Bioquimica i Biologia Molecular, Facultat de Quimica and CERQT at Parc Cientic de Barcelona Barcelona, Catalunya, Spain
²A.N.Belozersky Institute of Physico-Chemical Biology, Moscow State University 199899 Moscow, Russia
³School of Molecular and Microbial Biosciences, University of Sydney NSW, Australia
⁴Department of Pediatrics, Harbor-UCLA Medical Center, Research and Education Institute Torrance, CA 90502, USA

^*To whom correspondence should be addressed.

	Abstract

TOP Abstract INTRODUCTION SYSTEMS AND METHODS ALGORITHM IMPLEMENTATION AND DISCUSSION CONCLUSIONS REFERENCES

 

Motivation: Addition of labeled substrates and the measurement of the subsequent distribution of the labels in isotopomers in reaction networks provide a unique method for assessing metabolic fluxes in whole cells. However, owing to insufficiency of information, attempts to quantify the fluxes often yield multiple possible sets of solutions that are consistent with a given experimental pattern of isotopomers. In the study of the pentose phosphate pathways, the need to consider isotope exchange reactions of transketolase (TK) and transaldolase (TA) (which in past analyses have often been ignored) magnifies this problem; but accounting for the interrelation between the fluxes known from biochemical studies and kinetic modeling solves it. The mathematical relationships between kinetic and equilibrium constants restrict the domain of estimated fluxes to the ones compatible not only with a given set of experimental data, but also with other biochemical information.

Method: We present software that integrates kinetic modeling with isotopomer distribution analysis. It solves the ordinary differential equations for total concentrations (accounting for the kinetic mechanisms) as well as for all isotopomers in glycolysis and the pentose phosphate pathway (PPP). In the PPP the fluxes created in the TK and TA reactions are expressed through unitary rate constants. The algorithms that account for all the kinetic and equilbrium constant constraints are integrated with the previously developed algorithms, which have been further optimized. The most time-consuming calculations were programmed directly in assembly language; this gave an order of magnitude decrease in the computation time, thus allowing analysis of more complex systems. The software was developed as C-code linked to a program written in Mathematica (Wolfram Research, Champaign, IL), and also as a C++ program independent from Mathematica.

Results: Implementing constraints imposed by kinetic and equilibrium constants in the isotopomer distribution analysis in the data from the cancer cells eliminated estimates of fluxes that were inconsistent with the kinetic mechanisms of TK and TA. Fluxes measured experimentally in cells can be used to estimate better the kinetics of TK and TA as they operate in situ. Thus, our approach of integrating various methods for in situ flux analysis opens up the possibility of designing new types of experiments to probe metabolic interrelationships, including the incorporation of additional biochemical information.

Availability: Software is available freely at: http://www.bq.ub.es/bioqint/selivanov.htm

Contact: martacascante{at}ub.edu

	INTRODUCTION

TOP Abstract INTRODUCTION SYSTEMS AND METHODS ALGORITHM IMPLEMENTATION AND DISCUSSION CONCLUSIONS REFERENCES

 
‘Systems analysis’ of biochemical networks using isotopomer-distribution data encounters many problems, not only experimental, but also theoretical. In our previous publication (Selivanov et al., 2004) the problem of computation time for automatic construction and solving of huge systems of ordinary differential equations (ODEs) was addressed, when each metabolite consisting of n carbon atoms is represented by 2^{n 13}C carbon isotopomers. The previously presented method was faster than that of ‘isotopomer-mapping matrices’ (Schmidt et al., 1997) and in fact addressed the problem of stiffness of the resulting ODEs.

Another problem of metabolic flux analysis is reliability of the obtained set of fluxes and the main objective of the present work was to address it. Cellular metabolism of ¹³C labeled substrates introduces ¹³C into many metabolic intermediates which can be detected and quantified by NMR as distribution of individual isotopomers or by gas chromatography/mass spectrometry (GC/MS) as distribution of mass isotopomers. This distribution reflects respective metabolic fluxes; however, the obtained pattern in a particular experiment is often insufficient for reconstructing the profile of metabolic fluxes. The present work shows that without additional information, analysis of the isotopomer distribution evaluates large range of estimated fluxes compatible with a specific isotopomer pattern. Implementing the detailed kinetic mechanisms of transketolase (TK, EC 2.2.1.1 [EC] ) and transaldolase (TA, EC 2.2.1.2 [EC] ) provides additional information to restricting the range of flux profiles by the ones that are compatible not only with a given isotopomer pattern, but also with the data of previous biochemical studies. On the contrary, using the classical information makes available its validation for in vivo conditions of the label-distribution experiments.

TK and TA mediate multiple isotope-exchange reactions, which despite contributing to the isotopomer patterns, in the past have largely been ignored in metabolic flux analysis; exceptions are the works of Berthon et al. (1993) and van Winden et al. (2001). TK is a thiamine diphosphate (ThDP)-dependent enzyme that catalyzes cleavage of a carbon–carbon bond and reversibly transfers a two-carbon fragment (glycolaldehyde) from a ketose-phosphate donor to an aldose-phosphate acceptor. TA catalysis operates by a slightly different mechanism, involving a Schiff base formed directly between the enzyme and the substrate, and a three-carbon fragment (dihydroxyacetone) is transferred from a ketose-phosphate donor to an aldose-phosphate acceptor.

Reversible scission and reformation of the same sugar phosphate via TK or TA brings about isotope exchange without net formation of a particular sugar phosphate; these isotopomer exchange reactions are listed in Table 1. Since these reactions lead to the redistribution of labeled atoms, omitting them from a simulation scheme would result in the incorrect estimation of all fluxes in the network. In our previously developed software, this isotope exchange was independent, as is commonly accepted in flux analysis (Wiechert and de Graaf, 1996). In this way, it increased the number of kinetic parameters thus magnifying the problem of ambiguity of estimates that were consistent with the experimental isotopomer distributions. As shown here, using the general information on the reaction mechanism restricts the space of possible solutions. This constitutes the first step in developing learning software that will analyze new experiments using previously accumulated data and models.

View this table: [in this window] [in a new window]   Table 1 Simulation of 13C distribution measured in Rib5P of RNA from HT29 cells utilizing 100% [1,2-13C]D-glucose incubated in the absence (Control) or in the presence of oxithiamine (OT), which induced 75% inhibition of TK and 50% of accumulation of the cells during 72 h of incubation

 
The new software was used to estimate fluxes in glucose metabolism in HT29 cells that were incubated with [1,2-¹³C]D-glucose; the cells were studied in the absence and presence of oxythiamine (OT), an irreversible inhibitor of TK. As it was found experimentally, cell proliferation is highly dependent on TK activity (Comin-Anduix et al., 2001); since the activity of this enzyme is associated with central metabolism, the objective of the present example of software application was to reveal the changes in metabolic profile of central metabolism induced by TK inhibition by incubating HT29 cancer cells with OT. This study was assumed to establish the link between TK inhibition and cell proliferation.

	SYSTEMS AND METHODS

TOP Abstract INTRODUCTION SYSTEMS AND METHODS ALGORITHM IMPLEMENTATION AND DISCUSSION CONCLUSIONS REFERENCES

 
The six-step sequential BiBi kinetic mechanisms of the TK and TA reactions, expressed through the unitary rate constants using the method of King and Altman (e.g. Mulquiney and Kuchel, 2003; Cornish-Bowden, 2004) was integrated in the software for the analysis of isotopomer distributions of the metabolites of the glycolytic and pentose phosphate pathways. Two versions of the software were developed. One is based as before on Mathematica with a program written in C linked to it; this link enables high-speed calculations, further optimized and reprogrammed in assembly language. The second program was written entirely in C++ with some parts in assembly language. This program is presented as Windows and also console application with explicit-dialogue capability; it solves the differential equations for overall metabolite concentrations using the C++ class library ‘ODE++’ of Milde (2005) (http://www.minet.uni-jena.de/www/fakultaet/iam/ode++/ode++_e.html) using a series of Runge–Kutta and backward difference methods, and calculates isotopomer concentrations using the same algorithms as the Mathematica program.¹ The procedure for parameter optimization adapted to the present analysis and estimation of the confidence intervals was based on ² minimization using coordinate descent.

	ALGORITHM

TOP Abstract INTRODUCTION SYSTEMS AND METHODS ALGORITHM IMPLEMENTATION AND DISCUSSION CONCLUSIONS REFERENCES

 
Only the enhancements relating to the main objective, integration of kinetic modeling with isotopomer analysis, are given here. A detailed description of the numerical algorithm for isotopomer distribution analysis has already been presented (Selivanov et al., 2004).

Overall concentrations and fluxes. The total concentration of each metabolite is given by the solution of the ODEs that relate to the scheme presented in Figure 1, where the rate of change of each concentration is described as the sum of production rate of the given substance minus the rates of its consumption. The reaction rates are described mainly as the Michaelis–Menten equations, with the exception of the TK reactions 13–15 and TA reaction 16. These were described according to the reaction mechanism, shown in the Figure 2A, through the unitary rate constants, denoted as v_i and expressed according to the law of mass action as a product of the constant k_i, the respective enzyme species, and the substrate, where applicable, i.e. v₁ = k₁ x E x c₁. The overall reaction flux was described as the difference between the forward and reverse unitary rates for the same reaction step (v_f – v_r), which is the same for all of the six steps shown in the scheme (i.e. v_net = v₁ – v_–1) assuming that the enzyme is always under the condition of a quasi steady-state. Expressions for the unitary rates require a description of the concentrations of particular enzyme species, which are the product of the total enzyme concentration (e₀) and the fraction of the species. The fraction was expressed through the unitary rate constants by using the method of King and Altman (Cornish-Bowden, 2004; implemented in Mathematica in Mulquiney and Kuchel, 2003).

View larger version (19K): [in this window] [in a new window]   Fig. 1 Scheme of the metabolic reactions simulated in the model. The nodes represent metabolites, solid lines are reactions and numbers are the reference numbers of respective reaction modifier in the array of parameters (NV) used by the program. Several new reactions are added to the pathway considered in Selivanov et al. (2004): exchange between intracellular Glc6P and extracellular glucose, dependent on the parameters (1) and (0); exchange between the intracellular trioses and outside lactate, dependent on the parameters (7) and (8); efflux of GraP to the Krebs cycle (10), glycogen production (9) and influx of the medium fructose through the pool of trioses (11). The number of reactions corresponds to the number of parameter designed for fitting respective Vmax. The rate equations for all the reactions could be found in the Mathematica notebook, presented in the website.

 

View larger version (27K): [in this window] [in a new window]   Fig. 2 (A) scheme of the TK reaction, which consists of a series of reversible unitary steps, can be divided into two parts. The steps of first half-reaction are: binding of donor substrate (c1) and formation (k1, k–1) of a covalent enzyme–substrate complex (Ec1); splitting (k2, k–2) of donor substrate and formation of a covalently bound intermediate (the -carbanion of ,ß-dihydroxyethyl-ThDP—so-called ‘active glycolaldehyde’) and an aldose; both are localized in the active site of the enzyme (EGa1). This complex dissociates (k3, k–3) to the complex of the enzyme with active glycolaldehyde (EG) and the first product—free aldose (a1). In the second half-reaction, active glycolaldehyde interacts with the other aldose added to the reaction mixture (k4, k–4) and the new ketose is released after passing the same reaction steps in reverse order (k5, k–5 and k6, k–6). (B) CD spectra of holotransketolase isolated from baker's yeast according to the method described earlier (Wikner et al., 1994). The spectra were registered in 50 mM glycyl-glycine buffer, pH 7.6, containing 2.5 mM CaCl2, 30 µM ThDP, 4.2 µM TK. (1) TK without substrates, (2) TK after addition of 2 mM hydroxypyruvate and (3) TK after a series of simultaneous addition of equal amounts of Xul5P and GraP up to 180 µM with 15 µM concentration intervals. All additions after the first one did not change the spectral curve. (C) Evaluation of the dissociation constant of enzyme complex with Xul5P. Initial concentrations are: Rib5P, 800 µM, Xul5P, 25 µM (solid square) and 37 µM (open square). Calculated lines coincide with experimental data at K1 = 67.5 µM ([Xul5P]t=0 = 24 µM) and K1 = 69.5 µM ([Xul5P]t=0 = 37 µM). K6 = 2.7 mM for both curves that coincide with the experimental data. Catalytic constants for all curves are k2 = k–2 = k5 = 5700/min–1, k–5 = 1770/min. The product accumulation was measured spectrophotometrically at a wavelength of 340 nm on a Aminco DW-2000 spectrophotometer at 25°C by the rate of NAD+ reduction in a coupled system with D-glyceraldehyde-3-phosphate dehydrogenase (Kochetov, 1982). The reaction mixture (final volume 1 ml) contained: 50 mM glycyl-glycine, pH 7.6, 1 mM sodium arsenate, 0.37 mM NAD+, 3U D-glyceraldehyde-3-phosphate dehydrogenase, 3.2 mM dithiothreitol, 2.5 mM MgCl2, 0.08 mM ThDP, Xul5P, Rib5P and 0.02 U transketolase. (D) Evaluation of dissociation constant of the enzyme complex with Rib5P. Initial concentrations are: Xul5P, 200 µM, Rib5P, 10 µM (solid square), and Xul5P, 35 µM, Rib5P, 28 µM (open square). Measurement and other conditions are as in (C). Transketolase was homogenous by SDS–Na-electrophoresis; its specific activity was 20 U/mg. (E) A table summarizing the determined ratios of unitary rate constants.

 
For example, the fraction of the total enzyme, i.e. Ec1 is given by

(1)

where D_Ec1 is the tree of the Ec₁-node:

(2)

and is the sum of the trees, constructed analogously to D_Ec1, for all intermediate species.

Numerical solving of the relevant ODEs (presented as function ‘SS’ in the file ‘rib-ot.nb’) gives time courses of relaxation from initial values to a steady state for the metabolites and fluxes. The overall steady state values that were attained then allowed the automatic formulation of a huge system of differential equations for all the isotopomers [as was described in our previous work (Selivanov et al., 2004)].

Isotope exchange fluxes for isotopomer analysis. The ODEs for isotopomers contain similar flux terms as the respective equations for overall concentrations; the fluxes at this step correspond to influxes and effluxes from separate isotopomers. Moreover, these equations account separately for forward and reverse fluxes and isotope exchange, which do not lead to the change of overall concentrations (Table 1). Overall values of these fluxes aredescribed next.

As could be deduced from the reaction scheme for TK (Fig. 2A), the enzyme catalyzes two types of fluxes that differently contribute to isotope exchange. First, forward flux of the last n – 2 atoms of substrates (ketose; denoted c for Spanish ‘cetosa’) c₁ to (aldose) a₁ (e.g. Xul5P GraP) implies delivering the atoms through the three steps (c₁ + E Ec₁ EGa₁ EG + a₁). The intermediates Ec₁ and EGa₁ contain aldose fragments derived from two sources, c₁ and a₁. Quantification of the unidirectional steady-state flux of c₁ into a₁ via Ec₁ and EGa₁ requires the determination of the origin of transferred fragments in these two populations of enzyme complex. These fractions are fundamentally determined by the relative values of the respective unitary rates (described above). Specifically, the rate of delivery of c₁ atoms into a₁ is a part of the rate v₃; it is proportional to the concentration of c₁ atoms in EGa₁. The proportionality constant or fraction is denoted by , where the superscript ca denotes the ketose-to-aldose reaction direction. In turn, depends on the fraction of c₁ atoms in Ec₁, that partly consists also of a₁ atoms that enter via the reactions whose rates are v_–2 and v_–3; thus it is expressed as a ratio of the input of the donor atoms from Ec₁ () to the total input to EGa₁:

(3)

The proportion of atoms in Ec₁ that originated from c₁ () in turn is given by the ratio of influx of this kind of atom to the total influx into the compound Ec₁ at steady state:

(4)

Solving Equations (3) and (4) yields the expression:

(5)

The rate of the flux of c₁ atoms into a₁, v_if, where the subscript if denotes the ‘invisible forward’ reaction, is given by

(6)

Equation (6) gives the rate of forward delivery of the last n – 2 atoms in c₁ to a₁, expressed by using the unitary rates. Since these atoms can originate only from either c₁ or a₁, the fraction of atoms originating from a₁ is expressed as . Therefore, the reverse flux of the last n – 2 atoms of the aldose (a₁) to the ketose (c₁) can be described similarly to Equation (6) as

(7)

The exchange of atoms between c₂ and a₂ was simulated analogously.

Second, forward flux (v_f) of the first two atoms of c₁ to a second ketose/donor, c₂, implies delivering the atoms through six steps (c₁ Ec₁ EGa₁ EG EGa₂ Ec₂ c₂). This is a part of the rate of c₂ production (v₆) and it is proportional to the fraction of c₁ in Ec₂, namely , where the superscript c1 denotes its origin from c₁. This proportion is determined similarly to that described above, by solving the five equations for the fractions of atoms that originated from c₁ in all the species (analogously to the Equations (3) and (4)).

The reverse flux (v_r) of the first two atoms of c₂ to c₁ was simulated similarly. The details of these calculations are given in the Mathematica notebook ‘TK-scheme.nb’ on our website.

Since the main objective of the present work was to implement kinetic information (both from models and experiments) in isotopomer analysis, an example of experiments that can decrease the number of unknown parameters is described next.

Estimating values of unitary rate constants
The ratio of unitary rate constants in TK (k₂/k_–2 or k_–5/k₅) can be measured by using circular dichroism (CD) spectra of mixtures of the enzyme and substrate (for experimental details see the caption of Fig. 2). Specifically, Curve 1 in Figure 2B characterizes the negative band in the wavelength range 300–380 nm in the spectrum of the protein without bound substrates (E in Figure 2A). This band is sensitive to the presence of bound donor substrates (Heinrich et al., 1971; Kochetov et al., 1973; Kochetov, 1982) and Curve 2 characterizes the spectrum of EG, i.e. the only enzyme form that exists under conditions when hydroxypyruvate is used as a substrate; it is the only one that makes TK reaction to be irreversible. A similar spectrum is recorded if the reversibly converted ketose substrate (c₁ in Fig. 2B) is added and aldose product (a₁) is removed from the reaction mixture.

In the presence of ketose (Xul5P) and aldose (GraP) substrates, the CD band (Curve 3, Fig. 2B) occupies an intermediate position between those for E (Curve 1) and EG (Curve 2, Figure 2B). In this case, the enzyme is distributed among the four intermediate species (E, Ec₁, EGa₁, and EG; Fig. 2A). Increase of c₁ concentration in the reaction mixture induces transformation of the E species into Ec₁, whereas an increase of a₁ concentration induces conversion of the EG species into EGa₁.

It was found that a simultaneous increase of both Xul5P and GraP concentrations from 15 µM to 180 µM did not produce any change in the CD spectrum of the enzyme (Curve 3, Fig. 2B). The absence of change in the CD spectrum, despite the transition from E to Ec₁ and from EG to EGa₁, signified that E had the same spectrum as the Ec₁ and the EG had the same spectrum as the EGa₁. Therefore, the spectrum of Ec₁ coincided with the Curve 1 and EGa₁ with Curve 2. Under the conditions of excess of both substrates, when only two enzyme species (Ec₁ and EGa₁) would be present, they would be distributed according to the ratio of forward and reverse catalytic constants (i.e. k₂/k_–2) and the position of the spectral line would characterize this ratio. Specifically, at 330 nm (extremum of the negative band as shown in Fig. 2B)

(8)

where OD₃₃₀(X) is the OD of the test enzyme solution (in our case the respective point in Curve 3 in Fig. 2B), OD₃₃₀(Ec₁) corresponds to the point in the Curve 1, and OD₃₃₀(EGa₁) the corresponding point in Curve 2.

Since in the case shown in Figure 2B this ratio was close to one, k₂ = k_–2 was used in all subsequent modeling. Similar analysis of the TK spectra for the other substrate pair, Sed7P and Rib5P (data not shown) gave the constraint k_–5 = k₅/3. These defined ratios, when used in a regression procedure, decreased the number of unknown parameters of the system.

Considering the complete reaction, for instance, when Xul5P and Rib5P are present as substrates, both under conditions of excess, the specific activity of the enzyme depends only on the catalytic break-down constants, k₂ and k₅. Estimates of these values were obtained as follows: we set the ratio k₅/k_–2 = 1 as an initial value in regression analysis, and used the known values of k₂/k_–2 and k₅/k_–5. Then the experimentally measured rate of production of Gra3P (20 U/mg) gave the value k₂ = k_–2 = k₅ = 5300/min and k_–5 = 1770/min.

The enzyme–substrate dissociation constants were estimated by fitting the time course of product accumulation as it occurred under conditions when it was most sensitive to the parameter(s) in question. The data fitting was performed using the model of TK given in Figure 2A; the system of differential equations which was used is given on our website (‘TK-substrate.nb’). With a saturating concentration of Rib5P, and an initial concentration of Xul5P that was lower than its K_m such that Xul5P binding was the rate-limiting step, the initial rate of the TK reaction depended only on the ratio k_–1/k₁. Figure 2C shows the best fit for two experiments that involved product accumulation observed at the initial concentrations of Xul5P of 25 and 37 µM; the estimated value of k_–1/k₁ was 55 ± 7 µM.

In the above estimation we exploited the fact that the initial rate of the TK reaction, with Xul5P to be the only added ketose, does not depend on the dissociation constant of Sed7P (c₂) because its initial concentration is zero. As the product accumulates, the reaction of Sed7P binding becomes dependent on the value of k₆/k_–6. The simulated time course shown in Figure 2C coincided well with the experimental one at almost any value of the dissociation constant when it exceeded 300 µM. So it was concluded that the actual value was >300 µM.

When the initial concentration of Rib5P was lower than its K_m value and Xul5P was saturating, the kinetics of product accumulation depended only on the dissociation constant of Rib5P. Fitting these experimental data gave an estimate of k_–4/k₄ of 270 ± 30 µM (Fig. 2D).

The value of the ratio k₃/k_–3 was obtained by using the following relationship between the overall rate constant of the metaboliccycle, i.e.

(9)

where the equilibrium constant K_eq = [Rib5P][Xul5P]/([Sed7P][GraP]) = 0.48 (Casassa and Veech, 1986). It gave k₃/k_–3 < 34 µM (taking k₆/k_–6 > 300 µM).

The estimated values leave only five rate constants and the enzyme concentrations to be determined by fitting the isotopomer distribution and restrict a series of TK reactions, because they share common steps, as shown in the Table 1. The same values of unitary rate constants are used in the various different overall reactions. For example, reaction 14 has the same half-reaction Xul5P GraP as the reaction 13; and reaction 15 has the same steps Sed7P Rib5P, as reaction 13. The equilibrium constant of reaction 15, [Rib5P][Fru6P]/([Sed7P][Ery4P]) can be defined as a product of the equilibrium constants of the other two reactions [Rib5P][Xul5P]/([Sed7P][GraP]) and [Fru6P][GraP]/([Ery4P][Xul5P]) and one of the constants could be determinedfrom the equation for reaction 14 or 15 as is done in Equation (9).

Implementing these parameters can be found in the in the Mathematica notebook containing the main program ‘rib-ot.nb’ on our website.

	IMPLEMENTATION AND DISCUSSION

TOP Abstract INTRODUCTION SYSTEMS AND METHODS ALGORITHM IMPLEMENTATION AND DISCUSSION CONCLUSIONS REFERENCES

 
To show the advantages of implementing the kinetic mechanisms of TK and TA reactions in the way that is described here, we start from a reanalysis of the data from Selivanov et al. (2004), which is followed by the analysis of new experimental data. The new software was used to evaluate whether an association of cell proliferation and TK activity is based on changes in metabolic fluxes obtained here as the best-fit estimates. Mass isotopomer data were obtained from the cancer cells HT29 under control conditions and after TK inhibition with OT. The results and experimental conditions are given in Table 1.

Since a good fit to the experimental data was obtained previously, the present analysis was begun with the earlier fluxes. However, implementing the dependency between fluxes showed an inconsistency in the earlier results. In particular, the data revealed that the fluxes in both of the ‘invisible’ TA reactions (Fru6P Gra3P and Sed7P Ery4P) were almost two orders of magnitude greater than the forward and reverse fluxes in the corresponding two-substrate–two-product reactions (Fru6P Sed7P). The scheme shown in the Figure 2A makes it clear that the two ‘invisible’ exchange reactions taken together pass through the same unitary reactions as the full two-substrate–two-product reaction does; therefore such a high rate of ‘invisible’ reactions is inconsistent with such a low rate in the full reaction.

For TK reactions the opposite problem arose; the previously estimated ‘invisible’ fluxes were very low in comparison with the corresponding fluxes in the full two-substrate–two-product reaction; thus no realistic or consistent set of unitary rate constants was found.

Next we show that using kinetic and equilibrium constant constraints in the TK and TA reactions, as is implemented in the present software, filters out the inconsistent solutions.

Restricted by the interdependency of the fluxes via the TK and TA reactions, the program failed to fit the ribose-isotopomer data. It was impossible to increase m₁ and decrease m₂ in Rib5P, as the experimental data demanded, without a simultaneous decrease in m₀. The first numbers in the column denoted ‘Control’ in Table 1 show the fit of the model-simulation to the data. Although the isotopomer profile was similar it did not reproduce well the overall experimental data. Thus, instead of many solutions being compatible with the experimental data, accounting for the kinetic and equilibrium constant restrictions caused the opposite situation; there was no appropriate solution.

However, the computed profile of labeling was qualitatively similar to that obtained experimentally. This similarity with the difference referred mainly to the amount of unlabeled isotopomers, led to the assumption that an initial unlabeled RNA-ribose pool was persisting during the 72 h of incubation. This could be because of reutilization of the initially unlabeled nucleotides in subsequent generations of cells. Assuming this reutilization, we used the program to fit this portion of the unlabeled part and it gave a good fit with ² = 0.139. The best-fit value (Table 1, the column ‘Control’, third values) of 17% non-exchangeable and non-labeled ribose in the final amount of RNA is close to the expected amount contained in the initial cell population. Note that the cells underwent 3 cell cycles during the incubation. Overall, the analysis supported the hypothesis that practically all the nucleotides, once synthesized, are not degraded but are reused in RNA synthesis.

Partial inhibition of TK (upto 25% of that in the controls) involved incubating HT29 cells with OT and resulted in the mass isotopomer distribution given in Table 1 (column ‘OT’). To simulate this experiment, respective changes of TK activity in the model was made; the best-fit of the persistent unlabeled part was in this case 2 times higher than in the control, 29%, which corresponded to the observed decrease of cell accumulation. Moreover, this change simulated the measured profile of isotopomers, which also validated the model. This was in contrast to data fitting that was achieved without accounting for the rate-constant restrictions in the TK and TA reactions.

The fluxes, giving the shown label distributions in both conditions are also present in the Table 1. Despite the fact that the maximal TK activity was decreased for OT-inhibited cells proportionally for all TK-reactions, the profile of changes in real fluxes was not uniform; some fluxes changed proportionally, others remained practically the same or even increased owing to respective changes in the metabolite concentrations. Therefore, complex changes in metabolic fluxes in response to a change in one of the enzymes justify application of a complex analytical tool; it would definitely be impossible to predict the effect of inhibition without such an analysis.

The most natural way to connect the experimentally observed TK-dependent inhibition of cell proliferation with the main TK function in the PPP would be through the Rib5P concentration. R5P is produced in this pathway and is used for RNA and DNA synthesis. However, as is shown in the last row of the Table 1, inhibition of TK leads to an increase, and not a decrease in Rib5P concentration. Thus, if the analysis is correct, it suggests that biochemists should search for other ways in which TK participates in cellproliferation.

Though determined as described above, the K_d values for TK substrates were used simply as initial values and were free to change in the fitting procedure. The fitting procedure left them practically unchanged, giving them the following 99% confidence intervals. However, the fitted values were different from those defined from the analysis of human erythrocyte enzymes (Mulquiney and Kuchel, 2003) indicated in parentheses: Gra3P, 0.0374 ± 0.0136 (0.218); Xyl5P, 0.0495 ± 0.011 (0.176); Rib5P, 0.27 ± 0.027 (0.532), and K_d for Sed7P 0.248 (0.893) mM, obtained using the values of other constants and the equilibrium constant Rib5P*Xyl5P/(Sed7P*Gra3P) = 0.48 (Casassa and Veech, 1986). The algorithm also did not change the K_d values for Fru6P 1.71 ± 0.19 and Ery4P, 0.088 ± 0.00264 taken from study of human erytrocyte enzymes (Mulquiney and Kuchel, 2003). However the fitting procedure changed the starting value for K_d of TA, initially also taken from the analysis of the enzyme of human erytrocytes; Ery4P, initially 0.016, finally 0.288 ± 0.0064; Fru6P from 0.215 to 0.645 ± 0.0172; Gra3P from 0.12 to 0.0672 ± 0.0018; Sed7P from 0.077 to 0.048 ± 0.0019. These changes were critical for reproducing the experimental pattern. One could assume that different isoforms of TK represent majority in normal and transformed cells. This assumption is based on the discovery of new TK-related gene by Coy et al. (1996). who also suggested association between some pathologies and possible isoformsof TK.

	CONCLUSIONS

TOP Abstract INTRODUCTION SYSTEMS AND METHODS ALGORITHM IMPLEMENTATION AND DISCUSSION CONCLUSIONS REFERENCES

 
The software described here constitutes a new tool for the analysis of isotopomer distributions in metabolic systems measured using GC/MS or NMR. Its novelty resides in the integration of ‘classical’ steady-state enzyme kinetic data with isotopomer analysis, which made possible to connect all the isotope-exchange reactions catalyzed by the same enzyme, to determine the K_d values from isotopomer analysis as they apply in vivo, to restrict the possible ambiguity in evaluated fluxes. The example presented here entailed implementing the kinetic and equilibrium constant restrictions to filtering out the incorrect solutions, which fit the particular experimental pattern, but inconsistent with the enzyme kinetic mechanisms. The accepted solution, in addition to precisely reproducing the experimental pattern of isotopomers, proved its reliability by a correct description of the experimental behavior under TK inhibition. The analysis revealed unexpected information about reutilization of nucleotides, which also was consistent for both conditions of cell proliferation activity. After general validation, the model revealed the reliable metabolic profiles, which would be difficult to predict. In addition to flux analysis, it allowed the comparison of the determined in vitro kinetic constants with the best-fit values for in vivo conditions. This could provide understanding of the enzymological conditions in living cells and of the extent to which biochemical data obtained in vitro can be used to interpret the situation in vivo. Thus the earlier biochemical data may be functionally integrated with new information from whole cells and hence be of interest in SystemsBiology.

The present analytic tool can be easily expanded to enable analysis of more complex metabolic systems in living cells. For example, the first step in such an expansion has already been taken; the C++ program for the analysis of glycolysis, the PPP with the addition of substrate competition and the Krebs cycle with its anaplerotic reactions can be found on our website.

	Acknowledgments

 
This work was supported by the grants: Fundation la Caixa (ONO3-70-0), the Ministerio de Ciencia y Tecnología of Spanish Government (SAF2002-02785 and PPQ2003-06602-C04-04); NIH DK56090-A1 to W.N.P.L.; Generalitat de Catalunya (ABM/acs/PIV2002-32) to V.A.S; the Russian Foundation for Basic Research 03-04-49025 to G.A.K. 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. P.W.K.'s contribution was supported by a Discovery Grant for the Australian Research Council.

Conflict of Interest: none declared.

	Footnotes

 
¹In fact, this C++ version is more advanced; it considers the scheme with competition between the three TK reactions and in addition, the Krebs cycle and its anaplerotic reactions. The brief user-guides ‘readme console.doc’ and ‘readme windows.doc’ for this program can also be found on the website in the respective archive files containing executable programs. Only the Mathematica version ‘rib-ot.nb’ with linked ‘rrib.exe’ is described here.

Received on April 22, 2005; revised on June 3, 2005; accepted on July 5, 2005

	REFERENCES

TOP Abstract INTRODUCTION SYSTEMS AND METHODS ALGORITHM IMPLEMENTATION AND DISCUSSION CONCLUSIONS REFERENCES

 

Berthon, H.A., et al. (1993) ¹³C n.m.r. isotopomer and computer-simulation studies of the non-oxidative pentose phosphate pathway of human erythrocytes. Biochem. J., 296, 379–387.

Casassa, J.P. and Veech, R.L. (1986) The interdependence of glycolytic and pentose cycle intermediates in ad libitum fed rats. J. Biol. Chem., 261, 690–698[Abstract/Free Full Text].

Comin-Anduix, B., et al. (2001) The effect of thiamine supplementation on tumor proliferation. Eur. J. Biochem., 268, 4177–4182[ISI][Medline].

Cornish-Bowden, A. Fundamentals of Enzyme Kinetics, (2004) 3rd edn , London Portland Press.

Coy, J.F., et al. (1996) Molecular cloning of tissue-specific transcripts of a transketolase-related gene: implications for the evolution of new vertebrate genes. Genomics, 32, 309–316[CrossRef][Medline].

Heinrich, C.P., et al. (1971) A circular dichroism study of transketolase from baker's yeast. Biochem. Biophys. Res. Commun., 44, 275–279[CrossRef][Medline].

Kochetov, G.A. (1982) Determination of transketolase activity via ferricyanidereduction. Meth. Enzymol, 90, 209–223.

Kochetov, G.A., et al. (1973) The role of the charge transfer complex in the transketolase catalyzed reaction. Biochem. Biophys. Res. Commun, 54, 1619–1626[CrossRef][Medline].

Mulquiney, P.J. and Kuchel, P.W. Modelling Metabolism with Mathematica, (2003) , Boca Raton, FL CRC Press.

Schmidt, K., et al. (1997) Modeling isotopomer distributions in biochemical networks using isotopomer mapping matrices. Biotechnol. Bioeng, 55, 831–840[CrossRef].

Selivanov, V.A., et al. (2004) An optimized algorithm for flux estimation from isotopomer distribution in glucose metabolites. Bioinformatics, 20, 3387–3397[Abstract/Free Full Text].

Usmanov, R.A. and Kochetov, G.A. (1983) Interaction of baker's yeast transketolase modified by 2,3-butanedione with anionic and nonanionic substrates. Biochem. Int., 6, 673–683[Medline].

van Winden, W.A., et al. (2001) Possible pitfalls of flux calculations based on ¹³C-labeling. Metabol. Eng, 3, 151–162.

Wikner, C., et al. J. Biol. Chem, (1994) 269, 32144–32150[Abstract/Free Full Text].

Wiechert, W. and de Graaf, A.A. (1996) In vivo stationary flux analysis by ¹³C labeling experiments. Adv. Biochem. Eng. Biotechnol, 54, 109–154[Medline].

CiteULike    Connotea    Del.icio.us    What's this?

This article has been cited by other articles:

				R. J. Gillies, I. Robey, and R. A. Gatenby Causes and Consequences of Increased Glucose Metabolism of Cancers J. Nucl. Med., June 1, 2008; 49(Suppl_2): 24S - 42S. [Abstract] [Full Text] [PDF]

				V. A. Selivanov, S. Krause, J. Roca, and M. Cascante Modeling of Spatial Metabolite Distributions in the Cardiac Sarcomere Biophys. J., May 15, 2007; 92(10): 3492 - 3500. [Abstract] [Full Text] [PDF]

				V. A. Selivanov, S. Marin, P. W. N. Lee, and M. Cascante Software for dynamic analysis of tracer-based metabolomic data: estimation of metabolic fluxes and their statistical analysis Bioinformatics, November 15, 2006; 22(22): 2806 - 2812. [Abstract] [Full Text] [PDF]

This Article Abstract FREE Full Text (Print PDF) All Versions of this Article: 21/17/3558    most recent bti573v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Search for citing articles in: ISI Web of Science (4) Request Permissions Google Scholar Articles by Selivanov, V. A. Articles by Cascante, M. Search for Related Content PubMed PubMed Citation Articles by Selivanov, V. A. Articles by Cascante, M. Social Bookmarking          What's this?

Online ISSN 1460-2059 - Print ISSN 1367-4803

Copyright © 2008 Oxford University Press

Oxford Journals Oxford University Press

• Site Map
• Privacy Policy
• Frequently Asked Questions

Other Oxford University Press sites: Oxford University Press Oxford Journals Japan American National Biography Booksellers' Information Service Children's Fiction and Poetry Children's Reference Corporate & Special Sales Dictionaries Dictionary of National Biography Digital Reference English Language Teaching Higher Education Textbooks Humanities International Education Unit Law Medicine Music Online Products Oxford English Dictionary Reference Rights and Permissions Science School Books Social Sciences Very Short Introductions World's Classics