3236

Parameter estimation with matrix-based signal models using VARPRO for transient- and steady-state imaging
Nam Gyun Lee1 and Krishna S. Nayak2
1Biomedical Engineering, University of Southern California, Los Angeles, CA, United States, 2Electrical and Computer Engineering, University of Southern California, Los Angeles, CA, United States

Synopsis

VARPRO-based parameter estimation is extensively used in MR for its improved accuracy, precision, and convergence behavior compared to general nonlinear least-squares algorithms. This study investigates the feasibility of using matrix-based signal models with VARPRO instead of conventional analytic signal equations. Simulations and in-vivo study show that VARPRO with matrix-based signal models is identical to VARPRO with analytic signal equations, and both VARPRO approaches provide enhanced precision and accuracy in relaxometry maps compared to the conventional DESPOT1/2 methods from variable flip angle SPGR and bSSFP measurements.

Introduction

The recently introduced matrix-based signal modeling approach1 is attractive because closed-form expressions based on the Bloch equations can model off-resonance, finite duration RF pulses, and magnetization transfer effects. However, parameters were estimated by a general nonlinear least-squares algorithm with many initial starting points, which may cause convergence issues. The variable projection (VARPRO) algorithm2,3 has been proposed to provide a better convergence for "separable" least-squares problems. The VARPRO approach has been used extensively in MR4,5,6 because of its accuracy and precision for a parametric fitting of the data with analytic signal models. Here, we investigate VARPRO-based parameter estimation using matrix-based signal models for transient- and steady-state quantitative imaging.

Theory

Jaynes' matrix formalism11: We assume a single-pool model. We assume that an RF pulse is applied on the x-axis with instantaneous pulse approximation. The temporal evolution of the magnetization M=[MxMyMz]T is governed by:
dM(t)dt=(Ω(t)+Λ)M(t)+CM0,
where
Ω(t)=[00γB1,y(t)00γB1,x(t)γB1,y(t)γB1,x(t)0],Λ=[R22πΔf02πΔfR2000R1],C=[00R1].
The matrix Ω(t) describes RF rotation, and Λ and C describe evolution in the absence of RF. A general solution to the Bloch equations in the absence of RF is described by M(t+τ)=SM(t)+(SI)Λ1CM0 where S=exp(Λτ)=E(τ)P(τ). The relaxation E(τ) and free precession P(τ) operators are defined as:
E(τ)=[eτR2000eτR2000eτR1],P(τ)=[cos(2πΔfτ)sin(2πΔfτ)0sin(2πΔfτ)cos(2πΔfτ)0001].

Steady-state imaging: In the steady state, we have RΨM(t+τ)=M(t)=MSS where R=exp(Ωτ) and Ψ accounts for perfect spoiling (SPGR) or π phase cycling (bSSFP). Using the steady-state condition, we can express the steady-state magnetiztion immediately after an RF pulse as: MSS=(IRΨS)1RΨ(SI)Λ1CM0.The matrix-based signal expressions for SPGR and bSSFP (at TE) are expressed as:
MSPGR=(IRΨSPGRE1P1)1RΨSPGR(E1P1I)Λ1CM0
MbSSFP=(E2P2(IRΨbSSFPE1P1)1RΨbSSFP(E1P1I)Λ1C+(E2P2I)Λ1C)M0,
where P1=P(TR),P2=P(TE),E1=E(TR),E2=E(TE), ΨSPGR=diag([0 0 1]), and ΨbSSFP=diag([-1 -1 1]).

VARPRO formulation: The VARPRO method efficiently solves nonlinear least-squares problems when a general nonlinear model can be formulated as
minαSα,cScyΦ(α)c22,
where Φ(α) is the real-valued system matrix parameterized by nonlinear parameters α and c real-valued vector of linear parameters. For given α, we estimate c=Φ(α)y and replace c in the original formulation as:
minαSαyΦ(α)c22=yΦ(α)Φ(α)y22=(IΦ(α)Φ(α))y22.
The subproblem with nonlinear parameters is solved with a nonlinear least-squares algorithm (e.g., Levenberg-Marquardt) with the Frechet derivative of the (IΦ(α)Φ(α))y. The differentiation of pseudo-inverses Φ(α) with respect to α is well described by Golub and Pereyra2. The VARPRO formulation updates (c(i),α(i)) until a stopping criteria is met. The MATLAB implementation of VARPRO8 (varpro.m) requires two inputs: Φ(α) and Φα. We describe the system matrices for the variable-flip-angle (VFA) SPGR, bSSFP, SPGR & bSSFP, and SPGR & bSSFP with off-resonance as:
Φ(α1)SPGR=[Ωy(MSPGR(θ1))/M0Ωy(MSPGR(θm1)/M0],Φ(α2)bSSFP=[Ωy(MbSSFP(θ1))/M0Ωy(MbSSFP(θm2))/M0],Φ(α3)S&b=[Ωy(MSPGR(θ1))/M0Ωy(MSPGR(θm1))/M0Ωy(MbSSFP(θ1))/M0Ωy(MbSSFP(θm2))/M0],Ψ(α4)S&b,Δf=[Ωy(MSPGR(θ1))/M0Ωy(MSPGR(θm1))/M0Ωx(MbSSFP(θ1))/M0Ωx(MbSSFP(θm2))/M0Ωy(MbSSFP(θ1))/M0Ωy(MbSSFP(θm2))/M0],
where nonlinear parameters for each case are α1=R1,α2=R2,α3=[R1R2], and α4=[R1R2Δf], Ωx() and Ωy() select the x and y components of the magnetization, respectively, and m1 and m2 is the number of variable flip angles for SPGR and bSSFP, respectively. Note that c=M0 for all cases and when model parameters were not included in α, fixed parameters were assumed (R1=R1,true and Δf=0 for ΨbSSFP(α)). The derivatives of Φ(α) with respect to each element of α is defined as:
Φα=[Φα1Φα2Φαp]Rm×p.
Each column can be calculated using the product rule and the fact that ddα[A(α)1]=A(α)1[ddαA(α)]A(α)1 and [ddαA(α)]ij=ddαAij(α).

Methods

For both simulations and in-vivo experiments, excitation flip angles for SPGR were: 1.5°, 2.2°, 3.23°, 4.74°, 6.96°, 10.22°, 15° and for bSSFP were: 15° to 65° (steps of 10°).
Noiseless 1D signal simulation: To test the convergence (accuracy) of the proposed method, we evaluated the proposed method on noiseless simulated datasets. Four test datasets were generated: 1) VFA SPGR, 2) VFA bSSFP, 3) VFA SPGR & bSSFP, and 4) VFA SPGR & bSSFP with off-resonance (30Hz). The datasets 1-3 and dataset 4 were generated with the Ernst and Freman-Hill equations, and a matrix-based signal model, respectively.
Noisy Brainweb simulation: We performed noisy simulation studies using a Brainweb "fuzzy" (partial volume) dataset (T1,T2,M0): We generated 1) VFA SPGR, 2) VFA SPGR & bSSFP.
Human VFA SPGR: A 3D whole-brain VFA SPGR dataset with TR = 5ms was acquired at 3T. Bloch-Siegert-based B1 maps were acquired for B1+ compensation9.
Parameter estimation: Datasets were fitted with VARPRO with matrix-based signal models (eigenvalue), VARPRO with analytic signal equations (analytic), and the conventional DESPOT1/2 method10. For both VARPRO approaches, the same initial starting points (1,1) and stopping criteria were used throughout different voxels.

Results and Discussion

Figure 1 shows results from noiseless simulated datasets; The eigenvalue approach converged to the ground truth within a few iterations (<6) for all test cases (See legend for estimated parameters). All methods fit the data well with near-zero residual, however, only the eigenvalue VARPRO approach had no bias in all cases. Note also that the eigenvalue VARPRO approach can simultaneously estimate off-resonance (Figure 1d). Figures 2 and 3 show results from noisy simulated VFA SPGR and VFA SPGR & bSSFP Brainweb datasets, respectively; Both VARPRO approaches provide T1/T2/M0 estimates with superior accuracy and precision than DESPOT1/2. Figure 4 shows estimated T1/M0 maps from one healthy human VFA SPGR dataset; T1 histograms have lower coefficient of variation with the proposed approach suggesting that they may be more accurate. A comparison with a reference method is future work.

Conclusion

We demonstrate and validate an eigenvalue VARPRO approach for parameter estimation, specifically T1/T2/M0 mapping. We show improved accuracy and precision compared to DESPOT1/2 and robustness to off-resonance compared to analytic VARPRO.

Acknowledgements

We gratefully acknowledge funding support from NIH R01-HL130494.

References

1. Malik SJ, A.G. Teixeira RP, West DJ, Wood TC, Hajnal JV. Steady-state imaging with inhomogeneous magnetization transfer contrast using multiband radiofrequency pulses. Magn. Reson. Med. 2019;00:1-15. doi:10/1002/mrm.27984

2. Golub GH, Pereyra V. The differentiation of pseudoinverses and nonlinear least squares problems whose variables separate. SIAM Journal on Numerical Analysis. 1973;10:413–432.

3. Golub GH, Pereyra V. Separable nonlinear least squares: The variable projection method and its applications. Inverse Prob. 2003;19:R1–R26.

4. Haldar JP, Anderson J, Sun SW. Maximum likelihood estimation of T1 relaxation parameters using VARRO. In Proceedings of the 15thAnnual Meeting of ISMRM, Berlin, Germany, 2007 (Abstract 41).

5. Hernando D, Haldar JP, Sutton BP, Ma J, Kellman P, Liang ZP. Joint estimation of water/fat images and field inhomogeneity map. Magn Reson Med. 2008;59:571–580.

6. Zhao B, Lu W, Hitchens TK, Lam F, Ho C, Liang ZP. Accelerated MR parameter mapping with low-rank and sparsity constraints. Magn Reson Med. 2015;74:489–498.

7. A.G. Teixeira RP, Malik SJ, Hajnal JV. Fast quantitative MRI using controlled saturation magnetization transfer. Magn. Reson. Med. 2018;1–14. doi:10.1002/mrm.27442

8. O’Leary DP, Rust BW. Variable projection for nonlinear least squares problems. Comput. Optim. Appl. 2013;54(3):579-593.

9. Sacolick LI, Wiesinger F, Hancu I, Vogel MW. B1 mapping byBloch–Siegert shift. Magn Reson Med 2010;63:1315–1322.

10. Deoni S, Rutt BK, Peters TM. Rapid combined T1 and T2 mapping using gradient recalled acquisition in the steady state. Magn Reson Med. 2003;49:515‐526.

11. Jaynes ET. Matrix treatment of nuclear induction. Phys. Rev. 1955;98:1099–1105 doi:10.1103/PhysRev.98.1099.

12. Gloor M, Scheffler K, Bieri O. Quantitative magnetization transfer imaging using balanced SSFP. Magn Reson Med 2008;60:691–700.

Figures

Figure 1. Evaluation of eigenvalue, analytic, and DESPOT1/2 on noiseless simulated datasets. Simulation parameters are taken from Gloor et al12 and shown in legend (Ground truth). (a) VFA SPGR. (b) VFA bSSFP. (c) VFA SPGR and bSSFP. (d) VFA SPGR and bSSFP with off-resonance (30 Hz). The eigenvalue approach is identical to the analytic approach when signal models are derived under the same conditions (no off-resonance and instantaneous RF). Notice that while all methods have near-zero residual, the analytic and DESPOT1/2 methods produce up to 15% error in the noiseless case.

Figure 2. (a) T1 maps and M0 maps estimated by VARPRO (eigenvalue), VARPRO (analytic), and DESPOT1 from noisy simulated VFA SPGR images. Normalized voxelwise error maps are also shown. The overall normalized root-mean-square error (NRMSE) calculated within the brain is indicated in red. (b) White matter (WM) and gray matter (GM) T1 histograms. (c) WM and GM M0 histograms. The analytic approach resulted in the same parameter estimates as the eigenvalue approach and hence omitted.


Figure 3. (a) T1/T2/M0 maps estimated by VARPRO (eigenvalue), VARPRO (analytic), and DESPOT1/2 from noisy simulated VFA SPGR & bSSFP images. (b) WM and GM T1 histograms. (c) WM and GM T2 histograms. (d) WM and GM M0 histograms. The eigenvalue approach resulted in the same parameter estimates as the analytic approach. Both VARPRO approaches utilize VFA SPGR and bSSFP data simultaneously for relaxometry estimation and drastically improves accuracy and precision.

Figure 4. (a) T1 maps and M0 maps estimated by VARPRO (eigenvalue), VARPRO (analytic), and DESPOT1 from human VFA SPGR images. (b) WM, GM, and CSF T1 histograms. (c) WM, GM, and CSF M0 histograms. The eigenvalue approach resulted in the same parameter estimates as the analytic approach. The coefficient of variation (CV) for WM, GM, and CSF is indicated in legend.

Proc. Intl. Soc. Mag. Reson. Med. 28 (2020)
3236