Adaptive slice‐specific z‐shimming for 2D spoiled gradient‐echo sequences

Purpose To reduce the misbalance between compensation gradients and macroscopic field gradients, we introduce an adaptive slice‐specific z‐shimming approach for 2D spoiled multi‐echo gradient‐echoe sequences in combination with modeling of the signal decay. Methods Macroscopic field gradients were estimated for each slice from a fast prescan (15 seconds) and then used to calculate slice‐specific compensation moments along the echo train. The coverage of the compensated field gradients was increased by applying three positive and three negative moments. With a forward model, which considered the effect of the slice profile, the z‐shim moment, and the field gradient, R2∗ maps were estimated. The method was evaluated in phantom and in vivo measurements at 3 T and compared with a spoiled multi‐echo gradient‐echo and a global z‐shimming approach without slice‐specific compensation. Results The proposed method yielded higher SNR in R2∗ maps due to a broader range of compensated macroscopic field gradients compared with global z‐shimming. In global white matter, the mean interquartile range, proxy for SNR, could be decreased to 3.06 s−1 with the proposed approach, compared with 3.37 s−1 for global z‐shimming and 3.52 s−1 for uncompensated multi‐echo gradient‐echo. Conclusion Adaptive slice‐specific compensation gradients between echoes substantially improved the SNR of R2∗ maps, and the signal could also be rephased in anatomical areas, where it has already been completely dephased.


| 819
SOELLRADL Et AL. effective transverse relaxation rate R * 2 serves as a proxy for iron concentration and has been used to study inflammatory and degenerative diseases in the brain, 12,13 iron overload in the liver, 14 and myocardial iron overload. 15 However, obtaining quantitative R * 2 values from gradientechoes is typically subject to quantification errors due to faster intravoxel dephasing caused by macroscopic field variations 16 (such as near air/tissue interfaces). In 2D spoiled GRE imaging, dephasing is more pronounced along the slice direction because the slice thickness is usually much larger than the in-plane resolution. Consequently, an effective approach to reduce signal dephasing is to decrease the slice thickness, which is accompanied by reduced SNR and prolonged acquisition time. 17 Various postprocessing methods for reducing these dephasing effects by considering the slice profile and macroscopic field variations have been proposed, [18][19][20][21][22][23] but at very strong field gradients the correction of the fast signal decay might not be feasible.
Alternatively, z-shimming approaches allow compensating signal dephasing due to a certain macroscopic field gradient G z by variation of the compensation gradient moment in the slice-selective direction. Starting from the basic principle of changing the amplitude of the slice-selective refocusing gradient demonstrated by Frahm et al, 24 different methods have evolved. To minimize the effects of G z on T ′ 2 quantification, Ordidge et al proposed the acquisition of images with different refocusing gradients. 25 Similarly, the GRE slice excitation profile imaging method acquires several images with an equidistant spacing of compensation gradients to estimate the k-space shift in z-direction, which allows reconstruction of a magnitude image with minimal contribution of G z . 26 The method was applied to R * 2 mapping, 27 and by adding compensation gradients between the echo acquisitions, a more efficient sampling is possible. 28 Nonetheless, a drawback of these approaches is that several additional images need to be acquired, which limits their application in clinical routine because of the prolonged scan time.
Similarly, Meng et al started with one strong compensation gradient before acquiring the echo train and successively rephased it with small opposed inter-echo gradients. 29 To reduce acquisition time, Wild et al proposed a single-scan method by applying repetitively a triple of compensation gradients between echo acquisition of a spoiled multi-echo GRE sequence (mGRE). 30 However, all of these methods assume an ideal slice profile, which would give rise to a sinc-shaped signal decay in the presence of G z . 16 Addressing variations of slice profile, Nam et al proposed a single-scan z-shim approach that includes the slice profile and compensates for a positive and negative G z by applying compensation gradients that are linearly scaled with TE. 31 To avoid signal crushing of the first echoes, Lee et al started the z-shim after the fifth echo for myelin water signal mapping. 32 A common limitation of the aforementioned approaches is that the compensation gradients are fixed for the entire FOV (global z-shim), and, consequently, a misbalance with the actual field gradient results in incomplete rephasing or even spoiling of the signal.
We therefore propose an adaptive slice-specific z-shimming approach to address spatial variations of G z . The corresponding slice-specific compensation gradients are estimated for each slice individually from a fast prescan. Additionally, a more effective z-shim pattern is introduced, in which six G z values are successively compensated between echo acquisitions. By adapting a signal modeling approach for 2D spoiled mGRE sequences, 23 we compare this novel approach, in terms of R * 2 mapping, with a global z-shim approach with linearly increasing moments 31,32 and a conventional mGRE sequence without z-shim gradients. Furthermore, to highlight the importance of adequate signal modeling in the presence of G z , R * 2 is also estimated from the conventional mGRE data with a more widespread used mono-exponential signal model.

| Signal modeling
Signal dephasing due to a field gradient G z can be compensated at a TE by applying a short compensation gradient with duration t c and amplitude G c , which results in a compensation moment m c = G c t c = −G z TE. In the case of a train of k compensation gradients, each with the amplitude G c [k] and identical t c , the accumulative moment M c [i] for the ith echo at TE i is given by: The sum of all applied compensation moments m c [k] until TE i is equal to a single theoretical mean compensation gradient G c [i] applied over the entire duration TE i . This allows us to superimpose G z and G c [i] for signal modeling independent of the shape and duration of the applied compensation gradients. Assuming a mono-exponential signal decay with R * 2 , the signal S TE i of the spoiled gradient echo is given by integrating the complex transverse magnetization M (z) weighted with the phase dispersion induced by both gradients along the z-direction: Figure 1 shows a 2D spoiled mGRE sequence ( Figure 1A) and a combination of the global z-shim patterns proposed by Nam et al and Lee et al 31,32 ( Figure 1B) along with the proposed slice-specific pattern presented in this work ( Figure 1C). In addition, Table 1 lists the corresponding compensation gradients G c [i] for the z-shim approaches for each echo.
The compensation moments for the global z-shim method ( Figure 1B) are calculated for a single positive G + c and negative G − c value. The first applied gradient moment after the fourth echo (m c [5] compensates for the effects of negative G − z followed by nulling the accumulative moment by inverting m c [5] (m c [6] = −m c [5]). This step is repeated for a positive G + z by applying a negative compensation moment (m c [7] . To avoid crushing of the signal in the first echoes, z-shim gradients are not applied for the first echoes as proposed by Lee et al. 32 Our work extends the compensation pattern in Figure 1B  [n] are compensated for three consecutive echoes, which is followed by a nulling of the total moment for the subsequent echo. To give an example, the moments m c [n, 5] to m c [n, 7] in the proposed pattern ( Figure 1C) are calculated as follows, assuming equal echo spacing ΔTE: Moreover, to allow a more effective rephasing, the nonzero value is split into 1 5 , 2 5 , 3 5 , 4 5 , 5 5   is zero. In addition to the inserted z-shim gradients, for all variants in Figure 1 a navigator echo is acquired after the last echo to compensate for physiologically induced field variations. 33

estimation
For all measurements, the complex-valued raw data were first corrected with the phase of the navigator echo as described by Wen et al, 34 followed by a coil combination using the method proposed by Luo et al. 35 Then, F z−shim TE i was calculated as described in Soellradl et al 23 for the model F 4 (t). In this model, M (z) is estimated for a certain RF pulse shape and G with a numerical Bloch solver. 36 Additionally, two potential factors that might affect M (z) were included: first, the nominal flip angle deviations due to the transmit RF field B + 1 and second, G z is superimposed with G , which leads to a change of the spatial encoding from z to z � = z with = G G z +G . 37 Thus, depending on the sign and amplitude of G z , the nominal slice thickness Δz is changed to Δz � , which is given by Δz � = Δz .
After the estimation of M (z), F z−shim TE i was calculated for each echo by substituting and S 0 were estimated by nonlinear fitting of the reconstructed magnitude data to Equation (2) using the lsqnonlin() function in MATLAB (MathWorks, Natick, MA).

| Sequence and model evaluation
The differences between the investigated sequences and the proposed signal modeling were assessed by calculating four different R * 2 maps: From the measured data of all three sequences, R * 2 was estimated with the signal model described previously. Additionally, R * 2 maps were calculated by fitting the standard spoiled mGRE data to a mono-exponential signal decay S TE i = S 0 e −R * 2 TE i .

| Phantom experiments
All experiments were carried out on a whole-body 3T MRI system (Magnetom Prisma; Siemens, Erlangen, Germany) using an eight-channel knee coil. To evaluate the proposed z-shim pattern, a homogenous phantom (5 g/L agar doped with 110 µmol/L Magnevist to shorten T 1 ) was built. Measurements with a spoiled 2D mGRE ( Figure 1A), a global z-shim pattern ( Figure 1B), and the proposed slice-specific z-shimming approach ( Figure 1C) were performed. To allow a comparison between the acquisition methods for the estimation of R * 2 , all sequence parameters were set identically-except the Global z-shim and G − c [n] are estimated for each slice n.
To increase coverage of compensated field gradients, these values are divided into three fractions or five fractions in case of G is being equal to zero. amplitudes of the z-shim gradients. A sinc-Hanning windowed RF excitation pulse (pulse duration T pulse = 2 ms, time-bandwidth product = 2.7) with flip angle α = 60° was used. In total, 20 echoes with a monopolar readout and a bandwidth = 500 Hz/Px were acquired. The echo spacing was 3.4 ms for the first four echoes without z-shim gradients, starting with TE 1 = 2.8 ms up to TE 4 = 12.9 ms. For the subsequent echoes with z-shim gradients (t c = 2 ms), the echo spacing was increased to 5.4 ms (TE 5 = 18.2 ms to TE 20 = 98.8 ms). After the 20th echo, phase encoding was rewound to acquire a navigator echo at TE navi = 103.4 ms. A total of 26 slices with a spatial resolution of 1 × 1 × 4 mm 3 (FOV = 128 × 128 mm 2 ) were acquired in an interleaved slice-acquisition scheme with a TR = 3 seconds. For all z-shimming approaches, z-shim gradients were applied with t c = 2 ms starting after the fourth echo. For the measurements with the global z-shim pattern ( Figure 1B), G +∕− c was set to ±100 T∕m. This value was approximately half of the maximum magnitude of the observed field gradients G z in the phantom. In addition to the mGRE sequences, a B 1 map was acquired using a Bloch-Siegert approach. 38

| Prescan to estimate
For the proposed z-shim approach, a prescan was done to estimate G +∕− c [n]. The prescan acquisition was performed with the same slice thickness (4 mm), an in-plane resolution of 2 × 2 mm 2 (FOV = 64 × 64 mm 2 ), three echoes with TE = 2.7 ms, 4.8 ms and 6.9 ms, and GRAPPA acceleration of 2. The phase data of the prescan was then automatically processed to estimate the positive G + c [n] and negative G − c [n] for each slice as follows: The phase data were unwrapped using PRELUDE, 39 and the field map was estimated by dividing the phase difference by the TE difference between the third and first echo. For evaluation, a mask was created by thresholding the magnitude image, which then was eroded with disk elements (radius of 5 pixels) to eliminate outliers close to the border. To estimate the field gradient map G z,pre , the gradient in z-direction of the field map was calculated in regions within the mask and smoothed with a 3D Gaussian filter (kernel size of 1). Then, the G z,pre map was quantized with an interval of 10 µT/m. For . Before scanning with the proposed z-shimming approach, the specific interecho compensation moments were calculated based on the pattern listed in Table 1.

| In vivo experiments
The proposed slice-specific z-shimming approach ( Figure 1C) was evaluated for in vivo R * 2 mapping by comparing the results with the global approach ( Figure 1B) and the approach without z-shimming ( Figure 1A). In total, 3 subjects were scanned on the same 3T MRI system using a 20-channel head coil. The study was approved by the local ethics committee, and all subjects gave written informed consent. For all experiments, the same RF excitation pulse as in the phantom measurements was used. Sixteen echoes and one navigator echo were acquired with TE 1 = 3 ms to TE 4 = 9.7 ms (without z-shim gradients, echo spacing = 2.2 ms), TE 5 = 13.9 ms to TE 16  [n] were estimated from a prescan as described for the phantom measurements, except that the mask was generated with the brain extraction tool BET, part of FSL. 39 Sequence parameters of the prescan were 35 slices with a voxel size of 2.3 × 2.3 × 3 mm 3 (FOV = 96 × 78 mm 2 ), three echoes with TE = 2.7 ms, 4.8 ms and 6.9 ms, and a GRAPPA acceleration factor of 3 with 20 reference lines, TR = 344 ms, α = 20°.
The acquisition times were 15 seconds for the prescan and 7 minutes 20 seconds for each of the three GRE sequences. In addition to the mGRE sequences, an MPRAGE sequence with 1-mm 3 isotropic resolution was acquired for anatomical segmentation. Furthermore, B 1 mapping was performed with a highly accelerated approach based on the Bloch-Siegert shift. 40 The different methods were compared by calculating the median and interquartile range of R * 2 values in global white-matter and gray-matter masks. The global whitematter masks were obtained from MPRAGE images using SIENAX, 41 part of FSL, 39 and subcortical gray-matter masks using FSL FIRST. 42 Regional R * 2 evaluation (median; interquartile range) was performed after affine registration to mGRE space with FSL FLIRT. 43,44 3 | RESULTS Figure 2 shows the signal decay of the three investigated pulse sequences within one slice. To demonstrate effects of varying G z , three regions of interest (ROIs) (Figure 2B) with different G z intervals were defined, and their normalized averaged signal decay S ( Figure 2C) and averaged F z−shim ( Figure 2D) were plotted. The standard spoiled mGRE sequence reveals a faster decay of S with increasing magnitude of G z , whereas for the z-shim approaches, S is differently rephased or dephased. For the global z-shim,  ). Depending on the G z interval of each ROI, the best compensation varies with TE for the proposed approach.

| Phantom
In Figure 3, S and F z−shim are plotted as a function of the TE for three different slices. In each slice, the values were averaged within ROIs of different G z interval. Similar to Figure 2, with the global approach, the best signal rephasing is achieved when G − c ≈ −G z ≈ −100 T∕m ( Figure 3B). In contrast, with the proposed approach the signal is gradually rephased for all slices for each block of compensation gradients (G within the ROIs. Therefore, after each fifth compensation gradient, the signal is nearly ideally compensated in each block. This is indicated when comparing S of the echoes TE 9 = 39.7ms and TE 15 = 71.9ms with F z−shim . Here, the dephasing functions F z−shim ≈ 1 suggesting an ideal compensation of G z . Furthermore, when comparing S between the slices, S is approximately equal for these echoes independent of G z . Figure 4 shows the estimated G z map ( Figure 4A) and the obtained R * 2 maps ( Figure 4B-F). The R * 2 map from the monoexponential fit of the standard spoiled mGRE ( Figure 4B) reveals a strong overestimation proportional to |G z |, which can be drastically decreased by accounting for G z in the signal model ( Figure 4C). Nonetheless, compared with the R * 2 value of 6.4 s −1 in the center of the phantom (G z is close to zero), R * 2 becomes underestimated with increasing |G z |. Applying a global z-shim (|G +∕− c | = 100 T∕m) improves the results, especially in areas with |G z | around 100 T∕m ( Figure 4D, slices 5 and 20). Figure 4E demonstrates that the proposed slice-specific approach yields more homogenous R * 2 maps over a wide range of G z values (eg, slices 2 and 23).  Figure 5 shows the averaged R * 2 values of the phantom with a bin size of G z = 10 T m as a function of G z , and demonstrates the difference between the proposed approach and the global z-shimming. Although the global z-shim approach (|G +∕− c | =100 T m ) corrects R * 2 values at about |G z | = 100 T m to the expected value of 6.4 s −1 (R * 2 value at G z ≈ 0 T m ), the proposed approach yields constant R * 2 values over a broad range of G z from −150 T m to 125 T m . Furthermore, the results from the mono-exponential fit of the standard spoiled mGRE data clearly show the strong increase of R * 2 with |G z |. Figure 6 shows representative mGRE images for the three investigated sequences (12th to 16th echo). For the spoiled mGRE sequence ( Figure 6A), a faster signal decay in areas with strong G z , for example, close to the nasal cavities, can be observed. For all sequences, the 12th echo images as well as the 16th echo images are equal because of a zero net moment (M c,12 = 0 and M c,16 = 0). Between these two echoes, the signal in various brain areas is differently rephased and dephased, depending on the z-shim approach and G z . The global z-shim pattern with |G +∕− c | = 220 T∕m shows that negative G z values and positive G z values are rephased at the 13th and 15th, respectively ( Figure 6B). Instead of single positive and negative G z , a larger range of G z values can be covered by the proposed approach ( Figure 6C) (red arrows).

| In vivo
The R * 2 maps in Figure 7 demonstrate improvements in areas with strong G z from the global z-shim pattern using constant |G +∕− c | = 220 T∕m ( Figure 7C) over the spoiled mGRE ( Figure 7A,B), which are most pronounced in the temporal is zero for all sequences. With the proposed approach, the signal can also be rephased in areas where it has already been completely dephased (arrows). The complete series of the echoes with z-shim gradients is illustrated in Supporting Information Figure S2 (A) lobe and cerebellum (slice 3) or close to the sinuses (slice 9). Further improvements and additionally increased SNR are observed in the R * 2 maps obtained with the proposed adaptive z-shim ( Figure 7D). Table 2 summarizes global and regional R * 2 values for all subjects. In line with the visual assessment in Figure 7, the interquartile ranges are smallest for the proposed approach in white matter, followed by the global z-shim method. The largest interquartile ranges were obtained without z-shimming and when assuming a mono-exponential signal model.

| DISCUSSION
We have introduced an adaptive slice-specific z-shimming approach that allows one to minimize effects of macroscopic field gradients in the slice-selection direction in 2D mGRE sequences. For each slice n, a maximum positive and negative compensation gradient G +∕− c [n] is obtained from a fast prescan. To increase the coverage of compensated G z values, is split into three fractions: .

F I G U R E 7
Axial views of estimated in vivo R * 2 maps (left), with detailed views of the blue rectangular regions (right). A, The R * 2 maps were directly calculated from the spoiled mGRE data by assuming a mono-exponential signal model neglecting G z (F z−shim = 1). The other R * 2 maps were calculated using the proposed signal model for the spoiled mGRE (B), the global z-shim (|G +∕− c | = 220 T∕m) (C), and the proposed slice-specific approach (D) data. An increase in SNR can be observed from (C) to (D) due to higher signal recovery Based on these gradient values, a pattern of compensation moments between the echoes is calculated ( Figure 1C).
Our novel adaptive slice-specific z-shimming was compared with a conventional spoiled mGRE sequence and a global z-shimming approach that applies a positive and negative G +∕− c ( Figure 1B) independent of the slice position. 31,32 In contrast to modeling of the standard spoiled mGRE, the global z-shim enables us to recover R * 2 values in areas with strong G z , which is in line with the results of Nam et al. 31 By performing slice-specific z-shimming with more compensated G z values, the proposed approach results in SNR improvements (Figure 7). Quantitatively, the measured values are within the range of reported values in the literature at 3 T. The z-shim approach by Nam et al yielded an R * 2 of 20.77 s −1 for the putamen and 34.22 s −1 for the globus pallidus, 31 which is close to the mean values of our 3 subjects with 24.08 s −1 and 35.05 s −1 . When considering the age of the subjects, our R * 2 values are in good agreement with a study reporting different age ranges. 45 Subjects' regional R * 2 values in the caudate nucleus, thalamus, and brainstem are within the 95% confidence interval of this study. 45 For subjects 1 and 3, the R * 2 values in the globus pallidus are slightly above the 95% confidence interval as well as in the putamen for subject 3. [n] from the prescan field gradient map G z,pre [n], splitting of the compensation gradients into different magnitudes was performed. When using a single value (eg, maximum and minimum of positive and negative G z,pre [n]), improvements were only observed in areas with G z values close to the specific compensation gradient. To demonstrate this relation, additional measurements with a slice-specific approach but with a single G +∕− c [n] were performed. As shown in Supporting Information Figure S1, splitting G +∕− c [n] led to a more robust compensation over a wide range of G z values. A further refinement of our approach could be made by passing the desired compensation gradient for each echo G +∕− c n, TE i to the sequence. This comes with the advantage that the compensation gradients can be individually selected, based on the distribution of G z values in each slice.
Z-shim approaches primarily aim to avoid signal dephasing in areas with large G z . In this context, a rather unexpected finding was that also areas with relatively low field gradients (|G z | < 50 T∕m) yielded higher SNR in R * 2 maps by applying small compensation gradients compared with postprocessing-only methods (Figure 7, slice 24). This SNR increase might be especially promising for combined applications with acceleration methods such as parallel imaging 46,4748 The proposed approach has some limitations. First, a prescan with a duration of 15 seconds is necessary to estimate G z . However, this additional scan time is minimal compared with the fully sampled z-shim acquisition (7 minutes 20 seconds) itself, and the increase in SNR compensates for the prolonged scan time. Another issue, especially in vivo, is the estimation of a reliable field gradient G z | < 50 T∕m map from the prescan, which is used to define R * 2 . Here, we focused on a robust implementation and avoided potential gradient errors due to missing field-map values in the skull by eroding the G +∕− c [n] map. Nevertheless, it might result in nonoptimal compensation gradients in these areas. An alternative might be to match the slice position to a template G z map instead of performing a prescan. 49 This work focuses on z-shimming because the signal dephasing is major along the slice-selective (z-)direction compared with the orthogonal directions. In addition, strong in-plane field gradients can be considered by calculating additional compensation moments in in-plane directions or, as proposed by Yablonskiy et al, 50 by modeling the signal dephasing with the voxel spread function.
We have recently introduced a signal modeling approach for an arbitrary excitation pulse and Fz−shim , 23 which has been adapted in the current work to describe signal dephasing F z−shim due to G z and the compensation gradient G c . Because R * 2 is estimated from the measured data by nonlinear fitting of Equation (2), any modeling error in F z−shim will propagate into the R * 2 estimate. Here, B + 1 and have been considered for modeling, but additionally, the ratio TR∕T 1 can affect F z−shim . If the assumption TR ≫ T 1 is not fulfilled, M (z) changes according to the steady-state equation for spoiled gradient-echo sequences 51 and might bias F z−shim . To better assess the contributions of B + 1 , , and TR∕T 1 to F z−shim , additional simulations were carried out for different G z values (Supporting Information Figure S3). For a ratio of TR∕T 1 = 5, T 1 effects are negligibly small, while errors due to B + 1 increase with . Compared with B + 1 , the estimated errors caused by are similar for each . In contrast, for TR∕T 1 = 2, a T 1 bias can be observed, which is small compared with the B + 1 error. To investigate the influence of B + 1 and in vivo, Supporting Information Table S1 lists the results without considering B + 1 and . It reveals that the greatest relative change of R * 2 for the proposed approach was 2.7% for subject 3 in the brain stem (Supporting Information Table S2). These small changes in R * 2 suggest that B 1 mapping might not be necessary for the regions evaluated. However, when increasing or when evaluation regions with stronger G z , accounting for B + 1 might be beneficial. Based on the simulation results, a potential small T 1 effect cannot be excluded with the TR = 2.5 seconds used in vivo.

SOELLRADL Et AL.
Other sources for model deviations in F z−shim are the input parameters G z and G c . Similar as for the prescan, G z estimation is challenging if the field map values from adjacent slices are missing. For G c , it is assumed that it is ideally characterized by the actual applied gradient moment of the MRI system. Thus, errors might occur in case of gradient imperfections or when a different MRI system is used. In our experiments, a good correspondence between the predicted signal dephasing F z−shim and the measured signal S (Figures 2 and 3) was observed, indicating a reasonable accurate G c for the proposed approach.

| CONCLUSIONS
A new adaptive slice-specific z-shim approach in combination with signal modeling for 2D mGRE data was introduced to minimize the effects of macroscopic field gradients. The proposed approach allows a more robust correction of R *