Multi-echo quantitative susceptibility mapping: how to combine echoes for accuracy and precision at 3 Tesla.

PURPOSE: To compare different multi-echo combination methods for MRI QSM. Given the current lack of consensus, we aimed to elucidate how to optimally combine multi-echo gradient-recalled echo signal phase information, either before or after applying Laplacian-base methods (LBMs) for phase unwrapping or background field removal.
METHODS: Multi-echo gradient-recalled echo data were simulated in a numerical head phantom, and multi-echo gradient-recalled echo images were acquired at 3 Tesla in 10 healthy volunteers. To enable image-based estimation of gradient-recalled echo signal noise, 5 volunteers were scanned twice in the same session without repositioning. Five QSM processing pipelines were designed: 1 applied nonlinear phase fitting over TEs before LBMs; 2 applied LBMs to the TE-dependent phase and then combined multiple TEs via either TE-weighted or SNR-weighted averaging; and 2 calculated TE-dependent susceptibility maps via either multi-step or single-step QSM and then combined multiple TEs via magnitude-weighted averaging. Results from different pipelines were compared using visual inspection; summary statistics of susceptibility in deep gray matter, white matter, and venous regions; phase noise maps (error propagation theory); and, in the healthy volunteers, regional fixed bias analysis (Bland-Altman) and regional differences between the means (nonparametric tests).
RESULTS: Nonlinearly fitting the multi-echo phase over TEs before applying LBMs provided the highest regional accuracy of χ $$ \chi $$ and the lowest phase noise propagation compared to averaging the LBM-processed TE-dependent phase. This result was especially pertinent in high-susceptibility venous regions.
CONCLUSION: For multi-echo QSM, we recommend combining the signal phase by nonlinear fitting before applying LBMs.

INTRODUCTION

MRI QSM aims to determine the underlying spatial distribution of tissue magnetic susceptibility () from gradient‐recalled echo (GRE) phase data (): where is a vector of image space coordinates, the proton gyromagnetic ratio, TE the echo time, the ‐induced total field perturbation along the scanner's z axis, and the TE‐independent phase offset at a nominal TE = 0 ms.

For QSM, the acquired phase must be spatially (single‐echo data) or spatiotemporally (multi‐echo data) unwrapped to resolve aliasing. The unwrapped is proportional to (Equation 1), which is a combination of background () and local field contributions (): are induced by the global geometry, air–tissue interfaces, and any field inhomogeneities. reflect the tissue inside the region of interest (ROI), for example, the brain. For QSM, must be removed from . The resulting map is in the following relationship with : where is the magnetic dipole and denotes a spatially dependent convolution. Based on Equation 3, the local distribution of tissue , that is, the QSM map, is calculated by solving an ill‐posed ‐to‐ problem.

Recently, the QSM community critically reviewed how to best perform phase unwrapping and removal. Moreover, it promoted 2 challenges to compare algorithms for ‐to‐ inversion, but a consensus has yet to be reached. , A further open question toward QSM standardization is how and at which stage of the processing pipeline multi‐echo data from different echoes should be combined. This question is relevant because phase unwrapping, removal, or both, are often performed using Laplacian‐based methods (LBMs). ,

Laplacian phase unwrapping aims to calculate the ‐aliasing‐free phase as : with and the forward and inverse Laplace operators, and the aliased phase. Laplacian removal relies on the harmonicity of inside the ROI (i.e., ) and aims to solve this Equation :

The inverse discrete Laplace operator is not well defined and requires regularization, which is equivalent to spatially high‐pass filtering the phase or local field. However, the known varying frequency content at different TEs , could lead to different degrees of LBM‐induced high‐pass filtering at different echoes, alter the linearity of Equation 1, and thus introduce inaccuracies in the estimated and maps. To investigate the issue, this study aimed to compare existing strategies for combining the multi‐echo signal phase (see the Theory section for further details) when using LBMs for phase unwrapping or removal.

Five processing pipelines for QSM were designed incorporating LBMs for both phase unwrapping and removal and combining the signal from different TEs by fitting or averaging before or after applying LBMs. These pipelines were applied to both numerically simulated data and images acquired in vivo. Results from each pipeline were compared qualitatively by visual inspection and quantitatively via analysis of the regional bias and precision, as well as noise propagation.

THEORY

Multi‐echo combination

Previous studies employing multi‐step reconstruction pipelines for QSM at 3 Tesla have combined the signal from multiple echoes by either averaging , , , or fitting , , , before or after applying LBMs. For multi‐echo combination in QSM, this study focuses on approaches based on weighted averaging , , , or complex nonlinear fitting (NLFit) , that outperform approaches based on unweighted averaging or linear fitting.

Fitting

Multi‐echo combination via nonlinear fitting (NLFit) has formulated the temporal evolution of the complex signal as a nonlinear least squares problem : where denotes the acquired complex signal, the signal magnitude, the signal phase (Equation 1), and the i‐th echo time. This approach aims to mitigate noise in by correctly modeling as normally distributed the noise in the real and imaginary parts of the complex signal. Unlike weighted‐averaging–based approaches, NLFit enables estimating . Notably, the input phase to nonlinear fitting is minimally processed because it only requires to be temporally unwrapped, thus avoiding the application of LBMs before multi‐echo combination.

Weighted averaging

Multi‐echo combination (with echoes) via weighted averaging has been performed using either TE‐based weighting factors , or phase SNR‐based weighting factors.

TE‐based weighted averaging (TE‐wAvg) , accounts for the phase at shorter TEs being affected by larger noise levels than at longer TEs. TE‐wAvg calculates a combined as with weights equal to:

TE‐wAvg requires temporal unwrapping of the input multi‐echo phase, resulting in a combined , which still contains contributions.

SNR‐based weighted averaging (SNR‐wAvg) accounts for different tissue types reaching optimal SNR at different TEs and calculates a combined local field map () as: with weights equal to:

In Equation 10, denotes the map of voxel‐wise transverse relaxation rates, which is related to the signal magnitude by: with the initial transverse magnetization. SNR‐wAvg requires performing temporal and spatial unwrapping as well as background field removal on the input multi‐echo phase, resulting in a background field–free field map.

Alternatively, a distinct map has been calculated at each TE, and multi‐echo combination of over time has been performed via weighted averaging using magnitude‐based weighting factors (Susc‐wAvg). Based on the inverse proportionality of the phase noise and magnitude SNR, this method aims to improve the SNR of the combined map as: with weights equal to:

Noise propagation

Previous studies using fitting or SNR‐wAvg have calculated expressions for the noise in the total field map (). For TE‐wAvg or SNR‐wAvg, expressions for have not been calculated and were therefore derived here. For each multi‐echo combination method, expressions for noise propagation from to the corresponding and images were also derived here.

Based on error propagation, the noise in calculated using a linear least squares fitting approach is (see Supporting Information, Section 1):

Equation 14 also corresponds to the a priori noise estimate found in nonlinear least squares fitting, thus .

Based on Equations (7), (8), (9), (10), calculated using TE‐wAvg and calculated using SNR‐wAvg, respectively, have variances equal to: Assuming that noise in the single‐echo phase is temporally uncorrelated and based on error propagation, the noise in / calculated using TE‐wAvg/SNR‐wAvg is, respectively, equal to: Equations 17 and 18 omit the dependency because they combine the multi‐echo phase voxel by voxel.

Noise in the local field map

Based on error propagation and the orthogonality of and the in the ROI, for example, the brain: where denotes an “if and only if” relationship, and is the worst‐case scenario.

Noise in the susceptibility map

Due to the circular convolution theorem, the deconvolution operation in Equation 3 can be performed by point‐wise division in the Fourier domain. Thus, if the regularized inverse dipole kernel in k‐space can be analytically derived independent of the Fourier transforms of or , can be calculated as (see Supporting Information, Section 2): Where and , respectively, denote the direct and inverse Fourier transforms, and denotes k‐space coordinates. Deriving an analytical expression for is possible, for example, when considering thresholded k‐space division or the Tikhonov‐regularized minimal norm solution. , ,

For weighted averaging of TE‐dependent (Equations 12 and 13), based on phase error propagation over time and Equation 20, equals Based on Equations 12, 13, 20, and 21 (see Supporting Information, Section 3), is the square root of: Notably, Equation 22 requires analytically describing (Equation 20).

METHODS

Where not otherwise stated, image analysis was performed on a 64‐bit Microsoft (Redmond, WA) Windows 11 Pro operating system (Intel(R) (Santa Clara, CA) Core(TM) i5‐9400 CPU@2.90GHz processor; 16 GB RAM) using MatLab (R2021b, MathWorks, Natick, MA). Preliminary versions of this study were presented at the 2016 and 2018 annual meetings of the International Society for Magnetic Resonance in Medicine. ,

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 in 2 centers (University College London Hospital and Queen Square Multiple Sclerosis Centre, University College London) equipped with the same 3 Tesla MRI system (Philips Achieva, Philips Healthcare, Best, NL; 32‐channel head coil). Five subjects were acquired in each center. All volunteers provided written informed consent, and the local research ethics committees approved the experimental sessions. Images were acquired using a transverse orientation, FOV = 240×240×144 mm3, voxel size = 1‐mm isotropic, flip angle = 200, TR = 29 ms, 5 evenly spaced echoes (TE1/TE spacing = 3/5.4 ms), bandwidth = 270 Hz/pixel, SENSE factors = 2/1.5, flyback gradients = on, no flow compensating gradients, total scan duration = 04:37 min:s. Five subjects were scanned twice within the same session without repositioning to enable image‐based calculation of magnitude and phase SNRs.

Data simulation from a numerical head phantom

To ensure the availability of ground‐truth values against which to test the accuracy of QSM pipelines, a Zubal numerical head phantom was used with the following ROIs: the caudate nucleus (CN), globus pallidus (GP), putamen (PU), thalamus (TH), superior sagittal sinus, gray and white matter (GM and WM), and CSF. To match the acquisitions in vivo, the original 1.5‐mm isotropic phantom was resampled to a 1‐mm isotropic resolution with matrix size = 384×384×192 voxels. Compared to our previous study, the numerical phantom was updated to achieve realistic regional means ± SDs for , the proton density (), and the transverse relaxation rate () (see Supporting Information, Section 4). Simulated multi‐echo complex data were generated based on these ground‐truth spatially variable , , and distributions, as: with and , respectively, described by Equations 1 and 11, and TEs matched to the in vivo acquisitions. Random zero‐mean Gaussian noise with a SD = 0.07 was added to the real and imaginary parts of the noise‐free signal independently. , The random noise matrix was regenerated at each TE.

Data preprocessing

A brain mask was calculated for each subject by applying FSL brain extraction tool , with robust brain center estimation (threshold = 0.3) to the magnitude image at the longest TE. This choice of TE accounted for the greater amount of signal dropout near regions of high‐ gradients compared to shorter TEs.

A whole‐brain mask for the Zubal phantom was calculated by applying FSL brain extraction tool with robust brain center estimation (threshold = 0.5) to the map of the numerical phantom.

Processing pipelines for QSM

Five distinct processing pipelines (Figure 1) were applied to both the numerically simulated and the healthy volunteer data, and the time required to run each pipeline was measured using MatLab's stopwatch timer (MathWorks). Three of these pipelines (NLFit, TE‐wAvg, and SNR‐wAvg) combined the phase across TEs at different stages before performing the ‐to‐ inversion. Two other pipelines (Susc‐wAvg and Susc‐total generalized variation (TGV)‐wAvg) first calculated a distinct map at each TE and then combined the maps. The following paragraphs describe each processing pipeline in detail.

FIGURE 1

Processing pipelines for multi‐echo QSM. For each multi‐echo combination method (NLFit, TE‐wAvg, SNR‐wAvg, Susc‐wAvg, and Susc‐TGV‐wAvg), the processing steps are described as separate processing streams. NLFit, nonlinear phase fitting; SNR‐wAvg, SNR‐weighted phase averaging; Susc‐wAvg and Susc‐TGV‐wAvg, magnitude‐weighted susceptibility averaging; TE‐wAvg, TE‐weighted phase averaging

The NLFit pipeline first combined the complex GRE signal by nonlinear fitting over TEs using the Cornell QSM software package's Fit_ppm_complex function. It then applied simultaneous spatial phase unwrapping and removal using sophisticated harmonic artifact reduction for phase data (SHARP), a direct solver of Equation 5.

SHARP was chosen because it has been widely used in the literature on QSM and is both robust and numerically efficient. Moreover, in our recent study comparing multi‐echo and TE‐dependent QSM, a multi‐echo pipeline incorporating SHARP gave highly accurate multi‐echo QSM values. SHARP was applied using the minimum‐size 3‐voxel isotropic 3D Laplacian kernel, a threshold for truncated singular value decomposition equal to 0.05, and a brain mask eroded by 5 voxels.

The TE‐wAvg processing pipeline first applied Laplacian unwrapping to the phase at each TE using a threshold for truncated singular value decomposition equal to 10−10 (i.e., the default value in Ref. 29). Second, it calculated by averaging the unwrapped phase according to Equations 7 and 8. , Then, it calculated by applying SHARP with the same parameters as in NLFit.

The SNR‐wAvg processing pipeline first applied simultaneous phase unwrapping and removal to the phase at each TE using SHARP with the same parameters as in NLFit. An map was calculated by voxel‐wise fitting Equation 11 using MatLab's nlinfit function, where initial values for and were calculated by linearly fitting the log‐linearized version of Equation 11. The SNR‐wAvg pipeline then calculated by averaging the unwrapped and background field‐free phase according to Equations 9 and 10.

In the NLFit, TE‐wAvg, and SNR‐wAvg pipelines, ‐to‐ inversion was performed using Tikhonov regularization with correction for susceptibility underestimation and using the L‐curve method to determine the optimal value for the regularization parameter. , , This inversion method was chosen because it is computationally efficient and substantially reduces streaking artifacts relative to thresholded k‐space division.

The Susc‐wAvg processing pipeline calculated a distinct map at each TE by applying simultaneous phase unwrapping and removal using SHARP as in NLFit and performing the ‐to‐ inversion using Tikhonov regularization as in NLFit, TE‐wAvg, and SNR‐wAvg. This pipeline then calculated a combined map according to Equations 12 and 13.

The Susc‐TGV‐wAvg processing pipeline applied 1‐step TGV to the phase at each TE and then calculated a combined map as in Susc‐wAvg. The TGV method was tested because it avoids stair‐casing artifacts in the resulting map while correctly preserving structural borders. Moreover, in our recent study comparing multi‐echo and TE‐dependent QSM, TGV provided highly accurate TE‐dependent QSM images. TGV (v1.0.0_20210629) was run in Neurodesk (v20220302, https://neurodesk.github.io/) with the default parameter values , which are optimal for medical imaging applications.

ROI segmentation in the healthy volunteer images

Regional values were compared within the simulated and in vivo data, similarly to our previous study. ROIs similar to those in the numerical phantom were segmented in vivo: the CN, GP, PU, TH, posterior corona radiata (PCR) as a WM ROI, and the straight sinus (StrS) as a venous ROI. Briefly, for each subject, the CN, GP, PU, TH and PCR were segmented based on the Eve atlas, for which the GRE magnitude image was aligned to each subject's fifth‐echo magnitude image using NiftyReg (NiftK v14.11.0) , (TEEve/TE5 = 24/24.6 ms). The quality of ROI alignment was assessed by visual inspection. The ITK‐SNAP (active contour segmentation tool was used to segment the StrS over several slices based on the fifth‐echo magnitude image, which presented the best contrast between the StrS and the surrounding brain tissue.

Quantitative evaluation of the measured

In the numerical phantom simulations, each QSM pipeline's performance relative to the ground truth was visually assessed by calculating a difference image between the corresponding map and . Here, referred to the ground‐truth local field calculated using the reference scan method, , and to the ground‐truth magnetic susceptibility distribution with realistic regional means ± SDs of (Supporting Information Figure S1B). Means and SDs of were calculated for each pipeline in each ROI with , that is, the CN, GP, PU, TH, superior sagittal sinus, and WM. The RMSEs of both and relative to were also calculated throughout the brain volume. Notably, RMSEs for could not be calculated for the 1‐step Susc‐TGV‐wAvg pipeline.

In the volunteers, due to to the lack of a ground truth, representative susceptibility difference images were calculated relative to because the NLFit pipeline performed multi‐echo combination at the earliest possible stage; and relative to , because the TE‐wAvg pipeline had the lowest local field RMSE in the numerical phantom simulations (see the Results). Regional means and SDs of were calculated for each processing pipeline and compared against values in subjects of a similar age from the QSM literature. RMSEs could not be calculated because of the lack of a ground truth. For visualization purposes, the pooled averages and SDs were calculated after verifying that all intrasubject SDs of were larger than the intersubject SD of .

Noise propagation maps

Only the healthy volunteers scanned twice were considered for this analysis. To enable image SNR calculation, in 1 healthy volunteer, five 20×20‐voxel ROIs were drawn on a sagittal slice of the first‐echo magnitude image, including both the GM and WM and excluding regions with artifacts induced by SENSE, motion, or flow. All 5 ROIs were applied across the other 4 volunteers by using rigid alignment transforms (NiftyReg ) between the first‐echo magnitude images.

In each subject, ROI‐based magnitude (), magnitude noise (), and phase noise values () were calculated based on the SNR difference method : Here, and , respectively, denote the magnitude/phase images from the first and second scan, and is the 2D SENSE factor calculated by multiplying the SENSE factors applied along the 2 phase encoding directions. The values calculated based on Equations (24), (25), (26) were averaged across the 5 ROIs to calculate summary values of magnitude, magnitude noise, and phase noise at each TE.

For both the numerical phantom simulations and the healthy volunteers, a phase noise map was calculated at each TE as , : where was a constant equal to 1 (by definition) in the numerical phantom simulations and equal to: in the healthy volunteers. Notably, calculating phase noise analytically as in Equation 27 enabled the direct comparison of noise propagation between numerical simulations and data acquired in vivo.

For the NLFit, TE‐wAvg, and SNR‐wAvg pipelines, the map was calculated based on the multi‐echo maps according to Equations 14 and 17–19. Then, the map was calculated according to Equation 20 with Tesla and the Tikhonov‐regularized inverse magnetic dipole kernel with regularization parameter : and correction for underestimation. For the Susc‐wAvg pipeline, the map was calculated based on Equation 22. maps were not calculated for the Susc‐TGV‐wAvg pipeline because TGV estimates iteratively and an analytical expression for could not be derived.

To compare maps across pipelines, a line profile was traced in the same location of all images, and the value of each voxel along this profile was plotted. To compare the noise intensity and its variability between processing pipelines, the mean and SD of this representative line profile were also calculated.

Statistical analysis

Statistical analyses were performed based on the healthy subject data. For each ROI and each pair of pipelines, Bland–Altman analysis of the average was used to assess if pairs of pipelines systematically produced different results. For each ROI, statistically significant differences between pipelines 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. All statistical tests were 2‐tailed, and an uncorrected P value < 0.05 was considered significant.

RESULTS

Pooling of measurements

For each ROI and each processing pipeline, all intrasubject SDs of were larger than the intersubject SD. Thus, pooled means and SDs were calculated.

Performance of pipelines for multi‐echo QSM

The Susc‐wAvg and Susc‐TGV‐wAvg processing pipelines were the longest to run (Supporting Information Table S1) because they calculated a QSM map at each TE. The longer processing times required for the numerical phantom data were linked to the larger matrix size (384×384×192) compared to acquisitions in vivo (240×240×144).

In the numerical phantom, Figure 2 shows the ground‐truth (A, G), the QSM images calculated by each processing pipeline (B–L), their difference relative to the ground truth (M–V), and the RMSEs of throughout the brain volume (bottom row). Analogous results for are shown in Supporting Information Figure S2. In the numerical phantom simulations, the map calculated using the NLFit pipeline had the largest RMSE (109.4%) followed, in decreasing order, by the SNR‐wAvg (94.3%), TE‐wAvg (93.5%), Susc‐wAvg (91.5%), and Susc‐TGV‐wAvg pipelines (80.4%) (Figure 2, bottom row). The map calculated using the Susc‐wAvg pipeline had the largest RMSE (average across TEs: 85.0%), followed in decreasing order by the NLFit (81.2%), TE‐wAvg, and Susc‐wAvg pipelines (both 71.8%).

FIGURE 2

maps calculated using distinct multi‐echo combination methods in the numerical phantom simulations. The same transverse and sagittal slices are shown for the ground‐truth susceptibility map (A, G), and for the susceptibility maps calculated using NLFit (B, H), TE‐wAvg (C, I), SNR‐wAvg (D, J), Susc‐wAvg (E, K), and Susc‐TGV‐wAvg (F, L). The figure also shows the difference between each susceptibility map and the ground truth (M–V). The bottom row shows the RMSEs of for each pipeline. In all the sagittal images (G–L, R–V), the yellow and blue arrowheads, respectively, point at the same posterior and anterior locations in the superior sagittal sinus

In the numerical phantom and for each processing pipeline, Figure 3A shows the regional means and SDs of . The error between the calculated and ground‐truth appeared similar for all processing pipelines, although slightly larger SDs were always observed for the NLFit pipeline (Figures 2M–V and 3A). The superior sagittal sinus, which was the structure with the largest , showed the largest susceptibility errors for all processing pipelines (arrowheads in Figures 2G–L, R–V, and 3A).

FIGURE 3

Means and SDs of in the phantom and healthy volunteer ROIs. The means and SDs (error bars) of are shown in each ROI of the numerically simulated (A) and pooled healthy volunteer data (B) for each processing pipeline. In the numerical phantom, the ground‐truth is also shown. In the healthy volunteers, significant differences between pairs of pipelines are denoted using the symbols * ( P value <0.05) and ** ( P value <0.01) . CN, caudate nucleus; GP, globus pallidus; PCR, posterior corona radiata; PU, putamen; StrS, straight sinus; TH, thalamus; WM, white matter; SSS, superior sagittal sinus

For one representative volunteer, Figures 4A–J show the susceptibility images calculated using each processing pipeline. Additionally, differences images are shown relative to the (K‐R) and the maps (S‐Z). Here, susceptibility differences between processing pipelines were most prominent in the StrS (arrowheads in Figures 4F–J, O–R, W–Z). In the healthy volunteers, Figure 3B shows the pooled regional means and SDs of calculated by each processing pipeline. The average measured in the deep‐GM ROIs and in the PCR had values within the ranges reported by previous studies , , , , , : 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. In the StrS, only had an average value close to the previously reported range for venous blood, namely, 0.17–0.58 ppm , , , (Figure 3B).

FIGURE 4

maps calculated using distinct multi‐echo combination methods in a representative healthy volunteer. The same transverse and sagittal slices are shown for the susceptibility maps calculated using NLFit (A, F), TE‐wAvg (B, G), SNR‐wAvg (C, H), Susc‐wAvg (D, I), and Susc‐TGV‐wAvg (E, J). The figure also shows the differences between the TE‐wAvg, SNR‐wAvg, Susc‐wAvg, and Susc‐TGV‐wAvg maps and the NLFit map (K‐R); and the differences between the NLFit, SNR‐wAvg, Susc‐wAvg, and Susc‐TGV‐wAvg maps and the TE‐wAvg map (S‐Z). In all the sagittal images (F‐J, O‐R, W‐Z), the yellow arrowheads point at the same location in the StrS

In the CN, GP, and venous ROIs, the regional measured in vivo had a relative accuracy across pipelines similar to the numerical phantom simulations: the NLFit and Susc‐TGV‐wAvg pipelines, respectively, provided the highest and lowest means of , whereas the TE‐wAvg, SNR‐wAvg, and Susc‐wAvg pipelines provided intermediate values (Figure 3). Slightly different trends were observed in the PU, TH, and WM ROIs, where, in vivo, the NLFit and Susc‐TGV‐wAvg pipelines provided the lowest (absolute) means of ; and the TE‐wAvg, SNR‐wAvg, and Susc‐wAvg pipelines provided higher values (Figure 3B). In the numerical phantom simulations, the NLFit pipeline always had the largest SD of (Figure 3A). In contrast, in vivo, the NLFit pipeline had the smallest SD of in the CN, PU, TH, and PCR, and SDs of comparable to other pipelines in the GP and StrS (Figure 3B).

Phase noise propagation into the maps

Figures 5 and 6 show the estimated maps and profiles in the numerical phantom simulations and a representativehealthy subject, respectively. All images contained some degree of streaking artifacts, especially in vivo, because these are a manifestation of error propagation from to caused by the dipole kernel null space. In the phantom, the NLFit and Susc‐wAvg pipelines had the line profiles with the lowest means and SDs (Figures 5A,D,E,H,I). The NLFit and SNR‐wAvg pipelines always resulted in the lowest streaking artifacts burden (Figures 5A,E,C,G,I and 6A,E,C,G,I). Streaking artifacts were more severe in the TE‐wAvg and Susc‐wAvg pipelines, especially near high‐ venous structures.

FIGURE 5

Susceptibility noise profiles in the numerical phantom simulations. The same sagittal slice is shown for susceptibility noise () images calculated using the NLFit, TE‐wAvg, SNR‐wAvg, and Susc‐wAvg pipelines (A–D). The susceptibility noise is plotted for a line profile traced on the images (E–H), including the high‐ superior sagittal sinus. All line profiles are also shown combined in the same plot (I). The mean and the SD of are shown for each line profile

FIGURE 6

Susceptibility noise profiles in a representative healthy volunteer. The same sagittal slice is shown for susceptibility noise () images calculated using the NLFit, TE‐wAvg, SNR‐wAvg, and Susc‐wAvg pipelines (A–D). The susceptibility noise is plotted for a line profile traced on the images (E–H) including the high‐ StrS. All line profiles are also shown combined in the same plot (I). The mean and SD of are shown for each line profile

In the healthy volunteers, for all processing pipelines and ROIs, the Shapiro–Wilk test always rejected the hypothesis of normally distributed paired differences of . Therefore, pairwise comparisons between pipelines were always evaluated using the nonparametric sign test. Significant differencesbetween pipelines are shown in Figure 3B, whereas the between‐pipeline biases are shown in Figure 7. A lower threshold equal to |0.01| ppm (|·| denotes the absolute value) was chosen because, below this level, interpipeline differences cannot be disentangled from the intrapipeline variability (i.e., the regional SD).

FIGURE 7

Bias between multi‐echo pipelines for QSM in each ROI. The mean and SDs (error bars) of the bias are shown in each healthy volunteer ROI for all pairs of multi‐echo processing pipelines. 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 pair of QSM pipelines was considered negligible. ROI, region of interest

A bias greater than |0.01| ppm was observed for the NLFit pipeline relative to all other pipelines in the CN and StrS (Figure 7). A bias greater than |0.01| ppm was observed for the Susc‐TGV‐wAvg pipeline relative to all other pipelines in the GP, and relative to the TE‐wAvg and Susc‐wAvg pipelines in the StrS (Figure 7). These results suggests that, for accurate quantification, some of the significant differences detected by the sign test, for example, between the TE‐wAvg, SNR‐wAvg, and Susc‐wAvg pipelines, may be negligible (Figures 3B and 7).

DISCUSSION

Aiming to elucidate the optimal strategy for multi‐echo combination for QSM, this study compared multi‐echo combination methods applied at different stages of the QSM processing pipeline before or after LBMs for phase unwrapping or background field removal. Each pipeline was applied to numerically simulated data and images from healthy volunteers.

The higher relative accuracy of the NLFit pipeline in the numerical phantom simulations and in vivo suggests that, for QSM, combining the temporally unwrapped multi‐echo phase before applying LBMs for spatial phase unwrapping or removal is preferable to averaging the TE‐dependent LBM‐processed phase or . This suggestion appears to conflict with the higher RMSEs associated with the NLFit pipeline compared to other pipelines in the phantom simulations (Figures 2 and Supporting Information Figure S2). However, the RMSE jointly reflects systematic and random errors because it measures the bias between the estimated and true value and also reflects the variability of the estimated relative to its average value. Thus, the RMSE must always be interpreted in combination with complementary measurements of bias and precision. Furthermore, RMSEs of are difficult to interpret. In contrast with RMSEs of , they allow comparison of different pipelines without the effect of ‐to‐ inversion. However, they are voxel‐based measures based on a signal that is intrinsically nonlocal because variations extend beyond the anatomical region of shift that generated them. Thus, the best set of metrics for comparing images generated by a processing pipeline for QSM is still an active area of research. ,

There are several potential explanations as to why LBMs applied before multi‐echo combination reduce the overall accuracy of QSM. Firstly, in contrast with path‐based or region‐growing–based phase unwrapping, Laplacian phase unwrapping usually removes some components from the input phase image. Thus, the consistently lower accuracy of calculated using the TE‐wAvg processing pipeline was probably driven by the incorrect assumption that the TE‐dependent LBM‐unwrapped phase corresponded to the true unwrapped phase. LBMs for removal also applied truncated singular value decomposition (with a larger truncation threshold), but here the high‐pass filtering effect was expected because background fields are slowly varying. Finally, the similar accuracy and values of calculated using the TE‐wAvg, SNR‐wAvg, and Susc‐wAvg pipelines suggests a negligible difference between averaging the Laplacian unwrapped phase over TEs before (TE‐wAvg pipeline) or after TE‐dependent Laplacian removal (SNR‐wAvg and Susc‐wAvg pipelines).

In the numerical phantom simulations, all processing pipelines resulted in higher SDs of compared to the ground truth, suggesting that noise in the multi‐echo signal phase was amplified by all pipelines. This result is in line with the known noise amplification of ill‐posed inverse problems. However, the estimated SDs of varied across pipelines. In the simulations, the Susc‐TGV‐wAvg pipeline had the smallest SDs of (Figure 3A). However, in vivo, the NLFit pipeline generally had the smallest SDs of (Figure 3B). Both the TGV reconstruction pipeline and shortcomings of the numerically simulated data could explain these discrepancies in the performance of the Susc‐TGV‐wAvg pipeline. The numerically simulated data were generated based on a digital phantom which, despite varying regional values in a realistic fashion (see Supporting Information, Section 4 and Figure S1), ultimately still appeared as a smooth piece‐wise constant model (Figure 2A,G). As previously observed, piece‐wise constant geometries allow good recovery of the underlying distribution using TGV‐based algorithms because the piece‐wise constant constraints exactly match the underlying distribution. However, in regions with flow, anisotropic distributions, or microstructure, these numerical models are likely to depart from a realistic representation of the tissue . To overcome the limitations of this assumption, future studies could exploit a newly developed realistic head phantom for QSM, which does not have a piecewise constant distribution and incorporates microstructural effects.

In both numerical simulations and healthy volunteers, based on the line profiles traced on the noise maps and streaking artifact reduction, the NLFit pipeline had better noise mitigation (Figures 5 and 6). This result suggests that combining the temporally unwrapped multi‐echo phase by nonlinear complex fitting, designed to account for noise in the complex signal, results in better noise management than combining the multi‐echo phase by averaging. As previously shown, errors in the combined field map mainly result from both noise in the signal phase and phase unwrapping errors near high‐ regions (e.g., the veins). Both sources of error were successfully managed by nonlinear complex fitting, as shown by the dramaticreduction of streaking artifacts in Figure 6A. In line with previous observations, this result also suggests that the regularization strategy employed by the ‐to step can mitigate artifactual streaking errors only when major sources of error in the field map have been tightly constrained. At visual inspection in vivo, the SNR‐wAvg was the second‐best pipeline for the mitigation of streaking artifacts (Figure 6). Thus, if applying nonlinear complex fitting is not possible, averaging using the SNR‐weighting–based method could offer the best alternative for noise reduction.

All these indications do not necessarily apply to estimation in WM tissue. Indeed, due to WM's ordered microstructure, a comprehensive estimation of in WM requires acquiring gradient‐recalled echo images at multiple head orientations and modeling as a tensor. , Thus, further studies are needed to evaluate the applicability of these results to WM tissue.

In the present study, all experiments were limited to one field strength (i.e., 3 Tesla). Because tissue relaxation times (e.g., ) shorten with increasing field strength, but the signal phase at a given TE increases, further work is needed to assess the relevance of these results at ultrahigh fields. Finally, it must be noted that in the numerical phantom simulations, all processing pipelines underestimated (Figure 3). However, in QSM some degree of underestimation is always expected due to the ill‐posed nature of the ‐to‐ inverse problem.

CONCLUSION

The higher accuracy of regional values and better noise management of the NLFit pipeline suggest that, for QSM, combining the multi‐echo phase by nonlinearly fitting over TEs before applying LBMs is preferable to combining the TE‐dependent LBM‐processed phase or by averaging.

CONFLICT OF INTEREST

The authors declare no conflicts: AstraZeneca was not involved in any aspect concerning this work and has not influenced its content or the decision to submit it for publication.

FUNDING INFORMATION

Supported by the UK Engineering and Physical Sciences Research Council (EPSRC), award number: 1489882 (e.b.); by the EPSRC‐funded UCL Centre for Doctoral Training in Medical Imaging, grant EP/L016478/1 (a.k.), and the Department of Health's National Institute for Health Research funded Biomedical Research Centre at University College London Hospitals (a.k.); by the UCL Leonard Wolfson Experimental Neurology Centre, grant PR/ylr/18575 (d.l.t)

The Queen Square MS Centre, where part of the MRI scans for this work were performed, is supported by grants from the UK MS Society and by the National Institute for Health Research University College London Hospitals Biomedical Research Centre (UCLH/BRC). F. Grussu was supported by PREdICT, a study at the Vall d'Hebron Institute of Oncology in Barcelona funded by AstraZeneca (f.g.), and funding from the postdoctoral fellowships program Beatriu de Pinós (2020 BP 00117), funded by the Secretary of Universities and Research, Government of Catalonia (f.g.)

Section 1. Noise in total field map calculated by fitting

Section 2. Noise propagation from the local field to the susceptibility map

Section 3. Noise propagation in the Susc‐wAvg pipeline

Section 4. Numerical phantom simulations

Figure S1. Properties of the numerical phantom. in arbitrary units (a. u.), in ms and in parts per million (ppm) assigned to various ROIs in the numerical phantom are shown in (A). The location of these ROIs is shown in the (B) and maps (C) of the numerical phantom

Table S1 . Total image processing time. The table shows the time required to run each pipeline for the numerical phantom simulations and data acquired in vivo. The time reported for the SNR‐wAvg pipeline does not include the time required for mapping, as optimizing this step was outside the scope of the present study. For the Susc‐TGV‐wAvg pipeline, approximate timings are reported because image reconstruction was performed using first Neurodesk (QSM calculation at each TE) and then MatLab (multi‐echo combination)

Figure S2 . maps calculated using distinct multi‐echo combination methods in the numerical phantom simulations. The same transverse and sagittal slices are shown for the ground‐truth local field map (A, G), and for the local field maps calculated using NLFit (B, H), TE‐wAvg (C, I), SNR‐wAvg (D, J), and Susc‐wAvg at each TE (E, O). The figure also shows the difference between each local field map and the ground truth (P‐E2). The bottom row shows the RMSEs of for each pipeline

Click here for additional data file.