Robust 3D Bloch-Siegert based B1+ mapping using Multi-Echo General Linear Modelling

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 Shift (BSS) of two acquisitions with opposite off-resonance frequency RF pulses. Interleaving these RF pulses ensures robustness to motion and scanner drifts, however, here we demonstrate that doing so also introduces a bias in the B1+ estimates. Methods It is shown here via simulation and experiments that the amplitude of the error depends on MR pulse sequence parameters, such as TR and RF 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 (GLM) to estimate the correct BSS-induced phase is presented. Results In line with simulations, phantom and in-vivo experiments confirmed that the GLM-based method removed the dependency on tissue properties and pulse sequence settings. Conclusion The GLM-based method is recommended as a more accurate approach to BSS-based B1+ mapping.


Introduction
Knowledge of the spatial distribution of the radiofrequency transmit field (B1 + ) is crucial to many MRI applications. Moderate accuracy may suffice when setting transmitter gain 1 or calibrating multi-channel systems 2 . However, very high accuracy and precision are required for many quantitative MRI applications, e.g. mapping the longitudinal relaxation rate to characterise cortical myelination 3,4 .
Phase-based methods may be preferred as they are theoretically insensitive to T1 relaxation effects which often bias magnitude-based methods, especially at short TR. In the Bloch-Siegert (BS) 5,6 approach, an off-resonance RF pulse leads to the Bloch-Siegert frequency shift (BSS), and an associated phase accumulation, which is proportional to the square of the pulse amplitude thereby encoding the B1 + field. This technique performed favourably in a recent review of the accuracy, precision and practicality of a range of prominent B1 + mapping techniques 7 , and has been shown to be less sensitive to 0 inhomogeneities and chemical shifts 8 than other phase-based methods.
By subtracting the two phase images, the BSS effect is enhanced and unrelated phase components are removed, e.g. phase accumulated across echo time (TE) due to B0 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 B0 inhomogeneity on the BSS, up to first order 9 . These two 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 fashion disturbs the steady-state and introduces an additional phase difference between the two acquisitions, especially at short TR. We show that this additional phase difference leads to bias in the B1 + map that depends on the relaxation parameters of the studied tissue, the specifics of the RF spoiling regime and the actual B1 + amplitude. We additionally propose and validate a modified BSSbased 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 modelling 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 SAR or reducing the sensitivity.

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

The classic method: isolating the BSS phase by subtraction
The classic approach to BSS-based B1 + mapping consists of acquiring two datasets with BS pulses of opposite off-resonance frequency (i.e. + and − ). The phase difference between these is: Since the first order terms that depend on Δ 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 B0 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 B1 + estimate will be erroneous.

The GLM method: isolating the BSS phase by modelling a multi-echo acquisition
We propose an alternative to the classic BSS approach that computes accurate B1 + maps even if conditions vary between the two off-resonance frequency acquisitions. This approach relies on two novel features: 1. A dual-offset multi-echo sequence: In the modified BSS-based B1 + mapping sequence (Fig.1), multiple echoes are acquired after one excitation pulse. Two echoes, after the BS pulse, have previously been used 18 to concurrently compute the B0 inhomogeneity, whereas here 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.

A general linear model:
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.

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 B1 + 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 .
The impact of the RF spoiling increment was investigated by changing 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 was alternated between successive BS pulses. For sequential acquisitions, the sign was switched after half the total number of pulses ( 2 ), at which point RF spoiling was also reset. The estimated B1 + amplitude was calculated as follows having isolated the first term (only) of Eq.1: For the Classic method, Φ was taken to be half the difference in phase accumulated after the BS pulse of the acquisitions with opposite off-resonance frequency, and termed Φ (see Supporting Information).
For the GLM method, Φ was taken to be the mean absolute phase accrued during the BS pulses of these, and termed Φ (see Supporting Information).

MRI measurements  MR pulse sequence
Measurements were performed at 3T (Siemens Prisma) using a body coil for transmission and a 32channel head coil for signal reception using an in-house MR pulse sequence (Fig.1). A Fermi pulse of duration T=2ms imparted the Bloch-Siegert shift after the second echo. Acquiring two echoes before the BS pulse served to minimize the correlation of the 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 prior to the BS pulse, to ensure the same dephasing state for the magnetisation 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 PE1), to crush any undesired onresonance excitation and to minimise 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 non-ideal conditions is reduced by using a high crusher moment since 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 T1 of 1000ms, was chosen to maximise the precision of the B1 + 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 off-resonance 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.

 B1 + map estimation
All data, including B1 + 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 B1 + 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: Φ = Φ /2. 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 (Fig.1). The BSS phase Φ was captured by the first regressor of the model matrix such For both methods, 1 was computed on a voxel-wise basis using Eq.4.

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

Phantom Experiment 1: Comparing processing approaches
This experiment probed the impact of the processing method as a function of the RF spoiling increment, Two B1 + maps were estimated for each scan using the Classic and GLM methods respectively. The percent difference in B1 + was calculated with respect to the reference map for the same processing method. To compare processing methods, the percent difference between the two B1 + maps derived from the reference acquisition was also computed, with respect to that obtained with the Classic method. Table 1 lists sequence parameters of all experiments.

Effect of acquisition order and RF spoiling
For sequential acquisition ordering, once steady-state was reached for each off-resonance frequency blocks (Fig.2 a), the phase difference between these was zero before the BS pulse (Fig.2 c). It was nonzero after the BS pulse (Fig.2 e) since 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 (Fig.2 b). This temporal variance caused a non-zero phase difference between the interleaved acquisitions with opposite BS pulse frequencies, both before (

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 (Fig.3). The greatest errors were predicted for phase increments of 0° (equivalent to no RF spoiling) and 180°. Large errors were also predicted for phase increments of 60° and 120° whereas no error was predicted at 90°.

Effect of TR, B1 + , T2 and T1.
Increasing TR greatly reduced the predicted error in Φ (Fig.3 a). Longer T1 times led to larger predicted errors, but had less impact than TR, (Fig.3 c). Similarly, longer T2 times predicted larger errors (Fig.3 d), 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 (Fig.3 b). Note that in this case 1 + also increases, so while the predicted relative error decreased, the absolute error is actually increased.

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

Phantom experiments
Phantom Experiment 1: Comparing processing approaches With RF spoiling increments of 0 and 180°, the B1 + estimates were very different depending on the processing method used (Fig.4), but these estimates agreed when the RF spoiling increment was 90°.
The biases observed without RF spoiling and with = 120 ∘ relative to the sequential case were removed when the GLM method was used to process the same interleaved data (Fig.5 b; 0.21% (0.36%) for no RF spoiling and 0.13% (0.44%) for RF spoiling with = 120 ∘ ). Figure 5: Histograms of the difference in B1 + measured in a phantom with an interleaved acquisition scheme relative to the reference B1 + map acquired with a sequential acquisition and processed using the Classic method. The B1 + maps were calculated using either the Classic (a) or GLM (b) method. The data were acquired with a short TR (35ms) without (purple) or with RF spoiling ( = 120 ∘ ( ) = 90 ∘ ( )), or a long TR (100ms) and RF spoiling ( = 90 ∘ ) ( ). The reference B1 + map, with respect to which the error was calculated, was acquired with sequential ordering and RF spoiling ( = 90 ∘ ) and processed using the Classic method.

In vivo experiments
The reference acquisition, which used interleaved ordering, a longer TR of 100ms and an RF spoiling increment of = 120 ∘ , produced consistent B1 + maps with the Classic and GLM methods (Fig.7 c, f, i). The median (IQR) 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 B1 + values were observed in the ventricles of Participant 2 (Fig.8) of the map computed with the classic method. Consistent with this observation, a tail in the histogram was also present for Participant 2 (Fig.7 f).
High variability in B1 + bias was observed when the sequential acquisition ordering was used (Fig.7, blue curves) and artefacts were visible in the B1 + and difference maps (Fig.8 a-d). This was the case regardless of the processing approach.

Discussion
Efficient methods for mapping the B1 + 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 B1 + estimates may underlie inter-site differences in relaxation rates 24 while uncertainty in the estimates will lower reproducibility 4 . B1 + mapping based on the phase accrued due to the Bloch-Siegert shift has been reported to be an efficient technique for accurately estimating the spatial distribution of B1 + when compared to 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 B1 + field. However, in agreement with previous reports 17,10 , 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 B1 + maps had visible artefacts and large biases (Fig.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 steadystate but rather a pseudo-steady-state that alternated between two conditions depending on the frequency of the preceding off-resonance 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 (Fig.2). The additional phase accrued biased the estimated BSS phase (Fig.3) and therefore the B1 + estimates in phantom and in vivo experiments. We have shown that the bias depends on intrinsic tissue properties (T1, T2) 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 R1 estimate can be far greater.
For example, the Variable Flip Angle (VFA) technique, widely used for whole brain R1 mapping 27-30 is highly sensitive to B1 + inhomogeneity and therefore requires correction. In this case, the accuracy of the B1 + estimate is crucial since it can be shown that a given bias in B1 + will lead to a bias in the R1 estimate that is at least twice as large and increases by even more for acquisitions with high flip angles or large error in the B1 + estimate. In fact under certain conditions, such as in CSF where the T1 and T2 are long, the bias can reach 8% (Fig.8) when no RF spoiling is used, which would lead to a minimum of 16% bias in R1 with the VFA technique. Hence even small errors must be accounted for if accurate and robust R1 estimates are to be obtained.
Here, we have proposed and validated a novel acquisition and processing scheme for interleaved BSSbased B1 + 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, B0 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 two 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 (Fig.2 d). The use of two distinct regressors ( + and − ) in the model matrix of the GLM allows the two interleaves to differ, even prior to the BS pulse. This removes the bias in B1 + that would otherwise be present. Numerical simulations (Fig.3) and phantom experiments (Fig.4) 23 used. However, the latter estimation did not incorporate any correction for stimulated echoes 31 , and the T2 may be as short as 50 ms. Simulations using the same sequence parameters as the phantom experiments, with a T1 of 550 ms, and a T2 of 70 ms, predicted an error of 4.9% (Fig. 3d) for the case of no RF spoiling. With a shorter T2 of 50ms, 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 B1 + . Furthermore, while it is the phase component that is key to estimating B1 + , it is also worth noting that our simulations predicted a difference in the magnitude of the magnetisation 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 (Fig.8a) and experiment (Fig.8b). No such difference was predicted for sequential ordering of the off-resonance frequencies. Also in agreement with the numerical simulations, the proposed GLM method removed the dependence of the B1 + estimates on the RF spoiling increment, the TR and the acquisition mode in both phantom (Fig.5) and in vivo (Fig.6 and Fig.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 (SNR). 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 B1 + 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 penalise the Classic approach from a precision perspective since the TE is longer than necessary. Nonetheless, theoretical analysis would also predict improved precision via the GLM when compared to the Classic approach even with an optimal, shorter, TE. Determining the sequence settings that maximise 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 inhomogeneity and chemical shift. However this could be corrected with the information captured by the second regressor 0 . However, since the BSS phase is imparted only once, the BS pulse flip angle would need to be doubled to achieve the same phase-to-noise ratio. This is problematic from a SAR perspective, and the benefit of a single off-resonance frequency acquisition would be negated if the TR were also doubled to address it. Besides, removing the second acquisition with opposite off-resonance frequency prevents the isolation of the BS phase of interest because any phase caused by the eddy currents of the crushers, for example, would also be captured by the same regressor making the problem ill-posed. A workaround consisting of adding further gradients after the 4 th echo, to distinctly capture the effects of eddy currents, has been proposed and tested but has proven to be effective only in phantom experiments 32 .

Limitations:
Given that the GLM method relies on a multi-echo sequence, additional pre-processing steps are required compared to the Classic approach. Phase unwrapping across echoes is necessary. It has been necessary to spatially unwrap the phase difference between successive echoes to deal with large phase accumulation between successive echoes, then cumulatively add these to the first echo.
In the proposed method, multiple echoes are used to estimate the BSS phase, some of the echoes may suffer from dropout and potentially introduce noise into the estimate. To minimize this effect a weighted least-square (WLS) approach had been used to estimate the parameters of the GLM, down-weighting echoes with lower magnitude.
Conventionally, the BS pulse is applied just after the excitation pulse. Here, two echoes preceded the BS pulse and the difference of the third echoes, from the different off-resonance frequency acquisitions, was used to estimate B1 + using the Classic method. This increases the minimum TE (by ~4ms) and therefore lowers the SNR relative to the single-echo method. However, while this might reduce precision, it would not be expected to introduce bias.
This study focused on short TR 3D acquisitions. For 2D acquisitions, the shot-to-shot inconsistencies may be less problematic since the TR will be longer, although a bias was still observed in long T1 regions with a TR of 100 ms (Fig.8-First column). Besides,, the two acquisitions of one slice will be more separated in time which may result in additional phase differences due to motion, similar to the problem affecting sequential acquisitions.
Although more efficient pulses have been proposed 18,36 , only the commonly used Fermi shape for the BS pulse was investigated here. However it can be shown that for the same imparted Bloch-Siegert phase and the same off-resonance frequency, the bias introduced by the interleaved acquisition order is equivalent for a Fermi pulse and the more optimised pulse design suggested by Duan et al 18 .
While we have shown how to more accurately estimate the BSS phase, the conversion to 1 may still be a source of inaccuracy if any assumptions underlying Eq.4 are violated 18 . For the particular conditions we have explored (a Fermi pulse with 2ms duration, 1 / = 0.23 and 1 = 11 the error from this approximation is estimated from simulation to be less than 1%. Regardless of how the phase is converted to a B1 + value, it is imperative that the bias caused by interleaving the off-resonance pulses be removed.
The precision and accuracy of the GLM technique and the Classic method with ϕ = 90 ∘ relative to other B1 + mapping methods remain to be investigated. However, determining absolute accuracy will always be challenging since every method will have its own limitations.

Conclusion
Interleaved acquisitions are recommended for Bloch-Siegert based B1 + mapping to increase robustness to motion and scanner drift. However we have shown that with the Classic estimation method, this can introduce error into the B1 + estimates that will depend on tissue properties and sequence settings. In theory, one could use an RF spoiling increment of 90° to be immune to this error. However, we have also proposed and validated a multi-echo sequence design, combined with a GLM framework, to robustly isolate the BSS-induced phase regardless of the sequence parameters used. This allows bias free, low error estimates of the B1 + efficiency that do not depend on tissue properties, sequence settings and would furthermore be immune to reproducible hardware imperfections. Importantly, the proposed modifications do not extend acquisition time, reduce sensitivity or increase SAR. The latter is particularly important since SAR is a limiting factor at higher field strengths.

Supporting Information
Bloch-Siegert phase shift expression The first-order Taylor  In a static magnetic field, 0 , the precessional frequency is the Larmor frequency given by 0 = 0 .
The BSS approach applies an off-resonance RF pulse with a rotating frequency of 2 , such that = 0 − 2 ≠ 0. The off-resonance pulse (termed the BS pulse) modifies the precessional frequency during the time of its application. It has an amplitude 1 such that 1 = 1 .
In the frame of reference rotating at 0 , the precessional frequency during this pulse changes to: If the off-resonance frequency of the BS pulse is much larger than its strength in frequency units, i.e. ≫ 1 , the expression reduces to: Given that the precessional frequency is no longer equal to 0 during the BS pulse, an additional phase is accumulated. Since it is proportional to 1 2 this phase can be used to map the spatial distribution of B1 + .
Allowing for local 0 inhomogeneity, the local precessional frequency will be 0 + Δ 0 . In the rotating frame of reference and under the similar condition as before, ( + Δ 0 ) ≫ 1 , the precessional frequency during the BS pulse becomes: The first term on the right hand side of this expression, representing the frequency shift due to the BS pulse, will be called . Note that this effect adds to the phase accrued due to local 0 inhomogeneity.
The phase specifically accumulated due to the BS pulse, of arbitrary shape and duration T, is given by: where 1 is the peak amplitude of the BS pulse, and 1 is its normalized shape. Assuming |Δ 0 | ≪ | |, the first order Taylor expansion of this expression is:

Numerical simulations
Numerical simulations were carried out to evaluate the phase difference between two acquisitions with opposite off-resonance frequencies both before and after the BS pulse. A typical gradient echo acquisition was simulated in MATLAB (The MathWorks, Inc., Natick, MA) by a series of matrix operations as described below: , was then simulated as a series of small rotations indexed by ∈ [1: ]. Each small rotation k is a rotation of an angle | 1 ( * )| * about an axis defined by an azimuthal angle of * * + .
5-A rotation of the resulting magnetization vector about the z-axis by an angle -Ω 1 is then applied to simulate the crusher of opposite polarity. The resulting phase at this time is called + if the frequency of the preceding BS pulse was positive, or − if the frequency was negative. For simplicity, the pulse and crushers were assumed to be instantaneous from a relaxation standpoint such that Φ and Φ have the same TE. Note that the relaxation is instead accounted for in step 7. 6-Following the second crusher gradient, a spoiler gradient was applied such that the magnetization vector from step 5 was rotated around the z-axis by an angle of Ω 2 . 7-T1 and T2 relaxation were then simulated for the remainder of the TR period, i.e. TR-TE. 8-When RF spoiling was simulated, the phase of the RF pulses, , was incremented such that = + and = ϕ + ϕ .
These steps were repeated times, by updating the magnetization of Step 1 with the magnetization vector resulting from Step 7.
The simulation was performed for a 2D grid of spin ensembles ( * ) with varying values for Ω 1 and Ω 2 ranging from 0 ∘ to Ω 1 and Ω 2 , respectively. The values of the crusher dephasing Ω 1 varied along one dimension of the grid, whereas the values of the spoiler dephasing Ω 2 varied in the other direction to simulate crushers and spoiler along orthogonal axes. The integral of these spins gave the magnetisation in a single voxel.
For interleaved ordering, the sign of was switched before each BS pulse. For sequential ordering the sign was only switched after half the total number of pulses ( 2 ). If RF spoiling was included in the simulation, the RF spoiling phase (ϕ and ) were reset to 0 when the sign of was switched in the sequential case. In the case of an interleaved acquisition, the phase increment occurred for each repetition, including when the off-resonance frequency also changes.
The numerical values of the parameters used in the simulations are listed in the Supporting Information Table S1 Supporting Information Table S1: Parameters used in the numerical simulations The BS phase estimated by the Classic approach is simulated as: The BS phase estimated by the GLM approach is simulated as: