Correction of Main and Transmit Magnetic Field (B₀ and B₁) Inhomogeneity Effects in Multicomponent-Driven Equilibrium Single-Pulse Observation of T₁ and T₂

General mcDESPOT Theory and Modeling the Two-Component SPGR and bSSFP Signal

Solutions to the general Bloch-McConnell magnetization expressions (15) may be derived for any pulse sequence using discrete-time formalism, where the pulse sequence is modeled as a series of matrix expressions acting upon the magnetization vector. For example, the two-component bSSFP signal may be modeled as repeated RF pulse, relaxation, and longitudinal recovery “blocks” as shown in Fig. 2. The steady-state signal is determined through iterative evaluation of these blocks until the transverse magnetization components reach a constant value.

In the case of some pulse sequences a closed form solution for the magnetization and signal can also be found by equating the magnetization immediately before and following an RF pulse, relaxation, and longitudinal recovery block and solving symbolically. The closed form expressions for the two component SPGR and bSSFP magnetization evaluated in this manner are presented in Appendix A and comprise 11 free parameters, including: the free-water T₁ and T₂ (T_1,F, T_2,F); the myelin-associated water T₁ and T₂ (T_1,M, T_2,M); the relative myelin and free-water volume fractions (f_M and f_F); the myelin and free-water proton residence times (τ_M and τ_F); the myelin and free water off-resonance (with respect to the central water peak) values (δω_M, δω_F); and a factor related to the total equilibrium longitudinal magnetization and containing additional factors including scanner amplifier gains, RF coil receive biases, etc. (ρ).

The total number of free-parameters in these models may be reduced to six through four steps/assumptions: (1) normalizing all data with respect to the mean (thereby eliminating ρ); (2) Assuming both water pools are on-resonance (i.e., δω_M = δω_F = 0); (3) assuming only two water pools (f_M + f_F = 1); and (4) assuming the water pools are in exchange equilibrium (f_M/τ_M = f_F/τ_F). Presuming correct knowledge of the flip angle (α), Eqs. (A5) and (A6), can be fit to acquired SPGR and bSSFP data to derive estimates of T_1,M, T_2,M, T_1,F, T_2,F, τ_M, and f_M. However, errors are introduced into these estimates when a deviates from the prescribed nominal value, and when the transverse phase precession (δω_M and δω_F) is nonzero. To correct for these effects, we propose to calibrate α using the DESPOT1-HIFI B₁ mapping approach (14), and account for off-resonance affects by including an additional term, δω, where δω = δω_M = δω_F, in the multiparameter fitting; coupled with acquisition of bSSFP data with more than one phase-cycling pattern.

DESPOT1-HIFI and Modeling the Two-Component Inversion-Prepared-SPGR Signal

As described in (14), DESPOT1-HIFI involves the combined acquisition of multiangle SPGR data and at least one inversion-prepared (IR)-SPGR image. From this data, and assuming single-component T₁ relaxation, ρ, T₁ and the transmitted flip angle may be calculated. Here, the transmitted flip angle (α_T) is related to the actual “prescribed” value as α_T = κα_P, where κ is the spatially varying B₁ “calibration” field and includes tissue dielectric and RF pulse profile effects.

Two potential limitations of the DESPOT1-HIFI method may, however, prevent its use alongside mcDESPOT. The first is the assumed linear relationship between the prescribed and transmitted flip angle values (α_T = κα_P). In the presence of the frequency variations associated with slab-select gradient applied concomitant with the RF pulse, this linear relationship may be valid for a large flip angle range, or at all locations throughout the excited volume. The second is the single-component T₁ assumption in DESPOT1-HIFI. It is unclear if estimates of κ made under this assumption will be correct in multicomponent systems. In answer to this second point, if the exchange rate between the water pools is fast with respect to T₁ (i.e., the system is in the fast exchange regime), as is generally assumed, the longitudinal relaxation should be faithfully characterized by a single relaxation rate. In these systems, DESPOT1-HIFI should yield correct estimates of κ despite its simplified model.

Numerical simulations were performed to assess each of these limitations. To examine the linear flip angle relationship, a Bloch simulation was performed to evaluate the magnetization across a 14 cm imaging volume (using the same method as described in Ref. 16). T₁ and T₂ were set to 800 and 65 msec, respectively, for the simulation. The 800 µsec linear phase SLR RF pulse and gradient waveforms used for the in vivo imaging (discussed below) were used for the simulation. Simulation were performed separately for prescribed flip angle values of {2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, and 70}°. For each simulation, the transmitted flip angle and κ were determined at 256 points across the slab select direction. The percentage change between κ at 2° and 70°, for each of the 256 points, was then determined. Under ideal conditions of perfect linearity, the percent change should be zero.

To assess the accuracy of DESPOT1-HIFI κ estimates derived from two-component data, numerical simulations using theoretical data generated with two-component variants of the SPGR and IR-SPGR signal models were performed. Like the SPGR and bSSFP signals above, the two-component IR-SPGR signal may be modeled using an iterative discrete-time solution to the Bloch-McConnell expressions. The model comprises repeated “blocks” representing the inversion pulse; inversion time (TI) delay with T₁ and T₂ relaxation; and a train of N_RF readout pulses separated by TR (Fig. 2b). The transverse components are set to zero before each RF pulse to model ideal spoiling. These blocks are iteratively evaluated until the longitudinal magnetization following the inversion pulse, TI delay, and first readout pulse reach equilibrium (generally five iterations), and this value is taken as the IR-SPGR signal.

For each simulation, a numerical phantom with T_1,M = 465 msec, T_1,F = 970 msec, T_2,M = 15 msec, T_2,F = 86 msec, τ_M = 150 msec, and variable f_M was used for data generation, and consisted of 65,000 realizations with complex gaussian-distributed (0 mean, variable SD) noise added to the theoretical values. Accuracy was defined as the absolute percent difference between the mean of the 65,000 estimates and the “true” κ value, and precision was evaluated as the estimate mean divided by the SD.

DESPOT1-HIFI data were generated using the following parameters: SPGR: TR = 5.4 msec, α_P = {2, 4, 6, 8, 10, 12, 14, 16, and 18}°; IR-SPGR: TR = 5.4 msec, TI = 450 msec, α_P = 5°, N_RF = 64, as f_M was varied from 0.04 to 0.2 in increments of 0.02; as κ was varied from 0.8 to 1.2 in increments of 0.02; and for τ_M values of 150 and 350 msec. For each f_M, κ and τ_M combination, 65,000 estimates of κ were evaluated.

Building on this first simulation, the accuracy and precision of subsequently derived multicomponent mcDESPOT parameters (T_1,F, T_2,F, T_1,M T_2,M, τ_M and f_M, with δω = 0) was investigated. Two-component SPGR, IR-SPGR, and bSSFP data were generated as f_M was varied from 0.04 to 0.2 in increments of 0.02; and for κ values of 0.8, 1.00, and 1.20. Modeled “acquisition” parameters were: SPGR: TR = 5.4 msec, α_P = {2, 4, 6, 8, 10, 12, 14, 16, and 18}°; IR-SPGR: TR = 5.4 msec, TI = 450 msec, α_P = 5°, N_RF = 64; bSSFP: TR = 3.5 msec, α_P = {8, 16, 24, 32, 40, 48, 56, 64, 72}°. For the analysis κ was first estimated from the IR-SPGR and SPGR data and then used in the multicomponent model fitting. As above, 65,000 estimates of κ, T_1,F, T_2,F, T_1,M, T_2,M, τ_M and f_M were evaluated at each f_M and κ combination. For comparison, the simulation was also performed with κ set to one (i.e., no assumed B₁ inhomogeneity). The nonlinear mcDESPOT fitting was performed using the stochastic region contraction approach described in Appendix B and illustrated in Fig. 3. As effects of B₀ inhomogeneity were not included in the simulation, only a six-parameter fit for T_1,F, T_2,F, T_1,M, T_2,M, τ_M, and f_M was performed.

Accounting for Off-Resonance Affects of RF Phase-Cycling on bSSFP

The general two-component expression for the bSSFP signal is presented in Eq. (A5) with off-resonance effects denoted by (σ_RF + δω). To account for B₀ inhomogeneity and off-resonance we assume, as a first approximation, that within each imaging voxel the free water and myelin associated water share similar off-resonance characteristics, i.e., δω = δω_M = δω_F. This additional parameter may then be included in the nonlinear fitting (Appendix B) to model these effects.

For some values of δω, however, the resultant bSSFP signal will be very near zero for all values of α. To ensure sufficient data for robust fitting, bSSFP may be acquired with two or more phase-cycling patterns. Phase-cycling refers to incrementing the phase increment applied to each pulse along the RF pulse train [denoted by σ_RF in Eq. A3)]. Adjusting σ_RF shifts the spatial location of the signal null (13), thus ensuring at least one set of bSSFP data will have nonzero signal.

The accuracy and robustness of this phase-cycling approach was first examined through numerical simulations. Using the previously detailed numerical phantom (with f_M = 0.14) two-component SPGR and bSSFP data were generated while δω was varied from 0 to 450 Hz in increments of 25 Hz. Parameters used for the SPGR data were: TR = 5.4 msec and α = {2, 4, 6, 8, 10, 12, 14, 16, and 18}°; and for the bSSFP data: TR = 3.5 msec, α = {8, 16, 24, 32, 40, 48, 56, 64, and 72}° and phase-cycling increments of ω_rF = 0° and π°. Sixty-five thousand T_1,M, T_1,F, T_2,M, T_2,F, f_M, τ_M and δω estimates were calculated for each value of δω. For comparison, 65,000 estimates of T_1,M, T_1,F, T_2,M, T_2,F, f_M, and |_M were calculated using a six-parameter fit with δω = 0 using only the σ_RF = π° data. The accuracy and precision of each derived parameter were calculated at each δω increment and compared between the corrected and noncorrected methods. For this and the following simulation, no B₁ error was assumed.

Accuracy and precision of the approach were next evaluated over a range of myelin water volume fraction and off-resonance conditions. Theoretical data were generated for f_M from 0.02 to 0.20 (in increments of 0.04) and 800 from 0 to 450 Hz (in increments of 50 Hz). Sixty-five thousand T_1,M, T_1,f, T_2,M, T_2,F, f_M, τ_M and δω estimates were calculated for each f_M and δω combination and the accuracy and precision examined. Estimates were also calculated assuming δω = 0 and using only the σ_RF = π° data for comparison.

Putting It All Together—In Vivo Imaging

Building on the numerical simulations, a series of in vivo human brain imaging experiments were performed to evaluate the combined DESPOT1-HIFI + phase-cycled bSSFP mcDESPOT approach. All imaging experiments were performed on a GE Signa Excite 1.5 T HDx scanner with an eight-channel receive-only head array RF coil. Informed consent was obtained from each healthy volunteer before scanning and the study was performed with ethics approval from the host institute.

SPGR, IR-SPGR, and phase-cycled bSSFP data were acquired of four healthy volunteers, 28–42 years of age. A common 22 cm × 22 cm × 14 cm sagittal field-of-view (FOV), 128 × 128 × 80 acquisition matrix (with 0.75 partial Fourier acquisition), and the following sequence-specific acquisition parameters were used:

SPGR: TE/TR = 2.3 msec/5.4 msec, α = {2, 4, 6, 8, 10, 12, 14, 16, and 18}°, and receiver BW = ±19.23 kHz; IR-SPGR: TE/TR/TI/Tr = 2.3 msec/5.4 msec/350 msec/566 msec, α = 5°, BW = ±19.23 kHz; bSSFP: TE/TR = 1.8 msec/3.5 msec, α = {8, 16, 24, 32, 40, 48, 56, 64, and 72}°, σ_RF = 0° and π°, and receiver BW = ±50 kHz.

To achieve the various low-angle readout flip angles, the B₁ amplitude was scaled with duration held constant (800 µsec). The same RF pulse profile was used for all acquisitions, save the adiabatic IR-SPGR inversion pulse.

To exemplify the RF pulse profile effects, data were also acquired in the coronal orientation using the same acquisition parameters. Acquisition time per volunteer was ~16 min for each orientation (32 min in total). Including participant preparation, localizer scans, and participant removal, the entire protocol was ~40 min in length.

Alongside the 1.5 T data acquisitions, a proof-of-concept acquisition was also performed at 3 T. A 22 cm × 22 cm × 16 cm sagittal FOV, 128 × 128 × 92 acquisition matrix (with 0.6 partial Fourier acquisition) and sequence-specific acquisition parameters of:

SPGR: TE/TR = 2.1 msec/4.8 msec, α = {2, 4, 6, 8, 10, 12, 14, 16, and 18}°, and receiver BW = ±31.25 kHz; IR-SPGR: TE/TR/TI/Tr = 2.1 msec/4.3 msec/450 msec/656 msec, α = 5°, BW = ±31.25 kHz; bSSFP: TE/TR = 1.8 msec/3.6 msec, α = {6, 12, 18, 24, 30, 36, 42, 48, and 56}°, σ_RF = 0° and π°, and receiver BW = ±62.5 kHz. Total protocol (preparation, localizer scans, and removal) time for this acquisition was less than 20 min.

To minimize pulse profile effects, Dowell and Tofts (17) have suggested simply turning off the slab select gradient during the RF pulse. A limitation of this approach, however, is that it is restricted to applications in which the field of view covers the entire object of interest to avoid aliasing. Thus, imaging of the spinal cord or non-sagittal whole-brain acquisitions is problematic without the use of surface RF receiver coils. To investigate the utility of this approach, 3 T data were acquired of a healthy volunteer as described above with and without the slab-seclect gradient and rewinders turned on.

A potential pitfall of complex models is “over-fitting” the data, i.e., extracting information from models not supported by the data. To avoid this, one can use any number of parsimony metrics or statistical tests to justify the use of more complex models. Using the acquired 3 T data, we compared the conventional single-component DESPOT1 and DESPOT2 fits with the more complex two component mcDESPOT fit. Residuals to the single- and two-component fits were calculated and then compared using the Bayesian information criterion (BIC), (18) to determine which model was better supported by the data. The BIC compares the magnitude of the residuals while penalizing for the number of estimated parameters.

Following data acquisition, each volunteer’s data was linearly coregistered (19) and nonbrain signal was removed from each image (20). Three sets of mcDESPOT multicomponent parameter maps were calculated for each of the four volunteers: (1) No B₀ or B₁ correction, using only the SPGR and bSSFP (σ_RF = π°) data; (2) B₀ correction, but no B₁ correction; and (3) both B₀ and B₁ correction. After calculation of the maps, the coronal and sagitally oriented data were coregistered and mean frontal white matter T_1,M, T_1,F, T_2M, T_2,F, f_M, and τ_M were compared. Calculation time for the fully corrected whole brain maps was ~9 h on an eight-core (2.8 GHz) MacPro with in-house C code.

Simulation Results

To confirm the linear relationship between the prescribed and transmitted flip angles (α_T = κα_P) across the excited volume over the range of flip angles used, a Bloch simulation was performed examining this relationship for flip angles between 2° and 70°. Results of this simulation are shown in Fig. 4. Figure 4a – c contains results of the simulated pulse profiles for 2°, 10°, and 20° flip angles. Figure 4d contains the calculated κ values for each flip angle across the excited slab. Figures 4e and 4f are plots of the percent difference between κ calculated at 2° and 70° and 2° and 50°, respectively. Results show that while ideal linearity is not achieved across the full flip angle range, the deviation from linearity is small and occurs at flip angles greater than 50°. Between 2° and 70°, the maximum percent difference in κ was under 4%, and between 2° and 50°, the maximum percent difference was less than 0.6% toward the periphery of the volume.

Results of our investigation of the accuracy of DESPOT1-HIFI κ estimates derived from two-relaxation component data assuming a single-component T₁ model are shown in Fig. 5. For both residence times investigated (150 and 350 msec) the method provides accurate κ estimates, with less than 1% error, over the full range of κ and f_M investigated despite assuming a simplified T₁ relaxation model.

The accuracy of subsequent multicomponent estimates calculated with and without the DESPOT1-HIFI κ values is shown in Fig. 6. Incorrect knowledge of the transmitted flip angle results in substantial deviations in each of T_1,M, T_1,F, T_2,M, T_2,F, f_M, and τ_M. This result appears to differ from a prior investigation, which showed little influence of B₁ influence in mcDESPOT [9]. Further examination revealed a nonlinear response of f_M in the presence of flip angle errors. For 0.9 < κ <1.1, little effect is noted in the f_M estimates. For κ outside this range, however, the inaccuracy of the f_M estimates increases rapidly. In comparison, calibration of the transmitted flip angles using DESPOT1-HIFI results in substantial improvement in mcDESPOT parameter accuracy across the full range of myelin volume fraction investigated.

Results of the RF phase-cycling simulations are shown in Figs. 7 and 8, which show corrected and uncorrected T_1,M, T_1,F, T_2,M, T_2,F, f_M, and t_M estimates over a range of myelin water volume fraction and off-resonance values. The sensitivity of mcDESPOT to off-resonance is visually demonstrated in Figs. 7 and 8. The uncorrected estimates show significant deviation from the true theoretical values at even subtle δω. The inclusion of off-resonance modeling significantly improves the accuracy of the derived parameter estimates, with less than 5% absolute error observed in each of the parameters across the investigated myelin fraction and off-resonance values.

In Vivo Imaging Results

Representative imaging results are shown in Figures 9– 12. Figure 9 contains comparative sagittal and coronally acquired myelin water volume fraction’s maps, reformatted as axial slices at approximately the same level, calculated without B₀ and B₁ correction; with only B₀ correction; and with both B₀ and B₁ correction. These data clearly demonstrate the influence of flip angle and off-resonance errors in the derived maps. Errors related to flip angle inhomogeneity are greatest along the slab-select direction at the edges of excited volume and correspond to the profile of the RF pulse. Errors resulting from off-resonance conditions are maximal near the sinuses, inner ear, and brain stems.

Figure 10 displays maps of each of the six mcDESPOT parameters (T_1,M, T_1,F, T_2,M, T_2,F, f_M, and τ_M) calculated without B₀ or B₁ correction; with only B₀ correction; and with both B₀ and B₁ correction. As observed in Fig. 8, these data exemplify the B₀ and B₁-related errors in each of the multicomponent parameters and demonstrate the corrective ability of the presented flip angle and off-resonance modeling techniques. In general, it is noted that flip angle errors appear to predominately affect the T₁ measures, whereas off-resonance principally influence the T₂ measures. However, both have significant influence on the derived myelin water fraction estimates.

In Fig. 11, we show axially reformatted images through the uncorrected; B₀ only; and fully B₀ and B₁ corrected myelin water volume fractions maps at different levels throughout the brain. Again, these images demonstrate not only the severity of artifact in the uncorrected maps but also their effective removal and elimination when B₀ and B₁ are modeled and appropriately accounted for.

As noted throughout the numerical simulations, areas of off-resonance translate into over-estimations of myelin water fraction. Similarly, incorrect knowledge of the transmitted flip angle, such as at the edges of the excited volume, also translate into myelin water fraction over-estimation. A comparison of mean myelin water fraction values obtained from frontal white matter is shown in Table 1 and supports the qualitative results observed in Figs. 8–11. Averaging the corrected results from both orientation acquisitions across all four volunteers, mean frontal region white matter values were: T_1,M = 387 ± 26 msec, T_1,F = 946 ± 53 msec, T_2,M = 10.4 ± 3.2 msec, T_2,F = 108 ± 12 msec, f_M = 0.232 ± 0.0017, and τ_M = 98 ± 26 msec. No significant difference was noted between the values derived from the sagittal or coronally acquired data. These values compare favorably with previously reported values at 1.5 T of T_1,M = 350–460 msec, T_1,F = 970 msec, T_2,M = 10–55 msec, T_2,F = > 55 msec, and f_M = 0.12–0.32, respectively (12,13).

As proof-of-concept, demonstration of the technique at the higher 3 T field strength is shown in Fig. 12. As with the results shown previously at 1.5 T, the inclusion of B₁ and B₀ modeling removes these effects from the calculated mcDESPOT parameter maps.

As has been suggested previously, one approach to eliminate pulse profile affects is to eliminate the slab-selection gradient. Results of similar acquisitions with and without the slab-selection gradient are shown in Fig. 13 and show similar results regardless of the presence or absence of the gradient. This lends further credence to our simulation findings, which showed only a slight deviation in the linear relationship between the prescribed and transmitted flip angles even in the presence of off-resonance effects associated with the slab-selection gradient.

Finally, results of the single- and two-component residuals and BIC analysis are shown in Fig. 14. Despite the complexity of the two component mcDESPOT model, this analysis shows the mcDESPOT model is justified by the data throughout the majority of the image.