Multisite Kinetic Modeling of 13C Metabolic MR Using [1-13C]Pyruvate

Abstract

Hyperpolarized 13C imaging allows real-time in vivo measurements of metabolite levels. Quantification of metabolite conversion between [1-13C]pyruvate and downstream metabolites [1-13C]alanine, [1-13C]lactate, and [13C]bicarbonate can be achieved through kinetic modeling. Since pyruvate interacts dynamically and simultaneously with its downstream metabolites, the purpose of this work is the determination of parameter values through a multisite, dynamic model involving possible biochemical pathways present in MR spectroscopy. Kinetic modeling parameters were determined by fitting the multisite model to time-domain dynamic metabolite data. The results for different pyruvate doses were compared with those of different two-site models to evaluate the hypothesis that for identical data the uncertainty of a model and the signal-to-noise ratio determine the sensitivity in detecting small physiological differences in the target metabolism. In comparison to the two-site exchange models, the multisite model yielded metabolic conversion rates with smaller bias and smaller standard deviation, as demonstrated in simulations with different signal-to-noise ratio. Pyruvate dose effects observed previously were confirmed and quantified through metabolic conversion rate values. Parameter interdependency allowed an accurate quantification and can therefore be useful for monitoring metabolic activity in different tissues.

Pedro A. Gomez Damian 1, 2, 3 and Jonathan I. Sperl 1 and Martin A. Janich 1, 4, 5 and Oleksandr Khegai 1, 5 and Florian Wiesinger 1 and Steffen J. Glaser 5 and Axel Haase 3 and Markus Schwaiger 4 and Rolf F. Schulte 1 and Marion I. Menzel 1

Academic Editor:David Maintz

1, GE Global Research, 85748 Garching bei München, Germany
2, Medical Engineering, Tecnologico de Monterrey, 64849 Monterrey, NL, Mexico
3, Medical Engineering, Technische Universitat München, 85748 Garching bei München, Germany
4, Nuclear Medicine, Technische Universitat München, 81675 Munich, Germany
5, Chemistry, Technische Universitat München, 85748 Garching bei München, Germany

Received 30 August 2014; Revised 6 November 2014; Accepted 13 November 2014; 8 December 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

While ¹³ C magnetic resonance spectroscopy (MRS) has been utilized for in vivo imaging and spectroscopy of metabolism [1] for a long time, only the development of dynamic nuclear polarization (DNP) helped to overcome the inherent sensitivity limit; as through hyperpolarization using DNP followed by rapid dissolution, the ¹³ C MR signal can be amplified by more than 10,000-fold [2].

One of the most common and viable agents for in vivo use is [figure omitted; refer to PDF] 1-¹³ C]pyruvate (PYR) [3]. After intravenous injection, it is transported to the observed tissue or organ under observation, where it is enzymatically metabolized to its downstream metabolites [figure omitted; refer to PDF] 1-¹³ C]alanine (ALA) by alanine transaminase (ALT), [figure omitted; refer to PDF] 1-¹³ C]lactate (LAC) by lactate dehydrogenase (LDH), and [figure omitted; refer to PDF] ¹³ C]bicarbonate (BC) by pyruvate dehydrogenase (PDH) to varying extent, depending on tissue type and predominant metabolic activity. At the same time PYR is in chemical exchange with [figure omitted; refer to PDF] 1-¹³ C]pyruvate-hydrate (PYRH). As part of gluconeogenesis, PYR may also be carboxylated to oxaloacetate [4].

In order to quantify the metabolic exchange between PYR and its downstream metabolites, MRS data acquired over a certain time period after injection first require assignment of spectral peaks [5] in the spectral domain and second require quantification of these peaks over time. Several different methods have been used for this time-domain analysis, and among these the most simple and robust method is the determination of metabolite signal ratios. These ratios are usually obtained from the peak metabolite signals [6] or through integrating over time [5]. The latter approach has been employed in our previous study, conducted by Janich et al. [5], where hyperpolarized PYR spectra were quantified for different PYR doses and subsequently used to determine the dose effects on Wistar rats based on time integrated metabolite signal ratios.

Although the approach of obtaining relative metabolite signal ratios, LAC to PYR or ALA to PYR, is straightforward and robust, independently if obtained from peak signal or time integrals, the results suffer from an increasingly strong [figure omitted; refer to PDF] weighting of the integral, which skews the resulting ratios. Furthermore, although time-domain visualization and signal ratio determination is an effective tool for assessing the effect of different PYR doses, it provides no quantitative kinetic data of metabolic exchange.

In order to achieve this quantification, different methods for kinetic modeling of hyperpolarized ¹³ C MR data have been reported. Most approaches, derived from the modified Bloch equations, represent a two-site interaction between PYR and one specific downstream metabolite, for example, either LAC or ALA [7-14]. Modeling can be extended to include more sites (intra- and extracellular) or more metabolites [9, 12] (for a comprehensive comparison, see [15]). Even so, presumably for robustness, previous work focuses primarily on fitting data with just one downstream metabolite, keeping most parameters fixed, or even model free, based on signal ratios [5, 16, 17]. When PYR is injected and the corresponding metabolic reactions begin to take place, PYR is not metabolized exclusively into ALA (or LAC), but it changes dynamically into all of the aforementioned downstream metabolites [18]. There is furthermore some skepticism, if the implicit assumption of rate constant stability holds in all applications [17] and there are few analyses on model parameter dependence on SNR [19]. In particular, metabolic conversion in the heart predominantly follows the PDH path producing BC [6, 20]. We therefore hypothesize that the simultaneous consideration of various metabolic pathways is necessary to obtain an accurate evaluation of in vivo metabolic conversion rates. On this basis, we propose using a mathematical framework for multisite modeling (similar to [8, 21, 22]) by simultaneously fitting different possible ¹³ C metabolic pathways for PYR, which can typically be observed after injection of pyruvate labeled in the [figure omitted; refer to PDF] 1-¹³ C] position.

Additionally, although our prior work [5] evaluates quantification of spectra and employed a semiquantitative method to investigate metabolic conversion under different PYR doses (based on metabolite to PYR ratios), it does not provide fully quantitative kinetic data. Therefore, in this subsequent work we employ the experimental data obtained in [5] and implement the proposed multisite, dynamic model to determine metabolic conversion and signal decay rates for full quantification of the kinetics of metabolic conversion. Furthermore, the proposed model gives access to effective longitudinal relaxation times ( [figure omitted; refer to PDF] ), both for PYR and for the downstream metabolites.

Using the identical biological data, the kinetic parameters estimated by the multisite model are then compared to the parameters obtained using the two-site models proposed both in [8] and in [23]. The estimated parameters of all models are also compared between the three different doses utilized in [5], that is, 20, 40, and 80 mM (corresponding to 0.1, 0.2, and 0.4 mmol/kg bodyweight) of PYR, in order to evaluate the capability of the model for the assessment of dose response. As identical data is used, the evaluation allows for direct assessment of kinetic model stability and quality. Ideally, a successful kinetic model would allow the reduction of data variability due to modeling to a minimum, allowing the visualization of biological variability (i.e., as a response to dose treatment, etc.). In addition, using simulated metabolic data based on exemplary conversion rates, we assessed the variability and stability of the kinetic models under the influence of noise. Here, the expectation towards a model is that both systematic bias and standard deviation of the resulting metabolic conversion rates should be as low as possible over a large range of signal-to-noise ratio (SNR).

2. Theory

In our previous study [5], MRS spectral data after injection of pyruvate was acquired and analyzed utilizing time-domain fitting with AMARES [24], resulting in a time course of metabolite levels. To quantify the metabolic conversion, this previous study employed integrated metabolite signal ratios. In the following paragraphs, we will compare this simple integrative approach to kinetic modeling using three different approaches, which are two-site exchange differential model, two-site exchange integral model, and multisite exchange integral model.

2.1. Two-Site Exchange Differential Model

Using a two-site exchange differential model (2SDM) allows computing metabolic exchange rates [figure omitted; refer to PDF] and the respective metabolite's effective signal decay rates [figure omitted; refer to PDF] by solving a system of linear equations given in differential form [figure omitted; refer to PDF] The effective metabolite signal decay rate [figure omitted; refer to PDF] is dominated by [figure omitted; refer to PDF] relaxation, the respective backward metabolic exchange rate [figure omitted; refer to PDF] , and a flip angle (FA) term, which also depends on the repetition time (TR), accounting for the irreversible consumption of signal after successive excitations: [figure omitted; refer to PDF] with [figure omitted; refer to PDF] Hence, [figure omitted; refer to PDF] results in a single, inseparable term of signal decay. However, FA and TR are known from experimentation and can be corrected for. In case the backward exchange rate [figure omitted; refer to PDF] is assumed to be negligible, true [figure omitted; refer to PDF] relaxation times can be quantified; however, it remains unclear whether this assumption holds true in all physiological states of the animal.

2SDM does not assume a PYR input function and for that reason the first order differential equation (1) can be solved as a linear system. This approach is independent of the time course of PYR administration and is therefore straightforward to apply.

2.2. Two-Site Exchange Integral Model

Another approach in kinetic modeling, the two-site exchange integral model (2SIM), assumes a PYR input function that represents the PYR signal in time ( [figure omitted; refer to PDF] ). In Zierhut et al. [8] a series of piecewise defined exponential equations were presented: [figure omitted; refer to PDF] The first part of the equation takes into account PYR signal loss due to [figure omitted; refer to PDF] and the injection of PYR with a constant rate [figure omitted; refer to PDF] from the arrival time [figure omitted; refer to PDF] until [figure omitted; refer to PDF] . It nevertheless assumes that no conversion of PYR takes place during injection. The second part, for all time measurements later than [figure omitted; refer to PDF] , is characterized only by the PYR signal loss. In a similar manner, an assumption on the initial PYR concentration can be made instead of an assumption on the input function, leading to the modeling of only the exponential decay, as shown in [25]. Explicit modeling of [figure omitted; refer to PDF] allows for (1) to be solved yielding metabolite signals in time [8]: [figure omitted; refer to PDF] Alongside the parameters characterizing the PYR input function, these equations contain the same parameters ( [figure omitted; refer to PDF] and [figure omitted; refer to PDF] ) that were solved for using 2SDM.

2SIM can be considered as a two-step approach. First, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are determined by fitting (4) to the measured PYR signal. [figure omitted; refer to PDF] is simply calculated by summing [figure omitted; refer to PDF] and the known injection duration. These parameters are then utilized to fit (5) to the LAC and ALA signals. In [6], this model is also utilized to fit the BC signal. Finally the computed metabolic exchange rates [figure omitted; refer to PDF] , the decay rate [figure omitted; refer to PDF] , and the flip angle correction (3) can be used to estimate apparent [figure omitted; refer to PDF] relaxation of PYR.

2.3. Multisite Exchange Integral Model

As described above, the metabolite signal decay rate [figure omitted; refer to PDF] depends on [figure omitted; refer to PDF] relaxation, backward metabolic exchange rates [figure omitted; refer to PDF] , and signal loss from flip angle variations. On the other hand, the PYR signal decay [figure omitted; refer to PDF] does not depend on backward metabolic exchange, but on forward metabolic exchange rates [figure omitted; refer to PDF] . This signifies that the rate of PYR decay is also proportional to the rate of PYR downstream conversion.

Hence, when passing from 2SIM to a multisite exchange integral model (MSIM), the PYR input function (4)-represented in its differential form-needs to include all of the metabolic exchange rates: [figure omitted; refer to PDF] Note that both the PYR signal decay rate [figure omitted; refer to PDF] and the sum of all of the metabolic exchange rates [figure omitted; refer to PDF] are multiplied by the same term [figure omitted; refer to PDF] and can therefore be grouped into a total PYR signal decay rate: [figure omitted; refer to PDF] By replacing (7) in (6), the integral form of the new PYR input function reads [figure omitted; refer to PDF] The representation of the total PYR relaxation rate [figure omitted; refer to PDF] as the sum of the PYR relaxation rate and the metabolic conversion rates allows for a simultaneous fitting process, where the conversion rates are taken into account also in the PYR input function, creating dependent curves and a parameter interdependency. In addition, the estimation of [figure omitted; refer to PDF] values for PYR can be achieved directly using [figure omitted; refer to PDF]

Utilizing the same [figure omitted; refer to PDF] term for the metabolite signals, (5) becomes [figure omitted; refer to PDF] As seen in (2), the backward exchange rates are inseparably confounded with [figure omitted; refer to PDF] in the respective signal decay rate [figure omitted; refer to PDF] of each metabolite. A nonnegligible backward reaction thus leads to an overestimation of the true [figure omitted; refer to PDF] values for all of the downstream metabolites. For LAC, the overestimation might be considered negligible since the backward reaction was reported to have only a very small effect on kinetics [26], although earlier work indicates upregulated gluconeogenesis in liver-metabolism of tumor-bearing rats [27]. The assumption of negligible backward reactions might also not hold for ALA. There is no need to apply a backward exchange to BC; however, depending on pH, it is breathed out as ¹³ CO₂ and this could lead to an apparent shortening in [figure omitted; refer to PDF] . This signifies that the [figure omitted; refer to PDF] values for ALA and BC obtained utilizing this model can only be considered bounds for the true value.

3. Methods

3.1. Experimental Data

The experimental data was obtained from healthy male Wistar rats through the acquisition of slice-selective FID signals in heart, liver, and kidney tissue. Three different hyperpolarized PYR concentrations (20, 40, and 80 mM, which correspond to an injected dose of 0.1, 0.2, and 0.4 mmol/kg bodyweight) were utilized to measure a total of 15 animals. Each dose was injected into five different animals twice, resulting in a total of 10 measurements for each dose. A flip angle of 5° was utilized and TR was triggered to animal breathing yielding an average value of ~1 s. SNR was calculated by dividing the maximum PYR signal by the average noise for all time steps. More experimental details can be directly found in [5].

Further exemplary data to evaluate modeling performance at presence of pathology were obtained from adult female Fischer 344 rats (Charles River, Sulzfeld, Germany) beating subcutaneous mammary adenocarcinomas. The animals' anesthesia was maintained with 1-3% isoflurane in oxygen starting about 1 h before the first injection. During the experiment, the heart rate, temperature, and breathing signal were monitored using an animal monitoring system (SA Instruments, Stony Brook, NY, USA). All ¹³ C animal experiments were approved by the regional governmental commission for animal protection (Regierung von Oberbayern, Munich, Germany). Two injections were performed using an 80 mM concentration, allowing for direct comparison. For this set of experiments, a flip angle of 10° was utilized and TR was fixed to 1 s.

3.2. Data Processing

The experimental data [figure omitted; refer to PDF] with [figure omitted; refer to PDF] acquired at time steps [figure omitted; refer to PDF] was fitted to MSIM in a constrained least-squares sense; that is, [figure omitted; refer to PDF] with cost function [figure omitted; refer to PDF] parameters [figure omitted; refer to PDF] , and lower and upper bounds lb and ub, respectively. While [figure omitted; refer to PDF] was fixed to the time when the PYR signal reached 10% of its maximum peak value, [figure omitted; refer to PDF] was set as a fitting parameter accounting for various injection times. On the contrary, the implementation in [8] kept [figure omitted; refer to PDF] fixed while fitting for [figure omitted; refer to PDF] . Even though the duration of the injection was known, fixing [figure omitted; refer to PDF] in function of its peak value and calculating [figure omitted; refer to PDF] as a parameter allowed for different delivery and perfusion times. Delivery, perfusion, and export are however not implicitly included in the model. To improve the convergence properties of the optimization, the gradient of the cost function was calculated analytically. The optimization was carried out using the MATLAB function fmincon (MathWorks, Natick, MA, USA) employing the Trust Region Reflective Algorithm and a function tolerance of [figure omitted; refer to PDF] . The utilized bound constraints were set to physically relevant limits: upper bounds of 0.1 s^-1 for metabolic conversion rates [figure omitted; refer to PDF] , since they have been reported to be of a smaller order [8, 23], and of 0.005 s^-1 for the decay rates [figure omitted; refer to PDF] (equivalent to a 200 s inverse effective signal decay rate) and lower bounds establishing the positivity of all parameters. Note that the optimization always converged far away from the bounds and they were only implemented for numerical improvement. After optimization, [figure omitted; refer to PDF] values were estimated for all metabolites from the effective signal decay rate (see (2) and (9)). Initial conditions were fixed to expected normal parameters; however, randomizing the starting guess in between bounds and performing various iterations yielded comparable results.

Pyruvate-hydrate (PYRH), which is also present in spectroscopy, was not included in the minimization process. The reason for this is that conversion between PYR and PYRH is not enzymatic and we are interested in quantifying metabolic rates that lead to a better understanding of enzymatic activity. Additionally, since chemical exchange with PYRH is instantaneous and almost in equilibrium, including PYR would require adding three extra parameters to the minimization without providing additional information regarding metabolic activity. In fact, if PYRH were to be included, the immediate conversion of PYR to PYRH would lead to an overestimation of the apparent metabolic rate, which in turn would decrease all other parameters intrinsic in [figure omitted; refer to PDF] leading to an overestimation of [figure omitted; refer to PDF] values for PYR.

The same reasoning holds for the exclusion of additional pools. Although the MSIM model can be further extended to include multiple pools [15, 22], including them only adds variables to the minimization with no direct benefit to the determination of enzymatic conversion rates.

4. Results

4.1. Convergence and Quality of Fit

Parameter fitting with MSIM was shown to converge to an optimal point for every set of experimental data. Figures 1(a)-1(c) show the fitted curves of all metabolites for all models. The residuals for every metabolite and every measurement in the time domain were analyzed (Figures 1(d)-1(f)), and the error of the fitted curves and computed parameters was determined based on the parameter covariance matrix [28]. This error was utilized to determine 95% confidence intervals on the fitted data (see Figures 1(a)-1(c)).

Example of metabolic data acquired for a 40 mM (0.2 mmol/kg) dose in kidney predominant tissue, fitted curves (solid lines) using (a) MSIM, (b) 2SIM, and (c) 2SDM and 95% confidence intervals (dotted lines). (d-f) Residuals of fit.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

Note that for both MSIM and 2SIM the residuals have a distinct pattern. The pattern indicates that a linear injection rate does not fully model biological activity. In [9], the input function is modeled as a trapezoidal instead of a linear input, but the authors provide no residual analysis. On the other hand, assuming no input function by establishing a fixed initial PYR concentration [25] or solving the differential linear system may not fully account for the entire kinetic time course of the measured signals. In any case, this should be considered as a limitation for both models.

4.2. Model Comparison

For all of the experimental data, parameters were obtained utilizing the 2SDM, the 2SIM, and the MSIM. While a single implementation of MSIM brought forth parameter values for all downstream metabolites, an independent implementation for LAC, ALA, and BC was necessary in the two-site models. Since all three models were applied on exactly the same experimental data, the comparison between them and to the results obtained for the integrated metabolite signal ratios obtained from Janich et al. [5] directly allows assessing model accuracy separated from biological variability and experiment related inaccuracies like low SNR levels. Results from one exemplary minimization are shown in Table 1; Table 2 displays mean estimated [figure omitted; refer to PDF] values for all experiments and their respective SNR levels; and Figure 2 details the obtained metabolic conversion rates for all three models.

Table 1: Exemplary parameter estimates (± standard error) obtained from three different kinetic modeling methods for a 40 mM (0.2 mmol/kg) dose of kidney predominant tissue.

Model	MSIM	2SIM	2SDM
[figure omitted; refer to PDF] [s-1 ]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF] [s-1 ]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF] [s-1 ]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF] [s]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF] [s]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF] [s]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF] [s]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	N/A*

[figure omitted; refer to PDF] According to (1), 2SDM only fits for [figure omitted; refer to PDF] exchange rates and the corresponding [figure omitted; refer to PDF] values.

Table 2: [figure omitted; refer to PDF] calculated for MSIM and 2SIM and corresponding SNR levels for all concentrations and slices (mean ± standard deviation).

	[figure omitted; refer to PDF] (MSIM)	[figure omitted; refer to PDF] (2SIM)	SNR
20 mMol
Heart	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Liver	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Kidney	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
40 mMol
Heart	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Liver	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Kidney	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
80 mMol
Heart	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Liver	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
Kidney	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]

Metabolic conversion rates of LAC (left), ALA (center), and BC (right) obtained for heart, kidney, and liver predominant slices at 20, 40, and 80 mM concentrations (0.1, 0.2, and 0.4 mmol/kg doses) for 2SDM (top), 2SIM (center), and MSIM (bottom). Every box plot displays minima, 25th percentiles, medians, 75th percentiles, maxima, and outliers.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

(g) [figure omitted; refer to PDF]

(h) [figure omitted; refer to PDF]

(i) [figure omitted; refer to PDF]

Conversion rates and [figure omitted; refer to PDF] values calculated with MSIM tended to be lower than those of 2SIM and these in turn are lower than 2SDM (see Tables 1 and 2). Although performance is very similar for all models, reduced data spread can be observed in PYR to LAC conversion in kidney predominant tissue (Figure 2). Since MSIM fits up to nine parameters simultaneously, estimated error from the parameter covariance matrix was usually higher for MSIM.

Additionally, for an exemplary dataset, a noise analysis of all three models was implemented by adding Gaussian noise to different extent. Parameters were first obtained from an exemplary minimization with MSIM and were then subsequently used for time curve simulation. Every model was then fit 1,000 times with different initial parameters to this simulated time curve to create a model specific ground truth. Finally, based once again on 1,000 iterations, the simulated dataset was corrupted with random Gaussian noise and minimized with each model. Figure 3 displays mean and standard deviation of [figure omitted; refer to PDF] values up to a 10% noise level.

Figure 3: Noise level analysis for exemplary simulated data. Error bars show mean ± standard deviation.

[figure omitted; refer to PDF]

Figure 3 illustrates that although all models yield the same results in noise-free data, with increasing noise both bias and standard deviation of the two-site models 2SIM and 2SDM increase. As a consequence, the resulting metabolic conversion rates obtained from these two-site models increasingly suffer from systematic under- or overestimation. In contrast, the simulation demonstrates that the MSIM model remains bias-free, even with increased noise level, while exhibiting the smallest standard deviation compared to the two-site models.

From experimental results, it is clear that SNR increases with higher concentrations of injected PYR and that 20 mMol injections in liver and kidney predominant tissue had the lowest SNR (with corresponding noise levels of nearly 10%), whereas SNR in heart was generally higher but had a larger standard deviation (Table 2). According to noise simulations, it is precisely in low SNR regions that MSIM is expected to perform with lower deviations. Standard deviations for [figure omitted; refer to PDF] values and reduced data spread in 20 mMol [figure omitted; refer to PDF] quantification, especially in kidney predominant tissue, are indications that this holds.

4.3. Pyruvate Dose Assessment

The effects of PYR dose on Wistar rats were examined through the injection of solutions with concentrations of 20, 40, and 80 mM (doses of 0.1, 0.2, and 0.4 mmol/kg) hyperpolarized PYR. Kinetic data was obtained for all downstream metabolites and visualized with the same box plots used in [5]. With this approach, a direct comparison between the results previously obtained and the results obtained with kinetic modeling could be made, using median values as a distance dimension between the results obtained by the different models, rather than as confirmatory values (see Figure 2). As in [5], all median values suggest saturation effects. A more detailed assessment of the PYR dose effects on metabolism and its biological interpretation can be found in [5].

4.4. Tumor Evaluation

In tumor cells, it is well known that conversion from PYR to LAC is elevated even in the presence of oxygen [29, 30]. Additionally, some tumors show changes in alanine transaminase activity, leading to suppression of conversion of PYR to ALA [31-34]. Both effects were quantified by comparing experimental data obtained from a healthy rat and a rat with mammary carcinoma and using MSIM to obtain conversion rate parameters (see Figure 4). It can be seen that, for the same dose, the [figure omitted; refer to PDF] conversion rate was more than four times larger in tumor cells than healthy cells and the [figure omitted; refer to PDF] rate was more than 18 times larger in healthy cells than tumor cells. Therefore, obtained conversion rates provide a quantitative metric of metabolic differences between healthy and tumor cells.

Comparison of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] conversion rates between a healthy rat (from an 80 mM dose in kidney predominant tissue) and a rat with mammary carcinoma.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

5. Discussion and Conclusion

Three different kinetic modeling methods were implemented and investigated for the quantification of time-dependent metabolite levels. The two-site exchange differential model (2SDM) and two-site exchange integral model (2SIM) assume a two-site interaction between pyruvate (PYR) and one specific metabolite. The proposed multisite exchange integral model (MSIM) takes into account various downstream metabolites in one system and allows fitting in a one-step process. That is, all of the parameters are generated in a single minimization, avoiding the need for separate implementations for every specific metabolite and resulting in a robust, optimal convergence far from the imposed constraints.

The three models were compared by taking median values as a distance dimension and, using exemplary simulated data, performing a noise analysis. In this analysis, metabolic exchange rate values obtained with 2SDM and 2SIM showed a bias with increasing noise levels. On the other hand, MSIM showed almost no bias, maintaining the average computed value close to the ground truth even at high noise levels, with smaller standard deviations than 2SDM and 2SIM.

Using the experimental data of [5], all kinetic models were compared between different PYR concentrations to assess the effect of increased PYR doses on in vivo metabolism. Results obtained from all three kinetic models were very similar; however, MSIM yielded smaller data spread for metabolic conversion in low SNR experiments and more accurate effective [figure omitted; refer to PDF] values for PYR as downstream metabolite rates are taken into account during the optimization, while effective [figure omitted; refer to PDF] -estimation in 2SIM requires postprocessing corrections.

MSIM was then further utilized to evaluate model performance in disease. Obtained conversion rates from MSIM showed significant differences in healthy cells in comparison to tumor cells, where LAC conversion was elevated and ALA conversion, on the other hand, was suppressed.

Extending two-site models into a multisite model yields both biological and numerical insight. Biologically, it has been shown that calculated rates give proof of the saturation effects studied in [5] and can be used to quantify metabolic differences between normal and tumor cells. Numerically, a one-step fitting process with parameter interdependency performs marginally better than other fitting methods, particularly in regions with low SNR. Further work with the MSIM model will focus on pixelwise metabolic mapping of cellular activity and its application to different metabolic systems.

Abbreviations

[figure omitted; refer to PDF] :

Index of downstream metabolites: lactate, alanine, bicarbonate, and pyruvate-hydrate

[figure omitted; refer to PDF] :

Metabolic conversion rate from pyruvate to [figure omitted; refer to PDF]

LDH:

Lactate dehydrogenase

ALT:

Alanine transaminase

PDH:

Pyruvate dehydrogenase

CA:

Carbonic anhydrase

2SDM:

Two-site differential model

2SIM:

Two-site integral model

MSIM:

Multisite integral model

[figure omitted; refer to PDF] :

Time dependent signal for metabolite [figure omitted; refer to PDF]

[figure omitted; refer to PDF] :

Time dependent pyruvate signal

[figure omitted; refer to PDF] :

Flip angle function

[figure omitted; refer to PDF] :

Time of pyruvate arrival

[figure omitted; refer to PDF] :

Time at which pyruvate is no longer injected

[figure omitted; refer to PDF] :

Pyruvate injection rate

[figure omitted; refer to PDF] :

Metabolite signal decay

[figure omitted; refer to PDF] :

Pyruvate signal decay rate without metabolic conversion rates

[figure omitted; refer to PDF] :

Pyruvate signal decay rate including metabolic conversion rates

lb:

Vector of lower bounds

ub:

Vector of upper bounds

[figure omitted; refer to PDF] :

Index of all metabolites (lactate, alanine, pyruvate, and bicarbonate)

[figure omitted; refer to PDF] :

Sampling times

[figure omitted; refer to PDF] :

Vector of optimization parameters

[figure omitted; refer to PDF] :

Measured data point for metabolite [figure omitted; refer to PDF] and time step [figure omitted; refer to PDF]

[figure omitted; refer to PDF] :

Time dependent signal of metabolite [figure omitted; refer to PDF] as a function of parameters [figure omitted; refer to PDF] .

Acknowledgments

This work was funded by BMBF MOBITUM Grant no. 01EZ0826/7 and BMBF Grant no. 13EZ1114. The authors take responsibility for the content of the paper. Pedro A. Gomez Damian acknowledges the support from the German Academic Exchange Service (DAAD) and the Modality of Professional Experience (MEP) from the Tecnologico de Monterrey.

Conflict of Interests

Marion I. Menzel, Jonathan I. Sperl, Martin A. Janich, Florian Wiesinger, and Rolf F. Schulte are employed by GE Global Research. All other authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] K. M. Brindle, "NMR methods for measuring enzyme kinetics in vivo ," Progress in Nuclear Magnetic Resonance Spectroscopy , vol. 20, no. 3, pp. 257-293, 1988.

[2] J. H. Ardenkjær-Larsen, B. Fridlund, A. Gram, G. Hansson, L. Hansson, M. H. Lerche, R. Servin, M. Thaning, K. Golman, "Increase in signal-to-noise ratio of >10,000 times in liquid-state NMR," Proceedings of the National Academy of Sciences of the United States of America , vol. 100, no. 18, pp. 10158-10163, 2003.

[3] K. Golman, R. In 't Zandt, M. Thaning, "Real-time metabolic imaging," Proceedings of the National Academy of Sciences of the United States of America , vol. 103, no. 30, pp. 11270-11275, 2006.

[4] M. E. Bizeau, C. Short, J. S. Thresher, S. R. Commerford, W. T. Willis, M. J. Pagliassotti, "Increased pyruvate flux capacities account for diet-induced increases in gluconeogenesis in vitro," American Journal of Physiology , vol. 281, no. 2, pp. R427-R433, 2001.

[5] M. A. Janich, M. I. Menzel, F. Wiesinger, "Effects of pyruvate dose on in vivo metabolism and quantification of hyperpolarized ¹³ C spectra," NMR in Biomedicine , vol. 25, no. 1, pp. 142-151, 2012.

[6] H. J. Atherton, M. A. Schroeder, M. S. Dodd, L. C. Heather, E. E. Carter, L. E. Cochlin, S. Nagel, N. R. Sibson, G. K. Radda, K. Clarke, D. J. Tyler, "Validation of the in vivo assessment of pyruvate dehydrogenase activity using hyperpolarised ¹³ C MRS," NMR in Biomedicine , vol. 24, no. 2, pp. 201-208, 2011.

[7] S. E. Day, M. I. Kettunen, F. A. Gallagher, "Detecting tumor response to treatment using hyperpolarized 13C magnetic resonance imaging and spectroscopy," Nature Medicine , vol. 13, no. 11, pp. 1382-1387, 2007.

[8] M. L. Zierhut, Y.-F. Yen, A. P. Chen, R. Bok, M. J. Albers, V. Zhang, J. Tropp, I. Park, D. B. Vigneron, J. Kurhanewicz, R. E. Hurd, S. J. Nelson, "Kinetic modeling of hyperpolarized ¹³ C₁ -pyruvate metabolism in normal rats and TRAMP mice," Journal of Magnetic Resonance , vol. 202, no. 1, pp. 85-92, 2010.

[9] D. M. Spielman, D. Mayer, Y.-F. Yen, J. Tropp, R. E. Hurd, A. Pfefferbaum, "In vivo measurement of ethanol metabolism in the rat liver using magnetic resonance spectroscopy of hyperpolarized [1-¹³ C]pyruvate," Magnetic Resonance in Medicine , vol. 62, no. 2, pp. 307-313, 2009.

[10] O. Khegai, R. F. Schulte, M. A. Janich, "Apparent rate constant mapping using hyperpolarized [1-¹³ C]pyruvate," NMR in Biomedicine , vol. 27, no. 10, pp. 1256-1265, 2014.

[11] J. M. Park, S. Josan, T. Jang, M. Merchant, Y.-F. Yen, R. E. Hurd, L. Recht, D. M. Spielman, D. Mayer, "Metabolite kinetics in C6 rat glioma model using magnetic resonance spectroscopic imaging of hyperpolarized [1-¹³ C]pyruvate," Magnetic Resonance in Medicine , vol. 68, no. 6, pp. 1886-1893, 2012.

[12] S. Josan, D. Spielman, Y.-F. Yen, R. Hurd, A. Pfefferbaum, D. Mayer, "Fast volumetric imaging of ethanol metabolism in rat liver with hyperpolarized [1-¹³ C]pyruvate," NMR in Biomedicine , vol. 25, no. 8, pp. 993-999, 2012.

[13] M. I. Kettunen, D.-E. Hu, T. H. Witney, R. McLaughlin, F. A. Gallagher, S. E. Bohndiek, S. E. Day, K. M. Brindle, "Magnetization transfer Measurements of exchange between hyperpolarized [1-¹³ C]pyruvate and [1-¹³ C]lactate in a murine lymphoma," Magnetic Resonance in Medicine , vol. 63, no. 4, pp. 872-880, 2010.

[14] P. E. Z. Larson, A. B. Kerr, C. Leon Swisher, J. M. Pauly, D. B. Vigneron, "A rapid method for direct detection of metabolic conversion and magnetization exchange with application to hyperpolarized substrates," Journal of Magnetic Resonance , vol. 225, pp. 71-80, 2012.

[15] C. Harrison, C. Yang, A. Jindal, "Comparison of kinetic models for analysis of pyruvate-to-lactate exchange by hyperpolarized ¹³ C NMR," NMR in Biomedicine , vol. 25, no. 11, pp. 1286-1294, 2012.

[16] D. K. Hill, M. R. Orton, E. Mariotti, J. K. R. Boult, R. Panek, M. Jafar, H. G. Parkes, Y. Jamin, M. F. Miniotis, N. M. S. Al-Saffar, M. Beloueche-Babari, S. P. Robinson, M. O. Leach, Y.-L. Chung, T. R. Eykyn, "Model free approach to kinetic analysis of real-time hyperpolarized ¹³ C magnetic resonance spectroscopy data," PLoS ONE , vol. 8, no. 9, 2013.

[17] L. Z. Li, S. Kadlececk, H. N. Xu, D. Daye, B. Pullinger, H. Profka, L. Chodosh, R. Rizi, "Ratiometric analysis in hyperpolarized NMR (I): test of the two-site exchange model and the quantification of reaction rate constants," NMR in Biomedicine , vol. 26, no. 10, pp. 1308-1320, 2013.

[18] F. A. Gallagher, M. I. Kettunen, K. M. Brindle, "Biomedical applications of hyperpolarized ¹³ C magnetic resonance imaging," Progress in Nuclear Magnetic Resonance Spectroscopy , vol. 55, no. 4, pp. 285-295, 2009.

[19] M. F. Santarelli, V. Positano, G. Giovannetti, F. Frijia, L. Menichetti, J.-H. Ardenkjaer-Larsen, D. De Marchi, V. Lionetti, G. Aquaro, M. Lombardi, L. Landini, "How the signal-to-noise ratio influences hyperpolarized ¹³ C dynamic MRS data fitting and parameter estimation," NMR in Biomedicine , vol. 25, no. 7, pp. 925-934, 2012.

[20] M. E. Merritt, C. Harrison, C. Storey, F. M. Jeffrey, A. D. Sherry, C. R. Malloy, "Hyperpolarized ¹³ C allows a direct measure of flux through a single enzyme-catalyzed step by NMR," Proceedings of the National Academy of Sciences of the United States of America , vol. 104, no. 50, pp. 19772-19777, 2007.

[21] S. M. Kazan, S. Reynolds, A. Kennerley, E. Wholey, J. E. Bluff, J. Berwick, V. J. Cunningham, M. N. Paley, G. M. Tozer, "Kinetic modeling of hyperpolarized 13C pyruvate metabolism in tumors using a measured arterial input function," Magnetic Resonance in Medicine , vol. 70, no. 4, pp. 943-953, 2013.

[22] C. Yang, C. Harrison, E. S. Jin, D. T. Chuang, A. D. Sherry, C. R. Malloy, M. E. Merritt, R. J. DeBerardinis, "Simultaneous steady-state and dynamic ¹³ C NMR can differentiate alternative routes of pyruvate metabolism in living cancer cells," The Journal of Biological Chemistry , vol. 289, no. 9, pp. 6212-6224, 2014.

[23] F. Wiesinger, I. Miederer, M. I. Menzel, E. Weidl, M. Janich, J.-H. Ardenkjaer-Larsen, M. Schwaiger, R. F. Schult, "Metabolic rate constant mapping of hyperpolarized 13C pyruvate," ISMRM , no. 3282, 2010.

[24] L. Vanhamme, A. van den Boogaart, S. van Huffel, "Improved method for accurate and efficient quantification of mrs data with use of prior knowledge," Journal of Magnetic Resonance , vol. 129, no. 1, pp. 35-43, 1997.

[25] T. Harris, G. Eliyahu, L. Frydman, H. Degani, "Kinetics of hyperpolarized ¹³ C₁ -pyruvate transport and metabolism in living human breast cancer cells," Proceedings of the National Academy of Sciences of the United States of America , vol. 106, no. 43, pp. 18131-18136, 2009.

[26] T. Xu, D. Mayer, M. Gu, Y.-F. Yen, S. Josan, J. Tropp, A. Pfefferbaum, R. Hurd, D. Spielman, "Quantification of in vivo metabolic kinetics of hyperpolarized pyruvate in rat kidneys using dynamic ¹³ C MRSI," NMR in Biomedicine , vol. 24, no. 8, pp. 997-1005, 2011.

[27] J. D. Shearer, G. P. Buzby, J. L. Mullen, E. Miller, M. D. Caldwell, "Alteration in pyruvate metabolism in the liver of tumor-bearing rats," Cancer Research , vol. 44, no. 10, pp. 4443-4446, 1984.

[28] D. M. Bates, D. G. Watts Nonlinear Regression Analysis and Its Applications , John Wiley & Sons, New York, NY, USA, 2008.

[29] O. Warburg, "On the origin of cancer cells," Science , vol. 123, no. 3191, pp. 309-314, 1956.

[30] H. Lu, R. A. Forbes, A. Verma, "Hypoxia-inducible factor 1 activation by aerobic glycolysis implicates the Warburg effect in carcinogenesis," The Journal of Biological Chemistry , vol. 277, no. 26, pp. 23111-23115, 2002.

[31] W. Droge, H.-P. Eck, H. Kriegbaum, S. Mihm, "Release of L-alanine by tumor cells," The Journal of Immunology , vol. 137, no. 4, pp. 1383-1386, 1986.

[32] L. Brennan, C. Hewage, J. P. G. Malthouse, G. J. McBean, "Gliotoxins disrupt alanine metabolism and glutathione production in C6 glioma cells: a ¹³ C NMR spectroscopic study," Neurochemistry International , vol. 45, no. 8, pp. 1155-1165, 2004.

[33] H. R. Harding, F. Rosen, C. A. Nichol, "Depression of Alanine Transaminase Activity in the Liver of Rats Bearing Walker Carcinoma 256," Cancer Research , vol. 24, pp. 1318-1323, 1964.

[34] S. M. Ronen, A. Volk, J. Mispelter, "Comparative NMR study of a differentiated rat hepatoma and its dedifferentiated subclone cultured as spheroids and as implanted tumors," NMR in Biomedicine , vol. 7, no. 6, pp. 278-286, 1994.