Investigation of the influence of B0 drift on the performance of the PLANET method and an algorithm for drift correction.

PURPOSE: The PLANET method was designed to simultaneously reconstruct maps of T1 and T2 , the off-resonance, the RF phase, and the banding free signal magnitude. The method requires a stationary B0 field over the course of a phase-cycled balanced SSFP acquisition. In this work we investigated the influence of B0 drift on the performance of the PLANET method for single-component and two-component signal models, and we propose a strategy for drift correction.
METHODS: The complex phase-cycled balanced SSFP signal was modeled with and without frequency drift. The behavior of the signal influenced by drift was mathematically interpreted as a sum of drift-dependent displacement of the data points along an ellipse and drift-dependent rotation around the origin. The influence of drift on parameter estimates was investigated experimentally on a phantom and on the brain of healthy volunteers and was verified by numerical simulations. A drift correction algorithm was proposed and tested on a phantom and in vivo.
RESULTS: Drift can be assumed to be linear over the typical duration of a PLANET acquisition. In a phantom (a single-component signal model), drift induced errors of 4% and 8% in the estimated T1 and T2 values. In the brain, where multiple components are present, drift only had a minor effect. For both single-component and two-component signal models, drift-induced errors were successfully corrected by applying the proposed drift correction algorithm.
CONCLUSION: We have demonstrated theoretically and experimentally the sensitivity of the PLANET method to B0 drift and have proposed a drift correction method.

INTRODUCTION

Quantitative MRI is used widely to obtain quantitative characteristics of tissues related to their biological and physiological properties, based on which tissues can be differentiated and associated with specific diseases. The measurement of the relaxation times of tissues (or quantitative relaxometry) is particularly important for clinical applications in oncology and regenerative medicine.1 Many different techniques exist for quantitative relaxometry, such as standard inversion recovery and multiecho spin‐echo‐based approaches,2, 3, 4 many rapid SSFP approaches like inversion‐recovery TrueFISP,5, 6 the variable flip angle approach or DESPOT‐1 and DESPOT‐2,7, 8, 9 the triple‐echo steady‐state approach,10 the MR fingerprinting approach,11 and many others.

We recently introduced a method called PLANET to simultaneously reconstruct maps of the relaxation times T1 and T2, the local off‐resonance ∆f0, the RF phase, and the banding free signal magnitude, using phase‐cycled balanced SSFP (bSSFP) data.12 The method is based on linear least‐squares fitting of an ellipse to phase‐cycled bSSFP data in the complex signal plane and subsequent analytical parameter estimation from the fitting results.

A bSSFP signal is strongly dependent on local resonant frequency, and the use of RF phase cycling shifts the off‐resonance profile of the signal dependent on the RF phase increment. The main requirement of the PLANET model is a stationary main magnetic field (B0) over the course of the acquisition, which usually consists of 8‐10 dynamics and takes around 10 minutes for full brain coverage with FOV of 220 × 220 × 100 mm3 and voxel size of 1 × 1 × 4 mm3 (without any acceleration technique). In this case, accurate and precise parameter estimation can be achieved for a single‐component voxel, as we showed in a previous study,13 whereas systematic errors in parameter estimates are expected when multiple signal components with different relaxation times and frequencies are present within a voxel.13

Due to intensive gradient activity, the requirement of a stationary B0 field can be difficult to meet, and as a result, B0 drift can occur, which might result in errors in the estimated parameters. The severity of drift effect depends on the field strength, history of gradient activity and heating of metallic components of the scanner, the acquisition time, the used gradient mode, shimming, and more, which vary among different systems and over time.

The purpose of this work was to investigate the effects of B0 drift and to assess the influence of drift on the quantitative parameters estimated using the PLANET method. We first derived a geometrical interpretation of the influence of drift on a single‐component phase‐cycled bSSFP signal based on a mathematical model. Subsequently, based on this geometrical interpretation, we developed a strategy for drift correction. Next, we experimentally showed the influence of drift on the parameter estimates for a single‐component model in a phantom and for a two‐component model of white matter (WM) in the human brain. We assessed the effects of drift for both single‐component and two‐component signal models and evaluated the performance of the drift correction algorithm in both cases by looking at drift‐induced errors in the quantitative parameter estimates. Finally, we performed numerical simulations for both single‐component and two‐component signal models to verify the experimental results.

METHODS

How drift influences a single‐component phase‐cycled bSSFP signal

For a single‐component model with monoexponential transverse and longitudinal relaxation, the complex phase‐cycled bSSFP signal can be represented as an ellipse in the complex plane14, 15 as follows:where M eff, a, and b are parametric functions of T1, T2, TR, and flip angle (FA); ∆θ is the user‐controlled RF phase increment (rad); is the rotation angle of the ellipse around the origin with respect to its vertical form14; ; Δf 0 is the local off‐resonance (Hz); φ is the combined RF transmit and receive phase; and δ is the chemical shift of the species (Hz) with respect to the water peak.

After substitution, Equation 1 can be rewritten as

A graphical representation of the ellipse described by this equation in the complex plane is shown in Figure 1A.

Figure 1

Schematic representation of the influence of B0 drift on the data points in the complex signal plane. A, The ellipse without drift. B, Time‐dependent displacement of the data points along the ellipse due to B0 drift. The nondrifted data points (not influenced by drift) are blue; the displaced data points are red. C, Time‐dependent rotation of the data points around the origin. The rotated data points are green. D, The ellipse without drift is blue, and the ellipse fitted to the drifted data points (influenced by drift) is green. E, The vertical conic forms of the ellipse without drift is blue, and the ellipse fitted to the drifted data points is green

The frequency drift is modeled as , where Δf drift(t) is the time‐dependent frequency drift during PLANET acquisition and is equal to γΔB 0(t), where γ is the gyromagnetic ratio equal to 42.58 MHz/T. Then Equation 2 becomes

The first part of Equation 3 multiplied by the first exponential represents the elliptical equation in Equation 2, but with a modified time‐dependent RF phase‐increment scheme: . This corresponds to a drift‐dependent displacement of all data points along the ellipse as illustrated in Figure 1B. Note that if only this effect of drift is taken into account, the data points remain on the “nondrifted” ellipse (i.e., the ellipse fitted to the data points not influenced by drift).

There is, however, another effect of the drift, caused by the last exponential factor in Equation 3. Using this factor, drift leads to an additional rotation of the data points around the origin, as illustrated in Figure 1C,D. As drift is time‐dependent, the rotation angle differs per data point.

These two effects together relocate the data points in the complex plane away from the nondrifted ellipse and result in a nonelliptical distribution of the points. Ignoring the effects of B0 drift, fitting an ellipse to the drifted data points (i.e., data points relocated by drift from their nondrifted positions) would lead to different fit results compared with the fit for the nondrifted case (Figure 1E). After performing PLANET postprocessing,12 this would result in errors in the parameter estimates.

We propose a drift correction method that aims to relocate each data point back to the position in the complex plane that it would have in the case without drift.

Drift correction method

Based on this analysis, we propose a 3‐step drift correction algorithm:

Calculation of the spatio‐temporal B0 drift during the phase‐cycled PLANET acquisition Δf drift, (i, j)(t), where n is the number of the dynamic acquisition, t is the time, and (i, j) are the spatial indices of the voxel. One phase‐cycled PLANET acquisition consists of N acquisitions. Assuming temporally linear drift over the duration of the phase‐cycled PLANET acquisition, the frequency drift over the nth phase‐cycled acquisition is estimated by

where the total drift over the phase‐cycled PLANET acquisition Δf total drift(i, j) is calculated by subtracting 2 reference B0 maps acquired right before and right after PLANET acquisition, and t = n∆t is the time point within the PLANET acquisition scheme corresponding to nth dynamic acquisition, where the dynamic acquisitions in the phase‐cycled acquisition each have a duration ∆t.

Correction of M eff, T1, and T2 by multiplying the experimental complex data by , the geometrical equivalent of which is the rotation of each drifted data point around the origin back to the nondrifted ellipse.

Correction of Δf 0 and φ by defining , which geometrically moves the drifted data points along the ellipse back to their nondrifted positions, where Δθ is the user‐controlled RF phase increment , n = {1, 2…N}, covering a full cycle of 2π.

Temporal drift model

As we observed experimentally, B0 drift on a long‐time scale can be represented by an exponential function. In the proposed drift correction algorithm, we assumed the temporal evolution of the drift to be linear over the duration of one PLANET acquisition.

Here we compared two temporal drift models: where Δf drift, (i, j)(t) is the frequency drift over time t, A drift and b drift are parameters describing the global spatial drift characteristics, and (i, j) are the spatial indices of the voxel.

A linear model described by Equation 4; and

An exponential model described by

Accuracy and precision in the estimated parameters and drift correction performance

The accuracy of the method was assessed by calculating relative errors (ε) in T1, T2, Δf 0, and φ RF estimates before and after drift correction, as follows:

The precision of the method was assessed by calculating the relative SD of T1, T2, Δf 0, and φ RF estimates before and after drift correction, as follows:where refers to the average of the values X affected by drift; refers to the average of the values estimated after drift correction, assuming a true value of X for parameters T1, T2, ∆f 0, and φ RF; i is an index for the voxels in a region of interest (ROI) (in experiments) or the current number of the simulation (in numerical simulations); and Z is the total number of voxels in an ROI (in experiments) or the total number of simulations (in numerical simulations).

To quantify the drift correction on T1, T2, Δf 0, and φ RF estimates, Δcor were determined as

Experiments

Phantom experiments

To investigate the effects of drift on a single‐component phase‐cycled bSSFP signal model, and to test the drift correction algorithm, MRI experiments on a phantom (1.5‐L plastic bottle filled with an aqueous solution of MnCl2·4H2O [concentration of approximately 55‐60 mg/L]) were performed on a clinical 1.5T MR scanner (Ingenia; Philips, Best, Netherlands). A 15‐channel head coil was used as a receiver. The experimental design is shown in Figure 2A.

Figure 2

Experimental drift measurements in the phantom. A, Experimental design. B, Reference B0 maps (obtained in 3D and shown only for 1 axial slice of the phantom). C, Calculated B0 drift maps. D, Total drift map. E, Drift over 65‐minute interval for 1 voxel in the center of the slice: Exponential temporal drift curve (red) fitted to the experimental data points (blue dots); lines connecting the experimental data points (blue). F, The example of original (measured) data points and drift‐corrected data points with corresponding elliptical fits for 1 voxel in the middle of the slice. G, Example of original vertical ellipse and drift‐corrected vertical ellipse for 1 voxel in the center of the slice

To compare two temporal drift models, after the first reference B0 mapping acquisition we repeated 5 times the PLANET acquisition over the course of a 65‐minute experiment. Each PLANET acquisition was followed by the reference B0 mapping acquisition. This allowed us to assess the performance of the linear drift correction method in the presence of more pronounced long‐term drift, which is expected to be nonlinear. Six reference B0 maps were acquired over the course of a 65‐minute experiment, alternated with phase‐cycled PLANET acquisitions. The reference B0 maps were obtained using a dual‐echo approach. The B0 drift over the duration of each PLANET acquisition was calculated by subtracting the two reference B0 maps acquired just before and after the PLANET acquisition concerned.

For the PLANET acquisition with more severe drift, the T1, T2, Δf 0, and φ RF maps were reconstructed. Both linear and exponential temporal drift models were used to correct the drift over this PLANET acquisition; T1, T 1 , T 2 , Δf 0 , and φ RF maps were recalculated by applying the drift correction. Because T1 estimates depend on FA (see Equation 10 in Shcherbakova12), a B1 mapping sequence was acquired, and voxel‐wise B1 correction was performed while calculating the T1 maps. The B1 maps were calculated using a dual‐TR actual FA imaging technique.16 The reference T1 and T2 values of the phantom were measured using a simultaneous spin‐echo and inversion‐recovery method (2D MIXED).3 Relevant protocol parameter settings are presented in Table 1.

Table 1

Protocol parameter settings

Phantom experiment at 1.5 T
PLANET: 3D phase‐cycled bSSFP
FOV (m3)		Voxel size (mm3)		Acq. matrix			Rec. matrix				TR (ms)			TE (ms)				Flip angle				Number of RF increment steps				NSA			Readout direction					Dummy pulses					Total scan time (minutes)
160 × 160 × 159		1.1 × 1.1 × 3		144 × 145 × 53			160 × 160 × 53				10			5				30				10				1			AP					6 seconds for each dynamic					13:46
Reference B1 map (3D dual‐TR SPGR)
FOV (m3)		Voxel size (mm3)		Acq. matrix			Rec. matrix				TR (ms)				TE (ms)				Flip angle					NSA				Readout direction					Parallel imaging			Total scan time (minutes)
160 × 160 × 159		2.5 × 4 × 3		64 × 40 × 53			160 × 160 × 53				[30; 150]				1.82				60					1				AP					SENSE 2 in RL			03:12
Reference off‐resonance map (3D dual‐echo SPGR)
FOV (m3)		Voxel size (mm3)		Acq. matrix			Rec. matrix				TR (ms)			TE (ms)					Flip angle					NSA				Readout direction					Parallel imaging			Total scan time (minutes)
160 × 160 × 159		2.5 × 4 × 3		64 × 40 × 53			160 × 160 × 53				30			[4.6; 9.2]					60					1				AP					SENSE 2 in RL			01:04
Reference T1 and T2 map (2D MIXED)
FOV (m3)	Voxel size (mm3)		Acq. matrix		Rec. matrix			TR SE (ms)				TR IR (ms)					IR delay (ms)				TE (ms)								NSA		Readout direction				Total scan time (minutes)
160 × 160 × 3	2 × 2 × 3		80 × 80 × 1		160 × 160 × 1			1500				2000					500				[30; 60; 90; 120; 150; 180]								1		AP				04:47
In vivo experiments at 1.5 T and 3 T
PLANET: 3D phase‐cycled bSSFP
FOV (m3)		Voxel size (mm3)		Acq. matrix		Rec. matrix				TR (ms)			TE (ms)			Flip angle				Number of RF increment steps					NSA			Readout direction				Dummy pulses							Total scan time (minutes)
220 × 220 × 100		0.98 × 0.98 × 4		220 × 220 × 25		224 × 224 × 25				10			5			20				10					1			AP				10 seconds for each dynamic							11:00
Reference B1 map (3D dual‐TR SPGR)
FOV (m3)		Voxel size (mm3)		Acq. matrix			Rec. matrix				TR (ms)				TE (ms)				Flip angle				NSA				Readout direction					Parallel imaging						Total scan time (minutes)
220 × 220 × 100		3.44 × 4 × 4		64 × 55 × 25			224 × 224 × 25				[30; 150]				1.82				60				1				AP					SENSE 1.5 in RL						02:43
Reference off‐resonance map (3D dual‐echo SPGR)
FOV (m3)		Voxel size (mm3)		Acq. matrix			Rec. matrix				TR (ms)			TE (ms)					Flip angle					NSA						Readout direction			Parallel imaging			Total scan time (minutes)
220 × 220 × 100		3.44 × 4 × 4		64 × 40 × 53			224 × 224 × 25				30			[4.6; 9.2]					60					1						AP			No			01:21
Reference T1 and T2 map (2D MIXED)
FOV (m3)	Voxel size (mm3)		Acq. matrix		Rec. matrix				TR SE (ms)			TR IR (ms)					IR delay (ms)				TE (ms)										NSA			Readout direction			Total scan time (minutes)
220 × 220 × 4	2 × 2 × 4		80 × 80 × 1		112 × 110 × 1				2500			5000					500				[30; 60; 90; 120; 150; 180; 210; 240]										1			AP			14:00

Abbreviations: AP, anterior–posterior; bSSFP, balanced SSFP; IR, inversion recovery; NSA, number of signals averaged; RL, right–left; SE, spin echo; and SPGR, spoiled gradient echo.

The accuracy and precision in the parameter estimates, before and after linear drift correction, were assessed using Equations 6 and 7. Deviations quantifying the drift correction performed on T1, T2, Δf 0 , and φ RF estimates were calculated using Equation 8. The ROI analysis was performed on the quantitative T1 and T2 maps calculated for the phantom: The ROI (approximately 2000 voxels) was placed in the center of the phantom on the selected slice.

In vivo experiments

To investigate the effects of drift for a tissue in which multiple components are present, and to test the drift correction algorithm, experiments on the brain of healthy volunteers were performed on clinical 1.5T and 3T MR scanners. Both protocols included B1 mapping acquisition, one PLANET acquisition in between two reference B0 mapping acquisitions, and the reference T1 and T2 mapping acquisition with the protocol parameter settings given in Table 1.

A 2.5‐ms‐long RF excitation pulse was used in each PLANET acquisition to minimize magnetization transfer effects.17 Image registration (rigid) and Gibbs ringing filtering18 were applied to the brain data before performing the PLANET reconstruction. The B0 drift maps were filtered (using a circular averaging filter with radius of 15) before applying the drift correction algorithm.

The T1, T2, Δf 0 , and φ RF maps were calculated before and after drift correction. The B1 correction was performed voxel‐wise while calculating the T1 maps. Deviations quantifying the drift correction performed on T1, T2 , Δf 0 , and φ RF estimates were calculated using Equation 8. The ROI analyses were performed on the quantitative T1 and T2 maps for both 1.5T and 3T data. The ROIs were manually delineated in WM on the selected slice in the area where the drift was the most pronounced (each ROI was approximately 100‐150 voxels). The precision of the T1 and T2 measurements was evaluated by calculating SDs on T1 and T2 maps over the ROIs.

Numerical simulations

Drift‐induced errors and drift correction for a single‐component signal model

To investigate the errors caused by B0 drift for a single‐component tissue model, numerical simulations were performed with relaxation times equal to those of the phantom material: T1 = 430 ms and T2 = 50 ms. The following parameter settings were used in the simulations: FA in the range of 0º‐45º, TR in the range of 0‐20 ms, 10 RF phase‐increment values , n = {0, 1,.., 9}, M eff = 10 000, single peak with δ = 0, Δf 0 = 5 Hz, and φ RF = −0.2 rad (these values were obtained experimentally in the phantom). The B0 drift was assumed to be linear over time and spatially independent (Δf drift = [1 2 3 4 5 6 7 8 9 10] Hz), as we found in the experimental results in the phantom. Gaussian noise was added independently to the real and imaginary data, resulting in an SNR of about 230, which corresponds to the experimentally measured SNR in the phantom. The number of performed Monte Carlo simulations was 10 000.

The accuracy and precision in the T1, T2, Δf 0 , and φ RF estimates were assessed using Equations 6 and 7.

Drift‐induced errors and drift correction for a two‐component signal model

To investigate the errors in the parameter estimates caused by B0 drift in the case in which two components are present in the signal, numerical simulations were performed for WM tissue at 3 T, which is known to be a two‐component tissue.19, 20 Two single peaks were used in the simulations: the on‐resonant dominant component and the smaller component with an average frequency shift of Δf = 20 Hz.20 The dominant component has T1D = 1000 ms and T2D = 80 ms, with a volume fraction of 0.88; the smaller component has T1S = 400 ms and T2S = 10 ms, with a myelin water fraction of 0.12. The off‐resonance Δf 0 = 10 Hz was used, and the RF phase offset φ RF = −0.15 rad was used. Gaussian noise was added independently to the real and imaginary data, resulting in an SNR ranging from 30 to 150. The number of performed Monte Carlo simulations was 10 000.

The simulations were performed using the complex phase‐cycled bSSFP signal described by Equations 7 and 8 in our previous study13 for 3 cases:

No B0 drift;

Linearly increasing over time and spatially independent frequency drift Δf drift = [1 2 3 4 5 6 7 8 9 10] Hz; and

Linearly increasing over time and spatially independent frequency drift Δf drift = [1 2 3 4 5 6 7 8 9 10] Hz with subsequently applied proposed drift correction algorithm.

The accuracy and precision in the T1, T2, Δf 0 , and φ RF estimates were assessed using Equations 6 and 7, where the true parameter values were taken for the dominant WM component.

All simulations and calculations were performed in MATLAB R2015a (The MathWorks, Natick, MA).

RESULTS

Experimental results in the phantom

Experimental results in the phantom are shown in Figure 2. Six reference B0 maps acquired before and after each of 5 PLANET acquisitions and the corresponding calculated B0 drift maps are presented in Figure 2B,C. A total drift of 28 Hz over a 65‐minute scanning session was observed (Figure 2D). The temporal drift was analyzed voxel‐wise, and the example of the experimental data for one voxel (in the center of the phantom) is shown in Figure 2E. Over a 65‐minute scanning time the temporal drift can be considered as an exponential function. Over the 11‐minute duration of the PLANET acquisition, the drift with an average value of 10 Hz can be very well approximated with a linear function.

As an example, the initial data points and the data points after drift correction for one voxel are shown in Figure 2F, with corresponding elliptical fits. The conic vertical forms of these ellipses are shown in Figure 2G. The ellipses are different, as expected due to the drift.

The estimated T1, T2, Δf 0 , and φ RF maps of the phantom before and after linear drift correction, as well as the reference T1, T2, and Δf 0 maps, are shown in Figure 3A‐C. The drift correction was performed for the first PLANET acquisition, where the drift was more severe. The performance of linear and exponential drift correction was very similar; therefore, we did not include the maps of T1, T2, Δf 0 , and φ RF after exponential drift correction in Figure 3. A reference RF phase map was not acquired and therefore is not shown. Deviations quantifying the amount of linear drift correction performed on all quantitative parameters are shown in Figure 3D,E. The drift correction decreased the T1 values by about 4%, increased the T2 values by about 8%, decreased the Δf 0 values by about 120%, and increased the φ RF values by about 3%. The magnitude image with white vertical and horizontal lines used for T1 and T2 profiles, and T1 and T2 profiles on estimated, corrected, and reference maps, are shown in Figure 3F,G. The T2 estimates are more sensitive to the drift than T1 estimates.

Figure 3

Experimental results in the phantom. A, Reference T1, T2, and Δf0 maps. B, The T1, T2, Δf0, and φ RF maps before drift correction. C, The T1, T2, Δf0, and φ RF maps after linear drift correction. D, Maps of absolute Δcor quantifying the drift correction performed on T1, T2, Δf0, and φ RF. E, Maps of relative Δcor quantifying the drift correction performed on T1 and T2. F, Magnitude image with white vertical and horizontal lines in the center of the slice, used for T1 and T2 profiles. G, The T1 and T2 profiles (T1 and T2 values representing single voxels along the selected lines on estimated, corrected, and reference maps). The values were averaged voxel‐wise between the horizontal and vertical selected lines

The quantitative ROI analysis for parameters T1, T2, Δf 0, and the relative errors in these parameters before and after drift correction, are presented in Supporting Information Table S1. The T1 values were overestimated due to drift by around 5% compared with the reference values, and the corrected T1 values were in agreement with the reference values with an accuracy of 1%. The T2 values were underestimated due to drift by about 10% compared with the reference values, and after drift correction they were in agreement with the reference values with an accuracy of 2%. The Δf 0 values estimated by means of PLANET were about 80% overestimated due to drift, and after drift correction they became similar to the reference Δf 0 acquired right before PLANET acquisition.

Despite the fact that there are no directly visible B0 drift‐related artifacts in quantitative parameter maps, there are B0 drift‐related errors in the quantitative T1, T2, and Δf 0 maps.

Experimental results in the brain

1.5 T

The results of the experiment in the brain of a healthy volunteer at 1.5 T are shown in Figure 4. The results are presented for 1 central slice. Spatially homogeneous drift was observed over 11‐minute PLANET acquisition (Figure 4A) with an average value of 9 Hz. The T1, T2, Δf 0, and φ RF maps before and after linear drift correction are shown in Figure 4C,D. The banding‐free magnitude and the reference T1 and T2 maps are shown in Figure 4B. Deviations quantifying the amount of linear drift correction performed on all parameters are shown in Figure 4E. The mean T1 and T2 values were calculated for WM. The results of the ROI analysis are given in Table 2 for the estimated, drift‐corrected, reference, and literature‐published T1 and T2 values.

Figure 4

Experimental results obtained in the brain at 1.5 T. A, Reference B0 maps before (1) and after (2) PLANET acquisition, the corresponding drift map (1,2), and the drift map filtered using a circular averaging filter. B, Banding‐free magnitude image, the reference T1 and T2 maps. C, The T1, T2, Δf0, and φ RF maps before drift correction. D, The T1, T2, Δf0, and φ RF maps after linear drift correction. E, Maps of absolute Δcor quantifying the drift correction performed on Δf0, T1, T2, and φ RF. F, Maps of relative Δcor quantifying the drift correction performed on T1 and T2

Table 2

Quantitative results from the experiments in the brain at 1.5 T and 3 T: estimated, drift‐corrected, and reference T1 and T2 values in white matter

1.5 T	Estimated values		Drift‐corrected values		Reference values		Literature‐published valuesa
ROI	T1	T2	T1	T2	T1	T2	T1	T2
1	508 ± 34	55 ± 3	501 ± 26	61 ± 2	596 ± 21	76 ± 2	621 ± 61 (9)	58 ± 4 (9)
2	460 ± 33	55 ± 3	458 ± 20	57 ± 3	596 ± 19	77 ± 3	561 ± 12 (21)	73 ± 2 (21)
3	475 ± 30	56 ± 4	475 ± 32	58 ± 4	595 ± 27	84 ± 4
4	495 ± 33	57 ± 5	487 ± 38	58 ± 5	629 ± 24	85 ± 5
Mean	485 ± 33	56 ± 4	480 ± 30	59 ± 4	604 ± 23	81 ± 4
3T	Estimated values		Drift‐corrected values		Reference values		Literature‐published values
ROI	T1	T2	T1	T2	T1	T2	T1	T2
1	678 ± 33	51 ± 2	664 ± 35	53 ± 2	771 ± 19	68 ± 2	832 ± 1 (22)	80 ± 1 (22)
2	660 ± 42	50 ± 3	642 ± 39	53 ± 3	781 ± 18	70 ± 3	1084 ± 45 (23)	69 ± 3 (23)
3	636 ± 30	50 ± 2	624 ± 29	53 ± 2	771 ± 16	69 ± 2	781 ± 61 (24)	65 ± 6 (24)
Mean	658 ± 36	50 ± 2	643 ± 35	53 ± 2	774 ± 18	69 ± 2

The mean T1 and T2 values at 1.5 T were calculated for 1 slice of the brain by averaging over 4 regions of interest (ROIs) (each around 150 voxels) in white matter on corresponding T1 and T2 maps. The mean T1 and T2 values at 3 T were calculated for 1 slice of the brain by averaging over 3 ROIs (each around 100 voxels) in white matter on corresponding T1 and T2 maps.

Numbers in parentheses are reference citations.

After drift correction, T1 values decreased by about 1% compared with the uncorrected values, and T2 values increased by about 5% compared with the uncorrected values (Figure 4E,F and Table 2). The B0 values decreased by about 50%, and the corrected B0 map resembles the reference B0 map acquired right before the PLANET acquisition. The RF phase maps almost did not change after drift correction.

3 T

A spatially inhomogeneous drift was observed over the same 11‐minute PLANET acquisition in the brain of another healthy volunteer at 3 T, with a maximum value of 10 Hz for the selected slice (Figure 5A). The T1, T2, Δf 0, and φ RF maps before and after linear drift correction, the banding‐free magnitude, and the reference T1 and T2 maps, and deviations quantifying the amount of linear drift correction performed on all parameters are shown. The results of the ROI analysis are given in Table 2. Similar to the results at 1.5 T, the T1 values after drift correction did not change much; they locally decreased by about 2% compared with the uncorrected values in the area with more pronounced drift. The T2 values were more sensitive to drift, and after drift correction they increased by about 5%‐6% compared with the uncorrected values in the area with more pronounced drift (Figure 5E,F and Table 2). The B0 values decreased by about 50%, and the corrected B0 map resembles the reference B0 map acquired just before the PLANET acquisition. The RF phase maps did not change much after drift correction.

Figure 5

Experimental results obtained in the brain at 3 T. A, Reference B0 maps before (1) and after (2) PLANET acquisition, the corresponding drift map (1,2), and drift map filtered using a circular averaging filter. B, Banding‐free magnitude image, the reference T1 and T2 maps. C, The T1, T2, Δf0, and φ RF maps before drift correction. D, The T1, T2, Δf0, and φ RF maps after linear drift correction. E, Maps of absolute Δcor quantifying the drift correction performed on Δf0, T1, T2, and φ RF. F, Maps of relative Δcor quantifying the drift correction performed on T1 and T2

A remaining underestimation of about 20% in T1 values and about 30% in T2 values compared with the reference values are found in Figures 4 and 5 and Table 2.

Simulation results

Single‐component phase‐cycled bSSFP signal model of the phantom

Relative errors and SDs in T1, T2, Δf 0, and φ RF estimates for a single‐component signal model of the phantom are presented in Figure 6, affected by linear drift (Figure 6A,B) and after applying drift correction (Figure 6C,D). As shown, drift induced errors depend on the choice of FA and TR. For the combination of FA = 30° and TR = 10 ms, which was used in the experimental setup, the quantitative analysis of the errors is presented in Table 3. The T1 values are overestimated due to drift by about 4% compared with the true values; T2 values are underestimated due to drift by about 8% compared with the true values; Δf 0 values are overestimated due to drift by about 100%; and φ RF values are underestimated due to drift by about 5%. After applying the proposed drift‐correction algorithm, relative errors in all estimated parameters are almost zero, which demonstrates an accurate performance of drift correction. The SDs in all estimated parameters are not affected by drift correction much and are below 5%.

Figure 6

Simulation results for a single‐component signal model. A,B, Relative errors (ε) and relative SDs in T1, T2, Δf 0, and φ RF estimates (in percent) compared with their true values in the presence of linear over time and spatially independent drift: Δf drift = [1 2 3 4 5 6 7 8 9 10] Hz. C,D, Relative errors (ε) and relative SDs in T1, T2, Δf 0, and φ RF estimates (in percent) compared with their true values after applying linear drift‐correction algorithm. The initial settings: T1 = 430 ms, T2 = 50 ms, Δf 0 = 5 Hz, φ RF = −0.2 rad, and 10 RF phase increments

Table 3

Quantitative results of simulations for a single‐component signal model of the phantom and a two‐component signal model of WM at 3 T: the accuracy and precision in T1, T2, Δf 0, and φ RF estimates without drift, with drift, and after drift correction

Parameter	Single‐component model of the phantom, TR = 10 ms, FA = 30°		Two‐component signal model of WM at 3 T, TR = 10 ms, FA = 20°
	Relative error ε, %	SD, %	Relative error ε, %	SD, %
No drift
T1	0.07	1.7	−30.2	2.1
T2	0.04	1.0	−34.8	1.4
Δf 0	−0.01	2.8	13.7	1.5
φ RF	−0.01	4.7	20.0	5.6
Drift
T1	4.2	1.9	−29.5	2.0
T2	−7.5	1.1	−38.8	1.5
Δf 0	97.8	1.2	58.3	1.4
φ RF	−4.7	4.8	24.7	6.9
Drift‐corrected
T1	0.06	1.8	−30.6	1.9
T2	0.03	1.1	−34.3	1.5
Δf 0	0.01	2.7	10.1	1.9
φ RF	0.01	4.5	23.2	6.7

FA, flip angle; and WM, white matter.

These results are in agreement with the experimental results for the phantom shown previously: The simulated expected errors due to drift match the calculated errors in the estimated parameters.

Two‐component phase‐cycled bSSFP signal model

The simulation results for a two‐component signal model of WM are shown in Figure 7. Relative errors in T1, T2, Δf 0, and φ RF are shown for 3 cases: no drift, linear drift, and after applying the drift‐correction algorithm. The errors in the estimated parameters depend on the choice of FA and TR. As we showed in a previous study,13 in WM brain tissue the PLANET postprocessing results in systematic errors in estimated T1, T2, and Δf 0 values due to the presence of a second myelin‐related component in WM. Here we can observe similar behavior for the case without drift.

Figure 7

Simulation results for a two‐component signal model of WM at 3 T (T1D = 1000 ms, T2D = 80 ms, T1S = 400 ms, T2S = 20 ms, ∆f = 20 Hz, and myelin water fraction = 0.12). Value of Δf 0 = 10 Hz, φ RF = −0.15 rad, and 10 RF phase increments. Relative errors (ε) caused by drift in T1, T2, Δf 0, and φ RF estimates (in percent) compared with their true values (of the dominant component). A, Without drift. B, In the presence of linear over time and spatially independent drift: Δf drift = [1 2 3 4 5 6 7 8 9 10] Hz. C, After linear drift‐correction algorithm

For the combination of FA = 20º and TR = 10 ms, which was used in the experimental setup, the quantitative analysis of the errors is presented in Table 3. The T1 values are underestimated by 30% without drift, underestimated by 29.5% in the presence of drift, and underestimated by 30.5% after drift correction. The T2 values are underestimated by 35% without drift, underestimated by 39% in the presence of drift, and underestimated by 34.5% after drift correction. The Δf 0 values are overestimated by 14% without drift, overestimated by 58% in the presence of drift, and overestimated by 10% after drift correction. The φ RF values are overestimated by 20% without drift, overestimated by 25% in the presence of drift, and overestimated by 23% after drift correction.

The drift correction performed on all estimated parameters predicted by the simulations for the combination of TR = 10 ms and FA = 20º is similar to the drift correction performed experimentally in the brain: After drift correction, T1 values decreased by about 1% compared with the uncorrected values; T2 values increased by about 5% compared with the uncorrected values; Δf 0 values decreased by about 48%; and φ RF values decreased slightly by about 1.5%. In all cases, drift‐induced errors were corrected.

DISCUSSION

The PLANET method requires a stationary main magnetic field over the course of the acquisition. This requirement, however, can be difficult to meet, and as a consequence, B0 drift can occur. In this work, we investigated the sensitivity of the PLANET method to B0 drift and assessed the errors that drift can cause in the estimated T1, T2, Δf 0, and φ RF parameters.

We presented a mathematical interpretation of the influence of drift on the elliptical phase‐cycled bSSFP single‐component signal behavior and proposed a general strategy for drift correction. We demonstrated how drift influences the performance of the PLANET method experimentally in a phantom and in the brain of healthy volunteers. Consequently, we verified the effects of the drift by performing numerical simulations for the same phantom and in vivo setups.

The experimental results in the phantom showed that the drift of about 10 Hz, which occurred over the 11‐minute duration of the PLANET acquisition, induced the errors in estimated quantitative parameters: The T1 values were overestimated due to drift by about 5%; the T2 values were underestimated due to drift by about 10%; and the Δf 0 values were overestimated due to drift by about 80% compared with the corresponding reference values. The variance in the estimated parameters only slightly changed after drift correction. We demonstrated that both linear and exponential correction algorithms performed identically. The linear model for temporal evolution of the drift on a short time scale (0‐15 minutes) may be a fair approximation of the exponential drift in the experiments reported in this paper. Drift‐induced errors in T1, T2, Δf 0, and φ RF estimates in a phantom were successfully corrected by applying the drift‐correction algorithm. These results obtained experimentally were verified by numerical simulations for a similar setup: The bias and variance in all parameter estimates predicted by simulations matched the ones calculated using the experimental data of the phantom.

The investigation of the drift effects in the human brain showed that similar drift of about 10 Hz over the 11‐minute duration of the PLANET acquisition had a significant effect only on the estimated Δf 0 values: An overestimation of about 50% in Δf 0 values was caused by drift. The other quantitative parameters were only affected slightly: The drift induced an overestimation of about 1% in T1 estimates, an underestimation of about 5% in T2 estimates, and an overestimation of about 5% in φ RF estimates. The errors in the quantitative parameters calculated in the brain were in agreement with errors predicted by simulations for a similar experimental setup. The proposed drift‐correction algorithm performed well and corrected the errors caused by drift. However, the remaining underestimation by about 20%‐30% in T1 and T2 values compared with the reference and literature published values, which can be found in Figures 4 and 5 and Table 2, is not caused by B0 drift. It is caused by the effect that in WM tissue where multiple components are present, a single‐component PLANET model is not valid, as we already pointed out in a previous study.13 Obviously, such underlying errors were not and cannot be corrected by the drift‐correction algorithm. Keep in mind that any other techniques that assume a single‐component relaxation model will fail in this case as well.

The severity of drift effect depends on the field strength, history of gradient activity and heating of metallic components of the scanner, PLANET acquisition time, the used gradient mode, shimming, and more, which vary among different systems and over time. Even though the errors in estimated quantitative parameters caused by drift in human brain are small (1%‐5%) compared with the errors caused by the presence of multiple components (about 30% underestimation), as we have shown in this study, they cannot be predicted and can potentially affect reproducibility of the results, as drift effects are generally not reproducible. We have now shown that the drift‐induced errors can be successfully corrected by applying the proposed drift‐correction algorithm. Acquiring two quick low‐resolution reference B0 maps before and after the PLANET acquisition is generally a simple direct way to correct for drift and improve the quantitative parameter estimation using the PLANET method.

CONCLUSIONS

We have demonstrated that the PLANET method is sensitive to B0 drift. Although there may be no directly visible B0 drift‐related artifacts on the estimated parameter maps, drift can induce errors in these parameters. In the phantom, which can be described with a single‐component signal model, drift induced significant errors in the estimated parameters. However, in the human brain, where multiple components are present, drift had only a minor effect. We have now shown that the drift‐induced errors can be successfully corrected by applying the proposed drift‐correction algorithm for both cases of a single‐component and two‐component signal model.

TABLE S1 Quantitative results from the phantom experiment (the reference, estimated, and drift‐corrected T1, T2, and Δf 0 values) and the relative errors in estimated and drift‐corrected T1, T2, and Δf 0 values

Click here for additional data file.