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

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 (R 1), pure transverse (R 2) and free induction decay (R2*) 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.


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 (R 1 ), transverse relaxation rate (R 2 ) and proton density (PD) of brain tissues [2][3][4][5][6][7]. Other relaxation parameters whose MR signals can be usefully exploited to discriminate the microstructure of brain tissues are the free induction decay (FID) rate (R Ã 2 ) [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) B 1 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 R 1 , R 2 , PD, R Ã 2 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 R 1 , R 2 , PD and R Ã 2 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 R 1 , R 2 , PD, R Ã 2 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.

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 where K is the complex coil sensitivity, M 0 is the equilibrium magnetization, θ is the flip angle, E 1 exp(−T R Á R 1 ) (T R being the repetition time), T E 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 where E 2 exp(−T R Á R 2 ), φ = γ n ΔBT R + Δφ 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 φ $ ¼ g n DBT R .

PD-map
Once both R 1 and R Ã 2 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 jKjM 0 from (Eq 19), and averaging the results over the acquired echoes:

χ-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 B 0 is given by For the linear magnetic materials, (Eq 26) is completed by the magnetostatics constitutive equation: which reduces to for diamagnetic and paramagnetic substances whose susceptibility jχj ( 1. As in the MRI context B 0 ' B 0ẑ , theẑ component of the field variation due to susceptibility inhomogeneities is therefore given by where As we acquire a spoiled GRE sequence, the phase associated to susceptibility (see (Eq 1)) is therefore given by In order to extract ϕ χ (r) from the original wrapped phase ψ(r) (which also depends on background magnetic field inhomogeneities, the imaginary coil sensitivity I(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 (r 2 (ϕ−ϕ χ ) = 0) [14]: The inversion of (Eq 32) requires a regularization of the inverse of g(k) FT[G(r)] (see [8,15]), whence we finally obtain

R 2 -map
Once the R 1 -and the PD-maps are known, (Eq 6) makes it possible to derive the R 2 -map. From the following combination of Eqs (3), (5) and (6) if the acquired bSSFP scans are set with T E = T R /2, it follows that From the standard theory of the quartic equations it follows that the physically meaningful solution of (Eq 37) is given by where q a 1 a 4 ; ð43Þ The transverse relaxation rate is then found from the expression 3 Materials and Methods

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. 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 B 1 sensitivity profile, we describe the method developed to account for non-ideality of the B AE 1 profiles within the Field of View (FOV). Finally, we detail the way devised to estimate the overall reproducibility of the results.

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 mm 3 ).
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.
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.

Image restoration
In order to control the noise propagation in the relaxometry maps derived by inversion of the Bloch equations, the FC-FLASH and S 0 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-   (ρ and B), that act as follows: ; ð52Þ is the similarity radius, B * 1 is an adimensional constant to be manually tuned that determines the filter strength and s m ðxÞ is the standard deviation of noise of the m th image component atx 2 O (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.

-inhomogeneity correction
The derivation of R 1 , PD and R 2 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 B þ 1 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 B À 1 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 B þ 1 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 r þ 0 , modulated by a slab profile term: (z is the offset from the slab center in the slab select direction). As the Hessian tensor H + is symmetric, r þ 0 is a vector and k, ξ 2 R, a total of 11 free parameters need to be estimated in (Eq  53).
Similarly, the B À 1 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): Finally, the parameters are estimated by minimizing, with the Nelder-Mead method, the entropy of the 3D joint histogram related to R 1 , PD and R 2 maps in low gradient regions of the maps.
The effectiveness of the obtained B 1 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 R 1 and R 2 distributions reconstructed with and without the B AE 1 -inhomogeneity correction.

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 2 4 (spoiled GRE) × 3 (S 0 ) 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.
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 R 1maps (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.

Effectiveness of B AE 1 -inhomogeneity correction
The test of the effectiveness of the proposed B AE 1 -inhomogeneity correction on R 1 and R 2 is reported in Table 2.  Thanks to the acquired voxel size, the reconstructed datasets allow for reasonably good multiplanar reconstructions.

Map assessment
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. In both HC and MS patient the R 1 -map was rated as the most representative of the brain anatomy, with excellent differentiation between grey (GM) and white (WM) matter. The R Ã 2 -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 R 1 -and R 2 -maps, even in critical areas, such as the periventricular WM (Fig 6).

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

Discussion
We have described a 3D acquisition protocol that yields multiple quantitative parametric maps (R 1 , R 2 , R Ã 2 , 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,[20][21][22][23][24], 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 R 1 and R 2 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 R 1 -and R Ã 2 -maps of the brain [6,23], but not the quantitation of pure transverse relaxation rate (R 2 ) or susceptibility (χ). These may provide useful information when R 1 and R Ã 2 fail to provide adequate information or contrast to reveal the pathophysiological conditions of specific areas (e.g. Table 2. Percentual FWHM of the R 1 and R 2 distribution obtained without (pre) and with (post) B AE 1 -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.     in presence of high susceptibility gradients near petrous bones, nasal sinuses, etc. that cause an R Ã 2 enhancement uncorrelated to any tissue alteration).
Recently, some 3D approaches have been proposed to derive R 1 -and R 2 -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 T 2 , thus making the brain study unfeasible, due to CSF pulsation, unless switching to 2D [28].
In this respect, the family of DESPO methods [2,[29][30][31], while neglecting R Ã 2 -and χ-maps, do provide a full 3D high SNR R 1 -and R 2 -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 B 0 inhomogeneities sensitivity of bSSFP. However, some drawbacks of DESPOT2 may raise crucial issues even if R Ã 2 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 R 2 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 T R 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 B 0 inhomogeneity by means of a fully analytical solution (thus speeding up the computation step). As for the B AE 1 inhomogeneity dependence, one can completely remove it using a measured B 1 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.
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 B AE 1 -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.