Strategies for assessing diffusion anisotropy on the basis of magnetic resonance images: comparison of systematic errors.

Diffusion weighted imaging uses the signal loss associated with the random thermal motion of water molecules in the presence of magnetic field gradients to derive a number of parameters that reflect the translational mobility of the water molecules in tissues. With a suitable experimental set-up, it is possible to calculate all the elements of the local diffusion tensor (DT) and derived parameters describing the behavior of the water molecules in each voxel. One of the emerging applications of the information obtained is an interpretation of the diffusion anisotropy in terms of the architecture of the underlying tissue. These interpretations can only be made provided the experimental data which are sufficiently accurate. However, the DT results are susceptible to two systematic error sources: On one hand, the presence of signal noise can lead to artificial divergence of the diffusivities. In contrast, the use of a simplified model for the interaction of the protons with the diffusion weighting and imaging field gradients (b matrix calculation), common in the clinical setting, also leads to deviation in the derived diffusion characteristics. In this paper, we study the importance of these two sources of error on the basis of experimental data obtained on a clinical magnetic resonance imaging system for an isotropic phantom using a state of the art single-shot echo planar imaging sequence. Our results show that optimal diffusion imaging require combining a correct calculation of the b-matrix and a sufficiently large signal to noise ratio.

INTRODUCTION

Diffusion-weighted (DW) magnetic resonance imaging (MRI) has become a valuable technique for the detection and follow-up of cerebral and cardiac ischemia.[12345] In addition, many of studies have been performed on the anisotropic character of diffusion in brain,[67] muscle,[8910] bone.[1112131415] These studies were very informative about the dependence of water diffusion on the architecture of the tissues.[71216171819] Diffusion in an isotropic medium strongly influences contrast in DW images acquired by introducing diffusion weighting gradients before and after the 180° pulse of a spin-echo MRI sequence. Combining images corresponding to a fixed diffusion gradient direction, but variable strengths, an apparent diffusion coefficient map or ADC map can be calculated.

“Apparent” because the presence of cellular structures influences the random motion of the water molecules and leads to a diffusion coefficient that is smaller than that of pure water.

For an isotropic medium, the measured ADC is expected to remain the same if we change the direction of the diffusion weighting gradients.[2021] In order to characterize the diffusion behavior in an anisotropic medium, it is necessary to determine six independent elements of the symmetric 3 × 3-diffusion tensor (DT) D.[1721] This means that the experiment described earlier needs to be repeated for six independent choices for the direction of the diffusion weighting gradients. The derived DT components are rotationally variants, i.e., they are dependent on the directions of the applied diffusion gradients (i.e., the orientation in the laboratory frame[1721]). From the DT components; however, we are able to derive rotationally invariant (RI) tissue characteristics.[171821] These quantities include the three sorted eigenvalues and trace of the 3 × 3 symmetric DT D, the anisotropy ratio index and the volume ratio (VR) index.[18]

The accuracy of both the rotationally variant and rotationally invariant quantities is limited by the propagation of signal noise in the course of the calculations of the required quantities starting from the initial DW image data, i.e., the DW magnetic resonance images.[18] Another source of inaccuracy is the method used for expressing the impact of the sensitizing and imaging gradients on the signal measured.[12161718]

The goal of this paper is to compare the relative importance of these two factors in the determination of the DT D and the derived rotationally invariant quantities. Two schemes for the calculation of the b-matrix elements for a single shot echo-planar MRI sequence have been considered:

Only the diffusion gradients with their ramp-up and ramp-down times are taken into account in the calculation of the b-matrix elements, which we denote b0. This is the approach that is most often used in the clinical applications[45678910131415]

All diffusion and imaging gradients are taken into account in MRI sequence, including the sinusoidal readout gradients, the ramp-up and ramp-down times of all the gradients and the existence of cross terms. We denote the corresponding b-matrix b1.

We have calculated b-matrix elements using both methods and compared the diffusion results derived from data with varying signal to noise ratio (SNR) for each case. The experiments were carried out on an isotropic phantom, which should in principle lead to a DT with a single independent element.

DT imaging is considered as a relevant means for assessing the healthy and diseased tissues such as the alterations occurring in the pediatric brain myelination process, demyelization disorder, neoplasm involving white matter tract and stroke. It was found in studies involving patients with focal ischemia that the anisotropy ratio for healthy and diseased tissue was not significantly different, or slightly elevated during the 1st day post-occlusion, then decreased during the latter chronic phase, when diffusion is less impeded through cell membrane degeneration. However, a statistically significant increase was found in several rotationally variant indices of diffusion anisotropy measured in regions of small cortical and lacuna strokes in eight patients presenting 1-5 days post-occlusion.[222324]

THEORETICAL BACKGROUND

Background of the DT Theory

The diffusion phenomenon is a physical process involving the random motion of water molecules in a given medium while being driven by kinetic energy theory. The statistical physics described the phenomenon by the first law of Fick:[671617181920212223]

J = −D∇c      (1)

This law supports that a difference in concentration of a given moving molecules is creating a net flux from regions with a high concentration into the direction of regions with low concentration. Where J is the full flux, ∇c is the gradient of concentration and D is the diffusion coefficient. This, the diffusion coefficient D dependents on factors including dimensions of the diffusing molecule and the microstructure characteristics constituting the diffusion medium and the temperature of the medium. The diffusion coefficient is given by the Einstein equation:

Where k is the Boltzmann constant, T and η are the temperature measured in Kelvin and viscosity of the diffusion medium respectively, r is the radius of the diffusing molecule.

The net flux vanishes when concentration gradient become null, however molecules do continue moving to support the thermodynamic equilibrium.

This phenomenon is known as Brownian motion and was not been described as a stochastic motion mode. Molecules movement could be represented by a Gaussian distribution displacement model depending on D.

The probability for a given molecule to travel a distance x in the time t is defined by:

Indeed, the mean square displacement during the time Δt is given by:

〈x2〉 = 2DΔt      (4)

The x distance was travelled by the diffusing molecule during Δt. When there are not any boundaries for molecules motion, diffusion is described as free and follows the random Brownian motion model. This diffusion model is called isotropic.

For an isotropic medium, the DT becomes:

D = D0.I      (5)

Where Do is the scalar diffusion coefficient and I is the unity matrix given by:

A normalized DT D/D0 = Dnorm should reflect the specific constraint for an isotropic medium, namely that the three diagonal elements must be equal to each other (D/D0 = 1) and the off-diagonal elements should cancel to zero. Significant deviations from this property are to be ascribed to systematic differences in sensitivity of the various MRI sequences to diffusion phenomena. Indeed, connecting the observed anisotropy for isotropic media to errors in calibration, misalignment and cross-talk between the imaging and diffusion gradients has to be assessed.[671617181920212223]

Nevertheless, when molecules movement is hindered by obstacles, the diffusion will demonstrate preferential direction that is not associating barrier and the diffusion is said to be anisotropic. This last case is particularly important when studying the diffusion of water within human biological tissue, since composed by organized structures that are guiding the water movement which is not diffusing in free mode. Often these structures are not visible with conventional MRI techniques since their dimensions magnitudes are smaller than the resolution of a magnetic resonance image.

These aspects are important when interested by brain tissue magnetic resonance images. Thus, various brain tissues are highly characterized by their diffusion properties including white matter, gray matter and cerebrospinal fluid.[671617181920212223]

Expression for the b Matrix and DT

From the signal intensity in the DW magnetic resonance images, we can compute the DT elements:[2526]

The b matrix can be calculated as follows:

where

and γ is the gyromagnetic ratio, G (t) = [G (t), G, G] the column vector representing the gradient pulses in the diffusion imaging sequence, TE the echo time and ξ (t) the Heaviside function with:

The DT D is expressed in the laboratory frame x, y, z as follows:

In the same frame of reference expression[5] can be rewritten under the form:[25]

I (b) is the signal intensity for particular diffusion gradient sensitization and I (0) is the signal intensity measured in the absence of the diffusion gradient sensitization.

The b0 matrix elements are calculated using the following approximated expression:[25]

δ is the duration of the applied trapezoidal diffusion gradients, Δ the delay between the start times of the symmetric diffusion gradients before and after the 180° pulses, ε represents the fixed ramp-up and ramp-down time of the trapezoidal diffusion gradients.

The elements of the b matrix b1, are computed from Eqs. (7) and (8), with a G (t) that accurately reflects all gradients switched on at any time during the run of the DT MRI sequence.

Eigenvalues Sorted and the Trace

In the previous section DT, D, was expressed in the laboratory frame of reference. After diagonalization of the DT and sorting of the eigenvalues, we obtain a diagonal diffusivity tensor, λ, expressed in the following:

For an isotropic medium Eq. (11) becomes λ=λ0I, where I is the unit tensor and

where we have introduced the trace of the DT:

Anisotropy Indices

In order to express quantitatively the degree of anisotropy in the tissue studied, a set of indices derived from the sorted eigenvalues can be used:[16171825]

The anisotropy ratio, AR, compares the largest diffusivity to the mean of the two others:

The anisotropy ratio index is equal to 1 for a fully isotropic medium.

The VR, is the volume of an ellipsoid the main axes of which are the three sorted eigenvalues of the DT D, divided by the volume of a sphere with radius equal to the mean diffusivity.[18]

VR varies from 0 for strong anisotropy to 1 for complete isotropy.

MATERIALS AND METHODS

Numerical Calculation of the b Matrix Elements

Mathcad Software (Math Soft International, UK) running on a personal computer, was used to perform the calculation of the b-matrix elements for the two schemes considered. b0 was computed from the approximated Eq. (10), whereas b1 was compute shot DW echo planar imaging sequence (EPI) following a Stejskd from the integral expression of the b matrix, using the numerical integration for Eq. (7).[27]

Data Acquisition of the DW MRI

The DW MRI were obtained using a 1.5 Tesla Magnetom Vision whole body imager (Siemens, Erlangen, Germany). The system was equipped with an actively shielded gradient system with maximum gradient strength 25 mT/m. An isotropic phantom composed by water doped with NiSO4 × 6H2O. All experiments were carried out at a fixed temperature of ~ 22°C. The isotropic phantom considered was positioned in a manufacturer supplied head coil used for routine head clinical examinations. The DW MRI was performed using a single al-Tanner approach.[2829] The pulse sequence is illustrate in Figure 1.

Figure 1

Schematic diagram of the diffusion-weighted echo-planar imaging sequence

The MRI parameters were, echo time TE = 123 ms, field of view = 240 mm, square matrix size = 128 × 128, slice thickness = 5 mm. The duration of each diffusion gradient δ = 26 ms, the time between the starts of each of the two symmetrical consecutive diffusion gradients Δ = 59.7 ms. The rise times of the diffusion gradients were fixed to ε = 700 ms. The trapezoidal diffusion gradients strengths were 0, 11 and 22 mT/m. The single slice orientation was transverse, i.e. perpendicular to the external magnetic field. The basic measurement lasted 38 s, yielding a series of 13 images. Each series consisted of one image without diffusion weighting and one weakly and one strongly DW image for each of six non-collinear magnetic field gradient directions: {(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), . In those experiments we have increased the SNR of the DW-EPI images by averaging over multiple acquisitions, which is commonly done in the clinical practice. This consisted of repeating the entire imaging experiment number of acquisitions times and averaging the signal over these number of acquisitions measurements. As it is difficult to attach a single meaningful SNR to the 13 images comprising a DT measurement, we have preferred to report the results as function of the number of acquisitions, which is the same for all these images.

Data Post-processing

The echo-planar diffusion magnetic resonance images data was transferred from the scanner to a personal computer running Linux for off-line post-processing, using C-programs developed in our group for this purpose (30). From the 13 DW magnetic resonance images measured for each slice considered, we derived the six independent elements of the DT.[30] After diagonalization of the DT, a sorting routine produced the sorted eigenvalues. Finally, we calculated the trace of the DT and two anisotropy indices: The anisotropy ratio AR and the VR. All calculations were carried out on a pixel by pixel basis. We evaluated the different parameters in a fixed Region of Interest (ROI) on the phantom images.

RESULTS

in the following, all quantities derived from an analysis using b matrix models b0 and b1 will carry the superscripts 0 and 1 respectively. All errors correspond to one standard deviation.

Numerical Calculations of the b Matrix Elements

As you can be seen in Table 1 the b matrix elements calculated for three values of the diffusion gradients for six independent choices of the gradient direction of the MRI sequence using both algorithms are compared. The systematic deviation between the elements of b0 and those of b1 varies from less than 1-7% and depends strongly on the gradient directions.

Table 1

Percentage of deviation between both sets of b-matrix elements, b0 and b1. The elements in the table marked by (0) correspond to a value that is <1%

The Diagonal Elements of the DT

As shown in Figure 2a and b, the number of acquisitions dependence of the diagonal elements of the DT for an isotropic phantom in the laboratory frame derived using both b matrix models is shown and compared to λ0.

Figure 2

Diagonal elements of the diffusion tensor (mean ± standard deviation) × 10–3 mm2/s

With low number of acquisitions, the values of Dii0 and Dii1 (i: x, y and z) deviate significantly from the reference value λ0 = 2.17 × 10-3 mm2/s, defined as the asymptote of the arithmetic mean of the λk1 in function of number of acquisitions (k = 1, 2 and 3). Although the elements Dii1 are found to converge to λ0 with increasing number of acquisitions, this is not the case for Dii0, even at the highest number of acquisitions.

The Sorted Eigenvalues

Figure 3a and b describe corresponding results for the sorted eigenvalues, as it was the case for the diagonal elements of the DT. At low number of acquisitions, the values of λk0 and λk1 (k = 1, 2 and 3) deviate significantly from the reference value λ0. The deviations are found to be similar for λk0 and λk1. Again, while the elements λk1 are found to converge with increasing number of acquisitions to the known common value of reference λ0, this is not the case for λk0 even at the highest number of acquisitions. The deviation of λk0 from λ0 at highest number of acquisitions is stronger in comparison to Dii0.

Figure 3

Sorted eigenvalues (mean ± standard deviation) × 10–3 mm2/s

The Trace and the Anisotropy Ratio Index

As shown in Figure 4a, at high number of acquisitions, Trace1 has the same values as λ0, as expected, whereas it deviates from λ0 by less than 2% at the lowest number of acquisitions. Trace0 deviates from λ0 by an extra 1% with respect to Trace1 at all values of number of acquisitions.

Figure 4

Mean ± standard deviation of the trace, anisotropy ratio and volume ratio as function of number of acquisitions in an isotropic phantom

Figure 4b shows the AR dependence on number of acquisitions. At low number of acquisitions the values of AR0 and AR1 deviate from the isotropic value by 10% and 5% respectively, while their deviations are 5% and 2% at higher number of acquisitions. However, AR0 and AR1 deviate from each other by 5% at the lowest number of acquisitions while this deviation remains stable (3%) at number of acquisitions >20.

Figure 4c shows that the variation of VR0 and VR1 is within 1% for almost all values of number of acquisitions. This variation is practically negligible.

DT Maps and Indices

Figure 5 demonstrates maps of the diagonal element of the DT (Dxx, Dyy, Dzz) calculated using 20 averages in the first row and one acquisition in the second row (a). In (b) the eigenvalues of the DT (λ1, λ2, λ3) are derived from DT calculated using 20 averages in the first row and one acquisition in the second row. All maps were calculated from DW images of the same phantom. Indeed the maps acquired using higher number of acquisitions demonstrates converging contrast in this isotropic phantom. This reflects reduced systematic errors as function of number of acquisitions. Similarly, Figure 6 demonstrates Trace of the DT, Anisotropy ration (AR) and VR indices derived from DT calculated using 20 averages in the first row and one acquisition in the second row. All maps were calculated from DW images of the same phantom. Indeed the maps of indices acquired using higher number of acquisitions demonstrate converging contrast in this isotropic phantom. This reflects reduced systematic errors as function of number of acquisitions.

Figure 5

Maps of the diagonal element of the diffusion tensor (Dxx, Dyy, Dzz) calculated using 20 averages in the first row and one acquisition in the second row (a). In (b) the eigenvalues of the diffusion tensor (λ1, λ2, λ3) are derived from diffusion tensor calculated using 20 averages in the first row and one acquisition in the second row. All maps were calculated from diffusion weighted images of the same phantom

Figure 6

Trace of the diffusion tensor, anisotropy ration index and volume ratio index derived from diffusion tensor calculated using 20 averages in the first row and one acquisition in the second row. All maps were calculated from diffusion weighted images of the same phantom

DISCUSSIONS

When producing normative data for characterizing the anisotropy and fiber orientation in specific anatomic regions, accurate diffusion anisotropy maps are essential in achieving significant results.[3119] Any error in the calculated diffusion quantities reflecting the anisotropy could make an isotropic medium appear to be anisotropic and vice versa. Our data reflects the impact of the two main types of systematic error on the DT elements and derived indices.

Firstly, systematic errors those are due to neglecting the imaging gradients of the pulse sequence when calculating the b-matrix elements. Eqs (1) and (4) show that for a given number of acquisitions significant deviations between the diagonal elements of the DT for an isotropic phantom could be a consequence of the b-matrix being incorrectly calculated. The error in the calculations of the different diagonal and off-diagonal elements of the b-matrix is expected to vary because of the varying role played by the diffusion and imaging gradients in each of the various laboratory frame directions for diffusion weighting. Furthermore, the interaction between imaging and diffusion gradients is very complicated, especially when measuring diffusion along directions that do not coincide with the principal axes of the laboratory frame (e.g. xy, xz and yz). On the basis of computer simulations, Mattiello et al.,[26] have been reported that a given percentage of error in b-matrix elements produces the same percentage of error but of opposite sign in corresponding elements of the statistically estimated DT elements.

Secondly, systematic errors could also be caused by other factors such as signal noise, eddy currents, susceptibility effects and the calculation process that starts from the initial DW-imaging data and finishes at the desired diffusion parameter. In practice, diagonalizing the DT and sorting the eigenvalues is susceptible to accumulation of errors. Our data reflects error in DT elements due to the propagation of noise present in the acquired DW images.[3233] Using Monte Carlo simulations, Pierpaoli and Basser[18] and Bastin et al.,[34] have predicted that the effect of noise on the diffusion quantities (diagonal elements of the DT, eigenvalues sorted and the reported anisotropy indices) could be drastically reduced when measuring at high number of acquisitions. Our experimental results corroborate this finding.

The importance of the deviations of the diffusion parameters from λ0 is not easy to assess, as the random errors on the mean values of the parameters is smaller than the standard deviations plotted in Figures 2–4 by a factor (Number of pixels in the ROI)-1/2.

A comparison with the uncertainty in diffusion quantities due to the temperature fluctuation with 1°C shows that the average values of the ADCs and diffusivities obtained using the full b-matrix deviate by less than this uncertainty for number of acquisitions ≥5. When using the approximated b-matrix, this observation remains valid for the ADCs while the diffusivities exhibit larger deviation, even at the largest number of acquisitions.

The difference between the use of b0 and b1 for calculating the DT elements and derived quantities is most obvious at high number of acquisitions where only the effect of ignoring imaging gradients remains and the impact of noise is minimized. Even at the highest number of acquisitions considered, deviations up to 5% from λ0 are caused by using b0 instead of b1.

When using b1, the only systematic deviations that remain are due to noise [Figures 2b, 3b and 4]. In this case, at highest number of acquisitions, the Dii and λk deviate from λ0 by at most 2%. The effects of all systematic errors are less visible on Trace and VR than for AR. This means that the errors are not transmitted in the same way to these parameters. Figure 4 shows also that an estimate of the diffusion anisotropy can be made within the same margin of deviation as for Dii and λk.

Finally, DTI results obtained by combining number of acquisitions ≥5 and the use of the full b-matrix (b1) lead to the same level of deviation from λ0 as temperature fluctuations by 1°C.

CONCLUSIONS

To perform a measurement of diffusion anisotropy in non-homogeneous tissues like brain white matter, in which the fibers are oriented in complex way, is a challenging task. It requires performing the DW imaging with very high accuracy.

Due to image noise and b-matrix approximations systematic errors can produce erroneous conclusions of anisotropy. With correct calculation of b-matrix and sufficiently large number of acquisitions spurious anisotropy can be avoided.

BIOGRAPHY

Saïd Boujraf MedSc, PhD is graduated Medical Researcher from ‘Free University of Brussels’ (VUB), Belgium He is currently a Professor of Biophysics and Clinical MRI Methods, Department of Biophysics and Clinical MRI Methods, Faculty of Medicine and Pharmacy, University of Fez, Morocco Director of the Clinical Neuroscience Laboratory, Faculty of Medicine and Pharmacy, University of Fez, Morocco.

E-mail: sboujraf@gmail.com