Whole brain multiparametric mapping in two minutes using a dual-flip-angle stack-of-stars blipped multi-gradient-echo acquisition

.


Introduction
Multiparametric MRI of the brain can be used to improve the assessment of neurological diseases (Seiler et al., 2021).For example, myelin water fraction (MWF) mapping can be used to assess myelination in early development of the brain (Deoni et al., 2012) and demyelination in diseases, such as multiple sclerosis (Du et al., 2007;MacKay et al., 1994).T 1 mapping can be used to study Parkinson's disease (Baudrexel et al., 2010).Proton density (PD) mapping has been used to measure tissue water content in brain tumors (Neeb et al., 2006).R 2 * mapping can be used to study iron content in deep gray matter (Haacke et al., 2010).Quantitative susceptibility mapping (QSM) provides useful information for studying iron distribution, calcification, and metabolic oxygen consumption (Wang and Liu, 2015).Multiparametric mapping provides more information about brain tissue from different aspects in comparison with single parameter mapping and has applications in macrostructural quantitative MRI, microstructural tissue modeling, and synthetic MRI (Jara et al., 2022).However, the need for a long scan time hinders its clinical applications.
Several techniques have recently been developed for simultaneous multiparametric mapping with a clinical acceptable scan time.STAGE (Chen et al., 2018) uses a dual-flip-angle (DFA) dual-gradient-echo acquisition for the mapping of T 1 , PD, R 2 *, and QSM with a scan time of 5 min and a resolution of 0.67 mm × 1.33 mm × 2.0 mm.MR Multitasking (Cao et al., 2022) uses a hybrid T 2 -prepared inversion recovery pulse and multi-gradient-echo (mGRE) readouts for the mapping of T 1 , T 2 , T 2 *, and QSM with a scan time of 9.1 min and a resolution of 0.7 mm × 1.4 mm × 2.0 mm.VFA-EPTI (Dong et al., 2021) uses a variable-flip-angle (VFA) gradient-echo sequence and echo-planar readouts for the mapping of MWF, PD, T 1 , and T 2 * with a scan time of 5 min and a 1.5 mm isotropic resolution.The above Cartesian-based techniques demonstrate the evolving landscape of fast multiparametric mapping, and the field continues to explore other innovative approaches to shorten the scan time.Stack-of-stars (SOS) sampling represents one such direction.
SOS sampling has incoherent undersampling and motion robustness properties (Feng, 2022) in comparison with Cartesian sampling.The SOS trajectory is a three-dimensional (3D) hybrid trajectory with Radial sampling in the k x -k y plane and Cartesian sampling along the k z direction (Feng, 2022).K-space data can be undersampled along k x , k y and k z directions with an SOS trajectory.On the other hand, a Cartesian trajectory only allows undersampling along the phase-encoding direction.The SOS acquisition is, therefore, expected to achieve a higher acceleration factor (R) at comparable image quality than the Cartesian acquisition, as demonstrated in a water-fat MRI study (Li et al., 2022).Besides, subspace-based reconstruction (He et al., 2016;Tamir et al., 2017) with a locally low-rank constraint (Zhang et al., 2015b) can achieve a high acceleration factor with less partial volume effect (Dong et al., 2021) by exploring the global and local similarities of multi-contrast images.While the scan time can be greatly reduced using fast acquisition sequences and advanced reconstruction algorithms, corrections should be made to ensure the accuracy of the measurement.
MWF has important applications for understanding brain development and neurodegenerative diseases (MacKay and Laule, 2016).However, it would be overestimated if the T 1 saturation effect is ignored in 3D mGRE acquisition with a short repetition time (TR), because intracellular/extracellular water components with longer T 1 values are more saturated than myelin water components (Chan and Marques, 2020).On the other hand, inhomogeneities of B 1 transmit field (B 1 + ) will induce stimulated echoes in the multi-echo spin echo sequence, which should be corrected to obtain accurate MWF (Kumar et al., 2018;Prasloski et al., 2012).As for the quantification of MWF using the mGRE sequence, B 1 + and B 1 receive field (B 1 − ) can also introduce errors.
Therefore, the corrections for the T 1 saturation effect and B 1 homogeneities are necessary for accurate quantification of MWF.Besides, the effect induced by B 1 + /B 1 − inhomogeneities need to be corrected to obtain accurate T 1 (Jang et al., 2023;Preibisch and Deichmann, 2009) and PD maps (Volz et al., 2012b).
In this study, we developed a new imaging technique for fast simultaneous 3D multiparametric mapping of whole brain MWF, T 1 , PD, R 2 *, magnetic susceptibility, and B 1

+
. First, a DFA blipped mGRE (DFA-mGRE) sequence with an SOS trajectory was developed to accelerate the acquisition.Secondly, a novel joint-sparsity-constrained multicomponent T 2 *-T 1 spectrum estimation (JMSE) algorithm was proposed to improve the quantification of MWF by correcting for the T 1 saturation effect and B 1 + /B 1 − inhomogeneities.Finally, 3D whole brain multiparametric maps were obtained with a two-minute scan by integrating the DFA-mGRE SOS sequence, the JMSE algorithm, the tissue-priorbased B 1 + estimation algorithm, and the subspace-based reconstruction.

Multi-exponential decay model
MWF has commonly been estimated based on the analysis of T 2 or T 2 * (Chen et al., 2021;Du et al., 2007;Song et al., 2020) spectrum of different brain water components.To obtain the T 2 * spectrum, mGRE datasets are fitted to a multi-exponential decay model using the nonnegative least square algorithm (MacKay et al., 1994;Whittall and MacKay, 1989).MWF can then be obtained by calculating the ratio of short T 2 * components to the sum of all components (Nam et al., 2015a).The multi-exponential decay model can be described by the integral equation: where y(TE i ) is the collected mGRE datasets of a voxel at echo time TE i , s ( T *

2
) denotes the fraction of the water component with a known T * 2 value.After discretizing T * 2 into several different values T * 2p that corre-spond to different water components, Eq. ( 1) can be simplified into a linear model: where y is a vector representing the mGRE datasets of a voxel, A is the T 2 * bases matrix, the i-th row and p-th column element of A is exp , x is a vector representing the fractions of different water components to be solved, ε represents the measurement noise.

Multi-exponential decay model with T 1 relaxation
However, the correction for the T 1 saturation effect is necessary for accurate quantification of MWF in 3D acquisition with a short TR.In this study, a commonly used single-parameter (T 2 *) multi-exponential decay model is extended into a dual-parameter (T 2 *-T 1 ) model to account for the T 1 saturation effect.The DFA-mGRE sequence is used for data acquisition.
The extended model is given by the following equation: where y ( is the collected DFA-mGRE datasets of a voxel at echo time TE i with flip angle α j , s ) denotes the fraction of the water component with a known T * 2 value and a known T 1 value.Suppose m is the number of acquired echoes and d is the number of flip angles, DFA-mGRE datasets of a voxel can form a matrix Q of size m × d with the element Q(i, j) = y ( . First, Q is vectorized by stacking the columns to obtain y.Secondly, we discretize ) , which correspond to different water components.The T 2 *-T 1 bases matrix A of size md × n can be described by: where i ranges from 1 to m, j ranges from 1 to d, p ranges from 1 to n.
Finally, the extended model described by Eq. (3) can also be simplified into a linear model described by Eq. (2).

Joint-sparsity-constrained multicomponent T 2 *-T 1 spectrum estimation (JMSE)
In the extended T 2 *-T 1 model, the number of T 2 *-T 1 bases are much larger than the number of T 2 * bases in the T 2 * model.The inverse problem becomes more underdetermined and harder to solve if the T 2 *-T 1 spectrum is not constrained.Instead of estimating the spectrum voxel-by-voxel, we propose the JMSE algorithm and incorporate a joint sparsity constraint of the T 2 *-T 1 spectrum of all voxels.The joint sparsity constraint can reduce the degrees of freedom of the solution of an underdetermined linear system, leading to better performance in applications like image reconstruction (Wang et al., 2018;Lai et al., 2018).Specifically, JMSE uses the L 2,1 -norm regularization to simultaneously maintain group sparsity and restrict the fraction within a group.The optimization problem can be described by: min where X ∈ R n * l represents the joint sparse solution of all voxels, l is the number of voxels.‖ X‖ 2,1 is the L 2,1 -norm of X. X p ∈ R l is the p-th row vector of X, and regarded as an independent group representing a } is used for the slowly-relaxing intracellular/extracellular water components (IEW) as follows.First, T 2 * is always shorter than T 1 .Secondly, the commonly used dividing point of T 2 * between MW and IEW was 25 ms (Hwang et al., 2010).Thirdly, the mean T 1 of fast-relaxing/slowly-relaxing components in white matter is 226 ms (Du et al., 2014)/1.06s (Labadie et al., 2014).The index sets of MW, IEW, and all components (ALL) are given by: } , The optimization process of JMSE was demonstrated in the Appendix.

Data acquisition
The DFA-mGRE sequence with an SOS trajectory is shown in Fig. 1 (a).The spokes of the same echo are uniformly distributed in the k x -k y plane.The incremental angle between adjacent echoes is 2 • by inserting a small blip gradient between readout gradients to obtain incoherent sampling (Li et al., 2022).This sequence is repeated with a different flip angle.A small rotation angle of 1.25 • is inserted between these two trajectories of the acquisitions with different flip angles to further reduce the sampling coherence.Gradient spoiling and RF spoiling with W. Feng et al. quadratic phase increment (Zur et al., 1991) were used for the disruption of transverse coherences in the DFA-mGRE sequence.
Datasets were acquired on a 3T scanner (uMR790, United Imaging Healthcare, Ltd., Shanghai, China) with a 24-channel head coil using the DFA-mGRE sequence with Cartesian, fully-sampled SOS (full-SOS), and prospectively under-sampled SOS (pr-SOS) sampling.Cartesian measurements with a scan time of 15.4 min were treated as the reference.
A total of 377 spokes were acquired for full-SOS with a scan time of 24.3 min.For retrospectively under-sampled SOS (re-SOS) with R = 6/ 12/18, 63/32/21 spokes were uniformly selected from full-SOS.A total of 31 spokes were acquired for pr-SOS with a scan time of 2 min (R = 12).The pr-SOS was repeated twice (pr1-SOS, pr2-SOS) for the testretest analysis.
The protocol of the in vivo study (N = 9, age 25.2 ± 1.7) was approved by the Institutional Human Ethics Committee.Written consents were obtained before each scan.The FOV was 240 mm × 240 mm × 132 mm with the voxel size of 1 mm × 1 mm × 3 mm.The flip angles were 5 • and 26 • .The TR was 43.8 ms.Twenty echoes were acquired with the echo time (TE) of 2 ms for the first echo.The echo space was 2 ms.
The phantom (System Phantom Model 130, High Precision Devices, Inc., Boulder, Colorado) with three plates of T 1 , PD, and T 2 , was scanned in the phantom study.Each plate contains 14 deionized-water filled spheres.The voxel size was 1 mm × 1 mm × 2 mm, and the TR was 43.9 ms.Other acquisition parameters were the same as the in vivo study.

Image reconstruction
Gradient delay error of SOS was corrected using k-space center shifts (Peters et al., 2003) estimated from the data acquired in four TRs at the beginning of the scan.Respiration, heart beats, and bulk motion could introduce inter-shot phase error, leading to artifacts in multi-shot acquisitions (Glover and Lai, 1998).Phase correction was performed by aligning k-space center phase with the phase of reference shot (Glover and Lai, 1998).As shown in Fig. 1(c), the k-space center phase of Stack 2 ~ N are aligned with Stack 1 by multiplying the phase difference: where K cor ijmn is the phase-corrected k-space data of the i-th stack, j-th coil, m-th echo, n-th flip angle, K ijmn is the acquired k-space data, φ i=1, jmn is the reference k-space center phase of Stack 1, φ ijmn is the k-space center phase of K ijmn .
Images acquired by Cartesian were reconstructed using the Fourier transform (FFT).Images acquired by full-SOS were reconstructed using the Nonuniform Fourier transform (NUFFT) after delay and phase correction.The subspace-based reconstruction with a locally low-rank constraint was used for re-SOS and pr-SOS, as shown in Fig. 1(c).First, phase-corrected DFA-mGRE k-space data was obtained using Eq. ( 7).Secondly, to obtain the navigator of SOS, central projections of different stacks were extracted and averaged after performing FFT on the k-space data along k z .Thirdly, the temporal bases were estimated from the Casorati matrix formed by the navigator (Feng, 2023) using principal component analysis (PCA).The maximum-normalized eigenvalues (λ) decreased rapidly, and the first twenty bases were used in the reconstruction after parameter tuning.Finally, the spatial coefficient G with respect to bases Φ, and the reconstructed multi-contrast images I were estimated by: where nuF represents the NUFFT operator, C is the coil sensitivity map estimated by ESPIRiT (Uecker et al., 2014(Uecker et al., , 2015)), K cor is the phase-corrected k-space data, λ is the regularization coefficient (0.0001 used in this work after parameter tuning), R i (G) is the transform that extracts patch i of G and reshapes it into a Casorati matrix, ‖ ⋅‖ * denotes the nuclear norm.

Quantification of parameters
Fig. 2 shows the parametric quantification pipeline.The R 2 *, QSM, B 1 + , T 1 , B 1 − , PD, and MWF maps were estimated in order.The quantification of R 2 * and QSM are described in Section 3.3.1.The B 1 + and B 1corrected T 1 are estimated as described in Section 3.3.2.The B 1 − and B 1corrected PD are estimated as described in Section 3.3.3.B 1 correction was not performed in the phantom study, because the tissue-prior-based B 1 + estimation method is not appropriate for the phantom.The T 1 -B 1corrected MWF were estimated as described in Section 3.3.4.
QSM was estimated as follows.Phase images were processed using STISuite (https://chunleiliulab.github.io/software.html).First, the brain mask was obtained using FSL (Zhang et al., 2001).Secondly, phase images were unwrapped using a Laplacian based algorithm (Schofield and Zhu, 2003).Thirdly, the background phase was removed using VSHARP (Wu et al., 2012) to obtain the tissue phase.Fourthly, dipole inversion was performed using STAR-QSM (Wei et al., 2015).Finally, the susceptibility maps of different echoes were averaged for combination.Susceptibility weighted images (SWI) (Haacke et al., 2004) were also obtained and shown after minimum intensity projection.The slab thickness of the projection was 15 mm or 5 slices.

B 1
+ correction for T 1 In spoiled gradient echo imaging with TR, TE, and nominal flip angle α, the image intensity S is given by Volz et al. (2012a): where M 0 is the equilibrium magnetization, c is a constant scaling factor without spatial variation.The VFA-based T 1 mapping requires correction for B 1 + inhomogeneity (Preibisch and Deichmann, 2009), which can be explained by: where T 1app is the map without B 1 + correction.A 5 % deviation from nominal flip angle will result in a 10 % error for T 1 (Preibisch and Deichmann, 2009).
To obtain B 1 + from DFA datasets without additional scan, Eq. ( 11) can be written as: B 1 + in WM regions can be calculated using Eq. ( 12) with the prior information of the mean T 1 value in white matter (WM) regions, then extrapolated across the whole brain (Chen et al., 2018).In this study, three subjects were scanned using the DFA-mGRE sequence and dubbed dual refocusing echo acquisition mode (DREAM) (Nehrke and Börnert, 2012) sequence to obtain B 1 + map (B 1 + dream ) and B 1 -corrected T 1 map (T 1dream ).The parameters of DREAM were: FOV = 320 mm × 320 mm × 200 mm, resolution = 5 mm × 5 mm × 10 mm, scan time = 20 s.The mean T 1dream value in WM of these three subjects (1.19 s) was used for the tissue-prior-based B 1 + estimation.Specifically, Eq. ( 12) can be written as: As shown in Fig. 2(c), the B 1 + was estimated as follows.First, The apparent T 1 map was obtained by fitting DFA data without B 1 + map.
Secondly, B 1 + map in WM was calculated using Eq. ( 13) with the WM mask obtained by Freesurfer (Reuter et al., 2010).Finally, a fifth-order 3D polynomial fitting was used to extrapolate B 1 + across the whole brain (Dong et al., 2021).With the estimated B 1 + map, the B 1 + -corrected T 1 and M 0 maps were obtained, as shown in Fig. 2(d).

B 1 correction for PD
The effect induced by B 1 − inhomogeneity need to be corrected to obtain accurate PD map.The bias field can be used as B 1 − map in the B 1 − inhomogeneity correction for PD (Volz et al., 2012a).In this study, the T 1 weighted image (TE = 2 ms, α = 26 • ), the estimated B 1 + map, and the B 1

+
-corrected T 1 map were used to estimate the bias field with FSL (Zhang et al., 2001), as shown in Fig. 2(e).
The B 1 − estimation method was described in detail as follows (Supplementary Fig. S1).First, the original T 1 weighted image A1 was corrected using B 1 + and T 1 , to remove the B 1 + effect on the magnitude image.
The corrected T 1 weighted image A2 was given by: Secondly, the bias field estimation algorithm in FSL (Zhang et al., 2001)  To cancel the effect of constant scaling factor c, PD map was scaled to achieve an average PD value of 0.7 in WM regions.

T 1 -B 1 correction for MWF B 1
+ /B 1 − inhomogeneities also introduce errors in the quantification of MWF.Therefore, B 1 correction of DFA-mGRE datasets need to be carried out first.B 1 -corrected magnitude images of DFA-mGRE datasets can be obtained using the estimated T 1 , B 1

+
, and B 1 − maps.The image intensity Ŝ with B 1 correction can be obtained by: T 1 -B 1 -corrected MWF can then be obtained with JMSE algorithm, as shown in Fig. 2(g).

Quantitative analysis
In the phantom study, mean R 2 *, T 1app or M 0app values of 14 spheres on the corresponding plate were calculated.To evaluate the consistency between pr-SOS and Cartesian, linear regression was carried out, and intraclass correlation coefficient (ICC) was calculated.The test-retest analysis was examined for pr-SOS using Bland-Altman plot and paired t-test.Coefficient of variation (CV) was also calculated.
To evaluate the effect of B 1 correction on T 1 and PD quantification, the mean, standard deviation, and histograms of WM/gray matter (GM) after correction were compared with those without correction.The commonly accepted literature values were used as references.
To evaluate the effect of T 1 -B 1 correction on MWF quantification, the mean, standard deviation, and histograms of WM/GM after correction were compared with those without correction.Besides, the geometric mean T 2 * map (gmT 2 *) (Alonso-Ortiz et al., 2018b) and mean T 1 map (mT 1 ) were also calculated: where P is the index set.Here P can be P MW , P IEW , P ALL for MW, IEW, and ALL, respectively.The commonly accepted literature values were used as references.
In the retrospective undersampling experiment, normalized mean square error (NMSE) and structure similarity index measure (SSIM) were calculated to evaluate the quality of multiparametric maps obtained using re-SOS at different acceleration factors (R = 6, 12, 18).
In the prospective undersampling experiment, the DFA-mGRE datasets obtained using Cartesian, full-SOS, re-SOS, and pr1-SOS were registered to those obtained using pr2-SOS before parametric quantification.The T 1 -weighted magnitude image (TE = 2 ms, α = 26 • ) was segmented using Freesurfer (Reuter et al., 2010).After removing CSF regions, 155 Region of Interests (ROI) consisting of 65 WM regions, cortical regions, 25 subcortical regions, 2 cerebellum regions, were obtained to calculate the mean values for each subject.To evaluate the consistency between pr-SOS and Cartesian, Bland-Altman plot with a non-parametric distribution assumption and Wilcoxon signed rank test were performed.ICC and CV were also calculated.The CV of QSM was not used in the analysis, because the QSM values of different ROIs were distributed around zero.The test-retest analysis was also examined for pr-SOS.Besides, the mean values of WM/GM for nine healthy subjects were calculated for T 1 , PD, MWF, and R 2 *.The mean values of caudate nucleus (CN), globus pallidus (GP), and putamen (PUT) for nine healthy subjects were calculated for QSM.

Evaluation of B 1 correction for T 1 and PD
Fig. 4 shows the T 1 /PD maps and histograms with/without B correction with pr-SOS sampling.The estimated B 1 + and B 1 − maps are also shown.The uncorrected T 1app map and the estimated B 1 + map have similar spatial distributions: high at the center and low at the edge.The corrected T 1 map has a more uniform spatial distribution.The standard deviations of WM/GM are reduced from 0.229 s/0.314 s to 0.136 s/ 0.243 s after correction, which indicates the improved quantification of T 1 after B 1 correction.As for PD, the corrected PD map has a more uniform spatial distribution.The mean PD values of WM/GM after correction (0.700/0.787) are in the same range with literature values (0.709/0.812) (Neeb et al., 2006), (0.697/0.810) (Volz et al., 2012b), and the standard deviations of WM/GM are reduced after correction, which indicates the improved quantification of PD after B 1 correction.

Evaluation of T 1 -B 1 correction for MWF
Fig. 5(a) shows the uncorrected maps using the T 2 * model and singleflip-angle mGRE datasets (5 • or 26 • ), the T 1 -corrected maps using the extended T 2 *-T 1 model and uncorrected DFA-mGRE datasets, the T 1 -B 1corrected maps using the extended T 2 *-T 1 model and B 1 -corrected DFA-mGRE datasets with pr-SOS sampling.The uncorrected MWF is higher overall than the T 1 -B 1 -corrected MWF.The T 1 -corrected MWF and the estimated B 1 + map have similar spatial distributions.The T 1 -corrected MWF is lower than the T 1 -B 1 -corrected MWF in regions where B 1 + value is greater than one.The uncorrected and the T 1 -B 1 -corrected gmT 2 * MW maps are similar overall, while the T 1 -corrected gmT 2 * MW map is lower in regions where  ( Whittall et al., 1997), which indicates the improved quantification of MWF after T 1 -B 1 correction.

Retrospective undersampling experiment
Fig. 6 shows the multiparametric maps and corresponding error maps using re-SOS at different acceleration factors (R = 6, 12, 18).As the acceleration factor increases, the NMSE increases while the SSIM decreases.The errors of MWF, QSM, and R 2 * are higher in regions with severe B 0 inhomogeneity as indicated by the yellow arrows.The SSIM (0.9688 for T 1 , 0.9629 for PD, 0.9992 for B 1 + , 0.9762 for MWF, 0.9890 for QSM, 0.9621 for R 2 *) and NMSE (0.0174 for T 1 , 0.0027 for PD, 0.0002 for B 1 + , 0.0681 for MWF, 0.1204 for QSM, 0.0456 for R 2 *) indicate good quality of the multiparametric maps obtained using re-SOS at R = 12.The high quality of the reconstruction is maintained even at R = 18, which means that the scan time could be further reduced to 1.35 min. of Cartesian and pr2-SOS with a non-parametric distribution assumption.The ICC (0.986 for T 1 , 0.985 for PD, 0.994 for B 1 + , 0.967 for MWF, 0.973 for QSM, 0.920 for R 2 *) and CV (2.42 % for T 1 , 0.99 % for PD, 0.98 % for B 1 + , 17.2 % for MWF, 7.18 % for R 2 *) indicate good consistency between pr2-SOS and Cartesian reference.There is a small bias (P < 0.05) of − 0.013 s for T 1 , − 0.002 for PD, − 0.002 for B 1 + , − 0.001 for MWF, and 0.803 Hz for R 2 *.There is no significant (P > 0.05) bias for QSM (P = 0.10).
The mean values of multiparametric maps for nine healthy subjects in different ROIs using Cartesian, pr1-SOS, and pr2-SOS are listed in Table 1.The multiparametric maps and SWI images of nine healthy subjects with pr-SOS sampling acquired in two minutes are shown in Fig. 9.

Stop criterion and parameter selection of the JMSE algorithm
The problem to be solved by the JMSE algorithm involves the fitting of DFA-mGRE datasets using the T 2 *-T 1 bases, which is an underdetermined linear system.Fig. 10(a) shows T 1 -corrected and T 1 -B 1corrected MWF estimated by the JMSE algorithm at different iteration numbers.Fig. 10(b) shows the logarithm of the mean absolute fitting error (MAE) and standard deviation (SD) of MWF in WM plotted as a function of iteration number.The MAE decreases with increasing the iteration number until convergence is achieved.Besides, the MAE curve with two asymptotes exhibits a l-shape, suggesting that further increasing the iteration number has little contributions to the reduction of fitting error.On the other hand, the SD of MWF in WM curve exhibits a V-shape when the iteration number is larger than 150, which indicates that the JMSE algorithm falls into overfitting, as demonstrated in a multi-compartment relaxometry informed myelin water imaging study (Chan and Marques, 2020).To avoid underfitting and overfitting, the JMSE algorithm should be stopped appropriately.In this study, the intersection of the two asymptotes was used to determine the optimal iteration number T.
The algorithm was run with different β values to figure out the effect of the penalty parameter on the estimated MWF, as shown in Fig. 11(ad).The l-curve approach (Hansen, 1992) was used to determine the optimal penalty parameter.Referring to Eq. ( 5), the data consistency term is defined as ‖ AX − Y‖ F , and the regularization term is defined as ‖ X‖ 2,1 .For β = 1000, 100, the data consistency terms are relatively large, and MWF is overestimated.For β = 1, 0.1, 0.01, 0.001, and 0.0001, the regularization terms are relatively large, and the standard deviations of MWF in WM are relatively large.Therefore, β = 10 was used in this work.
The speed of convergence for the JMSE algorithm is determined by the step length γ.The upper bound of γ for theoretical convergence is 1.618 (Deng et al., 2013).The algorithm was run with different γ values to determine an appropriate step length.The fitting error curve with a small step length (γ = 0.001) drops slowly and does not start to converge when the algorithm stops, as shown in Fig. 11(e-f).The curve with γ = 0.01 drops monotonically and then converges without oscillations.The curve with γ = 0.1 drops quickly and then converges with small oscillations.The curve with a large step length (γ = 1) drops quickly and then

Discussion
This study has demonstrated the feasibility of a technique integrating the DFA-mGRE SOS sequence, the proposed JMSE algorithm, the tissueprior-based B 1 + estimation algorithm, and the subspace-based reconstruction for simultaneous 3D whole brain multiparametric mapping with a two-minute scan.The T 1 -B 1 correction method using the JMSE algorithm can improve the quantification of MWF in 3D acquisition.The B 1 correction method using the tissue-prior-based B 1 + estimation algorithm can improve the quantification of T 1 and PD.The phantom and prospective experiments have demonstrated good consistency between results of Cartesian reference and pr-SOS, as well as good reproducibility between two repeated prospective SOS scans.W. Feng et al. et al., 2016) and simulation-based methods (Dong et al., 2021;Tamir et al., 2017) are widely used for temporal bases estimation in subspace-based reconstruction.Navigator-based bases estimation, which does not require prior information, is adopted in this study.The navigator is formed by the central k-space lines, which takes advantage of the self-navigation property of SOS and saves time in acquiring the navigator data.The retrospective experiment has demonstrated good quality of the multiparametric maps obtained using re-SOS with the subspace-based reconstruction at R = 12 or 18 (scan time = 2 or 1.35 min).
In addition to the advantages mentioned above, SOS sampling also has several limitations (Feng, 2022).First, the overall reconstruction time of SOS sampling is longer than that of Cartesian sampling.Secondly, SOS sampling is sensitive to gradient delay.K-space samples will be shifted along the readout direction in the presence of gradient delay, leading to image inhomogeneity and blurring (Feng, 2022).Gradient delay correction should be performed carefully in the acquisition and reconstruction.Thirdly, SOS sampling is sensitive to eddy current.Strong eddy current will be produced by fast gradient switching.The magnitude of eddy current is related to the rotation angle between spokes and the echo space.Scanning with large rotation angle and small echo space requires fast gradient switching, which consequently leads to the generation of image artifacts associated with eddy current (Feng, 2022).It is important to choose the rotation angle, the echo space, and the spatial resolution appropriately.Instead of golden-angle rotation scheme, linear rotation scheme that produces small rotation angle, was In the normalization step, PD map was scaled to achieve an average PD value of 0.7 in WM regions.Therefore, the standard deviation of PD in WM across subjects was zero and was not listed.
Fig. 9. Multiparametric maps and SWI images of nine healthy subjects using prospectively under-sampled SOS acquired in two minutes.
used in this study to reduce the effect of eddy current.Finally, SOS sampling is sensitive to off-resonance (Feng, 2022) and B 0 inhomogeneity, leading to image artifacts and blurring.At present, the effect of B 0 inhomogeneity is unavoidable in the proposed technique.Typically, B 0 inhomogeneity significantly affects images with long TE, while its impact on images with short TE is relatively minor.The B 0 inhomogeneity effect on short TE images will be magnified using the subspace-based reconstruction with bases estimated from SOS navigator, because the SOS navigator are formed by data from different TEs (2 ~ 40 ms).The reconstructed DFA-mGRE images, including the images with short TE, will be affected by B 0 inhomogeneity related artifacts.The T 1 and PD maps estimated by DFA data of the first echo will also be affected, and the errors are relatively large at the edge of the brain with higher B0 inhomogeneity as shown in Fig. 6.To mitigate the B 0 inhomogeneity effect, the echo space and TE can be shortened while maintaining the spatial resolution by ramp sampling.On the other hand, more effective B 0 correction methods can be investigated.

Quantification of parameters
In the tissue-prior-based B 1 + estimation algorithm, there are three key points that require careful attention.First, the mean T 1 values in WM, as the prior information, should be determined based on the experiments using an independent B 1 + mapping technique.The mean T 1 values in WM reported by other study may not be directly applied to a different MRI scanner and sequence.The estimated B 1 + and B 1 -corrected T 1 will be biased if the T 1 prior information deviate from the actual or reference value.Secondly, Accuracy of the WM mask is important, for Eq. ( 13) only approximately holds in the WM region.Finally, a fast 3D global interpolation method is preferred than a 2D local interpolation method, to avoid slice discontinuity.Therefore, a fifth-order 3D polynomial fitting (Dong et al., 2021) was used in this study.Besides, The estimated B 1 + map was compared with the B 1 + map measured using DREAM (B 1 + dream ).The B 1 -corrected T 1 map using the estimated B 1 + was also compared with that corrected using B 1 + dream (T 1dream ), as shown in Supplementary Fig. S2.B 1 + and B 1 + dream have a similar spatial distribution: high at the center and low at the edge.T 1 and T 1dream have a more uniform spatial distribution than T 1app .The mean value of T 1 in WM/GM are close to the mean value of T 1dream , which indicate that the estimated B 1 + map can be used as an alternative to the B 1 + measured using DREAM.
As for the quantification of R 2 * and MWF, the B 0 inhomogeneity will result in signal loss in regions with severe B 0 inhomogeneity in the 3D mGRE-based imaging, and R 2 */MWF will be overestimated in these regions.With a linear approximation of the background field, the additional signal loss induced by the local field gradients can be corrected using a sinc function (Hwang et al., 2010;Yablonskiy, 1998) or the voxel spread function (Yablonskiy et al., 2013).As shown in Supplementary Fig. S3, the estimation of R 2 * and MWF in regions with moderate local field gradients can be improved after correction.On the other hand, the iron content has an effect on the estimation of MWF.In regions with high iron content, the signal attenuation occurs more rapidly, and the estimated MWF does not provide a realistic measure of myelination.Whereas, the effect of iron content can be considered negligible in regions with low iron content, and the estimated MWF can provide useful information for studying myelination in these regions.Moreover, iron and myelin could be separated based on the distinction that iron is paramagnetic and myelin is diamagnetic, which can be

Fig. 1 .
Fig. 1.(a) Diagram of the DFA-mGRE sequence with an SOS trajectory.Blip gradients represented by red triangles were inserted between readout gradients.(b) Spokes between adjacent echoes were rotated by 2 • .(c) Image reconstruction pipeline of re-SOS and pr-SOS.
was run three times.Finally, B 1 − was obtained by multiplying the three intermediate bias field maps.With the estimated B 1 − and R 2 * maps, PD was obtained after transverse relaxation correction and B 1 − inhomogeneity correction, as shown in Fig. 2(f).

W
. Feng et al.B 1 + value is greater than one.For gmT 2 * IEW and gmT 2 * ALL , the uncorrected, the T 1 -corrected, and the T 1 -B 1 -corrected maps are similar overall.The mean gmT 2 * IEW values of WM/GM with T 1 -B 1 correction (59.2 ms/61.2ms) are in the same range with literature values (WM: 61 ms, deep GM:53 ms, cortical GM: 66 ms) (Alonso-Ortiz et al., 2018b).The T 1 -B 1 -corrected mT 1MW map has a more uniform spatial distribution than the T 1 -corrected mT 1MW map.For mT 1IEW and mT 1ALL , the T 1 -corrected maps are higher than the T 1 -B 1 -corrected maps in regions where B 1 + value is greater than one.The mean mT 1MW (0.25 s) value of WM with T 1 -B 1 correction is in the same range with literature values (mean T 1MW of WM: 0.226 s) (Du et al., 2014).MWF histograms of WM/GM are shown in Fig. 5(b).The mean MWF values of WM/GM after T 1 -B 1 correction (0.099/0.033) are closer to literature values (0.100/0.028) (Chen et al., 2020), (0.113/0.031)

Fig. 3 .
Fig. 3. (a) Comparison of R 2 */T 1app /M 0app among Cartesian, full-SOS, re-SOS, pr1-SOS, and pr2-SOS.(b) Linear regression of Cartesian and pr2-SOS (14 ROIs).The red line represents the fitted line.The error bar represents the standard deviation.(c) Bland-Altman plots of pr1-SOS and pr2-SOS.The solid red line is the mean difference with the P value of paired t-test.The dotted lines refer to confidence limits of the mean difference ±1.96 standard deviation (SD).

Fig. 4 .
Fig. 4. (a) Comparison of uncorrected T 1app , M 0app maps and B 1 -corrected T 1 , PD maps.The estimated B 1 + and B 1 − maps are shown on the right.(b) T 1 (left) and PD (right) histograms of WM/GM.The number of bins is 100.The mean and standard deviation are shown using the corresponding color.

Fig. 5 .
Fig. 5. (a) Comparison of uncorrected, T 1 -corrected, and T 1 -B 1 -corrected MWF, gmT 2 *, and mT 1 maps.The mean and standard deviation of gmT 2 * and mT 1 of WM/ GM are shown in white color.The logarithm of the mean absolute fitting error (MAE) are plotted as a function of iteration number.The stop points (optimal iteration number T) of the JMSE algorithm are shown using the corresponding color.(b) MWF histograms of WM/GM.The number of bins is 100.The mean and standard deviation are shown using the corresponding color.

Fig. 6 .
Fig. 6.Evaluation of multiparametric maps obtained using re-SOS with R = 6, 12, 18.The multiparametric maps and corresponding 5 × error maps are shown on the left, the NMSE and SSIM are plotted on the right.

Fig. 7 .
Fig. 7. (a) Comparison of T 1 /PD/B 1 + among Cartesian, full-SOS, re-SOS, pr1-SOS, and pr2-SOS.(b) Bland-Altman plots of Cartesian and pr2-SOS.(c) Bland-Altman plots of pr1-SOS and pr2-SOS.The red solid line is the median difference with the P value of Wilcoxon signed rank test.The dotted lines refer to confidence limits of the median difference ±1.45 interquartile range (IQR).

Fig. 8 .
Fig. 8. (a) Comparison of MWF/QSM/R 2 * among Cartesian, full-SOS, re-SOS, pr1-SOS, and pr2-SOS.(b) Bland-Altman plots of Cartesian and pr2-SOS.(c) Bland-Altman plots of pr1-SOS and pr2-SOS.The red solid line is the median difference with the P value of Wilcoxon signed rank test.The dotted lines refer to confidence limits of the median difference ±1.45 interquartile range (IQR).ppm: parts per million; ppb: parts per billion.

Fig. 10 .
Fig. 10.(a) T 1 -corrected and T 1 -B 1 -corrected MWF estimated by the JMSE algorithm at the optimal iteration number and other iteration numbers.(b) The logarithm of the mean absolute fitting error (MAE) (blue dotted line) and standard deviation (SD) of MWF in WM (orange line) plotted as a function of iteration number.The blue line represents an asymptote.The coordinates of the intersection of two asymptotes are also shown.

Fig. 11 .
Fig. 11.(a-d) Penalty parameter (β or beta) selection experiments.(a) MWF estimated by the JMSE algorithm with different penalty parameters.(b) The logarithm of the MAE curves with different penalty parameters.(c) The mean (blue) and SD (orange) of MWF in WM plotted as a function of beta.(d) The logarithm plot of regularization term vs data consistency term (the l-curve approach).(e-f) Step length (γ or gamma) selection experiment.The logarithm of the MAE curves with different step lengths (e) and its magnified view (f).