Kinetic analysis of hyperpolarized data with minimum a priori knowledge: Hybrid maximum entropy and nonlinear least squares method (MEM/NLS)

To assess the feasibility of using a hybrid Maximum‐Entropy/Nonlinear Least Squares (MEM/NLS) method for analyzing the kinetics of hyperpolarized dynamic data with minimum a priori knowledge.


INTRODUCTION
Hyperpolarized metabolic imaging has enabled significant enhancements of the magnetic resonance (MR) signals of 13 C in a range of small molecules allowing the distinction and imaging of an injected molecule from its downstream metabolites. It has received much interest for the real-time study of kinetics of a range of metabolic processes in cancer (1), the heart (2), the liver (3), and in the brain (4). Metabolic conversion of hyperpolarized molecules can be studied in vitro in whole cells, ex vivo in perfused organs as well as in vivo. The acquisition of a time-series of hyperpolarized MR spectra leads to multi-exponential temporal dynamics that contain information on the underlying metabolic reactions, membrane transport rates, enzyme activities, availability of enzyme cofactors etc.
Dynamic curves can be fitted using different analysis methods to quantify the rates of the kinetic processes described in the data. A common approach to fit hyperpolarized data uses mathematical models characterized by several compartments that interconvert at characteristic rates. This method has been used to quantify the rate of conversion of hyperpolarized pyruvate to lactate in vitro in EL-4 mouse lymphoma cells (1) and in T47D human breast cells (5), as well as to characterize pyruvate metabolism in pig hearts in vivo (6). Compartmental modeling requires a knowledge of the chemical products of the reaction but also a priori assumption of the nature of the enzymatic reactions involved in the process. If the first point does not represent an issue because the number of metabolites can be identified by the number of peaks in the spectrum, the second point is more challenging particularly if the metabolic pathways are complex. A similar problem is encountered for the kinetic analysis of time-activity curves acquired using Positron Emission Tomography (PET). In this case it is often difficult to select a compartmental model among those possible that best approximates the interaction of the PET tracer with the region of interest (ROI). Spectral-based algorithms which require minimum a priori assumptions have been proposed as an alternative approach to compartmental modeling for the kinetic analysis of PET dynamic data (7,8). In this work we explore the possibility of using a hybrid maximum entropy/nonlinear least squares method (MEM/NLS) for the kinetic characterization of dynamic hyperpolarized 13 C data. This approach was originally proposed for the analysis of protein folding by Steinbach and colleagues and similarly to spectral-based algorithms allows the derivation of a spectrum of kinetic rates that characterize the biological phenomena described in the experimental data with minimum a priori assumptions (9).
Spectral-based algorithms, and similarly the MEM/NLS algorithm, can be applied only to single input systems, described by a set of first order linear differential equations (no input delay or input dispersion terms included). These assumptions, do not represent a restriction for their applicability because they are commonly satisfied by the systems assessed using hyperpolarized 13 C MR (10).
The performance of the MEM/NLS algorithm for the kinetic analysis of hyperpolarized in vitro 13 C timeseries was studied using Monte Carlo simulations over a range of known kinetics. The proposed method was validated by fitting the characteristic relaxation curves of hyperpolarized [1-13 C] pyruvate acquired at 9.4T using three different flip angles. Additionally, we used the MEM/NLS algorithm to assess the kinetics of time-series of hyperpolarized 13 C pyruvate-lactate exchange acquired in vitro in whole blood cells as well as to re-analyze the previously published in vitro reaction of hyperpolarized 15 N choline with choline kinase (11).

THEORY
The conversion of a hyperpolarized molecule to its downstream metabolites in vitro can be modeled as a multi-compartmental system described by the Bloch-McConnell equations that include both the effect of chemical reaction and relaxation of the hyperpolarized MR signal (Eq. [1]) (1).
M(t) is a time dependent vector of length N of the longitudinal magnetizations M z , K is a NÂN matrix whose elements k are the characteristic rates of the compartmental model chosen and R is the relaxation matrix containing the longitudinal relaxation rates r 1 ¼ 1/T 1 or accounting for the influence of flip angle (u) can be modified with effective relaxation rates r 1 ¼ 1/T 1 -ln(cosu)/TR for a given repetition time TR. The dimension of the system N depends on the number of metabolites detected after the injection of the hyperpolarized probe into the system of interest. Each element of the vector M(t) gives the evolution in time of one compartment.
The solution of the general system of differential equations presented in Eq. [1] to an impulse input can be formulated for a given metabolite as a function M(t) given by a discrete set of n exponentials: where A i are the characteristic amplitudes of the rates l i . The rates l i correspond to the eigenvalues of the matrix (K-R) in Eq. [1] and are, therefore, combinations of the rates of exchange between compartments and the relaxation rates. The value of n depends on the type of interaction between the compartments (i.e., uni-or bidirectional). For instance, the solution of the differential equation in Eq. [1] for a two compartment bi-directional in vitro system shown in Figure 1 is characterized by n ¼ 2 (bi-exponential kinetics) for both compartment A and B.
For a two-compartment, uni-directional system the dynamics of compartment A is characterized by n ¼ 1 (monoexponential) and compartment B by n ¼ 2 (bi-exponential).
In compartmental modeling the problem of finding the optimal value for A i , l i , and n is often formulated as an optimization problem and Eq. [2] is fit to the experimental data using a nonlinear least square (NLS) algorithm.
Instead of imposing a (discrete) compartmental model, experimental data can alternatively be fitted with a continuous distribution of exponential rates g(l) as follows (12): The distribution that best describes the data is obtained by applying the inverse Laplace transform to M(t) without any assumption on the structure of the data. The Laplace inversion of experimental data characterized by random noise is an ill-posed problem (i.e., an infinite number of solutions is available) and a regularization strategy is required to find the optimal distributions g(l). The maximum entropy method (MEM) finds the unique optimum distributions g(l) by maximizing the following cost function (13): S is the entropy term and x 2 is the normalized mean square error. a is the Lagrange multiplier varied during the optimization processes to maintain x 2 small while maximizing S. This approach gives a unique solution independent of initial guess of the model parameters. MEM, however, biases the results leading to nonrandom residuals (9).
MEM/NLS combines the advantages of NLS with those of MEM. Using this approach, the distribution of rates that best describe the experimental data is derived without any a priori assumption of the number of the kinetic components described in the time-series; it is unique and not biased. The algorithm is based on an iterative process that exploits a continuous distribution derived from the inversion of the experimental data through MEM to insert one exponential at a time into the NLS fit. Every time the MEM algorithm converges the number of peaks n in the continuous distribution is characterized and the NLS fit is performed with one exponential more than those used at the previous iteration. If the number of peaks in the continuous distribution is increased by more than one from the previous iteration, only the n þ 1 peaks characterized by the largest area are considered for the NLS fit in an iterative process until all the peaks that appear in the continuous distribution are accounted for in the discrete fit. The chosen optimal MEM distribution is found when x 2 has decreased by less than 1% with respect to the previous optimization. The optimal NLS fit is that characterized by several exponentials, n, for which x 2 has decreased by less than 1% from the discrete distribution characterized by n-1 kinetic rates (9). The final solution of the hybrid MEM/NLS is a discrete distribution of rates whose number n is that which best describes the nature of the data.
As for spectral-based algorithms in PET, it is possible, under certain hypotheses, to associate a compartmental model to the MEM/NLS solution and derive the characteristic rates of metabolic conversion (see Appendix for details). An advantage of this method is that no a priori assumption on the structure of the metabolic pathway has to be specified. Thus it represents a valid alternative to compartmental modeling when a certain degree of uncertainty on the metabolic pathway of the system is present.

Monte Carlo Simulations
The dynamics of in vitro uni-and bi-directional systems were simulated using the Bloch-McConnell equations (Eq. [1]) written for a two-compartment, bi-directional reaction system, in which a given chemical compound is converted to a single product ( Fig. 1) (Eqs. [5][6][7][8]).
k AB and k BA are the two apparent rates of the enzymatic conversion. r 1A ¼ 1 T1A À 1 TR lnðcosuÞ and r 1B ¼ 1 T1B À 1 TR lnð cosuÞ representing the effective relaxation rates of the hyperpolarized signal of A and B, respectively, that take into account the loss of signal due to T 1 and the application of the RF excitation pulse characterized by a flip angle u and time repetition TR (14). M 0 represents the intensity of the hyperpolarized signal of compartment A at t ¼ 0. The values used in the simulation for k AB and k BA are reported in Table 1. r 1A and r 1B were randomly generated within a chosen interval (0.022 s À1 < r 1A < 0.04 s À1 and 0.03 s À1 < r 1B < 0.05 s À1 ) for both uni-and bi-directional systems. All rate values were chosen within a physiological range characteristic of the in vitro conversion of pyruvate to lactate catalyzed by lactate dehydrogenase (LDH) and in agreement with previous studies (5,14,15). To assess the influence of random noise in the data, each dataset was simulated for a low and a high signal-to-noise ratio (SNR) (SNR ¼ 20 and SNR ¼ 90) in the range of those found in our experimental data. Usually the lowest SNR detected in vitro is associated with the hyperpolarized time-series for the product of the metabolic conversion, in this case represented by compartment B. The SNR was calculated as the ratio between the maximum value of the time-signal intensity curve of compartment B and the standard deviation of the noise of the same curve at thermal equilibrium (SNR B ) (16). Datasets were simulated with a time resolution Dt ¼ 2 s and repeated 100 times at each set of parameters. All simulated datasets were fitted with the MEM/NLS algorithm and the ability of the method to derive the correct number of known kinetic components in the simulated data was assessed.
The values of M 0 , r 1A , r 1B , k AB , and k BA were derived by associating a compatible compartmental model to the MEM/NLS solution and solving the system of algebraic equations shown in the Appendix (Eqs. [A3-A6] and Eqs. [A10-A14] for bi-and uni-directional systems, respectively). These values were compared with the simulated reference values as well as to the estimates derived using conventional compartmental modeling. The percent bias (%BIAS) for the kinetic rates k AB and k BA was calculated as a performance index: where p j and p TRUE are the estimated and true value of the indices p. Monte Carlo simulations were performed using a custom-made Matlab (MathWorks V R ) code, whereas the fitting was performed with MeMExp software (9) available online http://cmm.cit.nih.gov/memexp.

Experimental
Samples containing [1-13 C] pyruvic acid, 15 mM trityl radical and 1 mM Gadolinium (Dotarem) were hyperpolarized in a HyperSense V R (Oxford Instruments) DNP at 3.35T and 1.4K for $ 1 hr. Dissolution was performed using 4 mL of 40 mM TRIS buffer, 1 mM EDTA and neutralized to pH 7 with NaOH to yield a final solution of 50 mM hyperpolarized pyruvate. Experiments carried out on whole blood cells also included 50 mM unpolarized sodium lactate in the dissolution buffer. Polarizations were typically of the order P ¼ 20%.

Hyperpolarized [1-13 C] Pyruvate in Solution
Experiments were performed on a Bruker 9.4T Avance III spectrometer using a BBO broad-band observe probe at 298K. The hyperpolarized sample was contained in a 5mm tube with susceptibility matched plugs to restrict the sample volume to the active region of the coil. Hyperpolarized spectra were acquired from three identical hyperpolarized [1-13 C] pyruvate samples with a single acquisition using three different small flip angles: Hyperpolarized time-series were acquired with a temporal resolution Dt ¼ 2 s. The effective relaxation time constants of [1-13 C] pyruvate which include the influence of T 1 and flip angle (u) was derived by fitting the dynamic time-series of the hyperpolarized signal decay with both a mono-exponential function and the MEM/NLS algorithm and the results were compared.

Hyperpolarized [1-13 C] Pyruvate in Whole Blood
All procedures involving animals were carried out in accordance with the Home Office Guidance on the Operation of Animals (Scientific Procedures) Act 1986, HMSO (London). Blood samples were collected in heparinized MR tubes from male Wistar rats (250-300 g, n ¼ 4) under terminal anesthesia at the same time as the excision of the heart for a separate experiment (17). A total of 100 mL of the hyperpolarized [1-13 C] pyruvate containing unlabelled sodium lactate were vigorously mixed with 500 mL of whole blood in a 5-mm MR tube before being inserted immediately into the bore of a Bruker 9.4T Avance III spectrometer. We assume mixing of the hyperpolarized solution and the blood cells to be instantaneous. Time-series of hyperpolarized 13 C MR spectra (Dt ¼ 2 s) were acquired on a Bruker 9.4T at 310K and a small flip angle excitation (u ¼ 10 ). 13 C hyperpolarized MR spectra were integrated and the resulting time-series fitted with the hybrid MEM/NLS method (MeMExp software). The kinetic rates obtained using this approach were compared with those derived by fitting the data with a traditional compartmental model.

Hyperpolarized 15 N Choline
Hyperpolarized 15 N data analyzed in this work were previously published (11). 15 N labeled choline chloride mixed with 15 mM of trityl free radical (OX63) and DMSO/H 2 O was hyperpolarized at low temperature (1.4 K) for 2 h. After dissolution, a 20 mM choline hyperpolarized solution was collected in an MR tube containing purified human choline kinase and 10 mM ATP as cofactor for the enzyme reaction. A series of 15 N MR spectra (Dt ¼ 3.67 s) of hyperpolarized choline and its metabolic product phosphocholine was detected on a 9.4T MR spectrometer. Previously, the resulting timeseries were fitted with a two-compartment, uni-directional model (Fig. 1 with k BA ¼ 0). We re-analyzed these experimental data with the hybrid MEM/NLS algorithm and compared the results obtained from this method with those previously published (11).

RESULTS
The results from Monte Carlo simulations showed that the hybrid MEM/NLS algorithm is able to discern between bi-and uni-directional in vitro systems by esti-mating the number of exponentials that describe the data. In Figure 2 a representative result from the fitting of simulated data using the MEM/NLS method is shown for both bi-and uni-directional systems. The hybrid algorithm derives two rates for compartment B from the time-series of both systems simulated (black lines in Figures 2E and F), while two and one exponentials are identified from the dynamic curve of compartment A of a biand uni-directional system, respectively (blue lines in Figures 2E and F). This is in agreement with the solution of the Bloch-McConnell equations written for uni-and bi-directional systems.
The percentage error in estimating the number of components using the MEM/NLS algorithm is presented in Table 2 for both bi-and uni-directional systems. The ability of the hybrid method in deriving the correct number of rates from the time-series of compartment A was similar for the range of noise values simulated as well as for the simulated value of k AB and k BA. The percentage correctness (CE) was > 95%, with the exception of the bidirectional system simulated with the smallest k AB and the highest noise (CE ¼ 64%). The ability of the algorithm to derive the correct number of exponentials from compartment B decreased with increasing values of k AB for both simulated systems and was similar for the two SNR B simulated.
For bi-directional systems the maximum value of %Bias kAB and %Bias kBA (mean 6 SD) using the MEM/ NLS algorithm was reported at SNR B ¼ 20 equal to 1% 6 6% for k AB ¼ 0.06 s À1 and À18% 6 57% for k BA ¼ 0.002 s À1 , respectively. The minimum value of %Bias kAB and %Bias kBA was reported at SNR B ¼ 90 equal to À0.3% 6 2% for k AB ¼ 0.06 s À1 and to 22% 6 20% for k BA ¼ 0.008 s À1 . For unidirectional systems, the maximum %Bias value reported for k AB was 7% 6 18% at SNR B ¼ 20 and k AB ¼ 0.002 s À1 , whereas the minimum value was 4% 6 3% at k AB ¼ 0.006 s À1 and SNR B ¼ 90.
The MEM/NLS algorithm was used to fit the monoexponential relaxation decays of hyperpolarized [1-13 C] pyruvate acquired in solution at 9.4T using three different flip angles u ¼ 1 , u ¼ 5 , u ¼ 10 (Fig. 4A). The hybrid method correctly derived a single rate to describe the pyruvate decay curves corresponding to the relaxation rate of hyperpolarized [1-13 C] pyruvate which includes the effect of relaxation r 1p ¼ 1/T 1p and loss of polarization due to the flip angle (Fig. 4B). In Figure 4C the values of T 1p derived using MEM/NLS are plotted against those obtained fitting the experimental data with a mono-exponential function. The effective relaxation times reported using MEM/NLS were T 1p ¼ 54.5 s, 51.1 s, and 41.6 s for a nominal flip angle equal to 1 , 5 , and 10 , respectively. For very low flip angles the decay of the signal is determined only by T 1 (e.g., cos n (1 ) $ 1 for n < 100 pulses). We can therefore calculate the true flip   angles using the relation T À1 eff ¼ T À1 1 À TR À1 lnðcosuÞ and a true T 1 $ 54.5s to give actual flip angles of 4 and 8.6 which are less the nominal flip angles.
A representative hyperpolarized 13 C MR spectrum of pyruvate ($170 ppm) and its metabolite lactate ($183.5 ppm) in whole blood cells is shown in Figure 5A. The peak $179 ppm corresponds to pyruvate hydrate, a nonmetabolically active compound formed when pyruvate reacts with water. The hyperpolarized time-series of pyruvate (blue) and lactate (black) with overlaid MEM/NLS fits (red) are plotted as a function of time in Figure 5B. Fitting resid-uals are presented in Figure 5C. Figure 5D shows the results of the MEM/NLS fits to the hyperpolarized 13 C data. One discrete rate was derived through the hybrid algorithm from the pyruvate curve (black line), whereas two rates were obtained from the lactate time-series (blue lines). This solution is compatible with a twocompartment, uni-directional compartmental model and, therefore, Eqs. [A10-A14] presented in the Appendix were used to derive the rate of enzymatic conversion of pyruvate (compartment A) to lactate (compartment B) k AB as well as the relaxation rates r 1A and r 1B . The values estimated for k AB (k AB ¼ 0.0058 s À1 6 0.0020 s À1 , n¼ 4) are plotted against those derived using a two-compartment, unidirectional kinetic model (k AB ¼ 0.0059 s À1 6 0.0022 s À1 , n ¼ 4) in Figure 5E. The values reported for k AB derived using a two-compartment bi-directional kinetic model were found to be almost identical to those estimated with a unidirectional (k AB ¼ 0.0058 s À1 6 0.0024 s À1 ). The algorithm failed in estimating the values of k BA always returning the value of the lowest boundary condition set. The values reported using MEM/NLS for r 1A and r 1B were equal to 0.02 s À1 6 0.002 s À1 and 0.05 s À1 6 0.01 s À1 , respectively. Fitting the experimental data with a two-compartment, uni-directional kinetic model the values obtained for r 1A and r 1B were 0.03 s À1 6 0.004 s À1 and 0.04 s À1 6 0.005 s À1 , respectively. The values obtained for the rates of decay of the hyperpolarized signal of compartments A and B using a two-compartment, bi-directional model were r 1A ¼ 0.028 s À1 6 0.005 s À1 , r 1B ¼ 0.034 s À1 6 0.011 s À1 . The SNR B reported for these datasets was 35 6 10 (a.u.) in agreement with the values used in the Monte Carlo simulations.
In vitro hyperpolarized 15 N experimental data are presented in Figure 6A (colored dots) with overlaid fits (colored continuous lines) for a two-compartment, unidirectional, first-order kinetic model as previously published (11). The values reported for the rates of the metabolic conversion of choline to phosphocholine using this type of compartmental model are: r 1A ¼ 0.0041 s À1 , r 1B ¼ 0.024 s À1 , k AB ¼ 0.0015176 s À1 . The residuals of the choline and phosphocholine fits are shown in Figure 6C and display nonrandom structured residuals indicating that the chosen model is not adequate to fit the experi-mental data. In Figure 6E, the MEM/NLS shows the kinetic rates that characterize the fit curves of the experimental data obtained using a two-compartments, uni-directional compartmental model as previously published (l 1A ¼ 0.0057 s À1 , l 1B ¼ 0.0245 s À1 , l 2B ¼ 0.0057 s À1 ) (11). Figures 6B and F show the fits to the experimental data using the MEM/NLS algorithm and the kinetic rates derived, respectively. Two rates were identified for the time-signal curves of both choline (l 1A ¼ 0.0104s À1 , l 2A ¼ 0.0039s À1 , n¼ 1) and phosphocholine (l 1B ¼ 0.0138s À1 , l 2B ¼ 0.0086s À1 , n¼ 1). The SNR B values calculated for this dataset was 25 (a.u.) in agreement with the values used in the Monte Carlo simulations.

DISCUSSION
Results from the analysis of simulated and experimental data suggest that the MEM/NLS method is adequate for the kinetic characterization of hyperpolarized dynamic time-series. The proposed analysis method is able to derive the number of kinetic rates from the time-series of compartment A and B for both compartmental models simulated. The estimation of the values of A i and l i is affected by the presence of the noise in the data, nevertheless we showed that it is possible to derive the value of the characteristic kinetic rates (r 1A , r 1B , k AB and k BA ) of the system studied by solving Eqs. [A3-A6] and Eqs.
[A10-A14] (see Appendix) for bi-and uni-directional systems, respectively. The accuracy in estimating the value of k AB with the MEM/NLS algorithm was comparable with that of conventional compartmental modeling FIG. 5. a: Hyperpolarized 13 C spectrum of pyruvate and its metabolites lactate and pyruvate hydrate in whole blood cells. b: Shows representative fits to the data using the hybrid MEM/NLS method of pyruvate dynamics (blue) and lactate dynamics (black) in whole blood cells. c: Residuals to the fits are shown. d: Representative kinetic rates derived through the MEM/NLS analysis. Solid and dashed lines correspond to rates with positive and negative amplitude, respectively. e: The values of k AB representing the rate of conversion of pyruvate to lactate derived using the MEM/NLS approach are plotted against those obtained fitting the experimental data with a two compartments, uni-directional kinetic model.
for both types of systems, whereas the hybrid method overestimated the values of k BA at low values of the rate constant.
MEM/NLS fits of the relaxation curves of hyperpolarized [1-13 C] pyruvate in solution acquired using three different flip angles (Fig. 4) identified only one kinetic rate, consistent with mono-exponential decays. The values of the relaxation rates of pyruvate derived with the proposed method are in excellent agreement with those derived by fitting with a mono-exponential function and with those previously published (18). From our data we calculated the influence of the flip-angle correction on the decay of the hyperpolarized signal and found the actual flip-angle to be less than the nominal flip-angle. Flip-angle calibration was carried out on a thermally polarized pyruvate sample. We do not pursue this fur-ther here although possible reasons for this include nonlinearity of the RF amplifiers at low flip-angle or the influence of radiation damping in highly polarized samples.
We applied the hybrid MEM/NLS algorithm to in vitro hyperpolarized 13 C time-series of pyruvate-lactate exchange in whole blood cells (Fig. 5). MEM/NLS identified a single rate from the pyruvate curve, whereas two rates were detected from the lactate time-series. The conversion of pyruvate to lactate catalyzed by LDH is an exchange reaction that is close to equilibrium. In the absence of other information, it would be therefore reasonable to fit the experimental data to a twocompartment, bi-directional kinetic model. MEM/NLS however identifies only a single rate for the pyruvate curve suggesting that the experimental data are better FIG. 6. a: 15 N time-series of choline (blue) and phosphocholine (black) with a fit obtained assuming one way first order reaction kinetics (continuous lines). c: Shows the nonrandom distribution of the residuals due to the inability of the chosen kinetic model to accurately describe the experimental data. b: Experimental data are plotted with overlaid fits from the hybrid MEM/NLS method for choline (blue) and phosphocholine (black), respectively. d: The random residuals resulted from this fitting are shown. The spectra of discrete rates suggested by a two-compartment, uni-directional kinetic model are shown (e), whereas those identified by the MEM/NLS fit are plotted (f). Solid and dashed lines correspond to rates with positive and negative amplitude, respectively. described by a uni-directional model rather than a bidirectional system. Hyperpolarized pyruvate is transported across the cell membrane through specific monocarboxylate transporters. Once within the cytosol the intracellular pyruvate is in exchange with intracellular lactate through LDH. Extracellular pyruvate does not exchange directly with intracellular lactate. Modeling the metabolic conversion of pyruvate to lactate using a twocompartment, bi-directional kinetic model assumes that the hyperpolarized signal of pyruvate is from the intracellular rather than extracellular compartment. This is only true if the uptake of pyruvate across the cell membrane is faster than the rate of conversion of pyruvate to lactate. If the pyruvate uptake is slower than the rate of metabolic conversion then the hyperpolarized pyruvate signal will be dominated by the extracellular rather than intracellular compartment and transport may be rate limiting. In this case the system is better approximated by a twocompartment, uni-directional kinetic model where the compartment A represents the extracellular pyruvate, whereas compartment B represents the intracellular lactate. This information is lost if the experimental data are analyzed with a compartment model assuming a priori a bi-directional model. The values of k AB for the conversion of pyruvate to lactate derived using Eqs. [A10-A13] were found to be in very good agreement with those estimated using conventional compartmental modeling and with those previously published (15).
The utility of the hybrid method for the analysis of hyperpolarized dynamic data was also considered by reanalyzing hyperpolarized 15 N choline data previously published. Fitting with MEM/NLS demonstrated less structured residuals and improved model fit performance (Figs. 6C and D). Of interest, with MEM/NLS two rates were suggested to describe the choline (compartment A) curve and two for the phosphocholine (compartment B) (Fig. 6F) in contrast with the use of a uni-directional kinetic model (Fig. 6E). These were l 1A ¼ 0.0104 s À1 and l 2A ¼ 0.0039 s À1 , with the slowest corresponding to the r 1c ¼ 1/T 1c of choline. The two exponentials derived from the choline time-series could be compatible with a bidirectional enzymatic conversion. However, 31 P spectra acquired from the same solution before and after adding hyperpolarized choline demonstrated the presence of only ATP before the reaction and only ADP þ PCho after the reaction, respectively (data not shown) (11). These experimental findings suggest that there was no back enzymatic conversion of phosphocholine to choline as expected for the enzyme choline kinase. We speculate that the appearance of an additional rate in the choline curve suggests that the reaction stopped during the timecourse and this can be explained by second order kinetics of the enzymatic conversion. It has been previously shown that choline kinase activity is inhibited by ATP depletion and ADP accumulation as well as by product inhibition by phosphocholine (19). In this experiment a 20 mM hyperpolarized choline solution was added to purified human choline kinase containing 10 mM ATP and therefore the reaction was limited by the co-factor concentration, suggesting the importance of cofactor depletion during uni-directional reaction kinetics of hyperpolarized substrates.
We showed that the MEM/NLS algorithm is an adequate tool for the kinetic analysis of in vitro hyperpolarized dynamic data; however, our study has some limitations to be considered. We tested the sensitivity to noise of the hybrid method only for a range of SNR B s typical of in vitro experiments (20 < SNR B < 90). Hyperpolarized time-series acquired in vivo exhibit higher noise levels than those detected in vitro. The performance of the MEM/NLS algorithm at lower SNR B s should, therefore, be accurately assessed before extending the method to the analysis of in vivo data. This was beyond the scope of this work, which is a proof of concept analysis of the applicability of the hybrid method to hyperpolarized data. Extending the method to the analysis of in vivo data also requires further consideration on the input function. The measurement of the optimal input function during hyperpolarized in vivo experiments is still an unsolved problem. The hyperpolarized signal of the injected molecule in the region of interest (ROI) cannot be used as an input function because it contains information on the delivery of the molecule but also on its metabolic conversion. Recent studies have proposed methods to obtain an accurate measure of the delivery input function directly by measuring the hyperpolarized signal of the injected molecule in the arterial blood or indirectly by detecting the hyperpolarized signal of urea, a non metabolically active compound, injected into the system at the same time of the metabolite of interest (20,21). With an accurate knowledge of the delivery input function, the hyperpolarized time-series of detected metabolites can be deconvoluted before applying the MEM/NLS algorithm. From a mathematical point of view, MEM/NLS represents only one of the different possibilities to solve the problem of fitting a set of unknown exponential decaying functions (Eq. [2]). The solutions proposed for positron-emitting radiolabeled compounds (22)(23)(24)(25) offer an alternative for hyperpolarized data analysis but they need to be tested in the specific context.

CONCLUSIONS
This work is the proof of concept analysis of the feasibility of using the hybrid MEM/NLS method proposed by Steinbach J.P. and colleagues for the kinetic analysis of hyperpolarized in vitro data with minimum a priori assumption of the metabolic conversion studied. We demonstrated that this method is adequate to derive the number of kinetic rates that best describe the hyperpolarized 13 C MR time-series and how this can be used to inform on details of the enzymatic conversions that could otherwise be ignored if the experimental data were analyzed with traditional compartmental models where the type of reaction is assumed a priori. We also showed that it is possible to derive the value of the rates characteristic of the metabolic conversion by associating a compatible compartmental model to the MEM/NLS solution. However, results show that compartmental modeling is more accurate in estimating the values of the unknown kinetic parameters than MEM/NLS. The analysis approach presented here therefore will be useful to inform on the compartmental model that best approximate the biological system observed using hyperpolarized 13 C MR especially when the metabolic pathway assessed is complex or a new hyperpolarized probe is used.