Refined Analysis of Brain Energy Metabolism Using In Vivo Dynamic Enrichment of 13C Multiplets

Carbon-13 nuclear magnetic resonance spectroscopy in combination with the infusion of 13C-labeled precursors is a unique approach to study in vivo brain energy metabolism. Incorporating the maximum information available from in vivo localized 13C spectra is of importance to get broader knowledge on cerebral metabolic pathways. Metabolic rates can be quantitatively determined from the rate of 13C incorporation into amino acid neurotransmitters such as glutamate and glutamine using suitable mathematical models. The time course of multiplets arising from 13C-13C coupling between adjacent carbon atoms was expected to provide additional information for metabolic modeling leading to potential improvements in the estimation of metabolic parameters. The aim of the present study was to extend two-compartment neuronal/glial modeling to include dynamics of 13C isotopomers available from fine structure multiplets in 13C spectra of glutamate and glutamine measured in vivo in rats brain at 14.1 T, termed bonded cumomer approach. Incorporating the labeling time courses of 13C multiplets of glutamate and glutamine resulted in elevated precision of the estimated fluxes in rat brain as well as reduced correlations between them.

Metabolic modeling of Carbon labeling time courses of amino acids in brain allows quantitative measurements of metabolic fluxes in vivo (Gruetter et al., 2001;Rothman et al., 2003;Oz et al., 2004;Lanz et al., 2012). In the past decade, questions have been raised concerning the reliability of some estimated fluxes in brain (Shestov et al., 2007;Uffmann and Gruetter, 2007;Shen et al., 2009) and many efforts have been undertaken to improve their accuracy through adaptation of the metabolic models used to analyze the 13 C turnover curves. Monte-Carlo simulations showed that an increased number of measured time courses of 13 C positions of Glu and Gln as input for the metabolic modeling analysis resulted in an improved precision and reliability of brain metabolic fluxes (Shestov et al., 2007). However, the number of measured 13 C labeling time courses of metabolites relies on the NMR detection method, the chosen magnetic field strength and pulse sequence in different studies, and varied from two to seven 13 C uptake curves (Mason et al., 1995;Gruetter et al., 2001;Choi et al., 2002;de Graaf et al., 2004;Xin et al., 2010Xin et al., , 2015. Due to the difficulty of measuring 13 C multiplets of Glu and Gln in vivo with high temporal resolution, derivation of metabolic fluxes has been performed almost exclusively using dynamic positional enrichment (e.g., de Graaf et al., 2003;Gruetter et al., 2003;Henry et al., 2006;Duarte et al., 2011) which is the time course of total 13 C enrichment of each carbon position of each metabolite.
The measurement of labeling time courses of 13 C isotopomers noninvasively in the rat brain has been already reported by Henry et al. (2003a). A first analysis of in vivo measured 13 C dynamic labeling of isotopomers in a heart metabolic study led to improved precision of the estimated tricarboxylic acid (TCA) cycle flux (Jeffrey et al. 1999). In a study of brain metabolism by the same group, the 13 C time courses of the singlet and doublet of cerebral Glu labeled at position C4 was included in the metabolic analysis and resulted in a poorer determination of glutamate/alpha-ketoglutarate and oxaloacetate/aspartate exchange rate, called V x (Jeffrey et al., 2013), probably reflecting the incomplete description of the Glu C4 singlet turnover by the considered metabolic model.
Recently, a new 13 C metabolic modeling approach using the concept of bonded cumomers was proposed, taking full advantage of the available time courses of 13 C multiplets of Glu and Gln in the metabolic model (Shestov et al., 2012). In simulations, the bonded cumomer approach led to better precision than the positional approach in the determination of brain metabolic fluxes using artificial 13 C time courses of Glu and Gln multiplets generated with the two-compartmental description of brain energy metabolism (Gruetter et al., 2001) for different infused 13 C-labeled substrates.
Including the individual time courses of 13 C multiplets in the metabolic modeling process is expected to provide increased number and precision of the metabolic fluxes of interest. Therefore, the aim of this study was to incorporate the 13 C dynamic labeling of observable multiplets of Glu and Gln measured with in vivo 13 C MRS in the rat brain at 14.1 T to estimate brain metabolic fluxes in the case of a neuronal-glial compartmentalized metabolic network, under infusion of [1,[6][7][8][9][10][11][12][13] C 2 ] glucose.

Animal Preparation
All animal procedures were performed according to federal guidelines and were approved by the local ethics committee from the canton Vaud, Switzerland. The presented 13 C-glucose infusion studies were performed following a previously described protocol (Duarte et al., 2011). Briefly, four male Sprague-Dawley rats (276 AE 11 g, delivered from Charles River Laboratories, France) were fasted for 6 hr prior to the NMR experiment. Animals were intubated and ventilated with 2% isoflurane during surgery. Both femoral veins were catheterized for the administration of glucose and a-chloralose for anesthesia. Two arteries were cannulated for continuous monitoring of physiology (blood pressure and heart rate) and periodic blood sampling for blood gas, plasma lactate, and glucose concentration measurements as well as for further high resolution NMR analyses of 13 C enrichment of substrates. Body temperature was maintained between 37.0 and 37.5 C with a temperature-regulated circulating water bath. Following surgery, anesthesia was switched to intravenous a-chloralose administration (bolus of 80 mg/kg and continuous infusion rate of 28 mg/kg/hr). Animals were placed in a home-built holder and the head position was fixed using ear rods and a bite bar.
After adjustment of MRS parameters, an exponentially decaying 5-min bolus of 99%-enriched [1,[6][7][8][9][10][11][12][13] C 2 ] glucose (1.1 M in saline solution) was administered. The volume of this bolus was adapted to the basal glycemia in order to reach 70% plasma glucose fractional enrichment (FE) at the end of the 5 min. The bolus was followed by a continuous infusion of 70%-enriched glucose at a rate equivalent to the whole body glucose disposal rate of 33.2 mg/kg/min (Jucker et al., 2002) and adjusted based on periodically measured arterial plasma glucose concentrations in order to maintain a constant glycemia of 350 mg/dl throughout the experiment. This protocol results in a stable 70% plasma glucose FE reached within 5 min (Henry et al., 2002). At the end of experiment, rats were sacrificed and brain extracts prepared as previously described (Duarte et al., 2007) for further processing and analysis.

In Vivo and In Vitro NMR Spectroscopy
All in vivo spectra were acquired on a 14.1 T magnetic resonance imaging system interfaced to a 26-cm horizontal bore magnet (Magnex Scientific, Oxford, UK; Varian, Palo Alto, CA), equipped with 12-cm inner diameter actively shielded gradients reaching 400 mT/m in 120 ms. The coil assembly consisted of a home-built 1 H quadrature surface coil and an inner 13 C linearly polarized surface coil.
After initial setting, fast spin echo images (TR ¼ 5 s and TE ¼ 52 ms; eight echoes) were acquired to select a volume of interests of 320 ml in the brain. Magnetic field homogeneity was adjusted using FAST(EST)MAP (Gruetter and Tka´c, 2000). Localized 1 H NMR spectra were acquired using the spin-echo full-intensity acquired localized spectroscopy (SPECIAL) sequence (echo time of 2.8 ms and repetition time of 4 s). During the labeling experiment, Carbon-13 NMR ( 13 C NMR) spectra were acquired in vivo using the semi-adiabatic distortionless enhancement by polarization transfer (DEPT) technique combined with 3D-ISIS 1 H localization (Henry et al., 2003b). Water-soluble metabolites from brain extracts and plasma samples were quantified by 1 H and 13 C NMR spectroscopy with a 14.1 T DRX-600 spectrometer equipped with a 5-mm cryoprobe (BrukerBioSpin SA, Fa¨llanden, Switzerland).

Spectral Analysis
The signal-to-noise ratio (SNR) of acquired 13 C spectra is an important factor for the accuracy and precision in quantification of individual 13 C multiplets of Glu and Gln labeled at positions C4, C3, and C2. To achieve higher SNR, a possibility is to increase temporal averaging of 13 C spectra acquired with 5.3 min time resolution. After correction for phase and frequency drifts, consecutive MRS acquisition blocks were summed two by two resulting in dynamic spectral data with 10.6 min resolution. For every time point, the 10.6-min spectra from four animals were combined to a single data set after correction for frequency drifts. In the averaging process, the possible time shifts between different animals due to intermediate periods of shimming were taken into account. Summed 13 C NMR in vivo spectra were processed using LCModel (Stephen Provencher Inc., Oakville, Ontario, Canada). The basis sets for LCModel were generated using Matlab (TheMathWorks, Natick, MA) by simulating each 13 C isotopomer with the appropriate chemical-shift and J-coupling pattern as described by Henry et al. (2003b). The dynamically measured 13 C spectra were scaled based on the Glu pool size measured with 1 H MRS and the FE of GluC3, which was determined with the following formula: where C4 D34 indicates area of the C4 doublet resulting from double labeling of Glu at position C4 and C3, and C4 S the singlet area of Glu at position C4. Multiplets of GluC4 resonances were averaged over the last 30 min, assuming steady-state enrichment for GluC4 at this later stage of infusion.
Correction factors, found by in vitro 13 C NMR spectra from brain extracts and standard solutions containing the metabolites of interest (Duarte & Gruetter, 2013), were applied to account for the relative different efficiencies in signal enhancement by polarization transfer in DEPT for the considered carbon positions of Glu and Gln.

Metabolic Modeling
The 13 C turnover curves of Glu and Gln were analyzed using a metabolic model built based on the bonded cumomers concept (Shestov et al., 2012) in the case of a two-compartment neuronal-glial metabolic network description (adapted from Gruetter et al., 2001) shown in Figure 1. This model was adapted to include a nonzero concentration of aspartate in the glial compartment and a dilution factor at the level of glial acetyl-CoA as in our previous study (Duarte et al., 2011). This model introduces and assesses a new metabolic modeling approach that allows the analysis of individual 13 C multiplets time courses of in vivo NMR spectra of the rat brain.
The principles of bonded cumomers modeling were explained comprehensively by Shestov et al. (2012). In summary, a cumomer fraction, noted p M{i} , is the sum of isotopomer fractions for all isotopomers of the molecule M labeled at least at the set of positions i f g ¼ i 1 , i 2 , . . . , i n independently from the label at other positions. The size n of the set {i} will be referred to as the order of the p function. Bonded cumomers model consider those cumomers whose indexes refer to adjacent carbons because all possible metabolites isotopomers are not detectable in vivo by NMR. As four different splitting patterns are possible for any carbon with two direct carbon neighbors, they can be expressed in terms of the bonded ' bonded cumomers of order n 3.
The following equation shows the transformation for all observable Glu multiplets of the position C4 in a 13 C spectrum, where the indexes S, D, and Q refer to singlet, doublet, and doublet of doublet in the GluC4 signal, to the related cumomer fractions: Furthermore, if the couplings with the two neighboring carbons in the carbon chain of a given metabolite are identical, such as for C3 of Glu, the above equation simplifies to the following transformation, where the index T represents the triplet signal: In the equation above, intensities are normalized so that the total signal is the sum of all multiplet components. For example, in the case of GluC4, we have: The system of differential equations that describes labeling of all observable multiplets of Glu and Gln and the principal chemical intermediates implied in the generation of these multiplets was derived. The final system describing the incorporation of 13 C from labeled glucose into the multiplets of Glu and Gln consisted in a total of 133 differential equations (see Appendix).
The Cramer-Rao lower bounds (CRLBs) obtained from the LCModel spectral quantification (expressed as coefficient of variation [%standard deviation, SD]) were used as a confidence criterion to select the measurable individual 13 C multiplets of each metabolite. The threshold SD to include a multiplet in the metabolic analysis was set to 15% for every observable 13 C multiplet of Glu or Gln at labeling steady state.
The two-compartment metabolic model was first fitted to the 13 C enrichment curves of total concentration of Glu and Gln C4, C3, and C2 (i.e., the bonded cumomer model reduced to first-order cumomers, i.e., positional model) over time using the Levenberg-Marquardt algorithm for nonlinear regression, coupled to a Runge-Kutta method for nonstiff systems to obtain numerical solutions of the ordinary differential equations. The different 13 C uptake curves used in the regression were weighted according to square root of the inverse of the variance (square CRLBs) extracted from the LCModel quantification for every data point of Glu and Gln, to correct for the overweighting of noisy data in the multiple curve regression (Cobelli et al., 2000). The reliability of the estimated fluxes in brain was derived by Monte-Carlo simulation. Artificial data sets (500 series) were generated by adding random Gaussian noise with same variance as the experimental data to each model turnover curve obtained from the best fit of the experimental data. This resulted in 500 data sets with same characteristics as the experimental data but different noise realizations. All numerical procedures were performed in Matlab (MathWorks, Natick, MA).
In a second step, the bonded cumomer model was used for the fitting of the time courses of both the total 13 C enrichment of every carbon position and every NMRobservable 13 C multiplet of Glu and Gln which was detected experimentally with sufficient precision (<15% CRLB). The precision of the obtained metabolic fluxes was assessed with an equivalent Monte-Carlo simulation.
To investigate the particular changes in the correlations and probability distributions of V dil and V nt using both the positional model and bonded cumomer model, a sensitivity analysis was undertaken for the fluxes V dil and V nt . The mentioned fluxes were constrained to a range of chosen values and their effect on the 13 C labeling curves of Glu and Gln in the positional modeling and bonded cumomer modeling approaches were compared. For each constrained value of V dil or V nt , the metabolic model was fitted to the measured total 13 C enrichment time courses in positional modeling and total and multiplets 13 C enrichment time courses in bonded cumomer modeling. In the positional model, the time courses of 13 C multiplets of Glu and Gln were also simulated in every step of the sensitivity analysis to illustrate the separate effect of V dil and V nt on the simulated 13 C multiplets time courses and therefore the potential of adding those curves to the metabolic modeling process. Finally, the effect of removing the time course of total 13 C enrichment of Glu and Gln C4, C3, and C2 in the bonded cumomer model (referred in following as individual multiplet model) on the precision of the obtained metabolic fluxes was investigated. Statistical significance of differences in flux values found with the various modeling approaches was tested using two-sample unpaired t-test, corrected for multiple comparisons using Bonferroni correction when necessary.

Results
To increase the SNR enabling the measurement of isotopomers at low 13 C enrichment, temporal averaging to a time resolution of 10.6 min of the measured 13 C MRS data was used and 13 C spectra from different animals were combined prior to spectral quantification leading to spectra of higher quality and SNR (e.g., SNR of 24 after 2 hr [1,6-13 C 2 ] Glc infusion, based on the GluC4 S peak). The high quality of the spectra was judged from the well-resolved resonances and appearance of multiplet fine structures of Glu and Gln in spatially localized 1 H-decoupled 13 C NMR spectra ( Figure 2). The total metabolite concentrations measured by 1 H MRS and averaged over four animals were 8.4 AE 1.2 mmol/g for Glu, 4.3 AE 0.3 mmol/g for Gln, and 2.5 AE 0.3 mmol/g for Asp.
The overlap of the central peak of the triplet and the singlet resonances of Glu and Gln labeled at position C3 did not enable their individual spectral quantification, especially at low 13 C concentrations. A higher estimation of the triplet was consistently associated with an underestimation of the singlet and vice versa. Therefore, the sum of singlet and triplet resonances in GluC3 and GlnC3 time courses (GluC3 S þ T and GlnC3 S þ T) was considered as input for the cumomer model, taking advantage of the reliable spectral quantification of the summed multiplets GluC3 S þ T and GlnC3 S þ T, as judged from the respective CRLBs being less than 15% at steady state. The in vivo time courses of total 13 C labeling of Glu and Gln at carbon position C4, C3, and C2 and their observable 13 C multiplets quantified with SD lower than 15% were considered for metabolic modeling analysis during infusion of [1,[6][7][8][9][10][11][12][13] C 2 ] glucose ( Figure 3). The fitted time courses for the positional model included the total 13 C enrichment of Glu and Gln at positions C4, C3, and C2. The fitted 13 C turnover curves for the bonded cumomers approach included the following separate components: GluC4 total, GluC4 S, GluC4 D43, GluC3 total, GluC3 S þ T, GluC3 D, GluC2 total, GluC2 S, GluC2 D23, GlnC4 total, GlnC4 S, GlnC4 D43, GlnC3 total, GlnC3 S þ T, GlnC3 D, GlnC2 total, GlnC2 S (total 17 curves). The same notation as Henry et al., (2003a) and Lanz et al. (2013) has been used for the different 13 C labeling patterns of Glu and Gln.
To determine the metabolic rates upon infusion of 13 C-enriched glucose, either the positional or the bonded cumomer model was fitted to the aforementioned 13 C time courses. The resulting metabolic rates (Table 1) included neuronal TCA cycle ðV n pdh Þ, glial TCA cycle Figure 2. In vivo 13 C NMR spectra acquired three hours after starting [1,6-13 C 2 ] glucose infusion at 14.1 T from a 320 ml volume in the rat brain (summed spectra from four animals with averaged time of 10.6 min, no apodization applied). (b) The fine structure of Glu and Gln at position 4, 3, and 2 is depicted in details in the enlarged spectrum.
, the malate-aspartate shuttle activity (V x ), apparent neuro transmission flux (V nt ), and glial anaplerotic pyruvate carboxylation (V pc ). Assuming metabolic steady state for the metabolite pools, the cerebral metabolic rate of Glc oxidation (CMR glc(ox) ), the inflow of labeling from extracerebral lactate (V in ), and the Gln synthesis rate (V syn ) were calculated (Table 1)   effect of individual 13 C multiplets of metabolites as input information for metabolic modeling, brain metabolic fluxes were also estimated by considering all NMR-observable multiplets of Glu and Gln without total 13 C enrichment time courses (i.e., individual multiplet model) at each carbon position (Table 1). Overall, the estimated metabolic fluxes in all three approaches were precisely determined with SD lower than 20% (determined by Monte Carlo analysis) which owes to the good quality of input dynamic 13 C time courses with adapted time resolution and also the completeness of the used neuronal-glial metabolic model. The most significant differences in the value of estimated fluxes using the different mentioned approaches were for V nt and V dil , with differences of 55% between bonded cumomer model and positional model and 75% between individual multiplet model and positional model.
The probability distribution of metabolic flux values, calculated by Monte Carlo simulation, was fitted with a gamma distribution for every model (positional model, individual multiplet model, and bonded cumomer model) and compared with each other (Figure 4). The probability density for all estimated fluxes determined by Monte Carlo simulations approached a Gaussian distribution for the three modeling approaches, as a consequence of determined flux values clearly different from zero. In general, the distribution of most metabolic fluxes in every approach was narrow, reflecting the precision of the estimated fluxes and the fact that the models are numerically identifiable. Including 13 C multiplets time courses of Glu and Gln in the bonded cumomer model resulted in slightly narrower distributions for most metabolic fluxes. The most significant impact was on V dil distribution for which the full width at half maximum (FWHM) of the probability distribution decreased by a  factor of 5 in the bonded cumomer model compared with the positional model. Removing total enrichment in the individual multiplet model resulted in an approximately twofold increase of the FWHM of most fluxes probability distributions, as compared with the complete bonded cumomer model. The effect of the applied modeling approach on the correlation between metabolic fluxes ( Figure 5) was an overall decrease of correlations between most metabolic fluxes when applying bonded cumomer modeling. V out and V dil showed remarkably strong correlations in every approach, more than 80%, independently of the addition of the 13 C multiplets turnover curves. However, introducing information about time courses of 13 C multiplets of Glu and Gln had significant impact on reducing correlation between V nt and V dil in the bonded cumomer model and individual multiplet model, reduced from 80% to 45% and 30%, respectively.
On the other hand, removing information about total 13 C enrichment of Glu and Gln in the individual multiplet model increased slightly the correlation between several fluxes. For example, the correlation of V dil with V g increased from 13% in bonded cumomer model to 57% in individual multiplet model.
In addition, to visualize the improvements obtained with the bonded cumomer model and better understand the particular effect of it on the correlation of V nt and V dil , sensitivity analysis was undertaken for V nt and V dil . In this analysis, the remaining metabolic fluxes were adapted through the nonlinear regression process to get the best fit of the total 13 C enrichment curves with constrained V nt or V dil . Generally, the time courses of GluC4 and GlnC4 isotopomers showed the highest sensitivity to V nt and V dil , which is illustrated in Figure 6. Total 13 C enrichment time courses of GlnC4 and GluC4 showed low sensitivity to the constrained value of V nt either in the positional or bonded cumomer model. In the positional model, a combination of the remaining fluxes could compensate largely for the fixed V nt constraint, except for the steady-state level of GlnC4 total enrichment. On the contrary, the underlying 13 C multiplet time courses simulated for the positional modeling approach showed a higher sensitivity to changes in V nt in particular for GlnC4 doublet and singlet curves. This observation remained valid when including the 13 C multiplet curves in the fitting process (bonded cumomer model), which shows that no combination of the remaining adjusted parameters enabled to compensate for the fixed V nt constraint. 13 C-labeling curves for doublet and singlet of GlnC4 changed in opposite direction for increasing values of V nt . As a consequence, the total 13 C enrichment curve of GlnC4 was not significantly altered. Sensitivity analysis in both models showed that 13 C time courses of GluC4 isotopomers are not as sensitive as GlnC4 to the value of V dil . From this analysis, it also appeared that the sensitivity of total GlnC4 13 C enrichment to the value of V dil was a reflection of GlnC4 singlet sensitivity, while its doublet remained unchanged.

Discussion
This study demonstrates that using state of the art dynamic 13 C MRS data under 13 C-glucose infusion, the recently introduced concept of bonded cumomer metabolic modeling approach provides an extension to the estimation of neuroglial brain metabolic rates, reducing in particular the correlations between most estimated metabolic fluxes in the brain. This kind of approach was applied here for the first time to an extended range of isotopomers and total 13 C enrichment curves of Glu and Gln (total of 17 uptake curves), quantified from group-averaged in vivo 13 C spectra in rats. The advantages of the bonded cumomer modeling approach and its consistency with previously applied metabolic modeling approaches were evaluated on experimental data.

Preprocessing Approach for Reliable and Accurate Quantification of Multiplets
Previous simulation studies suggested that the inclusion of the time courses of individual 13 C multiplets from in vivo 13 C MRS of the brain in the metabolic modeling analysis can result in a significant improvement in the derivation of metabolic fluxes, determined with better precision (Shestov et al., 2012). Direct 13 C NMR spectroscopy detection at 14.1 T used here benefits from higher SNR and bigger chemical shift dispersion in the NMR signal, potentially allowing higher temporal resolution relative to lower fields (Duarte et al., 2011). The increased chemical shift dispersion of 13 C resonances allowed to better separate the GluC3 and GlnC3 resonances (Figure 2), which was not possible at lower field such as 9.4 T (Henry et al., 2003b). However, the detection and quantification of fine structures of 13 C multiplets by in vivo NMR spectroscopy are hampered by inherent effects of homonuclear spinspin J coupling, which is independent of the applied B 0 magnetic field strength. Therefore, although sensitivity and chemical shift dispersion increase linearly with the applied magnetic field, the spectral separation of the 13 C multiplets (in ppm) decreases. However, with the higher sensitivity allowed by surface coils, 1 H decoupling and efficient B 0 shimming using FAST(EST)MAP (Gruetter and Tka´c, 2000) method at 14.1 T provided enough homogeneity (less than 10 Hz linewidth) to clearly distinguish the 13 C multiplets in the spectral fitting.
The total 13 C enrichment of Glc, Glu, Gln, and Asp were well determined with CRLB lower than 15% at steady state when using a temporal resolution of 5.3 min and SNR of 9 (spectra not shown). However, to overcome the low signal amplitudes of some of the 13 C multiplets, temporal averaging was increased to 10.6 min and spectra were combined from four animals before spectral quantification, which altogether resulted in a threefold elevation of SNR and better defined 13 C spectral baseline. This allowed a more robust fitting using LCModel, which in turn resulted in more accurate and reliable quantification of 13 C multiplets with CRLB less than 15% at steady state ( Figure 3).
Significant overlap of the singlet with the center line of the triplet in GluC3 and GlnC3 resulted in high problem of cross correlation and quantification are same for Glu and Gln. Therefore, only the sum of singlet and triplet multiplets was considered for the quantification. The precision on the fit of the sum of singlet and triplet in GluC3 and GlnC3 was systematically better than the precision of the fit when using singlet and triplet components separately.

Positional Model
Nonlinear regression of the positional model to the total measured 13 C enrichment curves for Glu and Gln at position C4, C3 and C2 resulted in seven independently determined brain metabolic fluxes, without a priori constraints ( Table 1). The estimated fluxes and calculated CMR glc(ox) (0.28 AE 0.02 mmol/g/min) were overall comparable with previous 13 C NMR studies done under a-chloralose anesthesia (for CMR glc(ox) : Hyder et al., 1996;0.18 AE 0.3 mmol/g/min to 0.25 AE 0.03 mmol/g/min; Sibson et al., 1998;0.26 AE 0.04 mmol/g/min;Duarte et al., 2011;0.41 AE 0.02 mmol/g/min; Jeffrey et al., 2013; 0.35 AE0.11 mmol/g/min to 0.38 AE 0.12 mmol/g/min; Lanz et al., 2014;0.41 AE 0.05 mmol/g/min). When quantifying the 13 C spectra acquired from brain of four animals using the traditional preprocessing approach as in Duarte et al., 2011; namely spectral quantification for individual animals with time resolution of 5.3 min and averaging over the individually quantified data from four animals, the obtained brain metabolic fluxes using time courses of total 13 C enrichment of Glu and Gln for carbon positions C2, C3, and C4 showed a significantly different V nt (p < .05, corrected for multiple comparison) as the one estimated in this study (data not shown). The SNR of 13 C spectra acquired from every individual animal brain with time resolution of 5.3 min (SNR of 9 at steady state) was approximately three times lower than the one obtained by combining the 13 C spectra from four animals with temporal average of 10.6 min (SNR of 24 at steady state). Not surprisingly, the lower SNR resulted in 13 C multiplets determined with less precision, which likely reduced the precision of the averaged time courses of 13 C multiplets. These results illustrate the importance of the spectral preprocessing approach and spectral quality on the estimation of brain metabolic fluxes.
In this study, the exchange rate between cytosolic amino acids and mitochondrial TCA cycle intermediates, V x , was estimated more precisely (with a relative error lower than 5%) than in previous studies (Choi et al., 2002;Henry et al., 2002). The estimated V x was slightly higher than the pyruvate dehydrogenase rate, V n pdh and is consistent with malate-aspartate shuttle being the major mechanism for maintaining the cytosolic redox state under normoxic conditions, as reported before (Gruetter et al., 2001;de Graf et al,. 2004;Yang et al., 2009;Lanz et al., 2014).
In the present work, V nt was well determined with a relative SD below 10% using positional modeling under [1,[6][7][8][9][10][11][12][13] C 2 ] glucose infusion ( Table 1). The relative error of all fluxes except V dil was lower than 10%, ascribed to the high quality and SNR of the acquired 13 C spectra at 14.1 T and the refined preprocessing approach used for spectral quantification. Glial labeling dilution of the acetyl-CoA pool, V dil , had a slightly higher coefficient of variation (19%) and presented strong positive correlation to the apparent rate of neurotransmission, V nt . This correlation was attributed to the observed difference in 13 C labeling between Glu C4 and Gln C4. V dil creates a difference between the FEs of Glu and Gln while V nt is responsible for its dissipation by mixing the two pools.

Bonded Cumomer Model
The addition of dynamic time courses of 13 C multiplets of Glu and Gln to the input data for metabolic modeling yielded smaller SDs for most fitted metabolic rates compared with positional modeling (Table 1). Probability density functions for brain metabolic fluxes estimated by Monte Carlo simulations (Figure 4) indicated narrower distributions compared with the positional model, consistent with a higher precision of the determined fluxes. One of the largest differences was found in dilution of 13 C label at the level of glial acetyl-CoA, V dil , which was more precisely determined in the case of the bonded cumomer model approach.
Using the bonded cumomer model, that is, the inclusion of the dynamic 13 C labeling of Glu and Gln multiplets in metabolic modeling reduced the mathematical covariance between most fitted fluxes compared with the positional model ( Figure 5). The assessment of V dil using the positional model resulted in a large covariance between V dil and the apparent glutamatergic neurotransmission rate, V nt . The significantly reduced covariance between V dil and V nt in the bonded cumomer model is expected to reduce the influence of V dil on V nt and should allow a more precise and independent measurement of V dil and V nt . Therefore, the distributions of V dil and V nt determined with the bonded cumomer model are significantly different from ones of the positional model (p ¼ .0001 and .01, respectively). V out and V dil showed remarkably strong correlations in every approach, independently of the addition of the 13 C multiplet turnover curves of Glu and Gln, which is probably due of the fact that both fluxes contribute to dilute the 13 C enrichment at the level of acetyl-CoA.
Among the potential reasons for the higher precision and independency of brain metabolic fluxes is the simple increased number of 13 C time courses and experimental information, arising from dynamic 13 C enrichment of the individual multiplets of Glu and Gln, as input in the metabolic model, including partially redundant information. However, a more important part is the sensitivity of given 13 C isotopomers labeling to particular fluxes, which are then averaged out when summing all the isotopomers of the considered carbon position, which is evident in Figure 6 for GlnC4 sensitivity analysis to V nt and V dil .

Sensitivity Analysis of GluC4 and GlnC4 Multiplets to the Value of V nt and V dil
It is noteworthy that the 13 C enrichment time courses of total GlnC4 and total GluC4 13 C enrichments were less sensitive to V nt than the individual 13 C multiplets of GluC4 and GlnC4; when increasing V nt , the 13 C enrichment of the singlet of GluC4 decreased and the 13 C enrichment of the doublet of GluC4 increased, while the total 13 C enrichment time course of GluC4 did not change ( Figure 6). The remaining metabolic fluxes adapt to the constraints fixed on V nt to provide a similar fit quality for total GluC4 in the positional model. This is a manifestation of the higher correlation between V nt and most of the remaining fluxes found with the positional modeling approach ( Figure 5). For GlnC4, on the other hand, the total 13 C enrichment curve is sensitive to V nt resulting in differences in the 13 C turnover curve for different constraints on V nt . Nonetheless, the GlnC4 multiplet curves exhibited also a stronger sensitivity to V nt, in particular for GlnC4 S, this already at early time points. As a consequence, including the 13 C multiplet curves in the cost function of the regression led to a better-defined value of V nt in the bonded cumomer model compared with the positional model as well as a lower correlation between V nt and most other fluxes in the bonded cumomer model as summarized in Figure 5. Lower correlation of V nt with other fluxes in the bonded cumomer model is expected to result in a more independent estimation of V nt with less impact from other estimated fluxes.
The above argument was based on the flux V nt , which shows the strongest changes in terms of correlations with other fluxes between the positional model and the cumomer model. However, note that a general reduction of the correlation coefficients between fluxes is observed with the cumomer model ( Figure 5), leading to more precise flux measurements. Shestov et al. (2007) also showed lower sensitivity of simulated GlnC4 total 13 C enrichment turnover curve to the value of V nt in the positional model. They concluded that this effect was a consequence of the fast 13 C labeling of the smaller glial Gln pool compared with its indirect precursor, the large neuronal Glu pool. Glial Gln 13 C labeling therefore is largely influenced by the neuronal Glu pool. However, our results show that changing V nt has an effect on the 13 C turnover curves of GlnC4 multiplets in the bonded cumomer model even though the total 13 C enrichment of GlnC4 does not change significantly. This can be interpreted in terms of the different precursors of the essentially astrocytic Gln pool: The Gln pool reflects the 13 C enrichment of the small astrocytic Glu pool, which can be labeled either from the glial TCA cycle or from the neuronal TCA cycle via the neuronal Glu pool. However, the 13 C labeling patterns in the carbon positions 4 and 3 of Glu generated by the astrocytic or the neuronal TCA cycle are different due to the dilution of the 13 C labeling position C3 in the glial TCA cycle through pyruvate carboxylase. This can be seen by a lower total 13 C enrichment in the C3 of Gln as compared with Glu (see Figure 3). The position C4 is also diluted through metabolism of alternative substrates in the astrocytic TCA cycle. Therefore, the probability of generating a doubly labeled C4-C3 Glu molecule from the astrocytic TCA cycle is lower. Increasing the relative value of V nt will generate more Gln molecules with carbon backbones originating from the neuronal TCA cycle, with higher double C4-C3 labeling, increasing the doublet component of GlnC4.
For the same reason, the 13 C labeling time courses of multiplets for Glu and Gln in both models showed that the doublet of GlnC4 is not as sensitive as the singlet of GlnC4 to changes in astrocytic V dil , representing the glial oxidative metabolism of alternative, unlabeled substrates. Therefore, changes in the total GlnC4 enrichment reflect the sensitivity of GlnC4 S to the value of V dil while the 13 C enrichment curve of GlnC4 D43 does not change significantly.

Limitation of Modeling Based on 13 C Labeling Time Courses of Individual Multiplets
Removing total 13 C enrichment curves from the bonded cumomer model and considering just the 13 C labeling of individual multiplets had effect on the estimated brain metabolic fluxes in terms of value, precision, and covariance (Figures 4 and 5). The values found for the brain metabolic rates listed in Table 1 were in the range between those found with the positional model and the full-bonded cumomer model. One drawback of excluding the total 13 C enrichment curves of Glu and Gln in the bonded cumomer model is the reduction of the number of 13 C time courses as input to the metabolic model. Simulation and experimental studies verified that a higher number of experimental 13 C time courses of Glu and Gln carbon positions increased the precision of estimated brain metabolic fluxes (Jeffrey et al., 1999;Shestov et al., 2007). In addition, the spectral quantification of the 13 C multiplet time courses for Glu and Gln has a lower precision due to their lower concentrations when compared with the total 13 C enrichment.
Using [1,6-13 C 2 ] glucose as a substrate and acquiring 13 C spectra at a high field of 14.1 T in the present study led to high peak intensities in the 13 C multiplet resonances of Glu and Gln at position C2, C3, and C4, which made it easier to quantify 13 C multiplet signals with sufficient sensitivity. However, the 13 C multiplets are expected to be much more difficult to detect and quantify at lower magnetic fields, in particular when infusing [1-13 C] glucose and limiting the infusion experiment to 60 to 90 min, as typically undertaken in human studies for cost and practical reasons (Rothman et al., 2011). Therefore, the use of [1-13 C] glucose remains a limitation for a reliable and accurate quantification of time courses for individual 13 C multiplets of metabolites.

Conclusion
We conclude that incorporating the 13 C labeling time courses of multiplets of Glu and Gln measured by 13 C MRS at high magnetic field (14.1 T) in the neuronal-glial metabolic model brings improvement in the reliability and independency of estimated brain metabolic fluxes. An extended analysis of 13 C multiplet time courses of metabolites in vivo requires infusion of doubly-labeled Glc as substrate and strong B 0 magnetic field in order to get both enough 13 C enrichment and acceptable SNR in 13 C multiplets spectra. Therefore, this refined modeling approach may be only applicable for preclinical studies on small animals.

Author Contributions
M. M. D. undertook the essential part of the data analysis and interpretation, drafted the article, and undertook the implementation of the metabolic modeling; B. L. contributed to the study conception and design, supervised and contributed to the implementation of the metabolic modeling, helped with the data analysis and interpretation and revised the article for intellectual content; J. M. N. D. designed the animal infusion protocol, undertook the acquisition of the in vivo data and revised the article for intellectual content; N. K. participated to the analysis and quantification of the in vivo spectra and their interpretation and revised the article for intellectual content; R. G. contributed to the study conception and design, to the analysis and interpretation of the results, and revised the article for intellectual content. All mentioned authors gave their final approval of the version to be published and agreement to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Declaration of Conflicting Interests
In following equations, the indexes b and p specify metabolic pools in brain and plasma, respectively. Metabolic pools in the neuronal compartment are indicated with the index n and in the glial compartment with the index g.
The TCA cycle rate in the neuronal compartment, flux through pyruvate dehydrogenase: V n pdh The TCA cycle rate in the glial compartment: Flux through neuronal glutaminase: V nt The anaplerotic flux through pyruvate carboxylase: V pc Glutamate-alpha-ketoglutarate and oxaloacetateaspartate exchange rate: V x Outflow and inflow of labeling at the level of Lactate: V out and V in The dilution flux at the level of acetyl-CoA: V dil A-1. Glucose Transport Across the Blood-Brain Barrier Gln i, j f g : Asp i, j f g : Enrichment of glial amino acids.
Gln i, j f g : Gln i, j, k f g : Enrichment of single brain pyruvate pool by considering its fast equilibrium with lactate.

½ :
d dt Enrichment of glial TCA cycle intermediates.  OAA g ½ d dt