Accelerated free‐breathing liver fat and R2* quantification using multi‐echo stack‐of‐radial MRI with motion‐resolved multidimensional regularized reconstruction: Initial retrospective evaluation

To improve image quality, mitigate quantification biases and variations for free‐breathing liver proton density fat fraction (PDFF) and R2*$$ {\mathrm{R}}_2^{\ast } $$ quantification accelerated by radial k‐space undersampling.


INTRODUCTION
Quantitative measurement and evaluation of tissue biomarkers using MRI is of great clinical interest, because changes in the value of quantitative biomarkers can reflect the tissue property changes caused by diseases, and potentially at the early disease stage compared with qualitative diagnosis using conventional MR images. 1 Additionally, recording and comparing the longitudinal change of biomarker values are valuable and necessary for therapeutic evaluation and monitoring. 1,2at and iron deposition in tissue are indications of many diseases compared with normal conditions.][5][6][7][8][9][10][11][12][13][14][15][16][17][18] However, breath-hold MRI is challenging in children, elderly patients, and patients with breath-holding difficulties.Three-dimensional stack-of-radial MRI has been shown to provide accurate liver PDFF quantification in adult and pediatric patients using free-breathing acquisitions. 19,202][23] However, after self-gating, only a portion of the acquired radial views are used to reconstruct images and maps, potentially leading to residual radial undersampling artifacts and reduced SNR.Reduced SNR could limit the use of parameter maps, such as for accurate quantification in the low PDFF regime in potential liver donors. 248][29] However, previous studies reported that discrepancies could exist for quantitative results such as flow measured with CS accelerated methods compared to those with reference non-CS methods. 30,31[32][33][34] Recently, a stack-of-radial MRI with extra-dimensional golden-angle radial sparse parallel approach, 35 which regularizes the reconstruction in multiple dimensions, has demonstrated improved image quality in applications that use undersampled radial data such as dynamic contrast-enhanced liver imaging 35 and liver T 1 mapping. 36To our best knowledge, no previous work has investigated the feasibility and potential influence of CS reconstruction with multidimensional regularization in liver PDFF and R * 2 quantification yet.Two previous studies used model-based reconstruction to directly estimate water, fat, R * 2 , as well as B 0 and other parameters in the signal models. 37,38These methods use different signal models than that of previous studies, [21][22][23] which reconstruct images but not quantitative maps, therefore requiring careful initialization due to the nonlinear dependence on R * 2 and B 0 , and typically have very long reconstruction times (e.g., 1 h 35 min on a server with two Intel Xeon Gold 6238 central processing units [CPUs] 37 and 4 h on the Tesla V100 graphics processing unit [GPU] 38 ).
The purpose of this study was to develop a multi-echo stack-of-radial MRI method using motion-resolved reconstruction based on self-gating and CS with multidimensional regularization to improve image quality, mitigate quantification biases and variations for free-breathing liver PDFF, and R * 2 quantification accelerated by radial k-space undersampling.The proposed method was retrospectively and prospectively validated in motion phantoms with reference results from acquisitions without motion, retrospectively compared in vivo with a reference breath-hold 3D Cartesian method, and evaluated in terms of image quality, quantification bias, and variation on PDFF and R * 2 maps.

Multidimensional regularization
With the motion-state dimension resolved by self-gating, undersampled stack-of-radial data are reconstructed by CS with multidimensional regularization to minimize the following cost function 39 : where A is a catch-all matrix collectively applying all relevant operations consisting of multiplication with coil sensitivity maps and nonuniform Fourier transformation; x is the 3D image matrix (2D For a 1D discrete signal S(i), wavelet decomposition produces low (L) and high (H) frequency components through With the definition For a discrete 2D spatial + 1D motion-state signal, the wavelet decomposition of each permutation corresponding to I, J, T ∈ {L, H} can be generalized to Generally, the regularization term is a sum of the  1 -norm of all decompositions.In principle, each permutation of eight wavelet decompositions, (LLL), (LLH), (LHL), (LHH), (HLL), (HLH), (HHL) and (HHH), may have a different regularization strength.In the spirit of prior literatures, 39,40 the regularization weightings in the spatial dimensions are set with previous empirical values, and there is only a distinction in the low and high frequency decomposition in the motion state dimension.Consequently, two regularization parameters were considered: where  is the regularization weighting of the low frequency component and β is the relative change of the regularization weighting for the high frequency component.Circular boundary conditions are used for the spatial dimensions and open-boundary conditions for the motion-state dimension.The overall optimization was performed based on the FISTA algorithm 41 using the Chambolle-Pock algorithm 42 to evaluate the proximal operator.

Acquisition and reconstruction pipeline of the proposed method
A free-breathing multi-echo GRE stack-of-radial imaging sequence [19][20][21][22] with golden-angle ordering 43,44 was used for the data acquisition.The same radial view angle and partition encoding step were sampled for all echoes.The acquired k-space data underwent data-processing steps including k-space shift correction, [19][20][21][22][23]45 self-gating, [21][22][23] multidimensional regularization, 35,44 and quantification of PDFF and R * 2 . 14This pipeline was implemented in C++ and integrated into the vendor-provided research application sequence and image reconstruction framework (Siemens Healthcare, Erlangen, Germany), capable of both inline reconstruction after acquisitions and retrospective offline reconstruction using raw data files. The ilementation could use GPUs for faster reconstruction compared with CPU-only reconstruction.

Image-reconstruction methods for comparison
All radial raw k-space data in this study were acquired using a multi-echo GRE stack-of-radial research application sequence with golden-angle ordering and k-space shift calibration [19][20][21][22][23] with oversampled or fully sampled number of radial views 28 and saved for offline processing.Several methods were used to reconstruct images for comparison.First, termed as the "average" method, all radial views were reconstructed without any motion-compensation techniques.2][23] Alternately, self-gating was used to resolve multiple respiratory motion states or bins, termed as the "multi-bin" method.The results of the multi-bin method were regularized in the spatial dimension only ( = 0), termed as the "multi-bin with spatial regularization" method.Finally, results of the "proposed" method were obtained by performing CS reconstruction on multibin data with multidimensional regularization ( ≠ 0) (i.e., in both spatial and motion-state dimensions).

Motion phantom data and processing
From a previous study, 22 stack-of-radial MRI data of a custom phantom scanned at 3 T (MAGNETOM Prisma; Siemens Healthcare, Erlangen, Germany) was processed for new analysis in this study.The phantom had a saline bag, a gel holder, and four vials to achieve a range of R * 2 values.Details and plots are found in Zhong et al. 22 During the acquisitions, a programmable motion stage was used to move the holder and four vials according to a prerecorded respiratory motion waveform.The saline bag remained stationary throughout the experiment.Imaging parameters are listed in Table 1.Data from an acquisition performed with the same parameters and with the phantom being stationary were processed as a reference. 22riginally acquired with 768 radial views (oversampled, corresponding to a sampling factor of 1.5), raw data were undersampled to 512, 256, 128, and 64 radial views, equivalent to sampling factors 1.0, 0.5, 0.25, and 0.125.
Data with different radial views were reconstructed with different methods for comparison.
Mono-exponential fitting was used to calculate R * 2 .Using ImageJ 1.52n (National Institutes of Health, Bethesda, MD, USA), regions of interest (ROIs) were placed inside the vials, holder, and saline bag on each of two slices, one near the vial tips and the other at the center location of the vials.R * 2 values within the ROIs at the first motion state were reported as mean ± SD and compared to reference values without motion using Bland-Altman analysis with mean difference (MD) and limits of Imaging parameters and time saving of the sequences and protocols used in this study.

Free-breathing stack-of-radial
Field strength (T) Abbreviations: CAIPIRINHA, 2D parallel imaging (controlled aliasing in parallel imaging results in higher acceleration) 46 ; N/A, not applicable.
agreement (LoA = MD ± 1.96 × SD), which were reported as [MD; lower LoA, upper LoA].Furthermore, to assess the viability of the proposed method on future acquisitions and to incorporate the PDFF data into the phantom experiment, a new phantom (Calimetrix, Madison, WI, USA) was prospectively scanned on the same 3T scanner, which included seven vials and another five vials presenting varying ranges of PDFF and R * 2 values.The sequences and parameters were the same, except the image matrix was 256 × 256 and in-plane pixel size was 1.6 × 1.6 mm 2 .The prerecorded respiratory motion waveform differed yet exhibited similarities with the previous one. 22

In vivo data and processing
This study was Health Insurance Portability and Accountability Act-compliant and approved by the local institutional review board.Data included in a previous study 21 were used for new analysis in this study: free-breathing stack-of-radial and breath-hold Cartesian data from 5 healthy subjects (1 female, 34.6 ± 4.8 years, body mass index 22.6 ± 3.5 kg/m 2 ) and 6 subjects with nonalcoholic fatty liver disease (3 females, 58.5 ± 9.5 years, body mass index 28.9 ± 2.7 kg/m 2 ) acquired at 3 T (MAGNETOM Prisma or Skyra; Siemens Healthcare, Erlangen, Germany).
Originally acquired with fully sampled radial views (404 or 352, sampling factor 1.0), raw data were undersampled by factors of 0.5, 0.25, and 0.125 and reconstructed with different methods for comparison.Data from a breath-hold multi-echo Cartesian GRE acquisition were used as a reference. 21Imaging parameters are listed in Table 1.
PDFF and R * 2 values within 12 ROIs on different slices 21 were measured at the first motion state (the end expiration state) and compared with reference values measured by breath-hold Cartesian using Bland-Altman plots.SD values within the ROIs were displayed as error bars.
To offer supplementary assessment, additional measurement and comparison of PDFF and R * 2 values were performed within a larger ROI on the upper-mid-level slice covering the whole liver and avoiding large vessels.

Statistical analysis of linear mixed-effects model fitting
Linear mixed-effects models were used to fit the phantom and in vivo data.Mean and SD of PDFF and R * 2 were analyzed to evaluate significant biases and motion artifact-induced variations not corrected by the reconstruction method, respectively.The coefficient of variation (CV) on the first echo image was evaluated as a metric collectively considering SNR and undersampling/motion artifacts.More details are found in Supporting Information.

Parameter optimization of the proposed method
Our recommended value for the motion-state regularization parameter β was 80, and our recommended value for the number of bins was 4. Readers can refer to Supporting Information for details.

Phantom results
Example echo images, R * 2 maps, and Bland-Altman plots of the phantom using data with 128 radial views (sampling factor = 0.25) and different reconstruction methods are shown in Figure 1.Improved quality of echo images and R * 2 maps was evident by the proposed method compared with the current stable40 method.MD and LoA for the proposed method were [−0.6; −13.7, 12.6] s −1 for R * 2 , similar to those of the current stable40 method: [0.2; −13.4,13.7] s −1 .
For the linear mixed-effects model fitting results, among all sampling factors, the proposed method had the smallest CV, corresponding to the best image quality (Figure S7).The average method produced biased R * 2 values that were at least 17.8 s −1 greater than the reference results without motion (p = 0.013), whereas all other methods showed no significant differences (p > 0.131).For the sampling factor of 0.125, all methods had significantly increased SD values, such as significant residual undersampling/motion artifacts on the R * 2 maps (p < 0.033), whereas for the sampling factors of 1.0, 0.5, and 0.25, the proposed method consistently had insignificant residual artifacts (p > 0.281).More details are found in Supporting Information.
Bland-Altman plots illustrating outcomes for prospectively acquired data of PDFF and R * 2 vials with 128 radial views (sampling factor = 0.25) are shown in Figure S5.MD and LoA for the proposed method were [2.1; −29.8, 34.0] s −1 for R * 2 and [0.5%; −1.4%, 2.4%] for PDFF, exhibiting high consistency with the reference results, while both the average and stable40 methods exhibited large MD and LoA values for R * 2 .

F I G U R E 1
Example images, R * 2 maps, and Bland-Altman plots of the phantom using data with 128 radial views (sampling factor = 0.25).The reference results are from the data with 768 radial views acquired when the phantom was not moving (the first column).All other methods show results at the first motion state from the data acquired when the phantom was moving (the second to the fourth columns).All Bland-Altman plots show comparison between the corresponding method and the reference for the R * 2 values in the regions of interest (ROIs; six red circles).Red arrows indicate the motion-induced artifacts in the images and maps.The window/level settings are the same for all echo images and R * 2 maps, respectively.The echo images of stable40 (the third column of the first two rows) have different image scaling of the reconstruction pipeline from the other echo images.The R * 2 values are reported as mean ± SD.The relatively large error of R * 2 in V1 (14.9 ± 16.5 s −1 ) with the proposed method was caused by a bubble in that vial.

In vivo results
Example echo images and R * 2 and PDFF maps of the same subject in Figures S3 and S4 using data with 101 radial views (sampling factor = 0.25) and a clinical subject using data with 88 radial views (sampling factor = 0.25) are shown in Figures 2 and 3, respectively.Improved quality of images and maps using the proposed method compared with the current stable40 method was also observed.
Linear mixed-effects model fitting showed that the proposed method either had the best (lowest) CV or was not significantly different from the best method (Figure S10), and especially had significantly lower CV than all other methods for sampling factors of 0.125 and 0.25 (p < 0.001).For R * 2 SD values, except for the proposed method and the average method, all other methods showed significantly higher SD values for all sampling factors (p < 0.021).For PDFF SD values, for sampling factors of 0.25 and greater, the proposed method and the average method were not significantly different from each other (p > 0.812), but together were significantly lower than all other methods (p < 0.002).Supporting Information includes more detailed results.Bland-Altman analysis of PDFF and R * 2 values within the large whole-liver ROI on the upper-mid-level slice are shown in Figure S6.Similarly, MD and LoA for the proposed method using 101 radial views were [−3.9; −21.9, 14.1] s −1 for R * 2 and [−1.0%; −3.0%, 1.1%] for PDFF, closer to those of the stable 40 method using 404 radial views: [−0.8; −20.7, 19.1] s −1 for R * 2 and [−0.4%; −2.5%, 1.6%] for PDFF.The error bars also had similar appearances to the 12 ROI results.

Acquisition-time savings and reconstruction time of the proposed method
Total times of all original acquisitions, as well as the equivalent acquisition times of the proposed method with different sampling factors, are listed in Table 1.For phantom experiments, a sampling factor of 0.25 was the empirically preferred acceleration factor, leading to an equivalent acquisition time of 55 s compared with 156 s for the fully sampled baseline protocol.For in vivo experiments, depending on the slightly different protocols, the sampling factor of 0.25 was also preferred and led to acquisition time of 64 s or 59 s, compared with 171 s or 153 s.
The proposed method was implemented and tested on a desktop computer with an 8-Core i7-6700 CPU at

DISCUSSION
In this study, a multi-echo stack-of-radial MRI method using motion-resolved reconstruction with self-gating and CS reconstruction with multidimensional regularization for free-breathing liver PDFF and R * 2 quantification was developed.The method was validated in motion phantoms with reference results from acquisitions without motion, as well as in vivo with reference results from breath-hold 3D Cartesian acquisitions.
The regularization in this study used redundant Haar wavelets, which provide a uniform and flexible regularization processing in both spatial and motion state dimensions (Eqs.[2]- [4]) and an image representation allowing for the analysis of features at different scales (although this study used only one level).8][49] Moreover, the redundant Haar transform can be efficiently implemented in parallel computing architectures with multicore processors and GPUs, allowing for efficient image reconstruction.Finally, it is based on the same algorithm used for the commercially available GRASP application in previous studies. 39,40or phantom experiments, the proposed method supported to acquire the radial views with the recommended sampling factor of 0.25 (i.e., a 4-fold acceleration for the radial k-space acquisition), which saved 101 s in the equivalent total acquisition time, or 64.7% shorter compared with the fully sampled protocol with the sampling factor of 1.0.For in vivo experiments, acquiring the number of radial views with the recommended sampling factor of 0.25 could save 107 s (62.6% shorter) or 94 s (61.4% shorter), compared with the originally fully sampled protocols.This is a considerable time-savings for free-breathing stack-of-radial MRI (typically 2-5 min).
Our previous study used self-gating to accept data near end-expiration for reconstruction suffering from undersampling artifacts and influenced quantification accuracy if undersampled acquisitions are performed. 21Previous model-based reconstruction studies 37,38 directly estimated water, fat, R * 2 , and B 0 in signal models for optimization of CS reconstruction and regularization, leading to long reconstruction time (>1 h in Schneider et al. 37 and 4 h in Tan et al. 38 ).Our proposed method uses CS reconstruction and regularization to reconstruct echo images with improved image quality and reduced biases and variation compared with our previous study, 21 and allows the flexibility of combining an efficient image reconstruction with different PDFF and R * 2 quantification algorithms using echo images as the input.The performance of redundant Haar wavelet decomposition used in our proposed method is expected to be similar to the original extra-dimensional golden-angle radial sparse parallel method with total variation.The proposed method is suitable for both clinical and research settings, supporting online and offline reconstructions.

CONCLUSION
A multi-echo stack-of-radial MRI method using motion-resolved CS reconstruction and multidimensional regularization for free-breathing liver PDFF and R * 2 quantification was developed.Results demonstrated good agreement of PDFF and R * 2 measured by the proposed method compared with the reference results in motion phantoms and in vivo.Equivalent or less PDFF and R * 2 quantification bias and variations, as well as improved image quality by the proposed method, were observed compared with a self-gating method in previous studies.An accelerated in vivo imaging protocol with a radial view sampling factor of 0.25 and the proposed reconstruction method were validated retrospectively, which achieved an acceleration factor of 4 and an approximate acquisition time saving of 107 s (62.6% shorter) or 94 s (61.4% shorter).This proposed method may allow accelerated free-breathing liver PDFF and R * 2 mapping with reduced biases and variations.

F I G U R E 2
Example echo images and R * 2 and proton density fat fraction (PDFF) maps of the same healthy subject in Figures S3 and S4 using data with 101 radial views (sampling factor = 0.25).The reference results are from the breath-hold Cartesian acquisitions (the first column).All other methods show results at the first motion state (end-expiration state) from the data with free-breathing acquisitions (the second to the fourth columns).Red arrows indicate the motion-induced artifacts in the images and maps.The window/level settings are the same for all echo images and R * 2 and PDFF maps respectively.The echo images of stable40 (the third column of the first two rows) have different image scaling of the reconstruction pipeline from the other echo images.The images and maps are cropped for better visualization.The R * 2 and PDFF values within three regions of interest (ROIs; red, green, and blue circles) are reported as mean ± SD.

F I G U R E 3
Example echo images and R * 2 and proton density fat fraction (PDFF) maps of a clinical subject using data with 88 radial views (sampling factor = 0.25).The reference results are from the breath-hold Cartesian acquisitions (the first column).All other methods show results at the first motion state (end-expiration state) from the data with free-breathing acquisitions (the second to the fourth columns).Red arrows indicate the motion-induced artifacts in the images and maps.The window/level settings are the same for all echo images and R * 2 and PDFF maps, respectively.The echo images of stable40 (the third column of the first two rows) have different image scaling of the reconstruction pipeline from the other echo images.The images and maps are cropped for better visualization.The R * 2 and PDFF values within four regions of interest (ROIs; red, green, blue, and purple circles) are reported as mean ± SD.

F I G U R E 4
Bland-Altman plots comparing the R * 2 and proton density fat fraction (PDFF) results from different methods versus those from the reference breath-hold Cartesian acquisitions for data with 404 (A) and 101 radial views (B), respectively.The difference was obtained by results from free-breathing stack-of-radial acquisitions minus results from the reference breath-hold Cartesian acquisitions.Different colors represent the values within 12 regions of interest (ROIs) of different subjects.The error bars are the SD of the measured R * 2 or PDFF values within the ROI of each method.

3. 4
GHz (Intel, Santa Clara, CA, USA), 64 GB RAM, and a Quadro P2200 graphics card with 5 GB memory (NVIDIA, Santa Clara, CA, USA), which took 305 s to generate results for the in vivo free-breathing imaging protocol with 101 radial views.The CPU-only reconstruction required 493 s for the same data set on a vendor-specific simulator computer (Intel 16-Core Xeon Gold 6154 CPU at 3.0 GHz, 128 GB RAM, no GPU).With a 24-Core Intel Xeon E5-2680 CPU at 2.5 GHz, 256 GB RAM and two NVIDIA M60 GPUs each with 8 GB GPU memory on a scanner (MAGNETOM Vida; Siemens Healthcare, Erlangen, Germany), the inline reconstruction time of the same protocol was 248 s.