|Home | About | Journals | Submit | Contact Us | Français|
A key question in diffusion imaging is how many diffusion-weighted images suffice to provide adequate signal-to-noise ratio (SNR) for studies of fiber integrity. Motion, physiological effects, and scan duration all affect the achievable SNR in real brain images, making theoretical studies and simulations only partially useful. We therefore scanned 50 healthy adults with 105-gradient high-angular resolution diffusion imaging (HARDI) at 4 Tesla. From gradient image subsets of varying size (6≤N≤94) that optimized a spherical angular distribution energy, we created SNR plots (versus gradient numbers) for seven common diffusion anisotropy indices: fractional and relative anisotropy (FA, RA), mean diffusivity (MD), volume ratio (VR), geodesic anisotropy (GA), its hyperbolic tangent (tGA), and generalized fractional anisotropy (GFA). SNR, defined in a region of interest in the corpus callosum, was near-maximal with 58, 66 and 62 gradients for MD, FA and RA in respectively, and with about 55 gradients for GA and tGA. For VR and GFA, SNR increased rapidly with more gradients. SNR was optimized when the ratio of diffusion-sensitized to non-sensitized images was 9.13 for GA and tGA, 10.57 for FA, 9.17 for RA, and 26 for MD and VR. In orientation density functions modeling the HARDI signal as a continuous mixture of tensors, the diffusion profile reconstruction accuracy rose rapidly with additional gradients. These plots may help in making trade-off decisions when designing diffusion imaging protocols.
A common question in diffusion imaging is how many gradient images should be collected to provide adequate signal-to-noise ratio (SNR) for studies of fiber integrity in the brain. To maximize the information collected, several research groups have proposed high-angular resolution diffusion imaging (HARDI) - a powerful extension of MRI - that applies diffusion-sensitized gradients to the brain in large number of directions (Tuch, 2004). Diffusion-weighted images are collected by applying a magnetic field gradient in a specific orientation to quantify water diffusion in that direction. By studying anisotropy (directional dependency) in these diffusion profiles, one can gain exquisite insight into local fiber orientation and integrity in living tissues such as the brain and heart. Most early diffusion anisotropy studies used the diffusion tensor model (Basser et al., 1996), which describes the anisotropy of water diffusion in tissues by estimating the 3×3 covariance matrix of a Gaussian distribution from a set of K diffusion-sensitized images. Each voxel’s signal intensity in the k-th image is decreased by water diffusion according to the Stejskal-Tanner equation (Stejskal et al., 1965): Sk = S0 exp [-bgk T D gk], where S0 is the non-diffusion-weighted signal intensity, D is the 3×3 diffusion tensor, gk is the direction of the applied diffusion gradient and b is the instrumental scaling factor containing information on the pulse sequence, gradient strength, and physical constants (Le Bihan 1990a, 1990b). Seven gradients (6 diffusion-sensitized gradients and one diffusion non-sensitized gradient) are mathematically sufficient to determine the diffusion tensor, but MRI protocols with higher angular and radial resolution (e.g., HARDI) are commonly used, to resolve more complex diffusion geometries that a single tensor fails to model, e.g., fiber crossings and intermixing of tracts. The greater angular resolution is useful when examining fiber trajectories using tractography (Lenglet et al., 2008). Routine clinical studies of fiber integrity may also benefit from HARDI, as DTI-derived measures of anisotropy are biased (too low) in the ~40% of white matter voxels where more than one dominant fiber direction is detectable (Zhan et al., 2009a, Leow et al., 2009a).
Recent technical advances have made HARDI more practical. A 14-minute scan can typically sample over 100 angles (with 2 mm voxels at 4 Tesla). HARDI’s improved signal-to-noise ratio can be used to infer fiber pathways in the brain, identifying anatomical features, connections and disease biomarkers not seen with conventional MRI. If more angular detail is available, fiber orientation distribution functions (ODFs) may be reconstructed from the raw HARDI signal using several alternative methods: q-ball reconstruction based on the Funk-Radon transform (Tuch 2004; Hess et al., 2006), the Persistent Angular Structure technique and Maximum Entropy Spherical Deconvolution (MESD) (Jansons et al., 2003; Anderson, 2005), the DOT transform (Özarslan et al., 2006), deconvolution-based techniques (Tournier et al., 2004; Alexander, 2005; Kaden et al., 2007), or by fitting continuous mixtures of tensors (Leow et al., 2009; Jian et al., 2007). Recent work on stochastic tractography (Pichon et al., 2005; Jbabdi et al., 2007) has exploited the increased angular detail in HARDI, and fluid registration methods have been designed to align HARDI ODFs using specialized Riemannian metrics, such as the Fisher-Rao metric, to quantify the discrepancy between density functions (Chiang et al., 2008a).
The question “how many gradients are sufficient for diffusion imaging” has been studied theoretically and using computer-simulated models of signals. Papadakis et al. (2000) analyzed the performance of various DT-MRI angular sampling schemes in terms of the variance of different anisotropy measurements, and concluded that the minimum number of unique encoding directions N(bx) required for robust anisotropy estimation was around 18-21. Even so, Hasan et al. (2001) found that more than 6 sampling orientations provided little added benefit when computing tensor-derived measures, so long as the selected orientations were optimized to point to the vertices of an icosahedron. These findings received later support from the work of Batchelor (2002) but were challenged by Jones (2004). Jones found, in his simulation experiments, that at least 20 unique sampling orientations are necessary to estimate anisotropy robustly, and at least 30 unique sampling orientations are required to estimate tensor orientation and ADC robustly.
These guidelines are useful and highly cited, but many prior studies have been based on theoretical models or simulations of signal and noise, and may not reflect the achievable SNR in real human subjects. Initial studies by Landman et al. (2007) and Farrell et al. (2007) have begun to examine FA, MD, and principal eigenvector (PEV) accuracy, precision and reproducibility in real human data at 1.5T, by scanning a single subject with many alternative protocols (see Discussion for details). Here we extend these analyses to study a DTI protocol with more directions (94 rather than 30) and also in more subjects (50 rather than 1, although the one subject studied by Landman et al. (2007) and Farrell et al. (2007) was scanned using many protocols). In addition, we report SNR analysis for a HARDI-type index, GFA. GFA is derived from the full set of diffusion gradients and does not rely on fitting a diffusion tensor model to the signals. We also studied the recently introduced metric of geodesic anisotropy. Our study is the first, to our knowledge, to study DTI SNR at 4 Tesla, as prior studies were conducted at 1.5 T (Landman et al., 2007; Farrell et al., 2007). A very recent study of DTI SNR in 13 subjects (Zhu et al., 2009; also at 1.5T) appeared since this paper was submitted, and we now summarize that paper in the Discussion.
In real brain images, the achievable SNR is affected by motion, susceptibility-related distortions, physiological effects on the MR signal, and even scan duration. With these additional sources of error, the achievable SNR may asymptote more quickly than predicted, with respect to the gradient numbers, and may even depend on the subject. As such, a population-based study scanning a reasonably large number of subjects (here N=50), can help to assess whether real data follow theoretical predictions, whether SNR profiles depend on the subject scanned, and whether some DTI-derived measures are more robust than others when there is only enough time to collect limited numbers of gradients. Prior simulation studies focused on a few standard DTI-derived measures (FA, ADC, MD, and the principal eigenvector direction, PEV; Landman et al., 2007; Farrell et al., 2007). Here we also examined several other popular tensor-derived parameters, including RA, VR, GA and tGA, which have been somewhat widely used in clinical studies, and generalized fractional anisotropy (GFA). GFA is computed from the full gradient set (without fitting a tensor) and has been advocated as one possible HARDI analog of DTI-derived FA. We also assessed the optimal ratio of diffusion-sensitized to non-sensitized images, N(bx)/N(b0), to maximize SNR. Finally, to support our premise that the higher SNR protocols are also more accurate, we illustrated how the reconstruction accuracy of the diffusion orientation density function (ODF) depends on the applied gradient set.
Although the optimal trade-off between SNR and scan time depends on the time available and the need for angular precision, our goal was to provide useful graphical plots suggesting how to determine optimal gradient numbers for tensor-derived diffusion indices including mean diffusivity (MD), fractional anisotropy (FA), relative anisotropy (RA), volume ratio (VR), hyperbolic tangent of geodesic anisotropy (tGA), as well the generalized fractional anisotropy (GFA) derived from the ODF.
3D structural brain MRI scans and diffusion-weighted scans were acquired from 50 healthy adult participants (mean age: 22.6 years, S.D.: 2.1, 25 men/25 women) on a 4 Tesla Bruker Medspec MRI scanner using an optimized diffusion tensor sequence (Zhan et al., 2009a, 2009b). None of the subjects reported a history of significant head injury, a neurological or psychiatric illness, substance abuse or dependence, or had a first-degree relative with a psychiatric disorder. In addition, all subjects were screened, using a detailed neurocognitive evaluation (see de Zubicaray et al., 2008 for an overview), to exclude cases of pathology known to affect brain structure.
Imaging parameters were: TE/TR 92.3/8250 ms, 55 × 2mm contiguous slices, FOV = 23 cm. 105 gradient images were collected: 11 baseline (b0) images with no diffusion sensitization (i.e., T2-weighted images) and 94 diffusion-weighted images (b-value: 1159 s/mm2) in which gradient directions were evenly distributed on the hemisphere (Jones et al., 1999). The reconstruction matrix was 128 × 128, yielding a 1.8×1.8 mm2 in-plane resolution. Total scan time was 14.5 minutes.
Since each gradient image is applied in a direction that may be represented as one point i (1≤i≤94) on the surface of the unit sphere (Figure 1a), the optimality of a gradient directions set is typically defined using PDEs based on electrostatic repulsion, or angular distribution energy (see Hasan et al. (2001) for an excellent review). The parameterized normalized analytical set encoding gradient directions [gx, gy, gz] may be written as in Equation 1 (Wong et al., 1994)
In each of these encoding directions, diffusion-sensitized magnetic field gradient is applied. If the number of diffusion-sensitized gradients is large enough, the placement of points should be spaced nearly uniformly on the sphere. To evaluate the effect of varying the number of diffusion-sensitized gradients, N(bx), a series of subsets, of size 6≤N≤94, were chosen from among the original 94 diffusion-sensitized gradients, to maximize the total angular distribution energy. The angular distribution energy, Eij , of a pair of points, i and j, on the unit sphere may be defined as the inverse of the sum of the squares of (1) the least spherical distance between point i and point j, and (2) the least spherical distance between point i and point j’s antipodally symmetric point J. (Equation 2):
Here, i,j represent two different points on the spherical surface, J is the antipodally symmetrical point to j, and dist(i,j) is the least spherical distance between point i and point j. dist(i, J) is the least spherical distance between point i and point J (Figure 1b):
The total angular distribution energy, EL(N), for a specific subset of N diffusion-sensitized gradients may then be defined as the sum of the angular distribution energy of all pairs of unit gradient vectors on the sphere, using least spherical distances (Equation 3). The optimal sampled gradient subset N* can be achieved by selecting the set that maximizes E(N).
Gradients were added sequentially, and all the direction vectors present in the subset with N gradients were also present in the subsets with N+1, N+2, etc., gradients, so the sets form a nested chain (in the language of set theory). There are a couple of justifications for doing this. First, it is very efficient to search for the best gradient to add if the first N gradients are fixed. If the initial set of N gradients is fixed, there are only 94-N choices of gradients to add to the set, so the nested sets may be generated efficiently using a double loop where the associated energies are ranked and the best one retained. If we also allowed the first N gradients to vary, we would be faced with an N-dimensional combinatorial optimization with an initial search space of size 94CN for each successive N. Such an N-dimensional problem is tractable if the points on the sphere are each allowed to be continuously updated by gradient descent on the energy, and there are several such implementations of such an algorithm that run quickly. But because we wanted to avoid interpolation between the existing gradients, we wanted to choose from the available set of discrete directions, and it not clear how such a choice can be obtained by iterative optimization, unless non-global numerical optimization methods are used such as simulated annealing. Although there is some loss of optimality in our inductive approach, we found (see Figure 2) that the error is bounded below by 8% and is usually far less.
Thus, to avoid a time-consuming exhaustive search, we first chose one seed point from the original 94 diffusion sensitized gradient directions. Without loss of generality, we chose this to be (1, 0, 0), and found another point from the remaining 93 diffusion-sensitized gradient directions to maximize EL(2), then kept increasing the angular sampling one gradient at a time to maximize EL(3), EL(4), †,EL(93). EL(94) represents our original HARDI105 data. Gradient subsets with fewer than 6 independent directions will be ignored, as 6 diffusion-sensitized gradients are the minimum required for tensor estimation. In this initial test, we retained all 11 non-diffusion-sensitized gradients (b0); later, we assessed the effect of varying these as well, but for now they were kept constant.
In addition, it was fairer to use subsets of the original data rather than interpolate new data at new angular locations in between existing ones already collected, in case the interpolation itself influences the SNR. In addition, the selection of the first point may have a small influence on the obtained sampling: rotating the coordinates of the first point would lead to a rotated version of the optimal sampling. In specific brain regions, this may lead to a small local bias depending on the main orientation of the observed tissue, but this would likely be a very minor effect once a moderate number of gradients were available, and the effect would even out across the brain (as different fibers have a large range of different principal orientations).
Although increasing the number of diffusion-sensitized gradients samples the diffusion profile in more detail at each point, it is common to collect several baseline images, with no diffusion sensitization, which serve as a reference signal. Relative to these images, the attenuation of the MR signal due to diffusion is estimated. Collecting more of these images takes time but also better evaluates the baseline signal, so we also assessed the optimal ratio of the number of diffusion-sensitized images, N(bx), to the number of non-diffusion sensitized images, N(b0). To fairly compare scans that take the same amount of time, we made sure that N(b0)+N(bx) was held constant, and this constant was set to 21, 41, 61, and 81, in our experiments. By varying N(b0) from 1 to 11, we were able to associate each baseline set of images with a corresponding optimized subset of N(bx)=N-N(b0) diffusion-sensitized images, based on maximizing the total angular distribution energy, EL(N(bx)).
For each gradient subset, a diffusion tensor was fitted to the raw diffusion-weighted data at each voxel. The diffusion tensor was estimated using MedINRIA, which is a DTI analysis program developed by the INRIA research project Asclepios. (http://www-sop.inria.fr/asclepios/software/MedINRIA). (Pennec et al., 2004) Based on the estimated tensors, several commonly used scalar measures of diffusion anisotropy were computed from the tensors’ eigenvalues (λ1, λ2, λ3), including Fractional Anisotropy (FA), Relative Anisotropy (RA), Mean Diffusivity (MD), Volume Ratio (VR), Geodesic anisotropy (GA) and its hyperbolic tangent (tGA; it is common to apply this transform to GA so that the transformed values have the same ranges as the more commonly used measures, i.e., 0 to 1 for FA). All these scalar parameters were defined using the formulae listed in Table 1.
Orientation distribution functions (ODF) for water diffusion were also computed voxelwise from the HARDI signals using the Funk-Radon Transform (FRT) (Tuch, 2004). We used Descoteaux’s method (2007), which first expands the HARDI signals as a spherical harmonic (SH) series, simplifying the FRT to a linear matrix operation on the coefficients. To estimate the SH coefficients, we set the order of the SH series to 6. A generalized fractional anisotropy (GFA) map (Tuch, 2004) was constructed from the ODF (the formula is given in Table 1).
To assess the image SNR in a consistent region of the brain across all subjects, we used a region of interest (ROI) approach. As noted in Landman et al. (2007), the SNR for in vivo studies varies spatially, due to differences in coil sensitivity, and is tissue-dependent (depending on tissue T1, T2, and PD). To define a region that would be consistent with past studies (e.g., Landman et al., 2007) and would cover a region that is the target of many DTI studies, we chose a small region of interest in the corpus callosum, the main interhemispheric fiber tract. To define the ROI as consistently as possible across all 50 subjects, we first registered all 50 subjects’ b0 images to a common image template since these T2-weighted (b0) images are collected in perfect register with the diffusion tensor data. A single subject’s T2-weighted scan was randomly selected to serve as the registration target for all subjects. All other subjects’ T2-weighted (b0) images were aligned to this template using the publicly available FLIRT registration software (http://fsl.fmrib.ox.ac.uk/fsl/flirt), using 9-parameter linear registration and a correlation ratio cost function. Once linearly aligned, further registration was performed using a 3D Navier-Stokes-based fluid warping technique enforcing diffeomorphic mappings, using least squares intensity differences as a cost function (Lepore et al., 2008). Registration parameters, including the number of time steps, were optimized for the T2-weighted contrast images. Transformation matrices from the linear registration and deformation fields from the fluid registration were retained and were applied to all subjects’ diffusion anisotropy index maps. After registration, a small ROI (10×10×10 mm3) was manually traced in the middle of the corpus callosum (CC) on the common template using the “erosion with a 3D diamond” tool in the BrainSuite software package (Shattuck et al., 2002), to obtain a small and relatively homogeneous region. In this way, the biological variation in diffusion signals across the brain does not confound the assessment of SNR with respect to the number of gradients.
We computed SNR for all diffusion anisotropy indices map for all 50 subjects and for all angular sampling schemes EL(6) to EL(94). As simulation studies indicate that SNR tends to increase with the number of diffusion-sensitized gradients, we first plotted SNR against gradient numbers for all the parameter maps, by using all optimized gradient subsets. Although other definitions are possible, each map’s SNR was defined as the ratio of the mean voxel value to the standard deviation of the voxel values in the ROI. We preferred this method instead of estimating the noise from a background region of interest because it made more sense to estimate the signal and noise from the same region, and also the eigenvalues estimated from outside brain are more likely to be negative, which made no sense for the diffusion theory.
We also assessed ODF reconstruction accuracy for our different gradient subsets, using a computer simulation to define ground truth. This part of the study was intended to illustrate the effects of angular sampling on the ODF, and to help explain the slightly unexpected behavior of SNR versus gradient numbers for the GFA measures, which were derived from the full ODF and not from the single-tensor approximation.
We created 1000 different models of two-fiber systems, crossing at 90 degrees with equal volume fractions (w1 = w2 =0.5). Here we chose λ1=10×10−10 m2s−1 and λ2=2×10−10 m2s−1 as the eigenvalues for each individual tensor (FA=0.77, typical for white matter) to generate simulations using discrete mixtures of Gaussian distributions:
Using the gradient subsets mentioned in Section 2.2, we sub-sampled the simulated HARDI data with 94 values for S(q)/S(0), and the tensor distribution function (TDF) method (Leow et al., 2009) was used to analyze the sub-sampled data. The TDF method is essentially a deconvolution method, which we recently proposed, that estimates the angular distribution of the underlying fibers. We denote the space of symmetric positive definite three-by-three matrices by . The probabilistic ensemble of tensors, as represented by a tensor distribution function (TDF) P, is defined on the tensor space that best explains the observed diffusion-weighted images:
To solve for an optimal TDF P*, we use the multiple diffusion-sensitized gradient directions qi and arrive at P* using the least-squares principle:
These ODFs were rendered using 642 point samples, determined using a seventh-order icosahedral approximation of the unit sphere. To assess how accurately the diffusion profiles could be reconstructed in situations based on different angular sampling schemes, the symmetrized form of the Kullback-Leibler (KL) divergence, a commonly used measure from information theory, was used to measure the discrepancy between the reconstructed and ground truth ODFs. Reconstruction error was calculated from Eq.8, in which p(x) is the ODF derived from the sub-sampled schemes, and q(x) is the original ODF derived from the ground truth data:
In picking sequences with fewer directions than the total number acquired (94), we were faced with a dilemma: either we would have to re-scan people with a new energy-optimal direction set for each lower value of N (too time-consuming), or we would have to find a new energy-optimal direction set and interpolate the data, on the sphere, from the initially acquired directions. As the latter approach would require some interpolation, and therefore smoothing, which might affect the SNR relative to the data initially acquired, we wanted to take a simpler approach that inductively removed or added gradients in the initial set. As this set was very dense on the sphere, the approximation of the remaining subsets to an energy-optimal set is relatively close, given other considerations. Even so, it is possible that the subsampling of gradients from a given set would give some deviation from the theoretically optimal gradient set, and that the error might be different for each number selected, and reduced as the number increases (see Landman et al., 2007 for a discussion of this effect). This is an important consideration, as the SNR gains with increasing gradient numbers needs to be not entirely attributable to the deviations from optimal sampling. Figure 2 shows that even for gradient sets with very low numbers of directions (6), the energy deviation between the theoretically-optimal and generated subsets is less than 8%, and is only around 3-4% when there are more than 25 directions. As such we acknowledge that sub-optimal sampling may play a minor role in the SNR gains with respect to gradient numbers. Even so, the gradient sets were verified (in Figure 2) as being close to optimal, and as well-distributed as could be achieved without acquiring additional gradient sets in new directions.
We note that other subsets of gradients could be selected based on alternative criteria – there are numerical optimization equations (based on variational calculus) for uniformly sampling a sphere, such as Minimum Force (MF), Minimum Energy (ME), and Minimum Condition Number (MC), and minimum tensor variance (MV) (see Hasan, 2000; Hasan et al., 2001, for a review of these). The advantage of our optimization method is firstly that the angular distribution energy is monotonically increasing as the number of angular sampling points increases (Figure 2a). In addition, by including antipodally opposite points, we consider the distribution of the directions on the sphere rather than in just one hemisphere. Thus our optimization method may be regarded as related to the Minimum Energy (ME) method (Hasan, 2000). Also, to evaluate the associated errors for each subset generated from our optimization criteria, we plotted the total angular distribution energy for our inductively generated subsets (denoted by EL(N)) versus the total angular distribution energy for the corresponding optimal subsets (denoted by Eo(N)) calculated from Equation 1. The ratio of EL(N) to Eo(N) is plotted in Figure 2b. From Figure 2b, we can see that the error associated with our generated subsets relative to the theoretically optimal sampling is less than 8% even when gradient numbers are very small, and is less than around 3-4% when there are more than 25 directions.
In Figure 3, SNR curves are shown for MD, FA, RA, GA and tGA, as the number of diffusion-sensitized gradients, N(bx), increases from 6 to 94. For these experiments, the number of diffusion-sensitized images, N(b0), was held constant at 11 (this assumption is relaxed later). These plots may also be interpreted as plots of SNR versus acquisition time, as the scan time is proportional to the number of gradients, taking 14.5 minutes on our system for the full set of 105. In the plots, the mean SNR across all 50 subjects is shown, as is the SNR plot for the subject with the least SNR (blue curve) and for the subject with the highest SNR (red curve). SNR curves for the individuals resembled those for the population mean. As expected, there was a clear tendency for the SNR for all these five parameters to become saturated when the number of diffusion-sensitized gradients reached a particular level. Later, we define SNR convergence statistically, but for now we consider a simple measure of convergence, defined as the smallest number of gradients for which SNR comes within 1% of its maximal value. Table 2 shows this saturation number for each of these DTI-derived parameters.
As expected, SNR plots all rise initially as additional gradients collect more information for fitting the diffusion tensor. Rician noise in the gradient images leads to noisy estimates of the 3 tensor eigenvalues. SNR plots and saturation numbers differ for each DTI-derived parameter. This is because arithmetical operations on the eigenvalues such as addition and subtraction can cancel noise to some extent, whereas multiplication and division of eigenvalues tends to boost the noise level, and different types of arithmetic operations are used to calculate each DTI-derived parameter.
The saturation number was 58 for MD which is defined as the average of the 3 eigenvalues, whereas the FA and RA, which are defined as multiplication and division of tensor eigenvalues, both saturate with higher gradients (66 and 62 in respectively). Moreover, the GA (and its transform tGA), which are calculated from logarithm and/or exponentiation on the tensor eigenvalues, both saturate with fewer gradients (55~56). Even so, the hyperbolic tangent transformation can enhance SNR, because in general, mean SNR(tGA)>mean SNR(GA) (see Figure 3(d) and (e)). This supports the promise of tGA for clinical studies; our prior statistical maps comparing tGA with FA in a group of 90 twin subjects (Lee et al., 2008) found that tGA can provide slightly higher effect sizes than FA, perhaps because it accounts for the tensor manifold structure when computing differences and similarities between tensors.
Figure 4 shows the 50 subjects’ gains in SNR versus the number of diffusion-sensitized gradients for 2 additional measures of diffusion anisotropy: VR and GFA. The general tendency of mean SNR for these two parameters is to monotonically increase with the number of diffusion-sensitized gradients. SNR should eventually converge as the added value of new gradients diminishes, but convergence has not yet occurred even with 94 gradients - the highest angular sampling examined here - and there are gradual increases in SNR over the entire range. VR is computed from the eigenvalues of the fitted tensor, so it should eventually converge as N(bx) increases, but the number required to reach this plateau may be larger than 94. As seen in the formulae in Table 1, VR is expected to be noisier since it involves products of the 3 eigenvalues. Products of noisy variables tend to boost noise while the operations of addition and averaging tend to reduce it, and this can be confirmed from the curve of SNR(VR), which is more frustrating than other SNR plots. In general, our empirical findings are in line with theoretical predictions from error propagation analysis, in which measures relying on products of eigenvalues tensor to be noisier than summations where partially independent noise distributions can partially cancel each other out (Poonawalla et al., 2004; Koay et al., 2007). This can also be validated by the following simulation. Firstly, three positive eigenvalues were randomly produced between 0 and 1, and independent noise was added into the three eigenvalues, with σ(noise)=0.05. In comparing the DTI-derived indices before and after adding the noise, we confirmed that σ(MD) was the smallest (0.0781), while σ(VR) was the highest (4.5755), in line with the trajectories observed for the SNR.
Figure 4(b) shows that the mean SNR(GFA) has a generally increasing trend from 28 to 94 diffusion-sensitized gradients, but its SNR does not show evidence of saturation until N(bx) is close to 94. It is perhaps surprising how erratic this curve is compared to the others, especially when GFA has been advocated as a possible analog of FA for HARDI-based studies of fiber integrity. We note that GFA is the only index we examined that is computed from the orientation density function (ODF) rather than by first fitting a tensor. As such, it uses more of the information in the full set of diffusion gradients, and may only saturate when the number of gradients is very high (over 100). The GFA is calculated from the ODF, which is conventionally approximated by a spherical harmonic series (SH). In an SH series expansion of order k (here 6), the number of independent SH coefficients is (k+1)(k+2)/2, or 28 in our study. At least 28 independent diffusion-sensitized gradients are required to uniquely determine these 28 unknown SH coefficients, using least-squares fitting. As a large number of gradients are required to estimate the SH series stably, the many oscillations in the SNR plot for GFA may result from errors in the SH approximation. Intriguingly, in prior work using a 94-gradient HARDI sequence, we also found that FA provided better effect sizes than GFA did for correlations with IQ (intelligence quotient) in 92 normal subjects, although both measures were strongly correlated with IQ (Chiang et al., 2008, 2009). This is perhaps because the SNR is suboptimal for GFA unless more than 100 gradients are used. In addition, GFA is calculated from the ODFs, which are known to be noisy. To address this, additional ODF regularization and/or sharpening techniques have been proposed (Descoteaux et al., 2007, Goh et al., 2009a, 2009b) to smooth the signal and better detect the angular maxima. Although these post-processing operations may make it possible to smooth the erratic nature of the SNR (GFA) curve, here we just aim to give the reader a general impression of the SNR (GFA) curve without detailed subsequent post-processing.
To better understand whether SNR had converged, we performed a paired Student’s t test on SNR values, assessing the effects of adding one additional gradient from N-1 to N in all 50 subjects. If the t test result is declared significant (p<0.05) then this test confirms that adding one more gradient does indeed lead to better SNR. When this t test is not significant, there is no evidence that adding one more gradient to the acquisition protocol is helpful, so the signal may be said to be saturated. This result of this test clearly depends on the number of subjects (here 50), however, this is a reasonable and intuitive operational definition of saturation for practical purposes. We note that this test could be slightly improved by incorporating a multiple comparisons correction into the p-value, to control the false discovery rate, but we did not do so as the tests were intended as a heuristic to compare successive increments in gradient numbers. Figure 5 shows the paired t test results for MD, FA, RA, VR, GA, tGA and GFA.
Figure 5 shows the paired t test results for assessing the hypothesis that the SNR improves when adding one more gradient (i.e., SNRN>SNRN-1) for all the DTI-derived parameters. From Fig. 5(a)-(c) and (e)~(f), the SNR increase with increasing N(bx) is statistically significant, and then after some number is reached, the SNR oscillates and no longer shows a statistically significant improvement, consistent with Figure 3. The highest number of gradients for which the SNR improves when adding one more gradient is called the statistical saturation number, and this is shown in Table 3. We defined the meaning of this number to be that successively increasing the N(bx) always leads to statistically significant improvements in SNR until the statistical saturation number of gradients is reached.
Our saturation numbers are higher than comparable measures reported in other studies based on simulations, such as 18~21 by Papadakis et al. (2000), and 6 by Hassan et al.(2001). Fig. 5(d) confirm our previous results, shown in Fig. 4(a), that SNR values for VR increase with no obvious asymptote. As shown in Fig. 5(g), SNR(GFA) is so erratic for small numbers of gradients, that statistically significant improvements cannot be detected. This is because the SH approximation is extremely unstable when the number of gradients is low, making it difficult to make inductive arguments about SNR. Figure 6 shows some representative maps of FA, MD, RA, VR, GA and tGA based on HARDI17, HARDI38 and HARDI105 (where the numbers refer to the total number of gradient images, N(bx)+N(b0)).
Figure 7 shows the normalized SNR trends with the ratio of N(bx)/N(b0) in the condition of N(bx)+N(b0)=constant. We set this constant to be 21, 41, 61 and 81, and we have maximum N(b0)=11, so the ratio of N(bx)/N(b0) can vary from 0.9 to 80, and SNR was calculated for all parameters mentioned in the above section for different ratios. Table 4 shows the ratio that maximizes SNR for MD, FA, GA, tGA, RA and VR.
As noted earlier, the actual level of noise in the images depends on subject-specific factors (e.g., motion, physiology) and its effects are enhanced in measures based on products or ratios of eigenvalues, so the proportion of baseline images may differ when optimizing SNR for different DTI-derived measures. In real data, for all parameters except MD and VR, we found that SNR was optimized when using a slightly higher proportion of images with no diffusion sensitization. VR uses a higher-order product of eigenvalues, which tends to be noisier, so a higher proportion of diffusion sensitized gradients may be needed to overcome noise in the eigenvalue estimates. For MD, SNR value variation is within 1% of its maximum achievable value when the ratio is between 8 and 80, which indicates that MD may have more lenient requirements for the baseline signal. Considering its simple mathematical formula, which contains no products or ratios of eigenvalues, MD is relatively noise-resistant, requiring relatively few baseline and diffusion-sensitized images to estimate reliably.
This current paper is one of a few studies that have examined how SNR depends on gradient numbers in human data. Landman et al. (2007) studied the precision and accuracy of in vivo FA, MD, and principal eigenvector (PEV) estimates in the splenium of the corpus callosum and in the putamen. They scanned a single healthy 24-year-old male volunteer with a 30-direction diffusion-weighted protocol at 1.5 Tesla, using several optimized direction sets that optimized a potential energy measure derived by Jones et al. (1999), and models based on the Platonic solids. Some direction sets led to upward biases in the estimates of FA, and increased variability in MD; they also found that some of the sparser direction sets were more accurate in estimating FA for tensors in certain orientations: the orientations with improved FA accuracy were also the orientations with decreased PEV accuracy. These biases are of interest for multi-subject studies, as they imply that group analysis studies cannot achieve unbiased high-SNR results by averaging data across subjects, when each subject’s dataset is collected with low SNR (Farrell et al., 2007). A related study by Farrell et al. (2007) examined the same single-subject multi-protocol dataset, and noticed an upward bias in estimates of FA but no detectable bias in estimates of MD or the PEV direction, as SNR decreased. Higher SNR was needed to measure FA precisely in the gray matter versus the white matter. Farrell et al. (2007) also found that the inter-session reproducibility of the DTI-derived signal, over time, is much higher when more gradients are collected (either more directions, or repeated acquisitions, which both tend to increase SNR). This is a key point for longitudinal, interventional, or multi-site studies using DTI, where the changes in DTI signal may be of as much interest as the value of the DTI-derived parameters at any single time-point. A very recent paper by Zhu et al. (2009) scanned 13 healthy volunteers, at 1.5 Tesla, with protocols that contained 6, 21, and 31 unique gradient directions, with the latter two sequences being obtained using minimum potential energy computations for points on a sphere. They examined 15 ROIs in the white matter, and found that the uncertainty in the FA and MD measurements varied spatially across the brain, and varied systematically for different values of FA. They also examined, as we did here, the optimal ratio of diffusion-weighted to non-diffusion-weighted images (N(bx)/N(b0)), which they termed the DTIR (for DT image ratio). In simulations that examined DTIR values ranging from 2.89:1 to 61:1, they found that a DTIR value in the range of 5-6 gave well-balanced measurement quality for FA, MD, and PEV for most applications that have to operate under the constraints of a pre-defined scan time. This is largely in agreement with our findings. They also noted that the best DTIR for minimizing the uncertainty of FA was dependent on the FA value itself; in fact, a relatively high DTIR was marginally better for minimizing the uncertainty of FA, and mattered more for voxels with moderate FA (0.4-0.5) than for those with very high FA (0.8-0.9). Even so, DTIR as low as 3:1 was optimal for minimizing the uncertainty of MD, and had much more influence on MD than it did on FA. Specifically, they found that adding 7 more non-DW images to a 21-direction protocol took only 50 seconds but reduced measurement uncertainty for MD by 32%. This is in line with prior simulation studies that reported that for isotropic samples or for measuring a single ADC value, a DTIR of 3.6:1 is ideal (Eis and Hoehn-Berlage, 1995; Xing et al., 1997), but a ratio of 5.6:1 is ideal for isotropic samples (Jones et al., 1999). This led the authors of Zhu et al. (2009) to suggest that DTIR be set to 5-6, as higher DTIR values are preferable for better PEV estimation (Alexander and Barker, 2005; Kingsley, 2005). The only discrepancy between our work and that of Zhu et al. (2009) is that we found the DTIR to be 8-80 for MD, while they found a 3:1 ratio was optimal for estimating MD. One plausible reason for this discrepancy is that they were using sequences with fewer than 30 directions, whereas our HARDI sequence used 94 directions. When the number of directions is very high, it may not matter so much what the DTIR is, because the MD can be measured with good SNR. This is largely borne out by the asymptotes in the SNR plots that assessed the effect of DTIR.
Moreover, a related study by Skare et al. (2000) introduced the concept of the condition number to help in estimating upper and lower bounds for noise in tensor-valued signals. They noted that in DTI studies using a single b-value (as is the case here), the diffusion tensor is related by a transformation matrix to the measured apparent diffusion coefficients, and it is inferred via the Stejskal-Tanner equation from the diffusion-attenuated signal St in direction t (i.e., ADCi= ln(S0/St)/b). The transformation matrix depends only on the gradient directions. Using both a theoretical error propagation analysis and an empirical analysis of different direction sets, they found that the condition number of the transformation matrix was significantly correlated with the noise performances of the different data collection schemes. This is intuitively reasonable, as diffusion tensor estimation requires the inversion of a matrix that will tend to amplify noise in the gradient data when its condition number is high. In other words, when the condition number is large, the relative error in the estimated diffusion tensor elements is several times larger than the relative errors in measuring the ADC along each direction. In comparing various schemes (including tetrahedral, decahedral, tetraorthogonal, and Jones’ electrostatic repulsion arrangements), Skare et al. (2000) noted that, if other factors are equal, direction sets with smaller condition numbers should be favored where possible. They also proposed a new gradient optimization scheme (downhill simplex minimization; DSM) that optimized the condition number of gradient sets. In our current paper, the sequential optimization of energy for the nested gradient sets is intended to lead to gradient arrangements that tend to be uniformly distributed in 3D space, as far as possible given the constraints.
Eddy current corrections are sometimes necessary to correct for distortions between images with different diffusion gradient directions. Eddy currents can arise in DTI as the diffusion gradients change rapidly (Zhuang et al., 2006), and these can lead to distortion in the phase-encoding direction, slight blurring of GM/WM boundaries, and misregistration between individual diffusion-weighted (DW) images. In this study, we explicitly reduced eddy-current induced distortions in our scans by using single-shot echo planar imaging with a twice-refocused spin echo sequence (Reese et al., 2003), which was what we did in our experiment. In addition, we reported SNR from a small, relatively homogeneous region in the corpus callosum, so the effects of slight spatial shifts among the gradient images should be relatively low, at least compared with the overall spatial autocorrelation (coherence) in the diffusion signals in this relatively homogeneous region.
As in most DTI studies, and as in the other in vivo studies of SNR in DTI (Landman et al., 2007; Farrell et al., 2007; Zhu et al., 2009), we used only a single b-value of 1159 s/mm2. However, the restriction of DTI to a single b-value results in a bias toward a low bandwidth (in q-space) for the ODF. As noted by (Kuo et al., 2008), angular resolution is also a function of the q-space bandwidth. Our ongoing work using multiple b-value acquisitions at 7 Tesla (also known as diffusion spectrum imaging) has found that the ODF reconstruction error tends to be greater at higher b-values, but also falls with increasing angular resolution (Zhan et al., 2009).
To support the validity of our experiments, we also sought to confirm that the reconstruction error in resolving the diffusion profile was truly decreasing when using gradient sets with higher numbers of gradients. Figure 8 shows one example of this simulation. As expected, the sKL decreased with increasing SNR and when using more scanning directions (Figure 9). This result confirms the premise of this paper, namely that greater angular resolution leads to more accurate ODF recovery and less reconstruction error in HARDI. For these calculations, the ODF used is that derived analytically from the tensor distribution function (TDF) (Leow et al., 2009), rather than from the regularized ODF framework proposed by Descoteaux et al. (2007). Further work is necessary to fully understand how ODF accuracy depends on gradient numbers; recently, Cho et al. (2008) have examined the SNR needed to resolve fibers crossing at 45 degrees in phantoms scanned with different b-values.( Cho et al., 2008) As shown here, Cho et al. also noted that the accuracy and the angular resolution of q-ball imaging to define fiber orientations strongly depend on diffusion imaging parameters.
Optimizing diffusion imaging sequences for examining fiber integrity in the brain is a key issue in neuroscience and radiology research. Past studies based on mathematical theory and computer simulations suggested that around 18-21 gradients might suffice for accurate estimation of FA (Papadakis et al., 2000), but different studies disagree somewhat on the number of gradients needed to reliably estimate anisotropy indices. Here we improved upon past studies, by studying SNR profiles empirically in a relatively large number of normal subjects. We found that the SNR profiles and tendency to saturate depended on the DTI parameter of interest, each having a different profile of dependency on the scan time. Even so, a direct comparison of different studies is not possible, because the studies by Papadakis and Hasan used simulated data, for which noise is defined mathematically, and we used empirical data instead. There are also differences in the SNR definitions across studies, in that we included the biological variability of the signal over the region of interest, when taking into account the overall profile of variance.
The current study has three main conclusions. The number of diffusion-sensitized gradients affects the signal-to-noise ratio in different ways, for various DTI-derived parameters. For FA and RA, near-optimal SNR was achieved with 62~66 diffusion-sensitized gradients; for MD, this number was 58 and for GA and tGA, this number was about 55, while for VR, SNR increased well beyond 94 gradients. Our SNR plots, for a range of common anisotropy indices, may be useful in designing future acquisition protocols. Very high numbers of gradients were more beneficial when estimating parameters based on high-order products of tensor eigenvalues, which tend to be noisier. Differences in robustness for different DTI-derived measures derive from the arithmetic canceling or boosting of noise in the formulae for each index. The generalized fractional anisotropy (GFA), for example, required very large numbers of gradients to estimate reliably, and its SNR plot had not stabilized even with 100 gradients. Although GFA has been advocated as a good replacement for FA when large numbers of gradients are available, it requires the fitting of a spherical harmonic series expansion to the diffusion profile, which can be unstable unless very high angular sampling is available. Third, we confirmed that there is a special ratio of non-diffusion sensitized images, N(bx)/N(b0), that should also be collected to optimize SNR. Interestingly, this optimal ratio is also parameter-dependent, with values around 26 for VR, 10.57 for FA, 9.17 for RA, 9.13 for GA and tGA. This ratio was higher for some anisotropy indices (VR), lower for others (FA, RA, GA and tGA), and was not relevant at all for measuring mean diffusivity (MD), which is a relatively robust measure, which tends to cancel noise somewhat as the average of the diffusion tensor eigenvalues. Compared to other parameters, MD had relative higher noise resistance, needing fewer gradients to estimate, and fewer baseline signals.
The saturation number of each index should not be considered as a measure of the quality of each index for assessing biological characteristics of the white matter, such as fiber integrity. In our past studies, we found that FA outperformed GFA in terms of its correlation with IQ scores, even in 94-gradient scans that were expected to be sufficient to allow an accurate measure of GF (Chiang et al., 2008, 2009); we also found that genetic influences on fiber characteristics were marginally more evident with tGA than FA, perhaps partly because the geodesic anisotropy takes into account the manifold structure of the diffusion tensors, when comparing them (Lee et al., 2008). Ultimately, even if over 100 gradients are insufficient to reliably estimate GFA, it may not be any more useful for certain time-limited applications, even if it carries potentially different information not provided by other anisotropy indices.
Furthermore, there are several caveats regarding this study. Past studies relied mainly on computer simulations, which may not fully reflect the true sources of noise in a human subject’s scan. Even so, the achievable SNR in empirical data will also inevitably depend on the scanner, the magnetic field strength (here 4 Tesla), the spatial resolution, and circumstantial factors such as the amount of subject motion in the scanner. The issue of subject-specific variations was somewhat alleviated in the current study, as we examined 50 subjects’ scans, deriving plots for the minimum and maximum SNR achieved, as well as its mean profile, in a representative group of subjects. Even so, care may be needed in extrapolating our results to other scanners and field strengths. Second, if a study is designed to detect group differences in DTI (e.g., comparing patients with Alzheimer’s disease to groups of normal subjects), the image SNR is not the only relevant source of variation to consider – biological variability may be very high or relatively low, or the available sample sizes may be very high, so other considerations may supervene, and outweigh the need for very high SNR in each scan. Although our current study focuses on SNR improvement, whether or not this increased SNR translates into smaller minimal sample sizes to detect clinically relevant effects depends on the biological variation in these measures across subjects. Landman et al. (2007) also noted that the estimated orientation of the diffusion tensor (principal eigenvector) depends not just on the angular sampling, but also on patient motion, field inhomogeneity, and EPI-related distortions. Thus, further studies of scanner field strength, spatial resolution, tolerability, motion artifacts, test/re-test reliability, and clinical effect sizes are needed to evaluate the added benefit of HARDI’s SNR for radiologic and neuroscientific studies.
Finally, in considering the trade-off between SNR and scan time, alternative uses of the scan time must also be taken into account, such as acquiring other potentially interesting MRI sequences (e.g., spectroscopy or FLAIR scans), minimizing patient burden, and minimizing the risk of non-compliance with longer scan times (i.e., losing the scan altogether). High-angular resolution acquisitions are time-consuming to collect, but it must be conceded that they may provide new insight into fiber architecture and connectivity that cannot be achieved, even in principle, using smaller numbers of gradients. Tractography studies, for example, derive angular information on local fiber trajectories, and tend to benefit from higher angular resolution well beyond the sampling limits where the SNR profiles for scalar anisotropy measures are saturated. Even so, these calibration plots may be of interest in designing future DTI protocols for assessing fiber integrity in the living brain.
This work was supported by NIH grant R01 HD050735 and National Health and Medical Research Council, Australia grant 496682. Algorithm development was also supported by NIH grants EB007813, EB008281, EB008432, HD050735, and AG020098 (to P.T.).
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.