Computational estimation of tricarboxylic acid cycle fluxes using noisy NMR data from cardiac biopsies
BMC Systems Biology volume 7, Article number: 82 (2013)
The aerobic energy metabolism of cardiac muscle cells is of major importance for the contractile function of the heart. Because energy metabolism is very heterogeneously distributed in heart tissue, especially during coronary disease, a method to quantify metabolic fluxes in small tissue samples is desirable. Taking tissue biopsies after infusion of substrates labeled with stable carbon isotopes makes this possible in animal experiments. However, the appreciable noise level in NMR spectra of extracted tissue samples makes computational estimation of metabolic fluxes challenging and a good method to define confidence regions was not yet available.
Here we present a computational analysis method for nuclear magnetic resonance (NMR) measurements of tricarboxylic acid (TCA) cycle metabolites. The method was validated using measurements on extracts of single tissue biopsies taken from porcine heart in vivo. Isotopic enrichment of glutamate was measured by NMR spectroscopy in tissue samples taken at a single time point after the timed infusion of 13C labeled substrates for the TCA cycle. The NMR intensities for glutamate were analyzed with a computational model describing carbon transitions in the TCA cycle and carbon exchange with amino acids. The model dynamics depended on five flux parameters, which were optimized to fit the NMR measurements. To determine confidence regions for the estimated fluxes, we used the Metropolis-Hastings algorithm for Markov chain Monte Carlo (MCMC) sampling to generate extensive ensembles of feasible flux combinations that describe the data within measurement precision limits. To validate our method, we compared myocardial oxygen consumption calculated from the TCA cycle flux with in vivo blood gas measurements for 38 hearts under several experimental conditions, e.g. during coronary artery narrowing.
Despite the appreciable NMR noise level, the oxygen consumption in the tissue samples, estimated from the NMR spectra, correlates with blood-gas oxygen uptake measurements for the whole heart. The MCMC method provides confidence regions for the estimated metabolic fluxes in single cardiac biopsies, taking the quantified measurement noise level and the nonlinear dependencies between parameters fully into account.
Metabolic fluxes in animal tissues can be identified by measuring the incorporation of stable isotopes in intracellular metabolite pools. To quantify metabolic fluxes, isotope label incorporation is usually measured at several time points , among others in heart tissue –. Heterogeneity of metabolism inside the heart often confounds time series of small tissue samples, therefore a single time point protocol to quantify metabolic fluxes has been developed [5, 6]. Such single time point measurements in individual samples allow to define spatial profiles of metabolic fluxes in heterogeneous organs .
The incorporation of stable isotopes (e.g. 13C) in metabolic intermediates can be detected by nuclear magnetic resonance (NMR) spectroscopy or mass spectrometry (MS). The data is then analyzed with computational methods that require (i) detailed mathematical models of carbon transitions between the metabolites in the system and (ii) sophisticated optimization procedures for estimating the flux parameters. In the past, we have developed a bioinformatics method to estimate metabolic fluxes in aerobic metabolism from very noisy NMR measurements resulting from the Labelling with Isotope for a Pre-Steady-State Snapshot (LIPSSS) protocol . For LIPSSS, isotope labeled substrate for a metabolic pathway is infused for a short, definite period of time, and the metabolism is stopped before a steady state of label incorporation is reached. Finally, pathway metabolites are extracted and measured. Although the original computational analysis method  explores parameter space extensively to avoid local minima, only a rough estimate of parameter confidence regions was obtained by assuming local linearity. Here we introduce a Markov chain Monte Carlo (MCMC) parameter estimation strategy which allows a full description of the confidence regions of the estimated metabolic fluxes, including correlations and nonlinear dependencies between parameter estimates.
Brown et al.  and Gutenkunst et al.  sampled ensembles of parameter sets for systems biology models with MCMC. Correlations between model parameters were taken into account and confidence bounds for parameters and model predictions were defined [9, 10]. Monte Carlo methods have previously been applied to metabolic flux analysis (MFA) in order to handle inaccuracies in data and model . Sensitivity analysis by Monte Carlo sampling is also implemented in a 13C MFA analysis software package . In 13C MFA, MCMC sampling has been used for uncertainty analysis [13, 14], for flux estimation with noisy data , and for in silico experimental design to determine optimal substrate labeling protocols . Antoniewicz et al. proposed a different approach of determining confidence bounds on fluxes by calculating the agreement between model and experiment data as a function of the flux of interest .
We developed and applied an MCMC procedure to estimate the TCA cycle flux, carbon substrate uptake, and oxygen consumption from NMR spectra of 13C enriched glutamate sampled at a single time point. For the computational analysis, we expanded the R-package FluxEs . This analysis was applied to cardiac tissue biopsies flash-frozen 5.5 minutes following 13C acetate infusion in porcine hearts in vivo. The method was validated experimentally for a range of cardiac stress conditions. Our first goal was therefore to determine the uncertainty in the estimation of metabolic flux parameters based on the quantified uncertainty in the NMR measurements and in the prior knowledge. The second goal was to validate the computational estimations in experiments in vivo.
The study was approved by the Advisory Board for the Use of Experimental Animals of the Vrije Universiteit Amsterdam. The procedure is in accordance with the American Physiological Society “Guiding Principles in the Care and Use of Animals,” which state that muscle relaxants may be used in conjunction with drugs known to produce adequate anesthesia.
In this study the metabolic flux in the TCA cycle was measured in tissue biopsies taken from cardiac tissue via the LIPSSS experimental protocol which consists of a brief, timed infusion of 13C labeled acetate in the left anterior descending (LAD) coronary artery of anesthetized pigs . We began with unlabeled acetate which was infused for 30 minutes, in order to establish a stationary metabolic state, followed by [2-13C] acetate for 4 minutes and [1,2-13C] acetate for 1.5 minutes. After exactly 5.5 minutes of 13C enriched acetate infusion, metabolism was arrested by freeze-clamping part of the left ventricular wall of the heart before the isotopic steady state was reached. Biopsies from different regions of this part of the left ventricular wall were cut from the tissue slab after freeze-drying, and divided into approximately nine samples per heart with around 0.1 g dry mass per sample. After extraction with perchloric acid, the 13C NMR multiplets of glutamate were measured. 13C-NMR spectra were obtained at 100.62 MHz and analyzed with the MRUI/AMARES software package (more information about tissue preparation, NMR measurement and the package can be found in reference ).
Up to nine separate multiplet intensities were detected for glutamate. For independent testing of the LIPSSS method and the associated parameter estimation procedures, “gold standard” myocardial oxygen uptake was calculated from blood flow, hemoglobin content and blood-gas measurements taken before and during acetate infusion . Note that these classic oxygen uptake measurements are entirely independent of the LIPSSS method. We analyzed data from LIPSSS samples taken from N = 38 porcine hearts divided into 6 different experimental groups: (i) basal state of the heart (control group, n = 7), two groups with constriction (see below for method) of the coronary vessels to reduce blood flow ((ii) mild stenosis group, n = 7 and (iii) a moderate stenosis group, n = 6), (iv) peripheral venous infusion of dobutamine to induce cardiac stress (dobutamine group, n = 6) or (v) infusion of adenosine for cardiovascular dilatation (adenosine group, n = 4) and (vi) finally, a combination of stenosis and adenosine administration (stenosis + adenosine group, n = 8). In the mild and moderate stenosis groups, LAD blood pressure was adjusted with an occluder to amount to about 70 and 35 mmHg downstream of the occluder, respectively. In the adenosine and stenosis + adenosine groups, adenosine was infused into the LAD at a rate of 100 μg/kg/min. In the stenosis + adenosine group coronary blood pressure was reduced to about 45 mmHg. In the dobutamine group, dobutamine was infused at a rate of 10 μg/kg/min. Note that the dobutamine group initially contained 8 hearts from which two were excluded from further analysis, due to a low mean arterial blood pressure and insufficient NMR signal for parameter estimation (see ), respectively.
Anesthesia and animal experimental procedures
In all groups, sedation was performed with ketamine 15 mg/kg and midazolam 1 mg/kg intramuscularly, and anesthesia was maintained by continuous infusion of sufentanil (4 μg/kg/hr), midazolam (0.5 mg/kg/hr), and pancuronium (0.2 mg/kg/hr). The trachea was intubated and the lungs were ventilated with a mixture of 60% oxygen/40% air. Fluid-filled catheters were introduced and hemodynamic parameters collected as previously described (see ). A continuous infusion of lidocaine was started to help prevent cardiac arrhythmias (9 mg/kg/hr, with an initial bolus injection of 50 mg). Five cm H2O of positive end-expiratory pressure (PEEP) was applied before opening the thorax. The thorax was opened via a midsternal incision and the heart exposed by opening the pericardium. The left hemiazygos vein was tied off to prevent mixing of noncoronary venous blood with coronary venous blood. The LAD was dissected free over a distance of about 2 cm and was catheterized with a 24G catheter. In the stenosis and adenosine + stenosis groups, a custom-made adjustable aluminium occluder was placed around the artery, and LAD pressure was measured.
After finishing instrumentation the animal was allowed to stabilize for at least 15 minutes, the first batch of microspheres (labeled with 141Ce or 103Ru, in random order) was injected into the left atrium for baseline blood flow measurements. The intervention was performed and 30 minutes later a second batch of microspheres was injected for final blood flow measurements. Throughout the procedure hemodynamic data were recorded continuously.
The NMR measured enrichment of glutamate with isotopes was analyzed with a computational model. The model of carbon transitions in the TCA cycle used in this study was described previously in detail [5, 6, 8] and is therefore only described here in brief. The model contains ten metabolite pools, consisting of metabolites which contain 2–6 carbon atoms, and 50 transitions of carbon atoms between the metabolites (Figure 1). Isotopically labeled substrate enters the system via the acetate pool. Acetate is then converted into acetyl coenzyme A (acetyl-CoA), which then enters the TCA cycle. Since acetyl-CoA can also be formed from endogenous unlabeled substrates such as glucose, glycogen, or fatty acids, a diluting pool was introduced to account for dilution of the labeled acetate. The intermediates of the TCA cycle are represented by the 6-carbon metabolite pool labeled as citrate (which also comprises cis-aconitate and isocitrate), α-ketoglutarate, succinate (including succinyl-CoA) and oxaloacetate (representing a 4-carbon metabolite pool which also comprises malate and fumarate). Glutamate and aspartate are amino acids synthesized by transamination from α-ketoglutarate and oxaloacetate, respectively. The replenishment of TCA cycle intermediates was modeled by an anaplerotic influx connected to succinate. Malloy et al. have given detailed descriptions of the equations for anaplerosis [3, 19]. The metabolite concentrations were given as fixed parameters in the calculations: the glutamate pool size was measured by biochemical assay in each sample , because sensitivity analysis showed that results are sensitive to its value. However, the same sensitivity analysis showed that the metabolite pool concentrations of citrate, α-ketoglutarate, oxaloacetate and aspartate had little effect on the results and these concentrations were taken from previous studies [6, 18]. All metabolite concentration parameters were therefore fixed and all flux parameters estimated during the Markov chain Monte Carlo procedure (see below). More information about the model and a listing of all model equations can be found in Additional file 1.
The dynamic behavior of the model is affected by five system parameters (Figure 1). The flux parameters JTCA and Jexch, were expressed in μmol/(min*g dry weight [dw]) and represent reaction fluxes through the TCA cycle and exchange reactions with amino acids, respectively. The dynamics of incorporation of 13C label from acetate into the acetyl-CoA pool depends on transport in the blood vessels, permeation of the cell membrane, the flux of the conversion of acetate into acetyl-CoA, the flux of acetyl-CoA into the TCA cycle and the acetate and acetyl-CoA pool sizes. Fortunately, the time course of incorporation of 13C label into the acetate pool is almost mono-exponential  and can be represented by a single time constant which we term Ttrans. We incorporated this efficient way to represent acetyl-CoA dynamics into our model . The two parameters Pdil and Panap account for the degree of dilution of labeled acetate and the rate of anaplerosis relative to TCA cycle flux, respectively. Both are flux parameters which are expressed as fractions of JTCA. JTCA and Pdil describe energy and substrate turnover which are our targets to measure and are therefore labeled “primary parameters”. On the other hand, Jexch, Ttrans and Panap are constrained during parameter estimation by Bayesian priors (see below) and because they are not our primary target parameters they are termed “auxiliary parameters” which are allowed to vary to determine the uncertainty which they cause in the primary parameters, The LIPSSS estimate for myocardial oxygen consumption is calculated from the primary parameters, (see Eq. 5 below). Note that primary and auxiliary parameters are estimated together in the same procedure.
Matching model simulations to NMR measurements
The computational model described above accounts for all possible carbon isotope labeling states (isotopomers) of each of the metabolites. The system is described by 132 ordinary differential equations (ODEs) to calculate the rate of change of each isotopomer over time. For instance, the metabolite glutamate, which contains 5 carbons, is represented by 25 = 32 ODEs. The isotopomer composition is expressed as fractions of the metabolite concentration of the corresponding pool. Thus, at each time point, the sum of all isotopomer fractions in a pool is 1. All ODEs are then integrated over time to yield the simulated isotopomer fractions. For comparison with the 13C NMR measurements (mexp), simulated NMR multiplet intensities (msim) were calculated from the simulated isotopomer fractions for the time point at which the sample was taken in the experiment . To this end all isotopomers contributing to a particular NMR intensity were added. The simulated multiplet intensities are dependent on the values of the five model parameters. To quantify the agreement between model simulation and experimental data we define a least-squares cost function C as a function of the parameter vector θ, in which the squared residuals for all multiplets are weighted by their standard deviations and summed. Additionally, we include Bayesian prior terms in the cost function which reflect prior knowledge on auxiliary parameter values (see below):
The σi,exp represents the measurement error of the NMR intensity. This cost function is used for the optimization procedures. It is also used as the argument of the normal probability distribution used for the MCMC procedure (see below). The cost function integrates data measured directly in the experiment with literature information incorporated in the priors on parameter values.
Priors on parameter values
The main objective of this study was to estimate JTCA and Pdil, the two primary parameters which define aerobic and substrate metabolism and allow the calculation of oxygen consumption in the sample immediately before metabolic arrest. The three remaining parameters Ttrans, Panap, and Jexch are not our target parameters and cannot be determined with great precision. However, these auxiliary parameters are taken into account to evaluate their effect on the estimation of the primary parameters. To improve the estimation and to help avoid local minima in parameter space with physiologically implausible values of the auxiliary parameters, a priori information for such parameters (θi) can be directly entered into the cost function by adding a prior term to the cost function in Equation 1 for the deviation from a certain reference value θi*
where is the standard deviation for the auxiliary parameter in log-space. The advantage of logarithmic parameters is that the parameter values with a Gaussian prior distribution are positive and dimensionless. Note that the prior probability in Equation 2 does not include the normalization factor for the lognormal distribution of . Normalization was not necessary because our method applied the Metropolis-Hastings algorithm which uses the ratios of probabilities.
In previous studies, the value of Ttrans had been estimated to be 0.202 min which is compatible with the time constant of the enrichment of acetyl-CoA with radioactive label [6, 8]. We constrain Ttrans around this prior value with set to 0.336, a high value used in a previous study for unreported experimental errors . This is slightly higher than the value for the standard deviation of these parameters determined in simulations by Binsl et al. . The central 95% region of the prior for Ttrans lies between 0.202/1.96 and 0.202 *1.96 min, since = 0.336 = 1/4*(ln(θi * 1.96)-ln(θi/1.96)), (see ).
The accurate quantification of the exchange flux Jexch between α-amino and α-keto acids was found to be challenging [2, 22]. A previous analysis of the model used in this study revealed a low sensitivity of estimations of JTCA to variations of Jexch in the physiological range from 5–60 μmol/(min*gdw) . Reported values of exchange flux in the literature vary substantially. Some report high values for the exchange flux (e.g. 13-fold the flux of JTCA). Several other studies report Jexch to be approximately equal to JTCA[23, 24]. To address this issue, we set a prior on Jexch relative to the value of JTCA. Instead of calculating the prior cost directly from Jexch, it is therefore determined by entering the ratio θi = Jexch/ JTCA into Equation 2. The reference value θi* for the ratio is set to 1, based on values for Jexch/ JTCA reported by Nuutinen et al.  and Yu et al. .
Because of the large spread of values found in the literature (see above), we assumed a high standard deviation for the ratio Jexch/ JTCA and set to 1/4*(ln(θi * 15)-ln(θi/15)) = 1.345, with θi = 1. It is thereby ensured that Jexch lies with 95% probability between JTCA/15 and JTCA*15.
For the parameter Panap, the anaplerotic flux relative to the TCA cycle flux, most of the values found in literature were smaller than 1 and the highest experimental value found was reported to be 1 ± 0.3 [25, 26]. Hence the prior cost for Panap was set to be uniform for values of Panap between 0 and 1 combined with a half-normal distribution which had a standard deviation of 0.3 taken from Lloyd et al.  for the values above 1:
with and . The normalization constants c1 and c2 ensure that the probability density function of the prior is continuous and that its integral is equal to one. denotes the normal distribution. The probability density functions for prior(Ttrans), prior(Jexch), and prior(Panap) are shown in Figure 2 (solid lines).
Parameter estimation and sampling of parameter ensembles
In biological models, usually many different combinations of parameters can describe the experimental data . To address this, we decided to not merely rely on a single best-fit of the model parameters to the NMR data for fixed values of the auxiliary parameters, but instead, we systematically generated ensembles of model parameters that fit the data with reasonable precision. This approach clarifies how well the primary parameters are defined by the data despite uncertainty in the NMR intensities and auxiliary parameters. Through the use of an MCMC approach, confidence bounds can be set on the estimated parameter values. Sampling is based on Bayesian inference of a posterior parameter distribution
where Pr(D|θ) is the probability of a parameter vector θ to describe the given data D and Pr(θ) is the prior probability of the parameters (see above). The right-hand side of Equation 4 is equal to e− C(θ) where the cost function of Equation 1 (which includes the priors of Equations 2 and 3) is used. The probability functions were not all normalized because this was not necessary for the MCMC procedure which relies on the ratios of probabilities rather than absolute values. Note that the cost function (Equations 1, 2, 3) forms the basis of a probability function (Equation 4) that defined the ensemble of estimated parameter values.
In order to estimate the model parameters and to quantify the uncertainty of the estimated values, we sampled an ensemble of parameter sets which could describe the available NMR data by performing a random walk through the parameter space through the application of the Metropolis-Hastings algorithm. The starting point of the random walk was an optimized set of parameters, which had been obtained by a grid optimization strategy introduced by Binsl et al. . The grid optimization was designed to cope with a shallow basin shaped by the cost function in order to avoid optimization towards local minima. The procedure covered the 5-dimensional parameter space by a grid so as to find the best starting point for optimization. The second phase of optimization starting at this grid point was then performed using the Nelder & Mead simplex algorithm, and in the third phase we used the Metropolis-Hastings algorithm to sample a parameter ensemble with its probability density proportional to a probability function based on the cost function C(θ) of Equation 1 entered in Equation 4.
Quality criteria for flux estimations in NMR samples
In many of the available in vivo samples, NMR peak intensities are low and often below the threshold of observability, i.e. often six or seven of the nine multiplets of glutamate are not discernible from noise and were assigned an intensity of zero. In some of these low intensity samples, Monte Carlo sampling leads to very large ensemble standard deviations of the estimated primary parameters. We excluded such samples which did not yield reliable estimates for JTCA. The exclusion criterion was that the standard deviation of JTCA in the posterior parameter ensemble exceeded 10 μmol/(min*gdw).
Software package FluxEs
The analysis was performed using the R package FluxEs introduced by Binsl et al. . In order to process parameter ensembles, a Monte Carlo module was added to the software. This module uses the AMCMC algorithm implemented within the package spBayes [27, 28]. The AMCMC algorithm is a Metropolis-Hastings variant which automatically adapts the proposal step size for the sampled parameters in the random walk. This leads to quicker convergence to a posterior distribution. For the primary parameters, the time constant of the autocorrelation function of the sampled ensemble was calculated in order to inspect whether the algorithm converged to a stationary distribution. For samples with a high autocorrelation time in the primary parameters, we visually inspected the parameter trace.
A single model simulation run takes approximately 0.26 seconds on a computer with 2.26 GHz clock frequency. The grid optimization for a single sample took on average 115 minutes, the subsequent sampling with the adaptive Metropolis-Hastings algorithm took on average 540 minutes per sample.
The calculations for all samples were performed in parallel on the Lisa computer cluster system at SARA Computing and Networking Services (http://www.sara.nl). All code required for the analysis and part of the experimental data can be found in Additional file 2.
We estimated the TCA cycle flux from the NMR peaks of glutamate for 347 tissue samples from 38 hearts. Applying the exclusion criterion described above we removed 85 low-quality samples - leaving 262 samples for further analysis. For each sample, an ensemble of 35,000 parameter sets was generated with the Metropolis-Hastings algorithm. Although convergence was not the first criterion for sample rejection, all ensemble estimates with a high autocorrelation time constant were rejected according to the quality criterion.
An example of a parameter ensemble for one single sample of the control group is given in Figure 2. For Ttrans, Jexch, and Panap, the probability density functions of the prior distributions are plotted together with the histograms of the posterior distributions. The posterior distributions for these auxiliary parameters are very broad and relatively close to their corresponding prior distributions. In this way the MCMC ensemble method allowed defining the uncertainty in the primary parameters taking into account the large spread in auxiliary parameters. Despite the broad distribution of the auxiliary parameters, the estimates for JTCA and Pdil form relatively well-defined peaks and their standard deviations are relatively low.
For the primary parameters we can thus provide point estimates for each sample. To determine which measure best reflects the true value of a parameter, we conducted a simulation experiment in which multiple sets of artificial NMR multiplets were generated by model simulation and subsequent addition of Gaussian random measurement noise. The parameters were then re-estimated and we compared the estimates from the best fit after grid optimization (i.e. the fit with the lowest cost function value and therefore the highest likelihood, see Equation 4) and the mean, median, and mode of the Monte Carlo ensemble with the “true” parameter values from the initial simulation. Regarding the primary parameters, the best fit gave the most reliable point estimate. Below, we therefore report the best fit values for the primary parameters.
Validation by estimation of myocardial oxygen consumption
In order to validate our flux estimation method we compared the LIPSSS estimated myocardial oxygen consumption (MVO2, expressed in μmol/(min*gdw)) with independent “gold standard” measurements. The “gold standard” was determined by blood-gas oxygen and blood flow measurements and the LIPSSS estimated oxygen consumption was calculated from the parameter estimates of the model . The MVO2 for a single sample is determined from the primary LIPSSS flux parameters by stoichiometric biochemical relations and can be calculated as follows [8, 29]:
The MVO2 determined from blood-gas measurements reflects the oxygen consumption of the entire heart. When averaging the samples taken for LIPSSS measurements to estimate oxygen consumption for the entire heart (), individual sample sizes were taken into account. As in Binsl et al., the contributions of the individual samples were weighted by the dry weight wsample for each sample .
For all six experimental groups, the comparison of MVO2 estimated with the LIPSSS method (from the model parameters Pdil and JTCA) with the “gold standard” oxygen measurements is shown in Figure 3. One heart from the stenosis + adenosine group was excluded from the analysis since none of its samples satisfied the quality criterion.
For all groups, LIPSSS MVO2 correlated with blood-gas MVO2 relatively well. For the control group oxygen consumption measured by the two methods corresponded, but for the ischemic conditions (stenosis with and without adenosine), oxygen consumption tended to be lower for the LIPSSS method. We calculated Pearson correlation coefficients of 0.49 for control (n = 7, p = 0.26), 0.69 for mild stenosis (n = 7, p = 0.09), 0.66 for moderate stenosis (n = 6, p = 0.15), 0.99 for dobutamine (n = 6, p = 0.0003), 0.71 for adenosine (n = 4, p = 0.29), and 0.87 for the stenosis + adenosine group (n = 7, p = 0.01). The Pearson correlation for all groups combined was 0.85 (n = 37, p < 10-10). The dobutamine group showed higher oxygen consumption than the other groups reflecting the increased cardiac work load. It is important to note that the small tissue biopsies used in the LIPSSS experiment only covered a relatively small cardiac region, in contrast to the physiological blood-gas measurements which covered the entire left ventricle. Furthermore, the estimation of MVO2 from the parameters of the TCA cycle model only reflects myocardial oxygen consumption linked to the TCA cycle flux, disregarding other oxygen consuming reactions which were covered by the blood-gas measurements. The oxygen consumption measurements in a small ischemic region dependent on a constricted coronary artery would be very difficult to obtain with classic blood-gas measurements.
Estimation of TCA cycle fluxes
LIPSSS-based estimates for the primary model parameters under all experimental conditions are shown in Figure 4. Estimates for JTCA in the control group averaged 7.04 ± 0.79 (mean ± SEM) μmol/(min*gdw). For mild and moderate constriction of the coronary vessels, we estimated JTCA to be 4.12 ± 0.49 and 2.99 ± 0.36 μmol/(min*gdw), respectively.
Dobutamine infusion, which stimulates cardiac contraction, leads to a high average JTCA estimate of 11.18 ± 1.31 μmol/(min*gdw) of tissue. Estimations for the adenosine group show no difference with the control condition. The TCA cycle flux in the stenosis + adenosine group is in between the mild and moderate stenosis condition.
The relative contribution to the TCA cycle flux of substrates other than labeled acetate, i.e. Pdil is higher in all experimental groups compared with the baseline condition (see Figure 4). Low fractional acetate usage, i.e. high dilution has been previously documented in experiments with dobutamine .
The estimations of Ttrans for all the groups did not differ substantially from the prior value of 0.202 minutes (data not shown). Ensembles for the auxiliary parameters Jexch and Panap show large standard deviations. This indicates that these parameters cannot be estimated properly from the NMR data. Indeed, the experimental protocol was optimized to estimate the primary parameters, disregarding the auxiliary parameters. Nevertheless the effect of the potential spread in these auxiliary parameters on the uncertainty limits of the primary parameters was taken into account. Estimations of the auxiliary model parameters are described in supplemental file 3.
The fluxes of biochemical reactions linked to cardiac energy metabolism are of significant interest. Here we investigated a computational method to quantify fluxes in the TCA cycle using NMR data from 13C labeling experiments in porcine hearts. We took measurement error in the data and uncertainty of model parameters directly into account. To test the method, distinct 13C labeling patterns (isotopomers) in glutamate were measured under six different cardiac stress and control conditions. The data were analyzed with a detailed model of carbon transitions in the TCA cycle and two primary flux parameters of interest (reflecting total aerobic metabolism and uptake of the labeled substrate) were estimated. Possible variation in three auxiliary parameters, taken from experimental literature was included in the application of Bayesian priors. To define the uncertainty in estimated flux parameters from measurement error and uncertainty in prior knowledge, we used an MCMC method. As a result, we were able to derive estimates for the TCA cycle fluxes under various experimental conditions despite the high noise level in the available NMR data. For validation, we compared blood-gas measurements of myocardial oxygen consumption with oxygen consumption calculated from our own parameter estimates. The oxygen consumption estimated with our model correlated with the classic physiological measurements for the whole heart (Figure 3).
However, because the LIPSSS parameter estimates relied on small samples obtained from the heart while the blood gas measurements represented the oxygen consumption for the whole heart, the LIPSSS estimates are expected to deviate from the whole heart measurement. The deviation may have a random component because of the limited tissue sample size, and a systematic component because of functional differences between regions in the heart. The random component is expected because heterogeneity of blood flow and metabolism has been measured in heart muscle [18, 31]. A systematic component is expected especially in the stenosis groups, because the LIPSSS NMR measurements are taken from regions with lower oxygen consumption caused by low perfusion. However, it should be noted that this reasoning is incomplete because the blood gas estimation of oxygen consumption takes the local blood flow measured in the stenosed region into account. Nevertheless, systematic differences between the small region and the average for the whole heart may contribute to the deviation from the line of identity (see Figure 3) at low oxygen consumptions.
Additional physiological measurements of oxygen consumption and metabolic fluxes, independent from the stable isotope labeling experiments, are desirable for further validation of our method. Regional rates of oxygen consumption can be measured by measuring oxygen content in small veins  with a spectroscopical method in frozen tissue. The latter method is difficult and its validation has been criticized. A further method is the simultaneous determination of myocardial perfusion and oxygen content in small regions of the heart . Oxygen consumption can also be measured using PET and TCA cycle fluxes using in vivo NMR (e.g. ). However, these methods mostly have very limited spatial resolution  and were in turn subject to rather limited validation themselves. The difficulty in measuring local energy metabolic flux provided motivation to develop our present method in the first place. Despite the limited possibilities, further validation of the LIPSSS method in the future is desirable.
Part of the dataset used here, namely the control and dobutamine group, had been analyzed in a previous study . The estimates of Binsl et al.  relied on prior information on the model parameters Ttrans and Panap. The latter parameter, describing anaplerosis relative to the TCA cycle flux was constrained to 6 ± 3% of JTCA, based on information from literature studies on isolated hearts. The latter studies however, only accounted for anaplerosis from either propionate  or from pyruvate [35, 36]. It has been suggested that relative anaplerosis is often underestimated by conventional approaches, including isotopomer analysis or fractional enrichments of carbons in glutamate . Tracer experiments also exist using 13C labeled propionate that report the relative anaplerotic flux in rat hearts to be much higher than 6%, e.g. 16%  or 29% . Higher relative anaplerotic fluxes were reported during low flow ischemia, reaching 100%  and 35% . Higher values have also been reported for hypertrophy . Although our estimates for the parameter Panap in the present study have a relatively low precision, they suggest the possibility that in porcine heart anaplerotic flux in vivo is relatively high in contrast to low values often found in isolated hearts (see Additional file 3).
Since three different stenosis conditions were included in the present study, we chose a less constraining Bayesian prior on the parameter Panap which covered a broad range. It is important to note that the Bayesian priors were the same for the analysis of NMR data from all experimental conditions. Despite the use of different choices of priors on the parameters, and although a higher anaplerotic flux was estimated (see Additional file 3), our present estimates for fluxes in the control and dobutamine groups did not differ much from the previous estimates of Binsl et al. . Our estimates for cardiac ischemia induced by coronary stenosis show that the TCA cycle flux decreases whilst the relative anaplerosis increases (see Figure 4 and Additional file 3) which is compatible with existing literature (see references cited above).
Due to the high velocity of the exchange reactions between α-amino and α-keto acids, the determination of Jexch using tracer experiments is expected to be practically infeasible . Because of the uncertainty on Jexch, we decided to evaluate the effect of variation in Jexch. Values for Jexch/JTCA found in literature vary between 0.2 and 13 [20, 22], but are often around 1 in the heart , in contrast to the very high Jexch/JTCA reported for the human brain . Initial estimations of Jexch in our data showed that, particularly in samples with low NMR peak intensity, the simulated isotope enrichment was not very sensitive to Jexch. Rather than constraining Jexch around an absolute value, we chose to set a Bayesian prior relative to JTCA. The standard deviation of the prior was set to a very high value, reflecting the high variability of Jexch/JTCA measurements found in the literature. Jexch/JTCA estimated with our method ranged from 0.74 (median dobutamine group) to 1.75 (median control group). Weiss et al. reported a decreased absolute exchange flux compared with control conditions during post-ischemic reperfusion in rat hearts . A decrease in Jexch during stenosis was estimated in the present study (see Additional file 3).
Literature information on parameter values was incorporated into the analysis as Bayesian priors because of the high noise level in the NMR data. Without using prior information, flux parameters sometimes reach physiologically infeasible regions in parameter space. We investigated the sensitivity of our estimates of the primary parameters to the priors for the auxiliary parameters by re-performing the analysis with doubled prior standard deviations in equations 2 and 3. The estimate for parameter Pdil is rather insensitive to changes in the prior standard deviation (absolute difference in the estimated value averaged over all groups is 4.4±4.0%) while estimates of JTCA are more sensitive to alterations in the priors on auxiliary parameters (average absolute difference 20±21%). Especially in the moderate stenosis group, for which the NMR signals are on average very low, many estimates fail to meet the quality criteria if the standard deviation for all three priors simultaneously was made twice as large. This shows that the estimate of JTCA is sensitive to the prior. However, Bayesian priors are necessary to constrain the estimates within reasonable ranges. It is therefore important to emphasize that the prior values and their standard deviations are not arbitrarily chosen. The prior distributions of Panap and Jexch are based on experimental data [20, 23]– and were given large standard deviations. The prior on Ttrans is based on previous estimates [6, 8] and its standard deviation allows for a broad range. We therefore argue that although constraining parameters in this study was necessary due to the high noise in the data, our framework still allowed to define reasonable point estimates of flux parameters and additionally to define the variability in parameter estimates taking reasonable, sometimes deliberately high, values for the uncertainty of auxiliary parameters into account.
Parametric sensitivity analysis is commonly applied in systems biology . In this investigation, we chose an approach that explored the multidimensional space around a set of best-fit parameters using a random walk with the Metropolis-Hastings algorithm [9, 10, 21]. The advantage of this method is that it takes into account possible correlations and nonlinear dependencies between the model parameters. Antoniewiecz et al. approached the problem of defining confidence regions for flux estimates by minimizing a sum of squared residuals objective function as a function of the flux value . In their approach, the confidence interval for a flux of interest is derived by setting the flux constant while optimizing all remaining fluxes in the system. This step is repeated for a range of fixed flux values until the objective function value exceeds a predefined confidence limit. The advantage of the MCMC approach to determine confidence regions is that it takes all possible correlations between the fluxes into account, since no flux parameter is fixed during the MCMC sampling.
The challenge in analyzing the data in this study was the high noise level. Up to seven of the nine measured multiplet intensities could sometimes not be detected. Ensemble modeling proved to be a feasible method to separate samples with flux parameters that could be estimated from samples with poor information on the fluxes in the system. This ensemble approach made it possible to identify 262 out of 347 samples that gave useful estimates for the primary parameters. The quality selection of the samples allowed us to use the best-fit parameters from each sample as a point estimate for the primary parameters. The MCMC approach allowed us to define confidence bounds on all estimated parameter values taking their correlations into account. This is a significant advantage compared with previous approaches, where linearized or analytical methods were used to calculate errors on estimated model parameters [5, 6, 8].
Adding the Monte Carlo ensemble sampling to the LIPSSS framework enables us to estimate the confidence regions of flux parameters in a single sample. The small size of the tissue samples makes it feasible to identify the spatial variation of flux parameters expected because of the known heterogeneity in the tissue. The physiological meaning of our measurements of heterogeneity in metabolism in heart muscle will be addressed in future studies.
In this study we improved the LIPSSS method in order to quantify metabolic fluxes using stable isotope labeling integrated with mathematical models of carbon transitions: auxiliary information was taken into account in the form of Bayesian priors and emphasis was placed on the uncertainty analysis of the estimated flux parameters. The method was used to quantify TCA cycle fluxes from noisy NMR measurements in porcine hearts under different physiological conditions. Two important metabolic fluxes could be determined in single biopsies taken during animal experiments and confidence regions could be calculated for single samples.
Sauer U: Metabolic networks in motion: 13C-based flux analysis. Mol Syst Biol 2006, 2: 62.
Chance EM, Seeholzer SH, Kobayashi K, Williamson JR: Mathematical analysis of isotope labeling in the citric acid cycle with applications to 13C NMR studies in perfused rat hearts. J Biol Chem 1983,258(22):13785-13794.
Malloy CR, Sherry AD, Jeffrey FM: Analysis of tricarboxylic acid cycle of the heart using 13C isotope isomers. Am J Physiol Heart Circ Physiol 1990,259(3 Pt 2):H987-H995.
Schroeder MA, Atherton HJ, Ball DR, Cole MA, Heather LC, Griffin JL, Clarke K, Radda GK, Tyler DJ: Real-time assessment of Krebs cycle metabolism using hyperpolarized 13C magnetic resonance spectroscopy. FASEB J 2009,23(8):2529-2538. 10.1096/fj.09-129171
van Beek JH, Csont T, de Kanter FJ, Bussemaker J: Simple model analysis of 13C NMR spectra to measure oxygen consumption using frozen tissue samples. Adv Exp Med Biol 1998, 454: 475-485. 10.1007/978-1-4615-4863-8_58
van Beek JH, van Mil HG, King RB, de Kanter FJ, Alders DJ, Bussemaker J: A (13)C NMR double-labeling method to quantitate local myocardial O(2) consumption using frozen tissue samples. Am J Physiol 1999,277(4 Pt 2):H1630-H1640.
Alders DJC, Cornelussen RN, Prinzen FW, Specht PAC, Noble MIM, Drake-Holland AJ, de Kanter FJJ, van Beek JHGM: Regional sympathetic denervation affects the relation between canine local myocardial blood flow and oxygen consumption. Exp Physiol 2007,92(3):541-548. 10.1113/expphysiol.2006.036228
Binsl TW, Alders DJC, Heringa J, Groeneveld ABJ, van Beek JHGM: Computational quantification of metabolic fluxes from a single isotope snapshot: application to an animal biopsy. Bioinformatics 2010,26(5):653-660. 10.1093/bioinformatics/btq018
Brown KS, Hill CC, Calero GA, Myers CR, Lee KH, Sethna JP, Cerione RA: The statistical mechanics of complex signaling networks: nerve growth factor signaling. Phys Biol 2004,1(3–4):184-195.
Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP: Universally sloppy parameter sensitivities in systems biology models. PLoS Comput Biol 2007,3(10):1871-1878.
Schellenberger J, Palsson BØ: Use of randomized sampling for analysis of metabolic networks. J Biol Chem 2009,284(9):5457-5461.
Weitzel M, Nöh K, Dalman T, Niedenführ S, Stute B, Wiechert W: 13CFLUX2 - High-Performance Software Suite for 13C-Metabolic Flux Analysis. Bioinformatics 2013,29(1):43-45.
Kadirkamanathan V, Yang J, Billings SA, Wright PC: Markov Chain Monte Carlo Algorithm based metabolic flux distribution analysis on Corynebacterium glutamicum. Bioinformatics 2006,22(21):2681-2687. 10.1093/bioinformatics/btl445
Quek L-E, Wittmann C, Nielsen LK, Krömer JO: OpenFLUX: efficient modelling software for 13C-based metabolic flux analysis. Microb Cell Fact 2009, 8: 25. 10.1186/1475-2859-8-25
Yang J, Wongsa S, Kadirkamanathan V, Billings SA, Wright PC: Metabolic flux distribution analysis by 13C-tracer experiments using the Markov chain-Monte Carlo method. Biochem Soc Trans 2005,33(Pt 6):1421-1422.
Schellenberger J, Zielinski DC, Choi W, Madireddi S, Portnoy V, Scott DA, Reed JL, Osterman AL, Palsson BO: Predicting outcomes of steady-state 13C isotope tracing experiments with Monte Carlo sampling. BMC Syst Biol 2012, 6: 9. 10.1186/1752-0509-6-9
Antoniewicz MR, Kelleher JK, Stephanopoulos G: Determination of confidence intervals of metabolic fluxes estimated from stable isotope measurements. Metab Eng 2006,8(4):324-337. 10.1016/j.ymben.2006.01.004
Alders DJC, Groeneveld ABJ, de Kanter FJJ, van Beek JHGM: Myocardial O2 consumption in porcine left ventricle is heterogeneously distributed in parallel to heterogeneous O2 delivery. Am J Physiol Heart Circ Physiol 2004,287(3):H1353-H1361. 10.1152/ajpheart.00338.2003
Malloy CR, Jones JG, Jeffrey FM, Jessen ME, Sherry AD: Contribution of various substrates to total citric acid cycle flux and anaplerosis as determined by 13C isotopomer analysis and O2 consumption in the heart. Magma (New York, N.Y.) 1996, 4: 35-46.
Randle PJ, England PJ, Denton RM: Control of the tricarboxylate cycle and its interactions with glycolysis during acetate utilization in rat heart. Biochem J 1970,117(4):677-695.
Hettling H, van Beek JHGM: Analyzing the Functional Properties of the Creatine Kinase System with Multiscale “Sloppy” Modeling. PLoS Comput Biol 2011,7(8):e1002130. 10.1371/journal.pcbi.1002130
Jeffrey FMH, Reshetov A, Storey CJ, Carvalho RA, Sherry AD, Malloy CR: Use of a single 13C NMR resonance of glutamate for measuring oxygen consumption in tissue. Am J Physiol 1999,277(6 Pt 1):E1111-E1121.
Nuutinen EM, Peuhkurinen KJ, Pietiläinen EP, Hiltunen JK, Hassinen IE: Elimination and replenishment of tricarboxylic acid-cycle intermediates in myocardium. Biochem J 1981,194(3):867-875.
Yu X, White LT, Doumen C, Damico LA, LaNoue KF, Alpert NM, Lewandowski ED: Kinetic analysis of dynamic 13C NMR spectra: metabolic flux, regulation, and compartmentation in hearts. Biophys J 1995,69(5):2090-2102. 10.1016/S0006-3495(95)80080-9
Weiss RG, Chacko VP, Glickson JD, Gerstenblith G: Comparative 13C and 31P NMR assessment of altered metabolism during graded reductions in coronary flow in intact rat hearts. Proc Natl Acad Sci U S A 1989,86(16):6426-6430. 10.1073/pnas.86.16.6426
Lloyd SG, Wang P, Zeng H, Chatham JC: Impact of low-flow ischemia on substrate oxidation and glycolysis in the isolated perfused rat heart. Am J Physiol Heart Circ Physiol 2004,287(1):H351-H362. 10.1152/ajpheart.00983.2003
Rosenthal JS: AMCMC: An R interface for adaptive MCMC. Comput Stat Data Anal 2007, 51: 5467-5470. 10.1016/j.csda.2007.02.021
Finley AO, Banerjee S, Carlin BP: spBayes: An R Package for Univariate and Multivariate Hierarchical Point-referenced Spatial Models. J Stat Softw 2007,19(4):1-24.
Chatham JC, Forder JR, Glickson JD, Chance EM: Calculation of absolute metabolic flux and the elucidation of the pathways of glutamate labeling in perfused rat heart by 13C NMR spectroscopy and nonlinear least squares analysis. J Biol Chem 1995,270(14):7999-8008. 10.1074/jbc.270.14.7999
Robitaille PM, Rath DP, Abduljalil AM, O’Donnell JM, Jiang Z, Zhang H, Hamlin RL: Dynamic 13C NMR analysis of oxidative metabolism in the in vivo canine myocardium. J Biol Chem 1993,268(35):26296-26301.
Weiss HR, Sinha AK: Regional oxygen saturation of small arteries and veins in the canine myocardium. Circ Res 1978, 42: 119-126. 10.1161/01.RES.42.1.119
Reeder SB, Holmes AA, McVeigh ER, Forder JR: Simultaneous noninvasive determination of regional myocardial perfusion and oxygen content in rabbits: toward direct measurement of myocardial oxygen consumption at MR imaging. Radiology 1999, 212: 739-747.
Kofoed K, Hansen P, Holm S, Hove J, Chen K, Jin W, Jensen M, Iida H, Hesse B, Svendsen J: Regional myocardial oxygen consumption estimated by carbon-11 acetate and positron emission tomography before and after repetitive ischemia. J Nucl Cardiol 2000, 7: 228-234. 10.1016/S1071-3581(00)70011-0
Martini WZ, Stanley WC, Huang H, Rosiers CD, Hoppel CL, Brunengraber H: Quantitative assessment of anaplerosis from propionate in pig heart in vivo. Am J Physiol Endocrinol Metab 2003,284(2):E351-E356.
Panchal AR, Comte B, Huang H, Kerwin T, Darvish A, des Rosiers C, Brunengraber H, Stanley WC: Partitioning of pyruvate between oxidation and anaplerosis in swine hearts. Am J Physiol Heart Circ Physiol 2000,279(5):H2390-H2398.
Panchal AR, Comte B, Huang H, Dudar B, Roth B, Chandler M, Des Rosiers C, Brunengraber H, Stanley WC: Acute hibernation decreases myocardial pyruvate carboxylation and citrate release. Am J Physiol Heart Circ Physiol 2001,281(4):H1613-H1620.
Cohen DM, Bergman RN: Improved estimation of anaplerosis in heart using 13C NMR. Am J Physiol Endocrinol Metab 1997,273(6 Pt 1):E1228-E1242.
Kasumov T, Cendrowski AV, David F, Jobbins KA, Anderson VE, Brunengraber H: Mass isotopomer study of anaplerosis from propionate in the perfused rat heart. Arch Biochem Biophys 2007,463(1):110-117. 10.1016/j.abb.2007.02.022
Sorokina N, O’Donnell JM, McKinney RD, Pound KM, Woldegiorgis G, LaNoue KF, Ballal K, Taegtmeyer H, Buttrick PM, Lewandowski ED: Recruitment of compensatory pathways to sustain oxidative flux with reduced carnitine palmitoyltransferase I activity characterizes inefficiency in energy metabolism in hypertrophied hearts. Circulation 2007,115(15):2033-2041. 10.1161/CIRCULATIONAHA.106.668665
Mason GF, Gruetter R, Rothman DL, Behar KL, Shulman RG, Novotny EJ: Simultaneous determination of the rates of the TCA cycle, glucose utilization, alpha-ketoglutarate/glutamate exchange, and glutamine synthesis in human brain by NMR. J Cereb Blood Flow Metab 1995,15(1):12-25. 10.1038/jcbfm.1995.2
Weiss RG, Kalil-Filho R, Herskowitz A, Chacko VP, Litt M, Stern MD, Gerstenblith G: Tricarboxylic acid cycle activity in postischemic rat hearts. Circulation 1993,87(1):270-282. 10.1161/01.CIR.87.1.270
Perumal TM, Gunawan R: Understanding dynamics using sensitivity analysis: caveat and solution. BMC Syst Biol 2011, 5: 41. 10.1186/1752-0509-5-41
We thank SARA Computing and Networking Services (http://www.sara.nl) for their support in using the Lisa Computer Cluster. We thank Frans de Kanter for assistance with the NMR measurements. We also thank Erik van Dijk for scientific advice and Irisa Ono for editorial assistance.
Hannes Hettling was supported by the Centre for Medical Systems Biology (CMSB), the Netherlands Consortium for Systems Biology (NCSB) and the Netherlands Bioinformatics Centre (NBIC) which are centres of excellence supported by the Dutch government via the Netherlands Genomics Initiative. Parts of the open access costs for this article were paid by NWO, the Netherlands Organisation for Scientific Research.
The authors declare that they have no competing interests.
DJA, ABJG and JHGMvB designed and conducted the animal experiments. HH, TWB and JHGMvB designed and conceived the in silico experiments, HH and TWB performed the in silico experiments. HH, DJA, JH and JHGMvB analyzed the data. HH and JHGMvB wrote the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 1:Model equations. In this supplemental text, we give a detailed description of the computational model and list all model ODEs. (PDF 196 KB)
Additional file 2:Code and data. In this supplemental file we provide all R code and part of the experimental data used to produce the results of this study. The zip file contains a file README.txt which describes all code and data. (ZIP 59 KB)
Additional file 3:Estimation of auxiliary model parameters. In this supplemental text, the results of estimating the auxiliary model parameters Panap and Jexch are presented and discussed. (PDF 146 KB)
About this article
Cite this article
Hettling, H., Alders, D.J.C., Heringa, J. et al. Computational estimation of tricarboxylic acid cycle fluxes using noisy NMR data from cardiac biopsies. BMC Syst Biol 7, 82 (2013). https://doi.org/10.1186/1752-0509-7-82