Magnetization transfer and frequency distribution effects in the SSFP ellipse.

PURPOSE: To demonstrate that quantitative magnetization transfer (qMT) parameters can be extracted from steady-state free-precession (SSFP) data with no external T1 map or banding artifacts.
METHODS: SSFP images with multiple MT weightings were acquired and qMT parameters fitted with a two-stage elliptical signal model.
RESULTS: Monte Carlo simulations and data from a 3T scanner indicated that most qMT parameters could be recovered with reasonable accuracy. Systematic deviations from theory were observed in white matter, consistent with previous literature on frequency distribution effects.
CONCLUSIONS: qMT parameters can be extracted from SSFP data alone, in a manner robust to banding artifacts, despite several confounds.

INTRODUCTION

Quantitative magnetization transfer (qMT) imaging is of great interest for investigating multiple neuropathologies, particularly those relating to myelin.1, 2 The steady‐state free‐precession (SSFP) sequence is sensitive to MT and its short acquisition time is an advantage compared to more traditional gradient echo sequences.2, 4 qMT parameters can be found from SSFP data acquired at multiple flip angles,5 but this requires an additional T 1 map and assumes that all voxels are on‐resonance. This is valid at 1.5T but not at 3T or above where susceptibility gradients produce banding artifacts.6

The PLANET method can calculate T 1 & T 2 from complex‐valued SSFP data at a single flip angle but multiple phase increments, provided TR is sufficiently long for MT effects to be insignificant.7, 8 We show that by combining these methods qMT parameters can be calculated from SSFP data at multiple flip angles and phase increments. However, we demonstrate that frequency distribution effects in white matter, where multiple pools of water protons exist, present a significant confound.8, 9, 10

THEORY

We incorporated MT effects into the SSFP signal ellipse by following5 but including off‐resonance effects. We make the derivation in terms of the bound pool fraction f  = M 0/(M 0 + M 0) where the subscripts f & b denote the free and bound pools respectively. Hence the magnetization vector is as the bound pool is not visible during the MR experiment. We then define: R represents relaxation, with E 1 =  exp (−TR/T1), E 1 =  exp (−TR/T1) and E 2 =  exp (−TR/T2). X represents exchange, which is governed by f , f , the forward (bound‐to‐free) exchange rate k , the reverse (free‐to‐bound) exchange rate k and the exponential average of the exchange rates E  =  exp (−(k + k )TR). O describes accrued off‐resonance phase ψ = 2πΔf 0 TR.

P represents the RF pulse, including excitation of the free pool by the angle α, rotation by a phase increment ϕ , and saturation of the bound pool by the factor E  =  exp (−WT ). T is the length of the RF pulse, and W is the saturation rate: G 0 is the value of the lineshape at resonance, p 1 is the ratio of the pulse's mean B1 to the maximum B1, and p 2 is the ratio of the mean B12 to the maximum B12.11, 12 The factors p 1 and p 2 describe the shape of the RF pulse relative to a block pulse, and provide a simple way to scale the saturation rate for different flip angles and pulse lengths (the saturation rate of a block pulse simply scales with flip angle squared and the inverse square of pulse‐length). We assumed a fixed value of G 0 = 14 μs, which is the extrapolation of the Super‐Lorentzian from an off‐resonance value of 1 kHz to zero.5, 13 Strictly, the lineshape should be evaluated at the local off‐resonance value; however at 3T, the highest off‐resonance values in the human brain are approximately 100 Hz which is well below the 1 kHz value at which the Super‐Lorentzian is assumed to be flat. Finally, we define the following useful terms, similar to Equation 11 of Gloor et al5:We then define the magnetizations before and after the RF pulse asand then calculate the steady‐state solution at the echo time TE = TR/2:where the term represents relaxation and accrued off‐resonance between the RF pulse and TE. The complex‐valued steady‐state magnetization is then:where θ = ψ + ϕ . This can be converted to the elliptical form6 where

In the definition of M , we have included the decay of transverse magnetization at TE via the term, but off‐resonance in the same time period is incorporated into Equation 7 so that M remains a real quantity. Equation 7 also contains the B 1+ RF phase term ϕ , which was neglected earlier for simplicity.7

This equation has the same form as the elliptical expression of Xiang and Hoff, where M is the geometric solution, which is close to the geometric center of the ellipse.6 a, combined with M , determines the ellipse size, and is purely affected by T 2 and no MT parameters. b describes how flattened the ellipse is. In contrast to PLANET where M 0, T 1, and T 2 can be algebraically calculated from M , a, and b, it is not possible to extract M 0, T 1, f , and k from a single ellipse, principally because only two of the ellipse parameters (M and b) depend on MT. However, by acquiring multiple ellipses with different MT weightings, it is possible to recover the MT parameters.

METHODS

Model fitting

The following assumes multiple ellipses labeled i = 1 … m acquired with phase increments j = 1 … n. Each ellipse has a different α or T to change the MT weighting. Fitting proceeds in two steps. We first calculate the ellipse parameters for each ellipse separately, resulting in m values each of M , a and b . In the second step, the MT parameters are found from the ellipse parameters.

PLANET calculates M , a, and b algebraically for a single ellipse, but is difficult to use in the region near the Ernst angle because of a sign change in the formula for b.7 At 3T and above, due to B 1 + inhomogeneity, it cannot be guaranteed that α will be above the Ernst angle in all voxels. Also, although numerous algebraic ellipse methods exist, they have been developed for the general case where no a priori information is available,14, 15 but here the ellipse must be vertical (b < 2a/(1 + a 2)) and must not contain the origin.6 Hence in preference, we used a bounded nonlinear least‐squares fit of Equation 7 to find M , a, and b for each ellipse, implemented in our freely available toolbox.16 The cost function waswhere s is the acquired complex‐valued SSFP data for ϕ and m is the corresponding value calculated from Equation 7. The bounds were 0 < a <  exp (−TR/5000), 0 < b < 1, −2π < ψ < 2π, and −2π < ϕ  < 2π, and we forced the ellipse to remain vertical. The initial guess for M was |∑s /n|, while a and b were calculated assuming T 1/T 2 = 1000/50 ms for the specified TR and α. The initial guesses for the tightly coupled ψ and ϕ were more complex. Three guesses for ψ (−π, 0, +π) were tried, and for each ϕ was set to the phase of the complex mean of the data minus the guess for ψ. The guess with the lowest root‐mean‐square residual to the data was used as the initial guess for the full nonlinear least squares algorithm. A Huber Loss function was applied to the residuals to suppress the influence of outliers.17

From the set of a values, T 2 can then be found directly:To find the MT parameters, a nonlinear least squares fit was used with the cost‐function:where and are found from Equations (8a) and (8c), respectively. Values of M were scaled by their mean to make them comparable in magnitude to b. An initial guess of M 0 = 13, f  = 5%, k  = 2 s−1 and T 1 = 1000 ms was empirically determined to converge for a majority of voxels in the brain. T 1 was set equal to T 1.5

Simulations

Monte Carlo simulations were conducted to validate our method. The tissue parameters were M 0 = 100, T 1/T 2 = 1100/80 ms, k  = 4.45 s−1, and f  = 10%, and the sequence parameters were α  = 15, 15, 30, 30∘, T  = 256, 1024, 256, 1024 μs, and TR  = 7, 7.768, 7, 7.768 ms, producing four MT weightings. The sequence parameters correspond to our in vivo protocol which was dictated by system limits that are discussed below.

M , a , and b were then calculated with Equation 8 and combined with random values of ψ and ϕ drawn from uniform distributions between −π and π. m was calculated with Equation 7 and ϕ  = j60∘, j = 1 … 6. This gives 24 complex datapoints (four ellipses with six points each) to which complex‐valued Gaussian noise with σ = M 0/500 was added. This value corresponds to an SNR (defined as mean signal magnitude divided by noise) of approximately 70, which matches literature values for gradient echo sequences in white matter.184,096 realizations of noise were used.

In vivo measurements

Images from a healthy volunteers were acquired on a GE MR750 3T scanner equipped with 50 mTm−1 gradients and a 32‐channel head coil (Nova Medical) in accordance with local ethics procedures. The limiting factors for protocol design were the system SAR limit and eddy current artifacts,19 of which the latter were the dominant factor (see discussion). A 3D acquistion with 210 × 210 × 180 mm field‐of‐view and isotropic 1.5 mm voxel size (matrix size 140 × 140 × 120) was used. All phase increments for each ellipse were acquired sequentially before moving to the next ellipse, and 2 seconds of dummy TRs were added at the start of each volume acquisition to establish a steady state. The readout bandwidth was set to 25 kHz and the manufacturer's option to de‐rate the maximum gradient slew rate was enabled, resulting in TR = 7 ms for the 256 μs pulse and TR = 7.768 ms for the 1024 μs pulse. Partial k‐space acquisition with NEX = 0.75 (elliptical k‐space coverage) and parallel imaging (ASSET = 3) was used for a total acquisition time of 17 minutes. Other parameters were the same as for the simulations above. A map was acquired using the manufacturer's Bloch‐Siegert sequence to correct for RF inhomogeneity.20, 21 All images were motion corrected before further processing.22

RESULTS

Figure 1 shows the results of the Monte Carlo simulations for the ellipse parameters. Four histograms are present in each subfigure, one for each of the simulated ellipses, and the true values are shown by dotted lines. Parameters M , a, and b all have very low bias of less than 1% and coefficient of variation (CoV) of less than 2% for all the ellipses. Note that as the parameter a only depends on T 2 and TR, there are only two true values and the histograms for 256 μs and 1024 μs pulse‐lengths overlap.

Figure 1

Results of fitting the ellipse parameters M , a, b to the Monte Carlo simulations. The parameters for each individual ellipse are shown in the colors given in the legend, and true values are shown with dotted lines

Because the true values of θ 0 and ϕ were chosen randomly (to assess the robustness of the method to off‐resonance and effects), the true value modulus 2π was subtracted from the fitted value before plotting, and hence the histograms are centered on zero. Again, the fit is excellent with negligible bias and standard deviations of less than 2%. In the case of ϕ , a small number (<1%) of degenerate solutions were found at the +π/−π boundary. As ϕ is not a parameter of interest for this study, these fitting failures were considered inconsequential, and hence the range of the histogram is restricted to −2∘ to 2∘.

Figure 2 shows the results of fitting the MT parameters to the Monte Carlo simulations. The fits for M 0, f , T 1 and T 2 showed minimal bias (<1%) and reasonable CoV (<4% except for f which was 8%). k shows a skewed distribution, weighted toward higher values, and hence has a CoV of 22% and a bias of 4%.

Figure 2

Results of fitting the qMT parameters to the ellipse parameters. M 0, F, T 1,andT 2 show the expected normal distribution, but k shows a skewed distribution with a heavy tail toward high values, indicating that this parameter is difficult to fit correctly. True values are indicated with dotted lines

Experiment

Figure 3 shows the fitted ellipse parameters in the volunteer for the different MT weightings. The parameters M and b both show a reduction under increased MT weighting, while there is only a negligible change in a due to the different TR of the two settings. The off‐resonance frequency (calculated from θ 0) and B 1+ phase are also the same across the scans. While these parameters are not of direct interest to this study, the lack of change indicates excellent scanner stability.

Figure 3

Maps of the ellipse parameters for each MT weighting. The ellipse parameters M and b change across weightings, while a does not. Δf 0 and ϕ are not of direct interest to this study but are shown for completeness

We observed clear structure in ellipse root‐mean‐square error (RMSE), as shown in Figure 4 (first and third row). This clearly resembles white matter tracts. We plotted the ellipse in selected high‐ and low‐residual voxels (deliberately chosen to lie in opposing halves of the complex plane for clarity) to further illustrate the behavior (second and fourth row).

Figure 4

The residuals from fitting each of the four ellipses in two slices (first and third row). In each slice, a low‐ (blue) and high‐ (red) residual voxel are marked and the corresponding fitted ellipses shown (second and fourth rows). The large dots are acquired data (coloured by phase increment) and the line is the fitted ellipse. In low residual voxels, all acquired points on the ellipse, while in high‐residual voxels the ϕ  = 60∘ (light green), does not lie on the ellipse

Finally, Figure 5 shows maps of the MT parameters. The maps of M 0, f , T 1 and T 2 show good fitting quality. In contrast, the map of k displays a large number of fitting failures (bright voxels). A large bilateral ROI was drawn in WM and gave the following average values: f  = 11.3 ± 2.7%, T 1 = 794 ± 126 ms, T 2 = 54 ± 19 ms, and finally k  = 5.0 ± 2.4 s−1.

Figure 5

Maps of the fitted quantitative MT parameters through four slices of the brain. The exchange parameter k shows a large number of fitting failures where the value has reach the upper fitting bound; however, this does not appear to impact the fitting of the other parameters

DISCUSSION

The results above demonstrate that it is possible to extract T 1, T 2, proton density, bound‐pool fraction, and exchange rate information from SSFP data using an elliptical signal model, with no additional measurement of the apparent T 1. A B 1 + map is still required to calculate the local flip angle. Several confounds and difficulties present themselves with the current approach—in particular, the frequency distribution of white matter protons, the lengthy acquisition time, lack of precision in quantifying the exchange rate, and an inability to probe the bound‐pool lineshape.

The most significant issue with the current work is the increased residual for fitting the ellipses in white matter, which indicates that the above model does not fully explain the behavior of the SSFP signal in such regions. White matter is known to exhibit multiple water components with different T 2, which are generally attributed to free water and water trapped in myelin bilayers.18, 23 Indeed, Shcherbakova et al discussed how a two component system could be modeled with the ellipse formalism, but importantly only showed simulation results and did not discuss exchange effects between the pools.8 They suggested that within typical ranges of TR and flip angle the expected impact of the frequency distribution on PLANET is minimal. Similarly, our observed deviations from the ellipse appear small. Although the RMSE can approach 15% of M in WM, the majority of points exhibited an excellent fit and only a single outlier phase increment was observed in most voxels. The exact phase increment varied with off‐resonance frequency across the brain.

Previously, Miller et al discussed at length the asymmetries of the SSFP signal profile in white matter, but did not use the ellipse formalism.9, 10 We argue that the residual maps shown in Figure 3 are similar to Figure 3 of,10 and that if translated to the complex plane the asymmetry profiles would resemble the ellipses of Figure 4. Hence we conclude that the increased residuals displayed for fitting the ellipse in white matter are consistent with previous literature and can be attributed to the frequency distribution of protons in the myelin water pool. Quantifying the impact of the frequency distribution on the subsequent MT quantification would require modeling both MT and the frequency distribution, which is beyond the scope of this work.

The acquisition time and exchange‐rate issues are linked. There are two reasons for the relatively lengthy acquisition time described here, relative to the modest voxel size: the number of acquired images and the long TR. Our protocol required 24 separate SSFP images to be acquired with different flip angles, pulse lengths and phase increments, in order to use the two‐stage fitting process where we first fit the ellipse parameters and then the MT parameters. However, there are only seven parameters in the final set (M 0, f , k , T 1, T 2, off‐resonance, and B 1 + phase), and so the current protocol is highly redundant. In principle, fewer SSFP images would be required if the MT‐ellipse model were fitted directly to the data, which could yield a significant speed‐up. This was beyond the scope of this paper as we wished to illustrate the ellipse behavior in WM, which would be hidden by such a model. An additional benefit to a shorter protocol would be increased robustness against subject motion and frequency drift effects.24

We found that a relatively long TR was required to avoid eddy current artifacts that affected the different phase increment images in unpredictable ways if the full gradient performance of our system was exploited. These artifacts are not usually present in traditional SSFP images acquired with a 180∘ phase increment and standard Cartesian k‐space ordering because any eddy currents cancel out in successive TRs.19 However, in this work, images are deliberately acquired with phase increments that do not lead to such cancellation (see Supporting Information Figure S1). As noted in the methods section, for the purposes of this paper, we de‐rated the gradient slew rate to remove the artifacts, which almost doubled the TR.

Shortening the TR could also improve the fitting of k , which has been shown to contain clinically useful information.25 We found that the cost function, described by Equation 11, narrows in the exchange‐rate direction when TR is shortened, implying that k would then be easier to fit (see Supporting Information Figure S2). Due to our long TR, we did not consider the effects of finite RF pulse duration on our T 2 estimates, but this may become a significant effect with shorter TR.26

In this paper, a constant value of T 2 across the brain was assumed.5, 13 However, data from classic gradient‐echo qMT experiments indicate the apparent T 2 varies across brain when using a Super‐Lorentzian lineshape.27 However, due to the use of on‐resonance RF pulses in this method, the tails of the lineshape are not queried at all and hence it is not possible to fit for T 2. It may be possible to query the bound‐pool lineshape by acquiring additional images with a multi‐band excitation pulse,28 but this was beyond the scope of this article.

CONCLUSION

We have demonstrated that MT parameters can be extracted from the SSFP signal ellipse. This method does not require an additional T 1 map as this information comes directly from the model. The method is robust against SSFP banding artifacts because off‐resonance effects are directly incorporated into the model. The acquired data showed deviations from theory in white matter which can be attributed to frequency distribution effects.

FIGURE S1 Example raw images acquired at the minimum possible TR for our system. In addition to the SSFP bands (which are not artifacts in the context of this work), there are an additional zipper‐type eddy current artifacts present when ϕ ≠ 180∘ which are marked by arrows. These artifacts change position depending on the phase increment value, and only disappeared when TR was increased significantly

FIGURE S2 Contour plots of the MT fitting cost function in the f , k plane. The red dot marks the ground truth value. Reducing TR steepens the contours in the k direction, which would make fitting the true value easier, but does not significantly change the shape of the contours in the f direction

Click here for additional data file.