1. No category

YNIMG-12123; No. of pages: 8; 4C: 3, 4, 5, 6, 7
NeuroImage xxx (2015) xxx–xxx
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 xxxx
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 Published by Elsevier Inc.
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).
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
http://dx.doi.org/10.1016/j.neuroimage.2015.03.081
1053-8119/© 2015 Published by Elsevier Inc.
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081
2
Y. Nam et al. / NeuroImage xxx (2015) xxx–xxx
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.
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).
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081
Y. Nam et al. / NeuroImage xxx (2015) xxx–xxx
3
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.
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081
4
Y. Nam et al. / NeuroImage xxx (2015) xxx–xxx
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),
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.
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081
Y. Nam et al. / NeuroImage xxx (2015) xxx–xxx
5
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).
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081
6
Y. Nam et al. / NeuroImage xxx (2015) xxx–xxx
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.
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081
Y. Nam et al. / NeuroImage xxx (2015) xxx–xxx
7
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).
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081
8
Y. Nam et al. / NeuroImage xxx (2015) xxx–xxx
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.
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.
Please cite this article as: Nam, Y., et al., Improved estimation of myelin water fraction using complex model fitting, NeuroImage (2015), http://
dx.doi.org/10.1016/j.neuroimage.2015.03.081