PLANET: An ellipse fitting approach for simultaneous T1 and T2 mapping using phase-cycled balanced steady-state free precession.

PURPOSE: To demonstrate the feasibility of a novel, ellipse fitting approach, named PLANET, for simultaneous estimation of relaxation times T1 and T2 from a single 3D phase-cycled balanced steady-state free precession (bSSFP) sequence.
METHODS: A method is presented in which the elliptical signal model is used to describe the phase-cycled bSSFP steady-state signal. The fitting of the model to the acquired data is reformulated into a linear convex problem, which is solved directly by a linear least squares method, specific to ellipses. Subsequently, the relaxation times T1 and T2 , the banding free magnitude, and the off-resonance are calculated from the fitting results.
RESULTS: Maps of T1 and T2 , as well as an off-resonance and a banding free magnitude can be simultaneously, quickly, and robustly estimated from a single 3D phase-cycled bSSFP sequence. The feasibility of the method was demonstrated in a phantom and in the brain of healthy volunteers on a clinical MR scanner. The results were in good agreement for the phantom, but a systematic underestimation of T1 was observed in the brain.
CONCLUSION: The presented method allows for accurate mapping of relaxation times and off-resonance, and for the reconstruction of banding free magnitude images at realistic signal-to-noise ratios. Magn Reson Med 79:711-722, 2018.

INTRODUCTION

Quantitative MR imaging plays an important role in accurate tissue characterization for improving clinical diagnostic imaging and for planning, guidance and evaluation of image‐guided therapy. The mapping of longitudinal (T1) and transverse (T2) relaxation times is a particularly important tool for many clinical applications in oncology and regenerative medicine 1. Knowledge of T1 and T2 values allows optimizing the contrast‐to‐noise‐ratio between tissues by finding the optimal sequence parameter settings.

Various techniques are widely used for T1 and T2 relaxation time mapping. 2D inversion recovery spin echo (IR‐SE) and multi‐echo spin echo (ME‐SE) 2 are considered gold standard techniques, allowing accurate measurements of relaxation times. Unfortunately, scans based on these methods typically have a long acquisition time, which makes it challenging to use them in clinical practice. To speed up IR‐based T1 mapping, the Look‐Locker method was introduced 3. This approach is closely related to the IR‐SE, but instead of acquiring a single image for each inversion time, the Lock‐Locker method uses an inversion pulse followed by a train of low flip angle (FA) pulses, each followed by a read‐out, within each repetition time (TR). Although that considerably reduces the required scan time, it still is a time‐consuming 2D method, which results in a very long acquisition time to cover a complete 3D volume.

Another widely used method, which uses the variable FA (VFA) approach 4, 5, is DESPOT1 6. The method allows for rapid 3D high‐resolution T1 mapping and is easily implemented on clinical scanners. For this method, at least two acquisitions of spoiled gradient‐echo (SPGR) scans with different FAs are required. A signal model for the steady‐state is subsequently fitted to VFA data. A similar approach was developed by Deoni et al for T2 mapping under the name DESPOT2 and extended to the combined T1 and T2 mapping 7. DESPOT2 also requires at least two acquisitions of 3D balanced steady‐state free precession (bSSFP) with different FAs using prior knowledge of T1, often estimated using DESPOT1. Both methods have been shown error prone, which demands the optimization of the parameter settings. The combination of used FAs needs to be optimized for improved accuracy and precision 8, 9, 10. Furthermore, the influence of radiofrequency field inhomogeneity 11, off‐resonance effects 12, 13, and radiofrequency and gradient spoiling efficiency 14, 15 on the accuracy and precision of T1 and T2 measurements was investigated.

With the advent of stronger and faster gradient systems, bSSFP sequences have become widely used for rapid imaging with high signal‐to‐noise ratio (SNR) efficiency. Although the signal is a complex function of relaxation parameters T1 and T2, several bSSFP‐based approaches for relaxometry have been proposed, such as 2D inversion recovery TrueFISP 16, 3D triple echo steady‐state (TESS) 17. For instance, the TESS method allows simultaneous rapid 3D estimation of T1 and T2 within one single scan using two specific signal ratios between three echoes (SSFP‐FID and two SSFP‐echo modes) and an iterative golden section search algorithm.

Generally, bSSFP imaging has a high sensitivity to local magnetic field inhomogeneities, which results in banding artifacts. Radiofrequency (RF) phase‐cycled bSSFP imaging was introduced as a solution, and several algorithms were proposed to reduce banding artifacts 18.

Björk et al 19 introduced a parameter estimation algorithm to remove banding artifacts and simultaneously estimate relaxation times T1 and T2 from phase‐cycled bSSFP. In their work, a combination of linear least squares fitting followed by a subsequent nonlinear iterative fitting was used, called the two step LORE‐GN algorithm. They successfully reconstructed banding free magnitude images of a phantom and in vivo. Based only on simulations and numerical assessment of the Cramer‐Rao Bound (CRB), they concluded that simultaneous estimation of T1 and T2 from phase‐cycled bSSFP would be difficult at common SNR, because the CRB is high.

In this work, we introduce a novel approach, named PLANET, for simultaneous T1 and T2 estimation from phase‐cycled bSSFP 20. The elliptical signal model is used to describe the phase‐cycled bSSFP steady‐state signal 21. The fitting of the model to the acquired data is reformulated into a linear convex problem, which is solved directly by a linear least squares method, specific to ellipses 22. Subsequently, the relaxation times T1 and T2 are analytically calculated from the fitting results.

Our work shows that accurate mapping of the relaxation times T1, T2, the off‐resonance caused by local field deviations, and banding free magnitude is feasible for realistic SNRs and can be performed with a regular coil setup and scan protocol parameter settings.

METHODS

Elliptical Signal Model

The elliptical signal model for bSSFP was first used by Xiang and Hoff 23 to remove banding artifacts. The complex bSSFP signal right after the RF pulse (at echo time t = 0+) can be described as: where , , is the thermal equilibrium magnetization, α is the FA, TR is the repetition time, θ is the resonance offset angle (in radians), θ = θ ‐ Δθ, where , Δf0 is the off‐resonance caused by local field deviations, is the chemical shift of the species (in Hz) with respect to the water peak, Δθ is the user controlled RF phase increment (in radians). Parameters are all θ‐independent.

Parametric Equation [1] describes an ellipse in the complex signal plane. Each point on the ellipse represents real and imaginary components of the transverse magnetization, which are acquired with a certain RF phase increment, as illustrated in Figure 1. The total number of acquisitions N and the RF phase increment Δθ are user defined parameters.

Figure 1

The elliptical signal model of the bSSFP in the complex plane as a function of the resonance offset angle θ. In this particular case, and N = 10 acquisitions with different RF phase increments are shown.

Directly after the RF pulse (i.e., at t = 0+), the long axis of the ellipse is oriented vertically in the complex plane and the values for the parameters and are within the interval (0,1)21. The cross‐point is the geometric solution (GS), which is independent of θ and can be used to calculate a banding‐free magnitude image.

At echo time t = TE > 0 after the RF pulse, the signal is still a function of the resonance offset angle 21, but then the real and imaginary components of the signal are modulated by a factor and rotated around the origin by:

Taking into account the RF contribution, eddy current effects and B0 drift, the complex signal can be then described as: where is the effective magnetization, K is the magnitude of the combined receive field, is the RF phase offset, related to the combination of RF transmit and receive phases, is the extra phase errors due to eddy current effects, and is the extra phase errors due to B0 drift.

Reconstruction Method for Parametric Mapping

Essential to the PLANET method is a three‐step reconstruction algorithm to simultaneously estimate relaxation parameters T1 and T2, and an effective banding free magnitude from phase‐cycled bSSFP data. However, an additional step is required when the reconstruction of the off‐resonance map is also desired.

Step 1. Direct linear least squares ellipse fitting to phase‐cycled bSSFP data.

The first step consists of performing voxel‐wise direct linear least squares fitting of a general quadratic polynomial function to the data points in the complex plane 22: where x and y are real and imaginary components of transverse magnetization.

By minimizing the sum of squared algebraic distances of the ellipse to the data points under a proper scaling and an appropriate constraint specific to ellipses (discriminant – 4C1C3 = −1), we avoid the trivial solution C = 0 and exclude all nonelliptical fits, such as hyperbola and parabola. As a result, we find a unique set of coefficients C = [C1, C2, C3, C4, C5, C6] representing the ellipse. The fitting is based on a numerically stable version of the ellipse fit 24. This is a fast, direct, linear, non‐iterative ellipse fit. Because there are six unknowns, we need at least six data points x, which can be acquired by scanning with at least six different RF phase increment settings.

Step 2. Rotation of the ellipse to initial vertical conic form.

The rotation of data points of the ellipse to the initial vertical form, i.e., the orientation directly after RF excitation pulse, was performed by applying basic algebraic transformations to the polynomial representation of the ellipse found in the previous step (Eq. [6]). We found the rotation angle to be:

Because is defined within , we unwrapped it to cover the range by verifying that the ellipse of every voxel is vertical and that its center lies on the positive real axis.

After rotation, illustrated in Figure 2a, the conic equation for the vertical orientation can be used to describe the ellipse: where ( , 0) is the geometrical center of the ellipse, and are the semi‐axes of the ellipse.

Figure 2

a: Schematic representation of the ellipse at t = 0 + (red dashed line) and t = TE > 0 (blue solid line), which is rotated around the origin by φ. (xc, 0) = the geometrical center of the ellipse; A, B = semi‐axes of the ellipse. b: Geometrical determination of the angle using the locations of the data points (xn, yn) on the vertical ellipse.

Step 3. Analytical solution for parameters .

Parameters of the parametric form of the ellipse in Equation [4] are related to the geometric characteristics from the Cartesian form of ellipse in Equation [8] through a system of nonlinear equations (23):

We have solved this system for parameters analytically. The results are presented in the Appendix. The T1 and T2 estimates can be found from parameters a and b using the following equations:

Additional Step 4. Estimation of the local off‐resonance .

The rotation angle in Equation [7] includes the local off‐resonance and RF phase offset , which cannot be separated from Equation [5]. For simplicity, the chemical shift is ignored and the additional phase errors due to eddy current effects and B0 drift are assumed to be negligible:

The off‐resonance , however, can be estimated from the locations of the data points with different RF phase increment settings on the vertical ellipse: the precession angle during each TR depends only on the RF phase increment and the local off‐resonance , and not on the RF phase offset :

Using a Cartesian parametric equation of the ellipse: and after substitution of x and y from Equation [13] by the real and imaginary components of the signal in Equation [4], the relationship between parameters t and can be found as:

For each individual nth data point with n = {0,1,..N−1} Eq. [14] is equivalent to

The set of can be found from the set of using Equation [13]. The set of can be found geometrically from the data points on the vertical ellipse as illustrated in Figure 2b and consequently, the set of can be found from Equation [15] and can be represented by the sum of a sine function and a cosine function: .

Next, the coefficients and can be found by taking a linear least squares fitting approach and can be estimated from:

The off‐resonance can be found from Equations [12] and [16]. Because is defined within , we unwrapped it to cover the range , which results in a bandwidth ( ).

Sensitivity to FA Errors

To investigate how sensitive the method is to errors in the actual FA, simulations were performed for a range of nominal FAs between 1 ° and 90 ° and a range of deviation in actual FAs of ‐10% and +10%. The initial parameter settings were: KM0 = 1, T1 = 675 ms, T2 = 75 ms,  = 10 Hz, TR = 10 ms, TE = 5 ms,  = 0,  = 0, N = 10 phase cycles with phase increments , n = {0,1,..9}. The chosen T1 and T2 represent white matter at 1.5T. No Gaussian noise was added.

Experimental Validation

To investigate the performance of our method, both phantom and human volunteer experiments were performed on a clinical 1.5T MR scanner (Philips Ingenia, Best, The Netherlands). For all scans, a 16‐channel head coil (dS HeadSpine, Philips Ingenia, Best, The Netherlands) was used as a receive coil.

The phantom experiments were performed on a calibrated phantom consisting of gel tubes with known T1 and T2 values (TO5, Eurospin II test system, Scotland). Twelve tubes were chosen with T1, T2 combinations in the following ranges: T1 (220–1600 ms), T2 (50–360 ms).

First, the known temperature dependence of the relaxation times of the calibrated gels was used to assess the T1 and T2 values of the test tubes for the actual scanner room temperature. The temperature inside the phantom water was measured before and right after the experiment using a T‐type thermocouple. The difference between measured temperature values was below 0.5 ° and the average value was chosen for the correction.

The 3D phase‐cycled bSSFP sequence was performed with the protocol parameter settings, shown in Table 1. Complex‐valued data were acquired.

Table 1

Protocol Parameter Settings

Phantom experiments
3D phase‐cycled bSSFP
FOV (mm3)	Voxel size (mm3)	Acq. matrix	TR (ms)	TE (ms)	FA (˚)	Number of RF increment steps	NSA	Readout direction	Total scan time (min:s)
200x200x80	1.5x1.5x2.5	132x132x32	10	5	30	10	1	AP	06:51

NSA = Number of signal averages.

Reference T1 and T2 maps of the phantom were acquired using standard T1 and T2 mapping techniques. For the reference T1 mapping, a 2D turbo IR‐SE approach was used. For the reference T2 map, a 2D ME‐SE approach was used. The corresponding protocol parameter settings are shown in Table 1. A reference off‐resonance map was calculated using a dual echo SPGR method with the protocol parameter settings shown in Table 1.

Before voxel‐wise parameter estimation, all images were masked to exclude the borders of the tubes and the background from the analysis.

The reference T1 values were calculated voxel‐wise by performing the nonlinear fit of to multi TI IR‐SE data, with ρ and T1 as the fitting parameters. The reference T2 values were calculated voxel‐wise by performing the nonlinear fit of to ME‐SE data, with ρ and T2 as the fitting parameters.

To demonstrate the method in vivo, experiments were performed on the brain of three healthy volunteers on the same scanner. The protocol parameter settings for 3D phase‐cycled bSSFP are presented in Table 1. As shown in the Appendix, the FA should fulfill the condition . Thus, FA 30 ° was used, which should allow an accurate estimation from T1 > 100 ms onward for TR = 10 ms.

The reference T1 and T2 values of the brain were measured using a simultaneous (interleaved) spin echo and inversion recovery method (2D MIXED) 25 with the protocol parameter settings, shown in Table 1.

B1 map was calculated using a dual TR actual FA imaging 26 method with the protocol parameter settings, presented in Table 1. B1 correction was performed voxel‐wise for the calculated T1 maps.

To investigate an influence of magnetization transfer (MT) effects on the quantitative T1 and T2 mapping in vivo, experiments were performed on one volunteer using 3D phase‐cycled bSSFP with different RF excitation pulse durations, as suggested in work by Bieri and Scheffler 27. The protocol parameter settings were the same as shown in Table 1 (3D phase‐cycled bSSFP in vivo). The default pulse had a duration 0.84 ms. The long pulse optimized to minimize the MT effects had a duration 2.86 ms.

The signal ratio ΔS and MT ratio (MTR) were calculated as: , , where is the banding free magnitude measured with the default RF pulse, is the banding free magnitude measured with the long RF pulse (minimized MT effects).

A linear phase‐encoding profile order was used to minimize the eddy currents induced by changing phase‐encoding gradients 28.

The standard (fast channel combination) method, available on the scanner, was used for the combined phase reconstruction. Note that the RF phase offset remains the same for all dynamics with different RF phase increments settings.

To check the amount of B0 field drift, one additional dataset with RF phase increment was usually acquired at the end of the acquisition. In case of absence of B0 drift during the acquisition, the complex signals should be the same for data with and , otherwise the phase difference between these datasets is proportional to the amount of the drift.

The SNR for both phantom and in vivo data was calculated as defined in the work by Björk et al 19: where is the magnitude of nth phase‐cycled image, σ is the standard deviation of noise, is number of scans with different RF phase increment settings. The standard deviation of noise was calculated over the region of interest (ROI) on noise images (real and imaginary components), acquired dynamically using the same bSSFP sequence, without RF excitation and with no gradients applied.

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

RESULTS

Sensitivity to the Actual FA Errors

As can be observed from Figure 3, T1 estimates are highly sensitive to errors of the actual FA. The dependence is almost linear and the errors in T1 estimates increase with increasing FA. For example, 5% error (0.95) in actual FA results in 10% underestimation in T1 for the nominal FA of 30˚ and 12% underestimation in T1 for the nominal FA of 60 °. T2 estimates are not affected by the errors in actual FA.

Figure 3

Simulation results of sensitivity to the actual FA errors. The initial T1 = 675 ms, T2 = 75 ms, TR = 10 ms, TE = 5 ms. The horizontal line corresponds to the nominal FA = 10 °, which leads to collapsing of the ellipse to a line, as discussed in the Appendix.

Phantom Results

T1, T2, and Δf0 maps were first validated in a phantom magnitude images corresponding to different RF phase increment settings are presented in Figure 4a. The banding artifacts, the locations of which depend on the resonance offset angle, are shifted depending on the RF phase increment setting. The GS, representing the banding‐free effective magnitude image, was calculated for a set of acquired data and is presented in Figure 4b. The off‐resonance maps of the phantom were calculated using the PLANET method and using the reference method. The results are presented in Figures 4c,d. The two off‐resonance maps look almost similar. A minor deviation between the two calculated maps of [−2;+1] Hz was observed.

Figure 4

a: Magnitude images corresponding to different RF phase increments setting Δθ. b: The banding‐free effective magnitude. c,d: The off‐resonance maps calculated using the PLANET method and using the reference method.

T1 and T2 maps of the phantom were calculated using the PLANET method and using the reference methods. The results are presented in Figure 5. The processing time for the reconstruction of T1, T2, Δf0 and banding‐free effective magnitude for one slice was 6 s. We generally see a good quantitative agreement between reference and calculated from the PLANET method maps. However, there are some inhomogeneous regions inside some of the phantom tubes. The comparisons between the average T1 and T2 values for each of the phantom tubes are shown in Figure 6. Standard deviations in T1 and T2 were calculated for each tube. The estimated accuracy of tabulated T1 and T2 values of the test object, provided by the manufacturer, is ± 3%.

Figure 5

Experimental results from the phantom study: T1 and T2 maps calculated using the PLANET method and using the reference method.

Figure 6

Experimental results from the phantom study: comparison between average T1 (a) and T2 (b) values for the phantom tubes: blue = the tabulated values, green = calculated from the reference methods, red = calculated from the PLANET method. The mean T1 and T2 values of the gels were calculated for one slice in the center of the phantom by averaging over an ROI (around 250 voxels) inside each tube on estimated T1 and T2 maps. Precision of T1 and T2 measurement was evaluated by calculating standard deviations on estimated T1 and T2 maps over the same ROIs.

Results In Vivo

Figures 7a,b show the reference T1 and T2 maps of one axial slice of the brain of a healthy volunteer. The results of measurements in three different axial slices through the brain are shown in Figures 7c–f. The banding free effective magnitude is presented, as well as the T1 and T2 maps calculated using the PLANET method. The off‐resonance maps were calculated using the PLANET method and using the reference method. The minor observed deviation between the off‐resonance maps was [‐3;+3] Hz. The processing time for the reconstruction of T1, T2, Δf0, and for one slice was 7 s. On the T1 and T2 maps a good contrast between gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF) can be seen. The calculated mean T1 and T2 for WM and GM after B1 correction are presented in Table 2 in comparison with the reference values and those published in literature 7, 29.

Figure 7

Experimental results from the volunteer study: T1 map for one axial slice of the brain calculated using the reference method (a); T2 map for one axial slice of the brain calculated using the reference method (b); the banding free effective magnitude images (c); T1 maps calculated using the PLANET method (d); T2 maps calculated using the PLANET method (e); the off‐resonance maps calculated using the PLANET method (f); the off‐resonance maps calculated using the reference method (g). The position of the axial slice (a,b) is different from the positions of slices (c–g).

Table 2

Results from In Vivo Experiment: T1 and T2 Values Determined Using the PLANET Method and Using the Reference Method Compared With Published Valuesa

PLANET					Reference 2D MIXED method
	White matter		Gray matter			White matter		Gray matter
ROI #	T1 (ms)	T2 (ms)	T1 (ms)	T2 (ms)	ROI #	T1 (ms)	T2 (ms)	T1 (ms)	T2 (ms)
1	461 ± 19	62 ± 2	754 ± 47	80 ± 2	1	636 ± 15	75 ± 2	1016 ± 53	84 ± 3
2	466 ± 25	61 ± 2	749 ± 75	84 ± 8	2	602 ± 14	75 ± 2	1014 ± 30	82 ± 4
3	453 ± 15	62 ± 2	836 ± 70	83 ± 5	3	597 ± 13	73 ± 2	999 ± 52	84 ± 2
4	524 ± 19	64 ± 2	876 ± 71	82 ± 4	Mean	612 ± 14	74 ± 2	1010 ± 46	83 ± 3
5	512 ± 30	64 ± 2	837 ± 67	80 ± 4	Literature published valuesa
6	525 ± 18	63 ± 2	787 ± 24	98 ± 9
7	528 ± 22	64 ± 3	906 ± 48	90 ± 14	Refb	T1 (ms)	T2 (ms)	T1 (ms)	T2 (ms)
8			789 ± 24	84 ± 4	Ref IR;SE (7)	615 ± 12	69 ± 2	1002 ± 56	92 ± 3
9			787 ± 24	84 ± 4	Ref (7)	621 ± 61	58 ± 4	1060 ± 133	98 ± 7
Mean	496 ± 22	63 ± 2	813 ± 54	85 ± 5	Ref (29)	561 ± 12	73 ± 2	1048 ± 61	94 ± 6

The mean T1 and T2 values of WM were calculated for five slices of the brain by averaging over seven ROIs (each around 100 voxels) in WM on estimated T1 and T2 maps. The reference mean T1 and T2 values of WM were calculated for one slice of the brain by averaging over three ROIs (each around 100 voxels) in WM on the reference T1 and T2 maps. The mean T1 and T2 values of GM were calculated for five slices of the brain by averaging over nine ROIs (each around 30 voxels) in GM on estimated T1 and T2 maps. The reference mean T1 and T2 values of GM were calculated for one slice of the brain by averaging over three ROIs (each around 30 voxels) in GM on the reference T1 and T2 maps. The precision of T1 and T2 measurement was evaluated by calculating standard deviations on estimated T1 and T2 maps over the same ROIs.

Numbers in parentheses are reference citations.

An underestimation in T1 by approximately 15–20% in WM and GM compared with the reference values was observed. T2 values in WM were underestimated by approximately 10% compared with the reference values, T2 values in GM were determined very precisely compared with the reference values.

Examples of the acquired complex signals for WM, GM, and CSF with the corresponding elliptical fits are shown in Figure 8. The orientation of three ellipses is different, which is explained by difference in off‐resonance for those voxels.

Figure 8

Examples of the acquired complex signals for white matter, gray matter, and CSF with the corresponding elliptical fits from three voxels of the brain of a healthy volunteer.

T1 and T2 maps, calculated from datasets, acquired using the default and long RF excitation pulses, as well as signal ratio and MTR are presented in Supporting Figure S1, which is available online. The quantitative results from three ROIs placed in WM are shown in Supporting Table S1. The average relative signal loss due to MT effects in WM was found to be 13% and the average T1 shortening was 8%. The estimated SNR maps of the phantom and the brain are presented in Supporting Figure S2.

We did not observe any significant B0 drift between acquisitions with increments and , which were performed at the start and at the end of sequence. The maximum phase difference between these datasets was 0.04 rad for the phantom and 0.06 rad for the brain experiment.

DISCUSSION

The fitting of the elliptical signal model is a very important aspect of the proposed PLANET method. We reformulated the fitting procedure into a linear convex problem, which can be solved directly by using a linear least‐squares method. Prior knowledge of the elliptical trajectory in the complex plane allows to reduce the solution space to one unique solution by applying a proper scaling and an ellipse‐specific constraint. In combination with analytical solutions for parameters T1, T2, and , which take the fitting results as input, our approach becomes simple, robust, and fast. This is a clear advantage of our method compared with all iterative algorithms, which usually have longer reconstruction time and fitting problems related to local minima. The whole reconstruction time is very fast, which facilitates the adoption of the proposed method into clinical practice.

Compared with the work by Björk et al 19, who used a combination of linear fitting followed by subsequent nonlinear fitting and only four phase‐cycled acquisitions, the PLANET method requires at least six phase‐cycled acquisitions to directly fit the model to the experimental data. The inclusion of prior knowledge of the elliptical trajectory is essential and differs clearly from the methodology followed by Bjork et al. Contrary to their conclusions, which were based only on simulations, that the simultaneous T1 and T2 estimation using their algorithm is not feasible for realistic SNRs, we experimentally demonstrated that it is feasible to use the PLANET method at realistic SNRs, both in a phantom and in vivo, with a regular coil setup and protocol parameters settings.

The reported values of T1 and T2 in the phantom are in good agreement with the calibrated values and those calculated with reference methods. However, in some of tubes some inhomogeneous regions in the form of "ghosts" near the tube borders are observed in the resultant T1 and T2 maps of the phantom, which leads to an underestimation of the calculated T1 and T2 values for those tubes. This effect may have been caused by Gibbs ringing artifacts. The influence of these artifacts on estimated T1 maps needs to be further investigated and minimized.

T1 and T2 maps obtained in the brain of volunteers were generally in agreement with the reference maps and values in literature 7, 29; however, the T1 values were underestimated. B1 field inhomogeneities, resulting in errors in the actual FA, have shown a significant influence on the estimated T1 values. The errors caused by this effect depend on the used FA and were corrected using additionally acquired B1 maps. Unlike T1 estimates, T2 estimates were not affected by the errors in actual FA. MT effects were shown to have an influence on T1 quantification, particularly in WM. A partial mitigation strategy to minimize the impact of MT effects, as proposed by Bieri et al, was to use long RF excitation pulses in a combination with relatively low FA and long TR. Note that effects related to the presence of deoxyhemoglobin 30 and diffusion effects 31 were not included.

We believe that the observed underestimation of T1 and T2 in the human brain (particularly in WM) even after B1 correction might be caused by an inhomogeneous intravoxel frequency distribution and multicomponent relaxation 32, 33, 34, 35, 36, 37. The presence of different frequencies within a voxel results in asymmetries in the bSSFP signal profile, which have been found and comprehensively studied by Miller et al 38, 39. We also observed such asymmetries when we plotted the frequency responses for WM and GM. In the phantom, we did not observe asymmetries, because there are no structures with different frequency components, which can explain a good agreement of the found T1 and T2 values with the reference values.

Similar results were found in the work by Nguyen and Bieri 40. Their methodology, named MIRACLE, for T1 and T2 mapping is based on frequency‐shifted bSSFP scans with subsequent TESS processing for relaxometry. They used a similar experimental setup at 3T and showed a systematic underestimation of T1 values even after B1 correction in the brain, while the phantom results were in agreement with the reference. Particularly, they found a 40% underestimation in T1 for WM and a 20% underestimation of T1 for GM. They also believe that this is likely due to the asymmetric shape of the bSSFP frequency response in WM and GM due to presence of myelin 35, 36, 37. They investigated the effect by characterizing brain tissues with a two‐component relaxation parameter model, as proposed by Miller et al 39 and Deoni et al 41, in which the smaller myelin component had a lower combination of T1 and T2 compared with the dominant component. In their simulation, they observed a shift toward lower apparent T1 values which was in agreement with their experimental results and with the results which we presented in this study.

We believe that the performance and the results of the presented method in the brain deserves further examinations. In addition, further investigation of the in vivo protocol optimization will be the subject of our further research.

The relaxation times T1 and T2, the off‐resonance, and the banding free magnitude can be simultaneously and robustly estimated from one dynamic 3D phase‐cycled bSSFP sequence. This is an important difference and advantage compared with all existing bSSFP‐based techniques for relaxometry purposes.

Such quantitative mapping may be a useful addition to the common techniques for banding artifacts removal that rely on phase‐cycling 18. To accomplish this, just a few more additional bSSFP data sets with other RF phase increment settings are required. PLANET may be applied for investigating the local susceptibility and the electrical tissue properties: the off‐resonance maps can be used for quantitative susceptibility mapping (QSM) 42. RF phase offset maps, which can be in principle retrieved from Equation [11], could potentially be used for electric properties tomography 43.

Although bSSFP in general is a fast imaging technique with a high SNR efficiency, the disadvantage of using multiple phase‐cycled acquisitions is the increased scan duration. In this work, we used 10 steps, but theoretically, considering the number of fitting parameters, the minimum number of required steps is 6. Even though the scan duration is much shorter compared with the duration of the combined 2D gold standard IR‐SE and ME‐SE and 2D IR‐TrueFISP techniques, when 3D coverage is desired, and comparable to that of the combined 3D DESPOT1&2 or 3D TESS technique, we intend to further investigate ways to shorten the scan duration. Shorter acquisition times may be realized by minimizing the number of phase increment steps or by using acceleration techniques, such as compressed sensing 44 or dynamically phase‐cycled radial bSSFP 45.

We limited the model to the 3D acquisition mode, assuming a constant FA profile in the slice direction for each voxel. When volumetric coverage is not required, switching to the 2D acquisition mode would considerably decrease the acquisition time, but would lead to a nonideal FA profile over the slice which would compromise the required elliptical behavior of the integrated complex magnetization. An investigation of the feasibility of a 2D approach is subject of our further research.

In this study, we assumed that the chemical shift for all resonances was negligible (i.e., only water resonances present), which indeed was the case for the phantom and the brain experiments. For species with other chemical shifts, such as fat, the model should be adjusted to account for different initial conditions, corresponding to .

The method is sensitive to B0 drift, when it appears while acquiring acquisitions with different phase increment settings. This results in deviations from the single elliptical distribution in the complex plane and errors in the estimated parameters. In our experiments, we did not observe any significant B0 drift and did not compensate for it. If severe B0 drift appears, there is a need for correction.

The extension of this work will include a more detailed study of the precision and accuracy of the method in relation to the SNR.

CONCLUSIONS

We have presented a novel approach, named PLANET, for simultaneous estimation of relaxation times T1 and T2 from phase‐cycled bSSFP. Prior knowledge about the elliptical signal model was used to reformulate the fitting problem into a convex one, which can be solved directly using a linear least‐squares method. The unique ellipse‐specific solution of the fitting problem in combination with analytical solutions for T1 and T2 make our approach simple, robust and fast, additionally allowing for calculation of the off‐resonance and the banding‐free magnitude image from the same set of the acquired data.

We have demonstrated that accurate T1 and T2 mapping in a phantom as well as in the brain of healthy volunteers is feasible for realistic SNRs and can be performed with a regular coil setup and protocol parameter settings on a clinical MR scanner. We believe that the presented method may be applied in a wide range of applications.

Fig. S1. Experimental results from the MT study: (a) the banding free effective magnitude images, calculated for the default and long RF pulses; (b) T1 maps calculated for the default and long RF pulses; (c) T2 maps calculated for the default and long RF pulses; (d) the Signal Ratio and MTR calculated for the default RF pulse compared with the long RF pulse. The duration of the default RF pulse was 0.84 ms. The duration of long RF pulse was 2.86 ms.

Fig. S2. SNR maps calculated for one axial slice of the brain (a) and the phantom (b).

Fig. S3. Geometrical representation of the parameter space for a and b. The white area corresponds to the vertical ellipse for cases a > b and a < b, and the black area corresponds to the horizontal ellipse.

Table S1. The quantitative results from three ROIs placed in white matter for the default and long RF excitation pulses.

Click here for additional data file.