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 R2∗ and myelin water fraction (MWF) mapping. Results The influence of macroscopic field gradients on R2∗ 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 R2∗ 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.

(MWF) from analyzing the multi-exponential decay of GRE signal. 4 Methods exploiting the phase evolution of the signal are SWI 5 or QSM. 6 Quantitative susceptibility mapping and R * 2 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 R * 2 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 z as an approximation of the macroscopic field variation along z, the signal over time is sinc-weighted proportional to G z and slice thickness. 12 Based on this theory, Fernandez-Seara and Wehrli 19 estimated G z and R * 2 iteratively from the measured signal decay, which later was refined by initialization of G z using the field map 20 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 z . To account for various pulse shapes, Preibisch et al proposed a solution in which the signal modulation due to G z 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 T 1 to avoid saturation of the slice profile. 22 To achieve a smooth signal decay, they have used exponential RF pulses 23 and compared them with sincshaped and sinc-Gauss shaped pulses for R * 2 mapping. 24 Similar to R * 2 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 z 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 B + 1 and the effect of scaling of the slice profile due to superposition of G z and the slice-selection gradient G slice . We further demonstrate with phantom and in vivo measurements that, depending on the pulse shape for larger flip angles, the polarity of G slice 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 R * 2 maps and MWF maps acquired with arbitrary excitation pulses and flip angles.

SOELLRADL Et AL.
When TR is much larger than T 1 , Equation 2 simplifies and |M xy (z)| is obtained by the sine of α eff times the equilibrium magnetization S 0 : Here, α eff (z) is obtained for a certain slice-selection gradient G slice and the applied excitation pulse with a certain shape and amplitude. For small flip angles, the slice profile α eff (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 xy (z)| and φ xy (z).
Extending previous studies, the factors ξ and λ were added to describe 2 effects that affect α eff (z) and therefore signal dephasing. First, variations of the active transmit field (B + 1 ) cause a deviation from the nominal flip angle α, which can change the effective flip angle profile α eff (z), and therefore requires α to be scaled with ξ, obtained from the normalized B 1 map. Second, G z is superimposed with G slice , 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 z , 4 different models have been studied. Summarizing Equation 1 in a tissue-specific signal component S tissue (t) (e.g., S tissue (t) = S 0 e −R * 2 t ) and a component F i (t) describing the signal dephasing due to Δω(z), the model S i (t) can be written as S i (t) = S tissue (t)F i (t). The 4 models are defined as follows: The model S 1 (t) serves as an uncorrected reference without modeling M xy (z) and Δω(z). Then, for S 2 (t), only the magnitude along the slice |M xy (z)| was considered neglecting φ xy (z). In S 3 (t), φ xy (z) was included, and in S 4 (t) the model was extended by additionally incorporating B + 1 and λ variations.

| Numerical implementation
Signal dephasing due to G z was estimated numerically for F 2 to F 4 assuming E 1 = 0. In the first step, M xy was estimated for a certain RF excitation pulse and G slice 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 α eff (z = 0) = αξ in the center of the slice. Rather than estimating M xy for each voxel with ξ and λ, calculations were accelerated by estimating M xy 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 z was solved by numerical integration. The source code can be found at: https ://github.com/neuro imagi ng-mug.

| Simulations
To investigate the influence of the G slice polarity and flip angle α on F 3 , simulations for α = 30° and α = 90° with negative and positive polarity of G slice were performed. Based on the vendor's standard GRE pulse, a sinc-Hanning-windowed excitation pulse with a pulse duration T pulse of 2 ms and a bandwidth time (BWT) product of 2.7 was chosen for the experiments.
A G slice 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 xy | for α = 30°. Based on the observed field gradients in phantom measurements, G z 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 F 4 , a sensitivity analysis was performed for φ xy , B + 1 , and λ with the same sinc-Hanning-windowed excitation pulse. To estimate the relevance of φ xy , simulations with G z = 100 µT/m were carried out for F 4 with varying α from 10° to 90°, each with positive and negative G slice polarity. Results were compared with simulations for model F 2 considering only the magnitude |M xy | of the slice profile (φ xy = 0). For evaluation, the RMS error (RMSE) over time for each α between F 4 and F 2 was calculated.
The sensitivity for B + 1 was simulated by scaling B + 1 for each flip angle (α = 30° and α = 90°) with a factor ξ (ranging from 0.6 to 1.4) for G z = 100 µT/m. The results for F 4 obtained for different ξ were compared with those for ξ = 1 by plotting the RMSE. Same steps as for B + 1 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 xy | due to incomplete T 1 relaxation E 1 = e −TR∕T 1 ≠ 0 . Hence, the steady-state solution in Equation 2 was included to estimate signal dephasing F T1 in the presence of G z = 100 µT/m for different E 1 : For each TR/T 1 (ranging from 1 to 5), the Ernst-angle α Ernst was calculated and simulations with the sinc-Hanningwindowed RF pulse (T pulse = 2 ms and BWT = 2.7) were carried out by setting α = α Ernst , α = 0.8 α Ernst , and α = 0.6 α Ernst . Obtained results were compared by calculating the RMSE over time between F T1 and F 3 .

| Phantom experiments
To validate the results from the simulations of dephasing effects for different α and G slice , 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 T 1 .
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 slice . The same sinc-Hanning-windowed excitation pulse (T pulse = 2 ms and BWT = 2.7) as for the simulations was used, and |G slice | = 8.29 mT/m was used to achieve ∆z = 4 mm for α = 30°.
Other sequence parameters were as follows: FOV = 128 × 128 mm 2 , in-plane resolution = 1 × 1 mm 2 , 32 monopolar echoes with bandwidth = 500 Hz/px, TE 1 = 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 z 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 unwrapping 35 ). From the measured data, R * 2 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 slice amplitude slightly within the model, it was found that results could be further improved when using G slice = 8.5 mT/m for all analyses.

| Influence of TR/T 1
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, T 1 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 mm 2 in-plane resolution by changing the readout direction to the slice direction. Results were evaluated by estimating R * 2 with model S 4 for each TR and α.

and MWF experiments
To evaluate the proposed modeling for in vivo application, R * 2 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-mm 3 isotropic resolution for regional evaluation of R * 2 and MWF maps.
For R * 2 mapping, subjects were scanned twice with a mGRE sequence with alternating G slice polarity using a sinc-Hanning-windowed excitation pulse (T pulse = 2 ms and BWT = 2.7) with α = 85° (Ernst angle assuming T 1 = 1 second). Other sequence parameters were as follows: FOV = 256 × 208 mm 2 , in-plane resolution = 1 × 1 mm 2 , |G slice | = 11.05 mT/m to achieve ∆z = 3 mm, 17 monopolar echoes with bandwidth = 500 Hz/px, TE 1 = 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 TE navi = 65.4 ms, to correct for physiologically induced field variations. 36 Then for each channel, the nth phase-encoding line S n (k x , TE) was corrected as described by Wen et al 37 : where ϕ n 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 z was obtained from the difference of the unwrapped phase of the first and third echo divided by TE difference. From the data, R * 2 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 R * 2 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 T pulse = 1 ms and BWT = 2, α = 85°, G slice = 14.15 mT/m, FOV = 255 × 105 mm 2 , in-plane resolution = 1.14 × 1.14 mm 2 , ∆z = 4 mm, 27 bipolar echoes with bandwidth = 500 Hz/Px, TE 1 = 2.37 ms, ΔTE = 2.2 ms, TR = 2 seconds, TE navi = 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 FLIRT 45 before averaging. The MWF estimation was based on a multiexponential T * 2 relaxation times model 46 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 toolbox 47 and a cutoff for myelin water T * 2 my < 25 ms. 48 For S 3 and S 4 , the measured signal S was corrected with F 3 and F 4 , 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 (B 1 (t) = e − t 2σ ), and the other 3 were sinc-Hanning-windowed pulses with different BWT = 2, 2.7, and 8 and T pulse = 1 ms, 2 ms, and 4 ms. The value of G slice = 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 R * 2 mapping. The differences between the pulses were assessed by estimating R * 2 maps with S 4 . Figure 1 shows the simulation results for sinc-Hanning-windowed excitation pulse with positive and negative G slice polarity for α = 30° and α = 90°. It reveals that the polarity has no influence on |M xy (z)| of the slice profile ( Figure 1A,B), whereas φ xy (z) is inverted when flipping polarity ( Figure 1C,D). Consequently, F 3 depends on the polarity of G slice ( Figure 1E,F), an effect that is more strongly pronounced for α = 90°.

| Simulations
The sensitivities of the model parameters φ xy (z), B + 1 , λ, and TR/T 1 are illustrated in Figure 2. When neglecting φ xy (z) in Figure 2A, the RMSE substantially increases for α > 40° with larger RMSE for negative G slice . For α = 90°, the RMSE is 5.5% for negative polarity and 4.5% for positive polarity, respectively.
The influence of λ on the signal is relatively small compared with B + 1 and φ xy (z) with an RMSE of 0.8% for a strong G z with 500 µT/m and minimal dependency on α ( Figure 2C).
The simulated error due to neglecting T 1 for different TR/T 1 in Figure 2D shows an exponential decrease of the RMSE with increasing TR/T 1 for all simulated flip angles. For all TR/T 1 ratios, the highest RMSE was estimated when using the Ernst-angle α Ernst and declines nonlinearly for 0.8 α Ernst and 0.6 α Ernst . For example, for TR/T 1 = 1, the RMSE decreases from 2.8% to 1.8% to 1.2% for all simulated flip angles, whereas for TR/T 1 = 2 the RMSE reduces from 1.2% to 0.8% to 0.5%.
When comparing the simulated errors by neglecting φ xy in Figure 2A with T 1 effects in Figure 2D, the RMSE of φ xy becomes dominant with increasing TR/T 1 ratio. Given that TR/T 1 > 2, which results in α Ernst > 82°, the RMSE is smaller than 1.2%, whereas the RMSE due to neglecting φ xy is at least higher than 3.3% depending on the G slice polarity.

| Phantom experiments
The R * 2 values estimated with the mono-exponential model S 1 are plotted as a function of G z for α = 30° and α = 90° with positive and negative G slice polarity in Figure 3. The value of R * 2 increases proportional to G z for α = 30° ( Figure 3A) with up to 8-times higher R * 2 values for G z = 150 µT/m than for G z = 0 µT/m. For α = 30°, negligibly small differences between the polarity of G slice and the sign of G z were found, whereas for α = 90° (Figure 3B) Figure 3 shows the normalized averaged signal decay for |G z | = 100 µT/m plotted with positive and negative G z , explaining the difference in estimated R * 2 values. For α = 90° with positive G slice and G z > 0 (blue line), the signal decays faster than for G z < 0 (red) and vice versa when switching G slice polarity. Figure 4 compares the R * 2 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 slice polarity. In addition, the G z 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 slice polarity were found. When using only the magnitude |M xy | in model S 2 to estimate R * 2 for α = 90°, it was not possible to recover R * 2 without the influence of G z . The R * 2 values for G z > 0 were overestimated for positive G slice and underestimated for G z < 0, and switching to negative G slice polarity inverted the results. Extending the model S 2 by adding φ xy in S 3 yields better maps, which are influenced less by the G slice polarity. Additionally, including B + 1 and λ in S 4 substantially improves R * 2 maps, with minimal differences between G slice polarities. Further, estimated R * 2 maps using S 4 with α = 90° are comparable with maps estimated from α = 30° for both G slice polarities. Figure 5 illustrates the effects of neglecting T 1 for signal modeling. Estimated R * 2 maps with S 4 ( Figure 5A) indicate an overestimation of R * 2 , depending on TR and α in the presence of G z ( Figure 5B). For α = 30°, increased R * 2 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/T 1 ratio of 0.67 and 2.01 for the estimated T 1 = 740 ms. The origin for the R * 2 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 z . Between different TRs for α = 30°, the profiles show less variations compared with α = 90°, leading to different signal dephasing for the same G z . In addition to T 1 effects, for TR > 2 seconds, SNR benefits can be observed for maps acquired with α = 90° compared with α = 30°.

| In vivo experiments
In vivo results of R * 2 maps obtained with models S 1 and S 4 are illustrated in Figure 6 for both G slice polarities. When comparing S 1 (Figure 6A,B) with S 4 ( Figure 6D,E), much higher R * 2 values are observed in maps using S 1 compared with S 4 , thereby minimizing the effects of G z . In addition, the difference map between positive and negative G slice polarity for each model reveals strong variations of R * 2 values with up to 10 s −1 for S 1 in areas with strong G z ( Figure 6C). In contrast, maps estimated with S 4 substantially suppressed the effect of G slice polarity with difference values below 1 s −1 ( Figure 6F).
In Table 1 the regional evaluation of R * 2 values with the corresponding mean |G z | across all subjects is presented. Compared with the other models, the highest R * 2 values were obtained with S 1 in all anatomical regions. In addition, the difference between G slice polarities increases with the mean |G z | value in each region for S 1 . For example, in the caudate nucleus, where the smallest |G z | 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 z | of 89 µT/m. The R * 2 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 slice 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.
The difference between R * 2 estimation with S 4 and S 3 is shown in Figure 7, pointing out the effect of modeling B + 1 F I G U R E 5 Experimental evaluation of TR/T 1 dependency for R * 2 modeling in phantom measurements. A, Coronal R * 2 maps were estimated using S 4 for different TR and α. The minimum TR required for avoiding T 1 effects increases with the magnitude of G z (B) and α. The value of T 1 = 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 and λ in S 4 . When visually comparing the difference maps in Figure 7A,B, a strong correspondence between the magnitude of G z ( Figure 7C) and B + 1 ( Figure 7D) can be observed for both G slice polarities.
The R * 2 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 R * 2 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 R * 2 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 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 z , 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 F I G U R E 6 Comparison of coronal and axial R * 2 maps obtained from monoexponential model S 1 (A,B) with maps from the proposed numerical model S 4 (D,E) for positive and negative slice-selection gradient G slice . C,F, Difference map between G slice polarities for each model. The S 1 model shows R * 2 overestimation and substantial impact of the G slice polarity (C), which were mitigated using S 4 (F) T A B L E 1 R * 2 values (s −1 ) from models S 1 to S 4 in different brain regions for 10 subjects with the corresponding |G z | values for positive and negative G slice differences were found, indicating that B + 1 and λ have a neglectable small influence.
As shown in Figure 9, the MWF in the genu of the corpus callosum is underestimated with S 1 because of G z . 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 al 27 reported approximately 12% for their postprocessing approach, and Alonso-Ortiz et al 26 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 z | 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 z . In contrast to existing analytical solutions, our model is based on solving the Bloch equations numerically, which allows to F I G U R E 7 Difference between R * 2 maps estimated with S 4 (includes B + 1 and λ variations) and S 3 for positive (A) and negative slice-selection gradient G slice (B). Coronal (upper row) and axial (lower row) views are shown. C, B 1 map. D, G z map. Depending on G slice polarity, R * 2 varies in areas with higher B + 1 and G z variations F I G U R E 8 R * 2 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 R * 2 can be found in Supporting Information Table S1 estimate signal dephasing for any given flip angle α. We have shown that it is indispensable to consider the phase along the slice profile φ xy and the polarity of the slice-selection gradient G slice 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 z observed when using larger flip angles. As illustrated in Figure 1, this mismatch is explained by the phase variation φ xy 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 slice must be included for modeling, as switching polarity inverts φ xy and thus signal dephasing.
In addition to the polarity dependency of G slice , the effects of B + 1 variations and scaling of the slice profile with λ have been investigated in model S 4 . However, changes in R * 2 due to B + 1 and λ were relatively small compared with S 3 (Table 1). Evaluation has been performed under the assumption that with an ideal model the estimated R * 2 maps should be independent of G slice polarity. For the models S 1 and S 2 , strong differences between G slice polarities were found, primarily due to φ xy , 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 T 1 influence in the presence of G z . The experimental results in Figure 5A are in accordance with the simulation results in Figure 2D, where the error decreases with TR/T 1 , and the minimum TR/T 1 required enlarges with α. To gain SNR, it is desirable to use α Ernst , but care should be taken to prevent potential errors due to T 1 and B + 1 . By increasing TR/T 1 , both the α Ernst and the overall SNR increase; however, the errors due to B + 1 are magnified. For example, as illustrated in Figure 2, when TR/T 1 = 2, the error when neglecting T 1 is about 1.2% for α = α Ernst = 77°. By comparing errors caused by B + 1 variation, a deviation of ξ = 1.15 leads to errors in a similar range. Thus, without knowing T 1 , it is not possible to separate these effects, but it can be adjusted by the RF pulse shape. For instance, to estimate R * 2 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 B + 1 and TR/T 1 , but it would lead to stronger φ xy variations and zero crossings due to the sinc-shaped signal decay in the presence of G z . However, for MWF estimation, very short F I G U R E 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 z (e.g., in frontal areas) 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 decay 19,21 and models that use information from the pulse and field map to calculate F(t). 22-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 Wehrli 19 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 φ xy 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 al 51 and multiply the result with F 3 or F 4 , respectively. Given that G z is rather strong and that the signal dephasing is driven primarily by G z , a reliable parameter estimation is difficult to achieve due to the fast signal decay. To overcome this issue, for MWF and R * 2 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 B − 1 ) and receive coil-independent (e.g., B + 1 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 R * 2 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 z for larger flip angles requires the consideration of |M xy | and φ xy with correct G slice polarity. Furthermore, B + 1 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 z , which is highly relevant for accurate R * 2 and MWF mapping of the entire brain based on 2D mGRE.