0812

Flexible and Efficient Cramér-Rao Bound Optimization for Magnetic Resonance Fingerprinting using Automatic Differentiation of Bloch Simulations

Philip K. Lee^1,2, Lauren E. Watkins^2,3, Timothy I. Anderson¹, Guido Buonincontri⁴, and Brian A. Hargreaves^1,2,3

¹Electrical Engineering, Stanford University, Stanford, CA, United States, ²Radiology, Stanford University, Stanford, CA, United States, ³Bioengineering, Stanford University, Stanford, CA, United States, ⁴IRCCS Fondazione Stella Maris and Fondazione Imago7, Pisa, Italy

Synopsis

The Cramér-Rao Lower Bound (CRLB) is a metric for optimizing quantitative sequences that requires an analytical expression for the signal. The CRLB for Magnetic Resonance Fingerprinting (MRF) has a complex formulation that makes it difficult to account for system imperfections or relevant signal contributions such as diffusion. We apply automatic differentiation to Bloch simulations and choose flip angles and repetition times that optimize the CRLB of the MRF sequence without deriving an explicit analytical expression for the signal. This method is computationally efficient and can easily be extended to include B0 and B1 inhomogeneities. Results are validated with in-vivo measurements.

Introduction

The Cramér-Rao Lower Bound (CRLB) is a metric for optimizing quantitative sequences but requires an analytical expression for the signal. Magnetic Resonance Fingerprinting (MRF)¹ has a complex signal expression for the CRLB² making it difficult to account for B0 and B1 variations, or relevant signal contributions from diffusion and magnetization transfer. We apply Automatic Differentiation (AD) to the Bloch simulation and calculate the CRLB and its gradient with respect to sequence parameters for the MRF IR-FISP³ sequence in Figure 1. Reverse-mode AD (backpropagation) is applied when optimizing neural networks, which allows the gradient of the loss (a non-linear function of the weights), to be efficiently computed. AD allows us to perform the CRLB optimization with an asymptotic runtime of $$$O(N_{TR})$$$ and without deriving an expression for the signal.

Theory

The CRLB is the minimum achievable variance of a parameter estimate for an unbiased maximum likelihood estimator⁴. Assuming Gaussian noise, the CRLB of parameter $$$\theta_i$$$ is given by:

$$CRLB(\theta)_i=[V(θ)]_{ii}=[(\sum_{k=1}^K J_k^{}(\theta)J_k^T(\theta))^{-1}]_{ii},$$

where $$$V(\theta)$$$ is the CRLB matrix, and $$$J_k(\theta)$$$ is the Jacobian of echo $$$k$$$ with respect to the estimated parameters $$$\theta$$$. The relative CRLB (rCRLB) is defined as: $$rCRLB(\theta)_i=\frac{CRLB(\theta)_i}{\theta_i^2}.$$

The rCRLB objective represents the relative error of each parameter estimate and is non-convex requiring appropriate initialization.

AD applies the chain rule to calculate derivatives of multivariable functions⁵. A simple function for spin-echo imaging differentiated with AD is shown in Figure 2. Forward-mode AD is computationally efficient when computing the derivatives for a function $$$f\in R^n\to R^m,n<<m$$$. We employ forward-mode AD to calculate the Jacobian of the MRF Bloch simulation, which is efficient since $$$n=3$$$ for $$$\theta=[M0,T1,T2]$$$ and $$$m=N_{TR}$$$ (the signal at each echo). Reverse-mode AD is used to calculate the derivative of the rCRLBs with respect to sequence parameters for gradient-based optimization, illustrated in Figure 3.

Methods

Basic Bloch and EPG functions (available at http://web.stanford.edu/~pkllee) were implemented in autograd⁶. Flip angles (FAs) and TRs were chosen to optimize the rCRLB using the SciPy implementation of Sequential Quadratic Programming⁷. The sequence was simulated with Extended Phase Graphs (EPG) and optimized for grey matter at 3T (T1/T2=1330/80 ms) for $$$N_{TR}=400$$$. The optimization was initialized with the FA and TR schedule from Jiang et al.³ Constraints were: $$$10^{o}<\alpha<60^{o}$$$, $$$11<TR<16$$$ ms, with constraints on the FAs to promote smoothness in the magnetization evolution²: $$$|\alpha_i - \alpha_{i-1}|≤1^{o}$$$. A consenting volunteer was scanned at 3T under IRB approval using an 8-channel head coil. A fully-sampled MRF brain scan was performed once each for the rCRLB optimized sequence and the sequence used to initialize the optimization. A 10 second pause was added between interleaves to allow recovery to equilibrium. Scan parameters were: FOV$$$=$$$22$$$\times$$$22 cm, resolution$$$=$$$1.2 mm, 36 spiral interleaves, slice thickness$$$=$$$5 mm, duration 10:15. Reconstruction was done with a NUFFT and a dictionary T1=[200:4:2600], T2=[10:0.5:300] ms. Grey and white matter ROIs were manually chosen, and the relative errors defined as $$$\sigma_{ROI}/\mu_{ROI}$$$ were compared for the two acquisitions.

Results

Sequence parameters obtained by the rCRLB optimization and the grey matter fingerprint for the different schedules are plotted in Figure 4. The optimized FAs and TRs match the bang behavior observed by Zhao et al.²

Results from the in-vivo scan are shown in Figure 5. The relative error in the ROIs show that the T2 estimate is greatly improved while M0 and T1 variations remain similar in the rCRLB optimized sequence. This matches the prediction of the individual components in the rCRLB (Figure 4), where the largest improvement in rCRLB is for T2. Bias in the rCRLB optimized parameter maps may be due to unmodelled B1 slice variations which become significant at larger FAs⁸.

Discussion

Calculating the CRLB with AD is convenient since many MR phenomena are easy to simulate but cumbersome to differentiate analytically. This is exemplified by our use of EPG in the optimization, where the chain rule for applying spoiler gradients is easily determined with AD. Using AD to calculate the CRLB objective and gradient has $$$O(N_{TR})$$$ asymptotic runtime, which outperforms the $$$O(N_{TR}^2)$$$ scaling when approximating the gradient with finite differences. Finite differences has quadratic $$$N_{TR}$$$ scaling in MRF since $$$N_{TR}$$$ increases the number of variables to differentiate and the computation for each forward Bloch simulation. The AD framework has excellent flexibility and could be applied to MRF sequences with RF phase, a high number of TRs, or diffusion and magnetization transfer.

Conclusion

We have presented a framework for automatic differentiation of Bloch simulations and its application to MRF IR-FISP. We demonstrated that rCRLB optimization improves the accuracy of in-vivo parameter estimates, and that the optimization can be efficiently performed using automatic differentiation.

Acknowledgements

Grant support from GE Healthcare.

References

1. Ma D, et al. Magnetic resonance fingerprinting. Nature. 2013.

2. Zhao B, et al. Optimal Experiment Design for Magnetic Resonance Fingerprinting: Cramér-Rao Bound Meetse Spin Dynamics. arXiv. 2017.

3. Jiang Y, et al. MR fingerprinting using fast imaging with steady state precession (FISP) with spiral readout. Magnetic Resonance in Medicine. 2015.

4. Kay SM. Fundamentals of Statistical Signal Processing. Prentice Hall Signal Processing Series. 1993.

5. Baydin AG, et al. Automatic differentiation in machine learning: a survey. arXiv. 2015.

6. Maclaurin D. Modeling, Inference and Optimization with Composable Differentiable Procedures. PhD thesis, Harvard University. 2016.

7. Jones E, et al. Scipy: Open source scientific tools for Python. http://www.scipy.org/.

8. Buonincontri G, et al. MR fingerprinting with simultaneous B1 estimation. MRM. 2016.

Figures

Figure 1: Pulse diagrams for the MRF IR-FISP sequence. A single TR of the MRF IR-FISP acquisition is shown. The $$$\alpha_n$$$, and $$$TR_n$$$ are different for each TR which prevents a steady-state magnetization from forming. The rCRLB of the MRF sequence was optimized over $$$\alpha_n$$$ and $$$TR_n$$$ for $$$1 ≤ n ≤ N_{TR}$$$.

Figure 2: The derivative of a multiple-echo spin echo signal is calculated using backpropagation. The expression is split into elementary operations, each with a simple derivative that is applied during backpropagation. Performing backpropagation requires an evaluation of the graph in the forward direction. Evaluation of the forward direction is shown in the left table, creating temporary variables $$$U,V,W,X,Y$$$. Backpropagation is shown by green arrows in the graph. At each green arrow, the chain rule of the respective forward function is applied.

Figure 3: The automatic differentiation scheme used to obtain the CRLB and gradient of the CRLB with respect to sequence parameter is shown. a) A forward Bloch simulation is performed which generates the vector of echos $$$M_{echos}$$$ for a tissue $$$\theta$$$. b) Foward mode automatic differentiation (red) is used to calculate the Jacobian at each echo, which is assembled to obtain the CRLB cost function $$$C$$$. c) Reverse-mode automatic differentiation (blue) is performed through the graph operations generated in the forward mode calculation and the Bloch simulation to obtain the derivative of the cost function with respect to sequence parameters.

Figure 4: Top: Optimal FA and TR schedules from rCRLB minimization are compared to the initialization. The optimal schedules exhibit a bang nature taking on the minimum and maximum possible values. Inset: The loss curve shows that the optimization has converged and that the objective has been reduced by approximately 60% with the improvement being mainly from T2. The rCRLB optimization converged in 1.7 CPU hours which is significantly faster than the 19 CPU hours (per tissue) required in the implementation by Zhao et al. Bottom: The grey matter rCRLB optimized fingerprint has higher signal intensity compared to the initialization.

Figure 5: Top: Parameter maps show bias in the T1 and T2 estimate for the rCRLB optimized acquisition. This effect is most prominent in the white matter T2. Middle: Zoomed ROIs outlined in blue show that variation in the parameter estimate is reduced. Bottom: The rCRLB optimized sequence improves the relative error for parameter estimates in both grey and white matter. Even though grey matter (T1/T2 = 1330/80 ms) was the optimization target, improved parameter estimates in white matter (T1/T2 = 800/75 ms) are observed. ROI mean and standard deviations in the accompanying table demonstrate the bias between sequences.

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)

0812