- Open Access
Integration of enzyme kinetic models and isotopomer distribution analysis for studies of in situcell operation
BMC Neuroscience volume 7, Article number: S7 (2006)
A current trend in neuroscience research is the use of stable isotope tracers in order to address metabolic processes in vivo. The tracers produce a huge number of metabolite forms that differ according to the number and position of labeled isotopes in the carbon skeleton (isotopomers) and such a large variety makes the analysis of isotopomer data highly complex. On the other hand, this multiplicity of forms does provide sufficient information to address cell operation in vivo. By the end of last millennium, a number of tools have been developed for estimation of metabolic flux profile from any possible isotopomer distribution data. However, although well elaborated, these tools were limited to steady state analysis, and the obtained set of fluxes remained disconnected from their biochemical context. In this review we focus on a new numerical analytical approach that integrates kinetic and metabolic flux analysis. The related computational algorithm estimates the dynamic flux based on the time-dependent distribution of all possible isotopomers of metabolic pathway intermediates that are generated from a labeled substrate. The new algorithm connects specific tracer data with enzyme kinetic characteristics, thereby extending the amount of data available for analysis: it uses enzyme kinetic data to estimate the flux profile, and vice versa, for the kinetic analysis it uses in vivo tracer data to reveal the biochemical basis of the estimated metabolic fluxes.
Introduction: application of stable isotope tracer data in neuro -biology and -medicine
Metabolic networks of living cells produce the intricate redistribution of carbon skeleton atoms of substrates. If these substrates are artificially labeled by stable isotopes (such as 13C) at specific positions, the reorganization of carbon skeleton becomes measurable and its quantification provides insight to the respective metabolic reactions. Figure 1 shows an example of isotope exchange performed by one of the reaction catalyzed by transketolase (TK) in the non-oxidative pentose phosphate pathway (PPP). Interconnection of several isotope exchange reactions creates in each metabolite a variety of forms, which differ by the number and positions of 13C isotopes (13C isotopomers). A given set of metabolic fluxes produces a specific distribution of isotopomer fractions, and consequently, the isotopomer distribution indicates the underlying set of fluxes.
Two techniques for isotopomer detection were used to estimate metabolic fluxes in situ: nuclear magnetic resonance (NMR) and mass spectrometry (MS). The different kinds of data that can be obtained from NMR measurements were classified by Mollney et al . One-dimensional 1H NMR measurement can provide positional 13C enrichment, i.e. fraction of molecules with label in specific carbon position in a molecule [2, 3]; multiple peaks of one dimensional 13C NMR reveal groups of isotopomers . Two-dimensional 1H-13C NMR essentially allows quantification of individual isotopomer fractions [5, 6]. MS provides information on the fractions of different mass isotopomers [7, 8], i.e. groups of isotopomers with a fixed number of labels regardless their positions. It detects fewer isotopomers than NMR but is usually more sensitive.
Stable isotope detection is now used for the metabolic flux profile estimation in various areas of cell and tissue biology, and neuroscience is not an exception. Implementations of the in vivo NMR measurements to study brain function were described in the recent review article . Notably, 13C NMR studies recently showed the significance of glutamate release and recycling between neurons and glia. The neurotransmitter pool, which previously was considered to be small and metabolically inactive, appeared to be included in glutamine-glutamate recycling, the major neuronal metabolic pathway. The activity associated with glutamate neurotransmission was found to be linearly dependent upon glucose oxidation, and this supports the molecular model of stoichiometric coupling between glutamate neurotransmission and functional glucose oxidation. Thus, the considered here models of central glucose metabolism are clearly of great interest to the field of neuroscience.
Since 13C is known to be harmless, it has been used in human subjects , particularly for mapping the spatial localization of metabolites using NMR spectroscopic imaging [3, 11]. This method, designed for clinical use, has been applied as a diagnostic tool in pediatric and adult brain disorders [12, 13]. A non-invasive 13C magnetic resonance spectroscopy technique allowed the synthesis rate of N-acetyl-1-aspartate to be determined in patients with Canavan disease [14, 15]. It is also widely used in neuroscience to study brain metabolism in animals [16, 17]. Bixel et al  revealed an interesting application of this method by using it to define all products of leucine metabolism in cultured astroglial cells.
While most of the results described above were obtained by the qualitative comparison of NMR measurements, their quantitative analysis would offer much more profound insights into cell processes. A tool able to estimate metabolic fluxes would reveal the dynamic characteristics of cellular phenotype under specific conditions, thus complementing genomic and proteomic methods, which only reveal static characteristics of the cell. The specific tools for quantitative metabolic flux profile analysis from measured isotopomer distribution data are described next.
Current status of isotopomer distribution analysis
Very approximate estimation of metabolic fluxes could be provided by implementation even formulas derived using simplifications such as an assumption that the fluxes are unidirectional [19–22]. Such formulas were applied to the analysis of enrichment in specific carbon position or mass isotopomers. Although such analysis is easy to perform, it is incomplete, and as we have found (unpublished data), the simplifications could lead to unreliable estimates. More reliable tools that supplement data acquisition normally include the following: (i) a mathematical model, which simulates the distribution of isotopomers, (ii) algorithms used to fit respective experimental pattern and thus evaluate the matching flux profile, and (iii) statistical analysis of the best fit parameters.
Different levels of complexity could be implemented in mathematical models for the isotopomer distribution analysis. Katz and Wood  were the first who estimated the metabolic fluxes by solving the mass balance equations for the positional enrichment of metabolites. However, this approach was only appropriate for a small number of fluxes that, in addition, were unidirectional. Consideration of real bidirectional fluxes increased a number of unknown parameters rendering the systems of equations underdetermined. This disadvantage demanded more advanced analytical applications, which would be able to analyze more experimental data in order to unambiguously determine the parameters. Mason and Rothman  described the basic principles of simple model construction for the analysis of some selected positional enrichment together with executing a "classical" kinetic model that evaluated the time course of metabolite concentrations using the available kinetic information. The idea to complement kinetic models with isotopomer distribution data and, in the same time, to link isotopomer analysis with kinetic data increased the reliability and informative outcome of both approaches. However, the described method for isotopomer analysis cannot account for all the information provided by NMR experiments. For any given molecule containing n carbon atoms, it only filters a maximum of n enrichments at selected carbon positions, while the total number of all isotopomers is 2n and all of them could be important for the estimation of the metabolic flux profile. Thus a part of information contained in the NMR data is not used for the analysis. In fact, the NMR measurements, as it is mentioned in the introduction, could provide the relative concentrations of individual isotopomers or at least the sum of those that contain a specific label pattern.
The comprehensive analysis of mass isotopomer data obtained by MS also demands that all possible individual isotopomers are computed, because every mass isotopomer is the sum of several individual isotopomer concentrations and each component of the sum could be produced in the different way.
According to the total number of isotopomers, the calculation of their fractions requires solving 2n equations for the substance consisted of n carbons. In this way even for the model description of all isotopomers formation in glycolysis and PPP the algorithm must construct (it is practically impossible to write such a huge number of equations by hand) and solve several hundred equations.
Schmidt et al.  developed a simulation algorithm for the computation of all possible isotopomers. This algorithm of isotopomer mapping matrices, evolved from the earlier approach of atom mapping matrices [2, 26], allowed to construct and solve automatically hundreds of equations that describe all isotopomers transformation according to the specific reaction mechanisms. The iterative procedure, which was used for the solution, in fact provided the evaluation of time course for the relaxation of each isotopomer fraction to steady state. This method made it possible to simulate all isotopomers distribution for the given set of metabolic fluxes, while the best fit to the experimental distribution pointed to the set compatible with the measured distribution. However, the presence of large exchange fluxes caused severe instability of numerical solution or convergence problems , which restricted the application of this method.
Wiechert et al  found an elegant way to overcome the problem of instability by reformulating isotopomer balance equations into cumomer balance equations. The term cumomer (from cumulative isotopomers) designates a sum of isotopomers with fixed positions of label (a 1-cumomer fraction is a sum of the positional enrichments). Such reformulation of the equations in terms of cumomers simplified the task and made it possible to obtain the solution in one step based on matrix calculus. The equations formulated in terms of cumomer fractions can be solved explicitly as a cascade of linear systems, evaluating the cumomer fractions one by one starting from the 0-cumomer fraction. Then the cumomer fractions can be transformed back into isotopomer fractions.
This method was well elaborated and all the history of its development from analysis of positional enrichment experiments was well documented [27–30]. A computer program based on the cumomer-balance method performs the complete isotopomer analysis and finds the flux profile by experimental isotopomer distribution fitting . It is easy to work with this program; user just needs to write in symbolic form the reactions designed for analysis marking the respective re-formations of carbon skeleton. Using this input, the program automatically constructs isotopomer balances, transforms to cumomer analysis, solves the equations, and displays the results in the desired form. Creation of such a tool seemed to solve eventually the problems of isotopomer analysis.
However, there are at least three reasons to develop one more approach to isotopomer analysis. First, although the above described algorithm overcomes the problem of instability of iterative numerical solution, it is completely restricted by stationary flux analysis thus leaving without any examination the available time course of label distribution. In fact, the stability of solution could be controlled without losing the advantage of testing the time course of isotopomer accumulation, which could be more informative than steady state analysis.
Second, all the above approaches to the computation of isotopomer transformations consider fluxes as independent variables. This was noted as an advantage of the method because such an analysis did not need any assumptions regarding the biochemical basics of considered fluxes . However, this could be an advantage only if the experimental pattern of label distribution is sufficient for evaluation of a unique true set of fluxes. If the amount of information is insufficient, additional data are necessary for the unambiguous evaluation of flux profile. In this case the known kinetic characteristics of analyzed enzyme reactions and results of classical kinetic modeling could provide such additional necessary information.
Third, even if the analysis reveals the fluxes taking them as independent, the fluxes remain disconnected from the detailed mechanisms of the catalysis and regulation considered in kinetics models, thus the biochemical reasons for the observed behavior remain unclear. In this case, the use of kinetic modeling could also solve this problem. Long era of classical biochemistry developed a number of kinetic models of complex systems that use known characteristics of enzyme catalysis and regulation. These parameters could be employed in metabolic flux analysis, providing the necessary additional information. An excellent example is a model of erythrocyte central metabolism , which includes kinetic models of all involved enzymes with known regulatory and catalytic mechanisms and kinetic constants, verified by numerous experiments.
On the other hand, the kinetic models of complex systems, analyzing the experimentally observed cellular functions as a result of operation of many regulated processes, include many parameters and therefore, like the flux analysis, also suffer from insufficiency of experimental data. Moreover, the kinetic models normally include characteristics of enzymatic reactions obtained in vitro, and such data cannot always be used for the interpretation of in vivo experiments. The use of in vivo tracer data would animate the old classical district of kinetic studies.
In this situation integration of kinetic modeling with complete isotopomer analysis would provide:
for kinetic study of in vivo cell operation the new area of tracer data, which are necessary for understanding the organization and regulation of the processes in living cells and applicability of classical in vitro information;
for metabolic flux analysis the additional information to restrict the number of acceptable sets of metabolic fluxes by the ones that are compatible not only with a given pattern of isotopomer distribution, but also with the data of previous biochemical studies.
The way of such integration is described next.
Algorithms for integrated kinetic and metabolic flux analysis
As MFA is a commonly accepted acronym for metabolic flux analysis based on isotopomer data, the tool proposed by Selivanov et al [32, 33] designed for integrated kinetic and metabolic flux analysis is abbreviated here as IKMFA.
IKMFA was designed to possess the following characteristics:
To be compatible with any kinetic model of central metabolism, so that the MFA part accepts the total fluxes and metabolite concentrations predetermined by the kinetic model constituting the first part of the analytic software.
To use the total fluxes and metabolite concentrations, obtained from the kinetic model, for estimation of the time courses and distribution of all possible isotopomers.
To use fitting of the experimental time course, or/and steady state isotopomer distribution or/and global metabolite concentrations for the estimation of both fluxes and parameters of the kinetic model used.
These characteristics render IKMFA able to implement the detailed tissue-specific kinetic mechanisms of the enzyme reactions to describe the fluxes. In other words this approach builds a complete isotopomer analysis on the top of a kinetic model. Used for the analysis of isotopomer distribution data, the kinetic mechanisms, validated by all available kinetic models, could be examined for the in vivo conditions. This links the estimated metabolic fluxes with molecular and kinetic mechanisms of the metabolic pathways, i.e. with their biochemical basis. There are no limits in respect to the expressions for reaction rates, catalysis could be described according to the kinetic mechanism of any levels of complexity as it is accepted in kinetic modeling. An application of IKMFA for analysis of complex metabolic system that connects kinetic models of central metabolism with MFA has been considered elsewhere . This metabolic system is schematically represented in the Figure 2.
As indicated above, the first step of analysis is the solution of ordinary differential equations (ODEs), which describe the total concentration change as the sum of the production rate of the given metabolite minus the rates of its consumption.
In principle any ODE solver could serve in the kinetic part of the analysis. Since the software was written in programming language "C++", a compatible ODE solver was chosen. The existed source "C++" libraries, such as "Numerical recipes in C++", http://www.nr.com/ or the C++ class library 'ODE++' of Milde (2003), http://www.minet.uni-jena.de/www/fakultaet/iam/ode++/ode++_e.html were appropriate solvers. The metabolite concentration obtained from the ODEs solution become initial values of non-labeled isotopomers of internal metabolites in the following step of analysis. Other isotopomers are initially set to 0, except the outside metabolites, which have the initial distribution according to that added experimentally. Fluxes, also obtained from the solution of ODEs, are used in the second part to simulate the respective reactions between the isotopomers (Figure 2).
In spite of evident similarity between the first step of IKMFA and ordinary kinetic model execution, there is an essential difference between them. Kinetic models normally relate one enzyme reaction to the net flux as the difference between the forward and reverse fluxes, because only net fluxes define changes in the calculated total metabolite concentrations. In contrast, the first step of IKMFA needs to compute the forward and reverse fluxes separately and also some additional fluxes, which also serve as an input for the subsequent MFA. Definitions of such additional fluxes and a way to compute them are described next.
Definitions and algorithm for evaluation of all isotope-exchange fluxes
A classical example of a variety of fluxes is the non-oxidative PPP, the most problematic metabolic part related to numerous isotope-exchange reactions catalysed by TK and TA [34, 35]. TK, thiamine diphosphate (ThDP)-dependent enzyme, 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, forming new ketose- and aldose-phosphates. TA catalysis operates via another mechanism involving a Schiff base formed directly between the enzyme and the substrate, and, as a result, a three-carbon fragment (dihydroxyacetone) is transferred from a ketose-phosphate donor to an aldose-phosphate acceptor. However, a similar scheme could represent the steps of TA catalysis and TK catalysis is therefore considered here as representative of both enzymes. Although only two TK reactions are widely accepted to take place in PPP,
r5p + x5p <-> s7p + g3p (1)
e4p + x5p <-> f6p + g3p, (2)
this enzyme catalyzes in fact much more isotope-exchange reactions [34, 35]. Using as an example the first of two above reactions (also shown in Figure 1), Figure 3 illustrates all possible isotope exchange fluxes related with this TK reaction.
The reversible steps associated with the ping-pong mechanism of TK reaction involve exchange (i) between the ketose substrate and product of its cleavage, or, in the present example, between xylulose-5-phosphate (xu5p) and glyceraldehyde-3-phosphate (g3p) shown in Figure 3A and also between sedoheptulose-7-phosphate (s7p) and ribose-5-phosphate (r5p) shown in Figure 3B, and (ii) between two ketoses (xu5p and s7p, Figure 3C). Thus, six isotope exchange fluxes are associated with one TK reaction. As shown previously  all of them are essential for distribution of the label and should be taken into account for correct estimation of metabolic fluxes. These fluxes could be expressed through the rate constants of elementary steps as the legend to Figure 3 describes in details.
In fact, forward and reverse transfer of the first two carbons between the x5p and s7p are coupled with transfer of the three last carbons of x5p to g3p pool and back, and also the five last carbons of s7p to r5p pool and back. Therefore, transfer between ketose and product of its cleavage could be presented as consisted of two parts, one of which is coupled with transfer between two ketoses and some another, additional part. The legend to Figure 3 explains in details that these additional parts are equal in forward and reverse directions. Taking this into account, one should compute (in the way shown in the legend to Figure 3) four different isotope transfers related with just one TK reaction (1):
xu5p->s7p accompanied by the transfer xu5p->g3p and r5p->s7p
s7p->xu5p accompanied by the transfer s7p->r5p and g3p->xu5p
additional transfer xu5p<->g3p the same in both directions
additional transfer s7p<->r5p the same in both directions.
Presence of different substrates in the intracellular volume further complicates the situation. Figure 4 shows three reactions catalyzed by TK in vivo, which compete for the same enzyme. In this situation the number of isotope exchange fluxes related to the TK reactions increases, but all of them can be expressed through the elementary rate constants in the way similar to that indicated in the Figure 3 legend. The following isotope-exchange fluxes related to the TK reactions shown in the Figure 4 should be accounted for subsequent analysis in the way similar to that described in the legend of Figure 3 for one separated reaction:
xu5p->s7p accompanied by the transfer xu5p->g3p and r5p->s7p;
s7p->xu5p accompanied by the transfer s7p->r5p and g3p->xu5p;
xu5p->f6p accompanied by the transfer xu5p->g3p and e4p->f6p;
f6p->xu5p accompanied by the transfer f6p -> e4p and g3p -> xu5p;
s7p->f6p accompanied by the transfer s7p->r5p and e4p->f6p;
f6p->s7p accompanied by the transfer f6p->e4p and r5p->s7p;
additional transfer xu5p<->g3p the same in both directions;
additional transfer s7p<->r5p the same in both directions;
additional transfer e4p<->f6p the same in both directions;
Expressed through the elementary rate constants such isotope exchange fluxes are evaluated for all the enzymes in the course of executing the kinetic model. The obtained values are used in the second step of analysis, namely in simulation of isotopomer distribution, which is done in the way similar to that accepted in MFA (e.g. in ). The respective algorithm is described next.
Reactions between isotopomers
To interpret the result of isotope redistribution in metabolites of a metabolic pathway all the possible reactions between isotopomers are simulated at this step of analysis. An algorithm of simulation of isotopomer distribution, which is present below, is in principle similar to that described elsewhere , though it is implemented by different means. In contrast to the earlier algorithms , here the total fluxes are taken from the previous step of analysis, where they are expressed as functions of total concentrations of metabolites and effectors according to the reaction mechanisms. Moreover, the program calculates the real concentrations and time courses of isotopomers and then defines fractions corresponding to the analyzed experimental data.
For isotopomer designation we use a binary notation for the 13C and 12C atoms as it is helpful in optimization of references to a specific isotopomer. Since each carbon atom of a molecule can exist in one of the two states: labeled (marked as '1') or unlabeled ('0'), each metabolite in the model can be represented by an array of 2n of possible isotopomers, where n is the number of carbon atoms in the molecule. Each isotopomer in the model is represented as binary numbers; its digits correspond to the carbon atoms in a molecule (3 digits for trioses, 4 for erythrose, etc.). A '1' or a '0' in certain position in a string signifies that corresponding carbon atom is labeled or unlabeled. For instance, all isotopomers for glyceraldehyde-3-phosphate (g3p) are:
000, 001, 010, 01l, 100, 101, 110, 111
This representation of isotopomers as successive integer numbers is very convenient because the model uses it as references to the respective position in the existing array of concentrations of all isotopomers, different for each metabolite. This ordering of the isotopomers in the array allows optimization of referring the isotopomer products for any isotopomer substrate related to the considered reactions, as it is explained below. For instance, the reaction (xu5p + r5p ↔ s7p + g3p) between arbitrarily chosen isotopomers of x5p (10011) and r5p (10000) produces the following isotopomers of s7p and g3p:
10011(xu5p) + 10000(r5p) ↔ 1010000(s7p) + 011(g3p) (3)
This reaction could be easily simulated by manipulations with binary numbers, viz., the product of scission of the ketose substrate of n carbons will be a product of bitwise AND of the respective binary number with 2n-2-1, in the above example (10011 AND 111) = 011. To define the reference number of the produced new ketose, the reference number of initial ketose should be shifted n-2 positions right, and then (if m is number of carbons in the aldose substrate) m positions left. The bitwise OR operation performed for the obtained number and reference number of aldose substrate produces the resulted ketose isotopomer. In the above example, right shift of 10011 gives 10, its left shift gives 1000000 and (1000000 OR 10000) gives 1010000. These operations are fast and allow the reference number of isotopomer products to be defined for each pair of isotopomer-substrates.
Determination of the reference numbers of products for each pair of substrates is the basis of the optimization. Then the change of concentrations according to the succession of reaction is calculated. Assuming that all isotopomers have the equal affinity, the rate of reaction between pair of isotopomers is proportional to their concentrations, while the sum of reaction rates for all isotopomers gives the global metabolic flux calculated also at the previous, kinetic, step. If, for instance, Vf is the global forward flux for the TK reaction between X5P and R5P, the flux for this reaction between isotopomers i and j (V3fij) would be expressed as follows:
Vfij = Vf × [X5Pi] × [R5Pj]/([X5Ptot] × [R5Ptot]) (4)
Here, the indices i and j refer to the concentrations of the respective isotopomers and the index 'tot' refers to the total concentration of the metabolite, as calculated in the kinetic model. In this way, the first part, ODE solving, is linked to the second part that computes the label distribution: the kinetic model calculates global fluxes and concentrations and defines the values (as in the above example V3f/([X5Ptot] x [R5Ptot])), which are used to get the fluxes in reactions between isotopomers (as in Equation (4)). During the small time interval dt the TK reaction (3) considered as an example consumes the amount dt × V3fij of each of the isotopomers i and j and this value is subtracted from the concentrations of respective isotopomers.
[X5P i ]t+dt= [X5P i ] t - dt × Vf ij
[R5Pj]t+dt= [R5P j ] t - dt × Vf ij (5)
Here, the indices t and t + dt indicate the time of the simulated process. As usual for numerical solution of differential equations, this could be acceptable approximation if the time step dt is so small that consumption during one step is small compared to the amount of these isotopomers; practically, the value for dt is taken so that its further decrease does not affect the solution. The program adds the same amount to the concentrations of reaction products, which numbers (ra and rd) are defined from i and j as described above.
[GAP ra ]t+dt= [GAP ra ] t + dt × Vf ij
[S7P rd ]t+dt= [S7P rd ] t + dt × Vf ij (6)
The above algorithm is present in order to illustrate the main principles of TK reaction simulation between a pair of isotopomers according to the Equations 5 and 6. However this algorithm is not optimized in the best way. To simulate the reaction between n isotopomers of ketose and m isotopomers of aldose it is necessary to perform n × m cycles of such calculations. It could be optimized by the imaginary separation of the whole reaction into two steps corresponding to the aldose and ketose products formation. First scission of all ketose isotopomers could be simulated, producing respective isotopomers of aldose product according to the scheme x5p->g3p; this demands just n recalculations of x5p and g3p concentrations. Then the rest two-carbon fragment (for TK) with only 4 different possible combinations of label interacts with isotopomers of aldose substrate (C2+r5p->s7p). This demands 4 × m recalculations of r5p and s7p for one direction of the reaction. In this way the whole description of the TK reaction needs only 2×(n+4×m) recalculations of substrate-product pairs instead of n × m recalculations of four substances according to the above algorithm.
Analysis of experimental data starts from execution of the kinetic model simulating time course of metabolite concentrations and fluxes, which then are used in the second step of simulation of corresponding labeled isotopomer distribution. Experimental data fitting finds the flux profile and kinetic parameters compatible with the analyzed data. This is described next.
A combination of kinetic modeling with isotopomer distribution analysis allows combining also the respective experimental data, which can be analysed. Moreover, IKMFA expands the usual steady state isotopomer analysis to the non-steady state conditions. The following kinds of experimental data could be subjects of fitting:
– Measured rate of production of various metabolites in cell cultures under different conditions of incubation and intracellular concentration of some metabolites;
– Total metabolite concentrations;
– Distribution of labeled atoms such as 13C isotopes in metabolites (13C isotopomers), when the label was added with some of the substrates.
Optimization of a merit function in a multidimensional space of parameters is well elaborated and usually does not represent a problem [36, 37]. However, in our particular case, two kinds of computational problems arise. First, if the merit function (χ2, the sum of the squares of the vertical distances of the experimental points from the calculated curve divided by the standard deviation) is defined from the numerical solution of differential equations, abrupt change of parameters demanded by the applied fitting algorithm could render the system stiff and induce failure of ODE solver. Second, if the total metabolite concentrations are not measured, the merit function does not account for them, so that the best fit sometimes corresponds to evidently unreal concentrations.
The first problem was partly solved by our approach of two-step solution, when the isotopomer analysis comes after the differential equations that describe changes in total metabolite concentrations are solved. In the first step a small number of equations for total concentrations are to be solved and at this step the solution is usually robust. A number of numerical methods for ODE solution are implemented, in particular the Bulirsch-Stoer method  could be recommended. If the problem of stiffness nevertheless arises it could be solved using implicit or semi-implicit methods e.g. using the Bader and Deuflhard discretization of Bulirsch-Stoer method . The problem of stiffness often arises in the second step when much more equations for isotopomers have to be solved. The methods for stiff ODEs will be implemented for the solution of the huge systems for isotopomers, however, even the use of the steady state metabolite concentrations, known from overall the solution, as initial values for unlabeled isotopomers improves the robustness of the solution.
Another part of the stiffness problem in fact was solved when, to address the second problem (of unreal best-fit concentrations), we introduced the threshold for χ2 change, dependent on metabolite concentrations. According to the used Powel's algorithm for the local decrease of χ2 value, the program changed parameters one by one until the minimum was reached in each direction. To avoid long descent in almost flat valley, if a change in χ2 was lower than threshold, the program switched to another direction in the space of parameters. If some metabolite concentration has reached its critical value, the program automatically increased the threshold for χ2 change if the parameter change increased the critical concentration, but decreased the threshold to 0 if the parameter change induced also a decrease of critical concentration. These rules let to control the range of concentrations, and made the solution of ODEs more stable. After termination of descent in the space of parameters, the program make a step uphill by random change of a parameter as is supposed by Simulated Annealing algorithm, and then repeats the downhill descent. All the successful steps of parameter change are saved, so that if the program comes into the area of stiffness and the solver fails, the procedure could be started again from the last successful step after the necessary correction.
After the program finds the minimum of χ2 and completes the fitting procedure, it can determine the set of parameters most relevant to the fit using singular value decomposition of the second derivative matrix of χ2 with respect the parameters (Hessian matrix). Square roots of the diagonal elements of the inverse of Hessian matrix provide the standard deviation for the essential parameters. The present in our website instruction describes the way of using this feature of the program.
Thus, fitting algorithm accepts different kinds of data and estimates not only compatible set of metabolic fluxes but also parameters of enzyme reactions and their regulation, providing insight to the biochemical underground of the particular set of metabolic fluxes.
Prime results and perspective
The advantages of implementation of comprehensive enzyme kinetic mechanisms (e.g. as explored in the Figure 3 and shown in Figure 4) in metabolic flux analysis were examined in , where the two levels of profundity in flux analysis were compared. When the fluxes created in TK reaction were considered as independent variables, as was commonly accepted currently in metabolic flux analysis, they were compatible with the experimental r5p and lactate label distribution in cancer cell line HT29. However, the set of isotope-exchange fluxes obtained in this way was not unique and was found to be incompatible with the TK reaction mechanism. The implementation of dependency between all TK fluxes through the unitary rates, as was accepted in classical enzyme kinetic analysis, restricted the area of possible estimates thus making the analysis more reliable.
On the other hand, such integrated analysis brought a new area of isotopomer distribution data to the kinetic study of enzyme operation in vivo. Fitting isotopomer distribution directly accessing the kinetic constants provided insight to the enzyme kinetic mechanisms that operate in vivo. Some of the kinetic constants obtained from in vitro study were corrected based on the analysis of in vivo isotopomer distribution. The ability to accumulate the information from different sources and to check its compatibility makes this approach scalable into analytical tool which filters and orders available data, as well as evaluates metabolic fluxes compatible with maximum of available information.
The analysis of tracer data considered here focuses mainly on central metabolic pathways that, according to the recent review , are stoichiometrically coupled with glutamate neurotransmission. In this view a major limitation in the interpretation of tracer data is the poor definition of relationship between neuronal activity and neuroenergetic processes supported by glucose metabolism [38–40]. The term neuronal activity applies to a spectrum of energy requiring processes including action potential propagation, neurotransmitter release and uptake, vesicular recycling, and maintenance of membrane potentials . All of these processes are involved in short-term neuronal information encoding and the relative distribution of energy among them remains to be determined. The contribution of the different classes of neurons in a cortical region to the overall energy consumption remains also unclear. To address these questions the flux analysis based on tracer data should be coupled with the comprehensive models of central metabolism interactions with neuronal activity; this assumes further development of IKMFA towards higher complexity of the kinetic model underlying the flux analysis. Such a model could include different cellular compartments or even consist of partial differential equations describing the diffusion of components in intracellular volume. The requirement for multi-compartment analysis of label distribution data follows, in part, from a series of experiments where even qualitative analysis of metabolite labeling pointed out to the existence of different pools of pyruvate in neurons and astrocytes [42–44]. The existence of non-mixed pools of soluble substance in cytosole currently seems surprising. Using a multi-compartment model coupled with isotopomer analysis for the analysis of such data would allow other possible hypotheses to be excluded and would provide quantitative information regarding the relationship between the pools.
Although the data on the existence of different intracellular compartments not separated by membrane structures are not widely accepted, they appear in various areas of cell biology and deserve thorough analysis. Analysis of the behavior of ATP-sensitive K+ channel, located in sarcolemma revealed large differences in ATP levels at different distance from the sarcolemma [45, 46]. Similarly high diffusion limitations in the membrane vicinity were found for cAMP, as follows from the study of cAMP gated channels . Understanding real cellular processes requires integrated analysis of various data relating to the real multi-compartment intracellular space with restricted diffusion limitation for metabolites. The approach described here has an ability to develop into a tool that can be used in such an integrative analysis. An understanding of regulatory mechanisms as an addition to the reliable estimations of metabolic fluxes will increase outcome of stable isotope tracer analysis, which is already in clinical use.
pentose phosphate pathway
nuclear magnetic resonance
metabolic flux analysis.
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: 86-103. 10.1002/(SICI)1097-0290(1999)66:2<86::AID-BIT2>3.0.CO;2-A.
Marx A, de Graaf AA, Wiechert W, Eggeling L, Sahm H: Determination of the fluxes in the central metabolism of Corynebacterium glutamicum by nuclear magnetic resonance spectroscopy combined with metabolic balancing. Biotechnol Bioenerg. 1996, 49: 111-129. 10.1002/(SICI)1097-0290(19960120)49:2<111::AID-BIT1>3.0.CO;2-T.
Gruetter R, Novotny EJ, Boulware SD, Mason GF, Rothman DL, Shulman GI, Prichard JW, Shulman RG: Localized 13C NMR spectroscopy in the human brain of amino acid labeling fromD-[l-13C]glucose. J Neurochem. 1994, 63: 1377-1385.
Jeffrey FM, Rajagopal A, Malloy CR, Sherry AD: 13C-NMR: a simple yet comprehensive method for analysis of intermediary metabolism. Trends Biochem Sci. 1991, 16: 5-10. 10.1016/0968-0004(91)90004-F.
Szyperski T: Biosynthetically directed fractional 13C-labeling of proteinogenic amino acids – an efficient analytical tool to investigate intermediary metabolism. Eur J Biochem. 1995, 232: 433-448. 10.1111/j.1432-1033.1995.tb20829.x.
Henry PG, Oz G, Provencher S, Gruetter R: Toward dynamic isotopomer analysis in the rat brain in vivo: automatic quantitation of 13C NMR spectra using LCModel. NMR Biomed. 2003, 16: 400-412. 10.1002/nbm.840.
Lee W-NP, Byerley LO, Bergner EA, Edmond J: Mass isotopomer analysis: theoretical and practical considerations. BiolMass Spect. 1991, 20: 451-458. 10.1002/bms.1200200804.
Fischer E, Zamboni N, Sauer U: High-throughput metabolic flux analysis based on gas chromatography-mass spectrometry derived 13C constraints. Analytical Biochemistry. 2004, 325: 308-316. 10.1016/j.ab.2003.10.036.
Rothman DL, Behar KL, Hyder F, Shulman RG: In vivo NMR studies of the glutamate neurotransmitter flux and neuroenergetics: implications for brain function. Annu Rev Physiol. 2003, 65: 401-427. 10.1146/annurev.physiol.65.092101.142131.
Rothman DL, Novotny EJ, Shulman GI, Howseman AM, Petroff OA, Mason G, Nixon T, Hanstock CC, Prichard JW, Shulman RG: Prichard and R.G. Shulman, 1H-[13C] NMR measurements of [4-13C]glutamate turnover in human brain. Proc Natl Acad Sci USA. 1992, 89: 9603-9606. 10.1073/pnas.89.20.9603.
Beckmann N, Turkalj I, Seelig J, Keller U: 13C NMR for the assessment of human brain glucose metabolism in vivo. Biochemistry. 1991, 30: 6362-6366. 10.1021/bi00240a002.
Bluml S, Moreno A, Hwang JH, Ross BD: 1-(13)C glucose magnetic resonance spectroscopy of pediatric and adult brain disorders. NMR Biomed. 2001, 14: 19-32. 10.1002/nbm.679.
Bluml S, Moreno-Torres A, Ross BD: [1-13C]glucose MRS in chronic hepatic encephalopathy in man. Magn Reson Med. 2001, 45: 981-993. 10.1002/mrm.1131.
Moreno A, Bluml S, Hwang JH, Ross BD: Alternative 1-(13)C glucose infusion protocols for clinical (13)C MRS examinations of the brain. Magn Reson Med. 2001, 46: 39-48. 10.1002/mrm.1158.
Moreno A, Ross BD, Bluml S: Direct determination of the N -acetyl-L-aspartate synthesis rate in the human brain by (13)C MRS and [l-(13)C]glucose infusion. J Neurochem. 2001, 77: 347-350.
Fitzpatrick SM, Hetherington HP, Behar KL, Shulman RG: The flux from glucose to glutamate in the rat brain in vivo as determined by 1H-observed/13C-edited NMR spectroscopy. J Cereb Blood Flow Metab. 1990, 10: 170-179.
Künnecke B, Cerdán B, and Seelig J: Cerebal metabolism of [l,2-13C2]glucose and [U- 13C4]3-hydroxybutyrate in rat brain as detected by 13C NMR spectroscopy. NMR Biomed. 1993, 6: 264-277.
Bixel MG, Engelmann J, Willker W, Hamprecht B, Leibfritz D: Metabolism of [U- 13C]Leucine in Cultured Astroglial Cells. Neurochem Res. 2004, 29: 2057-2067. 10.1007/s11064-004-6879-8.
Jones JG, Solomon MA, Cole SM, Sherry AD, Malloy CR: An integrated (2)H and (13)C NMR study of gluconeogenesis and TCA cycle flux in humans. Am J Physiol Endocrinol Metab. 2001, 281: E848-E856.
Van Dijk TH, van der Sluijs FH, Wiegman CH, Baller JF, Gustafson LA, Burger HJ, Herling AW, Kuipers F, Meijer AJ, Reijngoud DJ: Acute inhibition of hepatic glucose-6-phosphatase does not affect gluconeogenesis but directs gluconeogenic flux toward glycogen in fasted rats. A pharmacological study with the chlorogenic acid derivative S4048. J Biol Chem. 2001, 276: 25727-25735. 10.1074/jbc.M101223200.
Perdigoto R, Rodrigues TB, Furtado AL, Porto A, Geraldes CF, Jones JG: Integration of [U- 13C]glucose and 2H2O for quantification of hepatic glucose production and gluconeogenesis. NMR Biomed. 2003, 16: 189-198. 10.1002/nbm.826.
Marin S, Lee WN, Bassilian S, Lim S, Boros LG, Centelles JJ, FernAndez-Novell JM, Guinovart JJ, Cascante M: Dynamic profiling of the glucose metabolic network in fasted rat hepatocytes using [1,2-13C2]glucose. Biochem J. 2004, 381: 287-294. 10.1042/BJ20031737.
Katz J, Wood HG: The use of glucose-14 C for the evaluation of pathway of glucose metabolism. J Biol Chem. 1960, 235: 2165-2177.
Mason GF, Rothman DL: Basic principles of metabolic modeling of NMR 13C isotopic turnover to determine rates of brain metabolism in vivo. Metab Eng. 2004, 6: 75-84. 10.1016/j.ymben.2003.10.003.
Schmidt K, Carlsen M, Nielsen J, Villadsen J: Modeling isotopomer distributions in biochemical networks using isotopomer mapping matrices. Biotechnol Bioeng. 1997, 55: 831-840. 10.1002/(SICI)1097-0290(19970920)55:6<831::AID-BIT2>3.0.CO;2-H.
Zupke C, Stephanopoulos G: Modeling of isotope distributions and intracellular fluxes in metabolic networks using atom mapping matrices. Biotechnol Prog. 1994, 10: 489-498. 10.1021/bp00029a006.
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: 69-85. 10.1002/(SICI)1097-0290(1999)66:2<69::AID-BIT1>3.0.CO;2-6.
Wiechert W, de Graaf AA: In vivo stationary flux analysis by 13C labeling experiments. Adv Biochem Eng Biotechnol. 1996, 54: 109-154.
Wiechert W, de Graaf AA: Bidirectional Reaction Steps in Metabolic Networks: I. Modeling and Simulation of Carbon Isotope Labeling Experiments. Biotechnol Bioeng. 1997, 55: 101-117. 10.1002/(SICI)1097-0290(19970705)55:1<101::AID-BIT12>3.0.CO;2-P.
Wiechert W, Siefke C, de Graaf AA, Marx A: Bidirectional Reaction Steps in Metabolic Networks: II. Flux Estimation and Statistical Analysis. Biotechnol Bioeng. 1997, 55: 118-135. 10.1002/(SICI)1097-0290(19970705)55:1<118::AID-BIT13>3.0.CO;2-I.
Mulquiney PJ, Kuchel PW: Modelling Metabolism with Mathematica. 2003, CRC Press, Boca Raton, FL
Selivanov VA, Puigjaner J, Sillero A, Centelles JJ, Ramos-Montoya A, Lee PW, Cascante M: An optimized algorithm for flux estimation from isotopomer distribution in glucose metabolites. Bioinformatics. 2004, 20: 3387-3397. 10.1093/bioinformatics/bth412.
Selivanov VA, Meshalkina LE, Solovjeva ON, Kuchel PW, Ramos-Montoya A, Kochetov GA, Lee PW, Cascante M: Rapid simulation and analysis of isotopomer distributions using constraints based on enzyme mechanisms: an example from HT29 cancer cells. Bioinformatics. 2005, 21: 3558-3564. 10.1093/bioinformatics/bti573.
van Winden W, Verheijen P, Heijnen S: Possible Pitfalls of Flux Calculations Based on 13C-Labeling. Metab Eng. 2001, 3: 151-62. 10.1006/mben.2000.0174.
Zupke C, Stephanopoulos G: Modeling of isotope distributions and intracellular fluxes in metabolic networks using atom mapping matrices. Biotechnol Prog. 1994, 10: 489-498. 10.1021/bp00029a006.
Kleijn RJ, van Winden WA, van Gulik WM, Heijnen JJ: Revisiting the C-label distribution of the non-oxidative branch of the pentose phosphate pathway based upon kinetic and genetic evidence. FEBS J. 2005, 272: 4970-4982. 10.1111/j.1742-4658.2005.04907.x.
Dennis JE, Schnabel RB: Numerical Methods for Unconstrained Optimization and Nonlinear Equations. 1983, Prentice-Hall, Englewood Cliffs, New Jersey, USA
Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C++: The Art of Scientific Computing. 2002, Cambridge
Shulman RG, Rothman DL: Interpreting functional imaging studies in terms of neurotransmitter cycling. Proc Natl Acad Sci USA. 1998, 95: 11993-11998. 10.1073/pnas.95.20.11993.
Posner MI, Raichle ME: The neuroimaging of human brain function. Proc Natl Acad Sci USA. 1998, 95: 763-764. 10.1073/pnas.95.3.763.
Fitzpatrick SM, Rothman D: New approaches to functional neuroenergetics. J Cogn Neurosci. 1999, 11: 467-471. 10.1162/089892999563454.
Shepherd GM: The Synaptic Organization of the Brain. 1994, New York: Oxford Univ Press
Zwingmann C, Richter-Landsberg C, Leibfritz D: 13C Isotopomer Analysis of Glucose and Alanine Metabolism Reveals Cytosolic Pyruvate Compartmentation as Part of Energy Metabolism in Astrocytes. Glia. 2001, 34: 200-212. 10.1002/glia.1054.
Cruz F, Villalba M, Garcia-Espinosa MA, Ballesteros P, Bogonez E, Satrustegui J, Cerdan S: Intracellular compartmentation of pyruvate in primary cultures of cortical neurons as detected by (13)C NMR spectroscopy with multiple (13)C labels. J Neurosci Res. 2001, 66: 771-781. 10.1002/jnr.10048.
Sonnewald U, Schousboe A, Qu H, Waagepetersen HS: Intracellular metabolic compartmentation assessed by 13C magnetic resonance spectroscopy. Neurochem Int. 2004, 45: 305-310. 10.1016/j.neuint.2003.10.010.
Abraham MR, Selivanov VA, Hodgson DM, Pucar D, Zingman LV, Wieringa B, Dzeja PP, Alekseev AE, Terzic A: Coupling Of Cell Energetics With Membrane Metabolic Sensing Integrative Signaling Through Creatine Kinase Phosphotransfer Disrupted By M-Ck Gene Knock-Out. J Biol Chem. 2002, 277: 24427-24434. 10.1074/jbc.M201777200.
Selivanov VA, Alekseev AE, Hodgson DM, Dzeja PP, Terzic A: Nucleotide-gated KATP channels integrated with creatine and adenylate kinases: Amplification, tuning and sensing of energetic signals in the compartmentalized cellular environment. Mol Cell Biochem. 2004, 256/257: 243-256. 10.1023/B:MCBI.0000009872.35940.7d.
Rich TC, Fagan KA, Nakata H, Schaack J, Cooper DM, Karpen JW: Cyclic nucleotide-gated channels colocalize with adenylyl cyclase in regions of restricted cAMP diffusion. J Gen Physiol. 2000, 116: 147-161. 10.1085/jgp.116.2.147.
Cornish-Bowden A: Fundamentals of Enzyme Kinetics. 2004, Portland Press, London, 3
This work was supported by the grants: Fundation la Caixa (ONO3-70-0), the Ministerio de Ciencia y Tecnologia of Spanish Government (SAF2005-01627 and PPQ2003-06602-C04-04); NIH DK56090-Al (to W.N.P.L); Generalitat de Catalunya (ABM/acs/PIV2002-32) to (V.A.S); Generalitat de Catalunya (2004 PIV2 14) to (T.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 authors also acknowledge the support of the Bioinformatic grant program of the Foundation BBVA and the Comissionat d'Universitats i Recerca de la Generalitat de Catalunya.
This article has been published as part of BMC Neuroscience Volume 7, Supplement 1, 2006: Problems and tools in the systems biology of the neuronal cell. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcneurosci/7?issue=S1.