Model Free Approach to Kinetic Analysis of Real-Time Hyperpolarized 13C Magnetic Resonance Spectroscopy Data

Real-time detection of the rates of metabolic flux, or exchange rates of endogenous enzymatic reactions, is now feasible in biological systems using Dynamic Nuclear Polarization Magnetic Resonance. Derivation of reaction rate kinetics from this technique typically requires multi-compartmental modeling of dynamic data, and results are therefore model-dependent and prone to misinterpretation. We present a model-free formulism based on the ratio of total areas under the curve (AUC) of the injected and product metabolite, for example pyruvate and lactate. A theoretical framework to support this novel analysis approach is described, and demonstrates that the AUC ratio is proportional to the forward rate constant k. We show that the model-free approach strongly correlates with k for whole cell in vitro experiments across a range of cancer cell lines, and detects response in cells treated with the pan-class I PI3K inhibitor GDC-0941 with comparable or greater sensitivity. The same result is seen in vivo with tumor xenograft-bearing mice, in control tumors and following drug treatment with dichloroacetate. An important finding is that the area under the curve is independent of both the input function and of any other metabolic pathways arising from the injected metabolite. This model-free approach provides a robust and clinically relevant alternative to kinetic model-based rate measurements in the clinical translation of hyperpolarized 13C metabolic imaging in humans, where measurement of the input function can be problematic.


Accounting for signal loss due to RF excitations
In this section we derive a general result that demonstrates that the signal loss incurred from each RF excitation can be accounted for by adding a constant term to the relaxation rate constants that depends on the flip angle and repetition time.
In the absence of any RF excitations, the continuous-time differential equation governing the evolution of the longitudinal magnetization of multiple metabolites is dM dt where the components of M(t) are the longitudinal magnetizations of each metabolite, K is a matrix of effective rate constants for the metabolic reactions and R is a diagonal matrix of T −1 1 values for each metabolite. With a repetition time of TR, RF excitations occur at times t n = n·TR, and with a flip-angle of θ we have M(t + n ) = M(t − n ) cos θ where t + n and t − n indicate times just before and just after t n . For t + n < t < t − n+1 , M(t) evolves according to (S1), that is Define a relaxation rate adjustment factor r θ = −TR −1 ln(cos θ), then the magnetization just before excitation n + 1 is where I is the identity matrix. The transverse magnetization resulting from excitation n is M(t − n ) sin θ, the components of which are proportional to the measured metabolite peak areas. This combined discrete-continuous system is therefore equivalent to a purely continuous-time model with overall rate matrix L = K − R − r θ I which is sampled at times n·TR. Effective longitudinal relaxation rate constants can be defined as for metabolite i, then L = K − diag(r 1 , r 2 , . . . , r N ) with N metabolites.

Laplace transforms for area under the curve ratio calculations
The Laplace transform of x(t) is defined as and so the area under the curve (AUC) of x(t) is obtained by setting s = 0 since This identity can be used to derive AUC expressions for systems involving multiple metabolites as follows.
With first-order exchange kinetics and longitudinal relaxation the system evolves according to the following differential equation where M(t) is a vector of the time curves for multiple metabolites, L is the overall rate matrix and M in (t) is a vector of the metabolite input curves. Taking Laplace transforms gives Since there are no hyperpolarized metabolites present at t = 0 we have M(0) = 0 and the above expression can be rearranged to give where I is the identity matrix. Setting s = 0 will give a vector of metabolite AUCs, In our application, hyperpolarized pyruvate is the only input metabolite so where P in (0) is the AUC of the pyruvate input curve. The AUC expression therefore simplifies to where v is the left-most column of L −1 . When AUC ratios are derived from the elements of M(0) the P in (0) terms cancel, and so we obtain the key general result that when the input term contains only one metabolite, AUC ratios do not depend on the time course of the input.
Example 1: Two-site plus-one Model The two-site model [1] considers pyruvate and lactate to be purely intra-cellular, and we augment these with a third metabolite X (e.g. alanine, bicarbonate, etc.) that is in exchange with pyruvate only. The metabolite vector and rate constant matrix are where r i is the effective longitudinal relaxation rate of metabolite i (see Eq. (S4)), and k ij is the reaction rate constant from i to j. For the above rate matrix we obtain where |L| is the determinant of L, and so The lactate to pyruvate AUC ratio is then simply which is proportional to the forwards rate constant k P L and is independent of both P in (0) and the rates associated with the additional metabolite X. An equivalent result for metabolite X also holds Example 2: Three-site plus-one Model The three-site model in [2] considers intra-and extra-cellular lactate separately while pyruvate remains purely intra-cellular. We augment this with a fourth metabolite X that is in exchange with pyruvate only. The metabolite vector and rate constant matrix are where k EL and k LE are the lactate transporter rates out of and into the cell respectively. Similar calculations as above give The total lactate to pyruvate AUC ratio is then This AUC ratio is also proportional to k P L and is independent of the input and additional metabolite rates. With no lactate transport from the cells k LE = 0, which gives γ = 1 and the AUC ratio simplifies to the two-site solution, as expected.

Modified cost function for kinetic model fitting
Although the standard deviation (STD) of the lactate residual is constant during each acquisition, the pyruvate residual starts with a STD around 20 times higher, which decays over around 100 seconds to the same level as the lactate STD. This has a biasing effect on the fitting, as the least-squares (LSQ) cost function is disproportionately influenced by the early pyruvate data with large residuals. We have developed a Bayesian noise model to counteract this bias that leads to a modified cost function that is approximately quadratic for pyruvate residuals of a similar magnitude to the lactate residuals, but rises more slowly than a quadratic for larger residuals.
The lactate signal model is where L n is the nth lactate datum measured at time t n , L(t n , x) is the lactate model function evaluated with rate constants in the variable x and n is a sample from a Gaussian distribution with (constant) variance σ 2 . Similarly, the pyruvate signal model is where ν n is a sample from a Gaussian distribution with variance σ 2 n , i.e. each datum has a distinct variance. These error models lead to data likelihood functions of the form where the distinct pyruvate variances are explicit. A prior distribution is used to model uncertainty in σ n conditional on σ, which is a truncated Jeffrey's prior, as suggested in [3]. This prior encodes the constraint that σ n > σ , i.e. the pyruvate residuals tend to be larger than the lactate residuals, but the prior is otherwise uninformative. Uncertainty in σ n can be accounted for by forming the joint probability of P n and σ n then marginalizing σ n via the following integral p(P n | x, σ) = where any constant terms have been neglected. This function is minimized with respect to x and σ to yield kinetic parameter estimates that are robust to the relatively large initial pyruvate residuals.