NeuroImage 116 (2015) 214–221

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Improved estimation of myelin water fraction using complex model fitting Yoonho Nam a,b, Jongho Lee b,⁎, Dosik Hwang a, Dong-Hyun Kim a,⁎⁎ a b

Department of Electrical and Electronic Engineering, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 120-749, Republic of Korea Department of Electrical and Computer Engineering, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 151-744, Republic of Korea

a r t i c l e

i n f o

Article history: Received 14 January 2015 Accepted 17 March 2015 Available online 6 April 2015 Keywords: Myelin water White matter T⁎2 relaxation Fiber orientation Frequency offset B0 field inhomogeneity

a b s t r a c t In gradient echo (GRE) imaging, three compartment water modeling (myelin water, axonal water and extracellular water) in white matter has been demonstrated to show different frequency shifts that depend on the relative orientation of fibers and the B0 field. This finding suggests that in GRE-based myelin water imaging, a signal model may need to incorporate frequency offset terms and become a complex-valued model. In the current study, three different signal models and fitting approaches (a magnitude model fitted to magnitude data, a complex model fitted to magnitude data, and a complex model fitted to complex data) were investigated to address the reliability of each model in the estimation of the myelin water signal. For the complex model fitted to complex data, a new fitting approach that does not require background phase removal was proposed. When the three models were compared, the results from the new complex model fitting showed the most stable parameter estimation. © 2015 Elsevier Inc. All rights reserved.

Introduction Myelin water imaging (MWI) is a noninvasive technique in exploring demyelinating diseases such as multiple sclerosis (Mackay et al., 1994). Conventional approaches use multi-echo spin echo (mSE) based T2 imaging in separating short T2 signals, which have been suggested to originate from myelin water (Mackay et al., 1994; Whittall et al., 1997), from long T2 signals. Recently, multi-echo gradient echo (mGRE) based T2⁎ imaging has been introduced as an alternative approach (Du et al., 2007; Hwang et al., 2011, 2010; Lenz et al., 2012). The mGRE MWI has several technical advantages such as large volume coverage, fast scan time, insensitivity to B1 inhomogeneity, and a low specific absorption rate over the mSE MWI. Moreover, mGRE can simultaneously generate various susceptibility related contrasts including T2⁎, phase and susceptibility maps (Duyn et al., 2007; Haacke et al., 2004; Oh et al., 2013b; Shmueli et al., 2009), which may provide additional information about myelin integrity (Lee et al., 2012; Li et al., 2009; Liu et al., 2011a).

⁎ Corresponding author. ⁎⁎ Corresponding author. Fax: +82 2 313 2879. E-mail addresses: [email protected] (J. Lee), [email protected] (D.-H. Kim).

http://dx.doi.org/10.1016/j.neuroimage.2015.03.081 1053-8119/© 2015 Elsevier Inc. All rights reserved.

In previous studies (Du et al., 2007; Hwang et al., 2010), mGRE MWI estimated the myelin water fraction (MWF) by fitting a signal decay model to mGRE data. The model had three components with different T2⁎ values. These three components have been suggested to originate from different water compartments (myelin water, axonal water, and extracellular water) in white matter. Among them, the myelin water has been associated with the shortest T2⁎ component due to the restricted mobility by myelin sheath. The other two long T2⁎ components have been suggested to originate from axonal and extracellular space waters (Du et al., 2007; Hwang et al., 2010; Lancaster et al., 2003). Recently, structured residual errors have been observed when the three-pool model was fitted to mGRE data from a region where the fiber orientation is perpendicular to B0 (van Gelderen et al., 2012). These fitting errors have been suggested to originate from the frequency offsets of the three water pools (van Gelderen et al., 2012). This observation has been further consolidated by other studies (Chen et al., 2013; Kim et al., 2014a; Sati et al., 2013; Sukstanskii and Yablonskiy, 2014; Wharton and Bowtell, 2012), demonstrating that the frequency offset and T2⁎ values of individual pools are dependent on the relative orientation of white matter fibers to the B0 field. In particular, the frequency offset between axonal water and myelin water becomes large when fiber orientation is perpendicular to the B0 field. These findings, in addition to the structured residual errors observed in the three-pool model, suggest that the frequency offsets need to be incorporated in the model in order to improve the accuracy of MWF estimation. When adding the frequency offset terms, one potential challenge is to generate reliable frequency estimation that

Y. Nam et al. / NeuroImage 116 (2015) 214–221

only includes frequency shifts from microstructural contribution (e.g. from myelin magnetic susceptibility) while excluding those from macroscopic (non-local) contribution (e.g. from air/tissue magnitude susceptibility difference). Several methods have been proposed for this purpose (Li and Leigh, 2001; Liu et al., 2011b; Neelavalli et al., 2009; Schweser et al., 2011c), but no solution completely separates the two. In this study, we introduce a new complex-valued three-pool model and a new data fitting approach that does not require the removal of the macroscopic frequency shifts. Using this method, we demonstrate that the new complex model provides a more stable MWF estimation than a magnitude model or a magnitude fitting of a complex model by showing reduced sensitivity to the number of echoes used for the fitting.

215

and fitted to the magnitude data. Similarly, here the three-pool complex model was constructed as follows:
Sðt Þ ¼jAmy e

−



1 þi2πΔ f my T2  my

Data acquisition

−

1 þi2πΔ f my−ex T2  my

¼ jAmy e

Signal models Three different models were compared to explain the signal decay characteristics in white matter. Model 1: three-pool magnitude model fitted to magnitude data The first model was the three-pool model suggested by Du et al. (2007) and Hwang et al. (2010). This model consisted of a myelin water pool (my), an axonal water pool (ax) and an extracellular water pool (ex) that had three different T2⁎ values. The model was formulated as follows:


 1 T2  my

t



1 þi2πΔ f ax T2  ax

t




−

þ Aax e


þ Aex e 

1 þi2πΔ f ax−ex T2  ax

t

−



1 þi2πΔ f ex T2  ex




−

þ Aex e

t



1 T2  ex

t

j

j

ð2Þ where Δfmy, Δfax and Δfex represent the frequency offsets of the three water pools. To reduce one fitting parameter, the relative frequency offsets (Δf my − ex and Δf ax − ex) with respect to the extracellular water pool were fitted in this model. The model was fitted to the magnitude signal of the mGRE data (van Gelderen et al., 2012).




Data from eleven volunteers (mean age = 24 ± 3 years) were collected. All scans were performed using a 3 Tesla clinical scanner (Tim Trio, Siemens, Erlangen, Germany) under approval of the local institutional review board. A 12-channel phased-array head coil was used for data reception. For MWI, whole brain 3D mGRE data were acquired using the following parameters: voxel size = 2 × 2 × 2 mm3, matrix size = 128 × 128 × 72, flip angle = 30°, TR = 120 ms, number of echoes = 32, TE = 2.1 ms to 61.93 ms with 1.93 ms echo spacing, unipolar readout, bandwidth per pixel = 1502 Hz, and data acquisition time = 18:43 min. To determine the fiber orientation, diffusion tensor imaging (DTI) was performed with the following parameters: voxel size = 2 × 2 × 2 mm 3, matrix size = 128 × 128 × 64, TR = 9700 ms, TE = 92 ms, diffusion directions = 64, b-value = 600 s/mm2, and data acquisition time = 5 min 30 s. The total scan time including localization and shimming was approximately 25 min. For data reconstruction, a Tukey window (parameter = 1/3) was applied in k-space data to reduce ringing artifacts. Prior to multi-channel signal combination, the global phase offsets of individual channels were corrected (Hammond et al., 2008). The acquired data were reconstructed and analyzed using MATLAB (MathWorks Inc., Natick, MA, USA).

−

þ Aax e 




−

Model 3: three-pool complex model fitted to complex data This new three-pool complex model included background phase term as follows:

Methods

Sðt Þ ¼ Amy e




t


t

−

þ Aax e

1 T2  ax




t

−

þ Aex e

1 T2  ex

 t

ð1Þ

where Amy, Aax and Aex represent the amplitude of the three water pools, and T2⁎ my, T 2⁎ ax , and T 2⁎ ex represent T2⁎ of the three water pools. The model was fitted to the magnitude signal of the mGRE data. Model 2: three-pool complex model fitted to magnitude data In van Gelderen et al.'s study (van Gelderen et al., 2012), the three-pool model was modified to include the frequency offsets

−

Sðt Þ ¼



1 þi2πΔ f my T2  my

Amy e


t

−

þ Aax e

 1 þi2πΔ f ax T2  ax

t


−

þ Aex e

1 þi2πΔ f ex T2  ex

! t

−i2πΔ f bg t−iK0

e

ð3Þ where Δfbg is a background frequency offset term that originates from the macroscopic (nonlocal) field inhomogeneity, and ø0 is an initial phase term that comes from B+ 1 phase offset (Kim et al., 2014b; Schweser et al., 2011a). In previous studies (Sati et al., 2013; Schweser et al., 2011b; Wharton and Bowtell, 2012), the background phase terms (i.e. Δfbg and ø0) were removed before three-pool fitting by background filtering techniques (Li and Leigh, 2001; Liu et al., 2011b; Neelavalli et al., 2009; Schweser et al., 2011c). In our new model, the background frequency offset term (Δfbg) was incorporated into each pool. Then a new frequency offset that included the background frequency offset term was estimated in the complex fitting (e.g. Δfmy + bg = Δfmy + Δfbg) in order to avoid errors in the background filtering process. The new model is written as follows:



 − T 1 þi2πðΔ f bg þΔ f ax Þ t ðΔ f bg þΔ f my Þ t ax 2 Sðt Þ ¼ Amy e þ Aax e
! − T 1 þi2πðΔ f bg þΔ f ex Þ t −iK 2 ex e 0 þAex e

 −

¼

1 þi2π T2  my

−

Amy e

−

þAex e




1 þi2πΔ f myþbg T2  my

1 þi2πΔ f exþbg T2  ex

t

þ Aax e ! t

e

−

1 þi2πΔ f axþbg T2  ax

t

ð4Þ

−iK0

where Δfmy + bg, Δfax + bg and Δfex + bg represent a sum of the background frequency offset and the frequency offset of each pool. In this model, Δfmy + bg − Δfex + bg and Δfax + bg − Δfex + bg correspond to Δfmy − ex and Δfax − ex in Model 2, respectively. The fitting parameters of all models are listed in Table 1 with the initial values and parameter search ranges. These parameters were estimated using an iterative non-linear curve-fitting algorithm (lsqnonlin function in MATLAB, TolX = 1e-5, TolFun = 1e-5). The MWF was defined as the ratio of the signal amplitude of the myelin water to that of the total water (i.e. MWF = Amy / (Amy + Aax + Aex)). ROI analysis Two regions of interest (ROIs) were manually drawn with the aid of a DTI fiber orientation map. Optic radiation regions where primary fiber orientation was perpendicular to B0 were chosen as perpendicular ROIs (Figs. 1a and b). Spinal-cortical tract regions where fibers were primarily aligned parallel to B0 were determined for parallel ROIs (Figs. 1f and g).

216

Y. Nam et al. / NeuroImage 116 (2015) 214–221

Table 1 Initial values and search ranges of the parameters for Model 1 (three-pool magnitude model fitted to magnitude data), Model 2 (three-pool complex model fitted to magnitude data) and n o N−1 ∑n¼1 Sn Snþ1 Model 3 (three-pool complex model fitted to complex data). S1 = S (TE1). Δ f bg;init ¼ ∠ : initial Δfbg (N = number of echoes used in fitting). 2πΔTE Models 1, 2, and 3

Initial value Lower bound Upper bound

Amy

Aax

Aex

T⁎2 my (ms)

T⁎2 ax (ms)

T⁎ 2 ex (ms)

0.1 × |S1| 0 2 × |S1|

0.6 × |S1| 0 2 × |S1|

0.3 × |S1| 0 2 × |S1|

10 3 25

64 25 150

48 25 150

Model 2

Initial value Lower bound Upper bound

Model 3

Δfmy − ex (Hz)

Δfax − ex (Hz)

Δfbg + my (Hz)

Δfbg + ax (Hz)

Δfbg + ex (Hz)

Ø0 (rad)

5 −75 75

0 −25 25

Δfbg,init Δfbg,init − 75 Δfbg,init + 75

Δfbg,init Δfbg,init − 25 Δfbg,init + 25

Δfbg,init Δfbg,init − 25 Δfbg,init + 25

∠S1 −π π

Each ROI region included 80–120 voxels (640–960 cm3) in mGRE images. The ROI data were processed as follows: First, magnitude ROI data were generated by averaging the magnitude of complex data over the ROI. For phase ROI data, the angles of a ROI-averaged complex data were taken. Then, the ROI-averaged signals, which were generated from the magnitude and phase ROI data, were fitted to each model. The residual errors (i.e. |data–model|) were calculated in each TE and averaged across all subjects. The residual errors between the models (Models 1, 2 and 3) as well as between the ROIs (i.e. parallel ROI vs. perpendicular ROI) were compared. To evaluate the reliability of each model, the ROI analysis was repeated for the data with a smaller number of early echoes (from the first 12 echoes to first 31 echoes). The resulting values of MWF and Δfmy − ex from each data were averaged across all subjects and compared. For each data, Student's t-test was performed with respect to the values estimated using all 32 echoes. The root mean squared errors (RMSE) were calculated for all models. Voxel-by-voxel analysis Whole brain MWF maps were generated for three different datasets (data with first 16 echoes, first 24 echoes, and all 32 echoes)

using the three models. Compared to the ROI-averaged signal fitting, voxel-by-voxel signal fitting was challenging due to a low SNR and other factors including macroscopic B0 field inhomogeneity. To reduce these effects, two adjacent slices were averaged, and weighted least-squares fitting was performed with the magnitude of each echo  2 n as a weight (i.e. minyb ∑i¼1 wi yi −ybι was performed where i is an ι

echo index, w i is magnitude of each echo, y i is the ith echo data, and ybι is the ith echo fitted value). The fitting was performed for each voxel using the aforementioned procedure. The estimated Δfmy − ex map was compared with the corresponding slice of DTI fiber orientation map to examine the reliability of the estimated Δfmy − ex map from Model 3. Furthermore, a myelin frequency offset map was established from DTI data using a hollow cylinder model suggested by Wharton and Bowtell (2012). Using the model, the frequency shift of myelin water can be written as follows: Δ f my ¼

γB0 2π

χI 2



2

cos θ−

 1 χ þEþ A 3 2

 !!! 1 1 3 r2i r 2 − 2 2 ln o − þ sin θ 3 4 2 ro −ri ri

ð5Þ where γ represents the gyromagnetic ratio and θ represents the angle between fiber bundles and B0. The following parameters were used to

Fig. 1. Population-averaged (n = 11) residual patterns for Model 1 (three-pool magnitude model fitted to magnitude data), Model 2 (three-pool complex model fitted to magnitude data) and Model 3 (three-pool complex model fitted to complex data) when all 32 echoes were used in the fitting. (a, b, f and g) Examples of perpendicular and parallel ROIs in the first echo GRE magnitude images and the corresponding DTI fiber orientation maps. (c, d and e) The residues (mean ± s.d.) of the perpendicular ROI. (h, i and j) The residues (mean ± s.d.) of the parallel ROI.

Y. Nam et al. / NeuroImage 116 (2015) 214–221

calculate the frequency offset of myelin water with a DTI fiber orientation map: isotropic susceptibility (χI) = −100 ppb; anisotropic susceptibility (χA) = −100 ppb; chemical exchange (E) = 20 ppb; gratio (ro/ri where ro: outer radius and ri: inner radius) = 0.8; and B0 = 3 T (Kim et al., 2014a). To evaluate the improvement of the proposed complex model fitting method (Model 3), a previous complex model fitting method, which removes background phase terms before the model fitting, was also performed for the same datasets with the first 16 echoes, and the estimated frequency offset map was compared with that from Model 3. For this purpose, phase data was high-pass filtered to remove the background frequency offset term (Δfbg) and the constant phase term (ø0). The high-pass filtering was performed by dividing the original complex images with the low-pass filtered complex images. The lowpass filtering was carried out using a Gaussian window. Three different window sizes (σ = 2, 4, and 8 mm) were tested to investigate the effects of the window size on the model parameters. Then, the following complex three-pool model (Sati et al., 2013) was fitted to the phase filtered complex data:
−

Sðt Þ ¼ Amy e

þ Aex e



1 þi2πΔ f my T2  my




−

1 þi2πΔ f ex T2  ex


−

t



þ Aax e

t

:

 1 þi2πΔ f ax T2  ax

t

ð6Þ

The initial values and parameter search ranges were identical to Model 3 as listed in Table 1 except the initial values of Δfmy, Δfax and Δfex which were set to zero. After estimating all parameters, a Δfmy map was compared with the Δfmy − ex map from Model 3. Results Fig. 1 shows the ROI locations (Figs. 1a and f), DTI fiber orientation maps (Figs. 1b and g), and residual errors from each model (Figs. 1c–e for the perpendicular ROI and Figs. 1h–j for the parallel ROI) averaged over eleven subjects. As demonstrated by van Gelderen et al. (2012),

217

the perpendicular ROI showed larger residues in Model 1 confirming systematic errors in the model. When Model 2 or 3 was fitted, the residues were substantially reduced. This result supports that the frequency offsets need to be included to explain the signal decay pattern of GRE data (Sati et al., 2013; van Gelderen et al., 2012; Wharton and Bowtell, 2012). In the parallel ROI, the residual pattern did not show a noticeable difference between the models (Figs. 1h, i, and j), indicating a smaller contribution of the frequency offsets in the signal when fibers were parallel to B0 (Sati et al., 2013; van Gelderen et al., 2012; Wharton and Bowtell, 2012). In Fig. 2, the MWFs estimated from six different numbers of echoes (from 12 echoes to 32 echoes) are plotted for the three models in both ROIs. When Model 1 was fitted to the perpendicular ROI (Fig. 2a), the MWFs showed large variation (Fig. 1a; minimum MWF is 14.2% with 32 echoes and maximum MWF is 19.1% with 20 echoes; asterisks are when p b 0.05 compared with MWF in the 32 echoes). On the other hand, the variation was substantially smaller when data were fitted with Model 2 (Fig. 2b; minimum MWF is 12.7% with 12 echoes and maximum MWF is 15.1% with 28 echoes) or Model 3 (Fig. 2c; the minimum MWF is 12.2% with 24 echoes and the maximum MWF is 14.4% with 16 echoes). These results demonstrate that when the frequency offsets were included in the model (Models 2 and 3), the MWF estimation became more stable. For the parallel ROI, the variations in the MWF were relatively small in all models (Figs. 2d, e, and f; no statistically significant difference). Similar results are quantitatively summarized in Table 2 where the MWFs, frequency offsets, and root mean squared errors (RMSE) of the two ROIs are averaged over all subjects for 16, 24 and 32 echoes. When Model 2 and Model 3 were compared, Model 3 showed smaller standard deviations (across the subjects) in Δfmy − ex estimation for the perpendicular ROI. This observation indicates a more stable estimation of Δfmy − ex in perpendicular fibers in Model 3. Fig. 3 shows the results from the voxel-by-voxel analysis. The magnitude images of last echoes used in the fitting (Fig. 3a), the MWF maps (Fig. 3b), the Δfmy − ex maps (Fig. 3c) are illustrated for the three different numbers of the data (first 16, first 24 and all 32 echoes)

Fig. 2. MWF estimated using various numbers of echoes (from 12 to 32 echoes). The means and standard deviations were calculated from the MWFs of the eleven subjects. (a) The estimated MWFs of the perpendicular ROI when Model 1 was fitted. (b) The estimated MWFs of the perpendicular ROI when Model 2 was fitted. (c) The estimated MWFs of the perpendicular ROI when Model 3 was fitted. (d) The estimated MWFs of the parallel ROI when Model 1 was fitted. (e) The estimated MWFs of the parallel ROI when Model 2 was fitted. (f) The estimated MWFs of the parallel ROI when Model 3 was fitted. Red asterisks indicate p b 0.05 in t-test with respect to the results estimated from all 32 echoes.

218

Y. Nam et al. / NeuroImage 116 (2015) 214–221

Table 2 ROI fitting results for Model 1 (three-pool magnitude model fitted to magnitude data), Model 2 (three-pool complex model fitted to magnitude data) and Model 3 (three-pool complex model fitted to complex data). The average (±standard deviation) MWF, Δfmy − ex, and root mean squared errors (RMSE) were estimated for the eleven subjects. Model 1

Model 2

Model 3

# of echoes used in fitting

MWF (%)

RMSE ×103

MWF (%)

Δfmy − ex (Hz)

RMSE ×103

MWF (%)

Δfmy − ex (Hz)

RMSE ×103

Perpendicular ROI 16 24 32

18.2 ± 3.1⁎ 19.1 ± 4.4⁎ 14.2 ± 2.2

1.3 ± 0.2 1.3 ± 0.3 2.0 ± 0.9

14.5 ± 2.9 14.6 ± 2.5 14.5 ± 3.0

8.1 ± 1.9 8.3 ± 3.7 8.6 ± 8.1

1.2 ± 0.2 1.1 ± 0.2 1.1 ± 0.2

14.4 ± 3.1 12.2 ± 1.4 13.5 ± 2.8

11.1 ± 1.4 13.3 ± 2.8 12.7 ± 1.7

1.1 ± 0.1 1.2 ± 0.2 1.3 ± 0.1

7.1 ± 1.4 5.7 ± 1.8 5.7 ± 2.0

1.1 ± 0.4 1.0 ± 0.3 1.1 ± 0.6

9.2 ± 1.7 8.5 ± 1.5 8.9 ± 1.7

2.4 ± 1.1 1.4 ± 1.7 1.3 ± 2.5

1.0 ± 0.3 0.9 ± 0.2 1.1 ± 0.6

8.5 ± 1.8 7.3 ± 1.9 7.3 ± 2.0

2.5 ± 0.9 2.4 ± 1.8 3.8 ± 3.1

1.1 ± 0.4 1.1 ± 0.4 1.3 ± 0.6

Parallel ROI 16 24 32 ⁎ Indicates p b 0.05.

used in the fitting. When the MWF maps were compared, Models 2 and 3 showed relatively consistent MWF maps across the different numbers of echoes. On the other hand, the MWF maps in Model 1 showed larger variations (for example, optic radiation regions; white arrows in Fig. 3b). When Model 2 and Model 3 were compared, Δfmy − ex maps from Model 3 (right column in Fig. 3c) showed less noisy and more consistent results over the different numbers of echoes than Model 2 (left column in Fig. 3c), agreeing with the ROI analysis results (Table 2). The Δfmy − ex maps from Model 3 also showed better agreement with a Δfmy map (Fig. 3e) generated by the hollow cylinder model (Wharton and Bowtell, 2012). Hence, Model 3 generated more stable maps (both MWF and Δfmy − ex) for the different numbers of echoes used in the data fitting. When the frontal lobe area was examined (red arrows in Fig. 3b), the MWF maps and the Δfmy − ex maps generated from the three different numbers of echoes showed large variations. This is probably due to the large macroscopic field inhomogeneity around the air–tissue interface (Yablonskiy, 1998; Yablonskiy and Haacke, 1994). The effects become more pronounced at late echoes. Hence, using a smaller number of early echoes (e.g. first 16 echoes) may improve the MWF estimation in the areas of large field inhomogeneity. The estimated parameter maps from Model 3 using the first 16 echoes are shown in Fig. 4. Additionally, a MWF map, a Δfmy − ex map, a color-coded DTI FA map, and a Δfmy map from the hollow cylinder model are included (Figs. 4d, e, i and j). The individual frequency offset maps (Figs. 4f, g, and h) included nonlocal background frequency shifts. These shifts, however, looked the same in all three frequency offset

maps and were removed in the difference map (Fig. 4i). The resulting Δfmy − ex map clearly showed that the characteristics of the myelin water pool, revealing large positive frequency shifts when the fibers were perpendicular to B0 (red and green regions in the DTI map) and slight negative to zero frequency shifts when the fibers were parallel to B0 (blue regions in the DTI map). The Δfmy − ex map also showed a similar pattern to the myelin water frequency map generated by the hollow cylinder model (Fig. 4j). When the axonal water frequency shift map (Δfax + bg, Fig. 4g) was examined, areas in which fibers were largely perpendicular to B0 yield slightly negative frequency shifts relative to neighboring areas agreeing with previous reports (Sati et al., 2013). Lastly, the extracellular water frequency shift map (Fig. 4h) showed almost no structure (i.e. Δfex ~ 0) as suggested in the previous study (Sati et al., 2013), supporting our approach of using Δfmy − ex as a substitute for Δfmy. The parameters estimated by Model 3 using the first 16 echoes in the two ROIs were summarized in Table S1. The results are in similar ranges to previous studies (Sati et al., 2013; Raven et al., 2014). Fig. 5 shows a 3D MWF map and a 3D Δfmy − ex map generated using the first 16 echo data in Model 3. Color-coded DTI FA maps of the corresponding slices are displayed as a reference. The MWF map showed regions with high MWF (e.g. optic radiation and corpus callosum) and the Δfmy − ex map highlighted regions that have fibers perpendicular to the B0 field. Fig. 6 demonstrates the advantage of our proposed complex fitting method. The method used an unfiltered phase image (Fig. 6a) and generated a myelin water frequency offset (Δfmy − ex) map without

Fig. 3. Estimated MWF, Δfmy − ex maps from 16, 24, and 32 echoes used in fitting and Δfmy map from DTI data and hollow cylinder model for Model 1, Model 2 and Model 3. (a) The magnitude images of the corresponding echo times. (b) The estimated MWF maps from different models (Models 1, 2 and 3). (c) The estimated Δfmy − ex maps and (d) The DTI FA map. (e) The calculated Δfmy map from DTI data using hollow cylinder model (HCM).

Y. Nam et al. / NeuroImage 116 (2015) 214–221

219

Fig. 4. Parameter maps from Model 3 (three-pool complex model fitted to complex data) fitted with first 16 echoes. The estimated fraction maps of (a) myelin water (Amy), (b) axonal water (Aax) and (c) mixed water (Aex) pools. (d) The estimated MWF map. (e) The corresponding DTI FA map. The estimated frequency offset (including background phase offset, Δfbg) maps of (f) myelin water (Δfbg + my), (g) axonal water (Δfbg + ax) and (h) extracellular water (Δfbg + ex) pools. (i) The calculated Δfmy − ex map and (j) the calculated Δfmy map from DTI data using hollow cylinder model (HCM).

noticeable artifacts even at the region with large B0 field inhomogeneity (e.g. dashed red circle in Fig. 6a). When the background frequency terms were removed before the model fitting, however, the myelin water fraction map and the frequency offset map became dependent on the filter size (Figs. 6f, g, h, j, k and l). Discussion and conclusion In this study, we introduced a new complex model and a new data fitting approach that explain the signal decay characteristics of GRE data. This method does not require any background field estimation and removal procedure and, as shown, provides more robust estimation of the model parameters as compared to conventional models.

One of important contribution of this new model and fitting approach is that any background field estimation and removal steps are not required. In previous studies, these steps were performed before complex model fitting (Sati et al., 2013; Schweser et al., 2011b; Wharton and Bowtell, 2012, 2013). As shown in Fig. 6, imperfect filtering leads to incorrect estimation of model parameters. Although a more sophisticated processing method (Li and Leigh, 2001; Liu et al., 2011b; Neelavalli et al., 2009; Sati et al., 2013; Schweser et al., 2011c) may improve the background field estimation, one cannot completely remove the background field. Hence, residues may potentially lead to significant errors during the data fitting process. One potential limitation of our approach is that the absolute frequency offsets of three pools cannot be estimated. However, since

Fig. 5. Seven representative slices of a healthy volunteer from (a) the estimated MWF map, (b) the estimated Δfmy − ex map and (c) DTI FA map. Model 3 (three-pool complex model fitted to complex data) was fitted to the first 16 echoes.

220

Y. Nam et al. / NeuroImage 116 (2015) 214–221

Fig. 6. Comparison of the estimated frequency offset difference between myelin water and extracellular water pools (Δfmy − ex) from the unfiltered phase data and the high-pass filtered phase data. The first 16 echoes were used for parameter estimation. (a) The 16th echo phase image. (b–d) The 16th echo filtered phase images for different window sizes of filter (σ = 2, 4, and 8 mm). (e) The estimated MWF map using the proposed complex fitting method (Model 3). (g–h) The estimated MWF maps from the high-pass filtered phase data for different window sizes. (i) The estimated Δfmy − ex map using the proposed complex fitting method (Model 3). (j–l) The estimated Δfmy maps from the high-pass filtered phase data for different window sizes.

the frequency offset value of the extracellular water pool has been demonstrated to be close to zero (Sati et al., 2013), Δfmy and Δfax can be approximated by Δfmy − ex and Δfax − ex respectively in our model. This is further supported by Δfex + bg map in Fig. 4 where the map shows almost no white matter structure. Shortening data acquisition duration has a practical value since it determines the minimum TR and, in turn, total scan time. Hence, using a smaller number of echoes helps to reduce the total scan time or increase spatial coverage. In addition, the use of late echoes in mGRE hampers the fitting results in the areas of large background field inhomogeneity because it progressively affects the signal over TE. Our results suggest that mGRE MWI may not need all 32 echoes. The use of 16 echoes (last TE = 31.1 ms) which reduces the effects of the macroscopic B0 field inhomogeneity and increases the spatial coverage appears to be a reasonable choice. Multiple methods have been proposed to compensate for the macroscopic B0 field inhomogeneity in mGRE (Han et al., 2015; Hernando et al., 2012; Nam et al., 2012; Oh et al., 2014; Yablonskiy et al., 2013). However, many of them require a longer minimum TE and may lose a significant portion of fast decaying myelin water signal. Hence, it may be challenging to apply them directly. In this work, the effects of T1 were not considered in the MWF estimation. It has been suggested that a short TR may overestimate MWF in mGRE (Du et al., 2007; Li et al., 2013; Nam et al., 2014; Oh et al., 2013a). In addition, the short TR may have introduced magnetization transfer effects in the data. Further research is necessary to investigate the effects of TR on MWF and to develop an approach to compensate for the overestimation. In perpendicular ROI, Model 3 showed more stable estimation of the frequency offsets than Model 2. In parallel ROI, however, Model 2 show

slightly smaller RMSE than Model 3. This is because the magnitude model has better estimation of the signal than the complex model when phase contrast does not exist (Lee et al., 2007). Since the difference in RMSE is small in parallel ROI and the improvement of Model 3 in perpendicular ROI are large, the use of Model 3 is overall beneficial. When compared with mSE MWI results (Laule et al., 2004; Oh et al., 2006; Prasloski et al., 2012; Whittall et al., 1997), our results overall showed similar MWF distribution showing higher MWF in corpus callosum, optic radiation, longitudinal fasciculus and internal capsule (see Fig. 6) although the MWF in internal capsule region is not as high as that in mSE MWI (Laule et al., 2004; Oh et al., 2006; Prasloski et al., 2012; Whittall et al., 1997). This may be related to a recent mSE MWI study suggesting an overestimated MWF in this area (Russell-Schulz et al., 2013). Since the new model includes phase signal, it may become more sensitive to physiological noise sources such as respiration that have more pronounced effects on phase (Noll and Schneider, 1994). Therefore, a proper correction for these factors may be necessary for an accurate estimation of the model parameters (Hu and Kim, 1994). Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.neuroimage.2015.03.081.

Acknowledgment This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (NRF2012R1A2A2A01009903).

Y. Nam et al. / NeuroImage 116 (2015) 214–221

References Chen, W.C., Foxley, S., Miller, K.L., 2013. Detecting microstructural properties of white matter based on compartmentalization of magnetic susceptibility. Neuroimage 70, 1–9. Du, Y.P., Chu, R., Hwang, D., Brown, M.S., Kleinschmidt‐DeMasters, B.K., Singel, D., Simon, J.H., 2007. Fast multislice mapping of the myelin water fraction using multicompartment analysis of T2* decay at 3 T: a preliminary postmortem study. Magn. Reson. Med. 58, 865–870. Duyn, J.H., van Gelderen, P., Li, T.-Q., de Zwart, J.A., Koretsky, A.P., Fukunaga, M., 2007. High-field MRI of brain cortical substructure based on signal phase. Proc. Natl. Acad. Sci. U. S. A. 104, 11796–11801. Haacke, E.M., Xu, Y., Cheng, Y.C.N., Reichenbach, J.R., 2004. Susceptibility weighted imaging (SWI). Magn. Reson. Med. 52, 612–618. Hammond, K.E., Lupo, J.M., Xu, D., Metcalf, M., Kelley, D.A., Pelletier, D., Chang, S.M., Mukherjee, P., Vigneron, D.B., Nelson, S.J., 2008. Development of a robust method for generating 7.0 T multichannel phase images of the brain with application to normal volunteers and patients with neurological diseases. Neuroimage 39, 1682–1692. Han, D., Nam, Y., Gho, S.-M., Kim, D.-H., 2015. Volumetric R2* mapping using z-shim multi-echo gradient echo imaging. Magn. Reson. Med. 73, 1164–1170. Hernando, D., Vigen, K.K., Shimakawa, A., Reeder, S.B., 2012. R2* mapping in the presence of macroscopic B0 field variations. Magn. Reson. Med. 68, 830–840. Hu, X., Kim, S.G., 1994. Reduction of signal fluctuation in functional MRI using navigator echoes. Magn. Reson. Med. 31, 495–503. Hwang, D., Kim, D.-H., Du, Y.P., 2010. In vivo multi-slice mapping of myelin water content using T2* decay. Neuroimage 52, 198–204. Hwang, D., Chung, H., Nam, Y., Du, Y.P., Jang, U., 2011. Robust mapping of the myelin water fraction in the presence of noise: synergic combination of anisotropic diffusion filter and spatially regularized nonnegative least squares algorithm. J. Magn. Reson. Imaging 34, 189–195. Kim, D., Lee, H.M., Oh, S.-H., Lee, J., 2014a. Probing signal phase in direct visualization of short transverse relaxation time component (ViSTa). Magn. Reson. Med. http://dx. doi.org/10.1002/mrm.25416. Kim, D.H., Choi, N., Gho, S.M., Shin, J., Liu, C., 2014b. Simultaneous imaging of in vivo conductivity and susceptibility. Magn. Reson. Med. 71, 1144–1150. Lancaster, J.L., Andrews, T., Hardies, L.J., Dodd, S., Fox, P.T., 2003. Three‐pool model of white matter. J. Magn. Reson. Imaging 17, 1–10. Laule, C., Vavasour, I., Moore, G., Oger, J., Li, D., Paty, D., MacKay, A., 2004. Water content and myelin water fraction in multiple sclerosis. J. Neurol. 251, 284–293. Lee, J., Shahram, M., Schwartzman, A., Pauly, J.M., 2007. Complex data analysis in high‐resolution SSFP fMRI. Magn. Reson. Med. 57, 905–917. Lee, J., Shmueli, K., Kang, B.-T., Yao, B., Fukunaga, M., van Gelderen, P., Palumbo, S., Bosetti, F., Silva, A.C., Duyn, J.H., 2012. The contribution of myelin to magnetic susceptibilityweighted contrasts in high-field MRI of the brain. Neuroimage 59, 3967–3975. Lenz, C., Klarhöfer, M., Scheffler, K., 2012. Feasibility of in vivo myelin water imaging using 3D multigradient-echo pulse sequences. Magn. Reson. Med. 68, 523–528. Li, L., Leigh, J.S., 2001. High-precision mapping of the magnetic field utilizing the harmonic function mean value property. J. Magn. Reson. 148, 442–448. Li, T.-Q., Yao, B., van Gelderen, P., Merkle, H., Dodd, S., Talagala, L., Koretsky, A.P., Duyn, J., 2009. Characterization of T2* heterogeneity in human brain white matter. Magn. Reson. Med. 62, 1652–1657. Li, W., Han, H., Guidon, A., Liu, C., 2013. Dependence of gradient echo phase contrast on the differential signal decay in subcellular compartments. Proc. Int. Soc. Magn. Reson. Med. p161. Liu, C., Li, W., Johnson, G.A., Wu, B., 2011a. High-field (9.4 T) MRI of brain dysmyelination by quantitative mapping of magnetic susceptibility. Neuroimage 56, 930–938. Liu, T., Khalidov, I., de Rochefort, L., Spincemaille, P., Liu, J., Tsiouris, A.J., Wang, Y., 2011b. A novel background field removal method for MRI using projection onto dipole fields (PDF). NMR Biomed. 24, 1129–1136. Mackay, A., Whittall, K., Adler, J., Li, D., Paty, D., Graeb, D., 1994. In vivo visualization of myelin water in brain by magnetic resonance. Magn. Reson. Med. 31, 673–677. Nam, Y., Han, D., Kim, D.-H., 2012. Single-scan R2⁎ measurement with macroscopic field inhomogeneity correction. Neuroimage 63, 1790–1799.

221

Nam, Y., Lee, J., Kim, D.-H., 2014. Dependence of short T2* fraction on TR in gradient echo imaging. 3rd International Workshop on MRI Phase Contrast & Quantitative Susceptibility Mapping. Neelavalli, J., Cheng, Y.C.N., Jiang, J., Haacke, E.M., 2009. Removing background phase variations in susceptibility‐weighted imaging using a fast, forward‐field calculation. J. Magn. Reson. Imaging 29, 937–948. Noll, D.C., Schneider, W., 1994. Theory, simulation, and compensation of physiological motion artifacts in functional MRI. Image Processing, 1994. Proceedings. ICIP-94., IEEE International Conference. IEEE, pp. 40–44. Oh, J., Han, E.T., Pelletier, D., Nelson, S.J., 2006. Measurement of in vivo multi-component T2 relaxation times for brain tissue using multi-slice T2 prep at 1.5 and 3 T. Magn. Reson. Imaging 24, 33–43. Oh, S.-H., Bilello, M., Schindler, M., Markowitz, C.E., Detre, J.A., Lee, J., 2013a. Direct visualization of short transverse relaxation time component (ViSTa). Neuroimage 83, 485–492. Oh, S.-H., Kim, Y.-B., Cho, Z.-H., Lee, J., 2013b. Origin of B0 orientation dependent R2* (= 1/T2*) in white matter. Neuroimage 73, 71–79. Oh, S.S., Oh, S.H., Nam, Y., Han, D., Stafford, R.B., Hwang, J., Kim, D.H., Park, H., Lee, J., 2014. Improved susceptibility weighted imaging method using multi-echo acquisition. Magn. Reson. Med. 72, 452–458. Prasloski, T., Rauscher, A., MacKay, A.L., Hodgson, M., Vavasour, I.M., Laule, C., Mädler, B., 2012. Rapid whole cerebrum myelin water imaging using a 3D GRASE sequence. Neuroimage 63, 533–539. Raven, E.P., van Gelderen, P., de Zwart, J., VanMeter, J., Duyn, J.H., 2014. In vivo assessment of age-related white matter differences using T2* relaxation. 3rd International Workshop on MRI Phase Contrast & Quantitative Susceptibility Mapping. Russell-Schulz, B., Laule, C., Li, D.K.B., MacKay, A.L., 2013. What causes the hyperintense T2-weighting and increased short T2 signal in the corticospinal tract? Magn. Reson. Imaging 31, 329–335. Sati, P., van Gelderen, P., Silva, A.C., Reich, D.S., Merkle, H., de Zwart, J.A., Duyn, J.H., 2013. Micro-compartment specific T2* relaxation in the brain. Neuroimage 77, 268–278. Schweser, F., Atterbury, M., Deistung, A., Lehr, B., Sommer, K., Reichenbach, J., 2011a. Harmonic phase subtraction methods are prone to B1 background components. Proc. Int. Soc. Magn. Reson. Med. 2657. Schweser, F., Deistung, A., Gullmar, D., Atterbury, M., Lehr, B.W., Sommer, K., Reichenbach, J.R., 2011b. Non-linear evolution of GRE phase as a means to investigate tissue microstructure. Proc. Int. Soc. Magn. Reson. Med. p4527. Schweser, F., Deistung, A., Lehr, B.W., Reichenbach, J.R., 2011c. Quantitative imaging of intrinsic magnetic tissue properties using MRI signal phase: an approach to in vivo brain iron metabolism? Neuroimage 54, 2789–2807. Shmueli, K., de Zwart, J.A., van Gelderen, P., Li, T.Q., Dodd, S.J., Duyn, J.H., 2009. Magnetic susceptibility mapping of brain tissue in vivo using MRI phase data. Magn. Reson. Med. 62, 1510–1522. Sukstanskii, A.L., Yablonskiy, D.A., 2014. On the role of neuronal magnetic susceptibility and structure symmetry on gradient echo MR signal formation. Magn. Reson. Med. 71, 345–353. van Gelderen, P., de Zwart, J., Lee, J., Sati, P., Reich, D., Duyn, J., 2012. Nonexponential T2 decay in white matter. Magn. Reson. Med. 67, 110. Wharton, S., Bowtell, R., 2012. Fiber orientation-dependent white matter contrast in gradient echo MRI. Proc. Natl. Acad. Sci. U. S. A. 109, 18559–18564. Wharton, S., Bowtell, R., 2013. Gradient echo based fiber orientation mapping using R2* and frequency difference measurements. Neuroimage 83, 1011–1023. Whittall, K.P., Mackay, A.L., Graeb, D.A., Nugent, R.A., Li, D.K.B., Paty, D.W., 1997. In vivo measurement of T2 distributions and water contents in normal human brain. Magn. Reson. Med. 37, 34–43. Yablonskiy, D.A., 1998. Quantitation of intrinsic magnetic susceptibility‐related effects in a tissue matrix. Phantom study. Magn. Reson. Med. 39, 417–428. Yablonskiy, D.A., Haacke, E.M., 1994. Theory of NMR signal behavior in magnetically inhomogeneous tissues: the static dephasing regime. Magn. Reson. Med. 32, 749–763. Yablonskiy, D.A., Sukstanskii, A.L., Luo, J., Wang, X., 2013. Voxel spread function method for correction of magnetic field inhomogeneity effects in quantitative gradient‐ echo‐based MRI. Magn. Reson. Med. 70, 1283–1292.