Abstract
This work investigates the possibilities of applying high-angular-resolution-diffusion-imaging- (HARDI-) based methods in a clinical setting by investigating the performance of non-Gaussian diffusion probability density function (PDF) estimation for a range of b-values and diffusion gradient direction tables. It does so at realistic SNR levels achievable in limited time on a high-performance 3T system for the whole human brain in vivo. We use both computational simulations and in vivo brain scans to quantify the angular resolution of two selected reconstruction methods: Q-ball imaging and the diffusion orientation transform. We propose a new analytical solution to the ODF derived from the DOT. Both techniques are analytical decomposition approaches that require identical acquisition and modest postprocessing times and, given the proposed modifications of the DOT, can be analyzed in a similar fashion. We find that an optimal HARDI protocol given a stringent time constraint (<10?min) combines a moderate b-value (around 2000?s/mm2) with a relatively low number of acquired directions (>48). Our findings generalize to other methods and additional improvements in MR acquisition techniques.
1. Introduction
Diffusion-weighted magnetic resonance imaging (DW-MRI) is a clinical medical imaging technique that provides a unique view on the structure of brain white matter in vivo. Right from its early stages, in the early 1990s, DW-MRI was perceived to have immediate value for the evaluation of neuropathologies such as acute ischemic stroke. Since then, numerous advents in diffusion imaging technology have greatly augmented image quality unraveling new clinical applications. Moreover, the debut of diffusion tensor imaging (DTI) and fiber tractography enabled a completely new, noninvasive view on white matter fibre bundles connecting gray matter neural populations, of increasing importance for cognitive neuroimaging applications. With DTI and fiber tractography, the understanding of several neurological and psychiatric disorders, such as schizophrenia, traumas, stroke, and edemas, has been increased, and they have also been applied clinically to aid presurgical planning before intracranial mass resections.
In DW-MRI, white matter fiber bundles are probed indirectly by measuring the directional specificity (anisotropy) of local water diffusion. During a small time interval (the “so-called” effective diffusion time ), at each location of the tissue, the diffusion causes a displacement of water molecules. Postprocessing of diffusion-weighted images is fundamentally aimed at calculating the probability density function (PDF) , for this displacement of water molecules in each imaging voxel, where is the 3D displacement vector. Since fiber bundles in tissue allow relatively free diffusion of water along the fibers but tend to obstruct the diffusion perpendicular to the fibers, one assumes that the reconstructed PDF, also known as the diffusion propagator, reflects the local fiber structure within a voxel. Diffusion tensor imaging (DTI) approximates the diffusion propagator by a 3D Gaussian function, which allows a relatively simple estimation of a small number of sampled diffusion directions [1]. However, it has recently become clear that a simple unimodal Gaussian approximation of the diffusion PDF cannot represent interesting complex fiber architecture such as single fiber bundles curving or fanning or multiple fiber bundles with different directions contained in a single imaging voxel [2–6]. Consequently, much of the modeling efforts have aimed at estimating multimodal (non-Gaussian) characterizations of the diffusion PDF that capture complex fiber architecture, particularly the presence of multiple fiber populations in single voxels. Although the full diffusion PDF provides the most extensive information about the structure of the material at microscopic scale, for many applications like fiber tracking, it is the orientation that is of most importance. Therefore, a simplified function is often derived, called the orientation distribution function (ODF), defined as the radial projection of the diffusion function. Crucially, a characterization of the ODF generally puts much lower requirements on the sampling of diffusion directions and, particularly, degrees of diffusion weighting than a characterization of the full diffusion PDF. Many recent modeling techniques, such as Q-ball [7–9], the DOT [10], and spherical deconvolution (SD) [11], focus on obtaining a low-variance unbiased estimate of the diffusion ODF (or the “fiber” ODF that represents probabilities of fiber orientations, rather than water diffusion positions) with only moderate acquisition requirements. The required acquisition is generally a high angular resolution diffusion imaging (HARDI) acquisition that acquires a large number of diffusion directions at a single level of diffusion weighting (a “moderate” q-space sampling on a single spherical shell).
The general convention for conducting HARDI acquisitions in presenting, extending, or testing non-Gaussian techniques is to sample as many gradient direction vectors, at as high as possible -values that the scanner and the subject stamina allow, while maintaining sufficiently high SNR. Tuch [4] reports -values of ?s/mm2 and number of gradients () of about 492, whereas Hagmann et al. [12] reports ?s/mm2 and , and Descoteaux [13] ?s/mm2 (although he states that ?s/mm2 is desirable) and . Compared to the modest q-space sampling for DTI (with total acquisition time of 3–6 minutes [12]), the HARDI techniques have 3-4 times longer acquisition times on average. Although there is some variation in the requirements between different techniques, generally a large number of diffusion directions (on the order of 50–300) at a considerable level of diffusion weighting (on the order of 1000–4000?s/mm2) with a high signal-to-noise ratio is needed. In addition, fiber tractography requires whole-brain coverage at a high spatial resolution (i.e., small isotropic voxel size). However, each of these requirements generally increase the acquisition time needed, potentially making HARDI unattractive for clinical applications or cognitive neuroscience investigations that want to add a limited time diffusion protocol to their fMRI measurements.
Therefore, an important combination of open questions about the measurement requirements of these HARDI modeling techniques arises. To what degree does the dense angular gradient sampling improve the accuracy of these techniques, particularly in less favorable SNR regimes? Do high -values contribute to the accuracy of the reconstructions, given the SNR penalty that needs to be payed? More generally, given the multiway trade-off between spatial resolution, angular resolution, -value, SNR, and measurement time, what are suitable settings for these parameters in a HARDI protocol for a modern 3T scanner and considerable time constraints (10–20 minutes) typical for a clinical or fMRI experiment setting? Finding the optimal parameters for the acquisition schemes in terms of both the direction sets and -value has been an open issue in HARDI acquisitions and a few investigations have partially tackled this issue. Crucially, all of these studies focussed on a situation where acquisition time constraints were very lenient (long diffusion-only protocols on experienced subjects) or absent (phantom and simulation studies). This allowed them to perform their investigation of optimal sampling schemes in a very favorable signal-to-noise ratio (SNR) regime.
This work investigates the possibilities of applying HARDI-based methods in a clinical setting by investigating the performance of non-Gaussian diffusion PDF estimation for a range of -values and diffusion gradient direction tables. It does so at realistic SNR levels achievable in limited time on a high performance 3T system for the whole human brain in vivo. We use both computational simulations and in vivo brain scans to quantify the angular resolution of two selected reconstruction methods: Q-ball imaging and the diffusion orientation transform (DOT) (and its derivations). Since the output of the DOT is the diffusion PDF, for fair comparison in this work, we derive the ODFs from the DOT and quantify the angular resolution in a similar fashion as for Q-ball. Moreover, we propose an analytical solution of the ODF derived from the DOT. Both techniques are analytical decomposition approaches that require identical acquisition and modest postprocessing times, and given the modifications of the DOT can be analyzed in a similar fashion. Our findings generalize to other methods and additional improvements in MR acquisition techniques.
2. Methods
2.1. Q-Ball Imaging
In DW-MRI, the obtained measurements do not give the diffusion propagator or its derivations, such as the ODF, directly. We review the transformation from measured signals to reconstructed directional profiles under a few of the most used modeling assumptions and conditions.
The simplest transformation that relates the measured signal and the full PDF is the 3D Fourier transform as shown by Callaghan et al. [14]: where is the normalized signal, is the unweighted or zero-weighted baseline signal obtained without any applied diffusion gradients, and is the diffusion 3D wavevector that characterizes the direction and strength of the diffusion weighting obtained from the nuclear spins in the presence of a diffusion sensitizing magnetic gradient. This relationship has been exploited in diffusion spectrum imaging (DSI) [15, 16], where the full PDF is reconstructed by a discrete Fourier transformation of a densely sampled (on a cartesian grid) signal . However, the requirements on measurement duration in DSI are considerable, and one might not be always interested in reconstructing the full and exact PDF. Therefore, depending on the application, simpler and less accurate methods are usually applied. For the purpose of simplification, many assumptions about the signal and/or the underlying PDF can be made, hence many different modeling techniques have been proposed in the literature among which are the follwing.(1)Assumption of a mono-exponential signal decay and Gaussian PDF: the simplest approach is DTI. In DTI one assumes a Gaussian diffusion propagator which results via Fourier transformation in a Gaussian signal , where is the acquisition parameter, is the effective diffusion time, and . The apparent diffusivity is modeled by a second order tensor , such that with gradient one assumes . DTI requires very modest gradient samplings on a q-shell [15] to determine its coefficients. However, this approach is limited, as it can resolve only one fiber bundle per voxel.(2)Assumption of a monoexponential signal decay : in HARDI, to increase the angular resolution, the apparent diffusivity is modeled in a less crude way by spherical harmonics (SH) or high order tensors (HOT) [3, 6, 17]. This enhances the characterization of tissue properties, but also increases the number of necessary samples, and gives incorrect results regarding the directions of the fiber crossings as shown by Özarslan et al. [10]. DOT maps the ADC into a probability function by solving the Fourier transform using the Rayleigh expansion of a plane wave in spherical coordinates. It is worth mentioning that DOT permits multi-exponential signal decay which allows the exploration of multiple q-shell samplings.(3)No assumption about signal decay or the underlying diffusion process: Q-ball imaging [4] simply maps the measurements on a q-shell via the Funk-Radon transform into another function on a sphere called orientation distribution function (ODF). This ODF is the integral in radial direction of the true diffusion propagator : where . One of the biggest recently developed advantages of Q-ball is its analytical solution of the Funk-Radon transform [8, 9, 18] utilizing a basis of spherical harmonics (SH), which significantly improves its performance.(4)Lately, spherical deconvolution (SD) methods proposed in the literature [11, 19] are becoming popular, especially for fiber tracking purposes due to sharper profiles of the reconstructed fiber orientation distribution (FOD) . One assumes, that the DWI signal can be modeled by a superposition of fiber bundle signals via the convolution . The FOD that one obtains by a spherical deconvolution of with reflects the angular distribution of fiber bundles. One has to keep in mind that this approach incorporates assumptions about the single fiber response . Recently, it has been shown that some of the bias in the fiber extraction from ODF and FOD (especially under small angles of crossings) can be improved using tensor decomposition approach as in the work of Schultz and Seidel [20].There are many other diffusion modeling techniques in the literature, but these exceed the scope of this paper. In our analysis we will include the analytical Q-ball proposed by Descoteaux et al. [8], since it is an analytical solution to the numerical ODF [15], and it has advantages as such. In our comparison, we excluded the Laplace-Beltrami regularization proposed by Descoteaux et al. [8], due to the increase of the complexity of the parameter space examined in this work as well as the difficulties that arise when trying to compare this method to the DOT (see below) that is proposed without any regularization. However, due to the assumption of the monoexponential decay for the signal attenuation in DOT, a smoothing in q-space can be assumed and has advantages for smaller -values where the monoexponential signal decay holds.
2.2. Diffusion Orientation Transform
The DOT [10] reconstructs the diffusion propagator at a given radius and a unit vector . In the DOT, the radius for reconstructing a single shell of the PDF has to be determined heuristically, and this process is not a trivial task, especially in the case of real data. is highly dependent on the scanning parameters, especially the -value, and possibly the tissue properties of the voxel that we try to resolve.
2.2.1. Numerical DOT-ODF
To avoid the selection of and inspired by definitions of the ODF from q-ball imaging [4] and the marginal ODF (mODF) from diffusion spectrum imaging (DSI) [21], we propose similar ODFs computed from the DOT as where is the PDF computed from the DOT [10], with as a unit vector defined by , where . In our calculations we truncate the integrals (summations in a discrete form) to an arbitrary high value, , and then perform a large number of summations from the full diffusion PDF. In Figure 1, we illustrate the diffusion profiles of the DOT and the proposed derivations. We observe that the DOT-mODF performs similarly as the DOT itself, and the Q-ball and DOT-ODF profiles exhibit resemblance as well.
2.2.2. Analytical DOT-ODF
Even though the sensitivity on the parameter for the numerical DOT-ODF is less than for the original formulation, an arbitrary high value for still must be specified. Additionally, a high number of integrations (or in discrete sense summations) have to be performed in each orientation of the full PDF to calculate the ODF value in that direction. Therefore, we propose the analytical derivation of the DOT-ODF explained in detail in Appendix A. The advantage of using the analytical DOT-ODF, is not only in reducing the computational expenses but also in avoiding the selection problem entirely. This way, we propose an ODF that can be easily compared to the ODFs proposed by Tuch [15] and later Descoteaux et al. [8]. Moreover, the analytical DOT-ODF method is less computationally expensive than the numerical method, since it does not require sampling of multiple shells of the PDF for summations. The required computational time is the same as for the analytical Q-ball.
2.3. Data
Previous studies have largely focussed on either software simulations or real MRI data. Where the analysis and validation of the simulated data can be absolutely quantitative, one must be cautious in generalizing results since simulated properties of tissue, acquisition, and noise parameters could be unrealistic for any simulation. In contradistinction, where full realism is evidently obtained in actual in vivo data, the analysis must be relatively qualitative because of the absence of an exact gold standard. To quantify the accuracy of the selected modeling techniques, we generate both synthetic data of two fiber crossings, and we extend our analysis to real data acquired in vivo on a human subject. In this way, we combine the advantage of a ground truth and full knowledge and control of all parameters in the synthetic data simulations, with the guaranteed realism of the tissue, acquisition, and noise parameters in the real data.
2.3.1. Synthetic Data
We generate synthetic software data according to the Söderman and Jönsson’s [22] model used in Özarslan et al. [10], Brampoutis et al. [23], and von dem Hagen and Henkelman [24].
We fix the parameters for this model similarly to Özarslan et al. [10] to fiber length , fiber radius , and free diffusion coefficient .
We calculate the effective diffusion time as where the correction is due to the diffusion which occurs during the time in which the gradients are on.
Using Söderman’s model we create two fiber crossings under angles of: 40°, 45°, 50°, 55°, 60°, and 90°. The choice for these angles was made given the expected accuracy for Q-ball imaging [8, 25] of around 50° (however, this is strongly dependent on the -value and applied regularization). To create different datasets with different acquisition parameters we vary the number of gradients and -values as in Tables 1 and 2, respectively.
Additional parameters that need to be specified for the Söderman’s model are the values for the gradient pulse duration and gradient spacing . We use the same values as from our real acquisition protocol where we optimize the echo time per -value.
This way we optimize the SNR per -value, and the resulting profiles will be less deteriorated from noise. However, one has to be careful here because the representations of the PDFs are for different measurement times. For all the data in this work, the same parameters listed in Table 2 are used in order to allow for a fair relation between the synthetic and in vivo data. To create noisy synthetic data we add Gaussian noise to the real and complex part of the signal, and we vary the standard deviation depending on the analysis we want to perform.
2.3.2. In Vivo Human Data
DW-MRI acquisitions were performed on a female subject using a twice refocused spin-echo, echo-planar imaging sequence on a Siemens Allegra 3T scanner (Siemens, Erlangen, Germany) equipped with a 40?mT/m head gradient set, an 8-channel cylindrical phased-array RF-coil. Informed consent was obtained prior to the measurement. Twenty axial slices (FOV?mm2, matrix ?mm3 voxel size) were positioned through the body of the corpus callosum and the centrum semiovale. EPI acquisition was performed with 3/4 partial Fourier and a readout bandwidth of 2840?Hz/Px. Datasets were acquired with diffusion gradient sampling schemes of 24, 36, 48, 60, and 72 unique directions obtained by a static repulsion algorithm, with the diffusion-weighted volumes interleaved with zero-weighted volumes every 12th scanned diffusion gradient direction. For each gradient scheme the scans were performed with -values of 1000, 1500, 2000, and ?s/mm2 (, resp., 77, 85, 91, and 100?ms, ?ms throughout), using the same gradient pulse duration and gradient spacing as given in Table 2. No averaging of repeated acquisition was performed, that is, there was only a single acquisition for each parameter set.
2.4. Analysis of Synthetic Data
To analyze the accuracy of the reconstructed HARDI models in the synthetic data sets, we report the detected angular error. We use a simple scheme for determining the angular error, that is, we find the local maxima of the reconstructed profiles on a mesh generated by tessellating icosahedron of certain order and then report the angular difference between these maxima and the simulated (true) fiber directions. To reduce the error that the discrete mesh imposes, in our analysis we use the 7th order of tessellation of the icosahedron, that imposes 1° spacing between tessellation points, thus, a maximum distance between any point and the nearest tessellation point of 0.5°. Depending on the orientation of the fibers in the space of the gradients, different errors might appear for the same simulated angle. For example, if we simulate an angle of 60° by two crossing fibers and then rotate the fibers in 3D space, one would assume that always the same angle of crossing is detected. However, in Söderman’s model, different angular errors are calculated for different orientations of the fibers. In general, these errors cause less than 5° angular error, and therefore while choosing our optimal parameters for certain angular configuration, we take these errors into account. For example, for an angle of 60°, the -value of 3000?s/mm2 might give an angular error of 4° and the -value of 4000?s/mm2 an angular error of 2.5°. However, due to the errors imposed by the tessellation of the mesh and the orientation of the fibers as well as our choice of tolerance for 5°, we report the b-value of 3000?s/mm2 as preferable, since both -values are within our range of tolerance, and we always report the smallest -value as optimal.
2.5. Analysis of Human Data
2.5.1. Criteria for In Vivo Data Analysis
For qualitative analysis of the real data, we select the centrum semiovale (CS) (see Figure 2) where crossings are to be expected. This is a challenging region for DW-MRI analysis techniques, since fibers of the corpus callosum (CC), corona radiata (CR), and superior longitudinal fasciculus (SLF) form a three-fold crossing. A region-of-interest (ROI) was defined on a coronal slice. Red rectangle indicates the position of the slices in the acquisition.
The qualitative criteria for comparison are the sharpness of the profiles in CC and lower parts of the CR, the presence of two ODF maxima for the CC and CR and the spread of the their crossing in the centrum semiovale, and the presence of the superior longitudinal fasciculus in the upper corners of the region. Furthermore, the presence of spurious peaks is an indication of poor SNR and unreliably detected crossings.
2.5.2. SNR in In Vivo Data
We calculate the SNR in the images via the “difference method” as described by Reeder, et al. [26]. For a set of two unweighted volumes, the SNR is computed via the difference method as where and are two unweighted images and and are the mean and the standard deviation respectively. The SNRs, computed in an ROI of 60 voxels in the middle of the corpus callosum, is 23.58, 20.75, 19.29, and 18.61 at -values of 1000, 1500, 2000, and 3000?s/mm2, respectively. The final reported SNR is calculated as the median of six SNRs, each computed from two unweighted images closest in time. Since unweighted volumes were interspersed between 12 diffusion weighted volumes, the SNRs are probably slightly under-estimated due to bulk motion and drift over the measurement time. Note that, due to the characteristics of the phased-array RF coil, the SNR in the CC in the middle of the head is a lower bound for the SNR in the rest of the scanned volume. SNR is expected to increase (quadratically) towards the periphery.
2.5.3. Scanning Time
To give an indication for the time requirements for the individual acquisition schemes, we estimate the total scanning time needed to obtain a full brain acquisition at the given parameters. Assuming at least 60 axial 2.5?mm slices would be needed for whole brain coverage and the shortest achievable TR would be used, whole brain measurement times are given in Table 3.
3. Results
In this section we will discuss the quantitative results from our synthetic data and give qualitative implications for the real data scanned under the same parameters.
3.1. Synthetic Data Analysis
To get an indication on the optimal parameters for each simulated angle, we run a series of tests, where we determine the minimal angular error that can be obtained with the minimal combination of -value and -order, under different gradient sampling schemes. The results are presented in Table 4.
We observe similar trends for both of the reconstruction techniques. Both Q-ball and DOT-ODF can accurately represent crossings of 90° for -values as low as 1000?s/mm2 and order of SH regardless of the choice for the gradient sampling scheme. The scenario changes for the angle of 60°, where the best choice for -value recommends 3000?s/mm2 and higher order of SH of . Trying to resolve smaller angles of crossings, with these analytical decomposition methods results in higher choices for -values and SH order even in noiseless cases. The SNRs in clinical scanners at high -values such as 4000?s/mm2 drop and therefore these -values are not recommended. The noise corrupts the reconstructed profiles, and modeling with high orders of SH will result in capturing the noise, which is dominant in these cases. Therefore, from the noise tables (see Supplementary Material available online at http://dx.doi.org/10.1155/2013/658583) we already get the indication that the lowest angle which can be captured within reasonably imposed angular error would be around 60°. To reduce the parameter space, we therefore select this interesting angle and SH order of (as Table 4 indicates for this angle) and look at some interesting influences of the -values and number of gradients, under different SNRs. Another motivation for this choice of reconstruction angle is the region that we select for illustrating the in vivo results. The centrum semiovale contains a large amount of voxels where the crossing angles are approximately between 60° and 90°. To investigate the effect of the gradient sampling schemes, we performed tests where for -values of 2000?s/mm2, we vary the number of gradients under different levels of added Rician noise. In Figure 3, we report the calculated angular error and the standard deviation of 100 noise realizations for simulations of 60° angle of crossing. We observe that for SNR of about 20 (similar to our case) of the recovered angle suffers from an angular error of less than 10°. For higher SNRs (20–40), we observe that around a number of 70 gradients the graph stabilizes and we do not gain much accuracy by having denser gradient sampling schemes. The results are similar to the ones from the DOT-ODF and for completeness are reported in Figures 8, 9 and 10.
In a similar fashion, we investigate the effect of the -value, on the accuracy of the reconstruction techniques under different SNRs. From the previous discussion, we concluded that the number of 70 gradient direction is sufficient for good detection of a 60° crossing fiber configuration, and therefore in this analysis we fix the number of gradients to 72 (Figure 4). For data with SNR 10, there is a clear tendency for an optimum -value of around 2000?s/mm2 as seen in the decrease in the angular error. Higher -values, produce worse results with high standard deviations. However, if the SNRs are higher, the angular error and the standard deviation decay and stabilize around -values between 2000?s/mm2 and 3000?s/mm2. Again, the results coincide with the ones from the analysis of analytical DOT-ODF, and for completeness are reported in Figures 8, 9 and 10.
The previous analysis emphasizes the importance of the SNRs for reliability in the reconstruction of analytical HARDI decomposition techniques. To address this important issue, we investigate the effect of the SNR under fixed parameters for -value and number of gradients and report in Figure 5. We observe that data with SNRs lower than 15 would produce very unreliable datasets, where the accuracy would drop to an angular error bigger than 10°. Smaller angles cannot be resolved accurately and most probably will not be detected as distinct fiber populations. It can also be seen that the effect of the number of gradient directions (at least at the level between 48 and 72) has much less of an effect on detection performance than the SNR level does. Take note that this is a simplified model of the diffusion process, where only two fiber populations are simulated. The scenario gets more complicated when more fiber populations interfere within a voxel. Increasing the number of gradients does not significantly improve the angular error in the case of SNR of about 20.
3.2. Real Data Analysis
For better appreciation of anatomical locations, we present the results by visualizing the recovered maxima, overlaid over an FA-map. To be able to observe the changes in the shape of the glyph profiles, in Figures 11, 12, 13, 14, 15, and 16 we additionally report the min, max normalized glyphs since the normalization enhances the profiles. However, in the noisy areas this considerably deforms the glyphs, and should be interpreted with care. For the color-coding we use standard RGB directional coloring. We present the results from Q-ball in all our figures since they are almost identical with the ones from the DOT-ODF (see Figures 1(a) and 1(e)).
Figure 6 illustrates the effect of the -order acquired under different parameters. The visualization is by recovered maxima. We observe that in all of the cases, the threefold crossing is better recovered at order 6, at the higher -value (3000?s/mm2), and higher number of gradients . At these parameter settings the structure of the SLF (yellow rectangle) becomes more evident. Therefore, we conclude that SH order is the optimal for recovering the crossing information in this region. This statement corroborates our findings from the synthetic data. Going to orders higher than results in many false maxima detected from the spurious peaks coming from the noise. Additionally, we observe that = ?s/mm2 performs better than ?s/mm2, since more threefold crossings are detected in the centrum semiovale.
The key investigation that answers the question of clinically realizable HARDI protocols is of the effect the -value (and thus the SNR level) and the number of gradients on the accuracy of the analytical decomposition techniques like Q-ball and DOT (and its previously defined derivations DOT-ODF and DOT-mODF). From the previous synthetic analysis, we observed that at about a number of 70 gradients the angular error stabilizes and does not decrease significantly when increasing the number of gradients. However, it is interesting to investigate how low the number of gradients can go at a given -value and SNR level, such that the data still has reasonable reconstructions and detected crossings. Therefore, we select gradient schemes with lower than 72 to investigate this question.
In Figure 7, we observe this influence of the -value and (for different -orders and glyph representations refer to Figures 11, 12, and 13). We observe that the crossings from CC and CR are already observed in the lowest parameter combination of ?s/mm2 and (for order already at ). However, the SLF structures are not that obvious and become more prominent at higher ?s/mm2 (see yellow marked rectangles in Figure 7). However at higher -values we observe many spurious peaks that disappear with the increase of the number of gradients (this can be better observed in Figure 11). At there is no appreciable difference between the profiles at ?s/mm2 and ?s/mm2. Given that the spurious peaks are slightly more prominent at ?s/mm2, and that lower -values are more attractive with respect to the total scanning time, we conclude that ?s/mm2 and would be preferable with a total scanning time of about 7–10 minutes. For datasets with SNR of about 20, the accuracy of the reconstructed profiles at angles of about 60° would be with angular error of about 8°.
4. Discussion
The benefits from applying HARDI modeling techniques that require modest q-space single shell sampling over the classical DTI model are well known. However, once the potentials are clear from the theoretical point of view, the questions about the employment of these techniques in every day practice (clinical or cognitive neuroscience research) arise. In the case of HARDI, not only the modeling techniques are complex, but there are many parameters from the acquisition to the postprocessing of the data that influence the accuracy in the estimation of the underlying diffusion profile. The main bottleneck for applying HARDI in clinical practice or as an added protocol in fMRI investigations is the acquisition time, given that decreasing it generally decreases the signal-to-noise ratio in the data thus inferring the accuracy in the reconstructions. The main parameters that influence the total acquisition time are the number of gradients, the spatial resolution (especially the number of slices), and the number of averages, many of which inevitably influence the SNR of the acquired data (please refer to Appendix B). In this study, we compared only analytical decomposition techniques such as Q-ball and DOT. Recent techniques with single q-space shell requirements that are known to have better angular resolution than Q-ball such as the solid angle consideration, Q-ball [27] and spherical deconvolution techniques [11, 28] could be used as well. However, we would expect the relative conclusion regarding empirically optimal short-time protocols to also hold for these techniques, possibly at a higher overall accuracy level.
Finding the optimal parameters for the acquisition schemes in terms of both the direction sets and -value has been an open issue in HARDI acquisitions and a few investigations have partially tackled this issue. Jones [29] investigated the effect of gradient sampling schemes on the correctness of estimating the diffusion tensor (DT) and scalar indices derived from it. However, this study was limited to DTI where the effect of the -value is not as significant as in HARDI-based modeling techniques, and therefore -values were not investigated. Moreover, it is based purely on computer simulated data. Alexander and Barker [30] investigated the effect of the -value on a correct estimation of fiber orientations. This Monte Carlo study was limited only to functions that are a mixture of Gaussian distributions and a reconstruction of algorithms that fit this same class or models to synthetic data. There have been a few recent studies that investigated the effect of the number of gradient direction on real experimental data. Zhan et al. [31] focused on a user study for investigating the effect of the number of gradients on the signal-to-noise ratio (SNR) of DTI derived scalar indices. Tournier et al. [32] focused on finding the minimum required number of gradients for successfully reconstructing HARDI diffusion profiles. Other studies have used a hardware diffusion phantom to investigate the effect of measurement parameters on the accuracy of directional diffusion estimation. Cho et al. [33] concentrated on finding the optimal -value for resolving small angles of crossing (in this case 45°) using Q-ball imaging. In Poupon et al. [34] new hardware diffusion phantoms are proposed to similarly investigate the implications on the optimal -value for reconstructing 45° of crossing with Q-ball imaging. Tournier et al. [35] report an extensive study on hardware phantoms comparing the angular resolutions of Q-ball Imaging, Spherical Deconvolution [36], and constrained spherical deconvolution (CSD) [11]. However, in this study the effect of gradient sampling was not examined. Crucially, all of these studies focussed on a situation where acquisition time constraints were very lenient (long diffusion-only protocols on experienced subjects) or absent (phantom and simulation studies). This allowed them to perform their investigation of optimal sampling schemes in a very favorable SNR regime.
In this work, we investigated the relations between above mentioned parameters. We addressed some important questions that relate the SNR of the acquired data with the accuracy of the reconstructions of some popular HARDI analytical decomposition techniques: Q-ball and DOT. To make both of the techniques comparable, we derived from the DOT two numerical ODFs: the DOT-ODF (comparable to the ODF derived from Q-ball) and the DOT-mODF. To utilize the performance of the DOT-ODF, we derive an analytical solution that is comparable to the speed of the DOT and Q-ball. Finally, with synthetic data simulation, and series of real data acquisitions, we gave implications that recommend a set of optimal acquisition parameters for these techniques. If the data is represented by order of spherical harmonics, an empirically optimal combination (given stringent time constraints) of -value and number of gradients that would resolve most of the angles 60° with error of less than 10° is ?s/mm2 and . However, even at a more modest combination of the acquisition parameters (?s/mm2 and ) most of the two fiber populations of larger angles will be resolved under the cost of larger angular error. Moreover, these parameters are recommended, given that the data has a minimal SNR (in the 0-weighted acquisitions) of about 20. Decreasing the SNR (increasing the spatial resolution or increasing acquisition speed by techniques such as parallel imaging) would produce unreliable results with angular errors of more than 30°.
It is important to remark that these conclusions to a very large degree hold even for radically different MRI hardware (e.g., high performance head gradient sets, phased-array RF coils with =32 channels, or even higher main field strengths than 3T) or new acquisition developments (e.g., highly accelerated parallel imaging) when two observations are taken into account. First, many techniques that are able to speed up the acquisition of a given protocol (e.g., parallel imaging) have an inevitable cost in the SNR level (see Appendix B). Therefore, employing them does not change the fundamental trade-off between spatial resolution, acquisition time, and SNR that define the playing field for the empirical optimization presented here. Second, where the available acquisition technology allows a fundamentally higher SNR or shorter acquisition times (e.g. phased-array coils, simultaneous multi-slice or multi-band acquisitions, and higher field strengths), it stands to reason that this margin would be used to increase the spatial resolution of the acquisition beyond the 2.5?mm or 2?mm isotropic voxel level, given that a satisfactory angular resolution and total acquisition time has been obtained. Particularly, when the final processing aim of a study is fiber tractography, there is much to gain in the veridicality of tracked fiber bundles by decreasing partial volume effects and increasing the ability to resolve finer white matter structures (see, e.g., the work of Roebroeck et al. [37]).
Appendices
A. Derivation of Analytical DOT-ODF
A.1. Diffusion Orientation Transform
The diffusion probability function is related to the measured MR diffusion signal by a Fourier integral as noted in (1). In the work of Özarslan et al. [10], the point-wise expansion of the plane wave in spherical coordinates is used in order to obtain the following expression: where where is the th order spherical Besel function. Normally, the radius is fixed, but one needs to find an optimal radius in order to have a sharp enough profile. In order to avoid this problem, we consider ODF, described in the following section.
A.2. ODF from DOT
In order to obtain the ODF, we need to integrate (2): According to the work of Özarslan et al. [10], , given by (A.2), can be computed analytically where is the confluent hypergeometric function of the first kind.
In order to compute , we use the integral representation of , and change the integration order. As a result, we obtain an analytical formula for the integral where The values for up to 8 are summarized in Table 5.
Therefore, if is given in spherical harmonic representation (A.3) and (A.6) give us Since tensor terms of order belong to the span of the spherical harmonics of the same order, we conclude that the same holds for the tensor representation where -coefficients of tensor expansion for .
B. The Basics of Signal and Noise in Pulse Gradient Spin echo–echo Planar Imaging (PGSE–EPI) Diffusion MR Sequences
Leaving static magnetic field B0 outside of consideration, the signal for a spin echo experiment is directly proportional to the available transverse magnetization, which is itself a function of , , , and : To maximize the available signal, must be maximized and minimized. Given the usual requirement in diffusion MR imaging for whole brain coverage using a 2D (slice selective) sequence, the is naturally long, usually longer than 4 times (which is about 800?ms for white matter at 3T), recovering at least 98% of the available signal. In contrast, the achievable in a Stejskal-Tanner PGSE sequences is typically long compared to T2 (which is about 84?ms for white matter at 3T), due to the diffusion gradient lobes. The amount of diffusion weighting, when gradient ramps are a small proportion of the diffusion-gradient lobes, is given by the -value: where is the 1H gyromagnetic ratio and , , and are the amplitude, time between onsets, and length of the two diffusion gradients, respectively. Since the -value has a cubic dependence on , quadratic on , and linear on , the most SNR-efficient way (i.e., minimum TE) to achieve any given -value is to play out diffusion gradients at maximum available amplitude and keep the as short as possible. Therefore, the SNR in diffusion imaging necessarily falls with increasing -value, through increasing TE. An estimate of the relative SNR for different -values can be made if one estimates a minimum achievable as (assuming , and 25?ms total for RF pulse lengths, pre-TE readout train length, and additional encoding and spoiling). The SNR relative to that for a equal to (assumed to be 70?ms) for different -values and two maximum achievable gradient strengths is given in the Table 6.
In addition, the SNR for each voxel is crucially dependent on the parameters of the imaging gradients and the readout train in EPI. The SNR in 2D Fourier-encoded MRI is proportional to with = , and are voxel sizes along readout and phase-encoded direction, respectively, is slice thickness, is the number of averaged repetitions of the entire sequence, and are the number of samples along readout and phase-encoded direction, respectively, is the readout sampling time, and is the readout bandwidth. This highlights the trade off between SNR, spatial resolution, spatial distortion, and dropout and total imaging time. In particular, for PGSE-EPI sequences, SNR can be traded off in the following ways.(i)Increasing voxel size increases the SNR at the cost of spatial resolution and increased susceptibility related distortions and dropout at air-tissue boundaries. Increasing voxel size from ?mm3 to ?mm3 or ?mm3 increases SNR with a factor 1.95 and 3.37, respectively (if other parameters, particularly and are held constant). (ii)Increasing the number of averaged repetitions of the entire sequence increases SNR at the cost of increased imaging time. Acquiring 2, 3, or 4 averages increases SNR with a factor 1.41, 1.73, and 2, respectively (if the patient does not move between acquisitions). (iii)Decreasing readout bandwidth increases the SNR at the cost of spatial resolution (due to -decay-related blurring) and increased susceptibility related distortions and dropout. In principle, however, the SNR gain can be partially compensated by increased with a longer readout train. (iv)Increasing the number of phase-encoding steps (equivalently: the number of readouts) increases the SNR when all else is held constant. However, this increases the readout train length and , which can compensate the SNR gain. In practice, is almost routinely decreased, by employing either partial Fourier sampling or accelerated (parallel) imaging along the phase encoded direction (or both), allowing a decrease in . The shorter readout trains also decrease susceptibility related distortions along the phase-encoded direction. In addition to the factor , parallel imaging further decreases the SNR by a spatially varying factor which depends on the coil array geometry and acceleration factor used.Finally, SNR crucially depends on the type and performance of the employed coil array for RF signal reception. A phased-array coil design can give considerable SNR advantages over a single circularly polarized birdcage design. However, the signal sensitivity for a phased array decreases quadratically with distance from the coil array, towards the center of the head, thus showing considerable spatial inhomogeneity of the SNR. A phased array of 8 coils can have a typical SNR advantage of 100% (6?dB) at the peripheral levels of the brain (gray matter) compared to a birdcage coil of similar dimensions but only 10% more signal-to-noise in the central white matter. Finally, a phased array coil will enable parallel imaging, with the additional SNR trade-off discussed above.
Acknowledgment
This work is part of the research project Multi-Field Medical Visualization (MFMV), Project no. 643.100.503, which is financed by The Netherlands Organisation for Scientific Research (NWO).
Supplementary Materials
Table 1: Results per b-value for different SNR, using Q-ball, Södermans model, with Ngrad=97, angle=60v and l=6. fix is the mean angular error,sis the standard deviation of the angular error and p is the fraction of samples for which we could retrieve two fiber populations. Table 2: Results per b-value for different SNR, using analytical DOT ODF, Södermans model, with Ngrad=97, angle=60v and l=6. fix is the mean angular error, sis the standard deviation of the angular error and p is the fraction of samples for which we could retrieve two fiber populations. Table 3: Results per b-value for different SNR, using DOT, Södermans model, with Ngrad=97, angle=60vand l=6. fix is the mean angular error s is the standard deviation of the angular error and p is the fraction of samples for which we could retrieve two fiber populations.