Effects of diffusion weighting schemes on the reproducibility of DTI-derived fractional anisotropy, mean diffusivity, and principal eigenvector measurements at 1.5T

Data acquisition

A healthy 24-year-old male volunteer participated in this study. Local institutional review board approval and written informed consent were obtained prior to examination. All data were acquired using a 1.5T MR scanner (Intera, Philips Medical Systems, Best, The Netherlands) with body coil excitation and a six-channel phased array SENSE head coil for reception. Three scanning sessions were performed over 2 days with the subject repositioned between each session. In each scanning session, 15 DTI datasets were acquired, yielding a total of 45 DTI scans. Each DTI dataset was acquired with the following imaging protocol. A multi-slice, single-shot EPI (SENSE factor = 2.0), spin echo sequence (TR/TE= 3632/100 ms) was used to obtain 25 transverse slices parallel to the line connecting the anterior and posterior commissures with no slice gap and 2.5 mm nominal isotropic resolution (FOV= 240 × 240, data matrix= 96 × 96, zero-filled and reconstructed to 256 × 256). Diffusion weighting was applied along 30 distinct directions (Jones et al., 1999; Skare et al., 2000) with a gradient strength of G = 19.5 mT/m and a b-factor of 1000 s/mm². Five minimally weighted images ($b_0\text{s}$) were acquired and averaged on the scanner as part of each DTI dataset. The total scan time to acquire one DTI dataset (30 DW images and 5 $b_0$ images) was 2 min 18 s. The total time, including image reconstruction, to acquire 15 DTI datasets and one anatomical MRI in a scan session was approximately 45 min.

Motion correction and SNR calculations

DTI datasets were processed offline using Matlab (Mathworks, Natick, MA) routines running on a Sun Fire V880 server (Sun Microsystems, Mountain View, CA). All data were coregistered with a purpose-constructed method (Landman et al., 2006) utilizing FLIRT (FMRIBTs Linear Image Registration Tool, Oxford, UK) (Jenkinson et al., 2002) to remove rigid body motion with six degrees of freedom registration. The gradient tables, which specify the direction of the magnetic field gradient corresponding to each particular DWI, were corrected to compensate for the rotational component of the coregistration procedure. Diffusion tensors were calculated using a multivariate log-linear fitting method in which each DW image and its corresponding vector from the gradient table entered the diffusion tensor calculation as unique entries—i.e., no averaging of DW images was performed. These analysis methods were encapsulated and publicly released in CATNAP (Coregistration, Adjustment and Tensor-solving, a Nicely Automated Program, http://iacl.ece.jhu.edu/~bennett/catnap/, Johns Hopkins University School of Medicine, Baltimore, Maryland, USA).

The log-linear fit was performed with no restrictions on the sign of the eigenvalues. To prevent taking the log of zero, a voxel with an intensity of zero in a DW image was set to half of the smallest non-zero voxel intensity recorded. The eigenvector associated with the largest eigenvalue (${\text{λ}}_1$), i.e., the principal eigenvector (PEV), is taken to represent the predominant fiber orientation in each voxel. In the color-coded PEV maps, the dominant fiber orientation in a voxel was represented by red, green, and blue colors which were assigned to right–left, anterior–posterior, and superior–inferior orientations, respectively (Jones et al., 1997; Pierpaoli, 1997).

There are several difficulties in calculating an experimental SNR value that can be precisely compared with the values resulting from simulation studies. Whereas the noise in simulation studies is a modeled parameter, the SNR in in vivo imaging is spatially varying (due to coil sensitivity) and is tissue dependent (via the interaction of imaging sequence parameters and local tissue characteristics, e.g., T1, T2, PD). Furthermore, MR images can be subject to chemical shift artifacts as well as artifacts due to field inhomogeneity, eddy currents, and EPI-based geometric distortions.

To compare results with the literature and ensure that subsequent simulation studies correspond to approximately the same SNR range as the experimental data, we used a manual region of interest calculation. For clarity, we report SNR in terms of power decibels $\left(20{\log }_{10}\frac{\text{signal amplitude}}{{\sigma }_{\text{noise}}}\right)$ and contrast ratio $\text{(signal amplitude:}{\sigma }_{\text{noise}}\text{)}$. The SNR in the splenium of the corpus callosum was approximately 25 dB (17.8:1) in the averaged $b_0$ image, which would indicate an approximate SNR of 18 dB (7.9:1) in a single $b_0$ image.

Partitioning and grouping of diffusion-weighted data

The DW data in this study were acquired with a PE optimized 30 direction scheme (Jones et al., 1999; Skare et al., 2000), hereafter referred to as the Jones30. To investigate the effects of different DW schemes, the Jones30 scheme was partitioned into optimal PE-minimized subsets of 6, 10, and 15 DW directions with a Monte Carlo pair-wise relaxation method (Landman et al., 2006). These subsets are denoted as PE₆, PE₁₀, and PE₁₅, respectively. The objective in this optimization step was to identify the closest approximation of a minimum PE scheme possible given the limitation of using data acquired with the Jones30 DW scheme. The Appendix Table A1 tabulates which directions from the Jones30 scheme were assigned to each of the optimized PE schemes.

The selection of optimal subsets was accomplished by Monte Carlo restarting of a local gradient descent algorithm. The algorithm was initialized by selecting N directions at random from the Jones30 scheme. This subset of size N was then optimized through an exhaustive search that involved a member–non-member exchange to minimize the previously described PE function (Jones et al., 1999). Fig. 1A illustrates the DW directions obtained by optimization for the PE₆, PE₁₀, and PE₁₅ schemes as compared to the standard Jones30. The PE of the partitioned DW schemes (99.0, 303.5, 726.7) was at most 1% higher than the unconstrained PE optima for the equivalent number of directions (98.3, 301.8, 719.5).

It is important to note the differences between the specified DW direction and the realized diffusion weighting direction due to patient motion. We use the term “specified” to refer to the set of DW directions that are given to the MR scanner at the time of acquisition. These directions are relative to either the slice layout (e.g., measurement, phase, and slice orientations) which was the case for this study or the fixed magnetic coordinate frame (e.g., the bore), depending on the scanner settings. When a heterogeneous sample, e.g., a human brain, is placed in the scanner, the diffusion weighting is also relative to the position and orientation of the anatomy. It is possible that volunteers will move or rotate during a sequence, and therefore the relative orientation between the scanner coordinate system and the anatomy changes. Image coregistration corrects for this in that it forces voxels in one anatomical location of one scan to correspond to the same anatomical location in a subsequent scan. However, motion correction also changes the relationship between the anatomical and scanner coordinate systems, and thus the diffusion weighting. When rotating the anatomy to establish spatial correspondence, the DW direction in the gradient table must also be rotated so that the vector that describes the diffusion weighting is correctly interpreted relative to the anatomy. We refer to the corrected DW directions, which are slightly different than the “specified” ones, as the “realized” DW directions. In an ongoing retrospective analysis of patient motion of 160 DTI scans from 52 subjects (30 ataxic patients, 22 controls), we observed an average of 1.9 ± 0.87 mm of translation and 0.53 ± 0.38° of rotation over a 3 min 57s time period (results presented in poster format Landman et al., 2006). Even with a healthy, well-behaved subject in this study, the specified and realized DW directions were different by 0.20 ± 0.07° within each session and a much larger difference occurred between sessions due to subject positioning (Fig. 1B).

To provide a benchmark for the differences between “poor” and “good” DW schemes, a direction-wise closest approximation to the tetrahedral 6 scheme was generated from the 30 direction PE scheme, denoted Tetra₆. The Appendix Table A2 lists the DW directions from the Jones30 scheme that were assigned to the Tetra₆ scheme. This scheme has been shown to have a high condition number (a measure of noise propagation in the estimation process), produce DTI results that exhibit prominent dependence on orientation, and lead to low SNR results (Skare et al., 2000).

To test if the partitioned subsets were a good approximation of the unconstrained optimal PE schemes with an equivalent number of directions, Monte Carlo simulation experiments were performed on diffusion tensor data derived from in vivo human brain observations. To minimize arbitrary effects of orientation, the unconstrained optimal PE schemes were first reoriented so that the directions were most similar (under mean squared error) to those of the partitioned subsets. One thousand tensors with FA> 0.25 were randomly selected from the DTI results of an in vivo human brain experiment. For each tensor, 100 simulated observations were created by applying the tensor model of diffusion to the set of DW directions for each DW scheme (Stejskal and Tanner, 1965). Gaussian noise was added in quadrature (e.g., in the frequency domain) to the DWIs and $b_0$ so that the resulting $b_0$ images had a SNR of 20 dB (10:1), which is approximately equivalent to the SNR observed in the experimental data. Diffusion tensors were then estimated from the simulated data for the partitioned subsets and unconstrained optimized schemes, respectively. Comparison of the mean FA and PEV using each DW scheme revealed mean differences between the unconstrained optimal PE schemes and the optimal PE partition schemes of (1) less than 0.005 in FA and (2) less than 0.5° in PEV orientation, which are both less than the differences observed in repeated acquisitions using a 30 direction DW scheme (Farrell et al., 2006). Thus, the optimal subsets (PE₆, PE₁₀, PE₁₅) are deemed appropriate surrogates for the unconstrained optimal schemes.

To ensure a scan time equivalent (STE) comparison (e.g., between DW schemes that could be achieved in the same amount of scan time), the acquired DTI data were grouped to form “composite acquisitions.” The composite acquisitions are groupings of DW volumes from various scans within a session, e.g., the PE₆ directions from the first 5 sets of scans in a session. The basis for comparison is 1 or 3 STE units which represent a total of 30 DW images+ 5 $b_0\text{s}$ or 90 DW images + 15 $b_0\text{s}$ respectively, no matter whether the DW images were unique or repeated observations. For example, the optimal 6 DW directions can be selected from scans 1 through 5, 6 through 10, and 11 through 15 respectively to yield three composite acquisitions with 6 DW directions at 1 STE. In each instance, 5 observations of each DW image are utilized. However, in order to construct a “composite acquisition” at 3 STE, the optimal 6 DW directions from all 15 scans must be combined. Similarly, this method produces 5 groupings of 10 DW directions at 1 STE and 1 grouping of 10 DW directions at 3 STE; 7 groupings of 15 DW directions at 1 STE and 2 groupings of 15 DW directions at 3 STE; and 15 groupings of 30 DW directions at 1 STE and 5 groupings of 30 DW directions at 3 STE.

Within session mean FA, MD, and PEV findings are defined as those obtained when all acquired in vivo data (15 acquisitions of 30 DW directions) from a session are utilized in a single diffusion tensor calculation. We note that directly combining raw DTI data from different scan sessions can be problematic due to changes in scanner calibration. Therefore, to further improve SNR and mitigate inter-session effects, the mean DTI contrasts from each session were averaged over the three combined analyses to produce “gold standard” contrasts.

Experiments that use acquired MR data only

Experiments that used acquired MR data and modeled noise

The effects of SNR on DTI contrasts have been demonstrated with numerical simulation and by combining several acquired in vivo datasets. These mutually exclusive techniques are difficult to reconcile, as in vivo datasets ensure a myriad of effects. In this study, a series of three simulations were conducted that added modeled noise to acquired MR data to explore how well the tensor model explained and replicated the findings obtained when acquired MR data are used and a data grouping method is implemented to increase SNR.

The following Monte Carlo framework was used to simulate imaging results on a set of baseline tensors for all simulations: (1) a theoretical noise-free intensity was computed according to the tensor model for each diffusion- and minimally weighted image specified by the simulated acquisition scheme (Stejskal and Tanner, 1965); (2) Gaussian noise was added to the acquired MR data in quadrature so that the composite, simulated minimally weighted acquisitions had a specific SNR; and (3) simulated observed tensors were estimated with a log-linear model. The baseline tensors were processed sequentially with a fixed number of Monte Carlo iterations per tensor.

Experiments that used acquired MR data only

PEV colormaps computed from PE₆, PE₁₀, PE₁₅, and Jones30 schemes at 1 STE are shown in Fig. 2, at the level of the lateral ventricles. Though these PEV colormaps appear to be quite similar and of comparable quality, visual inspection alone cannot discern the substantial differences due to the different DW schemes. For quantitative analysis to detect the differences, ROI and voxel-based analyses were performed.

Experiments that used acquired MR data and modeled noise

Selecting a diffusion weighting scheme

DTI results are systematically and statistically different when comparing studies using different, well-balanced DW schemes, a result that is in agreement with simulations. However, the magnitude of these differences is less than might be expected from other experimental parameters, such as changes in the number of DW images, number of averages, physiological noise, scanner hardware, or other factors that might influence the SNR. At 1 STE for example, using a DW scheme with 6 rather than 30 directions results in a differential orientation bias of approximately (1) 0.02 [FA] for FA, (2) 0.02 × 10⁻³ mm²/s for MD, and (3) 2° for PEV orientation. In order for these effects to be of clinical concern, they should be of the magnitude or greater than intra-session and inter-session reproducibility effects. In comparison for 1 STE with 30 directions, the intra-session test–retest variability of FA is ~ 0.1 [FA], MD is ~ 0.05 × 10⁻³ mm²/s, and PEV is ~ 5° (Farrell et al., 2006). In other words, compared to test–retest reproducibility and impact of SNR, the effects of the DW schemes are minor. We conclude that analyses from studies with different, but well-balanced DW scheme protocols are comparable.

When planning DTI studies, it is important to separate issues related to precision (reproducibility) and accuracy (bias) of each scheme. The optimized PE₆, PE₁₀, PE₁₅ and Jones30 schemes tested in this study have comparable precision (e.g., Δσ in Fig. 8). This means that they have comparable powers to discriminate normal from abnormal. For group studies (ROI-based or voxel-based), these schemes should perform equally well. The bias becomes an issue if one needs to compare results between different DW schemes. If datasets are simply compared with values reported in the literature, the bias reported in this study should be a minor influence. However, if results from two different schemes are statistically compared within one study, great care should be exercised. The difference due to the PE₆ and Jones30 DW schemes could be non-negligible when a group study is performed in brain regions with a specific fiber orientation, and the fiber is oriented in such a way that the diffusion weightings observed with the PE₆ and Jones30 difference are maximally different. For example, suppose that 10 subjects are scanned once with each of two protocols, once with Jones30 and again with five repetitions of PE₆. If one assumes good registration and anatomical correspondence, FA inter-session variability is ~ 0.1 [FA] (Farrell et al., 2006) with 1 STE (30 DWIs + 5 $b_0\text{s}$, approximately 3.5–4 min of scan time) and maximum orientation bias difference is ~ 0.03 [FA] (Fig. 7). A simple power analysis reveals that there is an 11% chance that a t-test would detect at least one significant difference due the bias (Altman, 1999). If 50 subjects are studied, the likelihood of detection increases to 54%. The amount of the bias decreases as the SNR increases (Fig. 11). Therefore, the difference could be harder to detect as the scan time increases. These conditions are met when the same subject is scanned on the same scanner. In reality, variability in group studies is dominated by variability among the subjects, scanner performance and image registration quality, which is not included in the above power analysis. Nonetheless, it is always advisable to use the same imaging parameters as much as possible for group studies. The choice of a diffusion weighting scheme should be recognized as one of the important parameters that should be considered.

In this study, cardiac gating was not applied. It is well known that cardiac pulsation can have a profound effect on DTI results (Golay et al., 2002; Jones and Pierpaoli, 2005), which could be a major contribution to the intra- and inter-session reproducibility noted in this study. Studies with cardiac gating should have better reproducibility and the relative significance of the PE₆ −Jones30 difference, described above, may increase in such studies. Unfortunately, cardiac gating has not been widely accepted in clinical studies because it lengthens the scanning time and subject motion conditions and arrhythmia can occasionally cause gating instability.

When a protocol is being designed for a new study, a DW scheme with a larger number of gradient orientations should be considered because of the decreased directional sensitivity (little orientation dependence for precision and accuracy) even though the improvement could be minor, as has been demonstrated in this and previous (Jones, 2004) studies. However, schemes with fewer DW directions may have practical benefits that are not considered in this paper such as more efficient diffusion weighting. Repeated measurements (e.g., five repetition of PE₆ as opposed to one Jones30 measurement) allow variability assessment among redundant datasets, which carry important information about image corruption and coregistration errors. Another example is that the observed signal can be increased by minimizing TE (Pierpaoli et al., 1996) or by applying a more efficient diffusion weighting scheme (Basser and Pierpaoli, 1996). For example, by changing from a PE₆ scheme to a Tetra₆ scheme, TE is reduced from 88 ms to 80 ms which corresponds to an 11% signal increase for white matter (assuming an approximate T2 of 80 ms). However, tradeoffs occur at the expense of increased condition number, hence noise propagation, and the benefit must be evaluated on a case by case basis (Batchelor et al., 2003; Skare et al., 2000). For example, the orientation variability for the PE₆ results at 1 STE is lower than it is for the Tetra₆ results at 3 STE, even though the SNR in the latter case is approximately $\sqrt{3}$ greater (Fig. 8). Such a large increase in SNR cannot be achieved by TE shortening, indicating that the PE₆ scheme will always yield better results than the Tetra₆ scheme.

A simple multivariate log-linear least-square fitting method was used in this study. Recently, other fitting methods such as robust tensor estimation methods (Chang et al., 2005; Cox and Glen, 2006; Koay et al., 2006) have been introduced. Differences due to robust and simple tensor estimation methods are systematic (Koay et al., 2006), so caution must be used when comparing results from studies using different methods. It would be an interesting area of future research to evaluate if robust tensor estimation methods have consistent impacts on bias and variability depending on tensor orientation or DW scheme.

The close agreement of experimental and simulation results reinforces the general purpose utility of single tensor simulation models for investigating the performance of DTI protocols. The correspondence between the acquired and simulated data was especially strong for the variability measures and the patterns of orientation dependence. There is a possible discrepancy between Figs. 8 and 11. Namely, while the simulated data (Fig. 11) indicate that the bias between PE₆ and Jones30 decreases as SNR goes up, the acquired data (1 STE versus 3 STE) indicate that the bias persists with a high SNR dataset (3 STE). We suspect that the bias at 3 STE (between the PE₆ and Jones30 schemes) may not be mitigated by SNR and is due to non-tensor behavior (e.g., crossing fibers) seen in acquired MR data. Recently, Behrens et al. (2006) demonstrated that approximately one third of brain voxels satisfy the criteria for having multiple fiber populations. If there are a significant number of voxels with the non-tensor behavior, it is reasonable to expect that differences between PE₆ and Jones30 results manifest as bias with respect to high SNR data because the gold standard was determined by using all data acquired with Jones30 scheme. This difference should be sensitive to b-value (the higher the b-value, the more pronounced the effect Frank, 2002). Nonetheless, the effect of non-tensor behavior is not considered in this article and is a fascinating future project.

Orientation analysis

In addition to conventional ROI-based analysis, we performed experiment (Fig. 7) and simulation-based (Fig. 9) orientation analysis. Even though many of the findings in this study simply confirmed results from previous simulation studies, we believe that it is an important effort to reconcile previous studies. FA, MD, and PEV showed consistent orientation patterns of bias and variability in both experiments and simulation data, as summarized in Table 1. FA precision and accuracy showed the greatest spatial dependence on the DW scheme, with a decrease in FA and increased variability roughly corresponding to DW directions. MD measurements were insensitive to DW scheme and thus to tensor orientation, while PEV measurements showed some dependence on DW scheme. The differences between the effect sizes and precise boundaries between the experimental and simulated data can be attributed to (1) spatially varying SNR, (2) differences between specified and realized DW directions, and (3) partial volume effects. As the sign, magnitude, and spatial patterns of the effects were highly consistent between simulation and acquired MR experiments, the data suggest that the single tensor model is useful for characterizing the behavior of differences between DW schemes.

Increasing the directional resolution of the DW scheme decreased the range of the orientation dependence for all contrasts (Fig. 8). This means that a DW scheme with more directions has a smooth precision and accuracy profile, which does not change profoundly as a function of fiber orientation. On average, low directional resolution showed increased FA and MD bias and variability. These effects are seen in Fig. 8 at both 1 and 3 STE ($b_0$ SNR approximately 18 dB (7.9:1) and 20 dB (10:1)).

Single tensor simulations (Fig. 10) provide intuition on how aggregate profiles develop by considering mean behavior over ROIs. Orientation dependence is dominated by moderate and high FA tensors, as these tensors exhibit the greatest degree of orientation variability. FA bias decreased and variability increased for tensors aligned with the DW directions for PE₆ relative to Jones30 (Fig. 10). Neither MD bias nor variability showed strong orientation dependence. The relationship between PEV, DW direction, and estimation reliability was complicated, however, in general AV tended to decrease and AB increase when the PEV was aligned with a DW direction.

The SNR tensor simulations demonstrate how these interactions generalize to other biological contexts or SNR levels (Fig. 11). At less than 15 dB (5.6:1), tensor estimation is exceptionally noisy, eigenvalues are notably biased, and there is substantial angular dependence introduced for the 6 DW direction scheme. At greater than 30 dB (31.6:1), neither minimal bias, orientation dependence, nor DW scheme dependence is apparent in either the FA or MD. Therefore, low SNR tends to exacerbate the orientation dependence of DTI contrast estimation and differences due to DW schemes.

Conceptualizing results

The results of this study are in agreement with Jones (2004) in that the error orientation profiles can be thought of “as a rubber sheet, and the sampling vectors as ‘fingers’” that serve to even out the surfaces. In other words, more independent sampling directions result in less orientation dependence in the precision and accuracy of derived quantities, when the total number of DW directions is held constant. Yet, these “fingers” are a bit more complex than weights on a surface serving to improve estimation.

The choice between sampling at independent directions versus repeated directions can be thought of as a tradeoff between shape (e.g., FA and MD) and orientation (e.g., PEV). The potential information about a tensor contained in a set of sampled DW directions is dependent on the tensor and noise. The results from this study (both simulation and in vivo measurements) show that estimation accuracy and precision of contrasts are not necessarily optimal for tensors aligned with the DW directions. This can be intuitively appreciated by considering a “diffusion peanut” for a prolate tensor (Jones and Basser, 2004). The physical diffusion process gives rise to a 3D diffusion profile, which is then sampled by DW imaging. The diffusion peanut is the conceptual representation of the diffusivity measured in each direction. In the case of a prolate tensor with high anisotropy, the derivatives of the spatial profiles of the diffusion peanut are relatively low at the poles and the equator. Sampling several locations (e.g., acquiring DW images with diffusion sensitization along these directions) in these regions would tell little about the orientation of the tensor, yet these would prove useful in determining anisotropy. Conversely, between the poles and equator, the profile of the diffusion peanut changes more rapidly, which well indicates orientation, but tells little about anisotropy. Furthermore, when noise is considered, different sampling directions may be determined with different accuracies due to varying signal attenuation and/or spatially varying SNR (e.g., due to parallel imaging techniques).

The results from this study agree with this analogy as we note that orientations of improved FA accuracy correspond to the orientations of decreased PEV accuracy. This effect can be thought of as a tradeoff between measuring the shape of the tensor and the orientation of the tensor. The high directional resolution scheme minimizes the variability for a tensor of unknown orientation, while a change to a low directional resolution scheme has varying impacts based on the orientation of the underlying tensor. If a DW scheme over samples slowly changing regions of the diffusion profile, the eigenvalues (hence, FA) will be relatively better determined. While if the rapidly changing regions of the diffusion profile are over sampled, the orientation (hence, PEV) will be relatively better determined. In light of this conceptual model, the experimental (Fig. 7) and simulation (Fig. 10) results are not surprising.

There clearly is inhomogeneity (orientation-dependent variation) in bias and variability. The inhomogeneity is not random and is related to sampling orientations. For FA, the maxima (or minima) correspond well to the sampling orientations, but for MD and PEV, the maxima or minima agree less well with the sampled orientations. In some cases, the perception of correspondence may be a function of the window-level of the figures. However, the phenomena are likely real, and may be explained by the contributions of all sampling directions and the secondary and tertiary eigenvectors.