Abstract
Background
Stimulus Response Experiments to unravel the regulatory properties of metabolic networks are becoming more and more popular. However, their ability to determine enzyme kinetic parameters has proven to be limited with the presently available data. In metabolic flux analysis, the use of ^{13}C labeled substrates together with isotopomer modeling solved the problem of underdetermined networks and increased the accuracy of flux estimations significantly.
Results
In this contribution, the idea of increasing the information content of the dynamic experiment by adding ^{13}C labeling is analyzed. For this purpose a small example network is studied by simulation and statistical methods. Different scenarios regarding available measurements are analyzed and compared to a nonlabeled reference experiment. Sensitivity analysis revealed a specific influence of the kinetic parameters on the labeling measurements. Statistical methods based on parameter sensitivities and different measurement models are applied to assess the information gain of the labeled stimulus response experiment.
Conclusion
It was found that the use of a (specifically) labeled substrate will significantly increase the parameter estimation accuracy. An overall information gain of about a factor of six is observed for the example network. The information gain is achieved from the specific influence of the kinetic parameters towards the labeling measurements. This also leads to a significant decrease in correlation of the kinetic parameters compared to an experiment without ^{13}Clabeled substrate.
Background
The recent developments in Metabolomics and Fluxomics open up new perspectives for detailed in vivo studies of the cellular metabolism. Various experimental approaches using GC or LCMS enable the measurement of intracellular concentrations and isotope enrichments with increasing accuracy. Nowadays, almost all metabolites of the central carbon metabolism can be detected [13]; although not all can be quantified exactly. Metabolic fluxes depend on the metabolite concentrations, enzyme concentrations (regulated by transcription [4]) as well as regulatory mechanisms (like enzyme modification [5] and allosteric inhibition or activation [6]). Different information about the in vivo mechanisms and fluxes can be collected from already established experimental approaches.
Figure 1 depicts the most prominent approaches in Fluxomics and Metabolomics and their information targets. In general, experimental approaches can be divided into metabolic stationary and metabolic nonstationary experiments. Whereas metabolic stationary approaches aim to quantify intracellular fluxes, nonstationary experiments are used to identify reaction kinetic properties like substrate affinities and allosteric inhibition. Fluxes at different metabolic states can only be calculated from a reaction kinetic model (with identified parameters). The following approaches are currently established:
Figure 1. Category of experiments. Overview of experiments used for the estimation of intracellular fluxes and enzyme kinetics. The new metabolic and isotopic nonstationary state is analyzed here within a simulation study.
1. Metabolic Flux Analysis (left state in Figure 1): At metabolic steady state it is assumed, that the intracellular concentrations and fluxes are not changing during the time of the analysis. A model based on mass balances is used to quantify the flux through the network [7].
2. ^{13}C Metabolic Flux Analysis (second left state in Figure 1): It is well known that solely stoichiometric balancing will not lead to a resolution of bidirectional or parallel fluxes [8]. Experiments with additional ^{13}C labeled substrates have been shown to increase the measurable information and to enable the quantification of exchange fluxes and parallel fluxes [9,10].
3. ^{13}C Metabolic Flux Analysis at isotopically nonstationary state (third state in Figure 1). Here the organism is kept at metabolic stationary state while a ^{13}C labeled substrate pulse is applied. Only recently, Nöh et al. [11] have evaluated an experiment of this type. It was shown that the duration of the labeling experiment can be drastically reduced from several hours to the order of minutes.
4. StimulusResponse Experiments (fourth state in Figure 1). At a metabolic nonstationary state the intracellular concentrations (and consequently also the fluxes) are not constant. These conditions target the identification of in vivo enzyme kinetics with one single experiment. The currently most prominent approach in this field is given by StimulusResponse Experiments [1215]. A culture is driven to a substrate limited state and then excited by a strong external stimulus, typically a substrate pulse. The metabolic response of the cells is tracked by rapid sampling on a subsecond timescale [1618]. Although lots of data is generated from a stimulus response experiment [16,19], the approach has shown to be limited: still only a part of the parameters can be identified. Even if the data from several experiments is gathered, it is not possible to determine all kinetic parameters [20].
These experiments are now common procedure. Yet missing is the experiment that results from a fusion of approach 3 and 4, i.e. labeling the substrate in a metabolically dynamic experiment. This type of experiment is theoretically analyzed to evaluate whether this type of labeling experiment can increase the parameter estimation accuracy of a reaction kinetic model. To this end, it is first shown how the experiments can be modeled based on the known concepts from ^{13}C Flux Analysis and reaction kinetic modeling. A series of papers [2123] are concerned with the integration of kinetic and isotopic information. In contrast to the present integration, the authors assumed most enzyme kinetic parameters to be known. Moreover the system was analyzed in a metabolic steady state and therefore resembles approaches 2 and 4 under the perquisite of mostly known kinetic parameters. The values for these parameters are taken from literature or additional experiments [23]. The approach presented here is based on a dynamic experiment rather than one or several metabolic steady states and all enzyme kinetic parameters are to be estimated from the experimental data.
The approach is based on a complete isotopomer model that is capable of exploiting all available isotopomer measurements from MS(/MS) or NMR devices. In order to demonstrate the general mathematical framework a small example reaction network is discussed the first time. Comparing the outcomes of an experiment with labeling to those of a reference without the addition of ^{13}C substrate the information surplus from the metabolic and isotopic transient state is elucidated. To cover different experimental settings, three scenarios of available measurements are introduced and compared.
Modeling metabolic and isotopic nonstationary systems
A general model structure for metabolic and isotopic nonstationary systems can be obtained by combining models from both worlds. However, this hybrid structure also inherits some standard assumptions from kinetic and isotopic models. Precisely, all upcoming equations are only valid if the following assumptions [24] are accepted:
1. Homogeneity: All intracellular metabolites are distributed uniformly.
2. No isotopic effects: The enzymatic reactions are not influenced by the labeling state of the metabolites.
3. Fast equilibrium of enzyme complexes: Intermediate enzyme complexes rapidly reach equilibrium such that the quasistationary assumption of MichaelisMenten [6] is valid.
Reaction kinetic model
With these assumptions, well known mechanistic or approximate enzymatic rate laws like MichaelisMenten [6] can be used to describe the fluxes v in dependence of substrate, product and effector concentrations, kinetic parameters α (like v_{max}, K_{M}) and external (experimental) parameters β (like substrate feed):
The concentrations of metabolites (c) do change depending on the in and effluxes v. The mass balances of the pools can be formulated using the well known stoichiometric matrix N [25]:
^{13}C Balance equations
Here, the labeling state is modeled using the isotopomer concept [24] to keep the contribution understandable on a basic level. Nevertheless, more evolved but mathematically equivalent concepts like cumomers [26], bondomers [27] or elementary metabolite units [28] can also be used. A metabolite with n carbon atoms has 2^{n }labeling states. Thus, 2^{n }isotopomer fractions are required for each metabolite to describe its labeling state. The overall isotopic state of the network is described by a vector x containing the isotopomer fractions of all metabolites. Isotopomer dynamics of a metabolite pool depends on the pool concentration c, the input labeling x^{inp }and the flux functions v. The overall system for the labeling state is given by [29]:
Here, the operator generates a diagonal matrix composed of concentrations c_{i }(with n_{i }repetitions which correspond to all isotopomer fractions, n_{i }being the number of isotopomers of metabolite i). The nonlinear function f includes the stoichiometry of the isotopomer transitions, the vector x describes the isotopomer distribution of the balanced intracellular pools (for an example take a look at Eq. (14) in the example section) and x^{inp }is defined by the isotopomer fractions of the input substrate.
Notice that unlike in stationary isotopomer balances [30], the vector v is not constant but changes over time in dependence of the metabolite concentrations. Combining both model equations (2) and (3), the metabolic and isotopic transients can now be completely described. Hence, the resulting model is a hybrid of a classical ^{13}C labeling and a reaction kinetic model. The functions f is constructed based on matrix calculus as described earlier by Wiechert et al [24] but replacing the vector v with a vector of kinetic rate terms. A clear difference appears at the point of solution; Instead of a numeric equation solver, numeric integration is applied to calculate the labeling transients.
Measurement model
The measurement model is a function that couples the dynamic isotopomer model with the available measurements. Within this section the measurement model will be described only very briefly in order to avoid rather unnecessary technical explanations (these can be found in Appendix A). Generally, at time point t_{k }measurements of the labeling and the concentrations are expected. Both quantities are gathered in the measurement vector y(t_{k}) that depends on the state of the system c(t_{k}) and x(t_{k}). Hence, we introduce a measurement function g:
Collecting all sampling time points into a vector ξ = (t_{1},...t_{N})^{T}, Eq. (4) compactly reads:
A concrete function g is presented within the example given in the appendix for massisotopomer measurements.
Parameter sensitivity
The final goal of the nonstationary experiment is the determination of the kinetic parameters α. All other quantities (fluxes v, concentrations c) follow from these parameters. Sensitivities quantify the influence of α on the model response, namely the concentration time course c(t) and the isotopomer distributions x(t). Implicit derivation of the differential equations (2) and (3) with respect to the parameters α leads to the sensitivity differential equations associated with the model (2)–(3):
Besides the influence of the parameters α on the system state, the influence on the expected measurements is now easily calculated using the chain rule:
As shown in the following sections, the sensitivities are essential for the statistical evaluation and comparison of different experimental configurations [10].
Nonlinear regression model
Regression analysis is the chosen method to estimate the kinetic parameters from the experimental data. Assuming that the model adequately describes reality and that experimental measurements y have some unknown error ε, Eq. (5) extends to:
Parameter estimation is performed by using weighted least square minimization. Expecting that the errors ε(ξ) are normally distributed with expectation E[ε(ξ)] = 0 and the given measurement covariances are collected in the measurement covariance matrix Σ, the weighted least squares functional for the calculation of the optimal parameters reads [31]:
Various algorithms can be applied to solve this minimization problem like derivationfree simplex methods [32], evolutionary strategies [33], or gradientbased methods [34].
Statistical evaluation
In order to obtain a statement about whether or not the information contained in the measurements is sufficient to identify the parameter values within a certain precision, it is essential to analyze the accuracy of the parameter estimation found by solving Eq. (8). Assuming that the estimate is close to the real solution α* application of linearized regression analysis is feasible. In experimental design studies usually a parameter set = α* is chosen. It is known, that linear statistics can produce imprecise estimates of parameter confidence regions [35]. Nevertheless, this tool still fulfils its purpose when it is used to compare different experiments. Here the relative improvement of the confidence regions is much more important than their precise knowledge. Moreover, the concept of correlation coefficients that characterizes the confidence regions is an inherently linear concept. It will turn out that parameter correlations are relevant because with ^{13}C labeled substrates parameter correlations strongly decrease that lead to smaller confidence regions.
The covariance of the parameters is calculated using the derivative and the covariance of the measurements Σ by [36]:
In the appendix section it is explained how the different types of measurements and additional knowledge are integrated.
Statistical measures to compare experiments
The parameter covariance matrix contains all information needed to judge the information gain of an experiment, namely parameter standard deviations, correlations and confidence regions. A drawback of using the covariance matrix is its size. Even for small models it is rather difficult to evaluate and compare experiments on its basis. Therefore statistical measures are used that reduce the complexity and facilitate an automated comparison of experiments. Usually, a reference experiment is specified to which all other experiments are compared.
A commonly used measure is the Dcriterion, which is calculated from the determinant of the covariance matrix. It comprehensively summarizes the available parameter information. To compare different experiments Möllney et al. [10] propose to use a root of the Dcriterion. The scaled ratio I_{D }between an experiment and the reference experiment (index ref) is defined by:
Roughly speaking the I_{D }value reflects an overall information gain based on the volume of the covariance ellipsoid (for details see [10]). The standard deviations and correlations of the judged experiment are – on an average – I_{D }times smaller than those from the reference experiment. However, this does not guarantee that every single standard deviation is improved.
Besides the Dcriterion, the Acriterion is frequently applied. The latter criterion reflects the mean of the diagonal elements of the covariance matrix. Again, a scaled scalar is introduced to compare two experiments:
In contrast to the Dcriterium, correlations between parameters are not taken into account; the Acriterium is based on the parameter variances only. Therefore, the square root has to be taken.
An exploratory study
Within this section the described framework is applied to an example network. As an example, some equations are explicitly shown to facilitate the understanding of the presented concepts. Some different conceivable scenarios regarding the available measurements are also introduced to cover different measurement setups and analyze the influence of, for example, a missing concentration measurement.
Example network
The network chosen for this study is shown in Figure 2. Although it is simple, it nonetheless reflects all typical reaction mechanisms of the central carbon metabolism. A linear reaction sequence (v_{upt}, v_{1}, v_{2}) with a feedback inhibited entry reaction (v_{upt}) as well as a cyclic reaction sequence (v_{3}, v_{4}) with some fillup reaction (v_{5}) and an exit (v_{6}) is included. Some reversible reaction steps (v_{1}, v_{5}) are incorporated as well as a bimolecular reaction step (v_{2}) with Hilltype kinetic mechanism and irreversible split reactions (v_{3}, v_{4}) are present (see Table 1 for a complete listing of reaction mechanisms and the kinetic expressions).
Figure 2. Example network structure including atom transitions. The reaction kinetics are given in Table 1.
Table 1. Kinetic rate laws used in the small example network.
The assumed metabolic regulation reflects some basic principles like product inhibition (reaction v_{upt }by A, v_{2 }by C) and allosteric inhibition (v_{4 }by C). Short chain molecules were chosen to keep the dimension of the isotopomer equations low. The substrate of the system contains only 2 carbon atoms (A_{feed}) and the molecule formed in the bimolecular reaction entering the cycle has 4 carbon atoms.
Isotopomer balance equations and simulation
Some assumptions about the experimental setup and the network will be needed to perform the simulation study of the nonstationary state: (1) the external influences β, (2) assumptions about the reaction mechanisms that will form the functions v and (3) estimates for the unknown kinetic parameters α.
The balance equations for the metabolite pools are not shown here because they are well documented in literature [7]; five pools are balanced for the example network. In total, 24 kinetic parameters are needed for the reaction kinetics. In order to model the reactor feeding profile, one input parameter () is introduced. The values used within this study are listed in Table 2. Having the balances in hand, the time course of the concentrations c(t) is calculated by integrating Eq. (2). Figure 3 shows the simulation result. Concentrations of the metabolite pools A and B show a short overshoot, while concentrations C, D and E show a monotonic increase after the pulse. Pools A and B rapidly reach a steady state after about 10 s, while the dynamics in C, D and E do not reach a new steady state in the given simulation time of 20 s.
Figure 3. Simulation of the example network and the assumed measurement standard deviations. The extracellular concentration (A_{ex}) rises from 0.05 mM to 1.5 mM by the pulse. The intracellular concentrations A and B show a short overshoot, C, D and E increase continuous and reach the steady state after about 25 s. The squares show the noise free measurements generated from the simulation. The standard deviations include the errors from sample processing and the MS measurements.
Table 2. Kinetic parameters for the simulation study.
The atom transitions have to be known to balance the isotopomer fractions (Eq. (3)). These are shown in Figure 2 and are also listed in Table 3. Forty isotopomer fractions are balanced in the example. Four input variables are needed to describe the labeling of A_{ex }(here set to x^{inp }= (0.01,0.01,0.01,0.97)^{T}). Exemplary, the balance equations for the isotopomer fractions a = (a_{00}, a_{01}, a_{10}, a_{11})^{T }of pool A are shown.
Table 3. Atom transitions for the small example network
Again it must be noted that the fluxes are functions of kinetic parameters (see Table 1). For example, the flux v_{1 }(A → B) is a reversible flux. Its rate is given by:
With known initial isotopomer fractions x_{0 }at time point t = 0 (usually determined by the natural isotope enrichments), the integration of Eq. (3) yields the time course x(t) of the isotopomer enrichments. In Figure 4 it can be seen that pools A and B are rapidly enriched with labeled carbon. Only at 2 s respectively at 5 s after the pulse they do show the labeling of the input substrate A_{ex}. Pools C, D and E exhibit a different behavior. Especially in C and D an interim labeling state can be observed, that results from the reaction of unlabeled E#00 with fully labeled B#11. This state is still seen in D while in E it is not as pronounced because E is fed from B and D. This already hints at additional information contents, because a more detailed behavior in comparison to the one shown Figure 3 is observed. In the following sections, it will be seen that the kinetic parameters do influence the isotopomer distribution differently to the concentration time course.
Figure 4. Simulation of the isotopomer enrichments. The extracellular substrate A_{ex }immediately switches from unlabeled to labeled. The enrichment of A and B is fast (steady state after <10 s), pool C shows intermediate labeling patterns before reaching the fully labeled steady state. These intermediate patterns can also be seen in D and E.
Scenarios
Since there is still much development in the analytical methods, scenarios with different measurements are investigated. In order to estimate the influence of increased or decreased availability of measurement information the following cases will be compared (summarized in Figure 5):
Figure 5. Summary of the compared measurement scenarios. In S_{all }it is assumed, that all concentrations and mass isotopomers can be measured. Scenario S_{Drel }is less informative, the concentration of metabolite D cannot be quantified. S_{C,Dfrag }assumes MS/MS measurements with fragmentation for the metabolites C and D.
• Scenario S_{all}: This is an optimistic scenario where all metabolite concentrations can be quantified and all mass isotopomer fractions are available. Here, the mass isotopomers of the entire carbon skeleton are measured without any fragmentation.
• Scenario S_{Drel}: Here, less information is available compared to S_{all}; this scenario is derived from findings that for some metabolites the concentration could not be determined [11]. Nevertheless, it was possible to measure the mass isotopomer ratios. As an example, the concentration of D cannot be determined, but its mass isotopomer ratios are measured.
• Scenario S_{C,Dfrag}: Compared to scenario S_{all }additional mass measurements of fragments are available. As an example, the labeling distribution of the pools C and D are measured in more detail because of a specific fragmentation in MS/MS measurement mode (see Figure 5). Metabolite D fragments into a C_{2 }and a C_{1 }skeleton. Here, the mass isotopomers of the C_{2 }body can be measured additionally. Instead of four mass isotopomers, in total six tandem mass isotopomer measurements are generated [37]. Metabolite C fragments into a C_{3 }and a C_{1 }body. Consequently, instead of five mass isotopomer ratios, eight tandem mass isotopomers are measured.
To summarize, scenario S_{all }is optimistic, S_{C,Dfrag }is even more optimistic, whereas in S_{Drel }some information is missing.
Sampling and quality of the measurements
Besides the measurements, the sampling frequency is also crucial for parameter identification [29]. For all scenarios, the following sampling setup is used:
• Three samples are taken every 1 s before the pulse (3...0 s),
• In the following 5 s samples are taken every 0.5 s (0.5...5 s),
• After 7 s samples are taken every 2 s (7...19 s)
The quality of the measurements depends on the accuracy of the measurement device and possible errors during sample preparation. Figure 6 depicts the sampling steps usually needed to determine intracellular concentrations. Essentially, the sampling steps lead to a dilution φ. The signal intensity depends on the concentration in the sample and the response factor γ that is determined from standard additions [38]. More details about the measurement model can be found in the appendix section. All sampling steps and possible errors have to be taken into account for the statistical analysis. Based on the current knowledge for the standard deviations the following factors are considered:
Figure 6. Scheme of the sample processing steps. The single sampling steps and their effects on concentration and labeling measurements.
• The response factors γ (set to 1 if determined) have a standard deviation σ_{γ }of 1%.
• The dilution factors φ (ξ) (also set to 1) have a standard deviation of 3%.
• The MS(/MS) measurement (area) itself is composed of a relative error which is assumed in the order of 1% and an absolute error that was set to 10 nM. The latter one appears rather high, but as the dilution usually will be around a factor 20 it is within a realistic range:
Results
Statistical analysis
The example system consists of 40 isotopomer and six concentration states with 24 parameters. Hence, the parameter sensitivity matrix (Eq. (6)) contains 46·24 = 1104 entries for each of the 19 sampling time points. Although interesting, a detailed analysis of these large matrices is extraordinarily timeconsuming. In fact, the influence on the measurable signals is more relevant. As it can be seen from equations, and the measurement equations are linear. Thus, the calculation of the output sensitivity is straight forward (Eq. (7)).
Figure 7 shows the course of the sensitivities with respect to kinetic parameters of the uptake reaction v_{upt }(K_{I}) and reaction v_{1 }(K_{mA}, v_{max }and K_{eq}). The sensitivities have to be compared on a comparable scale. Recall that the mass isotopomers are given as fractions of the pool. Multiplying the fractions with the pool concentration, the mass isotopomer concentrations are immediately obtained by c·x.
Figure 7. Parameter sensitivity plots. First row: Sensitivities of the concentrations to the parameters of reaction v_{1 }(K_{mA}, v_{max }and K_{eq}) and v_{upt }(K_{I}). Second to last row: Sensitivities of selected mass isotopomer concentrations to the same parameters.
These quantities can be directly compared to metabolite concentrations. Their parameter sensitivities can be simply obtained from Eq. (7) by using the chain rule. The sensitivity matrix dy/dα includes the sensitivities of the concentrations and of the mass isotopomers with respect to kinetic parameters. These sensitivities are plotted in Figure 7. It is observed that:
• After the pulse there is a strong increase in sensitivity of mass isotopomer concentrations of fully labeled pools (A+2, B+2, E+2, with a less extend C+4, D+3)
• The sensitivity of unlabeled concentrations (A+0, B+0, C+0, D+0, E+0) reaches zero, with some intermediate peaks.
• Some parameters show maxima in sensitivity, e.g. d(B·b_{+0})/dK_{mA}, d(C·c_{+2})/dK_{mA}, d(E·e_{+2})/dK_{mA}.
• Some parameters first show positive, then negative influence on the mass isotopomers (and vice versa), e.g. d(B·b_{+0})/dK_{mA}, d(B·b_{+0})/dK_{eq}, d(B·b_{+0})/dv_{max}.
• Most parameter sensitivities reach a new steady state within the given simulation interval of 20 s, e.g. d(A·a_{+2})/dK_{mA}, d(A·a_{+2})/dK_{eq}, d(A·a_{+2})/dv_{max}.
These findings clearly show that the kinetic parameters do not only influence the total metabolite concentrations, but also influence the labeling dynamics. Within the statistical analysis section it will be seen that this will help to achieve an identification of the enzyme parameters.
Parameter standard deviations
Figure 8 shows a cumulative plot of the estimated relative standard deviations of all parameters. Obviously, the majority of the parameters can hardly be determined from the (unlabeled) reference experiment. The best estimate is found to have about 46% relative error. The second step is seen at about 61%, which is the second best estimate. Three parameters are found within a standard deviation of about 137%. These high deviations are observed because the scenario includes 3% sample dilution error and additional noise from the MS device. Especially for low concentrated metabolites (like E, see Figure 3), high standard deviations occur due to the ground noise (10 nM) of the MS measurements. It has to be taken into account too, that there are more than the 24 kinetic parameters. Another 19 × 5 + 5 = 100 parameters (five φ 's for each measured metabolite and each sampling time point and five γ 's) are introduced for the measurement model. Although these parameters are measured directly from the sampling volume, added reagents and standard additions, these additional parameters do consume some information. It is seen from the statistical analysis that these parameters are also influenced and their estimation accuracy gets higher compared to the experimental errors.
Figure 8. Cumulative distribution of the estimated standard deviations. In the unlabeled reference experiment the most accurate parameter is determined with a relative standard deviation of about 46% (first step of the black line). Only two parameters can be determined with less than 70% standard deviation. In contrast, using a labeled substrate, 10 parameters (about 40% of the parameters) can be determined with a standard deviation of less than 70%.
A pronounced increase in estimation accuracy of the kinetic parameters is seen when comparing scenario S_{all }with the unlabeled reference (Figure 8). The line is shifted about a factor of 10 to the left. The increase is even more pronounced for the most accurate estimation (0.92% that is about a factor of 50). With up to 70% standard deviation, 10 parameters are found for the labeled experiment compared to only two for the reference.
When comparing the three labeled scenarios only little differences can be found. With some additional labeling information (scenario S_{C,Dfrag}) the estimates become slightly more accurate (red, dashdotted line). A missing concentration (scenario S_{Drel}) can obviously be compensated by the labeling information. The dashed line is shifted to the right slightly. The most accurate estimate has a standard deviation of 0.93%.
These findings are also reflected in the D and Acriterion (Table 4). The Acriterion is based on the variances and reflects an information increase that is comparable to observations from Figure 8. The increase in information is 11.01 (scenario S_{all}). The information gain from scenario S_{C,Dfrag }is about 1% higher; for S_{Drel }about 1% lower.
Table 4. Comparison of the D and Acrtiteria.
The information gain calculated from the Dcriterion is lower. Here a value of 6.91 is calculated for scenario S_{all}. With some more available labeling measurements (scenario S_{C,Dfrag}) a factor of 7.31 is found. The information gain for scenario S_{Drel }decreases only marginally compared to S_{all }(6.80).
Parameter correlations
So far, the improvement of the absolute variances of all parameters was investigated. In order to characterize the influence of the additional information from ^{13}C labeling on the relationship among different parameters, the correlation matrix (Figure 9) will now be analyzed. The correlations for the reference experiment and the labeled scenario S_{all }are shown. The parameters are grouped with respect to the reactions and a block structure becomes visible.
Figure 9. Color visualization of the correlation matrix. The parameter correlation matrix for the reference experiment (left) and the labeled experiment S_{all }(right). Deep blue represents a correlation of 1, deep red represents 1.
At first glance a strong decrease in correlation can be observed on the areas outside of the main diagonal blocks. There are blocks with only little changes close to the main diagonal. These blocks are formed from groups of parameters that belong to one single reaction step. For example, the kinetic parameters v_{max }and K_{mA }of the reaction v_{1 }are strongly positive correlated (bue). I.e. if one is estimated too high, the other will also be estimated too high. Also in the labeled scenario, these parameters cannot be determined separately. It is well known from enzyme kinetic studies that v_{max }and K_{mA }values can hardly be independently determined, if the experimental concentration range is not chosen correctly.
In the following it will be seen that the effect of improved parameter distinction between different reaction steps is based on indirect flux measurements. The parameters K_{mA }(of reaction v_{1}) and K_{I }(of v_{upt}) were chosen here as an example, because a strong decrease in correlation is observed. In Figure 7 the sensitivities of these parameters are plotted for the reference and the scenario S_{all}. Clearly, parameters with distinct influence on measurements will be less correlated. Within the reference the parameters exhibit a rather comparable influence on all measurements (though the magnitude is distinct). As a consequence, positive correlation is observed.
If labeling is introduced, the picture gets different. The influences on the measured mass isotopomers show some behavior that leads to a decoupling of parameters especially within the first 5 s after the pulse. The intermediate labeling patterns, particularly E+2, B+2 as well as E+0 and B+0 exhibit a distinct behavior for parameters K_{mA }and K_{I}. Whereas K_{mA }exhibits a strong influence on the time course of these mass isotopomers, the influence of K_{I }is relatively small. Additionally, the sensitivities of K_{mA }show a strong time dependency that increases the independent measurement of this parameter.
Conclusion
The exploratory study strongly suggests that performing a stimulus response experiment with ^{13}C labeled substrate and measuring the intracellular labeling states can significantly increase the accuracy of kinetic parameter estimation. A decrease in the variance of more than a factor of 10 was observed for the analyzed example network (Figure 8). With labeling, three parameters are determined with less than 10% relative standard deviation. More than 40% of the parameters (i.e. 10 of 24) can be determined with less than 70% standard deviation (Figure 8).
Additionally it was found that various parameter correlations are drastically reduced (Figure 9). When comparing the experiment with labeled substrate to its reference counterpart a seven fold information gain (based on the Dcriterion) was observed. Interestingly, the information loss due to the unknown concentration of pool D (scenario S_{all}) can be compensated by the measured mass isotopomer distributions. This scenario shows only slightly less accurate determined parameters (Figure 8 and Table 4)
It can clearly be seen from the sensitivities that the parameters influence the labeling time course more specifically than the concentration time course. When only measuring the concentrations, a concentration increase can be explained by an increased influx or a decreased efflux. The addition of isotopic transients to the concentration courses serves as additional information about the pool exchange. A higher influx will lead to a faster labeling enrichment, while a decreased efflux will slow down the labeling enrichment. Thus, the labeling transients include information about the pool exchange that delivers valuable constraints for parameter identification.
In the case of isotopically nonstationary experiments under metabolic steadystate conditions, a similar example network was analyzed to estimate the information gain compared to a classical ^{13}C labeling experiment at isotopically steadystate [39]. The conclusions from the simulation study have shown to be in accordance with findings from a complex E. coli network [11]. Therefore, the results obtained from the performed exploratory study should be transferable to real sized metabolic networks.
Clearly, an automated simulation tool is required to set up the equations based on the stoichiometry and the atom transitions. To solve these equations algorithms comparable to the tool of Nöh et al. [29] could be applied. Nöh et al. [29] have shown that the solution of approximately 4000 equations can be handled without problems. Compared to this number, the additional kinetic equations (1) can be neglected. Also, the sensitivities can be rapidly calculated using a cluster computer [29]. As further refinement, Antoniewicz et al. [28] demonstrated that the dimension of isotopomer balances can be further reduced using the concept of Elementary Metabolite Units (EMU).
Authors' contributions
SAW: Did the simulations and statistical evaluation, writing of the manuscript, KN: Wrote the section about statistical measures, correction of the manuscript, WW: Supervision of the presented work, correction of the manuscript.
All authors read and approved the manuscript.
Appendix: Details on the measurement model
Within scenario S_{all }all pool concentrations can be measured and the labeling state is observed by mass isotopomers from unfragmented molecules over time (Section "Scenarios"). For a metabolite with n carbon atoms, a total of n+1 different mass isotopomers will be measured. For example, the unlabeled pool A#00 will be found at a mass M(A)+0 (henceforth denoted as A+0). The signal A+1 is evoked by the singly labeled isotopomers A#01 and A#10. The fully labeled A#11 is found in the signal A+2. Usually some linear relationship between the signal intensity y_{A }and the concentration can be supposed. A scalable measurement model with some scaling factor ω is introduced by Möllney et al. [10]. Following this approach, the mass isotopomer measurements for pool A can be expressed by:
which can be shortly denoted by
The measurement equations with tandem mass isotopomers [37] are a little bit more evolved. The construction of the associated measurement matrix can be found in [37]. Its implementation is then done analogously to the example shown.
In contrast to an experiment without labeling, a metabolite pool is now distributed over various mass isotopomers depending on the labeling state. The total concentration has to be reconstructed from the single mass isotopomer measurements. Assuming that the labeling has no influence on the ionization, the calibration for the unlabeled metabolite can be used for all mass traces. Usually some linear relationship between the signal intensity y and the sample concentration is supposed, given by the response factor γ. The sample concentration originates from various preprocessing steps (Figure 6). Thus, the measured signal is obtained from the intracellular concentration c(t_{k}), the dilution factor φ(t_{k}) (sample preparation) and the response factor γ:
Some specialities of the different factors have to be pointed out:
1. The response factor γ is a device dependent factor that will usually be constant for all measurement time points, but different for each metabolite:
2. The dilution factor φ(t_{k}) reflects the sampling procedure and, thus, is time dependent (e.g. due to different amounts for neutralization). With regard to the metabolites, different scenarios could be taken into account. Current observations show, that during quenching metabolites leak to different extents [40]. Therefore in this study, a dilution factor for each metabolite is introduced that is independent. The vector of dilution factors thus reads:
Here a difference between the labeled and the nonlabeled experiment has to be noticed. The dilution factor φ(t_{k}) is a pool factor, the isotopomers of this pool are diluted to the same extent.
3. The scaling factor ω(t_{k}) is used for scaling the simulated labeling enrichments to the absolute scale of the measurement device. Consequently it can be replaced by the previously discussed scaling factors:
The errors of the single factors used for the simulation study have been mentioned in the main text. How these errors can be included to the regression model has not been discussed. The response factor γ is determined from a standard addition and is assumed to be constant over the measurement sequence. This factor is determined for each metabolite.
In contrast, φ(t_{k}) is needed for each measurement time point. Thus, the regression Eq. (8) has to be extended by the scaling factors γ and φ(t_{k}):
Construction of the measurement covariance matrix
The covariance matrix Σ contains the measurement variances that can be obtained from multiple measurements or more common from measurement validation. Its construction depends on the regression model (Eq. (21)). In case that g(x(ξ), c(ξ)) calculates the measurement vector (A_{+0}, A_{+1}, A_{+2}, B_{+0}, B_{+1}, B_{+2}, C_{+0}, C_{+1}, ... E_{+2}) at t_{1 }and in the following of t_{2 }until t_{N}. For sample preparation we have to take into account the dilution (Φ) of each metabolite at each time point. The accuracy of the concentrationresponse is determined for each metabolite. Thus, the matrix will have following structure:
Acknowledgements
We thank Dr. M.D. Haunschild for his support with the software tool MMT2.
References

Kopka J: Current challenges and developments in GCMS based metabolite profiling technology.
J Biotechnol 2006, 124:312322. PubMed Abstract  Publisher Full Text

Luo B, Grönke K, Takors R, Wandrey C, Oldiges M: Simultaneous determination of multiple intracellular metabolites in glycolysis, pentose phosphate pathway and tricarboxylic acid cycle by liquid chromatographymass spectrometry.
J Chromatogr A 2007, 1147:153164. PubMed Abstract  Publisher Full Text

VillasBoas SG, Mas S, Akesson M, Smedsgaard J, Nielsen J: Mass spectrometry in metabolome analysis.
Mass Spectrom Rev 2005, 24:613646. PubMed Abstract  Publisher Full Text

Caponigro GP R.: Mechanisms and control of mRNA turnover in Saccharomyces cerevisiae.
Microbiol Rev 1996, 60:233\&. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bendt AK, Burkovski A, Schäffer S, Bott M, Farwick M, Hermann T: Towards a phosphoproteome map of Corynebacterium glutamicum.
Proteomics 2003, 3:16371646. PubMed Abstract  Publisher Full Text

CornishBowden A: Fundamentals of Enzyme Kinetics. Edited by CornishBowden A. , Portland Press, London; 1995.

Stephanopoulos G, Aristidou A, Nielsen J: Metabolic Engineering: Principles and Methodogies. , Academic Press; 1998.

Wiechert W, Möllney M, Petersen S, de Graaf AA: A universal framework for C13 metabolic flux analysis.
Metab Eng 2001, 3:265283. PubMed Abstract  Publisher Full Text

Cohen SM, Glynn P, Shulman RG: 13C NMR study of gluconeogenesis from labeled alanine in hepatocytes from euthyroid and hyperthyroid rats.
Proc Natl Acad Sci U S A 1981, 78:6064. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Möllney 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:86103. PubMed Abstract  Publisher Full Text

Nöh K, Gronke K, Luo B, Takors R, Oldiges M, Wiechert W: Metabolic flux analysis at ultra short time scale: isotopically nonstationary 13C labeling experiments.
J Biotechnol 2007, 129:249267. PubMed Abstract  Publisher Full Text

Oldiges M, Takors R: Applying metabolic profiling techniques for stimulusresponse experiments: Chances and pitfalls.

Schäfer U, Boos W, Takors R, WeusterBotz D: Automated sampling device for monitoring intracellular metabolite dynamics.
Anal Biochem 1999, 270:8896. PubMed Abstract  Publisher Full Text

Theobald UM W.; Baltes, M.; Rizzi, M. & Reuss, M.: In vivo analysis of metabolic dynamics in saccharomyces cerevisiae: I. Experimental observations.
Biotechnol Bioeng 1997, 55:305316. Publisher Full Text

Visser D, van Zuylen GA, van Dam JC, Oudshoorn A, Eman MR, Ras C, van Gulik WM, Frank J, van Dedem GWK, Heijnen JJ: Rapid sampling for analysis of in vivo kinetics using the BioScope: A system for continuouspulse experiments.
Biotechnol Bioeng 2002, 79:674681. PubMed Abstract  Publisher Full Text

Buchholz A, Hurlebaus J, Wandrey C, Takors R: Metabolomics: quantification of intracellular metabolite dynamics.
Biomol Eng 2002, 19:515. PubMed Abstract  Publisher Full Text

Buziol S, Bashir I, Baumeister A, Claaßen W, NoisommitRizzi N, Mailinger W, Reuss M: New bioreactorcoupled rapid stoppedflow sampling technique for measurements of metabolite dynamics on a subsecond time scale.
Biotechnol Bioeng 2002, 80:632636. PubMed Abstract  Publisher Full Text

Lange HC, Eman M, van Zuijlen G, Visser D, van Dam JC, Frank J, de Mattos MJT, Heijnen JJ: Improved rapid sampling for in vivo kinetics of intracellular metabolites in Saccharomyces cerevisiae.
Biotechnol Bioeng 2001, 75:406415. PubMed Abstract  Publisher Full Text

Oldiges M, Kunze M, Degenring D, Sprenger GA, Takors R: Stimulation, monitoring, and analysis of pathway dynamics by metabolic profiling in the aromatic amino acid pathway.
Biotechnol Progr 2004, 20:16231633. PubMed Abstract  Publisher Full Text

Wahl SA, Haunschild MD, Oldiges M, Wiechert W: Unravelling the regulatory structure of biochemical networks using stimulus response experiments and largescale model selection.
IEE Proc Syst Biol 2006, 153:275285. Publisher Full Text

Selivanov VA, Puigjaner J, Sillero A, Centelles JJ, RamosMontoya A, Lee PW, Cascante M: An optimized algorithm for flux estimation from isotopomer distribution in glucose metabolites.
Bioinformatics 2004, 20:33873397. PubMed Abstract  Publisher Full Text

Selivanov VA, Meshalkina LE, Solovjeva ON, Kuchel PW, RamosMontoya 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:35583564. PubMed Abstract  Publisher Full Text

Selivanov VA, Marin S, Lee PWN, Cascante M: Software for dynamic analysis of tracerbased metabolomic data: estimation of metabolic fluxes and their statistical analysis.
Bioinformatics 2006, 22:28062812. PubMed Abstract  Publisher Full Text

Wiechert W, deGraaf AA: Bidirectional Reaction Steps in Metabolic Networks I. Modeling and Simulation of Carbon Isotope Labeling Experiments.
Biotechnol Bioeng 1997, 55:101117. Publisher Full Text

Schuster S, Fell DA, Dandekar T: A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks.
Nat Biotechnol 1999, 326333. PubMed Abstract  Publisher Full Text

Wiechert W, Möllney M, Isermann N, Wurzel W, de Graaf AA: Bidirectional reaction steps in metabolic networks: III. Explicit solution and analysis of isotopomer labeling systems.
Biotechnol Bioeng 1999, 66:6985. PubMed Abstract  Publisher Full Text

van Winden WA, Heijnen JJ, Verheijen PJT: Cumulative bondomers: A new concept in flux analysis from 2D [C13,H1] COSYNMR data.
Biotechnol Bioeng 2002, 80:731745. PubMed Abstract  Publisher Full Text

Antoniewicz MR, Kelleher JK, Stephanopoulos G: Elementary metabolite units (EMU): a novel framework for modeling isotopic distributions.
Metab Eng 2007, 9:6886. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Nöh K, Wahl A, Wiechert W: Computational tools for isotopically instationary C13 labeling experiments under metabolic steady state conditions.
Metab Eng 2006, 8:554577. PubMed Abstract  Publisher Full Text

Wiechert W: C13 metabolic flux analysis.
Metab Eng 2001, 3:195206. PubMed Abstract  Publisher Full Text

Bates DM, Watts DG: Nonlinear Regression Analysis and its Applications. New York, John Wiley & Sons; 1988.

Michalewicz ZF D. B.: How to Solve It: Modern Heuristics. 1st edition. Berlin, Springer; 2002.

Beyer HG: The Theory of Evolution Strategies. Berlin, Springer; 2001.

Deuflhard P, Bornemann F: Scientific Computing with Ordinary Differential Equations. In Texts in Applied Mathematics. Volume 42. Berlin, Springer; 2002.

Joshi M, SeidelMorgenstern A, Kremling A: Exploiting the bootstrap method for quantifying parameter confidence intervals in dynamical systems.
Metab Eng 2006, 8:447455. PubMed Abstract  Publisher Full Text

Chatterjee S, Hadi AS: Sensitivity Analysis in Linear Regression, Probability and Mathematical Statistics. New York, John Wiley & Sons; 1988.

Rantanen A, Rousu J, Kokkonen JT, Tarkiainen V, Ketola RA: Computing positional isotopomer distributions from tandem mass spectrometric data.
Metab Eng 2002, 4:285294. PubMed Abstract  Publisher Full Text

Bader A: A systematic approach to standard addition methods in instrumental analysis.

Nöh K, Wiechert W: Experimental design principles for isotopically instationary 13C labeling experiments.
Biotechnol Bioeng 2006, 94:234251. PubMed Abstract  Publisher Full Text

Bolten CJ, Kiefer P, Letisse F, Portais JC, Wittmann C: Sampling for metabolome analysis of microorganisms.
Anal Chem 2007, 79:38433849. PubMed Abstract  Publisher Full Text