Free-Breathing Water, Fat, R_2^{\star} and B_0 Field Mapping of the Liver Using Multi-Echo Radial FLASH and Regularized Model-based Reconstruction (MERLOT)
Zhengguo Tan, Sebastian Rosenzweig, Xiaoqing Wang, Nick Scholand, H. Christian M. Holme, Moritz Blumenthal, Dirk Voit, Jens Frahm, Martin Uecker
FFree-Breathing Water, Fat, R ∗ and B FieldMapping of the Liver Using Multi-Echo RadialFLASH and Regularized Model-basedReconstruction (MERLOT)
Zhengguo Tan , Sebastian Rosenzweig , Xiaoqing Wang , NickScholand , H Christian M Holme , Moritz Blumenthal , DirkVoit , Jens Frahm , and Martin Uecker Institute for Diagnostic and Interventional Radiology, University Medical CenterG¨ottingen, G¨ottingen, Germany German Center for Cardiovascular Research (DZHK), partner site G¨ottingen, G¨ottingen,Germany Biomedizinische NMR, Max-Planck-Institut f¨ur biophysikalische Chemie, G¨ottingen,Germany Cluster of Excellence “Multiscale Bioimaging: from Molecular Machines to Networks ofExcitable Cells“ (MBExC), University of G¨ottingen, G¨ottingen, Germany Campus Institute Data Science (CIDAS), University of G¨ottingen, G¨ottingen, Germany
January 11, 2021
AbstractPurpose:
To achieve free-breathing quantitative fat and R ∗ mapping ofthe liver using a generalized model-based iterative reconstruction, dubbedas MERLOT. Methods:
For acquisition, we use a multi-echo radial FLASH sequencethat acquires multiple echoes with different complementary radial spokeencodings. We investigate real-time single-slice and volumetric multi-echoradial FLASH acquisition. For the latter, the sampling scheme is extendedto a volumetric stack-of-stars acquisition. Model-based reconstructionbased on generalized nonlinear inversion is used to jointly estimate water,fat, R ∗ , B field inhomogeneity, and coil sensitivity maps from the multi-coil multi-echo radial spokes. Spatial smoothness regularization is appliedonto the B field and coil sensitivity maps, whereas joint sparsity regular-ization is employed for the other parameter maps. The method integratescalibration-less parallel imaging and compressed sensing and was imple-mented in BART. For the volumetric acquisition, the respiratory motionis resolved with self-gating using SSA-FARY. The quantitative accuracy Submitted to Magnetic Resonance in Medicine. a r X i v : . [ phy s i c s . m e d - ph ] J a n f the proposed method was validated via numerical simulation, the NISTphantom, a water/fat phantom, and in in-vivo liver studies. Results:
For real-time acquisition, the proposed model-based reconstruc-tion allowed acquisition of dynamic liver fat fraction and R ∗ maps at atemporal resolution of 0 . Conclusion:
The proposed multi-echo radial sampling sequence achievesfast k -space coverage and is robust to motion. The proposed model-basedreconstruction yields spatially and temporally resolved liver fat fraction, R ∗ and B field maps at high undersampling factor and with volume cov-erage. Keywords: model-based reconstruction, calibrationless parallel imaging,compressed sensing, water/fat separation, R ∗ relaxometry, multi-echo ra-dial sampling Introduction
Quantitative parameter mapping in the liver is of great interest in basic re-search and clinical practice. Specifically, quantitative proton density fat fraction(PDFF) and R ∗ maps have been shown to be non-invasive imaging biomark-ers for hepatic steatosis [1, 2] and iron overload [3, 4], respectively. Originat-ing from the 2-point chemical-shift-encoded Dixon method [5] for qualitativewater/fat separation, quantitative assessment of liver fat and iron decomposi-tion can be achieved via multi-point chemical-shift encoding (e.g. low flip anglemulti-gradient-echo acquisition) and physics-based modeling [6–13].For motion-robust and time-resolved quantitative mapping of the liver twoimportant techniques were proposed. First, radial sampling is motion robustand well suited for dynamic imaging [14–18]. Second, self-gating techniques havebeen introduced for liver fat and R ∗ quantification to correct for motion duringfree-breathing volumetric acquisition [19–21]. For a multi-echo stack-of-starsvolumetric acquisition, the long echo-train and volume coverage require longscan times of about 5 min. Therefore, real-time acquisition becomes infeasibleand self-gating is required to resolve respiration phases in free-breathing scans.These techniques, however, require about 800 radial spokes per echo and perpartition to balance volumetric spatial resolution and temporal resolution [19,21].To further accelerate data acquisition, undersampled multi-echo samplingcan be used with model-based reconstruction [20, 22]. First, the undersam-pling pattern differs amongst echoes and/or partitions in order to provide inco-herent sampling pattern as well as complementary k -space information in theparameter-encoding and/or partition directions. Second, model-based recon-struction integrates MR physics-based modeling into the forward model, directlyestimating quantitative parameter maps from the undersampled multi-echo k -space data [23–25]. Note that model-based reconstruction takes complementarysamples from all echoes into account, and thus enables high undersamplingfactors. Moreover, the combination of model-based reconstruction and parallel-imaging compressed-sensing reconstruction [26, 27] has achieved up to 7-fold ac-celerated whole liver fat and R ∗ mapping in a single breath hold (about 20 sec)[22]. To void the time constraints imposed by breath holding, Schneider etal. [20] used self-gating for free-breathing whole liver fat and R ∗ mapping with400 spokes per echo and per partition using model-based reconstruction. Thistechniques still requires pre-calibrated coil sensitivity maps and a pre-estimated B field map [28],Based on the real-time MRI framework [29], our previous work achievedthree-point qualitative water/fat separation via triple-echo radial FLASH ac-quisition and model-based reconstruction [30]. However, quantitative fat and R ∗ mapping was not yet achievable. Therefore, in this work we developed adedicated model-based reconstruction technique to jointly estimate water, fat, R ∗ , B field inhomogeneity maps as well as coil sensitivity maps from multi-echo radial FLASH sampled data. This offers a calibration-less technique which3 FGzGxGy
Figure 1: (Left) One representative TR block of the proposed multi-echo ra-dial FLASH sequence. (Right) The corresponding k -space trajectory. In thisexample, 9 echoes with different k -space spokes are acquired per RF excitation,whereas partial Fourier sampling is not used. The echoes are color coded, indi-cating the period when ADC is switched on, while the white solid lines indicateeither the ramp or the blip gradients. Note that spoiler gradients are ignoredin this illustrationcombines both parallel imaging and parameter mapping into a single nonlinearinverse problem. Exploiting the physics of multi-echo gradient-echo acquisition,four different signal models were implemented in BART [31] as nonlinear oper-ators, which allow (1) 3-point water/fat separation and B field mapping [30],(2) common R ∗ mapping between water and fat, (3) independent R ∗ mappingbetween water and fat, and (4) R ∗ mapping with only one dominant species.Each of these nonlinear operators can be combined with the nonlinear parallelimaging operator as described previously [29]. Auto differentiation implementedin BART then automatically yields the derivative required for the iterative non-linear optimization. Furthermore, this work combined the multi-echo radialsampling scheme with the stack-of-stars volumetric acquisition and the SSA-FARY self-gating technique [32] for whole liver coverage during free breathing. Data acquisition is based on a previously described triple-echo multi-spoke radialFLASH sequence [30]. Here, the echo train length (ETL) is extended up to 7gradient echoes per RF excitation, as depicted in Figure 1. The elongated ETLcaptures more of the T ∗ signal decay, and thus allows mapping of quantitative T ∗ values.The proposed multi-echo radial sampling sequence, in line with radial EPI[33], acquires multiple gradient echoes with different radial spoke encoding per4F excitation, thereby achieving rapid and efficient k -space coverage with im-proved temporal incoherence. Similar to our previous work [30], all spokeswithin one frame were uniformly distributed in k -space and the angle incrementbetween frames was the small Golden angle ( ≈ . o ). Multi-echo radial sam-pling is resistant to spatial distortions induced by B field inhomogeneities andimmune to motion artifacts (e.g. ghosting in Cartesian sampling). Parallel MRI [34–36] simultaneously receives signals from multiple receiver coils,and is extendable to include multiple echoes when using long echo-train MRIsequences, y j,m ( t ) = (cid:90) d (cid:126)r e − i π(cid:126)k ( t ) · (cid:126)r c j ( (cid:126)r ) ρ m ( (cid:126)r ) , (1)with c j and ρ m being the j th coil sensitivity map and the m th echo image,respectively. In the case of gradient echoes, ρ m is governed by ρ m = (cid:18) (cid:88) i I i · e i πf i TE m · e − R ∗ i TE m (cid:19) · e i πf B TE m , (2)where the first term sums up signals from all chemical species (indexed by i ),characterized by their corresponding proton density ( I i ), resonance frequency( f i ) and relaxation rate ( R ∗ i ). Here, the dependency on the spatial coordinates (cid:126)r is suppressed for simplicity. In addition, the echo signal is modulated by the B field inhomogeneity. TE m denotes the m th echo time.This generalized multi-species signal can be simplified to only two compart-ments [6, 7, 37], i.e. water (W) and fat (F), especially for liver and cardiacimaging, ρ m = (cid:18) W · e − R ∗ TE m + F · z m · e − R ∗ TE m (cid:19) · e i πf B TE m . (3)The chemical-shift phase modulation from fat is denoted as z m with the 6-peakfat spectrum [9, 38], while all fat peaks are assumed to have an equal R ∗ [39].W and F are complex-valued, while the other parameters in x are constrainedto be real.Noteworthy, in cases of relatively low fat fraction [40], the independent- R ∗ model in Equation (3) can be simplified as a single- R ∗ model proposed by Yuet al. [8] ρ m = (cid:18) W + F · z m (cid:19) · e − R ∗ TE m · e i πf B TE m , (4)which assumes common R ∗ between water and fat. Due to its numerical stabilityand simplicity, this single- R ∗ model has been used in numerous studies [9, 11,13, 19, 20, 41]. Moreover, in cases with only water protons, the model reducesto ρ m = ρ · e − R ∗ TE m · e i πf B TE m . (5)5ll the above-mentioned signal models, as well as the triple-echo water/fat sep-aration model [30], were implemented in BART [31].Given the above MR signal models, the nonlinear forward model in multi-coilmulti-echo acquisition can be written in the operator form y j,m = F j,m ( x ) := P m F M SB , (6)with x = (W , R ∗ , F , R ∗ , f B , c , · · · , c N ) T . j is the the coil index ( j ∈ [1 , N ]),and m the echo index ( m ∈ [1 , E ]). The nonlinear operator ( B ) calculates echoimages according to the parameter maps in x and the corresponding signalmodel as given in Equations (3) to (5). Every echo image is then point-wisemultiplied by a set of coil sensitivity maps in x , as denoted by the operator S . Afterwards, all multi-echo coil images are masked to be restricted to a givenFOV ( M ), Fourier-transformed ( F ), and sampled ( P ) at each echo. In addition,A k -space filter [42] was applied onto the sampling pattern P . The joint estimation of the unknown x , i.e. simultaneous and quantitative map-ping of water/fat-separated R ∗ , field f B , as well as the coil sensitivity maps,is a nonlinear inverse problem. A solution can be estimated using a regularizedleast-squares problem minimize (cid:107) y − F ( x ) (cid:107) + λR ( x )subject to R ∗ ≥ , R ∗ ≥ y is the measured multi-coil multi-echo k -space data, and R ( x ) is the reg-ularization applied onto x with strength λ ( λ > (cid:96) R ∗ , F, and R ∗ , a non-negativityconstraint onto the relaxation rates R ∗ and R ∗ , and spatial smoothness ofthe field inhomogeneity ( f B ) and coil sensitivity ( c i ) maps via a Sobolev-norm[29, 30]. In addition, (cid:96) x in Equa-tion (7) can be extended to include all temporal frames. In addition to theaforementioned spatial regularization, temporal total variation (TV) regular-ization [43] was applied onto all parameter maps except coil sensitivity maps.The objectve functional in Equation (7) was solved by the iteratively regular-ized Gauss-Newton method (IRGNM). In each Newton iteration, the nonlinearobjective function is linearized and solved by the alternating direction methodof multipliers (ADMM) [44]. For details about this algorithm please refer toAppendix. 6 Methods
The proposed reconstruction method was validated via a numerical phantomand phantom experiments at 3 T (Skyra, Siemens Healthineers, Erlangen, Ger-many).We used an analytical numerical phantom [45] implemented in BART tosimulate k -space data for ten T ∗ tubes with R ∗ values equally ranging from 5 . . − . Base resolution was 200. Seven echoes per excitation were usedwith both first TE and echo spacing set to 1 . − was added to the data. The proposed model-basedreconstruction was tested using 101 and 33 radial spokes per echo.The NIST phantom (System Standard Model 130, CaliberMRI, Boulder,Colorado, USA) was utilized for quantitative validation with a 20-channel headcoil. Detailed multi-echo radial FLASH acquisition parameters were 2D singleslice with field of view (FOV) 280 mm, base resolution 256, slice thickness 5 mm,bandwidth 890 Hz pixel − , and 7 echoes (TEs: 1 . , . , . , . , . , . , .
40 ms and TR 11 .
70 ms). For comparison, fully-sampled single-sliceCartesian multi-gradient-echo data using bipolar readouts was acquired withTE 1 . , . , . , . , . , . , .
16 ms and TR 22 ms. The otherparameters were the same as the radial acquisition.Furthermore, a simple water/fat phantom was constructed to validate esti-mation of fat fraction maps. Details of this phantom are provided in Support-ing Information Figure S1. Validation experiments were conducted with a 18-channel body matrix coil together with a spine coil. Fully-sampled single-sliceCartesian multi-gradient-echo data was acquired using the following parame-ters: FOV 192 mm, base resolution 192, bandwidth 840 Hz pixel − , flip angle 5 ° ,bipolar readouts with TEs 1 . , . , . , . , . , . , .
64 ms and TR293 ms. In addition, the proposed multi-echo radial FLASH sequence was usedwith TEs 1 . , . , . , . , . , . , .
90 ms and TR 12 .
30 ms, while theother parameters were the same as the Cartesian acquisition.The acquired Cartesian data was Fourier transformed and then coil combinedusing coil sensitivity maps estimated via ESPIRiT [46]. Afterwards, a pixel-wise fitting routine implemented in BART was used to reconstruct parametermaps using the model in Equation (5) for the NIST phantom, while for thewater/fat phantom the graph-cut reconstruction [28] was employed. In contrast,the acquired multi-echo radial data was reconstruction by the proposed model-based reconstruction method.
For in vivo scans, an 18-channel body matrix coil and a spine coil were employedfor parallel acquisition. The FOV was chosen as 320 ×
320 mm and voxel sizeas 1 . × . × , resulting in an image matrix size of 200 × ° and the bandwidth was 1090 Hz pixel − . The ETL was 7, with7Es 1 . , . , . , . , . , . , .
69 ms and a TR of 9 .
89 ms. Totest the repeatability of the proposed method, the scans were performed twice.Several subjects with no known illness participated in the development of thismethod. All volunteers gave written informed consent before MRI.
Single-Slice Acquisition
For real-time single-slice acquisition, one k -space frame consisted of 33 RF exci-tations. Thus, the acquisition time per frame was 326 ms for 231 radial spokes(33 excitation × R = 200 × ( π/ / ≈
10. Three slices were acquired for each scan of the2D acquisition.
Stack-of-Stars Volumetric Acquisition
For stack-of-stars volumetric acquisition [18], one TR block is repeated for everypartition until it changes readout gradients to sample spokes at different angularpositions. In this work, the total number of partitions was 32 with 12 . × N partition = 9 .
89 ms × ≈
356 ms). Forbreath-hold scans, 45 excitations per frame were used, leading to a total scantime of TA = T × N excitation = 16 s. For free-breathing acquisitions, a total of330 excitations were used, and thus the total scan time was 1 .
95 min.
The proposed regularized nonlinear model-based reconstruction has been imple-mented in BART [31], building on previous work on model-based T1 mapping[47, 48]. A brief description of the reconstruction procedure is given here.
Pre-Processing
The acquired multi-coil multi-echo data was compressed to 10 virtual coils viaprincipal component analysis [49]. The multi-echo sampling trajectory was cor-rected for gradient delays using the RING method [50]. Note that RING wasapplied onto every echo to estimate its corresponding gradient delay coefficients,which were then used to correct for trajectory.
Binning
As multi-echo data was continuously sampled without any physiological gating,a natural image reconstruction method was to reconstruct the dynamic datasequentially. That is, every frame was reconstructed by iteratively minimizingEquation (7). This non-gated reconstruction method resolved respiratorymotion without the need for gating techniques.8lthough the non-gated reconstruction can provide motion-resolved dynamicparameter maps, it is worthwhile to exploit the periodicity of physiological mo-tion using a self-gating reconstruction of the continuously acquired data.First, self gating dramatically reduces the amount of frames to be reconstructed.Second, it becomes necessary to apply self-gating techniques for volumetric free-breathing acquisition, because the temporal footprint of one frame becomes toolarge for a non-gated reconstruction. In this work, the recently developed SSA-FARY technique [32] was adapted, where the continuously acquired multi-echodata spanning multiple breathing cycles is binned into one cycle consisting ofseven respiration phases. As one echo-train readout was relatively short, onlythe first echo was extracted and used for SSA-FARY.
Initialization
Off-resonance phase modulation (refer to the second term in Equation (3))causes phase wrapping along the echoes, especially in cases of long echo trainand large B field inhomogeneity. Consequently, multiple local minimum couldoccur for the field map f B . To prevent this, W, F and f B maps were initial-ized by their estimates from a model-based 3-point water/fat separation [30]. R ∗ and coil sensitivity maps were initialized as 0.For 3D data acquired via multi-echo radial stack-of-stars FLASH, initializa-tion was conducted via applying temporal regularization [51] along the partitiondimension. Afterwards, the model-based reconstruction was performed individ-ually on every partition without temporal regularization. Iterative Reconstruction
For the model-based reconstruction, the regularization strength in Equation (7)is reduced along Newton iterations: λ n = 1 /D n − , with n being the n th Newtoniteration and the reduction factor D >
1. In this work, D = 3 was used. ForADMM, the maximum number of iterations in each Newton iteration were givenas: n maxiter = min( M, × − ln λ n ), where M can be specified by the user (set as100 by default). Consequently, the maximal iterations gradually increase alongNewton steps. The ADMM penalty parameter ρ was set as 0 . . (A) 101 spokes per echo(B) 33 spokes per echo B0 R2* Figure 2: Model-based reconstruction with the model in Equation (5) as wellas the quantitative analysis on a numerical phantom acquired by (A) 101 and(B) 33 spokes per echo, respectively. Displayed images are the ρ , R ∗ and B field maps. The right panel shows the mean and standard deviation values ofeach selected ROI (colored in the R ∗ map in (A)), with the y axis labeling thereference R ∗ values Numerical Simulation
Figure 2 shows the model-based reconstruction and quantitative analysis resultsof the numerical phantom. All displayed maps were reconstructed via jointlyestimating all parameters in Equation (5). The model-based reconstructionwas capable of recovering accurate quantitative R ∗ maps with R ∗ tubes equallyranging from 5 . . − . When reducing the number of radial spokes perecho from 101 to 33, the standard deviation of tubes with small R ∗ values wasslightly increased, but quantitative accuracy stayed consistent. NIST Phantom
Figure 3 compares the reconstructed R ∗ maps via pixel-wise fitting of the Carte-sian multi-gradient-echo and model-based reconstruction of the radial multi-echo data (MERLOT), respectively. Both reconstructions used the model inEquation (5). Due to the non-convexity of the signal model in Equation (5),only the absolute part of the complex echo images from Cartesian sampling wasused in the pixel-wise fitting. On the contrary, the spatial smoothness constrainton the B field inhomogeneity map in MERLOT prevents spatial discontinuity10 C a r t e s i an + m oba f i t M E R L O T R2* B0
Figure 3: Comparison of R ∗ maps from (top left) pixel-wise fitting of Cartesianmulti-gradient-echo images and (bottom left) MERLOT (i.e. multi-echo radialFLASH with model-based reconstruction). Both reconstructions used the modelin Equation (5). The right panel displays the correlation plot between thereference R ∗ (obtained via Cartesian + pixel-wise fitting) and the MERLOT R ∗ values based on the seven depicted region of interests (ROI)in the B map, thus the complex k -space data was used. Quantitative analysisof the seven selected ROIs shows good match between these two methods, whilethe MERLOT approach shows lower standard deviation. Water/Fat Phantom
Figure 4 shows the reconstructed water, fat and fat fraction maps from graph-cut reconstruction [28] of the multi-echo Cartesian data and model-based re-construction of the multi-echo radial data, respectively. Quantitative analysisof the computed fat fraction was provided in the right. Both reconstructionmethods provided good water and fat images. Fat fraction values showed goodmatch between the two reconstructions and to the product description.11 C a r t e s i an + g r aph c u t M E R L O T Figure 4: Comparison of fat fraction maps from (top left) graph-cut reconstruc-tion of Cartesian multi-gradient-echo images and (bottom left) MERLOT withthe independent- R ∗ signal model. The right panel displays the correlation plotbetween the reference fat fraction (obtained via Cartesian + graph cut) and theMERLOT fat fraction values based on the four depicted ROIs I ndependen t - R * S i ng l e - R * Water Fat Fat FractionB0R2* of waterR2* of fatR2*
Figure 5: Model-based reconstruction of the real-time free-breathing liver radialmulti-echo FLASH acquisition. Displayed images are results based on (top) thesingle- R ∗ signal model in Equation (4), and (bottom) the independent- R ∗ signalmodel in Equation (3) 12able 1: Quantitative analysis of liver R ∗ and fat fraction from different acqui-sition protocols. The ROIs shown in Figure 5 were used for analysis.2D FB
3D BH
3D FB WFR2S WF2R2SScan 1 R ∗ (s − ) ROI 1 50 . ± .
95 48 . ± .
88 59 . ± .
53 53 . ± . . ± .
13 45 . ± .
09 47 . ± .
01 39 . ± . . ± .
31 47 . ± .
18 47 . ± .
39 46 . ± . . ± .
51 7 . ± .
56 8 . ± .
84 8 . ± . . ± .
87 10 . ± .
91 13 . ± .
98 13 . ± . . ± .
38 10 . ± .
35 12 . ± .
99 13 . ± . R ∗ (s − ) ROI 1 48 . ± .
22 47 . ± .
20 46 . ± .
35 48 . ± . . ± .
17 42 . ± .
12 47 . ± .
17 53 . ± . . ± .
61 41 . ± .
74 37 . ± .
44 45 . ± . . ± .
43 7 . ± .
41 11 . ± .
32 8 . ± . . ± .
82 12 . ± .
89 14 . ± .
30 15 . ± . . ± .
14 11 . ± .
23 16 . ± .
51 12 . ± .
2D FB refers to single-slice free-breathing acquisition. The frame at end exhalation was selected for analysis.Quantitative analysis was done for both single- R ∗ (WFR2S) and independent- R ∗ (WF2R2S) model-basedreconstructions.
3D BH refers to stack-of-stars breath-hold acquisition. The 17th partition was selected for analysis, given itsspatial similarity to the 2D FB slice.
3D FB refers to stack-of-stars free-breathing acquisition. The 5th bin of the 17th partition was selected foranalysis.
Even with a high acceleration factor per echo of about 10, the proposed regu-larized model-based reconstruction technique was capable of jointly estimatingseparated water/fat images, R ∗ maps, B field inhomogenity maps, as well asa set of coil sensitivity maps with good quality, as shown in Figure 5. Thereconstruction results based on the single- R ∗ model do not differ much from theindependent- R ∗ signal model (refer to Table 1).Figure 6 shows reconstruction results with only spatial and with joint spa-tial and temporal sparsity regularization. When using the latter regularizationmethod quantitative parameter maps (especially the R ∗ map) have reducednoise and streaking artifacts compared to the reconstruction with only spatialregularization. The data shown for the volumetric acquisition was obtained from the same vol-unteer as in the single-slice example. Figure 7 shows the reconstruction resultsof three selected partitions from the 16 s breath-hold scan. Model-based recon-struction yields good water/fat images, R ∗ , and field inhomogengeity maps. Fur-13 S pa t i a l + T e m po r a l S pa t i a l S pa r s i t y Water Fat R2* Fat FractionB0
Figure 6: Comparison between regularization terms using the same data asshown in Figure 5. (Top) Model-based reconstruction with only spatial sparsityregularization (joint (cid:96) P a r t i t i on P a r t i t i on P a r t i t i on R2*Water Fat B0 Fat Fraction
Figure 7: Model-based reconstruction on stack-of-stars radial multi-echo liveracquisition with 16 s breath hold. The (top) 23rd, (middle) 17th, and (bottom)12th partitions are displayed 14 a t e r R * Bin 1 Bin 2 Bin 3 Bin 4 F a t f r a c . Bin 5 Bin 6 Bin 7 (C)
Reformat View P a r t i t i on P a r t i t i on P a r t i t i on Water R2* Fat Fraction (A)
Bin 5 (B)
Bin 7 P a r t i t i on P a r t i t i on P a r t i t i on Water R2* Fat Fraction
Figure 8: Model-based reconstruction of free-breathing stack-of-stars multi-echoradial acquisition. SSA-FARY was employed to extract seven respiratory phases.(A) and (B) display the reconstructed water, R ∗ and fat fraction maps fromthree selected partitions at the 5th bin (end exhalation) and the 7th bin (endinhalation), respectively. (C) Reformatted view of all respiration phases withdotted green lines indicating the respiratory motion15her, Figure 8 shows motion-resolved model-based reconstruction results withthe application of SSA-FARY as well as joint spatial and temporal regulariza-tion. In both cases, every respiratory bin contains about 45 spokes per echo,which enables (1) 16 s breath-hold scan and (2) less than 2 min free-breathingscan for whole liver coverage.Quantitative analysis of the model-based reconstruction on single-slice free-breathing as well as stack-of-stars breath-hold and free-breathing acquisition isprovided in Table 1. Due to the variance of breathing motion, we found it diffi-cult to identify the exact same slice to be analyzed, but overall the quantitative R ∗ and fat fraction values agree in different scan protocols. For the 2D scans,the 3rd slice of scan 1 and the 1st slice of scan 2 were chosen for quantitativecomparison because their location agreed best. Quantitative analysis in Table 1shows close similarity among the first and second scan for each volunteer. This work presented dynamic liver water/fat separated R ∗ and B field map-ping based on efficient multi-echo radial sampling and regularized nonlinearmodel-based reconstruction, building on our previous work combining non-linearmodel-based reconstruction with radial sampling [23, 53–55].The multi-echo radial sampling [30] was also integrated with stack-of-starsvolumetric acquisition, offering the possibility of whole liver water/fat sepa-ration and R ∗ mapping. Accelerated acquisition with reduced scan time wasachieved via regularized nonlinear model-based reconstruction, which requiredabout 70 s per frame on the GPU and 17 min on the CPUs.Both single- R ∗ and independent- R ∗ models were implemented in this work.Results in Figure 5 show no big difference between these two models. However,note that the latter model assures a non-biased modeling and is applicable topatients with relatively high fat fraction.In general, motion is a challenging obstacle in quantitative MRI. One solu-tion to this challenge is real-time imaging [51], which remains limited to two-dimensional imaging. Here, this work achieved dynamic free-breathing liverwater/fat separation and R ∗ mapping with 1 . × . × voxel size and326 ms per frame comprising 33 excitations and 7 echoes per excitation.For 3D volumetric acquisition, motion was then either avoided via breathholding or resolved with the SSA-FARY technique. The proposed regularizednonlinear model-based reconstruction enabled an accelerated acquisition of a 3Dvolume within a single breath hold. Further acceleration could be achievablewith non-aligned acquisition, which encodes complementary spatial informationalong the partition dimension [56].This work implemented the proposed regularized model-based reconstructionin BART [31]. Under the general nonlinear inversion reconstruction framework[29], the nonlinear physics-based signal model was linearized in every Newtonstep. This linear system is then minimized as a least square problem via ADMM,which allows the flexible use of generalized and multiple regularization terms.16lthough the nonlinear model-based reconstruction is computationally demand-ing, it allows modelling of the complex multi-gradient-echo signal Equations (3)to (5) using a minimal number of unknowns in the forward model.Linear subspace methods similar to the T2 shuffling technique [57] are moreefficient, but are not suitable for R ∗ mapping for two reasons. First, a shorterecho train length (i.e. 7 in this work) is used in R ∗ mapping. Second, B fieldinhomogeneity induced phase modulation requires more principal componentsto represent the complete signal characteristics (results not shown), renderingthis technique infeasible. Therefore, it is advantageous to directly perform thenonlinear model-based reconstruction, which allows for the joint regularizationand estimation of parameter maps.While Berglund et al. [58] proposed to regularize the second-order deriva-tive of the B field map in image space, this work employed the Sobolev-normregularization [29] on the B field map. Beside rendering spatially smooth fieldmaps, more importantly, the proposed model-based reconstruction used joint (cid:96) -Wavelet regularization on W , F and R ∗ maps to achieve high-resolution R ∗ mapping. A dynamic liver water/fat separation as well as R ∗ and B field mapping tech-nique was presented. The use of multi-echo radial or stack-of-stars samplingoffered rapid and efficient k -space coverage. Generalized nonlinear model-basedreconstruction with joint-sparsity constraints on parameter maps and smooth-ness constraints on the B field and coil sensitivity maps achieved scan timereduction and accurate reconstruction of R ∗ maps. This reconstruction methodformulated the integrated parallel imaging and parameter mapping as a general-ized regularized nonlinear inverse problem, and thus required neither calibrationscans nor pre-computation of coil sensitivity maps.17 unding Statement This work was supported by the DZHK (German Centre for CardiovascularResearch), by the Deutsche Forschungsgemeinschaft (DFG, German ResearchFoundation) under grant TA 1473/2-1 / UE 189/4-1 and under Germany’s Ex-cellence Strategy—EXC 2067/1-390729940, and funded in part by NIH undergrant U24EB029240.
Open Research
In the spirit of reproducible and open research the proposed regularized model-based nonlinear inverse reconstruction is made openly available as part of theBART toolbox [31] in https://github.com/mrirecon/bart/. Scripts to producethe experiments are available at https://github.com/mrirecon/multi-echo-liverwith data published as DOI:10.5281/zenodo.4359744.
Appendix A
Regularized Model-based Nonlinear Inverse Reconstruction
To solve for the unknowns in Equation (6), IRGNM linearizes the nonlinear for-ward model via Taylor expansion in each Newton step, thus the data-consistencyterm in Equation (7) becomes (cid:107) DF ( x n )d x − [ y − F ( x n )] (cid:107) (A.1)where the Jacobian DF ( x n ) denotes the derivative of the forward operator withregard to the k th-step estimate. Given that IRGNM starts with initial guess x [29], one can denote x = x n +1 − x = x n + d x − x . As a result, Equation (A.1)becomes (cid:107) DF ( x n )( x + x − x n ) − [ y − F ( x n )] (cid:107) ⇒ (cid:107) DF ( x n ) x − [ DF ( x n )( x n − x ) + y − F ( x n )] (cid:107) (A.2)whose minimum occurs when its derivative is set to 0, DF H ( x n ) (cid:8) DF ( x n ) x − [ DF ( x n )( x n − x ) + y − F ( x n )] (cid:9) = 0 ⇒ DF H ( x n ) DF ( x n ) x = DF H ( x n ) (cid:8) DF ( x n )( x n − x ) + y − F ( x n ) (cid:9) (A.3)which represents the normal equation of the linear least square problem. Denote A := DF H ( x n ) DF ( x n ) and b := DF H ( x n ) (cid:8) DF ( x n )( x n − x ) + y − F ( x n ) (cid:9) ,and given generalized (cid:96) (cid:107) Ax − b (cid:107) + λ (cid:107) z (cid:107) subject to T x − z = 0 (A.4)18he updates can be derived, x ( k +1) := ( A H A + 0 . ρT H T )[ A H b + 0 . ρT H ( z ( k ) − µ ( k ) )] z ( k +1) := T λ/ρ ( T x ( k +1) + µ ( k ) ) u ( k +1) := u ( k ) + T x ( k +1) − z ( k +1) (A.5)The x update is solved by the conjugate gradient method, and the z updateis computed via soft thresholding ( T λ/ρ ), where λ is passed from IRGNM anditeratively reducing along Newton steps, λ = 1 /D n − with D > n the n the Newton iteration. ρ is known as the penalty parameter in ADMM.With the alternating update scheme in Equation (A.5), the following trans-form ( T ) and proximal operators were used. First, the (cid:96) W , R ∗ , F, and R ∗ was realized by the Wavelet soft-thresholding prox-imal operator, i.e. perform Wavelet transform of the input maps, apply thesoft-thresholding prox on the transformed coefficients, and then inverse Wavelettransform of the thresholded coefficients to obtain the maps. Second, the non-negativity constraint on R ∗ maps was achieved via the projection prox. Third,the (cid:96) T was identity for these regularizations. Fourth,to apply temporal TV regularization onto parameter maps, first-order deriva-tives were computed along the temporal dimension, acting as the transformoperator. A soft-thresholding operator was then applied onto the transformeddata. Appendix B
The iterative solution to Equation (A.3) requires the computation of the Jabo-cian DF ( x ) and its corresponding adjoint DF H ( x ) operator, with the forwardoperator F ( x ) denoted in Equation (6). Note that the forward operator canbe split into two nonlinear operators: the parallel imaging operator ( P F S ) andmulti-echo signal model operator ( B ). Since the first one has already been im-plemented in BART for NLINV [29], only the second operator is required to beimplemented. Afterwards, the two nonlinear operators can be chained together.Therefore, only the operator B is explained in details here.As denoted in Equation (3), the nonlinear operator B presents the mappingfrom the parameter maps (W , R ∗ , F , R ∗ , f B ) to the multi-echo images ( ρ m ),thus B : C N × (cid:55)→ C N × E . (A.6)Here, N denotes the image size, 5 represents the number of parameter mapsand E the number of echoes. Therefore, its Jacobian matrix DB ∈ C N × E × .Denote B m as the operator output corresponding to the m th TE, its corre-19ponding Jacobian is DB m = ∂B m ∂ W ∂B m ∂R ∗ ∂B m ∂ F ∂B m ∂R ∗ ∂B m ∂f B T = e − R ∗ TE m · e i πf B TE m W · (cid:16) ( − TE m ) e − R ∗ TE m (cid:17) · e i πf B TE m z m · e − R ∗ TE m · e i πf B TE m F · z m · (cid:16) ( − TE m ) e − R ∗ TE m (cid:17) · e i πf B TE m (cid:16) W · e − R ∗ TE m + F · z m · e − R ∗ TE m (cid:17) · ( i π TE m ) · e i πf B TE m T . (A.7)The adjoint operator is then its complex conjugate transpose.20 eferences [1] Caussy C, Reeder SB, Sirlin CB, Loomba R. Noninvasive, quantitativeassessment of liver fat by MRI-PDFF as an endpoint in NASH trials. Hep-atology 2018;68:763–772.[2] Hu HH, Branca RT, Hernando D, Karampinos DC, Machann J, McKen-zie CA, et al. Magnetic resonance imaging of obesity and metabolic dis-orders: Summary from the 2019 ISMRM Workshop. Magn Reson Med2020;83:1565–1576.[3] Wood JC. Impact of iron assessment by MRI. Hematology 2011;2011:443–450.[4] Hernando D, Levin YS, Sirlin CB, Reeder SB. Quantification of liver ironwith MRI: State of the art and remaining challenge. J Magn Reson Imaging2014;40:1003–1021.[5] Dixon WT. Simple proton spectroscopic imaging. Radiology 1984;153:189–194.[6] Reeder SB, Wen Z, Yu H, Pineda AR, Gold GE, Markl M, et al. MulticoilDixon chemical species separation with an iterative least-squares estimationmethod. Magn Reson Med 2004;51:35–45.[7] Reeder SB, Pineda AR, Wen Z, Shimakawa A, Yu H, Brittain JH, et al.Iterative decomposition of water and fat with echo asymmetry and least-squares estimation (IDEAL): Application with fast spin-echo imaging.Magn Reson Med 2005;54:636–644.[8] Yu H, McKenzie CA, Shimakawa A, Vu AT, Brau ACS, Beatty PJ, et al.Multiecho reconstruction for simultaneous water-fat decomposition and T ∗ estimation. J Magn Reson Imaging 2007;26:1153–1161.[9] Yu H, Shimakawa A, McKenzie CA, Brodsky E, Brittain JH, Reeder SB.Multiecho water-fat separation and simultaneous R ∗ estimation with mul-tifrequency fat spectrum modeling. Magn Reson Med 2008;60:1122–1134.[10] Doneva M, B¨ornert P, Eggers H, Mertins A, Pauly J, Lustig M. Compressedsensing for chemical shift-based water-fat separation. Magn Reson Med2010;64:1749–1759.[11] Vasanawala SS, Yu H, Shimakawa A, Jeng M, Brittain JH. Estimation ofliver T ∗ in transfusion-related iron overload in patients with weighted leastsquares T ∗ IDEAL. Magn Reson Med 2012;67:183–190.[12] Hu HH, B¨ornert P, Hernando D, Kellman P, Ma J, Reeder SB, et al.ISMRM workshop on fat-water separation: Insights, applications andprogress in MRI. Magn Reson Med 2012;68:378–388.2113] Sharma SD, Hu HH, Nayak KS. Accelerated T ∗ -compensated fat frac-tion quantification using a joint parallel imaging and compressed sensingframework. J Magn Reson Imaging 2013;38:1267–1275.[14] Glover G, Pauly JM. Projection reconstruction techniques for reduction ofmotion effects in MRI. Magn Reson Med 1992;28:275–289.[15] Rasche V, de Boer RW, Holz D, Proksa R. Continuous radial data acqui-sition for dynamic MRI. Magn Reson Med 1995;34:754–761.[16] Peters D, Korosec FR, Grist TM, Block WF, Holden JE, Vigen KK, et al.Undersampled projection reconstruction applied to MR angiography. MagnReson Med 2000;43:91–101.[17] Zhang S, Block KT, Frahm J. Magnetic resonance imaging in real time:Advances using radial FLASH. J Magn Reson Imaging 2010;31:101–109.[18] Block KT, Chandarana H, Milla S, Bruno M, Mulholland T, FatterpekarG, et al. Towards routine clinical use of radial stack-of-stars 3D gradient-echo sequences for reducing motion sensitivity. J Korean Soc Magn ResonMed 2014;18:87–106.[19] Zhong X, Armstrong T, Nickel MD, Kannengiesser SAR, Pan L, Dale BM,et al. Effect of respiratory motion on free-breathing 3D stack-of-radial liver R ∗ relaxometry and improved quantification accuracy using self-gating.Magn Reson Med 2020;83:1964–1978.[20] Schneider M, Benkert T, Solomon E, Nickel D, Fenchel M, Kiefer B, et al.Free-breathing fat and R ∗ quantification in the liver using a stack-of-starsmulti-echo acquisition with respiratory-resolved model-based reconstruc-tion. Magn Reson Med 2020;84:2592–2605.[21] Zhong X, Hu HH, Armstrong T, Li X, Lee YH, Tsao TC, et al. Free-breathing volumetric liver R ∗ and proton density fat fraction quantificationin pediatric patients using stack-of-radial MRI With self-gating motioncompensation. J Magn Reson Imaging 2020;53:118–129.[22] Wiens CN, McCurdy CM, Willig-Onwuachi JD, McKenzie CA. R ∗ -Corrected water-fat imaging using compressed sensing and parallel imaging.Magn Reson Med 2014;71:608–616.[23] Block KT, Uecker M, Frahm J. Model-based iterative reconstruction forradial fast spin-echo MRI. IEEE Trans Med Imaging 2009;28:1759–1769.[24] Fessler JA. Model-based image reconstruction for MRI. IEEE Signal Pro-cessing Magazine 2010;27:81–89.[25] Wang X, Tan Z, Scholand N, Roeloffs V, Uecker M. Physics-based Reconstruction Methods for Magnetic Resonance Imaging. arXiv2020;2010.01403. 2226] Block KT, Uecker M, Frahm J. Undersampled radial MRI with multi-ple coils. Iterative image reconstruction using a total variation constraint.Magn Reson Med 2007;57:1186–1098.[27] Lustig M, Donoho D, Pauly JM. Sparse MRI: The application of com-pressed sensing for rapid MR imaging. Magn Reson Med 2007;58:1182–1195.[28] Hernando D, Kellman P, Haldar JP, Liang ZP. Robust water/fat separationin the presence of large field inhomogeneities using a graph cut algorithm.Magn Reson Med 2010;63:79–90.[29] Uecker M, Hohage T, Block KT, Frahm J. Image reconstruction by regu-larized nonlinear inversion – Joint estimation of coil sensitivities and imagecontent. Magn Reson Med 2008;60:674–682.[30] Tan Z, Voit D, Kollmeier JM, Uecker M, Frahm J. Dynamic water/fatseparation and B inhomogeneity mapping – Joint estimation using un-dersampled triple-echo multi-spoke radial FLASH. Magn Reson Med2019;82:1000–1011.[31] Uecker M, Ong F, Tamir JI, Bahri D, Virtue P, Cheng JY, et al. BerkeleyAdvanced Reconstruction Toolbox. In: Proceedings of the 23th AnnualMeeting of ISMRM, Toronto, CAN; 2015. p. 2486.[32] Rosenzweig S, Scholand N, Holme HCM, Uecker M. Cardiac and Respi-ratory Self-Gating in Radial MRI using an Adapted Singular SpectrumAnalysis (SSA-FARY). IEEE Trans Med Imaging 2020;39:3029–3041.[33] Silva AC, Barbier EL, Lowe IJ, Koretsky AP. Radial Echo-Planar Imaging.J Magn Reson 1998;135:242–247.[34] Roemer PB, Edelstein WA, Hayes CE, Souza SP, Mueller OM. The NMRphased array. Magn Reson Med 1990;16:192–225.[35] Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: Sensi-tivity encoding for fast MRI. Magn Reson Med 1999;42:952–962.[36] Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, WangJ, et al. Generalized autocalibrating partially parallel acquisitions(GRAPPA). Magn Reson Med 2002;47:1202–1210.[37] Chebrolu VV, Hines CDG, Yu H, Pineda AR, Shimakawa A, McKenzieCA, et al. Independent estimation of T ∗ for water and fat for improvedaccuracy of fat quantification. Magn Reson Med 2010;63:849–857.[38] Hamilton G, Yokoo T, Bydder M, Cruite I, Schroeder M, Sirlin CB, et al.In vivo characterization of the liver fat 1H MR spectrum. NMR Biomed2011;24:784–790. 2339] Reeder SB, Bice EK, Yu H, Hernando D, Pineda A. On the performanceof T ∗ correction methods for quantification of hepatic fat content. MagnReson Med 2012;67:389–404.[40] Yokoo T, Bydder M, Hamilton G, Middleton MS, Gamst AC, Wolfson T,et al. Nonalcoholic fatty liver disease: Diagnostic and fat-grading accuracyof low-flip-angle multiecho gradient-recalled-echo MR imaging at 1 . R ∗ relaxometry: Theory, optimization, and clinical validation. Magn ResonMed 2013;70:1319–1331.[42] Pruessmann KP, Weiger M, B¨ornert P, Boesiger P. Adcances in sen-sitivity encoding with arbitrary k-space trajectories. Magn Reson Med2001;46:638–651.[43] Feng L, Grimm R, Block KT, Chandarana H, Kim S, Xu J, et al. Golden-angle radial sparse parallel MRI: Combination of compressed sensing, paral-lel imaging, and golden-angle radial sampling for fast and flexible dynamicvolumetric MRI. Magn Reson Med 2014;72:707–717.[44] Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimizationand statistical learning via the alternating direction method of multipliers.Foundations and Trends in Machine Learning 2010;3:1–122.[45] Guerquin-Kern M, Lejeune L, Pruessmann K, Unser M. Realistic analyt-ical phantoms for parallel magnetic resonance imaging. IEEE Trans MedImaging 2012;31:626–636.[46] Uecker M, Lai P, Murphy MJ, Virtue P, Elad M, Pauly JM, et al. ESPIRiT– an eigenvalue approach to autocalibrating parallel MRI: Where SENSEmeets GRAPPA. Magn Reson Med 2014;71:990–1001.[47] Wang X, Kohler F, Unterberg-Buchwald C, Lotz J, Frahm J, Uecker M.Model-based myocardial T mapping with sparsity constraints using single-shot inversion-recovery radial FLASH cardiovascular magnetic resonance.J Cardiovasc Magn Reson 2019;21:1–11.[48] Wang X, Rosenzweig S, Scholand N, Holme HCM, Uecker M. Model-basedreconstruction for simultaneous multi-slice T mapping using single-shotinversion-recovery radial FLASH. Magn Reson Med 2020;85:1258–1271.[49] Huang F, Vijayakumar S, Li Y, Hertel S, Duensing GR. A software channelcompression technique for faster reconstruction with many channels. MagnReson Imaging 2008;26:133–141.[50] Rosenzweig S, Holme HCM, Uecker M. Simple auto-calibrated gradientdelay estimation from few spokes using Radial Intersections (RING). MagnReson Med 2019;81:1898–1906. 2451] Uecker M, Zhang S, Voit D, Karaus A, Merboldt KD, Frahm J. Real-timeMRI at a resolution of 20 ms. NMR Biomed 2010;23:986–994.[52] Berglund J, Johansson L, Ahlstr¨om H, Kullberg J. Three-point Dixonmethod enables whole-body water and fat imaging of obese subjects. MagnReson Med 2010;63:1659–1668.[53] Roeloffs V, Wang X, Sumpf TJ, Untenberger M, Voit D, Frahm J. Model-based reconstruction for T mapping using single-shot inversion-recoveryradial FLASH. Int J Imaging Syst Technol 2016;26:254–263.[54] Tan Z, Roeloffs V, Voit D, Joseph AA, Untenberger M, Merboldt KD,et al. Model-based reconstruction for real-time phase-contrast flow MRI:Improved spatiotemporal accuracy. Magn Reson Med 2017;77:1082–1093.[55] Wang X, Roeloffs V, Klosowski J, Tan Z, Voit D, Uecker M, et al. Model-based T mapping with sparsity constraints using single-shot inversion-recovery radial FLASH. Magn Reson Med 2018;79:730–740.[56] Breuer FA, Blaimer M, Heidemann RM, Mueller MF, Griswold MA, JakobPM. Controlled Aliasing in Parallel Imaging Results in Higher Acceleration(CAIPIRINHA) for Multi-Slice Imaging. Magn Reson Med 2005;53:684–691.[57] Tamir JI, Uecker M, Chen W, Lai P, Alley MT, Vasanawala SS, et al. T shuffling: Sharp, multicontrast, volumetric fast spin-echo imaging. MagnReson Med 2017;77:180–195.[58] Berglund J, Kullberg J. Three-dimensional water/fat separation and T ∗ estimation based on whole-image optimization – Application in breathholdliver imaging at 1.5 T. Magn Reson Med 2012;67:1684–1693.25 upporting Information
7% 15% at least 30% 92 g per 100 mL7% 15% at least 30% 92 g per 100 mL