Radiology Research and Practice

Volume 2014, Article ID 871619, 10 pages

http://dx.doi.org/10.1155/2014/871619

## Multisite Kinetic Modeling of ^{13}C Metabolic MR Using [1-^{13}C]Pyruvate

^{1}GE Global Research, 85748 Garching bei München, Germany^{2}Medical Engineering, Tecnológico de Monterrey, 64849 Monterrey, NL, Mexico^{3}Medical Engineering, Technische Universität München, 85748 Garching bei München, Germany^{4}Nuclear Medicine, Technische Universität München, 81675 Munich, Germany^{5}Chemistry, Technische Universität München, 85748 Garching bei München, Germany

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

Academic Editor: David Maintz

Copyright © 2014 Pedro A. Gómez Damián et al. 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.

#### Abstract

Hyperpolarized ^{13}C imaging allows real-time *in vivo* measurements of metabolite levels. Quantification of metabolite conversion between [1-^{13}C]pyruvate and downstream metabolites [1-^{13}C]alanine, [1-^{13}C]lactate, and [^{13}C]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.

#### 1. Introduction

While ^{13}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 ^{13}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 1-^{13}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 1-^{13}C]alanine (ALA) by alanine transaminase (ALT), 1-^{13}C]lactate (LAC) by lactate dehydrogenase (LDH), and ^{13}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 1-^{13}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 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 ^{13}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 ^{13}C metabolic pathways for PYR, which can typically be observed after injection of pyruvate labeled in the 1-^{13}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 (), 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 and the respective metabolite’s effective signal decay rates by solving a system of linear equations given in differential form The effective metabolite signal decay rate is dominated by relaxation, the respective backward metabolic exchange rate , and a flip angle (FA) term, which also depends on the repetition time (TR), accounting for the irreversible consumption of signal after successive excitations: with Hence, 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 is assumed to be negligible, true 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 (). In Zierhut et al. [8] a series of piecewise defined exponential equations were presented: The first part of the equation takes into account PYR signal loss due to and the injection of PYR with a constant rate from the arrival time until . It nevertheless assumes that no conversion of PYR takes place during injection. The second part, for all time measurements later than , 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 allows for (1) to be solved yielding metabolite signals in time [8]: Alongside the parameters characterizing the PYR input function, these equations contain the same parameters ( and ) that were solved for using 2SDM.

2SIM can be considered as a two-step approach. First, , , and are determined by fitting (4) to the measured PYR signal. is simply calculated by summing 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 , the decay rate , and the flip angle correction (3) can be used to estimate apparent relaxation of PYR.

##### 2.3. Multisite Exchange Integral Model

As described above, the metabolite signal decay rate depends on relaxation, backward metabolic exchange rates , and signal loss from flip angle variations. On the other hand, the PYR signal decay does not depend on backward metabolic exchange, but on forward metabolic exchange rates . 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: Note that both the PYR signal decay rate and the sum of all of the metabolic exchange rates are multiplied by the same term and can therefore be grouped into a total PYR signal decay rate: By replacing (7) in (6), the integral form of the new PYR input function reads The representation of the total PYR relaxation rate 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 values for PYR can be achieved directly using

Utilizing the same term for the metabolite signals, (5) becomes
As seen in (2), the backward exchange rates are inseparably confounded with in the respective signal decay rate of each metabolite. A nonnegligible backward reaction thus leads to an overestimation of the true 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 ^{13}CO_{2} and this could lead to an apparent shortening in . This signifies that the 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 ^{13}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 with acquired at time steps was fitted to MSIM in a constrained least-squares sense; that is,
with cost function
parameters , and lower and upper bounds lb and ub, respectively. While was fixed to the time when the PYR signal reached 10% of its maximum peak value, was set as a fitting parameter accounting for various injection times. On the contrary, the implementation in [8] kept fixed while fitting for . Even though the duration of the injection was known, fixing in function of its peak value and calculating 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 . The utilized bound constraints were set to physically relevant limits: upper bounds of 0.1 s^{−1} for metabolic conversion rates , since they have been reported to be of a smaller order [8, 23], and of 0.005 s^{−1} for the decay rates (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, 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 leading to an overestimation of 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)).