Mutual-information based optimal experimental design for hyperpolarized \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{13}$$\end{document}13C-pyruvate MRI

A key parameter of interest recovered from hyperpolarized (HP) MRI measurements is the apparent pyruvate-to-lactate exchange rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{{\textrm PL}}$$\end{document}kPL, for measuring tumor metabolism. This manuscript presents an information-theory-based optimal experimental design approach that minimizes the uncertainty in the rate parameter, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{{\textrm PL}}$$\end{document}kPL, recovered from HP-MRI measurements. Mutual information is employed to measure the information content of the HP measurements with respect to the first-order exchange kinetics of the pyruvate conversion to lactate. Flip angles of the pulse sequence acquisition are optimized with respect to the mutual information. A time-varying flip angle scheme leads to a higher parameter optimization that can further improve the quantitative value of mutual information over a constant flip angle scheme. However, the constant flip angle scheme, 35 and 28 degrees for pyruvate and lactate measurements, leads to an accuracy and precision comparable to the variable flip angle schemes obtained from our method. Combining the comparable performance and practical implementation, optimized pyruvate and lactate flip angles of 35 and 28 degrees, respectively, are recommended.


Mutual
A key parameter of interest recovered from hyperpolarized (HP) MRI measurements is the apparent pyruvate-to-lactate exchange rate, k PL , for measuring tumor metabolism.This manuscript presents an information-theory-based optimal experimental design approach that minimizes the uncertainty in the rate parameter, k PL , recovered from HP-MRI measurements.Mutual information is employed to measure the information content of the HP measurements with respect to the first-order exchange kinetics of the pyruvate conversion to lactate.Flip angles of the pulse sequence acquisition are optimized with respect to the mutual information.A time-varying flip angle scheme leads to a higher parameter optimization that can further improve the quantitative value of mutual information over a constant flip angle scheme.However, the constant flip angle scheme, 35 and 28 degrees for pyruvate and lactate measurements, leads to an accuracy and precision comparable to the variable flip angle schemes obtained from our method.Combining the comparable performance and practical implementation, optimized pyruvate and lactate flip angles of 35 and 28 degrees, respectively, are recommended.
The potential of hyperpolarized (HP) 13 C-Pyruvate magnetic resonance imaging (MRI) to characterize cancer biology, predict progression, and monitor early responses to treatment is well known (e.g., [1][2][3][4][5][6][7][8][9] ).Ongoing studies in prostate, brain, breast, liver, cervical, and ovarian cancer 1,3,4 have shown that HP 13 C-Pyruvate MRI is safe and effective.One of the central aspects of HP-MRI that make it appealing is the elevated chemical conversion of pyruvate to lactate, even in the presence of oxygen, via lactate dehydrogenase (LDH), also referred to as the Warburg effect 10,11 .The higher production of lactate has been shown to correlate with disease presence, the aggressiveness of the disease (e.g., cancer and inflation), and response to therapy.The rate of pyruvate-to-lactate exchange ( k PL ) is a crucial parameter of interest in locating aggressive disease and assessing the biological state of the tissue.HP-MRI presents a unique opportunity to observe tumor metabolism in vivo 2,3,5,6 and use this information to make inferences about tumor aggressiveness and response to therapy.However, a recent white paper 3 highlights the need for further development of spatial, temporal, and spectral encoding strategies that minimize uncertainty while maximizing the efficiency of HP-MRI.An example of uncertainty is the variability of the reported HP measurements in the literature [12][13][14] .In the present work, we develop an information-theorybased approach to determine the optimal MRI design parameters, such as flip angles, with a goal of reducing the uncertainty in the recovered rate parameter, k PL .
The physics of the dynamic nuclear polarization (DNP) method that enables MR imaging of HP 13 C-pyruvate is described in 3,5 .The time history of pyruvate and lactate magnetization within the imaging voxels constitutes the data of interest.Together with a pharmacokinetic HP-MRI model, pyruvate and lactate magnetization data are employed to recover the pyruvate-to-lactate apparent exchange rate, k PL ; e.g., 2, 15, 16 .Accuracy and uncertainty in the recovered rate parameter depend on the data's information content and the pharmacokinetic model's fidelity.This work aims to determine the MRI design parameters that increase the information content in the data and reduce the uncertainty in the rate parameter.The mutual information (MI) between the data and critical model

Hyperpolarized (HP) MRI model
Consider a tissue domain within which different constituents evolve depending on the local environment which includes, for example, extravascular (interstitial) and vascular.In this work, two spatial compartments, namely, extravascular and vascular, each containing HP pyruvate and lactate and complement of these two constituents, are considered.The model employed is based on 2,15 and accounts for T 1 relaxation loss, signal loss due to excita- tion at each scan, and pyruvate-to-lactate and lactate-to-pyruvate exchanges.We expect a dynamic, spatially, and spectrally localized data from the HP MR data acquisition such that the spatially invariant model may represent the data collection on a voxel-by-voxel basis.
Consider a tissue domain cell of volume | cell | in units of cm 3 .Assuming cell is small enough that the spatial variation of hyperpolarized agents in can be ignored, we let cP (t) and cL (t) denote the relative density of HP pyruvate and lactate, respectively, at time t (in units of s).Here, the relative density of species A ∈ {P, L} is defined as the ratio of the density of species ρ A and the total density ρ , i.e., cA = ρ A /ρ .Then the total mass of HP pyruvate and lactate are ρ|� cell |c P (t), ρ|� cell |c L (t) , respectively, ρ being the mass density (g/cm 3 ) of the continuum mixture.Discrete times at N scans are denoted by t k , 1 ≤ k ≤ N ; θ k P and θ k L are flip angles in k th scan whereas TR k = t k − t k−1 , k > 1 , are repetition times and TR 1 = 0 .With c = (c P , cL ) T , the HP pyruvate and lactate available for measurement at the (k + 1) th scan, k ≥ 1 , are given by 15 (1) where the matrix A accounts for T 1 relaxation losses and pyruvate-lactate exchanges Here T 1,P , T 1,L denote the T 1 relaxation times (s), k PL , k LP pyruvate-to-lactate and lactate-to-pyruvate exchange rates (s −1 ), k ve vascular-tissue exchange rate (s −1 ), and ν e the extravascular volume fraction.In (1), C k denotes the matrix that takes into account the loss of signal due to excitation at k th scan: Lastly, VIF = VIF(t) is the vascular input function (dimensionless).Empirically determined blood flow through nonbranching vessels have been shown to correspond to a gamma variate input function 30 .Thus, VIF is taken to be where σ P , α P , β P are constants, and γ denotes a gamma probability density function given by The constant t0 is associated with bolus arrival time and is treated as one of the uncertain model parameters.Constants σ P and α P are dimensionless while β P is in units of seconds.
In (1), the parameters can be gathered in two different classes: 1) model parameters that depend on the tissue domain and specific voxel and includes P = (T 1,P , T 1,L , k PL , k LP , k ve , ν e , t0 ) , and 2) HP-MRI design param- eters such as repetition times and flip angles ) .For simplicity, the parameters T 1,P , T 1,L , ν e , k LP are assumed to be known and fixed so that P = (k PL , k ve , t0 ) .The default values of model and design parameters are provided in Tables 1and 2, respectively.

The total signal
The magnetization intensity of constituent A , A ∈ {P, L} , at k th scan is assumed to be proportional to the mass of the constituent in cell , i.e., there is a constant C such that Cρ cP (t k )|� cell | and Cρ cL (t k )|� cell | are the total magnetization of pyruvate and lactate, respectively.Without loss of generality, C is such that Cρ|� cell | = 1 .In this case, cP (t k ), cL (t k ) are the total magnetizations, sin(θ k P )c P (t k ) and sin(θ k L )c L (t k ) are the measured signal, and cos(θ k P )c P (t k ) and cos(θ k L )c L (t k ) are the signals remaining for the next measurements, respectively.The signals measured at the k th scan are, see 15 , (2) www.nature.com/scientificreports/i.e., one only measures the sin(θ k P ) and sin(θ k L ) fractions of the magnetization leaving the cos(θ k P ) and cos(θ k L ) fractions for the next measurement.The total signal is the sum of the individual signals at different scans and is given by where the dependence of G on the design parameters K and model parameters P is clear from ( 1) and (6).Note that the total signal depends on the entire history of the magnetization time evolution.The magnetization at the current acquisition depends on the HP signal available from the previous acquisition and the total signal is the sum over each time point.In general, the total signal is also proportional to the magnetization weighted by B1 sensitivity maps of the receive coils used for the acquisitions.Here, we are considering the HP data fits on a per pixel basis and assume that B1 variations are small across a given pixel.

Synthetic data
To verify the uncertainty reduction in k PL using the optimal design parameters, HP-MRI experiment is syntheti- cally simulated using the model in (1) with the optimal design parameters and the control design parameters (default values listed in Table 2).Suppose K S is the design parameter associated with some scenario S ; S = default correspond to the default design parameters, S = OED SNR correspond to the optimized design parameters for the specific SNR value.
Using the pharmacokinetic model (1), "ground truth" (signals at N scans) is generated using K S design parameters and default model parameters listed in Table 1.The data (simulation results) is denoted by Y S and takes the form: where, s i P , s i L denote the pyruvate and lactate signals at i th scan, see (6).A total of 25 samples of noisy data for different values of SNR are considered.Noisy samples are computed as follows: where σ s (SNR data ) , given in (21), is the amount of noise in measured signals that depends on the assumed SNR, SNR data , and a j , b j ∼ N (0, 1) for each j = 1, 2, .., N .Since SNR in the data, SNR data , is expected to vary with pixelwise location, a range of SNR data for different optimal design parameters is considered to comprehensively evaluate the accuracy and precision of the k PL parameter recovery.The noise values correspond to the previous SNR range considered: σ s (SNR data ) for SNR data ∈ {2, 5, 10, 15, 20}.

Mutual information based optimization of MR scan parameters
A major goal of this study is to formulate an optimization problem to determine the design parameters ) such that the MRI measurements reduce uncertainty in the rate parameter, k PL .Treating total signal, G defined in (7), as the data, and model parameters, P = (k PL , k ve , t0 ) , and data as random variables, an optimization problem of maximizing the mutual information (MI) between the data and model parameters is proposed.It is shown that uncertainty in recovered k PL from synthetic noisy data is reduced when optimal design parameters are considered; see section "Results".
In what follows, after defining some preliminary quantities, the mutual information between data and the model parameters is defined.Let z ∈ Z = R , P ∈ ⊂ R 3 , and K ∈ D ⊂ R 3N−1 , where Z, , D are data, model parameter, and design parameter spaces, respectively.It is remarked that mutual information are intractable and suffer from the curse of dimensionality as the model complexity is increased to consider more variables.The curse of dimensionality appears from the nested integrals inherent to the mutual information calculation.Further, mutual information calculations in this work utilizes the spatially invariant model while incorporating spatial variations in the mutual information through the parameter variance.For example, spatial variations in the T1 relaxation times among biological compartments are expected and are represented by the variance in the T1 parameters of our model.Our Bayesian approach explicitly accounts for modeling uncertainty including differences among biological compartments through the parameter variance.
Prior, likelihood, and evidence Suppose p 0 (P) is the prior probability distribution function (PDF) of model parameters P and p(z) is the PDF of the data z .Then the joint PDF p(z, P) must satisfy ( 6) where p(z|P) is the conditional PDF of data conditioned on model parameters P (also referred to as the likeli- hood function) and p(P|z) the conditional PDF of model parameters P conditioned on data z (posterior of P ).The prior is taken as multi-variate Gaussian with mean µP and covariance matrix P : Here �µP − P� 2 . Within this Bayesian setting, P , is representative of biologi- cal variability.Within the scope of this manuscript, biological variability refers to potential patient specific differences in the pyruvate to lactate exchange rates, vascular tissue exchange rates, and bolus timing of the injected pyruvate.
To derive the expression for the likelihood function, first suppose that G = G(K, P) is the model prediction of data.Data and the model prediction are assumed to be related as follows: where an additive model of noise is assumed and the noise, ε , is taken as Gaussian with a zero mean and standard deviation σ z , i.e., ε ∼ N (0, σ 2 z ) .Here, we assume that the the phase of the real and imaginary MR data acquisition is constant over time such that the signal may be phase corrected.The MR data acquisition is, in general, complex valued with real and imaginary components that are independent and Gaussian.However, within our phase corrected approach, only the real component is non-zero and Gaussian noise for the real channel is appropriate.
Therefore, the likelihood function p(z|P) takes the Gaussian form: Here || • || denotes the Euclidean norm.Technically, p(z|P) is also conditioned on K , but, for simplicity in nota- tion, the dependence on K is suppressed.
With p 0 (P) and p(z|P) defined as above, the joint PDF p(z, P) can be computed using (9).Further, using ( 9), the PDF of data z ∈ R (evidence), p(z) , can be computed by marginalizing p(z, P) with respect to P: where is the parameter space.

Mutual information
Next, the statistical mutual information is defined and the optimization problem for design parameters K is posed.Given HP-MRI data, the accurate inference of pyruvate-to-lactate exchange rate parameters (and other parameters in P ) depends on the specific choice of control (design) parameters K as selection of K affects the information content in the measured data.The notion of mutual information 29 provides a way to quantify the information content about the model parameters P in the data z .The mutual information between the two random variables z and P with PDFs p(z) and p 0 (P) and the joint PDF p(z, P) is defined as Here, the mutual information I depends only on design parameters K and the forward model (1).
Using Bayes theorem, p(z, P) = p(z|P)p 0 (P) , it can easily be shown that or, where the second term in the above equation is identified as information entropy H(z; K) and the first term as negative of the cross-information entropy, H(z|P; K) .That is

Optimization problem
In order to maximize the reduction in the uncertainty in the model parameters (i.e. to have the most confident estimates of the parameters P ), we propose to maximize the mutual information between the observation data and parameters of interest: (9) p(z, P) = p(z|P)p 0 (P) = p(P|z)p(z), (10)   where K * is the design parameter for which I is maximum (assuming K * exist).
Given a Gaussian prior and conditional probability of the data with respect to the parameters, the entropy of the data conditioned on the model parameter, i.e., H(z|P) , is constant.Therefore, the optimization problem simplifies to

Approximation of mutual information and information entropy
Mutual information calculations are computationally demanding due to high-dimensional integration over the parameter and data spaces.Combinations of both quadrature and sampling methods have been employed for mutual information calculations, each of which is well-suited to certain function classes [32][33][34][35][36][37][38][39][40][41][42][43][44] .These methods include Monte Carlo and Quasi-Monte Carlo methods [33][34][35] , lattice rules 36 , adaptive subdivision 37,38 , neural network approximations 39 and numerical quadrature 40 .Here the problem structure is used to accelerate computations and facilitate tractable numerical integration.Following 45 , Gauss-Hermite quadrature is applied in each dimension of mutual information integrals defined in (16) to numerically integrate multi-variate Gaussian random variables.Quadrature order is increased until convergence is observed.Finite difference time stepping is used to solve the ODE-based model presented in Sect."The total signal".Automatic differentiation accelerated optimization for OED calculations.For constant TR optimization, the auto-differentiation functions of MATLAB were used to calculate gradients of (16).In particular, design parameters K and state variables c were considered as optimization variables to minimize the objective function (16)  with respect to the model constraints (1).Auto-differentiation provides the derivatives of the objective function and constraints with respect to this full space formulation.Given the derivatives in the full space formulation, the reduced space gradient of the objective function ( 16) with respect to the design parameters K may be calculated using an adjoint solve.For varying TR optimization, adjoint gradients were calculated by hand.
Inverse problem to recover rate parameter.A MATLAB routine fmincon is used to solve the inverse problem of recovering model parameters P in the model (1) from the data generated in Sect."Synthetic data" As an objective function for the inverse problem, square of the difference between data and model prediction of signals is used.Similarly, derivatives from the automatic differentiation feature in the MATLAB are utilized for numerical optimization.
The Cramér-Rao bound is computed to provide a reference for the uncertainty observed in the recovered k PL .
Here, J is the Fisher information matrix.Each time point of the HP signal evolution, Eq. ( 1), may be considered as an independent random variable, thus following 46 , the Cramér-Rao bound may be computed analytically with the Fisher information matrix given as follows.

Results
In this section, the main results of our analysis are presented.First, the optimal design parameters for different signal-to-noise ratios (SNRs) are shown.Optimal design parameters for both temporally constant and varying flip angles at each data acquisition are considered.Next, the reduction in uncertainty of k PL when using optimal design parameters generated synthetic data is shown.

Optimized design parameters
As mentioned in Sect."Mutual information based optimization of MR scan parameters", multi-variate Gaussian is taken as a prior for uncertain model parameters, P = (k PL , k ve , t0 ) .The mean and diagonal covariance matrix are fixed to For this choice of mean and variance, all quadrature points for a fifth order Gauss-Hermite quadrature approximation of numerical integration were positive.The remaining model parameters are fixed according to Table 1.
Next, to fix the likelihood function, the Gaussian noise distribution, i.e., ε ∼ N (0, σ z ) , is needed to be fixed.To consider reasonable values of σ z , first the reference peak pyruvate signal s Pref is calculated using the solution of the model with the default model and design parameters in Tables 1 and 2; it is found to be s Pref = 0.6173 .Then for different signal-to-noise ratios ( SNR ), the noise (standard deviation) in the individual signals, σ s , and the standard deviation of the total signal, σ z , are computed as follows, for SNR ∈ {2, 5, 10, 15, 20}, (20) µ P = (0.15, 0.05, 4), � P = diag(0.03,0.01, 1.3). ( where N is the total number of scans.For each SNR and corresponding σ z (SNR) in the above list, optimal design parameters are obtained by solving the optimization problem (19).For simplicity, let K OED SNR denote the optimized design parameter corresponding to SNR and σ z (SNR) .The rationale behind considering different SNR is that, in reality, data is expected to have varied signal-to-noise ratios and this is shown to impact the choice of optimal design parameters.The initial values of the design parameter K are listed in Table 2. Two cases of optimization as described below are considered: • Constant flip angle optimization.In this case, TR k is fixed to 3 seconds for all k, and the flip angles are assumed to be same for all scans.As a result, the optimization problem involved only two variables, θ P and θ L .• Variable flip angle and TR optimization.In this case, flip angles and TR values at all scans are optimized.
For the case of constant flip angle optimization, the optimized design parameters, K OED SNR , corresponding to the five SNRs are tabulated in Table 3. Figs. 1 and 2 represent optimal design parameters, K OED SNR , for SNR = 2 , SNR = 10 , and SNR = 20 for the two cases of optimization problems, respectively.Fig. 1 presents the optimal solution when considering a fixed repetition time of 3s and optimizing for pyruvate and lactate flip angles that are constant in time.Fig. 2 presents the optimal solution when allowing the repetition time and flip angles to vary at each acquisition for a fixed number of data acquisitions, N .The optimal values of design parameters are shown in Fig. 1(a-c) and Fig. 2(a-c).The time varying optimized design parameters are significantly different Table 3. Optimized design parameters considering constant flip angles throughout the scan.Repetition time is fixed to TR k = 3 s, for all k. from the constant value flip angle scheme.Fig. 1(d-f) and Fig. 2(d-f) show the transverse magnetizations for s k P , s k L .For the case of variable design parameters, optimized TR values are shown in Fig. 3.

Validation of the uncertainty reduction using optimal design parameters
The basic workflow in verification includes generating the samples of noisy data associated with different design parameters following Sect."Synthetic data" and then solving the inverse problem to recover uncertain model parameters for each sample of noisy data.From the recovered k PL for different samples of data, the mean and the standard deviation is computed.Specifically, the standard deviation is used as a measure of the uncertainty in the recovered k PL .and varying flip angle scheme; respectively.Fig. 4(b) and Fig. 4(e) correspond to K OED 10 for a constant and varying flip angle scheme; respectively.And finally Fig. 4(c) and Fig. 4(f) correspond to K OED 20 for a constant and varying flip angle scheme; respectively.The Cramér-Rao bound is provided as a reference for the lower bound of the variance in each.Fig. 5 provides a control for the accuracy and precision obtained for flip angles similar to those currently in use in our clinic 47,48 , θ k P = 20 and θ k L = 30 .Similarly, the Cramér-Rao bound is provided as a reference for the lower bound of the variance.The results show that the uncertainties in recovered parameter using K OED 2 are comparable to the current clinical pulse sequence implementation.In fact, K OED 2 demonstrates improved performance for all SNR data except SNR data = 20.
Generally, an improvement in the accuracy and precision of the recovered parameter is seen with increasing SNR data .A time-varying flip angle scheme leads to a higher parameter optimization that can further improve the quantitative value of mutual information over a constant flip angle scheme.However, the constant flip angle scheme, 35 and 28 degrees for pyruvate and lactate measurements, leads to the accuracy and precision comparable to the variable flip angle schemes obtained from our method.

Discussion
The evaluation of accuracy and precision as a function of pulse sequence design is effectively a bi-level optimization problem where the goal is to solve two nested optimization problems: (1) find the pulse sequence that produces the best accuracy and precision for the (2) best curve fit to the data.Direct numerical optimization of the bi-level cost functions(s) is challenging 49 .The mutual information objective function of this study as well as signal maximization and Fisher information in the literature 15,18,19 are effectively proposing a surrogate objective function as a numerically tractable approximation to the bi-level optimization problem of interest.Within the mutual information based optimal experimental design formulation, a time varying flip angle and repetition time scheme is seen to provide significant differences in the pulse sequence as compared to the case when excitation angles are fixed to a constant value over time with a fixed repetition time of 3 s.Indeed, the varying flip angle and repetition time scheme leads to a higher dimensional parameter optimization that provides more degrees of freedom to further improve the quantitative value of mutual information over the constant flip angle scheme.However, as seen in Fig. 4, the constant flip angle scheme leads to comparable accuracy and precision when considering the inference from noise-corrupted data.The time varying scheme is seen to be more sensitive to noise corruption of the expected signal and is generally seen to have the higher variance in the parameter recovery at lower SNR data .The mutual information calculations are generally seen to achieve the Cramér-Rao bound with higher SNR.
The purely theoretical nature of these results is a limitation of the study.However, the two-compartment model analyzed in this work is utilized in numerous studies [50][51][52][53][54][55][56] and, when parameterized with physically meaningful values of the parameters, good data agreement is shown within these studies.Results of this manuscript are thus relevant towards guiding future data collection.
The reduction in the recovered variance is seen to be correlated with the assumed noise value added to the data.Intuitively, less noise resulted in less variance in the parameter recovery.Less intuitively, the optimal MI solutions for flip angles are seen to vary with the noise value of the signal conditional probability model p(z|P) .The greatest reduction in measurement uncertainty is seen for the MI optimal solution corresponding to low SNR of the signal conditional probability model.Here, the lower flip angle is applied to the non-injected secondary substrate and higher pyruvate signal is maintained throughout the acquisition.This could be due to the system being so signal limited for the low SNR case that it is forced to leverage the pyruvate signal to extract any additional information about the metabolic exchange rate.The excitation angles optimized for the higher SNR condition reduce the pyruvate excitation angle to save magnetization for subsequent conversion while simultaneously increasing the lactate flip angle.Within the time varying optimization, the pyruvate excitation angle is reduced to zero after 20 s.For high SNR this suggests that the lactate signal is sufficient to accurately determine the metabolic exchange rate and measuring a large pyruvate signal after the initial bolus is less important.
This work considers uncertainty in the vascular-tissue exchange parameter, bolus arrival time, and rate constants modeled through a Gaussian prior.However, a more comprehensive evaluation of additional uncertain parameters would further evaluate the stability of our results.Additionally, the effect of alternative prior formulations such as uniform distributions for the prior parameters may also be investigated.The numerical computation in this work is also limited by the quadrature scheme for numerical integration of the mutual information integrals.Adding additional sources of uncertainty suffers from the well-known curse of dimensionality 57 and alternative integration schemes such as Markov chain Monte Carlo may be more effective.
Further, the current approach considers the real component of the readout and assumes SNR such that Gaussian statistics is an appropriate noise model for the signal acquisition.Rician statistics 58 is known to be more appropriate as the noise model for low SNR and the low SNR range is expected to be more important toward the end of the HP data acquisition as the signal decays.Rician statistics will be considered in future efforts to optimize acquisition parameters at low SNR or when considering both the real and imaginary components of the signal magnitude.
Alternative to the spatial-invariant model, high-fidelity models may also be considered to determine optimal design parameters and to recover model parameters from the data.However, for such an approach to work, a realistic high-fidelity model is needed keeping in mind the major factors in HP-MRI physics.Additional model fidelity may include permutations of lactate and pyruvate that are endogenous as well as hyperpolarized.Intravascular, extracellular, and intracellular species may also be considered.Nonlinear pyruvate-to-lactate conversion parameters such as from a Michaelis-Menten relationship may also be considered.Additional formulations may also consider the impact of blood flow in the simulations directly though Dirichlet boundary conditions, www.nature.com/scientificreports/as convective transport through porous media 59 , or as a sophisticated 3D-1D coupling with vasculature treated as 1D curvilinear segments 60,61 .In summary, our results suggest that the constant flip angle scheme corresponding to K OED 2 is the best choice in terms of accuracy and precision of the parameter recovery.Results at K OED 2 , θ k P = 35 and θ k L = 28 , are comparable to the current clinical pulse sequence implementations, θ k P = 20 and θ k L = 30 , and demonstrate an improved performance at low SNR data .Further, the constant flip angle scheme may represent a practical choice for implementation on the pulse sequence hardware.

Figure 1 .
Figure 1.Optimized design parameters along with the signals of pyruvate and lactate obtained from the solution of the forward model using optimal design parameters.In (a), for noise σ z (2) , i.e., σ z (SNR) for SNR = 2 , the optimized flip angle scheme is shown for constant flip angles throughout the acquisition (optimal angles are θ k P = 35 degrees and θ k L = 28 degrees, for all 1 ≤ k ≤ N ).The corresponding signal evolution of the transverse magnetizations () are shown in (d) for the constant flip angle case.Similarly, (b) and (c) present the optimized flip angles for SNR = 10 and SNR = 20 ; respectively.The corresponding signal evolution of the transverse magnetizations are shown in (e) and (f).

Figure 2 .
Figure 2. Optimized design parameters along with the solution of the forward model.In (a), for noise σ z (2) , i.e., σ z (SNR) for SNR = 2 , the optimizer considers jointly varying the flip angle and repetition time each acquisition.The corresponding signal evolution of the transverse magnetizations () are shown in (d).Similarly, (b) and (c) present the optimized flip angles and repetition time for SNR = 10 and SNR = 20 , respectively.The corresponding signal evolution of the transverse magnetizations are shown in (e) and (f).Note here that x-axis in all plots is for time and therefore the plots implicitly include the values of optimized TR k , 1 ≤ k ≤ N.

Table 1 .
. Default model parameters, P , and remaining fixed model parameters used in initialization, optimization, and verification steps.