Robust 3D Bloch‐Siegert based B1+ mapping using multi‐echo general linear modeling

Purpose Quantitative MRI applications, such as mapping the T1 time of tissue, puts high demands on the accuracy and precision of transmit field (B1+) estimation. A candidate approach to satisfy these requirements exploits the difference in phase induced by the Bloch‐Siegert frequency shift (BSS) of 2 acquisitions with opposite off‐resonance frequency radiofrequency pulses. Interleaving these radiofrequency pulses ensures robustness to motion and scanner drifts; however, here we demonstrate that doing so also introduces a bias in the B1+ estimates. Theory and Methods It is shown here by means of simulation and experiments that the amplitude of the error depends on MR pulse sequence parameters, such as repetition time and radiofrequency spoiling increment, but more problematically, on the intrinsic properties, T1 and T2, of the investigated tissue. To solve these problems, a new approach to BSS‐based B1+ estimation that uses a multi‐echo acquisition and a general linear model to estimate the correct BSS‐induced phase is presented. Results In line with simulations, phantom and in vivo experiments confirmed that the general linear model‐based method removed the dependency on tissue properties and pulse sequence settings. Conclusion The general linear model‐based method is recommended as a more accurate approach to BSS‐based B1+ mapping.

and an associated phase accumulation, which is proportional to the square of the pulse amplitude thereby encoding the B + 1 field. This technique performed favorably in a recent review of the accuracy, precision and practicality of a range of prominent B + 1 mapping techniques, 7 and has been shown to be less sensitive to B 0 inhomogeneities and chemical shifts 8 than other phase-based methods.
The BSS technique is flexible and can be integrated into a multitude of pulse sequences, such as 2D 9 or 3D gradient echo (GRE), 8,10,11 interleaved echo planar imaging, spiral GRE, 12 and spin echo [13][14][15][16] acquisitions. The BSS technique, as typically implemented, requires two acquisitions with opposite off-resonance frequencies. By subtracting the 2 phase images, the BSS effect is enhanced and unrelated phase components are removed, e.g., phase accumulated across echo time (TE) due to B 0 inhomogeneity or chemical shifts, due to eddy currents, or any initial phase due to the transmitting and receiving coils. This subtraction also has the advantage of removing the effect of B 0 inhomogeneity on the BSS, up to first order. 9 These 2 acquisitions can either be played out sequentially or by interleaving the opposite off-resonance frequencies. Previous reports have shown the interleaved approach to be more robust to motion 17 and magnetic field drift. 10 In this work, we focus on a 3D spoiled GRE implementation and demonstrate, through both simulations and experiments, that alternating the sign of the off-resonance frequency from shot to shot in an interleaved manner disturbs the steady-state and introduces an additional phase difference between the 2 acquisitions, especially at short TR. We show that this additional phase difference leads to bias in the B + 1 map that depends on the relaxation parameters of the studied tissue, the specifics of the RF spoiling regime and the actual B + 1 amplitude. We additionally propose and validate a modified BSS-based approach that removes these dependencies. The solution consists of a multi-echo acquisition in which several echoes are acquired before and after the BS pulse and modeling the phase evolution with a general linear model (GLM). We demonstrate that this GLM-based approach to isolating the BSS phase allows the interleaved approach to be used without introducing any error, extending the acquisition time, increasing the specific absorption rate (SAR) or reducing the sensitivity.

| THEORY
In line with Sacolick et al, 9 and as detailed in the Supporting Information, which is available online, the BS phase introduced by an RF pulse with peak amplitude B p 1 and normalized shape B norm 1 , is proportional to the square of the peak pulse amplitude: off is the off-resonance frequency of the pulse and Δ B 0 is the local field inhomogeneity, both in Hz.

| The classic method: Isolating the BSS phase by subtraction
The classic approach to BSS-based B + 1 mapping consists of acquiring 2 datasets with BS pulses of opposite off-resonance frequency (i.e., + off and − off ). The phase difference between these is: Because the first order terms that depend on Δ B 0 cancel, this expression simplifies to: Previously, this subtraction was assumed to also remove any phase accumulated from other sources, such as eddy-currents, transmit/receive-related phase offsets, chemical shifts, and local B 0 inhomogeneities. However, crucially, this is only true if the additional phase components are identical for each of the off-resonance frequencies. If this assumption is violated, the B + 1 estimate will be erroneous.

| The GLM method: Isolating the BSS phase by modeling a multi-echo acquisition
We propose an alternative to the classic BSS approach that computes accurate B + 1 maps even if conditions vary between the two off-resonance frequency acquisitions. This approach relies on two novel features: a dual-offset multi-echo sequence and a GLM.

| Dual-offset multi-echo sequence
In the modified BSS-based B + 1 mapping sequence (Figure 1), multiple echoes are acquired after one excitation pulse. Two echoes, after the BS pulse, have previously been used 18 to concurrently compute the B 0 inhomogeneity, whereas here (1) multiple echoes are acquired before and after the BS pulse. As in the classic method, a second acquisition is performed with the opposite off-resonance frequency, either sequentially or interleaved.

| GLM
The GLM approach models the phase data, Y, as multiple linear sources of phase accumulation over time: Each row of the model matrix X, corresponds to a single echo. The following explanatory variables model distinct sources of phase evolution in separate columns: • X BSS : Models the phase accumulated due to the BS pulse, specifically the first term of the sum in Equation 1 (only). This phase should only be present after the BS pulse, and should change sign with the off-resonance frequency. The regressor consists of zeros for echoes before the BS pulse and either 1 or −1 afterward depending on the BS pulse frequency.
• X Δ B 0 : Models the phase accumulated due to local B 0 inhomogeneities. Regressor values are TEs. • X Odd∕Even : Models the phase difference between odd and even echoes due to eddy currents generated by the bipolar readout gradients. X Odd∕Even is 1 for odd echoes and -1 for even echoes. • X Offset + and X Offset −: Models the initial phase offset of the acquisitions with positive and negative off-resonance frequencies, respectively. X Offset + is 1 for all echoes from a TR with a BS pulse with positive off-resonance frequency, and 0 for all echoes from a TR with a BS pulse with negative off-resonance frequency; vice versa for X Offset −. • X SameSign : Models phase consistently accumulated during the BS pulse and the crushers, regardless of the sign of the BS pulse's off-resonance frequency. This includes the second term in Equation 1 and, for example, any phase due to eddy currents generated by the crushers. X SameSign is 0 for echoes before the BS pulse and 1 for echoes after The sequence diagram of the modified 3D multi-echo GRE for BSS-based B + 1 mapping. Two echoes are acquired before, and 6 echoes after, the BS pulse, which is flanked by crushers in 1 phase-encoding direction (PE 1 ) to destroy any inadvertent on-resonance excitation and minimize dependence on excitation flip angle. The gradients on each axis are balanced before the BS pulse. In this example, 8 echoes are acquired for each off-resonance frequency, resulting in a total of 16 phase images from which the B + 1 efficiency is mapped. (B) A GLM is used to model the phase variation across TEs. (C) Typical maps of model coefficients obtained in vivo exemplify the phase accrued due to the BSS ( BSS ), B 0 field inhomogeneity, Δ B 0 , alternating readout (RO) polarity, Odd∕Even , initial phase offsets specific to the off-resonance frequency of the pulse, Offset + , Offset − , and any additional phase due to the presence of the block of crushers and the BS pulse that is independent of the sign of the off-resonance frequency of the pulse, SameSign the BS pulse, regardless of the pulse's off-resonance frequency.
The parameters constitute the regression coefficients and are estimated voxel-wise by means of the weighted leastsquare approach: ̂= X T WX −1 X T WY. W is a diagonal weighting matrix in which the weights are the magnitude of the echoes. ̂B SS , the parameter of interest, is considered an estimate of Φ BSS . ε is an error term.

| Numerical simulations
Simulations were used to evaluate the difference, both before and after the BS pulse, between two acquisitions with opposite off-resonance frequencies. Any difference present at this point would introduce an error in B + 1 estimated with the classic method. A typical GRE acquisition was simulated by a series of matrix operations (described in Supporting Information). Several configurations were investigated.

| Effect of RF Spoiling increment
RF spoiling modifies the phase of the excitation pulse and the BS pulse ( ) so that the phase difference between successive TR periods increases linearly by a constant amount BaseInc .
The impact of the RF spoiling increment was investigated by changing BaseInc from 0° to 180° in 10° steps. A phase increment of 0 corresponds to no RF spoiling (i.e., = 0 for all pulses).

| Sequential or interleaved
For an interleaved acquisition, the sign of off was alternated between successive BS pulses. For sequential acquisitions, the sign was switched after half the total number of pulses N exc 2 , at which point RF spoiling was also reset. For the Classic method, Φ BSS was taken to be half the difference in phase accumulated after the BS pulse of the acquisitions with opposite off-resonance frequency, and termed Φ Classic BSS (see Supporting Information). For the GLM method, Φ BSS was taken to be the mean absolute phase accrued during the BS pulses of the acquisitions with opposite off-resonance frequency, and termed Φ GLM BSS (see Supporting Information).

| MR pulse sequence
Measurements were performed at 3T (Siemens Prisma) using a body coil for transmission and a 32-channel head coil for signal reception using an in-house MR pulse sequence ( Figure 1). A Fermi pulse of duration T = 2 ms imparted the BSS after the second echo. Acquiring 2 echoes before the BS pulse served to minimize the correlation of the X BSS regressor with other regressors while maintaining a reasonable TE for the echo after the BS pulse. The encoding gradients on all axes were balanced immediately before the BS pulse, to ensure the same dephasing state for the magnetization across TRs, and played again just after. Crusher gradients were played out either side of the BS pulse, concurrently with the balancing/phase-encoding gradients (on PE 1 ), to crush any undesired on-resonance excitation and to minimize any dependence on the excitation flip angle. As demonstrated by Duan et al, 18 perfect dephasing of the transverse magnetization before the BS pulse is required to fully remove any such dependence. Sensitivity to nonideal conditions is reduced by using a high crusher moment because the greater the dephasing, the smaller the dependence on the excitation flip angle. Therefore, a relatively large crusher moment, 19 designed to generate a theoretical dephasing of 6 rad across a voxel, was used. An excitation flip angle of 15°, corresponding to the Ernst angle for a TR of 35 ms and a T 1 of 1000 ms, was chosen to maximize the precision of the B + 1 mapping. 18 Although this may be somewhat sub-optimal from a precision perspective for the phantom acquisitions and the long TR in vivo acquisitions, it was used for all acquisitions to ensure consistency.
For sequential acquisitions, the positive and negative offresonance frequency pulses were played out in consecutive blocks. In the interleaved case, the off-resonance frequency was alternated across successive TR periods. To achieve steady state, 200 dummy cycles were executed before each block in the sequential case, 400 cycles were used at the outset of the interleaved case (to match acquisition times), and in both cases a spoiler gradient, set to reach a dephasing of 6 at the end of the TR, was applied in the readout direction after the last echo. The RF spoiling was reset at the end of the first block in the sequential case. Table 1 lists sequence parameters of all experiments. (4)

map estimation
All data, including B + 1 maps, were reconstructed in real time using in-house code implemented in Gadgetron. 20 An apodizing filter was applied to the raw k-space data along each dimension to minimize ringing artefacts. After Fourier transformation the images were adaptively combined across coil elements. 21 Two B + 1 maps were computed for each dataset using the Classic and GLM methods, respectively. For the Classic method, the phase of the third echoes, which were acquired after the positive and negative off-resonance frequency BS pulses, were subtracted to estimate the BSS phase: For the GLM method, the phase images for each off-resonance frequency were temporally unwrapped, by spatially unwrapping the differences between successive echo pairs, 22 and cumulatively adding these to the phase of the first echo. The phase images were subsequently used to estimate the parameters of the GLM (Figure 1). The BSS phase Φ BSS was captured by the first regressor X BSS of the model matrix such that Φ BSS =̂B SS .
For both methods, B p 1 was computed on a voxel-wise basis using Equation 4.

| Phantom experiments
B + 1 maps were acquired on an FBIRN gel phantom. 23 Experiments 1 and 2 used both the Classic and GLM processing methods to construct B + 1 maps. B + 1 errors were quantified as the percent difference of the estimated B + 1 with respect to a reference B + 1 map..

Phantom experiment 1: Comparing processing approaches
This experiment probed the impact of the processing method as a function of the RF spoiling increment, B + 1 maps were created using both Classic and GLM methods and compared with the reference map. Histograms of the relative difference in B + 1 estimates were calculated for each RF spoiling condition and processing method. Two B + 1 maps were estimated for each scan using the Classic and GLM methods, respectively. The percent difference in B + 1 was calculated with respect to the reference map for the same processing method. To compare processing methods, the percent difference between the two B + 1 maps derived from the reference acquisition was also computed, with respect to that obtained with the Classic method.

RF spoiling
For sequential acquisition ordering, once steady-state was reached for each off-resonance frequency block (Figure 2A), the phase difference between these was zero before the BS pulse ( Figure 2C). It was nonzero after the BS pulse ( Figure 2E) because it contained the BSS phase. However, this BSS phase estimate remained constant over time and was independent of the RF spoiling condition. For interleaved acquisitions, the phase before the BS pulse varied from pulse to pulse regardless of pulse number ( Figure 2B). This temporal variance caused a nonzero phase difference between the interleaved acquisitions with opposite BS pulse frequencies, both before ( Figure 2D) and after ( Figure 2F) the BS pulse. Furthermore, the phase difference after the BS pulse, i.e., the estimate of the BSS phase, differed from the phase difference obtained with the sequential approach and depended strongly on the RF spoiling conditions.

| Effect of RF spoiling increment
The error in the BSS phase resulting from simulating an interleaved acquisition and using the Classic method depended strongly on the RF spoiling increment used (Figure 3). The greatest errors were predicted for phase increments of 0° (equivalent to no RF spoiling) and 180°. Large errors were F I G U R E 2 Numerical simulation results. Phase accrued before the BS pulse in case of sequential acquisitions (A) or interleaved acquisitions (B), in the presence (red and blue curve) of RF spoiling Φ BaseInc = 50 • and 117 • or without RF spoiling (green). Phase difference before the pulses of the two acquisitions with opposite frequencies (C,D) and phase difference after the BS pulse (E,F) also predicted for phase increments of 60° and 120°, whereas no error was predicted at 90°.

, T 2 , and T 1
Increasing TR greatly reduced the predicted error in Φ BSS ( Figure 3A). Longer T 1 times led to larger predicted errors, but had less impact than TR ( Figure 3C). Similarly, longer T 2 times predicted larger errors ( Figure 3D), especially with phase increments of 60 and 120° where the predicted error was already large.
The amplitude of the BS pulse also had a small impact whereby the relative error was predicted to be lower for higher amplitude pulses ( Figure 3B). Note that in this case B + 1 also increases, so while the predicted relative error decreased, the absolute error is actually increased.

BSS estimation
The numerical simulations predicted that these errors were removed by using the GLM approach. As a result, the derived B + 1 estimates were predicted to be stable and agree with the B + 1 estimated from the sequential case using the Classic method. This was the case regardless of the RF spoiling conditions, TR, T 1 , T 2 , or BS pulse amplitude.

| Phantom experiments
Phantom experiment 1: Comparing processing approaches With RF spoiling increments of 0 and 180°, the B + 1 estimates were very different depending on the processing method used (Figure 4), but these estimates agreed when the RF spoiling increment was 90°. Two further peaks in the discrepancy were observed with RF spoiling increments of 60° and 120°. The difference observed without RF spoiling (Φ BaseInc = 0 • ) was 3.0%.
The biases observed without RF spoiling and with BaseInc = 120 • relative to the sequential case were removed when the GLM method was used to process the same interleaved data ( Figure 5B; 0.21% (0.36%) for no RF spoiling and 0.13% (0.44%) for RF spoiling with BaseInc = 120 • ).

| In vivo experiments
The reference acquisition, which used interleaved ordering, a longer TR of 100 ms and an RF spoiling increment of BaseInc = 120 • , produced consistent B + 1 maps with the Classic and GLM methods (see Figure 7C,F,I). The median (interquartile range) differences relative to the GLM method were: 0.09% (0.72%), 0.30% (0.97%) and 0.11% (0.60%) for participants 1, 2, and 3, respectively. However, higher B + 1 values were observed in the ventricles of Participant 2 (see Figure 8) of the map computed with the classic method. Consistent with this observa tion, a tail in the histogram was also present for participant 2 (see Figure 7F).
Large bias was seen with respect to this reference when B + 1 maps were estimated from interleaved data, acquired without RF spoiling, using the Classic method (see Figure 7A,D,G). Qualitatively, the bias was highly visible in the B + 1 maps (see Figure 8A) following anatomical detail, and was greatest in the ventricles with long T 1 and T 2 . In the difference map, strong bias was visible along the cortical ribbon (see Figure 8B). Median (interquartile range) differences were 4.52% (2.92%), 3.99%

F I G U R E 5 Histograms of the difference in B +
1 measured in a phantom with an interleaved acquisition scheme relative to the reference B + 1 map acquired with a sequential acquisition and processed using the Classic method. The B + 1 maps were calculated using either the Classic (A) or GLM (B) method. The data were acquired with a short TR (35 ms) without (purple) or with RF spoiling (Φ BaseInc = 120 • (blue) and Φ BaseInc = 90 • (yellow) ), or a long TR (100 ms) and RF spoiling (Φ BaseInc = 120 • ) (red). The reference B + 1 map, with respect to which the error was calculated, was acquired with sequential ordering and RF spoiling (Φ BaseInc = 90 • ) and processed using the Classic method (2.48%), 3.86% (2.63%) for participants 1, 2, and 3, respectively. These biases were greatly reduced when the data were processed using the GLM method: −0.33% (1.78%), −0.10% (1.22%), and −0.35% (1.98%), respectively (see Figures 7 B,E,H and 8C,D).
High variability in B + 1 bias was observed when the sequential acquisition ordering was used (see Figure 7, blue curves) and artefacts were visible in the B + 1 and difference maps (see Figure 8A-D). This was the case regardless of the processing approach.

| DISCUSSION
Efficient methods for mapping the B + 1 transmit field with high accuracy and precision are prerequisite for demanding MRI applications, such as the quantification of the longitudinal relaxation rate. 3 Biases in B + 1 estimates may underlie intersite differences in relaxation rates, 24 while uncertainty in the estimates will lower reproducibility. 4 B + 1 mapping based on the phase accrued due to the BSS has been reported to be an efficient technique for accurately estimating the spatial distribution of B + 1 when compared with other magnitude or phase-based techniques. 7,25 Our numerical simulations indicate that the sequential approach for acquiring the necessary BSS data will deliver a bias-free estimate of the B + 1 field. However, in agreement with previous reports, 10,17 our in vivo experiments show that this approach is sensitive to phase perturbations over time such as those caused by motion and scanner drifts. The resulting B + 1 maps had visible artefacts and large biases (see Figure 8, right column). The greater robustness of the interleaved acquisition scheme has led to its adoption in more recent work using this technique. 11,26 In the interleaved case, however, our numerical simulations showed that the phase never reached steady-state but rather a pseudo-steady-state that alternated between two conditions depending on the frequency of the preceding offresonance BS pulse. As a result, the difference in phase between interleaves is not solely due to the BSS effect and, therefore, does not match the phase difference of the sequential ordering scheme (Figure 2). The additional phase accrued biased the estimated BSS phase (Figure 3) and, therefore, the B + 1 estimates in phantom and in vivo experiments. We have shown that the bias depends on intrinsic tissue properties (T 1 , T 2 ) as well as sequence parameters (TR, RF spoiling increment, amplitude of the BS pulse).
Although the biases observed here are relatively small, their impact on the R 1 estimate can be far greater. For example, the variable flip angle technique, widely used for whole brain R 1 mapping 27-30 is highly sensitive to B + 1 inhomogeneity and, therefore, requires correction. In this case, when estimating R 1 from 2 weighted images with different flip angles, the accuracy of the B + 1 estimate is crucial because it can be shown that a given bias in B + 1 will lead to a bias in the R 1 estimate that is at least twice as large and increases by even more for acquisitions with high flip angles or large error in the B + 1 estimate. In fact, under certain conditions, such as in cerebrospinal fluid where the T 1 and T 2 are long, the B + 1 bias can reach 8% (see Figure 8) when no RF spoiling is used, which would lead to a minimum of 16% bias in R 1 with the variable flip angle technique. Hence, even small errors must be accounted for if accurate and robust R 1 estimates are to be obtained.
Here, we have proposed and validated a novel acquisition and processing scheme for interleaved BSS-based B + 1 mapping that does not suffer from these biases. Crucially, multiple echoes are acquired either side of the BS pulse and a GLM framework is used to describe the phase evolution over time. The GLM models the effects of the BS pulse, B 0 inhomogeneity, eddy currents and phase offsets both common to, and specific to, the positive and negative off-resonance frequency interleaves. The bias observed with the Classic method results from the invalid assumption that the only difference between the 2 interleaves is the phase imparted by the BS pulse. Indeed, as demonstrated by the numerical simulations, a difference is already present before playing out the BS pulse ( Figure 2D). The use of 2 distinct regressors (X + offset and X − offset ) in the model matrix of the GLM allows the 2 interleaves to differ, even before the BS pulse. This removes the bias in B + 1 that would otherwise be present. Numerical simulations ( Figure 3) and phantom experiments (Figure 4) confirm this, with both showing peak differences for RF spoiling increments of 0° (equivalent to no RF spoiling), 60°, 120°, and 180°.
Inversion recovery and multi-echo spin echo experiments indicate T 1 and T 2 times of 550 ms and 70 ms for the FBIRN phantom 23 used. However, the latter estimation did not incorporate any correction for stimulated echoes, 31 and the T 2 may be as short as 50 ms. Simulations using the same sequence parameters as the phantom experiments, with a T 1 of 550 ms, and a T 2 of 70 ms, predicted an error of 4.9% ( Figure 3D) for the case of no RF spoiling. With a shorter T 2 of 50 ms, a lower error of 3.4% was predicted by the simulations (data not shown). These results are in broad agreement with the somewhat lower error of 3.0% observed experimentally for this case. Of note, incorporation of diffusion effects into the simulations 32,33 had little impact on the level of bias in the estimated B + 1 . Furthermore, while it is the phase component that is key to estimating B + 1 , it is also worth noting that our simulations predicted a difference in the magnitude of the magnetization between TRs with interleaved off-resonance frequencies, and that this difference would depend on the RF spoiling increment. Good agreement was again seen between prediction (see Figure 8A) and experiment (see Figure 8B). No such difference was predicted for sequential ordering of the offresonance frequencies. Also in agreement with the numerical simulations, the proposed GLM method removed the dependence of the B + 1 estimates on the RF spoiling increment, the TR, and the acquisition mode in both phantom ( Figure 5) and in vivo (Figures 6 and 7) experiments.
The robustness of the GLM method to the sequence parameters, and the RF spoiling increment in particular, makes this method more flexible, which can be exploited to optimize the signal-to-noise ratio. Given the dependence of the signal amplitude on the RF spoiling increment, 34 a small gain in reproducibility can be expected by choosing the optimal value. In fact, theoretical analysis of the variance of the B + 1 estimates can be used to show that the GLM should deliver higher precision. This has been verified empirically (data not shown) when using the same data for each processing method, as has been done for all of the experiments presented in this work. Although this does not affect the accuracy, it does penalize the Classic approach from a precision perspective because the TE is longer than necessary. Nonetheless, theoretical analysis would also predict improved precision with the GLM when compared with the Classic approach even with an optimal, shorter, TE. Determining the sequence settings that maximize the reproducibility and quantifying the full benefit that can be gained empirically will be the focus of future work.
In theory, the GLM method could use just a single off-resonance frequency with a reduced model matrix containing only half-length regressors for X BSS , X B 0 , X Even∕Odd and X +

Offset
. 35 In this case, X BSS would model all the phase imparted by the BS