Investigating the accuracy and precision of TE-dependent versus multi-echo QSM using Laplacian-based methods at 3 T

Purpose: Multi-echo gradient-recalled echo acquisitions for QSM enable optimiz-ing the SNR for several tissue types through multi-echo (TE) combination or investigating temporal variations in the susceptibility


| INTRODUCTION
Quantitative susceptibility mapping aims to determine the underlying spatial distribution of tissue magnetic susceptibility (χ) 1,2 from the phase ( ) of a gradient-recalled echo (GRE) MRI sequence: where TE is the echo time; r is the voxel position in the image; is the proton gyromagnetic ratio; ΔB Tot is the χ-induced total field perturbation along the scanner's z-axis; and 0 is a phase offset at TE = 0.
Processing pipelines for QSM usually require unwrapping to resolve spatiotemporal 2 aliasing, removing background field variations (ΔB Bg ) resulting from χ sources outside the brain, and solving a local field (ΔB Loc )-to-χ ill-posed inverse problem. 2 Both single-echo and multi-echo GRE acquisitions can be used for QSM. Acquiring multiple TEs is advantageous primarily because the phase SNR (SNR( )) can be optimized in multiple tissues simultaneously, 1-3 as choosing TE = T * 2 for a tissue maximizes SNR( ) in that tissue. 1,4 Thus, combining multiple echoes via fitting or averaging over TEs results in a TE-independent field map designed to optimize SNR( ) (see Equation 1). By acquiring multiecho GRE images and processing each TE separately, a method referred to as TE-dependent QSM, 5 some studies [5][6][7] have investigated the TE dependence of the QSM susceptibility resulting from nonlinear temporal variations of the phase. Acquiring multi-echo GRE images and combining all TEs improves the SNR of the resulting field map 3,7 or QSM image 1,7 compared to acquiring single-echo GRE images, but the regional accuracy and precision of χ calculated using TE-dependent versus multi-echo QSM has undergone limited investigation.
Regarding the accuracy of χ in TE-dependent QSM, χ measured using numerical simulations in a vein perpendicular to the external magnetic field was shown to be constant across a long range of TEs. 1 However, χ measured using numerical simulations 6 or real data 5,6 was also shown to be TEdependent, particularly at short TEs. In high-χ regions, such as microbleeds, such a TE dependence could arise from a failure of Laplacian phase unwrapping to recover the true phase, 6 or, motivated by the observation of a regional TE dependence of χ, it could be linked to intrinsic tissue properties. 5,6 Regarding the precision of χ, simulated 1,3 and in vivo data 7 have shown that, for the range of TE values such that the measured MR signal is above the noise floor, single-echo χ noise decreases at longer TEs, combining multiple TEs by averaging over echoes reduces the propagation of phase noise into ΔB Loc compared with using single echoes 3 , and fitting an increasing number of TEs progressively reduces the propagation of phase noise into the multi-echo χ map compared with using single echoes. 7 Processing pipelines for QSM often incorporate Laplacian-based methods (LBMs) for phase unwrapping 8 or ΔB Bg removal. 9 Laplacian phase unwrapping aims to identify the unwrapped phase whose local derivatives are the most similar to the derivatives of the wrapped phase. Laplacian ΔB Bg removal aims to eliminate ΔB Bg from the brain by exploiting their harmonicity (ie, ∇ 2 ΔB Bg = 0) inside the brain. Because LBMs perform nonlinear operations, they must be applied carefully to avoid introducing inaccuracies into the processed signal phase. Notably, the effect on the accuracy and precision of χ using Laplacian-based phase unwrapping and ΔB Bg removal in the same pipeline for multi-echo or TE-dependent QSM has not been investigated. In fact, previous studies on the noise in single-echo versus multi-echo QSM have not used LBMs 1,7 or have used LBMs to perform ΔB Bg removal but not phase unwrapping. 7 Moreover, they have not compared the accuracy and precision of ΔB Loc and χ calculated at each TE with the accuracy and precision of ΔB Loc and χ calculated by combining all TEs. 1,3,[5][6][7] Because LBMs for phase unwrapping and ΔB Bg removal are frequently applied within the same QSM pipeline, 5,6 it is important to systematically evaluate their effect on the accuracy and precision of TE-dependent versus multi-echo QSM. This could contribute to the standardization of image processing pipelines for TE-dependent QSM. Therefore, this study aimed to compare the effect of LBMs on the accuracy and precision of χ in TE-dependent versus multi-echo QSM. This was investigated using numerically simulated data and images acquired in vivo, incorporating LBMs for both phase unwrapping and ΔB Bg removal into the QSM processing pipelines.

| METHODS
When not otherwise stated, image analysis was performed using MATLAB (R2017b; The MathWorks, Natick, MA) magnetic resonance imaging, multi-echo QSM, quantitative susceptibility mapping, single-echo QSM, TE-dependent QSM and statistical analysis was performed using Stata (R15; StataCorp, College Station, TX).

| In vivo data acquisition
Multi-echo 3D GRE imaging of 10 healthy volunteers (average age/age range: 26/22-30 years, 5 females) was performed on a 3 T system (Philips Achieva [Philips Medical Systems, Best, NL]; 32-channel head coil). A preliminary version of this study was presented at the 2017 meeting of the International Society for Magnetic Resonance in Medicine. 10 All of the volunteers provided written informed consent, and the local research ethics committee approved the experimental sessions. Images were acquired using a transverse orientation (FOV = 240 × 240 × 144 mm 3 , isotropic voxel size = 1 mm 3 , flip angle = 20º, TR = 29 ms, five evenly spaced echoes [(TE 1 /TE spacing = 3/5.4 ms], bandwidth = 270 Hz/pixel, SENSE 11 along the first/second phase-encoding directions = 2/1.5, flyback gradients along the frequency-encoding direction, and without flow-compensating gradients).

| Data simulation from a numerical head phantom
To ensure the availability of ground-truth χ values against which to test the accuracy of the QSM pipelines, a piece-wise constant Zubal numerical head phantom 12 was used with a matrix size of 256 × 256 × 128 voxels 3  , as QSM has often been used to study iron deposition in these regions 4 ; one venous ROI in the superior sagittal sinus (SSS), as QSM can be used to investigate vein oxygenation in the healthy brain 13,14 and in cerebrovascular disease 15,16 ; and, to represent the other main brain tissue types, cortical gray matter (GM), white matter (WM), and CSF ROIs were also included (see Supporting Information Figure S2). Based on this phantom, the ground-truth total field perturbation (ΔB True Tot ) at 3T MRI was simulated using a Fourier-based forward model. 17 A noise-free multi-echo complex MRI signal C was simulated throughout the brain using Equation 1 with ΔB True Tot and a spatially variable 0 (see Supporting Information) as well as a mono-exponential T * 2 magnitude signal decay: where i is the unit imaginary number. To enable the comparison of in vivo results and numerical simulations, the TEs in Equation 2 were matched to those of the in vivo 3D GRE protocol. Random zero-mean Gaussian noise with a standard deviation (SD) equal to 0.07 was added to the real and imaginary parts of the noise-free signal independently. 18 This value of SD was calculated by fulfilling the following high-SNR condition at the longest TE in the numerical phantom ROI with the shortest T * 2 (ie, the SSS 19 ): where A and (A), respectively, denote the magnitude image and its noise. Therefore, realistic noise affected the numerical simulations at all TEs and in all ROIs except the region outside the head, which was masked out during the QSM calculation.

| Data preprocessing
A brain mask was calculated for each subject by applying the FSL Brain Extraction Tool (BET) 20,21 with robust brain center estimation (threshold = 0.3) to the magnitude image at the longest TE. The mask was based on the longest TE, to account for the greater amount of signal loss (signal dropout near regions of high-susceptibility gradients) compared to shorter TEs. A whole-brain mask for the Zubal phantom was calculated by applying FSL BET 20,21 with robust brain center estimation (threshold = 0.05) to the T * 2 map of the numerical phantom (see Supporting Information).

| Processing pipelines for QSM
For TE-dependent QSM, the signal phase from each TE was processed separately. For multi-echo QSM, all of the processing pipelines received as an input the same total field map, which was calculated by nonlinearly fitting (NLFit) the complex multi-echo signal over TEs 22 using the Cornell QSM software package's Fit_ppm_complex function (http:// weill.corne ll.edu/mri/pages /qsm.html).
Three distinct processing pipelines incorporating LBMs for phase unwrapping and ΔB Bg removal were implemented as follows (Figure 1). The first processing pipeline applied simultaneous phase unwrapping and ΔB Bg removal using SHARP (sophisticated harmonic artifact reduction for phase data), which is a direct solver of the discretized Poisson equation. 23 SHARP was chosen because it has been widely used in the literature on QSM and has been shown to be both robust and numerically efficient. 9 To minimize the number of voxels at the edge of the brain affected by convolution artifacts, SHARP was implemented using the minimum-size 3-voxel isotropic 3D Laplacian kernel. 23,24 SHARP was applied using a threshold for singular value decomposition equal to 0.05 and a brain mask with 2-voxel (numerical phantom) or 4-voxel erosion (healthy subjects). The threshold for singular value decomposition value used here was consistent with the values used in the original studies on SHARP 23,24 and led to the best ΔB Bg removal on visual inspection. Both the second and third processing pipelines first applied Laplacian phase unwrapping with a threshold for singular value decomposition equal to 1 -10 (ie, the default value in the Cornell QSM software package). Then, because larger kernel sizes have been shown to better recover small susceptibility variations, the second processing pipeline applied SHARP with variable kernel size (V-SHARP) 25 to the Laplacian-unwrapped phase with an initial kernel radius of approximately 7 mm and a 1-voxel step size/final kernel radius. The third strategy applied the Laplacian boundary value (LBV) method 26 to the Laplacian-unwrapped phase. The LBV was tested because, in a recent multicomparison study, 9 it quantitatively outperformed all other methods for ΔB Bg removal. The Cornell QSM software package's LBV function was used with default parameter values.
For TE-dependent QSM only, residual phase offsets were removed by applying V-SHARP with an initial kernel radius of approximately 40 mm and a 1-voxel step-size/final kernel radius. 27 In all these pipelines, ΔB Loc -to-χ inversion was performed using Tikhonov regularization 2,28 with correction for susceptibility underestimation 23 and using the L-curve method 29 to determine the optimal value for the regularization parameter. This inversion method was chosen because it is computationally efficient and has been shown to substantially reduce streaking artifacts relative to the truncated k-space division method. 15 For each processing pipeline in the healthy volunteers, the mean across subjects of the individually optimized Tikhonov regularization parameters was used.
Finally, to investigate a single-step processing pipeline incorporating LBMs, a fourth processing pipeline using total generalized variation (TGV) 30 was applied to the wrapped phase at each TE (TE-dependent QSM) or the NLFit field map (multi-echo QSM). The TGV method was chosen because it has been shown to avoid stair-casing artifacts in the resulting QSM image while correctly preserving structural borders. 30 It was implemented using Singularity (https:// github.com/CAIsr /qsm; Sylabs, Albany, CA) and the default Tot ) underwent exactly the same processing steps, following each processing stream described in the gray box parameter values 1 , 0 = (0.0015, 0.005) adopted in the original TGV study based on the criterion that the 1 : 0 = 3:1 ratio is optimal for medical imaging applications. 30

| Region-of-interest segmentation in the healthy subject images
To enable the comparison of regional χ values between the simulated and in vivo data, ROIs similar to those defined in the numerical phantom were segmented in vivo (ie, the CN, GP, PU, TH, posterior corona radiata [PCR] as a WM ROI, and the straight sinus [StrS] as a venous ROI). We chose to consider the PCR instead of the whole WM, to minimize WM microstructural effects in the regional χ analysis, as WM microstructure was not modeled in the numerical phantom simulations. We chose to consider the StrS instead of the SSS as the venous ROI, because QSM calculation in vivo required larger amounts of brain mask erosion, which left only a few SSS voxels available for regional χ analysis.
For each subject, the CN, GP, PU, TH, and PCR were segmented based on the Eve χ atlas 31 ( Figure 2), whose GRE magnitude image was aligned to each subject's fifth-echo magnitude image using NiftyReg. 32,33 The fifth-echo magnitude image was chosen because it had the closest TE (24.6 ms) to that of the Eve GRE magnitude image (24 ms). Moreover, because the fifth-echo magnitude image had the best contrast between the StrS and the surrounding brain parenchyma, it was used to segment the StrS of each volunteer using the ITK-SNAP active contour segmentation tool 34 ( Figure 2).

| Quantitative evaluation of the measured χ
All of the processing pipelines ( Figure 1) were applied to the noisy simulated images of the numerical phantom and the images of the 10 healthy volunteers. To avoid combining effects induced by TE dependence or multi-echo combination in a reference ROI with similar effects in the ROIs under investigation, the χ maps were not referenced (for example, to the CSF). Instead, it was assumed that ΔB Bg removal provided sufficient χ referencing for the purpose of this study. 9 In the numerical phantom simulations, the performance of each QSM pipeline relative to the ground truth was visually assessed by calculating a difference image between the corresponding χ map and True . In the volunteers, because of the lack of a ground truth, TE-dependent difference images were calculated relative to the corresponding multi-echo QSM image.
In the numerical phantom simulations, the means and SDs of χ were calculated for each pipeline in each phantom ROI with ≠ 0 (ie, the CN, GP, TH, PU, SSS, and WM) (Supporting Information Figure S2). The root mean squared errors (RMSEs) of χ relative to True were also calculated throughout the brain volume as where ̂ denotes the χ map calculated by the QSM pipeline, and ‖⋅‖ 2 is the Euclidean norm.
In the volunteers, the means and SDs of χ were calculated in each ROI and for each processing pipeline and compared against χ values in subjects of a similar age from the literature. In the volunteers, RMSEs could not be calculated because of the lack of a ground truth.
For visualization purposes in the volunteers, if pooling the regional means and SDs of χ across subjects was possible (ie, all intrasubject SDs of χ were larger than the intersubject SD of χ), the pooled χ averages (m p ) and SDs (SD p ) were respectively calculated as and where m i and SD i respectively denote the individual mean and SD of χ in the ith subject, and n i denotes the ROI size in voxels.

| Statistical analysis
Statistical analyses were only performed on the volunteer data, because they required several distinct measurements for the average χ in each ROI. The presence of a systematic bias in TE-dependent versus multi-echo QSM within the same processing pipeline (Figure 1) was investigated using Bland-Altman analysis of the average χ at each TE and for each ROI. Statistically significant differences in multi-echo versus TE-dependent QSM within the same processing pipeline and ROI and at each TE were tested by considering the corresponding distributions of average χ values across subjects. To assess whether to apply parametric paired t-tests or nonparametric sign tests, the normal distribution of the differences between paired χ values was assessed using the Shapiro-Wilk test.

| Pooling of χ measurements
For each ROI and each processing pipeline, all intrasubject SDs of χ were larger than the intersubject SD; thus, m p and SD p were calculated according to Equations 5 and 6.

| Performance of TE-dependent versus multi-echo QSM
Here, only the results of the SHARP-based and TGV-based processing pipelines are presented, because these pipelines best estimated True in the numerical phantom simulations. The results of the V-SHARP-based and LBV-based processing pipelines are available in the Supporting Information.
In the numerical phantom, Figures 3 and 4 show the multiecho or TE-dependent susceptibility images (B-G and I-N) For the numerical phantom simulations, the error between the calculated and ground-truth χ appeared fairly constant over TEs (Figures 3 and 4) except in some high-χ regions in which the error increased with TE (ie, the GP and SSS in the SHARP-based pipeline [ Figure 5A] and the SSS only in the TGV-based pipeline [ Figure 5B]).
The SSS, which was the highest χ structure ( True = 0.3 ppm), showed the largest susceptibility errors for both SHARP-based and TGV-based TE-dependent QSM ( Figure  5A,B). In the SHARP-based pipeline, large susceptibility errors were also observed in the GP, which was the second-highest χ structure ( True = 0.19 ppm).
For one representative volunteer, Figures 6A-L and 7A-L show the multi-echo and TE-dependent susceptibility images respectively calculated using the SHARP-based and TGV-based pipelines. In addition, Figures 6M-V and 7M-V respectively show TE-wise differences between TE-dependent and multi-echo QSM. For each ROI in the healthy volunteers, Figure 5 shows the pooled means and SDs of multi-echo and TE-dependent χ calculated using the SHARP-based (C) or TGV-based pipelines (D). In the volunteers, susceptibility differences between TE-dependent and multi-echo QSM were most prominent in the StrS (arrowheads in Figures 6G-L,R-V and 7G-L,R-V) and at the shortest TEs.
In general, the average susceptibilities measured in the deep-GM ROIs and in the PCR ( Figure 5C,D) had values within the ranges reported by previous studies 25,31,[35][36][37][38] : 0.01-0.13 ppm for the CN, 0.06-0.29 ppm for the GP, 0.02-0.14 ppm for the PU, −0.02-0.08 ppm for the TH, and −0.06-0.03 ppm for the PCR. An exception was observed in the CN, in which the mean of TE 1 was always smaller than the corresponding minimum value reported in the literature (ie, 0.01 ppm). In the StrS, only the mean of the multi-echo χ calculated using the SHARP-based pipeline had a value close to the previously reported range for venous blood, namely, 0.17-0.58 ppm 13,15,37,39 (Figure 5C,D).
In both the numerical phantom simulations and the healthy volunteers, the QSM images appeared less noisy with

| Statistical analysis
In the healthy volunteers, for all processing pipelines and ROIs, the Shapiro-Wilk test rejected the hypothesis of normally distributed paired differences of χ. Therefore, TEdependent versus multi-echo QSM values were always evaluated using the nonparametric sign test.
For the SHARP-based and TGV-based processing pipelines, significant differences between TE-dependent and multi-echo QSM are highlighted in Figure 5C,D, whereas the bias between TE-dependent and multi-echo χ is shown in Figure 8. The same information for the V-SHARP-based and LBV-based pipelines is available in Supporting Information Figures S6C,D and S9. Biases smaller than |0.01| ppm (|⋅| denotes the absolute value) were considered negligible.
For the SHARP-based processing pipeline, the bias of TEdependent versus multi-echo χ was greater than |0.01| ppm in the CN, GP, and StrS at all TEs ( Figure 8A). For the TGVbased processing pipeline, the bias of TE-dependent versus multi-echo χ was greater than |0.01| ppm in the CN at TE 1 and TE 2 , and in the StrS at all TEs ( Figure 8B). In the CN for both processing pipelines and in the GP for the SHARP-based pipeline, the bias monotonically decreased with increasing F I G U R E 5 The means and SDs (error bars) of χ are shown in each ROI of the numerically simulated (A,B) and pooled healthy volunteer data (C,D) for the TE-dependent and multi-echo SHARP-based and TGV-based processing pipelines for QSM. In the numerical phantom, the groundtruth χ is also shown. In the healthy volunteers, the TEs at which the TE-dependent χ differed significantly from the multi-echo χ are denoted using the symbol ★ (P-value < .05). Abbreviations: CN, caudate nucleus; GP, globus pallidus; PCR, posterior corona radiata; PU, putamen; StrS, straight sinus; TH, thalamus; WM, white matter TEs, whereas in the StrS the bias decreased from the first to the third TE, then slightly increased again from the third to the fifth TE ( Figure 8). In all the ROIs, the SD of the bias decreased with increasing TEs (Figure 8), reflecting the reduced SD of χ at longer TEs ( Figure 5C,D).
In general, the bias was greater than |0.01| ppm at the same TEs at which there was a significant difference between TE-dependent and multi-echo QSM ( Figures 5C,D  and 8).

| DISCUSSION
This study aimed to compare processing pipelines for QSM that combine the GRE signal acquired at multiple TEs (multiecho QSM) with those that process each TE individually (TEdependent QSM) while applying four alternative strategies for Laplacian phase unwrapping and ΔB Bg removal (Figure 1). Multi-echo and TE-dependent QSM were compared by applying each pipeline to data simulated using a head numerical phantom (Supporting Information Figures S1 and S2) and images acquired in 10 healthy volunteers. The resulting χ values depended on several factors. The first was the LBM used for ΔB Bg removal: simultaneous phase unwrapping and ΔB Bg removal using SHARP gave the most accurate multi-echo χ, and one-step TGV gave the most accurate TE-dependent χ, even at short TEs ( Figure 5 and Supporting Information Figure S6). The second factor was the (ground-truth) χ of the ROI: larger reconstruction errors were observed in high-χ regions, such as the GP and the veins ( Figure 5 and Supporting Information Figure S6). The third factor, for TE-dependent QSM only, was the TE: all pipelines except TGV gave less accurate χ values than multi-echo QSM, especially at shorter TEs ( Figure 5 and Supporting Information Figure S6).
In the numerical simulations, TE-dependent QSM always suffered greater susceptibility noise than multi-echo QSM in all ROIs except the veins (Figures 3-7). The imprecision of TE-dependent QSM was closely dependent on noise and decreased with increasing TE (difference images in Figure 3 and SDs in Figure 5A,B). This finding is in line with previously published results on noise in QSM over TEs, 1,7 and, combined with the generally smaller RMSEs of multi-echo χ, suggests that multi-echo QSM offers the best processing pipeline performance across the whole brain. In the numerical phantom simulations, the accuracy of TE-dependent χ was similar at all TEs but decreased with increasing TE in the GP ( SHARP ) and the SSS ( SHARP and TGV ) ( Figure 5A,B). High-χ structures such as the SSS also showed the largest χ errors in the difference images, where they were increasingly visible at longer TEs (arrows in  Figures 3U-Z and 4U-Z). This result is in line with previous studies on simulated data, 6 which suggest that this TE dependence can arise from the failure of Laplacian phase unwrapping to accurately recover the true phase, especially in high χ regions. To address this limitation, future studies will investigate the application of alternative phase unwrapping methods. 40 The numerical phantom used in this study was characterized by an ROI-wise constant χ distribution. Therefore, in the numerical phantom simulations, an average χ close to True and a small SD of χ were positive indicators of a pipeline's accuracy and precision, respectively. However, brain regions in vivo are characterized by a heterogeneous χ distribution within each ROI. Thus, without knowledge of the true SD of χ, a smaller SD does not necessarily indicate a pipeline's greater precision in the healthy volunteers. For a comprehensive evaluation of the effect of overregularization by particular pipelines, future work will involve implementing a numerical phantom with a known heterogeneous χ distribution (SD) in each ROI.
In the healthy volunteers, in contrast with the numerical phantom simulations, the accuracy of TE-dependent QSM varied nonmonotonically over TEs (see the CN, PU, TH, and StrS in Figure 5C,D). As is discussed below, distinct factors could motivate this discrepancy in the deep GM and the venous ROIs.
In the deep GM ROIs of healthy volunteers, the greater inaccuracy of shorter versus longer TEs ( Figure 5C,D and Supporting Information Figure S6C,D) was possibly caused by additional tissue microstructure in vivo, which might cause a nonlinear temporal phase accumulation unresolvable by conventional QSM algorithms. 5,6 In particular, short -T * 2 components at a subvoxel level, which are preserved when using LBMs for phase unwrapping, could manifest as apparent nonlinearities at shorter TEs, thus explaining the TE dependence in the Laplacian-processed phase. 6 Attempting to explain the TE dependence of in vivo QSM, previous studies 5, 6 have described the measured χ using multicompartment models, showing that three tissue compartments might accurately define the major χ sources in deep GM, 5 and that a geometric model with axon, myelin, and extra-axonal space microcompartments might accurately define the major χ sources in WM. 41 However, in WM the TE dependence of phase 41 was observed in fiber bundles with highly ordered microstructure, which was not the focus of our study.
In the StrS, the nonlinear variation of TE-dependent QSM and the lower accuracy of shorter versus longer TEs compared with the simulations could perhaps be explained by the lack of flow compensation in the acquisition sequence, 42 which was a limitation of this study and particularly affected shorter TEs. As the maximum velocity of blood flow in healthy cerebral veins is about 25 cm/s, the maximum displacement of spins flowing in the StrS at increasing TEs were estimated as x TE i = 250 mm∕s ⋅ TE i [s]: 0.8, 2.1, 3.5, 4.8, and 6.2 mm. The exact encoded location of such flowing spins would depend on the StrS orientation relative to the applied encoding gradients, and on the total phase accumulated by these spins at each TE. Over time, the zero-crossing points of the first moment of the frequency-encoding gradient were increasingly closer to the TE (Supporting Information Figure S3), indicating that the total phase accumulated by flowing spins along this encoding direction would decrease over TEs. Therefore, the flow-induced χ error in the StrS would be expected to decrease at increasing TEs. Instead, flow-induced artifacts did not appear prominent in the multi-echo χ maps ( Figures 6G  and 7G and Supporting Information Figures S7G and S8G), possibly because they were mitigated by the nonlinear fitting method used to combine multiple TEs. 22 By assigning smaller weights to the phase at shorter TEs (at each TE, weights are inversely proportional to the level of noise in the phase), any flow-induced phase error or spatial shifts in the vessel location at the two shortest TEs were probably mitigated in the combined field map.
In the StrS, the greater discrepancy between multi-echo χ and TE-dependent χ at shorter TEs always resulted in a large χ estimation bias (Figure 8) corresponding to statistically significant differences between the two methods as measured by the sign test ( Figure 5C,D). Although smaller, a bias corresponding to significant differences in the sign test was also found at longer TEs ( Figures 5C,D and 8).
In ROIs other than the StrS, the difference between TGVbased TE-dependent and multi-echo QSM was only significant at the first or second TEs ( Figure 5D), whereas the difference between SHARP-based TE-dependent and multiecho QSM was sometimes significant also at longer TEs ( Figure 5C).
The results of the statistical analyses suggest that in TGV at shorter TEs, the estimation bias was mainly driven by the lower SNR of TE-dependent versus multi-echo QSM (larger limits of agreement at shorter TEs in Figure 5D). In SHARP , significant differences were sometimes found at longer TEs (see the CN, GP, and TH in Figure 5C), suggesting a potential underestimation of the true χ in addition to the lower SNR at shorter TEs. In both cases, unresolved tissue microstructure components in the deep-GM ROIs and the lack of flow compensation in the StrS might have also contributed to the systematic bias between TE-dependent and multi-echo QSM.
In both the simulated and in vivo data, shorter-TE QSM images appeared noisier and were potentially affected by the fast decaying signal components linked to tissue microstructure. Nonetheless, including all TEs always resulted in multiecho QSM images with less noise or smaller RMSEs (in the numerical simulations) compared to excluding shorter TEs (Supporting Information Figures S10-S12).
The numerical phantom simulations did not include modeling of receive bandwidth-induced distortions along the frequency-encoding direction, as the effect of such distortions on the ROI-based χ calculation was found to be negligible in healthy volunteers. Indeed, here the largest bandwidthinduced distortions, calculated as the ratio between the unwrapped total field map and the pixel bandwidth, 43 were located near the air-tissue interfaces and were within the [−1, 1] voxel range throughout the multi-echo field map.
In this study, all experiments were limited to one field strength (ie, 3 T). Although scanning at higher field strengths (eg, 7 T) would result in a higher phase at a given TE, this benefit may be counterbalanced by the decrease in T * 2 in tissue compartments already subject to rapid decay at 3 T. Thus, even if nonlinearities in the phase over TEs have been observed at both 3 T 6 and 7 T, 5 further work is needed to assess the scalability of our results to ultrahigh fields.
Finally, it must be noted that in the numerical phantom simulations, all processing pipelines underestimated True ( Figure 5A,B and Supporting Information Figure S6A,B). In QSM, some degree of underestimation is always expected, because the lack of MRI signal in the image background and zeros in the dipole kernel make the local field-to-susceptibility problem ill-posed. 2

| CONCLUSIONS
For studies on the TE-dependent variation of χ, the LBMs in the QSM processing pipeline should be carefully chosen to minimize biasing the results with LBM-related temporal variations. To this aim, a processing pipeline based on onestep TGV could provide slightly more accurate results over different TEs compared with processing pipelines based on SHARP, V-SHARP or LBV, and thus enable a more accurate investigation of TE-dependent χ variations deriving from tissue microstructure. This does not necessarily apply to the quantification of χ in the WM. For maximum accuracy, WM χ should be modeled as a tensor rather than a scalar value, because of the ordered microstructure found in WM. 41,44 Because there is no direct relationship between susceptibility tensor solutions and the single-orientation solution, 45 our results cannot directly inform on the accuracy and precision of tensor-based χ calculated from all the TEs versus single TEs of a multi-echo GRE protocol. Studies aimed at investigating the TE-dependent evolution of χ in the veins should implement flow compensation to avoid different flow-induced phase contributions at different TEs. Finally, the observed nonlinearity in the signal phase at shorter TEs highlights the need for accurate modeling of short T * 2 tissue components. Such modeling could also inform the design of more accurate digital brain phantoms, which is useful in methodological studies on TE-dependent QSM. strong background fields on the inferior side of the brain. The original Zubal numerical phantom with each region uniquely labeled from 1 to 125 is shown in (A). Here, the arrow points at the unrealistically flat interface between the brain and the nasal cavity. The susceptibility and total field calculated from (A) are respectively shown in (B) and (C). Figure (D) shows the gross anatomy of the nasal cavity and highlights the presence of a bony structure, called the cribriform plate of the ethmoid bone, between the brain and the nasal cavity. The susceptibility and total field calculated from (A) with the nasal cavity reshaped are respectively shown in Figures (E)  The means and SDs of χ are shown in each ROI of the numerically simulated (A,B) and pooled healthy volunteer data (C,D) for the TE-dependent and multi-echo V-SHARP-based and LBV-based processing pipelines for QSM. In the numerical phantom, the ground-truth χ is also shown. In the healthy volunteers, the TEs at which the TEdependent χ significantly different from the multi-echo χ are denoted using the symbol * (P-value < .05) FIGURE S7 TE-dependent vs. multi-echo χ maps in a representative healthy subject calculated using the V-SHARP-based processing pipelines. The same transverse and sagittal slices are shown for the susceptibility maps calculated using the V-SHARP-based multi-echo (A,G) and TE-dependent pipelines (B-F, H-L). The figure also shows the difference between each TE-dependent susceptibility map (M-V) and the corresponding multi-echo image. In all the sagittal images (G-L, R-V), the yellow arrowheads point at the same location in the StrS FIGURE S8 TE-dependent vs. multi-echo χ maps in a representative healthy subject calculated using the LBV-based processing pipelines. The same transverse and sagittal slices are shown for the susceptibility maps calculated using the LBVbased multi-echo (A,G) and TE-dependent pipelines (B-F, H-L). The figure also shows the difference between each TEdependent susceptibility map (M-V) and the corresponding multi-echo image. In all the sagittal images (G-L, R-V), the yellow arrowheads point at the same location in the StrS FIGURE S9 Bias of TE-dependent vs. multi-echo QSM.

ACKNOWLEDGMENTS
The mean and SD of the bias of TE-dependent vs. multi-echo QSM is shown for the V-SHARP-based (A) and LBV-based (B) processing pipelines. In both the bar plots, the gray band denotes the [−0.01 -0.01] ppm interval. If the mean of the bias was within this interval, the difference between the corresponding TE-dependent QSM pipeline and the multi-echo QSM pipeline was considered negligible FIGURE S10 Five-TE QSM vs. four-TE QSM in the numerical phantom. In the numerical phantom simulations and for the two best performing processing pipelines (TGV-based and SHARP-based), the figure shows the QSM images calculated using all TEs (A,D) or all TEs except the shortest (B,E). For each processing pipeline, the figure also shows the difference between the two QSM images (C,F). The RMSE value is reported for each QSM image FIGURE S11 Five-TE QSM vs. four-TE QSM in vivo. For one representative healthy volunteer and for the TGV-based pipeline, the figure shows the QSM image calculated using all TEs (A) or using all TEs except the shortest (B). The figure also shows the difference between the two QSM images (C) FIGURE S12 SD of five-TE QSM vs. four-TE QSM in vivo. For one representative healthy volunteer and the TGV-based processing pipeline (see Figure S11), the standard deviation of χ is shown in the whole brain and in each ROI considered for the healthy volunteer analysis