A new algebraic method for quantitative proton density mapping using multi-channel coil data

Abstract

A difficult problem in quantitative MRI is the accurate determination of the proton density, which is an important quantity in measuring brain tissue organization. Recent progress in estimating proton density in vivo has been based on using the inverse linear relationship between the longitudinal relaxation rate T₁ and proton density. In this study, the same type of relationship is being used, however, in a more general framework by constructing 3D basis functions to model the receiver bias field. The novelty of this method is that the basis functions developed are suitable to cover an entire range of inverse linearities between T₁ and proton density. The method is applied by parcellating the human brain into small cubes with size 30mm x 30mm x 30mm. In each cube the optimal set of basis functions is determined to model the receiver coil sensitivities using multi-channel (32 element) coil data. For validation, we use arbitrary data from a numerical phantom where the data satisfy the conventional MR signal equations. Using added noise of different magnitude and realizations, we show that the proton densities obtained have a bias close to zero and also low noise sensitivity. The obtained root-mean-square-error rate is less than 0.2% for the estimated proton density in a realistic 3D simulation. As an application, the method is used in a small cohort of MS patients, and proton density values for specific brain structures are determined.

Keywords: Bias field; Proton density; Quantitative MRI; Receiver coil sensitivity; T(1); Transmission coil sensitivity.

Fig.1

Flow chart showing the proposed analysis steps. The novelty of our approach is how the transmission coil sensitivity is estimated, the spatial basis functions are created and used to model the receiver coil sensitivities, the spatial expansion coefficients of g are calculated, and the scaling factors of the proton density are determined simultaneously for all cubes in one step.

Fig.2

Top: Inverse linear relationship between proton density ρ and T₁ for a single parameter {A, B} (blue line) and for a family of parameters from which the novel basis functions are derived (green lines covering an entire area). Bottom: First nine spatial basis functions derived after PCA orthonormalization for a cube (30mm × 30mm × 30mm) in white matter. A 3D view of each basis function is shown as well as five 2D-slices of each basis function at indicated distances from the right front face. The first basis function represents a constant. All the other basis functions are orthonormal over the cube. For the nearest coil (coil 1), dominant contributions are from functions 1, 2, 3 and 4 whereas for the other more distant coils, main contributions are from functions 1, 2, 3. Functions 5 to 9 also contribute but have a weight that is about one magnitude smaller.

Fig.3

2D noiseless toy example involving 3 coils on a 64×64 pixel grid and error obtained in recovering the images from the signal equation (Eq.(3)) for three different methods. The coil sensitivities contain a mixture of polynomials up to 2^nd order using 6 functions {const, x, y, xy, x², y²} and an additional 2D Gaussian function (amplitude 30, centered at x = y = 32, standard deviation σ = 40). The proton density is a uniform random image with values 0.5 < ρ(x, y) < 1. A corresponding T₁ image satisfying the inverse linear relationship $\frac{1}{\rho }=0.879+\frac{503\mathit{\text{ms}}}{T_1}$ was also generated. The signal amplitudes for all coils were computed according to Eq.(3). A. Simulated coil sensitivities g^(1,2,3)(x,y) and proton density ρ(x,y). B. Error Δg^(1,2,3) = g^(1,2,3)(ground truth) − g^(1,2,3)(estimated) and Δρ = ρ(ground truth) − ρ(estimated) obtained using a polynomial basis up to 2^nd order {const, x, y, xy, x², y²} in the modeling of the coil sensitivities. C. Error Δg^(1,2,3) = g^(1,2,3)(ground truth) – g^(1,2,3)(estimated) and Δρ = ρ(ground truth) − ρ(estimated) obtained using a polynomial basis up to 3^rd order with 8 functions {const, x, y, xy, x², y², x³, y³}. D. Error Δg^(1,2,3) = g^(1,2,3)(ground truth) − g^(1,2,3)(estimated) and Δρ = ρ(ground truth) − ρ(estimated) obtained using the proposed optimized basis set similar to the ones in Fig.2 except for the 2D case.

Fig.4

Images of simulated (left) and estimated (right) MR quantities. From top to bottom: A. longitudinal relaxation rate T₁ (units are ms), B. observed transmission coil sensitivity m, C. receiver coil sensitivity g⁽ⁱ⁾ of one arbitrary chosen coil element, D. proton density ρ.

Fig.5

Estimation accuracy of simulation: A. longitudinal relaxation rate T₁ (in ms), B. observed transmission coil sensitivity m, C. receiver coil sensitivity g⁽ⁱ⁾ averaged over all coil elements, D. proton density ρ.

Fig.6

Root-mean-square error (RMSE) in % of the estimation accuracy for longitudinal relaxation rate T₁, transmission coil sensitivity m, signal amplitude M₀, receiver coil sensitivity g⁽ⁱ⁾ (averaged over all 32 coil elements) and proton density ρ, using simulated data with different noise fractions. A noise fraction of 0 indicates no noise added and a noise fraction of 1 indicates the same noise level as estimated from real MRI data. The overall error for all estimated quantities is less than 1 %. The small bias at zero noise fraction is due to partial volume effects introduced by the coregistration step of all generated images. Note that the noise sensitivity of ρ is by a factor of about 0.4 smaller than the noise sensitivity of T₁.

Fig.7

Comparison of proton density results using a numerical brain phantom. A. Ground truth images of proton density. B. Results obtained using the method by Mezer et al. (2016) rescaled so that the mean proton density is the same as in the ground truth images. C. Difference map ρ(B) – ρ(A) in percent. D. Percent difference map between the images obtained by the proposed method (this research) and the ground truth in A. Note the difference in scale between C and D.

Fig.8

Proton density results using public data for a normal subject from Mezer et al. (2016). A. Scaled proton density map results according to Mezer et al. (2016) B. Proton density map obtained using the proposed method (this study). The images in A were scaled so that the mean proton density is the same as in B. C. Difference map ρ(A) − ρ(B) in percent. D. Effect of CSF normalization using the final proton density images with no scaling involved for the images obtained by the method of Mezer et al. (2016) and the proposed method. Here, none of the proton density images were mean adjusted. The vertical axis shows the empirical cumulative density function (CDF) for the proton density in the lateral ventricles segmented by FreeSurfer. Note that 20% of voxels have a ρ > 1 using the method by Mezer et al. (blue curve) whereas for the proposed method max(ρ) = 1 (red curve).

Fig. 9

Estimated images of a representative subject (subject 1): A. longitudinal relaxation time T₁, B. observed transmission coil sensitivity m, C. kernel-density estimate of the distribution of the observed transmission coil sensitivity m, D. signal amplitude M₀, E. receiver coil sensitivity g⁽ⁱ⁾ of one arbitrary chosen coil element, F final proton density ρ after normalization.