2.1 Description of the Pulse Sequence

The ME-MP2RAGE pulse sequence, shown in Fig 1, is obtained from MP2RAGE by acquiring multiple gradient echoes for each phase-encoding step. The signal expression for MP2RAGE derived in [5] is also valid for ME-MP2RAGE. Briefly, the signals (at any echo time, T_E,i) can be written as: (1) and (2) where M₀ is the equilibrium magnetization, E₁ ≡ exp(−T_R,GRE/T₁), E_A ≡ exp(−T_A/T₁), E_C ≡ exp(−T_C/T₁), and n is the number of k-space lines in one GRE block. The intervals T_R,GRE, T_A, and T_C, as well as the flip angles α₁ and α₂ are defined in Fig 1; η ≡ (1/2)[1 − M_z(0⁺)/M_z(0⁻)] is the inversion efficiency of the adiabatic pulse, where M_z(0⁻) and M_z(0⁺) are the longitudinal magnetizations immediately before and after the application of the pulse, respectively [35]. is the steady-state longitudinal magnetization, which is obtained by solving the Bloch equations with appropriate boundary conditions (see Eq A1.4 in [5]). Finally, all scanner-dependent parameters required for converting magnetization into signal voltage are lumped together in the scaling constant ξ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Schematic diagram of the ME-MP2RAGE sequence.

After an adiabatic inversion pulse (INV), two GRE readout blocks are collected with excitation pulse flip angles, α₁ and α₂. Both GRE blocks consist of n acquisitions of k-space lines, each of duration T_R,GRE, stepping linearly through the second phase-encoding direction (as in standard MP2RAGE). To each T_R,seq corresponds a k-space line acquisition for the first phase-encoding direction. The center of k-space is acquired at times T_I,1 and T_I,2. Examples of the additional mono-polar gradient lobes of the ME readouts are indicated by blue color. The sequence repetition time, T_R,seq, is defined as the time between two successive inversion pulses. The total acquisition time is thus defined by T_R,seq multiplied by the number of steps in the second phase-encoding direction. T_A, T_B, and T_C denote, respectively, the durations from the initial inversion pulse to the onset of the first GRE block, between the two GRE blocks, and from the end of the second GRE block to the next inversion pulse. Note that the first and the second phase-encoding direction may be interchanged.

https://doi.org/10.1371/journal.pone.0169265.g001

The (complex) signals from the two inversion-contrast image volumes are combined by computing (3) where S* denotes the complex conjugate of S and returns the real part of x. The expression for ρ can be found analytically, and it is used to invert the equation for T₁ numerically for a specific range of T₁ values to calculate the maps. Several parameter combinations can yield equivalent T₁ maps, provided that the appropriate expression for ρ (or a suitable approximation) is used, and specific requirements on such parameters are met, as detailed in [5].

The same equations can also be used to show that the (absolute) signal level of the image volume recorded at T_I,2, being acquired later on the recovery curve, is generally larger than the one recorded at T_I,1, in agreement with experimental observations. Hence, the second inversion image volume was employed for obtaining and χ estimates throughout this work. Instead, the T₁ estimates were always obtained from the inversion contrast recorded with the first echo.

2.2 Simulation Studies

In order to evaluate the sensitivity of the ME-MP2RAGE signal to variations, we used simulations based on the Bloch equations. In particular, we plotted the signal expression ρ as a function of T₁ with additional consideration of ±20% and ±40% variations in (through the flip angles). The convergence of the varied curves to the original signal expression indicates insensitivity to inhomogeneity of the field.

Because of the well-known inhomogeneity at 7 T, the nominal flip angle, α_nom, may differ significantly from the effective flip angle, α_eff, generated experimentally in a tissue. For more realistic comparisons, the flip angle values of the simulations and the experiments were matched using the efficiency factor η_α ≡ α_nom / α_eff, where α_eff was determined experimentally using maps.

The (ME-)MP2RAGE sequence yields T₁ maps through an inversion-recovery acquisition where a mono-exponential recovery constant is estimated from two support points, T_I,(1,2), by inverting the signal expression for T₁. Similarly, maps are obtained via a fitting procedure, where a mono-exponential decay constant is interpolated on ME measurements.

It is possible to investigate, through Monte-Carlo-like simulations, with the addition of noise [36], how the Signal-to-Noise Ratio (SNR) and the number and distribution of the support points affect the quality of the relaxation constant estimates. In particular, n_T = 100 different relaxation time values from the intervals T₁/ ms ∈ [500, 3500] and were assumed, and up to N = 20000 simulations of the obtained signal levels were performed for each value with random variation of the noise contribution to compute standard deviations (SD), σ_T, of the relaxation-time estimates. Based on typical experimental results, the SNR was set to 25, 50, or 100, with T_I,1/ ms ∈ [500, 1500] and T_I,2/ ms ∈ [1500, 3500] (constrained to T_I,2 − T_I,1 ≥ 1000 ms) in simulations of T₁ mapping, and with T_E,i/ ms ∈ [2, 30] and n_E ∈ [3, 5] in simulations of mapping (n_E is the number of echoes). The mean and SD of the σ_T values for the full range of relaxation times, μ_σ and σ_σ, were then used to assess the robustness of the estimation.

Unless otherwise noted, simulations, data analyses (see below), and the visualization of the results were performed using the “Scientific Python” ecosystem [37–39].

2.3 Acquisition Parameter Considerations

This work targets at acquiring multiple quantitative maps of the full brain at sub-millimeter resolution with an acquisition time similar to current clinical practice, which implies that only a restricted set of acquisition parameter values are desirable or accessible. In order to achieve an advantage from the simultaneous mapping approach, the following conditions should be met: (i) the T₁ map obtained with ME-MP2RAGE should be as accurate as a corresponding map acquired with MP2RAGE, and (ii) the acquisition time required for ME-MP2RAGE should be shorter than the sum of the acquisition times of corresponding MP2RAGE and ME-FLASH scans.

We specifically considered the following influences from the most relevant acquisition parameters: (a) lower values of α_1,2 translate into improved insensitivity but also lower SNR, which directly affects the accuracy and precision; to maximize SNR, α₂ was set to the Ernst angle, assuming T₁ ≈ 1.6 s as the mean of the expected values for WM (≈ 1.2 s) and GM (≈ 2.0 s) [10]; (b) the number of k-space lines per GRE block affects the block duration and limits the range of possible inversion times; n increases with matrix size and Field Of View (FOV) and may be decreased through Parallel Acquisition Techniques (PAT); (c) the choice of T_I,(1,2) has an effect on the sensitivity and the accessible value range for T₁ estimates as well as its precision; (d) the repetition time of the sequence, T_R,seq, defines the total acquisition time and also has an (indirect) impact on the sensitivity to variation; (e) a short first T_E is desirable to maximize the SNR for T₁ mapping; however, a relatively long (later) T_E in the order of 15 ms is useful for χ mapping (at the cost of longer T_R,GRE with concomitant constrains for matrix size and FOV) as it increases the Contrast-to-Noise Ratio (CNR) of the phase images and avoids possible anisotropy bias in the WM caused by non-linear phase evolution at short T_E [25, 26]; finally, acquiring more echoes improves mapping but might require a larger bandwidth, Δν.

The relatively large number of parameters (and, hence, optimization criteria) and their complex relations makes it impossible to define an unbiased cost function that simultaneously maximizes the accuracy and precision of all maps.

For this reason, a manual optimization approach was adopted: initially, the FOV, matrix, and PAT parameters were defined to determine n; then, we accommodated the maximum number of echoes within a fixed T_R,GRE, long enough for reliable and χ mapping; lastly, we adjusted all other parameters to achieve sufficient insensitivity, and T₁ mapping accuracy with the minimum possible T_R,seq.

For the MP2RAGE acquisitions, a similar optimization approach was used, since the values suggested by [5] were not optimal for the desired resolution and FOV. In this case, the shorter T_R,GRE (due to n_E = 1) allowed to increase n, which was exploited to decrease the acquisition time. This was achieved both by shortening T_R,seq and by switching the first and second phase-encoding directions, so that the main (T_R,seq) loop is performed along the GRAPPA-accelerated direction. However, in this case, the obtained acquisition parameters cause a systematically lower and potentially ambiguous T₁ estimate in the CSF value range, which is not relevant in most applications. Note that this feature is intrinsic to the MP2RAGE-based T₁ estimation procedure (see also: Fig 2) and the acquisition parameters should be adjusted to avoid this issue in the interval of interest.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Estimated T₁ values as a function of (a) the ME-MP2RAGE and (b) the MP2RAGE signal intensity parameter ρ.

The green lines correspond to the effective acquisition parameters while red/blue lighter/darker lines indicate, respectively, ±20% and ±40% offsets of . The flip angle values are adjusted for the accuracy factor η_α. Note that the estimated T₁ is limited to an upper value of approximately 3 s for the acquisition parameters that were used for MP2RAGE, which would result in T₁ values exceeding this limit (e.g., in CSF) to be underestimated.

https://doi.org/10.1371/journal.pone.0169265.g002

Exploratory measurements (summarized in the Supporting Information) with a large coverage of different acquisition parameters were used to infer their effects, in relation to the corresponding simulations. The manually chosen parameters were then tested on a cohort of subjects to provide additional experimental evidence (allowing for group statistics).

2.4 MRI Experiments

Quantitative 3D MRI of the brain was performed in 19 healthy subjects (9 females, 23–31 years) on a Magnetom 7T scanner (Siemens Healthcare, Erlangen, Germany).

This study was approved by the Ethics Committee at the Medical Faculty of the University of Leipzig (application n.: 273-14-25082014). All participants had given informed written consent prior to the examination.

The ME-MP2RAGE pulse sequence was provided by Siemens. The acquisitions were performed using a monopolar readout (to avoid potential effects related to a constant offset in the gradients system), Generalized auto-calibRating Partially Parallel Acquisitions (GRAPPA) along the first phase-encoding (PE) direction, and 6/8 partial Fourier factors (both for the first and the second PE direction). The in vivo experiments were divided into two studies.

In Study 1, exploratory measurements were performed with a circularly-polarized transmit/24-channel receive head coil and a 32-channel receive version of the same coil design (Nova Medical, Wilmington, MA, USA) to evaluate the following range of acquisition parameters: FOV between 224×224×145.6 mm³ and 192×168×112 mm³ with nominal isotropic resolutions in the 0.6–0.9 mm range; GRAPPA acceleration factors, f = 2–3; T_R,seq = 4700–8000 ms; T_I,(1,2) = 750–1100, 2750–3500 ms; α_1,2 = 2–5, 3–10 deg; T_E = 2.32–23.05 ms; n_E = 3–6; and Δν = 280–500 Hz/px yielding total acquisition times, T_acq = 5:55–22:25 min. MP2RAGE and ME-FLASH data were also recorded for the validation of the results with acquisition parameters similar to previously published values [4–6] or in-house standard settings. The FOV, nominal resolution, and GRAPPA acceleration were matched to the specific ME-MP2RAGE acquisition; further parameters were: for MP2RAGE T_R,seq = 5000–8000 ms; α_1,2 = 4–5, 3–10 deg; T_I,(1,2) = 800, 2400 ms; T_E = 2.15 ms; Δν = 280 Hz/px; T_acq = 9:42 min; and for ME-FLASH T_R = 31–35 ms; α = 10–11 ms; T_E = 2.94–29.59 ms; n_E = 5–6; Δν = 280–500 Hz/px; T_acq = 6:04–12:40 min.

From the exploratory results, a preferred set of ME-MP2RAGE parameters (see Table 1) was derived to test, in Study 2, the intra-subject reproducibility with a FOV of 192 × 144 × 134.4 mm³ (axial orientation, first PE direction from right to left, second PE with 14% oversampling) and an acquisition matrix of 320 × 280 × 224 (i.e., 0.6 mm isotropic nominal resolution).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Details of the manually optimized acquisition parameters that were used for the simulations in Fig 2 and for the acquisitions of Study 2 with a nominal spatial resolution of 0.6 mm (isotropic).

https://doi.org/10.1371/journal.pone.0169265.t001

These scans were recorded with the 32-channel coil on 7 of the 19 subjects. To obtain reference data, further measurements were made in the same sessions with MP2RAGE and ME-FLASH (parameters included in Table 1) using the same FOV, matrix size, orientation, GRAPPA acceleration, and partial Fourier settings.

A typical session, thus, consisted of two MP2RAGE, two ME-FLASH and two ME-MP2RAGE acquisitions acquired in random order without repositioning.

In some sessions, additional maps were acquired with 3 mm isotropic nominal resolution employing an in-house modification (based on complex instead of magnitude images) of the Actual Flip-angle Imaging (AFI) technique [40]. Due to time restrictions, these maps were acquired in place of an ME-FLASH or MP2RAGE scan; thus, in such cases, ME-FLASH or MP2RAGE reproducibility results are not available.

2.5 Image Processing

T₁ maps were reconstructed on-line by algorithms integrated in the Image Calculation Environment (ICE) provided by the vendor. maps were reconstructed off-line by log-linear least-squares fitting to the signal magnitudes, using the standard polynomial fit approach based on singular value decomposition. The QSM analysis was implemented in C++ using the ODIN library [41]. Susceptibility mapping was based on the phase of the echo acquired at the longest T_E for ME-MP2RAGE (T_E = 14.77 ms), and a similar value (T_E = 15.00 ms) was used for the ME-FLASH reference scans in order to improve comparability. The phase images were unwrapped using a Fourier method [42], which offers the advantage of describing singularities as continuous functions, and high-pass filtered using the Sophisticated Harmonic Artifact Reduction for Phase data (SHARP) approach [43]. Filtered phase images were converted into units of ppb (division by 10⁻⁹ γ T_E B₀; γ is the gyromagnetic ratio of the proton), and quantitative maps were calculated with the Superfast Dipole Inversion (SDI) method [44].

Image segmentation and registration of brain tissues were accomplished using the FMRIB Software Library (FSL), Ver. 5.0 [45]. For each subject, a brain mask was generated (from the GRE readout at T_I,2 for MP2RAGE or ME-MP2RAGE) using the Brain Extraction Tool (BET2) (with the ‘robust brain center estimation flag’ [-R] active), a Gaussian filter on BET’s result (for smoother edges) with SD σ = 0.5 px, and thresholding for values above 5% of the maximum. The extracted masks were furthermore shrunk (to avoid border effects) by the application of a binary erosion filter with a “spherical” kernel of diameter 1 px, iterated five times to remove the outer layer of CSF. Maps were then linearly co-registered using the FMRIB Linear Image Registration Tool (FLIRT) with rigid-body transformations, the correlation ratio metric and the previously obtained mask as weights. Once the volumes were registered, the non-brain tissues were masked out.

2.6 Statistical Analyses

Results obtained with ME-MP2RAGE were compared against MP2RAGE as reference for T₁ maps and against ME-FLASH as reference for and χ maps. The equivalence of the maps was evaluated voxel-by-voxel with two-dimensional (2D) correlation histograms and difference images. For the quantification of correlations and mutual consistency of the maps, the following parameters were investigated:

1. The squared Pearson’s correlation coefficient, r², as an indication of how well the test measurements reproduce the reference measurements (neglecting the noise of the reference measurements).
2. The means and SDs of the image volumes’ difference, as an estimate of respectively accuracy and precision, (4) as well as their absolute difference, (5) where x_i and y_i are defined by the map intensity of a voxel in the test data and in the reference data, respectively, and N is the number of voxels after the masking step. Note that μ_|D| and σ_|D| are equivalent (except for a factor ) to the average and SD of the Euclidean distance from the identity line in the 2D histogram, and reflect the combined effects of deviation from the identity line and (random) spreading.

While r² is often used in correlation studies, it only reflects potential deviations from a linear relationship, but it is not sensitive to the slope of the linear component of the relationship. It is thus of limited relevance if the measured values should be, ideally, identical as in the current case. As r² is dimensionless, it more directly permits inter-modality comparisons. The accuracy and the precision of the measurement can be estimated by μ_D and σ_D, respectively (in units of the quantity being measured). Additionally, if the accuracy is good, that is when μ_D ≈ 0, then μ_|D| (along with σ_|D| as its error) may be used as an estimate of the overall reproducibility. A sensible reference for these parameters is the size of the value range of the obtained maps. Note that these parameters are quite sensitive to outliers and may vary considerably depending on the efficacy of some of the procedures used in this work (e.g., the brain tissue masking, the fits used to calculate the maps, the registration, the considered value range, etc.).

To investigate the effects of potential errors arising from the registration procedure, we performed the voxel-by-voxel analysis described earlier on a given map in comparison to the same map after application of a linear rigid transformation. More specifically, translations in axial direction were examined for voxel offset values between 0.0 and 1.0 px; rotations about the axial direction and centered in the geometrical center of the FOV were examined for rotation angles between 0.0 and 1.0°; and rotations followed by translations were examined within the same value ranges. A complementary estimate of the effects of the registration step was obtained by applying the usual voxel-by-voxel analysis to “self-registered” maps. More precisely, reference maps were compared with the same maps registered back onto themselves after the application of a significant rigid transformation (in this study, a 10° rotation). These analyses were performed for T₁ maps obtained with MP2RAGE and for and χ maps obtained with ME-FLASH.

Because of the relatively complex reconstruction technique required for QSM, it is not possible to use an equivalent simulation approach as for relaxometry to study the impact of the SNR. Nevertheless, it is possible to investigate the effects of SNR on the final QSM result by using an approach similar to what was done to investigate the effects of mis-registration. Particularly, we reconstructed several QSM maps obtained from the phase of the ME-FLASH complex images with the addition of Gaussian noise (with increasing SD in a range of 0–30 arb.units, corresponding to up to ≈ 100% decrease in SNR). The effect of the SNR on the χ maps was evaluated using the same voxel-by-voxel comparison analyses described earlier.

Additionally, we performed an analysis on a Region-Of-Interest (ROI) basis, where we calculated the group average parametric values across the different brain areas: frontal lobe (Fro.), temporal lobe (Tem.), parietal lone (Par.), occipital lobe (Occ.), insula (Ins.), putamen (Put.), caudate (Caud.), thalamus (Tha.), cerebellum (Cereb.) as defined by the MNI Structural Atlas. Additionally, white matter (WM) from the Harvard-Oxford Subcortical Structural Atlas was considered. Volumes from the second inversion images of the ME-MP2RAGE data were registered non-linearly to the MNI152 1 mm nominal resolution template (using default options) and such transformations were then used to bring the maps to the same space. Then, for each area (defined by the probabilistic masks thresholded for above 80%, which may incidentally leave out some portions of the WM), all the subjects from Study 2 were group-averaged, and the group SD was also calculated.

3.1 Simulations

The results obtained by relaxometry simulations (an overview of which is given in S3 and S4 Tables with additional details reported in S1 and S2 Figs) indicate that, for a given SNR, the number of echoes and their distribution is crucial for mapping, while the specific choice of inversion times (for the range examined in our study) is less critical for T₁ mapping. As a rule of thumb, choices of inversion times according to T_I,1/2 < T₁ < 2T_I,2 and echo times according to (where max indicates the longest T_E) yield robust estimates of T₁ and , respectively, unless the SNR is insufficient. Detailed results for the specific settings of T_I and T_E that were selected for Study 2 are shown in Fig 3. They indicate that the accuracy and precision (represented graphically by the distance from the identity line and by the size of the error bars, respectively) depend on the exact value of the relaxation times, and at this level of SNR are comparable between the test (ME-MP2RAGE) and the reference (MP2RAGE for T₁, ME-FLASH for ) sequences. Additional information on the SNR dependence can be found in the support material indicated above.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Simulations for the expected accuracy and precision ofT₁ and maps at SNR = 50.

Rows refer to T₁ (a, b) and (c, d) simulations, while columns indicate the choice of the simulation parameters reflecting the acquisition: ME-MP2RAGE (a, c), MP2RAGE (b) or ME-FLASH (d). Each panel contain a plot of the estimated relaxation time as a function of the exact value, with the error bars indicating the standard deviations, σ_T, for N = 20000 simulations. The distance from the identity line indicates the expected accuracy, while the size of the error bars reflect the expected precision. For T₁ / exponential recovery simulations, the range 0.5–3.5 s was probed, and the exact acquisition parameters simulated were T_I,(1,2) = 800, 2400 ms for ME-MP2RAGE, and T_I,(1,2) = 750, 2900 ms for MP2RAGE. For T₂ / exponential decay simulations, the range 2–60 ms was probed, and the exact acquisition parameters simulated were n_E = 4, T_E,1 = 2.5 ms, ΔT_E ≈ 4.2 ms for ME-MP2RAGE, and n_E = 5, T_E,1 = 3.0 ms, ΔT_E = 6.0 ms for ME-FLASH.

https://doi.org/10.1371/journal.pone.0169265.g003

Fig 4 shows combined histograms of the flip-angle accuracy factor measured in the brain with the 32-channel coil. The mean plus/minus SD were 0.8 ± 0.2. Hence, flip angle values used in the simulations were divided by 0.8 and rounded to the nearest integer to obtain corresponding settings for α_nom in the experimental protocols.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. inhomogeneity displayed as histograms of the flip-angle accuracy factor inside the brain of six healthy human volunteers.

Voxels outside the head or not containing brain tissue were masked out. The dashed lines indicate the mean (thick line) plus/minus the standard deviations (thin lines) of η_α across subjects.

https://doi.org/10.1371/journal.pone.0169265.g004

The variation of simulated intensity parameters ρ of the (ME-)MP2RAGE signal in dependence of T₁ is presented in Fig 2. The underlying acquisition parameters are those from Table 1. These graphs were used to visually inspect the sensitivity and the acquisition parameters were selected accordingly. The level of the steady-state magnetization is the major factor affecting sensitivity, while the accessible T₁ range is mostly influenced by the specific choice of α_1,2 and T_I,(1,2). Therefore, the number of k-space lines acquired in each GRE block, n, and T_R,GRE, which relate directly to the image resolution and SNR, have a stronger impact on the sensitivity when the inter-acquisition times T_A, T_B, and T_C are shorter. Therefore, a relatively long T_R,seq is often needed for higher accuracy and robustness.

3.2 MRI Experiments

The results from the exploratory Study 1 (summarized in the S5–S7 Tables) were used to investigate the impact of the acquisition scheme on the mapping. Further results obtained with a nominal resolution of 0.6 mm will be discussed in greater detail below. For higher resolutions, the statistical analyses did not indicate a strong impact from the particular imaging sequence. At lower resolutions, T₁ maps obtained with ME-MP2RAGE and MP2RAGE were also essentially identical, whereas and χ maps obtained with ME-FLASH seemed slightly more accurate than those obtained with ME-MP2RAGE. This counter-intuitive observation is explained by the different PE schemes used in ME-MP2RAGE acquisitions with lower (e.g., 0.9 mm with first PE direction from head to foot) or higher (e.g., 0.6 mm with first PE direction from left to right) nominal resolutions, which yielded stronger restrictions in the choice of T_E,max (and, hence, less accurate fits) for the lower-resolution scans.

Based on the combined results from the simulations and Study 1, the parameters in Table 1 were defined for an in-depth evaluation of data acquired with a high nominal resolution of 0.6 mm. An example of the data obtained with the ME-MP2RAGE sequence is shown briefly in Fig 5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Example of the images obtained from an ME-MP2RAGE acquisition.

Each row represents a different echo time. The columns show in order: first inversion magnitude (1st) and phase (2nd); second inversion magnitude (3rd) and phase (4th). Magnitude images are shown in arb.units, while phase image are in radians, both using a gray scale. Note that: (i) the phase images for the first inversion point show an abrupt change in their value corresponding to the zero crossing of the signal in the T₁ recovery curve; (ii) the phase images for the second inversion point present some coil combination pole artifacts resulting in a corresponding degradation of the QSM maps at these locations.

https://doi.org/10.1371/journal.pone.0169265.g005

A comparison of T₁ maps acquired with ME-MP2RAGE and with MP2RAGE is shown in Fig 6. Similar comparisons of and χ maps derived with ME-MP2RAGE and with ME-FLASH are shown in Figs 7 and 8, respectively. These results are based on single subject acquisitions. Note that the χ maps are displayed in gray scale to be consistent with the QSM literature, but a diverging color map with linear luminance could be a better option for representing signed data (see also S8 Fig). The group statistics results are reported in Table 2. Individual subject results are reported in S8–S10 Tables.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Single subject acquisition of T₁ maps acquired with (a) ME-MP2RAGE, (c) MP2RAGE and (e) their difference, and test-retest reproducibility evaluation with voxel-by-voxel correlation 2D histograms for (b) ME-MP2RAGE, (d) MP2RAGE, and (f) ME-MP2RAGE versus MP2RAGE.

Note how the CSF voxel are underestimated by the MP2RAGE as a result of the specific choice of the acquisition parameters.

https://doi.org/10.1371/journal.pone.0169265.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Single subject acquisition of maps acquired with (a) ME-MP2RAGE, (c) ME-FLASH and (e) their difference, and test-retest reproducibility evaluation with voxel-by-voxel correlation 2D histograms for (b) ME-MP2RAGE, (d) ME-FLASH, and (f) ME-MP2RAGE versus ME-FLASH.

https://doi.org/10.1371/journal.pone.0169265.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Single subject acquisition of χ maps acquired with (a) ME-MP2RAGE, (c) ME-FLASH and (e) their difference, and test-retest reproducibility evaluation with voxel-by-voxel correlation 2D histograms for (b) ME-MP2RAGE, (d) ME-FLASH, and (f) ME-MP2RAGE versus ME-FLASH.

https://doi.org/10.1371/journal.pone.0169265.g008

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of the group averages μ_g and SDs σ_g results of the correlation parameters for the T₁, and χ maps acquired during Study 2.

https://doi.org/10.1371/journal.pone.0169265.t002

The 2D histograms of voxel-by-voxel correlations of T₁ maps underline that the acquisition scheme produces very consistent results, both for within- as well as across-sequence comparisons. This is also supported by the mean and SDs of the image volumes’ differences, μ_D and σ_D, being close to zero (relative to the T₁ value range). However, some deviation is observed for very long T₁ values corresponding to CSF voxels, as anticipated, because the selected MP2RAGE acquisition parameters were expected to systematically underestimate these values (see also Figs 2 and 6 and §2.3)

For the and χ maps, the 2D histograms and voxel-by-voxel correlations indicate that both ME-MP2RAGE and ME-FLASH do not achieve a similar degree of global reproducibility as observed for the T₁ maps. While both and χ maps are essentially compatible across sequence types with μ_D values close to zero, the σ_D values are not negligible and, at least in the case of , of almost the same order as the value range. Also, the bias introduced in the QSM results by the pole artifacts in the phase images seems to be independent of the acquisition scheme (see also S9 Fig). However, this does not seem to indicate inaccuracy of a specific acquisition scheme as the reproducibility does not substantially change for the within-sequence comparison. These results were consistently observed for all subjects.

The inter-acquisition averages of μ_D, μ_|D|, σ_D, σ_|D| for the test ME-MP2RAGE acquisition, for the reference MP2RAGE and ME-FLASH acquisitions, and for the cross-comparisons are always within twice their respective SDs. Additionally, the r² parameters for different acquisition schemes of the same mapping are always within a SD of each other. This suggests that the accuracy, precision and reproducibility of the mappings are only marginally dependent on the tested acquisition schemes.

The results from the ROI-based analysis for all maps are presented in Table 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Group averages, μ_g, and SDs, σ_g, for the ROI-based analysis of ME-MP2RAGE acquisitions.

Abbreviations: WM = White Matter; Fro. = Frontal Lobe; Tem. = Temporal Lobe; Par. = Parietal Lobe; Occ. = Occipital Lobe; Ins. = Insula; Put. = Putamen; Caud. = Caudate; Tha. = Thalamus; Cereb. = Cerebellum. Susceptibility is indicated here by Δχ as a reminder of the values being referenced to an arbitrary offset, implying that the group average and SD values are expected to be biased by that and, possibly, by the pole artifacts of the phase images.

https://doi.org/10.1371/journal.pone.0169265.t003

3.3 Effects from Mis-registration and Noise

The effects of registration are reported in Fig 9, where we show the 2D correlation histograms both on artificially mis-registratered and on self-registered T₁, , and χ maps, respectively. The original maps were obtained with MP2RAGE for T₁ and with ME-FLASH for and χ in these cases. Additional results for correlation histograms and the statistical comparisons related to the mis-registration are provided in Simultaneous S3–S5 Figs and S11–S13 Tables, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. 2D correlation histograms showing the effects of registration.

Columns refer to T₁ (a, c), (b, d) and χ (e, f) maps, while rows show either mis-registration (a, b, c) or self-registration (d, e, f) effects, respectively. The mis-registration and the self-registration effects are here illustrated by comparing each map with itself: after the application of a small transformation without registration (mis-registration) or after the application of a large transformation followed by a registration step (self-registration). The x-axis indicate the unmodified map, while the y-axis indicate the transformed or self-registered map. The small tranformation considered for the misregitration is a 0.5px translation followed by a 0.5° rotation. The large transformation considered for the self-registration is a 10° rotation.

https://doi.org/10.1371/journal.pone.0169265.g009

The 2D histograms obtained from these analyses are comparable to those obtained from the comparison of different acquisitions, even for relatively small transformation parameters. These essentially capture the global effects of the considered transformations. Moreover, the correlation coefficients indicate different sensitivities for the different map types. Specifically, a relatively linear behavior of r² for the explored range of parameters is found for T₁ and χ maps, while there is a remarkably non-linear drop in r² for the same range of transformations in the case of (e.g., r² dropped from 0.952 to 0.918 for T₁ and from 0.913 to 0.850 for χ but from 0.615 to 0.401 for in analyses of subtle translations of 0.375 px and 0.500 px). The comparisons of the “self-registered” maps indicate that the effect of the registration step on the correlation analyses is at least of the same order of magnitude of a 0.2° rotation, which is alone sufficient to explain an r² of ≈ 0.99, 0.91 and 0.97 for T₁, and χ respectively, the major effect again being observed with maps. Note that self-registration results indicate only a lower limit on the effects of the registration step because of the negative effects of noise. Additionally, neither r² nor the other correlation parameters are linear with respect to SNR, mis-registration or their combination; therefore these results cannot be used to quantify the single effects.

The explored noise level range correspond to up to about half of the SNR of the original images. The values of the obtained coefficients indicate that such SNR reductions degrade the voxel-by-voxel correlations with an effect size that is comparable to that from a translation of approximately 0.5 px, a rotation by approximately 0.9°, or a 0.4° rotation followed by a 0.4 px translation. An example of the effects of an increased noise level on χ maps is reported in S6 Fig. A summary of the correlation analysis is provided in S14 Table.

4.1 Acquisition Parameters

In this work, simulations and exploratory MRI experiments were employed to optimize an ME-MP2RAGE sequence modification for high-resolution simultaneous mapping of T₁, , and χ. The proposed set of acquisition parameters is specifically tailored to the tested resolutions and FOV. The software tools developed for this work can effectively be used for optimizing the acquisition parameters for different use cases. However, their limitations should be properly considered.

Particularly, the solutions of the Bloch equations for evaluating the sensitivity are ignoring the experimentally available SNR. This is of particular relevance when considering that the T₁ maps are more insensitive to inhomogeneities at lower (relative) , which is typically achieved by lowering the flip angles α₁ and α₂ or increasing the T_R,seq. However, this might have a negative impact on the acquisition time or the SNR (depending on T_R,GRE, which in turn determines the Ernst angle).

The fitting procedures that were used to infer the influence of SNR and sampling points on the accuracy and precision of the estimated relaxation time constants (i.e., T₁ from T_I,(1,2) and from different T_E,i) are based on the implicit assumption that the signal is well modeled by a mono-exponential function. However, recent works suggest that more complex curves (e.g., multi-exponential) might be required for modeling both longitudinal [12] and transverse [25, 46, 47] magnetization behaviors. More elaborate acquisition schemes with many more support points would be required to further investigate this aspect, which is neglected at this level.

Additionally, while the physical experiment is fundamentally identical, the methods for obtaining T₁ are slightly different for the simulation and the measurements. While the simulations rely on the fitting of magnitudes images, the measurements use a combination of the complex images to obtain the MP2RAGE signal which is then converted to T₁ via a look-up table (see also: Fig 2). In fact, the MP2RAGE approach is more accurate [5], essentially because the phase images provide additional information, particularly for the first inversion image, where the magnitude-only images have very low intensity due to the proximity to the zero-crossing in the inversion recovery curve. Therefore, the simulation constitutes a lower limit (except for SNR considerations) to the accuracy of the T₁ estimates.

The proposed procedure for the choice of the acquisition parameters requires a manual optimization, which inherently implies a certain degree of arbitrariness. While a more formal approach would be possible, there are several aspects that should be considered. Firstly, the arbitrariness is inherent to the problem, because a particular desired feature (e.g., resolution, accuracy, speed) is obtained at the expenses of the others. For example, it would be possible to write a cost function for insensitivity and a cost function for the coverage of a specific range of T₁ values. However, there is no unbiased way of combining the two without introducing some degree of arbitrariness, and further weights would need to be integrated into the model in order to cover the most relevant aspects of the acquisition. Such weights might not have immediate physical links with the desired features, making it difficult to find a rationale for their definition. Secondly, given the structure of the manifold parameters to be explored, there exist several suboptimal combinations as already suggested by [5], and an optimization algorithm might not find the best solution. Although a metric could be defined to describe the marginal gain associated with slight changes in parameters (e.g., by taking the partial derivatives of “first order” cost functions with respect to the parameters optimized by the other cost functions), this would lead to a much more complex mathematical model and would introduce additional arbitrary weights. Lastly, given the limits to specify sequence parameters on the protocol level (e.g., the flip angle can only be set to integer values in our current implementation), an automated optimization procedure is likely to produce a set of parameters that is not accessible with the same exactness in an experimental situation. Ultimately, this does not justify the additional complexity of such an approach. For these reasons, we performed the simulations by choosing the parameters with a rationale based on provisions of the desired resolution and scan time. These results were then extensively tested experimentally, thus providing additional information on the impact of aspects that were neglected or only partially addressed by the simulations.

4.2 Simultaneous Parameter Mapping

Focusing on the voxel-by-voxel comparisons of results obtained with different or the same acquisition modalities, a major point to discuss is the validity of the method for assessing the impact of the acquisition technique on the map accuracy. Since there was no direct access to the “ground truth”, we assumed that a more-established acquisition technique can be used as a “gold standard” for a specific map. Based on this, we consider two different acquisition schemes to be experimentally equivalent if they have the same degree of reproducibility (e.g., similar μ_D and σ_D) as two separate repetitions with the same sequence. This rationale relies on the assumption that the assessment of the reproducibility is unbiased. However, the registration step has an important effect on the parameters used to quantify reproducibility. This bias is always present (for non-simultaneous acquisitions), arising from the numerical interpolation required to address geometrical transformations acting on a sub-voxel scale—even in ideal cases where other spatial inaccuracies due to movements of the subject during the acquisition or noise associated with physiological processes could be ignored. While errors arising during registration are not sufficient to entirely explain the observed level of reproducibility, our simulation results are compatible with the hypothesis that even subtle mis-registrations are able to largely explain the values of the statistical coefficients used to assess the reproducibility. This is especially true for the measured (and, partially, also for χ) maps, because the inter-voxel values fluctuations are observed at a much finer spatial scale compared to the measured T₁ maps. These differences are mostly localized at the boundaries of regions characterized by different average values (e.g. GM/WM, GM/CSF, etc.) and are sufficiently strong to be captured by the global coefficients and plots considered.

The lower reproducibility that is observed for the and χ (compared to T₁) may not entirely be explained by registration issues, and further investigation may be required to elucidate other potentially relevant sources of error, for example physiological noise or motion. This would improve the understanding of the link between the tissue microstructure and mapping at a finer level. To reduce the effect of noise in these maps, it would be possible to introduce in the fit algorithm a regularization term incorporating additional information, for example from neighboring voxels [12, 46] or from less noisy acquisitions like the simultaneously measured T₁.

Regardless of the acquisition scheme, our results also indicate a higher reproducibility of χ results compared to . However, this result would need a more robust validation. In fact, the correlation coefficients considered in this work are not appropriate for accurate quantification of this aspect. Additionally, the QSM results are obtained on a smaller volume (due to limitations of the processing pipeline) and, hence, a comparison of both contrasts might be biased by this difference. Note also that values were generally longer in the voxels excluded for QSM and might be estimated with an increased error because the acquisition parameters had been optimized for values up to ≈ 45 ms.

The SNR available for and χ estimates was inherently lower for ME-MP2RAGE as compared to ME-FLASH with the tested parameters. This is because substantially longer T_R values and higher flip-angles were used for ME-FLASH. Additionally, the SNR of the second GRE block in (ME-)MP2RAGE is inversely modulated by T₁ relaxation. Despite these differences in the SNR, the reproducibility coefficients are not appreciably different across acquisition schemes, suggesting that either the mapping or the comparison procedure (or both) are not sensitive enough for subtle SNR variations.

Of note, when choosing T_E for both ME-FLASH and ME-MP2RAGE, we ignored potential effects related to the phase difference between water and fat. This was not a significant issue because usually the fat contribution to the signal in brain tissues is negligible, since most lipid molecules are immobilized in membranes and will have a in the μs range, which is too short to be observed with the typical echo times of a standard ME-FLASH experiment. Otherwise, we would have expected an acquisition-related bias in the maps, which would have caused, for example, the μ_D parameter between different acquisition to appreciably deviate from zero, and this was not observed.

The T₁ values obtained from the ROI-based analysis indicate excellent agreement with [10] for WM, while GM values are consistently lower in our experiments. This might be explained by a slightly different definition of the ROIs and varying partial voluming effects, especially with CSF. values are generally in agreement with previously published results (e.g., [48, 49]), although it should be noted that results are always modulated—and therefore biased—by the B₀ shimming. The susceptibility results are relative to an arbitrary reference, which was chosen to obtain an average susceptibility value very close to zero. For this reason, the group-averages in this case may be of limited significance.

In its current form, the ME-MP2RAGE pulse sequence achieves simultaneous acquisition of T₁, , and χ maps as a trade-off between resolution, SNR, and scan time. Our set of parameters (Table 1) showed that a slightly reduced SNR to stay within the constrained acquisition time was still sufficient for relatively accurate simultaneous mapping. A more significant impact on the map accuracy, however, is expected for more aggressive reductions of the scan time (e.g., by further increasing the GRAPPA factor). Instead, if the resolution and FOV constraints are relaxed, it is possible to reach the point where the MP2RAGE and its ME variant have identical scan times, but the latter allows the simultaneous acquisition of T₁, and χ maps without a further penalty for the achieved accuracy. This situation might be more interesting for 3 T, where typical resolutions are lower. This aspect remains to be experimentally investigated, but the already mentioned simulation tools can be effectively used to predict the reliability of T₁ estimates.

Further sequence developments may allow to address some of the current limitations of ME-MP2RAGE. A finer control over timing and k-space coverage will improve the flexibility of the sequence, thus allowing for better trade-offs in resolution, scan time, and maps accuracy. This could be achieved for example by mixing the first and the second phase encoding directions in each GRE block (eventually integrating compressed sensing techniques) and/or by using a different number (and eventually timings) for the acquired echoes across the GRE blocks, but this would require both major modification to the sequence code and careful consideration on the effects on the T₁ point-spread function, and ultimately a separate study to evaluate its accuracy in estimating T₁ maps.

Such modifications may also lead to the opportunity of acquiring more sampling points in the inversion recovery curve. This is interesting because it might offer an effective and time saving sequence for the purpose of investigating non-mono-exponential behaviors, and hence partial volumes effects.

Another approach toward the acquisition of multiple images at different inversion times was recently proposed with the MPnRAGE sequence [50], where a radial readout scheme is employed instead of the cartesian one.