Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

A Novel Multiparametric Approach to 3D Quantitative MRI of the Brain

  • Giuseppe Palma ,

    giuseppe.palma@ibb.cnr.it

    Affiliation Institute of Biostructure and Bioimaging, National Research Council, Naples, Italy

  • Enrico Tedeschi,

    Affiliation Department of Advanced Biomedical Sciences, University of Naples Federico II, Naples, Italy

  • Pasquale Borrelli,

    Affiliation IRCCS SDN, Naples, Italy

  • Sirio Cocozza,

    Affiliation Department of Advanced Biomedical Sciences, University of Naples Federico II, Naples, Italy

  • Carmela Russo,

    Affiliation Department of Advanced Biomedical Sciences, University of Naples Federico II, Naples, Italy

  • Saifeng Liu,

    Affiliation The MRI Institute for Biomedical Research, Waterloo, ON, Canada

  • Yongquan Ye,

    Affiliation Department of Radiology, Wayne State University, Detroit, MI, United States of America

  • Marco Comerci,

    Affiliation Institute of Biostructure and Bioimaging, National Research Council, Naples, Italy

  • Bruno Alfano,

    Affiliation Institute of Biostructure and Bioimaging, National Research Council, Naples, Italy

  • Marco Salvatore,

    Affiliation IRCCS SDN, Naples, Italy

  • E. Mark Haacke,

    Affiliations Department of Radiology, Wayne State University, Detroit, MI, United States of America, The MRI Institute for Biomedical Research, Detroit, MI, United States of America

  • Marcello Mancini

    Affiliation Institute of Biostructure and Bioimaging, National Research Council, Naples, Italy

Abstract

Magnetic Resonance properties of tissues can be quantified in several respects: relaxation processes, density of imaged nuclei, magnetism of environmental molecules, etc. In this paper, we propose a new comprehensive approach to obtain 3D high resolution quantitative maps of arbitrary body districts, mainly focusing on the brain. The theory presented makes it possible to map longitudinal (R1), pure transverse (R2) and free induction decay () rates, along with proton density (PD) and magnetic susceptibility (χ), from a set of fast acquisition sequences in steady-state that are highly insensitive to flow phenomena. A novel denoising scheme is described and applied to the acquired datasets to enhance the signal to noise ratio of the derived maps and an information theory approach compensates for biases from radio frequency (RF) inhomogeneities, if no direct measure of the RF field is available. Finally, the results obtained on sample brain scans of healthy controls and multiple sclerosis patients are presented and discussed.

1 Introduction

Multi-parametric quantitative Magnetic Resonance Imaging (qMRI) based on the relaxometry properties of intracranial tissues has long been and still is an active field of research in medicine and physics [1]. Several approaches have been used, aiming at obtaining quantitative estimates of the longitudinal relaxation rate (R1), transverse relaxation rate (R2) and proton density (PD) of brain tissues [27]. Other relaxation parameters whose MR signals can be usefully exploited to discriminate the microstructure of brain tissues are the free induction decay (FID) rate () [4] and the magnetic susceptibility (χ) [8]. In this paper, we present a new acquisition and processing procedure that, starting from a set of conventional high resolution 3D Steady State sequences, makes it possible to achieve full brain coverage for absolute measurements of intracranial compartments. For the first set of tissue properties, Radio Frequency (RF) B1 inhomogeneities can create problems. As pointed out in [9], a 3D acquisition protocol is preferred to remove intra-voxel biases, known to arise from imperfect 2D RF-pulse profiles that cause different isochromat evolutions in response to different effective flip angles.

Here, we prove that a set of steady-state sequences, acquired with variable flip angles and different phase coherence, makes it possible to derive quantitative volumetric R1, R2, PD, and χ maps of the brain tissues. We also show that high Signal-to-Noise Ratio (SNR) full brain coverage with a sub-millimeter resolution may be obtained in a total acquisition time of 14 minutes. The proposed solution for the relaxation rates is fully analytic and allows for the inversion of the signal equations for R1, R2, PD and maps in a typical 3D dataset within a few seconds (5 s for a 320 × 270 × 128 voxel array on a 2.7 GHz Intel Core i5 processor).

The plan of the paper is as follows. In §2, the Bloch equations of the sequences are briefly reviewed; then, the proposed mathematical scheme to derive R1, R2, PD, and χ maps is presented in full detail. In §3, we provide a description of a sample experimental setup performed to prove the feasibility of the general theoretical scheme of §2. Finally, results are presented in §4 and discussed in §5.

2 Theory

2.1 Sequence Bloch equations

Within the steady-state sequence family, we selected the spoiled Gradient Echo (GRE) and the balanced Steady-State Free-Precession (bSSFP) sequences because of their low sensitivity to flow artifacts in both CSF in ventricles and blood within vessels.

The complex signal of a spoiled GRE sequence is (1) where K is the complex coil sensitivity, M0 is the equilibrium magnetization, θ is the flip angle, E1 ≡ exp(−TRR1) (TR being the repetition time), TE is the echo time, γn is the gyromagnetic ratio of the imaged nucleus, ΔB is the local magnetic field variation and ϕ0 is the phase-shift induced by the RF-pulse.

The complex signal of a bSSFP sequence is (2) where (3) (4) (5) E2 ≡ exp(−TRR2), φ = γnΔBTR + Δφ is the resonance offset angle and depends on how the phase-cycling Δφ is achieved: if the frequency of the Larmor reference frame is formally shifted, then ; otherwise, if a constant phase term is added to the previous RF phase each time a new RF pulse is applied, then .

2.2 Band-free bSSFP reconstruction

Due to the dependence of (Eq 2) on φ, bSSFP images are well known to be affected by severe banding artifacts, particularly in regions with rapidly varying susceptibility or main magnetic field [9]. To avoid a biased R2-map, in the following a new dataset S0 (independent of φ) is derived from a set of n bSSFPs (n ≥ 3) acquired with increasing resonance offset angles (phase cycling) as described in [10, 11].

The φ-independent contrast S0 we aim to derive corresponds to (6) For each phase-cycling j ∈ {1, …, n}, we can formally introduce (7) (8) (9) (10) Defining (11) (12) from Eqs (2)–(10) it can be shown that (13)

If n ≥ 3, the Moore-Penrose pseudoinverse of A provides the least-squares estimate of the unknown variables α, β and γ (cast into x form) as a function of the bSSFP datasets (cast into y form—each is known and corresponds to the acquired scan multiplied by a function of the acquisition parameters): (14)

From Eqs (8)–(9), it follows that (15) whence, from (Eq 8), we obtain (16)

2.3 R1-map

The acquisition of two or more spoiled GRE scans performed at different flip angles and fixed TR makes it possible to estimate R1 [9]. The magnitude of (Eq 1) satisfies (17) (Eq 17) represents the equation of a straight line with slope E1 in a plane whose coordinate pairs are (SFL/tan θ, SFL/sin θ). Therefore, if we write the signal intensity of the j-th echo (j ∈ {1, …, m}, m ≥ 2) of the i-th spoiled GRE acquisition (i ∈ {1, …, l}, l ≥ 2) as (18) where (19) R1 can be estimated via simple linear regression of (Eq 17) as (20) with (21)

2.4 -map

If we define (22) from (Eq 18) it follows that (23) whence turns out to be the linear coefficient between and TE,j. A simple linear regression yields therefore to (24)

2.5 PD-map

Once both R1 and maps are known, it is also possible to provide an estimate of proton density (in arbitrary units) by inverting (Eq 18) for each i, j with respect to |K|M0 from (Eq 19), and averaging the results over the acquired echoes: (25)

2.6 χ-map

An extensive description of the susceptibility quantitation techniques is reported in [8]. Here we briefly summarize the actually adopted strategy.

The net magnetic field B(r) resulting from the magnetization M induced within a matter distribution by an external magnetic field B0 is given by (26) For the linear magnetic materials, (Eq 26) is completed by the magnetostatics constitutive equation: (27) which reduces to (28) for diamagnetic and paramagnetic substances whose susceptibility |χ| ≪ 1.

As in the MRI context , the component of the field variation due to susceptibility inhomogeneities is therefore given by (29) (30) where (31)

As we acquire a spoiled GRE sequence, the phase associated to susceptibility (see (Eq 1)) is therefore given by (32)

In order to extract ϕχ(r) from the original wrapped phase ψ(r) (33) (which also depends on background magnetic field inhomogeneities, the imaginary coil sensitivity ℑ(K), RF-pulse design etc.), phase unwrapping was performed using the 3D-SRNCP algorithm [12], followed by the SHARP filter [13], which exploits the harmonicity of static magnetic fields within homogeneous media (∇2(ϕϕχ) = 0) [14]: (34)

The inversion of (Eq 32) requires a regularization of the inverse of g(k) ≡ FT[G(r)] (see [8, 15]), whence we finally obtain (35)

2.7 R2-map

Once the R1- and the PD-maps are known, (Eq 6) makes it possible to derive the R2-map. From the following combination of Eqs (3), (5) and (6) (36) if the acquired bSSFP scans are set with TE = TR/2, it follows that (37) where (38) (39) (40) (41)

From the standard theory of the quartic equations it follows that the physically meaningful solution of (Eq 37) is given by (42) where (43) (44) (45) (46) (47)

The transverse relaxation rate is then found from the expression (48)

3 Materials and Methods

3.1 Study design

We set up an experimental configuration to test the practical feasibility of the above general theoretical algorithm (schematically shown in Fig 1). It should be noted that, depending on the specific features of the available MR system, different strategies may be adopted.

thumbnail
Fig 1. Schematic flowchart of the implemented algorithm.

From FC-FLASH phase data the χ-map is directly derived (see §2.6). The TrueFISP collection is used to extract the band-free S0 dataset (see §2.2). Then, FC-FLASH and S0 series are denoised by the SVN-MNLM scheme (see §3.3). The -map is computed from resulting FC-FLASH series (see §2.4). The introduction of an estimate of the B1-map (actually measured or simply guessed) allows for the extraction of the R1- (see §2.3), PD- (see §2.5) and R2- (see §2.7) maps. Depending on the reliability of the B1-map estimate, a B1-correction scheme may be iteratively applied (dashed lines—see §3.4).

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

First, we describe the acquisition protocol; then we present the image restoration algorithm we applied to increase the Signal to Noise Ratio (SNR) of the datasets. Further, as on our scanner there is no protocol for mapping B1 sensitivity profile, we describe the method developed to account for non-ideality of the profiles within the Field of View (FOV). Finally, we detail the way devised to estimate the overall reproducibility of the results.

3.2 Acquisition protocol

MR data were collected on a 3 T scanner (Trio, Siemens Medical Systems, Erlangen, Germany) with a volume transmitter coil and 8-channel head receiver coil. The “Carlo Romano” ethics committee for biomedical activities of “Federico II” University of Naples (Italy) specifically approved the study and the written informed consent form, which was signed by the subjects (one Healthy Control (HC) and one Relapsing-Remitting Multiple Sclerosis (MS) patient) undergoing the MR scan.

The 3D Steady-State sequences were acquired with a pseudo-axial orientation (along the anterior commissure-posterior commissure line) and shared the same FOV covering the whole brain (230 mm in anteroposterior direction; 194 mm in laterolateral direction; 166 mm in craniocaudal direction) with an in-plane resolution of 0.65 mm (laterolateral phase encoding direction) and a slice thickness of 1.3 mm (voxel size ΔV = 0.55 mm3).

We acquired 2 fully flow-compensated double-echo spoiled GRE (FC-FLASH) sequences with flip angles of 3° and 20°, repetition time of 28 ms, echo times of 7.63 ms and 22.14 ms, a GRAPPA factor of 2 and a bandwidth of 190 Hz/pixel (acquisition time: 2 × (4’ 46”)). For each echo, a magnitude/phase reconstruction was enabled, thus obtaining a total of 4 complex volume datasets.

Then, a set of 4 phase-cycled bSSFP (TrueFISP) sequences was acquired with resonance offset angle steps of π/2 (i.e. Δφ ∈ {0, π/2, π, 3π/2}), flip angle of 10°, repetition time of 7.3 ms, echo time of 3.65 ms, a GRAPPA factor of 2 and a bandwidth of 252 Hz/pixel (acquisition time: 4 × (1’ 20”)). For each echo, a magnitude/phase reconstruction was enabled, thus obtaining a total of 4 complex volume datasets (i.e. one for each resonance offset angle).

It may be noted that, depending on the parameter subset of actual interest, this acquisition protocol may be shortened following the recipes shown in Table 1.

thumbnail
Table 1. Design guidelines for minimal acquisition protocols (second column) making it possible to reconstruct each parameter listed in the first column.

The corresponding total acquisition time and the list of other maps that can be obtained for free are reported in the third and fourth columns, respectively.

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

3.3 Image restoration

In order to control the noise propagation in the relaxometry maps derived by inversion of the Bloch equations, the FC-FLASH and S0 image series were pre-processed with a novel ad hoc version of the Non-Local Means (NLM) filter, adapted for parallel MRI (see [16]) and multi-spectral images (see [17]) with an improved handling of the arbitrary noise power in each image component.

Spoiled GRE ({Si,j}) and bSSFP (S0) datasets were considered as a discrete version of a general multi-component image X: ℝN → ℝM (N = 3 and M = 5) with a bounded support Ω ⊂ ℝN. In this sense, the Spatially-Varying-Noise Multi-component NLM (SVN-MNLM) is a class of endomorphisms of the image space, each identified by 2 parameters (ρ and ς), that act as follows: (49) where (50) (51) (52) is the similarity radius, ς ∼ 1 is an adimensional constant to be manually tuned that determines the filter strength and is the standard deviation of noise of the mth image component at (the noise power maps were estimated following [16]).

Due to the high computational complexity of the above scheme, a multi-Graphics Processing Unit (GPU) implementation of the NLM algorithm [18] was adapted to (Eq 49) and exploited for fast image processing.

3.4 -inhomogeneity correction

The derivation of R1, PD and R2 maps in §2.3, §2.5 and §2.7 critically depends on the flip angle sensed by each voxel all-over the FOV. However, transmission coil design issues, the non-ideality of 3D slab profiles and the dielectric resonances, which are more prominent on high field scanners, may appreciably alter the homogeneity, thus resulting in a low frequency variation of the derived relaxometry parameters as we move from internal to peripheral brain regions. Moreover, the PD map also depends on the map of the receiver coil, while the other maps do not depend on the receiver coil sensitivity, as this constant factor is eliminated in the ratios of signal intensities in Eqs (20), (24) and (36).

Taking advantage of the large scale nature of the variation pattern, we adopted an information theory approach to remove the biases due to FA inhomogeneities. Essentially, a useful approximation of the peripheral FA decay is given by its second order Taylor logarithmic expansion about its stationary point , modulated by a slab profile term: (53) (ζ is the offset from the slab center in the slab select direction). As the Hessian tensor H+ is symmetric, is a vector and k, ξ ∈ ℝ, a total of 11 free parameters need to be estimated in (Eq 53).

Similarly, the map can be estimated by a different second order Taylor logarithmic expansion (this time no slab profile affects the receiver coil sensitivity, thus reducing the number of free parameters to 9): (54)

Finally, the parameters are estimated by minimizing, with the Nelder-Mead method, the entropy of the 3D joint histogram related to R1, PD and R2 maps in low gradient regions of the maps.

The effectiveness of the obtained B1 inhomogeneity correction was evaluated by acquiring the full set of sequences on a cylindric homogeneous phantom (diameter: 18 cm; height: 30 cm) and computing the relative Full Width at Half Maximum (FWHM) of the R1 and R2 distributions reconstructed with and without the -inhomogeneity correction.

3.5 Reproducibility estimation

To estimate the confidence interval of the proposed maps (taking into account image uncorrelated noise, variations in signal amplification of the MR scanner, tissue temperature fluctuations, etc.), each FC-FLASH was acquired twice and each TrueFISP three times in the HC during the same exam session, without moving the subject. Given the redundant set of complex datasets, an ensemble of 24 (spoiled GRE) × 3 (S0) different complete relaxometry protocols was produced. Forty-eight different estimates of the relaxometry maps were thus derived and used to estimate the degree of reproducibility via mean and standard deviation maps.

3.6 Image processing and analysis

All data processing was performed using an in-house developed library for Matlab (MATLAB® Release 2012b, The MathWorks, Inc., Natick, MA, USA), partly described in previous works [1618], on a commercial workstation (Intel® Core i7-3820 CPU @ 3.6 GHz; RAM 16 GB) equipped with 2 GPU boards (NVIDIA GeForce® GTX 690).

The definition of the brain structures was qualitatively assessed in consensus by three experienced neuroradiologists, who were asked to indicate the specific diagnostic contribution provided by each map.

To further assess the capability of the quantitative parameters to distinguish between different intracranial structures, a set of bilateral ROIs (see Fig 2) was manually drawn on the R1-maps (and then automatically transferred to the other maps) of both the HC and MS subjects. ROIs were intentionally drawn well within the anatomical structures to avoid partial volume effects.

thumbnail
Fig 2. Axial brain slices at the 4 different levels showing the position of ROIs drawn for measurement of mean relaxometry values in selected brain structures: cortex, white matter and CSF (a); head of caudate and thalamus (b); globus pallidus and putamen (c); red nucleus and substantia nigra (d).

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

4 Results

4.1 Effectiveness of -inhomogeneity correction

The test of the effectiveness of the proposed -inhomogeneity correction on R1 and R2 is reported in Table 2.

thumbnail
Table 2. Percentual FWHM of the R1 and R2 distribution obtained without (pre) and with (post) -inhomogeneity correction.

For each case, we report the values for the entire phantom and for its eroded version (defined by a spherical structuring element with radius of 1 cm) simulating the brain without the skull.

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

4.2 Map assessment

R1, R2, PD, and χ maps were obtained for each subject in about 30 minutes after completion of the acquisition. Examples of the maps for the 43-year-old female HC are displayed in Fig 3. Thanks to the acquired voxel size, the reconstructed datasets allow for reasonably good multi-planar reconstructions.

thumbnail
Fig 3. Axial (left), coronal (middle) and sagittal (right) brain slices of R1 ([0 ∼ 2] s-1), R2 ([0 ∼ 50] s-1), ([0 ∼ 50] s-1), PD ([0 ∼ 1] arbitrary units) and χ ([−300 ∼ 300] ppb) maps (from top to bottom) in a 43-year-old female HC.

As coronal and sagittal images are derived by a multi-planar reconstruction of the original axial dataset, their in-plane resolution is 0.65 × 1.3 mm2, with a slice thickness of 0.65 mm. Small insets pointing out the position of the slices on a perpendicular plane are shown as anatomical reference in the upper row.

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

In Fig 4 we report the mean quantitative maps and their corresponding voxel-by-voxel standard deviation maps, obtained as described in §3.5. On the whole, the standard deviation ranged from 2% to 4% of the quantitative absolute values, showing an excellent reproducibility of the map estimates.

thumbnail
Fig 4. R1 ([0 ∼ 2] s-1), R2 ([0 ∼ 50] s-1), ([0 ∼ 50] s-1), PD ([0 ∼ 1] arbitrary units) and χ ([−300 ∼ 300] ppb) (from top to bottom) maps (left) and corresponding confidence interval (right) in a 43-year-old female HC, displayed with a markedly different (scale factor of 20) grayscale to highlight tiny differences in measure uncertainties.

A small inset pointing out the position of the slices on a sagittal plane is shown as anatomical reference in the upper row.

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

4.3 Morphological evaluation

Fig 5 is an example of the maps assessed by the neuroradiologists in the qualitative analysis. In both HC and MS patient the R1-map was rated as the most representative of the brain anatomy, with excellent differentiation between grey (GM) and white (WM) matter. The - and χ-maps clearly showed the subcortical structures known to induce significant local field inhomogeneities due to build-up of paramagnetic substances. In addition, in the MS patient, demyelinating lesions were clearly displayed on both R1- and R2-maps, even in critical areas, such as the periventricular WM (Fig 6).

thumbnail
Fig 5. Axial brain slices at the level of midbrain (left), basal ganglia (middle) and fronto-parietal convexity (right) of R1 ([0 ∼ 2] s-1), R2 ([0 ∼ 50] s-1), ([0 ∼ 50] s-1), PD ([0 ∼ 1] arbitrary units) and χ ([−300 ∼ 300] ppb) maps (from top to bottom) in a 43-year-old female HC.

Small insets pointing out the position of the slices on a sagittal plane are shown as anatomical reference in the upper row.

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

thumbnail
Fig 6. R1 (a), FLAIR (b), R2 (c), PD (d), (e) and χ (f) axial brain slices at the level of the lateral ventricles in a 45-year-old female MS patient (Expanded Disability Status Scale—EDSS: 3.5; disease duration: 7 years).

A small inset pointing out the position of the slices on a sagittal plane is shown as anatomical reference in the upper row.

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

The mean and standard deviation of the quantitative values of each structure and map are reported in Table 3.

thumbnail
Table 3. Relaxometry and susceptibility properties measured in selected brain locations.

In each cell, ROI mean and standard deviation are reported.

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

5 Discussion

We have described a 3D acquisition protocol that yields multiple quantitative parametric maps (R1, R2, , PD and χ) based on the relaxometric properties of brain tissues.

Image quantitation in MRI has been widely studied in the last decades, with particular regard to tissue relaxometry, and almost any subset of the quantitative parameters obtained with our approach has been the subject of at least one study [1, 8, 9, 19]. However, a great majority of the proposed schemes largely rely on 2D acquisition sequences [3, 2024], which, in this context, may be detrimental for at least two reasons. First, SAR issues limit the slice thickness to values considerably larger than the voxel thickness achievable in 3D acquisitions. Second, and most important, voxel intensities reconstructed from 2D acquisitions derive from the integration of isochromat signals over a non-ideal slice profile. Their Bloch equations in fact evolve according to flip angles that greatly vary over each voxel, unlike 3D sequences, in which intra-voxel flip angle variation is usually small compared to its mean value [9]. Therefore, 2D sequences usually lead to biases in R1 and R2 calculations, which are particularly hard to estimate if the RF-pulse is not completely known or for sequences leading to entangled relaxometry equations (e.g., any type of SSFP [25]).

Conversely, other groups have adopted 3D acquisition strategies to recover R1- and -maps of the brain [6, 23], but not the quantitation of pure transverse relaxation rate (R2) or susceptibility (χ). These may provide useful information when R1 and fail to provide adequate information or contrast to reveal the pathophysiological conditions of specific areas (e.g. in presence of high susceptibility gradients near petrous bones, nasal sinuses, etc. that cause an enhancement uncorrelated to any tissue alteration).

Recently, some 3D approaches have been proposed to derive R1- and R2-maps based on different modes of unbalanced SSFP signals [5, 7], with promising results in the study of joints, cartilage and, in general, of the musculoskeletal system [26, 27]. However, 3D unbalanced SSFP sequences are known to be highly sensitive to pulsatile flow of fluid with relatively long T2, thus making the brain study unfeasible, due to CSF pulsation, unless switching to 2D [28].

In this respect, the family of DESPO methods [2, 2931], while neglecting - and χ-maps, do provide a full 3D high SNR R1- and R2-maps [2], using sequences poorly prone to motion artifacts (spoiled GRE and bSSFP [32, 33]). Moreover, the numerical approach described in [31] is able to subtract the banding artifacts due to B0 inhomogeneities sensitivity of bSSFP. However, some drawbacks of DESPOT2 may raise crucial issues even if and susceptibility information are not of interest. In particular, even if bSSFP sequences are usually quite insensitive to motion, when flip angle gets close to 90 degrees, CSF pulsatile flow may cause some periventricular hyperintensities due to ghosting artifacts in phase encoding directions. This, in turn, leads to a biased estimate of R2 in the periventricular WM that may mimic, for example, demyelinating lesions. Moreover, the same high flip angle values prescribed for at least one bSSFP of the DESPOT2 protocol are likely to exceed SAR constraints, especially if a short TR is desired to shorten total acquisition time.

The scheme we propose is meant to solve the above mentioned issues. As DESPO, our approach entirely relies on widely available 3D acquisition sequences. Moreover, it allows for the quantitation of 5 independent parameters and gets rid of the sensitivity to B0 inhomogeneity by means of a fully analytical solution (thus speeding up the computation step). As for the inhomogeneity dependence, one can completely remove it using a measured B1 field map, if an ad hoc protocol is available on the scanner; otherwise, as shown in §3.4 and §4.1, an information theory approach can be exploited to largely recover the map homogeneity. Furthermore, a judicious use of the Bloch equations for the acquired MR signals enabled us to skip the acquisition of high flip angle bSSFPs, thus limiting the required acquisition time and avoiding both SAR issues and CSF pulsation artifacts.

As for the quantitation of brain tissue parameters, the results are encouraging. From a pure theoretical point of view, it is noteworthy that the well known physical constraints on relaxation rates () are satisfied for any structure/tissue analyzed in Table 3, despite the highly uncorrelated acquisition strategies we adopted to estimate the associated maps. A comparison of ours with the corresponding estimated values present in the literature also confirms the reliability of the ranges we found in the different compartments. In particular, a good agreement is found for longitudinal relaxation (GM: 0.609 ± 0.008 s-1 in [6] or 0.622 ± 0.043 s-1 in [34]; WM: 1.036 ± 0.036 s-1 in [6] or 1.190 ± 0.071 s-1 in [34]; caudate nucleus: 0.683 ± 0.022 s-1 in [6] or 0.719 ± 0.025 s-1 in [34]) and transverse FID (GM: 15.2 ± 0.4 s-1; WM: 21.0 ± 0.8 s-1; caudate nucleus: 18.2 ± 1.2 s-1 in [6]) rates measured at 3 T. As for R2 values, which are usually less sensitive to B0 magnitude than R1 and , 3D measures performed at 1.5 T substantially match our findings (GM: 10.20 ± 0.73 s-1; WM: 18.5 ± 1.4 s-1; caudate nucleus: 11.24 ± 0.76 s-1 in [29]), except for an unusually high value found in WM by [29], which significatively differs also from several other 2D measures [3, 9, 20, 35]. Also, our PD estimates are in line with those reported in the literature (GM/CSF: 0.796 ± 0.078 in [20] or 0.844 ± 0.019 in [6]; WM/CSF: 0.742 ± 0.070 in [20] or 0.683 ± 0.006 in [6]; caudate nucleus/CSF: 0.851 ± 0.084 in [20] or 0.827 ± 0.016 in [6]). Finally, our χ measures fit well with the results of previous studies (GM: 16 ± 10 parts-per-billion (ppb); WM: −18 ± 9 ppb; globus pallidus: 200 ± 40 ppb; putamen: 70 ± 20 ppb; substantia nigra: 200 ± 60 ppb; red nucleus 90 ± 30 ppb in [13]). However, it should be borne in mind that plenty of variable factors (such as inter-subject variability, B0, voxel size, ROI definitions, scanner-specific biases, 2D or 3D imaging modality, etc.) may alter the qMRI results and, as such, should recommend caution in interpreting differences or similarities [6].

Visual inspection of our parametric maps also showed the expected features of normal structures (Fig 5) and of pathological tissue changes, such as demyelinating lesions (Fig 6).

We are aware that the present study is not without limitations. Due to the absence of a protocol for mapping the RF-transmission profile and the receiver coil sensitivity, we were forced to enable the -inhomogeneity correction introduced in §3.4, which can greatly limit the map inhomogeneities, but not completely compensate for it. Moreover, the underlying assumption of the present version of our scheme is that both longitudinal and transverse relaxations follow a mono-exponential recovery of the equilibrium condition. Even if this is a common assumption in relaxometry studies, there is evidence that WM shows multiple relaxation components, each with its own relaxing water pool [30]. In this respect, further theoretical and clinical analysis will be devoted to extend the present scheme to a multicomponent generalization.

In conclusion, we found a fully analytical solution to a 5 parameter qMRI problem that enables the reconstruction of 3D brain datasets with submillimiter resolution in a clinically feasible acquisition time (less than 15 minutes in our experimental setup). Our maps do not suffer from many of the issues affecting previously reported strategies, and allow for a more accurate characterization of brain structures. We believe this will trigger a convenient access to unconventional ways of studying a large spectrum of pathophysiological conditions.

Acknowledgments

Antonietta Canna and Federica Riccitiello are gratefully acknowledged for their valuable assistance when drafting the manuscript.

Author Contributions

Conceived and designed the experiments: GP. Performed the experiments: SC CR ET. Analyzed the data: GP PB SL. Contributed reagents/materials/analysis tools: PB SL YY. Wrote the paper: GP ET PB SC CR EMH. Designed the software used in analysis: GP PB MC BA. Obtained permission for MRI acquisition: ET MS MM.

References

  1. 1. Tofts P. Quantitative MRI of the brain: measuring changes caused by disease. John Wiley & Sons; 2005.
  2. 2. Deoni SC, Rutt BK, Peters TM. Rapid combined T1 and T2 mapping using gradient recalled acquisition in the steady state. Magnetic Resonance in Medicine. 2003;49(3):515–526. pmid:12594755
  3. 3. Warntjes J, Leinhard OD, West J, Lundberg P. Rapid magnetic resonance quantification on the brain: Optimization for clinical usage. Magnetic Resonance in Medicine. 2008;60(2):320–329. pmid:18666127
  4. 4. Draganski B, Ashburner J, Hutton C, Kherif F, Frackowiak R, Helms G, et al. Regional specificity of MRI contrast parameter changes in normal ageing revealed by voxel-based quantification (VBQ). Neuroimage. 2011;55(4):1423–1434. pmid:21277375
  5. 5. Palma G, Greco D, Innocenti S, Alfano B. Method of generating 2D or 3D maps of MRI T1 and T2 relaxation times. Google Patents; 2012. US Patent App. 13/409,953.
  6. 6. Weiskopf N, Suckling J, Williams G, Correia MM, Inkster B, Tait R, et al. Quantitative multi-parameter mapping of R1, PD*, MT, and R2* at 3T: a multi-center validation. Frontiers in neuroscience. 2013;7:95. pmid:23772204
  7. 7. Heule R, Ganter C, Bieri O. Triple echo steady-state (TESS) relaxometry. Magnetic Resonance in Medicine. 2014;71(1):230–237. pmid:23553949
  8. 8. Haacke EM, Reichenbach JR. Susceptibility weighted imaging in MRI: basic concepts and clinical applications. John Wiley & Sons; 2011.
  9. 9. Brown RW, Cheng YCN, Haacke EM, Thompson MR, Venkatesan R. Magnetic resonance imaging: physical principles and sequence design. John Wiley & Sons; 2014.
  10. 10. Björk M, Gudmundson E, Barral JK, Stoica P. Signal processing algorithms for removing banding artifacts in MRI. In: Proceedings of the 19th European Signal Processing Conference (EUSIPCO-2011), Barcelona, Spain; 2011. p. 1000–1004.
  11. 11. Björk M, Reeve Ingle R, Gudmundson E, Stoica P, Nishimura DG, Barral JK. Parameter estimation approach to banding artifact reduction in balanced steady-state free precession. Magnetic Resonance in Medicine. 2013;72(3):880–892. pmid:24166591
  12. 12. Abdul-Rahman HS, Gdeisat MA, Burton DR, Lalor MJ, Lilley F, Moore CJ. Fast and robust three-dimensional best path phase unwrapping algorithm. Applied optics. 2007;46(26):6623–6635. pmid:17846656
  13. 13. Schweser F, Deistung A, Lehr BW, Reichenbach JR. Quantitative imaging of intrinsic magnetic tissue properties using MRI signal phase: an approach to in vivo brain iron metabolism? Neuroimage. 2011;54(4):2789–2807. pmid:21040794
  14. 14. Li L, Leigh JS. High-precision mapping of the magnetic field utilizing the harmonic function mean value property. Journal of magnetic resonance. 2001;148(2):442–448. pmid:11237651
  15. 15. Haacke EM, Liu S, Buch S, Zheng W, Wu D, Ye Y. Quantitative susceptibility mapping: current status and future directions. Magnetic resonance imaging. 2015;33(1):1–25. pmid:25267705
  16. 16. Borrelli P, Palma G, Comerci M, Alfano B. Unbiased noise estimation and denoising in parallel magnetic resonance imaging. In: Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on. IEEE; 2014. p. 1230–1234.
  17. 17. Borrelli P, Palma G, Tedeschi E, Cocozza S, Comerci M, Alfano B, et al. Improving Signal-to-Noise Ratio in Susceptibility Weighted Imaging: A Novel Multicomponent Non-Local Approach. PLoS ONE. 2015 06;10(6):e0126835. pmid:26030293
  18. 18. Palma G, Piccialli F, Comerci M, De Michele P, Borrelli P, Cuomo S, et al. 3D Non-Local Means denoising via multi-GPU. Proceedings of the 2013 Federated Conference on Computer Science and Information Systems. 2013;p. 495–498.
  19. 19. Deoni SC. Quantitative relaxometry of the brain. Topics in magnetic resonance imaging: TMRI. 2010;21(2):101. pmid:21613875
  20. 20. Warntjes JB, Engström M, Tisell A, Lundberg P. Brain Characterization Using Normalized Quantitative Magnetic Resonance Imaging. PloS one. 2013;8(8):e70864. pmid:23940653
  21. 21. Bonnier G, Roche A, Romascano D, Simioni S, Meskaldji D, Rotzinger D, et al. Advanced MRI unravels the nature of tissue alterations in early multiple sclerosis. Annals of clinical and translational neurology. 2014;1(6):423–432. pmid:25356412
  22. 22. Stanisz GJ, Odrobina EE, Pun J, Escaravage M, Graham SJ, Bronskill MJ, et al. T1, T2 relaxation and magnetization transfer in tissue at 3T. Magnetic Resonance in Medicine. 2005;54(3):507–512. pmid:16086319
  23. 23. Neeb H, Schenk J, Weber B. Multicentre absolute myelin water content mapping: Development of a whole brain atlas and application to low-grade multiple sclerosis. NeuroImage: Clinical. 2012;1(1):121–130.
  24. 24. Neema M, Goldberg-Zimring D, Guss ZD, Healy BC, Guttmann CR, Houtchens MK, et al. 3 T MRI relaxometry detects T2 prolongation in the cerebral normal-appearing white matter in multiple sclerosis. Neuroimage. 2009;46(3):633–641. pmid:19281850
  25. 25. Hänicke W, Vogel HU. An analytical solution for the SSFP signal in MRI. Magnetic resonance in medicine. 2003;49(4):771–775. pmid:12652550
  26. 26. Soscia E, Palma G, Sirignano C, Iodice D, Innocenti S, Greco D, et al. Relaxometric Maps: Sequence Development And Clinical Impact. Initial Observations. European Congress of Radiology. 2011.
  27. 27. Sirignano C, Soscia E, Iodice D, Innocenti S, Palma G, Greco D, et al. T2 & T2 Maps: Sequence Development And Clinical Impact On Joint Study. European Congress of Radiology. 2012.
  28. 28. Heule R, Bär P, Mirkes C, Scheffler K, Trattnig S, Bieri O. Triple-echo steady-state T2 relaxometry of the human brain at high to ultra-high fields. NMR in Biomedicine. 2014;27(9):1037–1045. pmid:24986791
  29. 29. Deoni SC, Peters TM, Rutt BK. High-resolution T1 and T2 mapping of the brain in a clinically acceptable time with DESPOT1 and DESPOT2. Magnetic resonance in medicine. 2005;53(1):237–241. pmid:15690526
  30. 30. Deoni SC, Rutt BK, Arun T, Pierpaoli C, Jones DK. Gleaning multicomponent T1 and T2 information from steady-state imaging data. Magnetic Resonance in Medicine. 2008;60(6):1372–1387. pmid:19025904
  31. 31. Deoni SC. Transverse relaxation time (T2) mapping in the brain with off-resonance correction using phase-cycled steady-state free precession imaging. Journal of Magnetic Resonance Imaging. 2009;30(2):411–417. pmid:19629970
  32. 32. Zur Y, Wood M, Neuringer L. Motion-insensitive, steady-state free precession imaging. Magnetic resonance in medicine. 1990;16(3):444–459. pmid:2077335
  33. 33. Chavhan GB, Babyn PS, Jankharia BG, Cheng HLM, Shroff MM. Steady-State MR Imaging Sequences: Physics, Classification, and Clinical Applications 1. Radiographics. 2008;28(4):1147–1160. pmid:18635634
  34. 34. Wright P, Mougin O, Totman J, Peters A, Brookes M, Coxon R, et al. Water proton T 1 measurements in brain tissue at 7, 3, and 1.5 T using IR-EPI, IR-TSE, and MPRAGE: results and optimization. Magnetic Resonance Materials in Physics, Biology and Medicine. 2008;21(1–2):121–130.
  35. 35. West J, Blystad I, Engström M, Warntjes J, Lundberg P. Application of quantitative MRI for brain tissue segmentation at 1.5 T and 3.0 T field strengths. PloS one. 2013;8(9):e74795. pmid:24066153