Assessment and correction of macroscopic field variations in 2D spoiled gradient-echo sequences.

PURPOSE: To model and correct the dephasing effects in the gradient-echo signal for arbitrary RF excitation pulses with large flip angles in the presence of macroscopic field variations.
METHODS: The dephasing of the spoiled 2D gradient-echo signal was modeled using a numerical solution of the Bloch equations to calculate the magnitude and phase of the transverse magnetization across the slice profile. Additionally, regional variations of the transmit RF field and slice profile scaling due to macroscopic field gradients were included. Simulations, phantom, and in vivo measurements at 3 T were conducted for R 2 ∗ and myelin water fraction (MWF) mapping.
RESULTS: The influence of macroscopic field gradients on R 2 ∗ and myelin water fraction estimation can be substantially reduced by applying the proposed model. Moreover, it was shown that the dephasing over time for flip angles of 60° or greater also depends on the polarity of the slice-selection gradient because of phase variation along the slice profile.
CONCLUSION: Substantial improvements in R 2 ∗ accuracy and myelin water fraction mapping coverage can be achieved using the proposed model if higher flip angles are required. In this context, we demonstrated that the phase along the slice profile and the polarity of the slice-selection gradient are essential for proper modeling of the gradient-echo signal in the presence of macroscopic field variations.

INTRODUCTION

The sensitivity of gradient‐echo (GRE) imaging to variations in magnetic susceptibility has led to a widespread range of applications in MRI. Approaches focusing on the signal decay are methods such as FMRI, in which the susceptibility difference between deoxyhemoglobin and oxyhemoglobin is measured,1, 2 perfusion MRI with gadolinium‐based contrast agent,3 mapping by acquiring multi‐gradient‐echo (mGRE), or the determination of the myelin water fraction (MWF) from analyzing the multi‐exponential decay of GRE signal.4 Methods exploiting the phase evolution of the signal are SWI5 or QSM.6 Quantitative susceptibility mapping and contrast in the brain have been used to study iron content in inflammatory and neurodegenerative diseases.7, 8, 9, 10, 11

A challenge with quantifying mGRE data are the macroscopic field variations that arise, for example, from air/tissue boundaries, which lead to a faster signal decay and consequently mask tissue‐relevant mesoscopic‐scale and microstructure‐scale effects.12 Approaches for correcting the influence of macroscopic field variations can be grouped into methods based on sequence design or postprocessing of conventional mGRE data. Methods requiring sequence adaptations aim to compensate intravoxel dephasing by varying the slice‐rephasing gradient and/or by applying compensation gradients in the slice‐selective direction.13, 14, 15, 16, 17 Beyond sequence programming, a powerful approach is to simply increase spatial resolution,18 which is not always feasible in a clinical setup with scan‐time constraints or limited SNR.

The focus of the present work is the correction of the influence of macroscopic field variations by postprocessing of mGRE data. Given that for 2D‐GRE acquisitions the slice thickness is usually much larger than the in‐plane resolution, signal dephasing is largely driven by macroscopic field variations along the slice‐selective direction z.13 Assuming an ideal slice profile and a constant field gradient G as an approximation of the macroscopic field variation along z, the signal over time is sinc‐weighted proportional to G and slice thickness.12 Based on this theory, Fernandez‐Seara and Wehrli19 estimated G and iteratively from the measured signal decay, which later was refined by initialization of G using the field map20 and further extended by Yang et al by modeling field variation as a quadratic function, a condition that has become more relevant at ultrahigh‐field MRI.21 Depending on the RF excitation pulse, deviations from the ideal slice profile cause different dephasing in the presence of G. To account for various pulse shapes, Preibisch et al proposed a solution in which the signal modulation due to G and the slice profile is described by the envelope of the RF pulse.22 This model allows to describe the signal decay for flip angles α of less than 60° and a much longer TR than the longitudinal relaxation time T1 to avoid saturation of the slice profile.22 To achieve a smooth signal decay, they have used exponential RF pulses23 and compared them with sinc‐shaped and sinc‐Gauss shaped pulses for mapping.24

Similar to mapping, in MWF mapping, macroscopic field variations need to be accounted for.25 Here, modeling approaches also assume an ideal slice profile.25, 26 In recent work, Lee et al combined z‐shimming with modeling of the magnitude of the slice profile to achieve a better modeling of the signal decay.27

In contrast to the analytical solution in Preibisch et al,22 which is limited by the small flip angle approximation,28 we here propose a numerical model for solving the signal dephasing in the presence of G for an arbitrary excitation pulse and flip angle. Consequently, the model allows to benefit from increased SNR in measurements with interleaved slice acquisition, especially for large flip angles (α > 60°). Extending the model from Hernando et al,29 we also investigate variations of the transmit RF field and the effect of scaling of the slice profile due to superposition of G and the slice‐selection gradient G. We further demonstrate with phantom and in vivo measurements that, depending on the pulse shape for larger flip angles, the polarity of G has to be considered, because through‐slice phase variations can severely affect signal dephasing. With the proposed model, it is possible to substantially improve the quality of maps and MWF maps acquired with arbitrary excitation pulses and flip angles.

METHODS

Theory

In the presence of macroscopic field variations Δω(z), the signal S(t) of a 2D spoiled GRE is proportional to the integral over the complex transverse magnetization weighted with Δω(z) along the slice‐selective direction z. Thus, depending on and Δω(z), additional signal dephasing is observed in contrast to theoretical mono‐exponential signal decay with . If Δω(z) is smooth and slowly varying in space, Δω(z) can be approximated with a linear function in each slice.12 By assuming the origin of z being in the center of the slice, the equation for S(t) reads as follows:where G denotes the field gradient and Δω 0 denotes the field offset. The value of depends on several factors (including ξ, λ, and E 1), discussed in detail subsequently. Depending on the ratio of the TR and the T1, which is included in the exponential term , and the effective flip angle α(z) along the slice, the solution for |M(z)| changes according to the steady‐state equation for spoiled GRE sequences30 as follows:

When TR is much larger than T1, Equation 2 simplifies and |M(z)| is obtained by the sine of α times the equilibrium magnetization S 0:

Here, α(z) is obtained for a certain slice‐selection gradient G and the applied excitation pulse with a certain shape and amplitude. For small flip angles, the slice profile α(z) can be estimated for an RF pulse envelope B 1(t) with the small flip angle approximation.28 However, larger flip angles require solving the Bloch equations for |M(z)| and φ(z).

Extending previous studies, the factors ξ and λ were added to describe 2 effects that affect α(z) and therefore signal dephasing. First, variations of the active transmit field () cause a deviation from the nominal flip angle α, which can change the effective flip angle profile α(z), and therefore requires α to be scaled with ξ, obtained from the normalized B 1 map. Second, G is superimposed with G, leading to either broadening or narrowing of the slice profile described by the factor λ31 as follows:

To investigate the effect of the described parameters on signal dephasing in the presence of G, 4 different models have been studied. Summarizing Equation 1 in a tissue‐specific signal component S tissue(t) (e.g., ) and a component F(t) describing the signal dephasing due to Δω(z), the model can be written as . The 4 models are defined as follows:

The model S 1(t) serves as an uncorrected reference without modeling (z) and Δω(z). Then, for S 2(t), only the magnitude along the slice |M(z)| was considered neglecting φ(z). In S 3(t), φ(z) was included, and in S 4(t) the model was extended by additionally incorporating and λ variations.

Numerical implementation

Signal dephasing due to G was estimated numerically for F2 to F4 assuming E 1 = 0. In the first step, was estimated for a certain RF excitation pulse and G with a freely available numerical Bloch solver using MATLAB (MathWorks, Natick, MA).32 Simulations were carried out with temporal resolution of 2 µs and spatial resolution of 80 µm with 2501 spatial points. The normalized envelope B 1(t) was scaled to achieve α(z = 0) = αξ in the center of the slice. Rather than estimating for each voxel with ξ and λ, calculations were accelerated by estimating in steps of Δξ = 0.05 followed by linear interpolation to Δξintp = 0.005. Variations of λ were incorporated by multiplying the sampling points along z with λ, to scale the thickness of the slice. In the last step, the integral along z for given G was solved by numerical integration. The source code can be found at: https://github.com/neuroimaging-mug.

Simulations

To investigate the influence of the G polarity and flip angle α on F3, simulations for α = 30° and α = 90° with negative and positive polarity of G were performed. Based on the vendor’s standard GRE pulse, a sinc‐Hanning‐windowed excitation pulse with a pulse duration Tpulse of 2 ms and a bandwidth time (BWT) product of 2.7 was chosen for the experiments. A G of 8.29 mT/m was determined with the Bloch solver to achieve a slice thickness ∆z of 4 mm, as defined by the FWHM of |M| for α = 30°. Based on the observed field gradients in phantom measurements, G was set to 100 µT/m for all simulations. In in vivo measurements of the brain, field gradients up to 300 µT/m have been reported in areas such as orbitofrontal cortex or inferior temporal lobe.33

Exploiting the relevance of individual parameters for modeling F4, a sensitivity analysis was performed for φ, , and λ with the same sinc‐Hanning‐windowed excitation pulse. To estimate the relevance of φ, simulations with G = 100 µT/m were carried out for F4 with varying α from 10° to 90°, each with positive and negative G polarity. Results were compared with simulations for model F2 considering only the magnitude |M| of the slice profile (φ = 0). For evaluation, the RMS error (RMSE) over time for each α between F4 and F2 was calculated.

The sensitivity for was simulated by scaling for each flip angle (α = 30° and α = 90°) with a factor ξ (ranging from 0.6 to 1.4) for G = 100 µT/m. The results for F4 obtained for different ξ were compared with those for ξ = 1 by plotting the RMSE. Same steps as for were carried out for λ by changing the value from 0.8 to 1.2.

A crucial assumption with the proposed models is that for a given α, TR is long enough to avoid changes of |M| due to incomplete T1 relaxation . Hence, the steady‐state solution in Equation 2 was included to estimate signal dephasing FT1 in the presence of G = 100 µT/m for different E 1:

For each TR/T1 (ranging from 1 to 5), the Ernst‐angle α was calculated and simulations with the sinc‐Hanning‐windowed RF pulse (Tpulse = 2 ms and BWT = 2.7) were carried out by setting α = α, α = 0.8 α, and α = 0.6 α. Obtained results were compared by calculating the RMSE over time between FT1 and F3.

Phantom experiments

To validate the results from the simulations of dephasing effects for different α and G, polarity phantom measurements were performed. For the phantom, a plastic cylinder (Ø = 12 cm and length = 20 cm) was filled with agarose gel (5 g/L), which was doped with 110 µmol/L MAGNEVIST to shorten the T1.

The phantom was scanned on a 3 T MRI system (Magnetom Prisma; Siemens, Erlangen, Germany) twice by a mGRE sequence with α = 30° and α = 90°, each with alternating polarity of G. The same sinc‐Hanning‐windowed excitation pulse (Tpulse = 2 ms and BWT = 2.7) as for the simulations was used, and |G| = 8.29 mT/m was used to achieve ∆z = 4 mm for α = 30°.

Other sequence parameters were as follows: FOV = 128 × 128 mm2, in‐plane resolution = 1 × 1 mm2, 32 monopolar echoes with bandwidth = 500 Hz/px, TE1 = 4 ms, ΔTE = 5 ms, TR = 3 seconds, 25 slices, with 0% interslice gap. For B 1 mapping, a Bloch‐Siegert sequence with the same resolution was used.34

The G map was obtained by using the central difference from the field map ΔB 0 to estimate the gradient in the ith slice:

Single side difference was used for the first (i = 1) and last slice (i = N). The value of ΔB 0 was estimated from a linear fit of the first 6 echoes of the unwrapped phase (PRELUDE unwrapping35). From the measured data, maps were estimated in MATLAB using the lsqnonlin() function with models S 1 to S 4.

As indicated in Supporting Information Figure S1, when varying G amplitude slightly within the model, it was found that results could be further improved when using G = 8.5 mT/m for all analyses.

Influence of TR/T1

Phantom measurements with different TRs (125 ms, 250 ms, 500 ms, 1 second, 1.5 seconds, 2 seconds, 3 seconds, and 5 seconds) and α (30°, 60°, and 90°) were carried out with the mGRE sequence to investigate steady‐state effects for modeling. A Bloch‐Siegert sequence was used for B 1 mapping.34 In addition, T1 was estimated with an inversion recovery sequence with 6 TIs (100 ms, 200 ms, 400 ms, 800 ms, 1.6 seconds, and 3.6 seconds), and the excited slice was measured with 0.5 × 0.5 mm2 in‐plane resolution by changing the readout direction to the slice direction. Results were evaluated by estimating with model S 4 for each TR and α.

In vivo and MWF experiments

To evaluate the proposed modeling for in vivo application, and MWF mapping experiments were performed on the same 3 T MRI system with 10 subjects (age range = 26‐50 years). The study was approved by the local ethics committee, and all subjects gave written informed consent. In addition, subjects were scanned with an anatomical MPRAGE with 1‐mm3 isotropic resolution for regional evaluation of and MWF maps.

For mapping, subjects were scanned twice with a mGRE sequence with alternating G polarity using a sinc‐Hanning‐windowed excitation pulse (Tpulse = 2 ms and BWT = 2.7) with α = 85° (Ernst angle assuming T1 = 1 second). Other sequence parameters were as follows: FOV = 256 × 208 mm2, in‐plane resolution = 1 × 1 mm2, |G| = 11.05 mT/m to achieve ∆z = 3 mm, 17 monopolar echoes with bandwidth = 500 Hz/px, TE1 = 2.87 ms, ΔTE = 3.59 ms, TR = 2.5 seconds, and 30 slices with 0% interslice gap. The last echo was a navigator echo at TEnavi = 65.4 ms, to correct for physiologically induced field variations.36 Then for each channel, the nth phase‐encoding line was corrected as described by Wen et al37:where ϕ and ϕ 1 are the mean phase values of the nth navigator echo and the reference phase of the first navigator echo, respectively. To account for phase accumulation after excitation, the estimated phase difference was scaled with the TE. Afterward, corrected k‐space data were combined with the method proposed in Luo et al.38

For B 1 mapping, a highly accelerated method based on the Bloch‐Siegert shift was used.39 The field map for calculating G was obtained from the difference of the unwrapped phase of the first and third echo divided by TE difference. From the data, maps were obtained using the models S 1, S 3, and S 4. The difference between the models was assessed regionally by calculating the mean and SD of in all subjects in gray matter and global white‐matter masks. Gray matter masks were segmented from the MPRAGE images with FSL FIRST,40 and the global white matter masks with SIENAX,41 part of FSL.42 All masks were affinely registered to mGRE space with FSL FLIRT.43, 44

For MWF mapping, all subjects were scanned with a slightly adapted mGRE sequence to account for the fast decaying myelin water component. Short echo spacing (ΔTE = 2.2 ms) was achieved with a bipolar readout gradient, which was inverted in a second acquisition to compensate for phase errors between even and odd echoes. Other sequence parameters were as follows: sinc‐Hanning‐windowed excitation pulse with Tpulse = 1 ms and BWT = 2, α = 85°, G = 14.15 mT/m, FOV = 255 × 105 mm2, in‐plane resolution = 1.14 × 1.14 mm2, ∆z = 4 mm, 27 bipolar echoes with bandwidth = 500 Hz/Px, TE1 = 2.37 ms, ΔTE = 2.2 ms, TR = 2 seconds, TEnavi = 63.8 ms, 25 interleaved slices with 0% interslice gap, and total scan time = 12 minutes. Again, a highly accelerated B 1 map was acquired.39

After correction of the data with the navigator echoes, the 2 mGRE images were registered using FSL FLIRT45 before averaging. The MWF estimation was based on a multi‐exponential relaxation times model46 with M = 200 water components:

Evaluation of data was performed by estimating MWF maps using models S 1, S 3, and S 4 with the nonnegative least squares algorithm of the MERA toolbox47 and a cutoff for myelin water my < 25 ms.48 For S 3 and S 4, the measured signal S was corrected with F3 and F4, respectively, before parameter estimation.

Regional evaluation of MWF maps was performed in white matter tracts with the JHU white‐matter atlas.49 The atlas was nonlinearly registered with FSL FNIRT to the MPRAGE images and transformed to the mGRE space using FSL FLIRT.43, 44 Before evaluation, masks were manually checked and adjusted with ITK‐SNAP.50

In a single scan session, 8 mGRE data sets were acquired from 1 subject (male, age = 29) using 4 different excitation pulses with α = 30° and 85° for each pulse. The first pulse was a 2‐ms‐long Gaussian pulse with σ = 280 µs (), and the other 3 were sinc‐Hanning‐windowed pulses with different BWT = 2, 2.7, and 8 and Tpulse = 1 ms, 2 ms, and 4 ms. The value of G = 10.56 mT/m, 18.87 mT/m, 11.05 mT/m, and 15.65 mT/m was estimated with the Bloch solver for ∆z = 3 mm and α = 30°. Other sequence parameters, as well as B 1 mapping, were as described for in vivo mapping. The differences between the pulses were assessed by estimating maps with S 4.

RESULTS

Figure 1 shows the simulation results for sinc‐Hanning‐windowed excitation pulse with positive and negative G polarity for α = 30° and α = 90°. It reveals that the polarity has no influence on |M(z)| of the slice profile (Figure 1A,B), whereas φ(z) is inverted when flipping polarity (Figure 1C,D). Consequently, F3 depends on the polarity of G (Figure 1E,F), an effect that is more strongly pronounced for α = 90°.

Figure 1

Simulation results for magnitude |M| (A,B) and phase φ (C,D) of the slice profile and the resulting dephasing F3 (E,F) with a macroscopic field gradient G = 100 µT/m for a sinc‐Hanning‐windowed excitation pulse (pulse duration Tpulse = 2 ms, bandwidth time product BWT = 2.7). For each α (top α = 30°, bottom α = 90°), simulations were performed with positive (red dotted line) and negative (solid blue line) G polarity. There is no difference in the magnitude (A,B) but the mirrored phase for α = 90° (D) causes different F3 (F)

The sensitivities of the model parameters φ(z), , λ, and TR/T1 are illustrated in Figure 2. When neglecting φ(z) in Figure 2A, the RMSE substantially increases for α > 40° with larger RMSE for negative G. For α = 90°, the RMSE is 5.5% for negative polarity and 4.5% for positive polarity, respectively.

Figure 2

Sensitivity analysis of the numerical model parameters. A, Comparison of the effect of including phase φxy in F4 versus a magnitude model F2 for positive and negative G polarity and G = 100 µT/m. B, Effects of variations in F4 with a macroscopic field gradient G = 100 µT/m. C, Influence of G on the slice encoding described with λ. D, The RMS error (RMSE) for neglecting T1 for different TR/T1 ratios is plotted assuming G = 100 µT/m. For each TR/T1, the RMSE was estimated between the F4 and FT1 for α, 0.8 α, and 0.6 α

The sensitivity for variations in Figure 2B depends strongly on the nominal flip angle α. For α = 30°, the RMSE was below 0.5% for all simulated values of ξ (α = α*ξ), with a moderate increase for α = 60° to 1% for ξ = 1.3. With 2.9%, the RMSE was 3 times higher for α = 90°.

The influence of λ on the signal is relatively small compared with and φ(z) with an RMSE of 0.8% for a strong G with 500 µT/m and minimal dependency on α (Figure 2C).

The simulated error due to neglecting T1 for different TR/T1 in Figure 2D shows an exponential decrease of the RMSE with increasing TR/T1 for all simulated flip angles. For all TR/T1 ratios, the highest RMSE was estimated when using the Ernst‐angle α and declines nonlinearly for 0.8 α and 0.6 α. For example, for TR/T1 = 1, the RMSE decreases from 2.8% to 1.8% to 1.2% for all simulated flip angles, whereas for TR/T1 = 2 the RMSE reduces from 1.2% to 0.8% to 0.5%.

When comparing the simulated errors by neglecting φ in Figure 2A with T1 effects in Figure 2D, the RMSE of φ becomes dominant with increasing TR/T1 ratio. Given that TR/T1 > 2, which results in α > 82°, the RMSE is smaller than 1.2%, whereas the RMSE due to neglecting φ is at least higher than 3.3% depending on the G polarity.

The values estimated with the mono‐exponential model S 1 are plotted as a function of G for α = 30° and α = 90° with positive and negative G polarity in Figure 3. The value of increases proportional to G for α = 30° (Figure 3A) with up to 8‐times higher values for G = 150 µT/m than for G = 0 µT/m. For α = 30°, negligibly small differences between the polarity of G and the sign of G were found, whereas for α = 90° (Figure 3B), positive and negative G yield different values and a dependency on the polarity of G. Moreover, Figure 3 shows the normalized averaged signal decay for |G| = 100 µT/m plotted with positive and negative G, explaining the difference in estimated values. For α = 90° with positive G and G > 0 (blue line), the signal decays faster than for G < 0 (red) and vice versa when switching G polarity.

Figure 3

Comparison of values estimated from the phantom experiments with the mono‐exponential model S 1 are plotted as function of G for α = 30° (A) and α = 90° (B) with positive and negative slice‐selection gradient G. Additionally, the averaged normalized signal decay is plotted for |G| = 100 µT/m. The dotted red line represents a positive G and the solid blue line represents a negative G. For α = 30°, no relevant differences between the polarity of G and G are observed, whereas for α = 90°, flipped G polarity substantially affects

Figure 4 compares the maps obtained from fits using models S 2, S 3, and S 4 for α = 30° (Figure 4A) and α = 90° (Figure 4B), each with positive and negative G polarity. In addition, the G map and B 1 map are illustrated in Figure 4C. Although results for α = 30° are comparable for all models, considerable differences for α = 90° between models and G polarity were found. When using only the magnitude |M| in model S 2 to estimate for α = 90°, it was not possible to recover without the influence of G. The values for G > 0 were overestimated for positive G and underestimated for G < 0, and switching to negative G polarity inverted the results. Extending the model S 2 by adding φ in S 3 yields better maps, which are influenced less by the G polarity. Additionally, including and λ in S 4 substantially improves maps, with minimal differences between G polarities. Further, estimated maps using S 4 with α = 90° are comparable with maps estimated from α = 30° for both G polarities.

Figure 4

Coronal and axial slices of estimated maps from the phantom measurements for different signal models (S2 ‐S4). Although all correction models yield relatively comparable values for α = 30° (A), the high flip angle results for α = 90° (B) highlight the effect of and λ correction. Full modeling with S 4 also eliminates the influence of the polarity of the slice‐selection gradient G at α = 90°. The corresponding G maps and B 1 maps are shown in (C)

Figure 5 illustrates the effects of neglecting T1 for signal modeling. Estimated maps with S 4 (Figure 5A) indicate an overestimation of , depending on TR and α in the presence of G (Figure 5B). For α = 30°, increased values are observable only up to a TR of 500 ms, whereas for α = 90° these effects extend up to a TR of 1.5 seconds. These TR values correspond to a TR/T1 ratio of 0.67 and 2.01 for the estimated T1 = 740 ms. The origin for the overestimation is shown in Figure 5C, where the averaged measured signal along the slice profile is plotted. Depending on α and TR, the steady‐state solution changes, causing a modeling error in the presence of G. Between different TRs for α = 30°, the profiles show less variations compared with α = 90°, leading to different signal dephasing for the same G. In addition to T1 effects, for TR > 2 seconds, SNR benefits can be observed for maps acquired with α = 90° compared with α = 30°.

Figure 5

Experimental evaluation of TR/T1 dependency for modeling in phantom measurements. A, Coronal maps were estimated using S 4 for different TR and α. The minimum TR required for avoiding T1 effects increases with the magnitude of G (B) and α. The value of T1 = 740 ± 86 ms was estimated with an inversion‐recovery sequence. C, The measured signal along the slice for each α and TR shows the different steady‐state solutions

In vivo experiments

In vivo results of maps obtained with models S 1 and S 4 are illustrated in Figure 6 for both G polarities. When comparing S 1 (Figure 6A,B) with S 4 (Figure 6D,E), much higher values are observed in maps using S 1 compared with S 4, thereby minimizing the effects of G. In addition, the difference map between positive and negative G polarity for each model reveals strong variations of values with up to 10 s−1 for S 1 in areas with strong G (Figure 6C). In contrast, maps estimated with S 4 substantially suppressed the effect of G polarity with difference values below 1 s−1 (Figure 6F).

Figure 6

Comparison of coronal and axial maps obtained from mono‐exponential model S 1 (A,B) with maps from the proposed numerical model S 4 (D,E) for positive and negative slice‐selection gradient G. C,F, Difference map between G polarities for each model. The S 1 model shows overestimation and substantial impact of the G polarity (C), which were mitigated using S 4 (F)

In Table 1 the regional evaluation of values with the corresponding mean |G| across all subjects is presented. Compared with the other models, the highest values were obtained with S 1 in all anatomical regions. In addition, the difference between G polarities increases with the mean |G| value in each region for S 1. For example, in the caudate nucleus, where the smallest |G| was observed with 20 µT/m, the difference between polarities is below 0.1 s−1, whereas in the brainstem it is 7.46 s−1 at a mean |G| of 89 µT/m. The values generally decrease when using S 2, but the difference between polarities slightly increases compared with S 1. Applying models S 3 and S 4 reduces the discrepancy between G polarities to a maximum of 2.01 s−1 and 1.25 s−1 in the brainstem. In all other regions the difference is much smaller, with values below 0.8 s−1. Between models S 3 and S 4, rather small changes can be observed generally.

Table 1

values (s−1) from models S 1 to S 4 in different brain regions for 10 subjects with the corresponding |G| values for positive and negative G

Region	Gslice	S 1	S 2	S 3	S 4	|G z| (µT/m)
Global WM	pos.	26.34 (1.16)	21.11 (0.61)	20.10 (0.58)	19.63 (0.62)	43.06 (8.81)
	neg.	23.72 (0.97)	18.17 (0.58)	19.89 (0.54)	20.17 (0.50)	43.68 (8.40)
Caudate Nucleus	pos.	23.36 (1.53)	21.46 (1.47)	21.81 (1.40)	21.82 (1.40)	20.47 (3.57)
	neg.	23.41 (1.61)	21.58 (1.32)	21.58 (1.28)	21.52 (1.27)	20.42 (3.22)
Pallidum	pos.	39.83 (2.78)	36.56 (2.58)	35.60 (2.65)	34.85 (2.71)	34.38 (8.75)
	neg.	36.86 (2.75)	33.50 (3.08)	35.07 (2.85)	35.61 (2.81)	34.03 (8.73)
Putamen	pos.	29.11 (2.14)	25.78 (1.70)	25.02 (1.76)	24.49 (1.80)	32.69 (5.84)
	neg.	26.97 (2.06)	23.52 (2.05)	24.91 (1.90)	25.28 (1.87)	32.90 (5.87)
Thalamus	pos.	25.84 (1.80)	22.34 (0.64)	21.33 (0.62)	20.34 (0.84)	33.65 (8.75)
	neg.	22.61 (0.97)	18.80 (1.24)	20.50 (0.87)	21.22 (0.74)	34.41 (8.83)
Brainstem	pos.	35.15 (7.97)	20.34 (2.07)	17.58 (1.77)	15.10 (1.60)	88.61 (35.73)
	neg.	27.70 (6.99)	11.21 (1.96)	15.08 (1.53)	16.45 (1.55)	89.90 (34.79)

The and |G| values are shown as mean (SD).

Abbreviations: neg., negative; pos., positive; and WM, white matter.

The difference between estimation with S 4 and S 3 is shown in Figure 7, pointing out the effect of modeling and λ in S 4. When visually comparing the difference maps in Figure 7A,B, a strong correspondence between the magnitude of G (Figure 7C) and (Figure 7D) can be observed for both G polarities.

Figure 7

Difference between maps estimated with S 4 (includes and λ variations) and S 3 for positive (A) and negative slice‐selection gradient G (B). Coronal (upper row) and axial (lower row) views are shown. C, B 1 map. D, G map. Depending on G polarity, varies in areas with higher and G variations

The maps from data acquired with 4 different excitation pulses and 2 different flip angles are shown in Figure 8. Visually, only minor differences among all maps are observable. Higher SNR can be observed in maps with α = 85° compared with α = 30°. Mean regional values are in good agreement after applying models S 3 and S 4 (Supporting Information Table S1). For example, in global white matter, the largest SD of between the acquisitions was found for S 1 with 1.59 s−1, due to the different pulses and flip angles. By using S 2, it decreases to 0.82 s−1, and for S 3 and S 4 the estimated values are 0.19 s−1 and 0.2 s−1, respectively.

Figure 8

maps estimated with model S 4 from multi‐gradient‐echo (mGRE) data acquired with 4 different excitation pulses (A‐D) for α = 30° (top row) and α = 85° (bottom row). Regional evaluation of can be found in Supporting Information Table S1

Figure 9 shows representative slices of MWF maps from 5 subjects obtained with models S 1, S 3, and S 4. It shows that with S 1, in areas with strong G, such as in the frontal and temporal lobe, the MWF estimation was not feasible, whereas the proposed approaches allowed a reconstruction in these areas. Between maps with models S 3 and S 4, no considerable differences were found, indicating that and λ have a neglectable small influence.

Figure 9

Representative myelin water fraction (MWF) maps from 5 subjects, obtained using models S 1 (A), S 3 (B), and S 4 (C). The proposed models S 3 and S 4 allow us to recover MWF values in areas strongly affected by the field gradient G (e.g., in frontal areas)

As shown in Figure 9, the MWF in the genu of the corpus callosum is underestimated with S 1 because of G. Using S 3 and S 4 enabled us to recover MWF values in these areas with a median of 12.09% and 12.66%, respectively. Our MWF results are within the range of reported values: For the genu of the corpus callosum, Lee et al27 reported approximately 12% for their postprocessing approach, and Alonso‐Ortiz et al26 reported approximately 16%. Furthermore, in the body of the corpus callosum, the proposed models yield to an increase of MWF from 3.7% with S 1 to 6.65% and 6.67% for S 3 and S 4, respectively. Interestingly, this analysis demonstrated that rather small |G| with around 10 µT/m in the body of the corpus callosum severely affects MWF estimation when using the simple model S 1. Supporting Information Table S2 summarizes the median MWF values in all 10 subjects in different white‐matter regions for models S 1, S 3, and S 4.

DISCUSSION

In this work we have introduced a numerical model for the signal dephasing of 2D mGRE sequences for arbitrary excitation pulses in the presence of a macroscopic field gradient G. In contrast to existing analytical solutions, our model is based on solving the Bloch equations numerically, which allows to estimate signal dephasing for any given flip angle α. We have shown that it is indispensable to consider the phase along the slice profile φ and the polarity of the slice‐selection gradient G for describing the signal dephasing for higher α. In our experiments, the threshold was approximately 60°, but this may also vary with the RF‐pulse shape.

Compared with existing models,19, 20, 22, 23, 24, 26, 27 which include the slice profile and assume linear varying macroscopic field variations, with the proposed model it is possible to explain different signal decays for different signs of G observed when using larger flip angles. As illustrated in Figure 1, this mismatch is explained by the phase variation φ along the slice profile, causing either a faster dephasing or a short period of rephasing followed by dephasing. Consequently, depending on the pulse shape and effective flip angle, the polarity of the gradient G must be included for modeling, as switching polarity inverts φ and thus signal dephasing.

In addition to the polarity dependency of G, the effects of variations and scaling of the slice profile with λ have been investigated in model S 4. However, changes in due to and λ were relatively small compared with S 3 (Table 1). Evaluation has been performed under the assumption that with an ideal model the estimated maps should be independent of G polarity. For the models S 1 and S 2, strong differences between G polarities were found, primarily due to φ, and by using S 3 it was substantially reduced, which indicates improved modeling. However, the main challenge for validation of the models was that, in vivo, no ground truth was available.

Another important aspect is the assumption that TR for a given α is sufficiently long to avoid T1 influence in the presence of G. The experimental results in Figure 5A are in accordance with the simulation results in Figure 2D, where the error decreases with TR/T1, and the minimum TR/T1 required enlarges with α. To gain SNR, it is desirable to use α Ernst, but care should be taken to prevent potential errors due to T1 and . By increasing TR/T1, both the α and the overall SNR increase; however, the errors due to are magnified. For example, as illustrated in Figure 2, when TR/T1 = 2, the error when neglecting T1 is about 1.2% for α = α = 77°. By comparing errors caused by variation, a deviation of ξ = 1.15 leads to errors in a similar range. Thus, without knowing T1, it is not possible to separate these effects, but it can be adjusted by the RF pulse shape. For instance, to estimate more accurately, longer RF pulses can be used to obtain a slice profile closer to the ideal, rectangular shape. This would have the advantage that signal dephasing is influenced less by and TR/T1, but it would lead to stronger φ variations and zero crossings due to the sinc‐shaped signal decay in the presence of G. However, for MWF estimation, very short pulses are needed, which will be more sensitive to these factors. Optimization of the RF pulses for specific applications was beyond of the scope of this work, but different pulses and their effects can be included and studied with the provided framework.

When comparing different modeling approaches, we can distinguish between models that fit parameters of F(t) from the signal decay19, 21 and models that use information from the pulse and field map to calculate F(t).22, 23, 24 Approaches that fit F(t) are more flexible in terms of model deviations from the ideal slice profile. For example, the sinc function used in the model approach by Fernandez‐Seara and Wehrli19 is well‐suited to model a variety of signal decays observed with different excitation pulses. Similarly, when modeling the macroscopic field as a quadratic function, the effects of a nonideal slice profile are inherently compensated for.21 However, in these models, the parameter estimation is often challenging due to the multiplication of F(t) with S tissue(t), thereby requiring the acquisition of many echoes. In contrast, with the analytical solution or our proposed numerical approach for F(t), only the parameters of the tissue model S tissue need to be estimated. Thus, if the properties of the RF pulse are available, a detailed description of F(t) is possible, favoring a closed or numerical solution. To select an appropriate model for a certain RF pulse and flip angle, the provided framework can be used to evaluate the expected error of different modeling approaches. If φ might be neglected for a specific RF pulse and flip angle, then an analytic solution yields a faster solution of F(t).

This work has similar limitations as other related postprocessing approaches.19, 20, 22, 24, 26, 29 The assumption of a linear varying magnetic field in slice direction might not hold in some areas with large susceptibility changes, which is especially pronounced at higher field strengths. However, as we have solved the dephasing along the slice direction by numerical integration, the model can also easily be adapted to describe the dephasing also for a quadratic varying magnetic field. Furthermore, in‐plane dephasing effects are neglected. In 2D acquisitions the slice thickness is usually much larger than the in‐plane resolution, but this might reduce accuracy in areas where the macroscopic in‐plane field variations are high. A possible solution to account for in‐plane dephasing could be to calculate the voxel spread function in‐plane as proposed by Yablonskiy et al51 and multiply the result with F3 or F4, respectively. Given that G is rather strong and that the signal dephasing is driven primarily by G, a reliable parameter estimation is difficult to achieve due to the fast signal decay. To overcome this issue, for MWF and it has been shown that z‐shim gradients between echoes can improve maps by rephasing the signal with appropriated compensation gradients.27, 52 Therefore, future work will focus on extending our model by including the moment of the z‐shim gradients in the modeling to describe the signal dephasing accordingly for every echo.

In addition to variations of the macroscopic field, variation of the phase offset φ 0 at TE = 0 could potentially influence signal dephasing. Contributions to φ 0 in phased array coils can be divided into receive coil–dependent (receive sensitivity ) and receive coil–independent (e.g., phase).53 To reconstruct the navigator‐corrected raw data, a multi‐echo approach was used to combine the individual coil data.38 In this approach, for each coil, images from all echoes are multiplied with the complex conjugate of the first echo, which removes inherently all components of φ 0 of the coil combined data. The development of the proposed models pointed out that the use of navigator echoes is highly recommended to compensate for phase errors arising from physiological fluctuations. As illustrated in Supporting Information Figure S2, depending on the subject’s reconstruction of parameter maps, not having the navigator echoes caused similar artifacts, as reported by Nam et al.54 If variations of φ 0 should be included, a ROEMER/SENSE reconstruction could be applied, as an example.55, 56

The scan time of the proposed applications is about 6 minutes for maps and 12 minutes for MWF maps. This is clinically acceptable for whole‐brain investigations, but further investigations will also focus on combination with accelerated imaging methods such as 2D CAIPIRINHA.57

CONCLUSIONS

Proper modeling of the signal dephasing in the presence of G for larger flip angles requires the consideration of |M| and φ with correct G polarity. Furthermore, and λ variations can potentially lead to a bias in the estimated model parameters, depending on the excitation pulse. Consequently, the proposed model allows to minimize the effects of G, which is highly relevant for accurate and MWF mapping of the entire brain based on 2D mGRE.

FIGURE S1 Coronal maps from the phantom measurements (α = 90°) estimated for a varying slice‐selection gradient G within the model. The most homogenous map was obtained with G = 8.5 mT/m

FIGURE S2 A,B, The MWF maps from 2 subjects. Maps are shown without and with correction of the raw data with the phase of the navigator echo

TABLE S1 Influence of pulse shape and flip angle for modeling . Note: The values (s-1) were estimated with models S 1 to S 4 from mGRE data acquired with 4 different pulses and α = 30° and α = 85°. It shows a flip angle and pulse shape dependency for S 1 in all regions. By applying S 2, differences decrease but values remain larger for α = 85° than for α = 30°. With S 3 and S 4, the flip angle dependency can be improved, leading to minimal differences of between the pulses. In the S 4 model, and λ have a small additional effect on estimation, compared with S 3

TABLE S2 Myelin water fraction values (%) with models S 1, S 3, and S 4 in different white matter regions for 10 subjects. Note: The MWF values are shown as median (interquartile range). The corresponding |G| values are listed as mean (SD)

Click here for additional data file.