Although attractive in its simplicity, the diffusion tensor model has been shown to be inadequate in the many regions of the brain that contain so-called “crossing fibers” (62
), whereby two or more differently oriented fiber bundles are co-located within the same voxel. The term “crossing fibers” is itself somewhat misleading, as it includes any situation where multiple fiber orientations contribute to the signal measured for the same imaging voxel. Therefore, this also applies to configurations that may not initially have been thought of as “crossing fibers,” for example, fiber bundles that “brush” past each other within the same imaging voxel, or even curving or “fanning” fibers (). Crossing fibers are endemic to DWI, due to its coarse resolution (~2 to 3mm) compared with the white matter structures of interest even the pyramidal tracts are only ~3-mm thick in subcortical regions (96
)]. Indeed, recent studies have shown that a significant proportion of the white matter contain crossing fibers (97
), with the most recent estimating that multiple fiber orientations can be detected in over 90% of white matter voxels (98
FIG. 12 Simulated configurations of complex fiber bundle architecture at the length scale of a single voxel (61). Note that any tract organization different from a single straight fiber population is typically referred to as “crossing fibers,” (more ...)
FIG. 13 The number of distinct fiber orientations detected within each voxel (98), overlaid on the corresponding anatomical T1-weighted image. Each voxel within the mask is colored according to the number of orientations detected (red: single orientation; green: (more ...)
These effects have an obvious impact on the diffusion tensor and any measures derived from it (101
). As the mean ADC is largely unaffected, anisotropy measures such as FA (102
) are particularly sensitive to the presence of crossing fibers (101
), as are the axial and radial diffusivities (86
) (). This has important consequences for their interpretation, as they are commonly regarded as surrogate markers of white matter “integrity.” Given the extent and profound impact of crossing fiber regions, such interpretations should only be made with extreme caution.
Crossing fibers are even more problematic for tensor-based tractography methods (62
): if one corrupt orientation estimate is encountered, the tracking algorithm may venture off course into an adjacent white matter structure, leading to both false-positive and false-negatives connections (97
). Moreover, the problem is far greater than might initially be expected: any given white matter tract of interest will traverse a large number of voxels, any of which might contain crossing fibers. It can readily be appreciated that the proportion of tracts traversing at least one affected voxel must be much greater than the proportion of affected voxels. If as much as 90% of white matter voxels are affected (98
), it is unlikely that any tracts will remain unaffected throughout their entire course.
In practice, the orientation produced by the diffusion tensor model is likely to be fairly close to the largest contributing fiber direction; this is why tensor-based tractography algorithms can often follow the corticospinal tract through known crossing fibers regions such as the pons and centrum semiovale (105
). The impact of these effects is, therefore, most severe for nondominant tracts. For example, the lateral projections of the corticospinal tract cross through regions where the superior longitudinal fasciculus is dominant. The commissural projections of the corpus callosum cross through the more dominant fibers of the corona radiata. These and other pathways, such as the acoustic radiations, cannot be delineated reliably without the ability to resolve crossing fibers (97
). Moreover, pathology can turn a normally dominant tract into a nondominant tract, as has been shown with Wallerian degeneration in the corticospinal tract (104
For these reasons, there is increasing interest in using higher order models to capture the information provided by DWI more fully. An overview of these models is provided here.
Methods based on q
-space provide an estimate of the spin propagator (or at least its angular dependence) by exploiting its Fourier relationship with the DW signal measured as a function of the q
-vector (related to the direction and intensity of the DW gradient pulse) (106
). The spin propagator corresponds to the spin displacement probability density function—in others words, it provides the probability that a randomly chosen water molecule within the volume of interest (i.e., a single voxel) will have a particular displacement over the diffusion time. As water molecules are more likely to move along fiber orientations, the spin propagator will have higher probability along these orientations. As there is no explicit need for a model of diffusion in white matter, these methods are often considered to be “model-free.”
A common criticism of practical implementations of q
-space based methods is the violation of the Narrow Pulse Approximation. The q
-space formalism is only strictly valid if the spins are approximately static during the application of each DW gradient pulse. For the in vivo case, this requires DW gradient pulse durations of the order of 1 ms or less (107
). Unfortunately, due to the limited gradient amplitudes available on current clinical systems, the required diffusion weighting cannot be obtained with such short DW pulse durations. However, it has been shown that in the presence of longer DW gradient pulses, the spin displacements obtained reflect the difference between the spin’s time-averaged positions during each DW gradient pulse (108
). In a restricted environment, this will cause an underestimation of quantitative measurements of displacement (109
), but importantly will not necessarily affect the estimated orientations (91
). In fact, recent studies suggest that using long DW gradient pulses may be beneficial for fiber orientation estimation by exaggerating the orientation dependence of the DW signal (110
Another concern with the q
-space approach is that the spin propagator does not necessarily reflect the true fiber orientation distribution (112
). Fiber orientations are typically extracted from the spin propagator by identifying the directions along which the probability of displacement is highest (91
). Unfortunately, due to the nature of diffusion, these high probability regions will be relatively broad and overlap significantly. While not necessarily a problem in itself, it does have consequences for the estimated orientations when using peaks in the spin propagator: (i) closely aligned fiber orientations will be “blurred” together and will thus be identified as a single orientation; and (ii) this can lead to a bias in the estimated fiber orientations (114
). These issues can be addressed by the introduction of a suitable model for diffusion in white matter (112
), although these methods could then no longer be called “model-free.”
With the exception of diffusion spectrum imaging (DSI), most q-space methods are based on the shorter and arguably more efficient high-angular resolution DW imaging (HARDI) acquisition (see “Data Acquisition” below). However, as data are collected on a spherical shell in q-space, it is not possible to perform directly the 3D Fourier transform required for q-space analysis. The problem is made tractable by assuming a particular functional form for the radial dependence of either the DW signal or of the spin propagator. The various methods based on both HARDI and q-space differ in the specific assumptions they make in this respect.
Diffusion Spectrum Imaging
DSI is the direct application of q
-space in 3D (91
). It requires data to be acquired on a 3D Cartesian grid in q
-space, from which it is trivial to perform the required 3D Fourier transform. Fiber orientations are identified by reducing the 3D spin propagator to its 2D radial projection, the diffusion orientation density function (ODF), and finding the peaks of this function. The largest obstacle to its routine use is the large amount of data required and the correspondingly lengthy acquisition (see “Data Acquisition” below).
Since the introduction of DSI, a number of variations have been proposed. For instance, the 3D Cartesian q
-space framework has been extended to spherical coordinates, to allow direct reconstruction from data sampled on multiple shells in q
-space (i.e., multiple q
-values per direction) (117
). Another approach, generalized DTI parameterizes the spin propagator using higher order tensors, also allowing the use of a multishell acquisition (118
). These variations may offer advantages over DSI, in terms of a simplified acquisition protocol and improved robustness of reconstruction.
Q-ball imaging (QBI) provides an estimate of the diffusion ODF using the significantly shorter HARDI acquisition (89
). It can be shown that the Funk-Radon transform of the DW signal (i.e., the integral over a great circle in q
-space) provides an approximation to the radial integral of the spin propagator (i.e., the diffusion ODF). As with DSI, fiber orientations can be extracted from QBI by finding the peaks in the diffusion ODF, and used to track through crossing fiber regions (121
). QBI has been shown to be both fast and capable of producing results similar to DSI with substantially reduced acquisition times. However, to ensure adequate SNR in the DW images, QBI is typically performed using low to intermediate q
-values; this will introduce significant blurring into the ODF, since the approximation inherent in QBI isvalid in the limit of large q
-values. Although blurring can be reduced using large q
-values, this can only be done at the expense of SNR and/or scan time.
Persistent Angular Structure MRI
Persistent angular structure MRI (PAS-MRI) provides an estimate of the so-called PAS of the spin propagator using HARDI data (94
). The underlying principle is that spins are assumed to diffuse by a fixed distance r, with an angular distribution given by the PAS. With this definition of the radial dependence of the spin propagator, it becomes possible to perform the 3D Fourier transform required for q
-space analysis. PAS-MRI is also combined with a maximum entropy constraint to improve the stability of the results. Unfortunately, the current implementation of PAS-MRI is computationally intensive, limiting its practical use. On the other hand, the entropy constraint allows PAS-MRI to operate on low b
-value data (e.g., b
= 1156 s/mm2
in Ref. 94
The Diffusion Orientation Transform
The diffusion orientation transform also operates on HARDI data and provides an estimate of the spin propagator evaluated at any given radius R0
). The 3D Fourier transform is made tractable by assuming a mono-exponential radial dependence for the DW signal. The diffusion orientation transform differs from most other q
-space methods in that the ODF provided is not a radial projection of the spin propagator, but corresponds to the amplitude of the spin propagator for a chosen displacement R0
. This has the advantage of providing increased separation between the various fiber orientations when using larger values of R0
Methods based on a mixture model rely on explicit models to provide an estimate of the DW signal arising from each distinct fiber population. The DW signal that would be measured for a particular combination of fiber orientations is assumed to be given as the weighted sum of each population’s contribution to the DW signal. Estimating the fiber orientations given the data then becomes a matter of fitting the model to the measured DW signal.
Mixture models assume that the DW signal measured from a given voxel is simply the linear sum of the DW signals for each fiber population present within that voxel. This allows the problem to be expressed as a linear combination, which can be solved using relatively simple and established methods. This assumption is met as long as there is negligible exchange of water molecules between the various fiber populations during the diffusion time (i.e., ~50 ms). Note that this is distinct from exchange between different cellular compartments (e.g., intracellular/extracellular). In this context, we are concerned with exchange between much larger scale bundles of fibers.
This assumption can be justified by the following arguments. First, the typical mean spin displacement across white matter fibers is of the order of 10 mm or less (85
). Significant exchange can, therefore, only occur over the ~10 mm immediately adjacent to the interface between fiber bundles. Therefore, exchange effects can only become significant if fibers from the different bundles interdigitate at the micron scale. Second, there is increasing evidence that water molecules within the intra-axonal space are to a good approximation restricted (85
), in which case there will be negligible exchange between adjacent axons, and by extension between different fiber bundles. Although the extra-axonal water may exchange, its impact will remain small as it occupies only ~20% by volume and is, moreover, generally assumed to diffuse more rapidly (and hence will be more strongly attenuated) than the intra-axonal water.
Mixture model methods also typically assume that fiber bundles in the brain share at least some, if not all of their DWI characteristics, which make it possible to reduce the complexity of the model and increase the stability of the reconstruction. This can range from simply assuming axial symmetry of the diffusion signal about the fiber axis (i.e., the diffusion tensor is prolate for each bundle) to stating that the diffusion signal is effectively identical for all fiber bundles (i.e., the diffusion tensor has fixed eigenvalues for each bundles). The latter assumption is commonly made as it leads to much better conditioned reconstructions and allows the problem to be expressed as a spherical deconvolution (see below). By extension, this also implies that any observed variations in white matter anisotropy (as measured by DTI) are due entirely to crossing fiber effects.
Although this assumption may seem naïve given the known variations in axonal diameters, packing density, and myelination levels, amongst other potential confounding factors (e.g., 125
), there are reasons to suggest that it is justified, at least in the context of fiber orientation estimation. First, these parameters have a relatively weak effect on the anisotropy: axonal membranes are sufficient to drive anisotropic diffusion (85
). Second, minor changes in parallel diffusivity can only weakly affect the axial DW signal as it is already strongly attenuated, while changes in perpendicular diffusivity are almost indistinguishable from changes in the volume fraction of the corresponding bundle (92
). In both cases, the estimated orientation will not be affected.
An additional advantage of these methods is the possibility of introducing constraints based on prior knowledge about the fiber orientation distribution. In particular, the constraint of positive (or at least non-negative) volume fractions is commonly included, either implicitly or explicitly (85
). A maximum entropy constraint can also be used, which favors a distribution of fiber orientations with few well-defined peaks (129
). Such constraints, where they can be applied, help to improve the conditioning of the problem, and can provide results that are much more robust to noise.
Multitensor fitting extends the diffusion tensor model to handle multiple fiber orientations and typically makes use of HARDI data. In this model, the DW signal is assumed to originate from a mixture of compartments, each characterized by its own diffusion tensor. To improve the stability of the method, the shape of each diffusion tensor is assumed to be prolate and axially symmetric. Its anisotropy is also typically fixed (90
), although some implementations allow the anisotropy to vary (127
). An additional isotropic compartment is sometimes included to account for CSF or gray matter contamination (97
). Importantly, these methods require an estimate of the number of fiber orientations to include in the model for each voxel, which is typically achieved by some form of model comparison. Some implementations are framed as Bayesian inference problems, providing the posterior distribution of fiber orientations given the noisy data. Therefore, they form a natural basis for probabilistic tractography methods that rely on the availability of such a distribution (97
Combined Hindered and Restricted Model of Diffusion
Combined hindered and restricted model of diffusion (CHARMED) is strongly related to multitensor fitting methods in that a discrete number of diffusing compartments are included in the model. In this case, however, the model consists of one extra-axonal compartment (characterized with a single diffusion tensor), and a number of intra-axonal compartments each corresponding to a distinct fiber population (characterized using a model of restricted diffusion within cylinders) (132
). This approach requires a more demanding 3D q
-space acquisition to discriminate between the hindered and restricted components of the model.
Spherical deconvolution forms the basis of a number of recently proposed methods. Fundamentally, it extends the multitensor concept by increasing the number of fiber populations to infinity. In this limit, the summation becomes an integral over the distribution of fiber orientations, allowing the problem to be expressed as a spherical convolution (92
). By assuming a particular convolution kernel (representing the DW signal for a single fiber orientation), the fiber orientation distribution can be estimated by performing the spherical deconvolution operation (92
). Spherical deconvolution results obtained using constrained spherical deconvolution (99
) are shown in .
FIG. 14 Fiber orientation distributions for each voxel, for a coronal section showing the lateral projections of the corpus callosum (left-right: red lobes) crossing through the fibers of the corona radiata (inferior–superior: blue lobes) and of the superior (more ...)
The various implementations differ in the assumed convolution kernel, with some methods assuming a diffusion tensor model (93
), and others measuring it directly from the data (92
). They also differ in the constraints placed on the solution, with many implementations introducing a non-negativity constraint (99
), and others including a maximum entropy term (129
Apparent Diffusion Coefficient Models
Some methods have been proposed to characterizing the angular dependence of the ADC in the presence of crossing fibers. These include fitting spherical harmonics to the ADC profile (95
), or fitting a generalized diffusion tensor series to the ADC profile (118
). However, ADC-based methods have been shown to be inconsistent in the presence of non-gaussian diffusion (136
) and are, hence, inappropriate to characterize crossing fibers. Moreover, these methods cannot be used to infer the fiber orientations themselves without significant further processing (124
More recently, diffusion kurtosis imaging has been proposed as a way of characterizing non-gaussian behavior in DWI (137
). In its simplest form, it extents the simple linear relationship between the logarithm of the DW signal and b-value by introducing a quadratic term. This can be extended to 3D with a rank-4 diffusional kurtosis tensor (138
). These approaches provide promising new measures of tissue microstructure that relate to non-gaussian diffusion in the brain.
Data processing and Quantification
Given the increased complexity of nontensor models, it is inevitable that the amount of postprocessing will increase. However, in many cases, computation times can be kept very short: many algorithms can be expressed as linear matrix operations, and performed in seconds on modern workstations. These include DSI, QBI, diffusion orientation transform, and spherical deconvolution. Other algorithms, while requiring more computation, still provide results within a timeframe of a few minutes for a whole-brain data set (e.g., constrained spherical deconvolution). In other cases, the amount of computation can become prohibitive (hours or days), and limit the potential of these methods for clinical use.
Most of the methods described above (with the exception of ADC-based methods) aim to provide improved estimates of fiber orientations. While these are crucial for fiber tracking, there is a need for scalar measures to quantify features of the white matter that are insensitive to crossing fibers. Measures of diffusion anisotropy based on the diffusion tensor model have been used as surro-gate markers of white matter integrity in countless studies, but these are profoundly affected by crossing fibers. Unfortunately, deriving such crossing-fiber invariant measures is extremely difficult, as most white matter changes have relatively subtle effects on the DW signal. Nonetheless, diffusion kurtosis imaging, while still in its early stages of development, may provide biologically useful measures (137
), although it is unclear whether they are entirely unaffected by crossing fibers.
Some of the algorithms described above required a 3D q-space acquisition, whereby DW images are acquired with the DW gradients applied over a range of orientations and amplitudes. DSI for instance operates on data acquired with q-vectors arranged in a Cartesian grid, to facilitate the application of the subsequent 3D Fourier transform. The CHARMED model requires data acquired at multiple q-values per DW orientation, along multiple orientations. Moreover, these algorithms will require relatively large maximum q-values, which can only be achieved on clinical systems by increasing the echo time. For these reasons, the scan times required to acquire data suitable for these types of analyses tend to be long.
On the other hand, most algorithms make use of data acquired using the HARDI strategy (90
), whereby a relatively large number of DW directions are used with a constant b
-value. This allows the acquisition to focus on the angular part of the DW signal and select the most appropriate diffusion weighting so as to maximize contrast-to-noise per unit of scan time. The HARDI sequence is, therefore, arguably more efficient than the 3D q
-space acquisition for the purposes of fiber orientation estimation. However, the optimal b
-value and number of directions are both still the subject of ongoing research. This is partly a consequence of the wide range of algorithms and associated parameters that can be used to analyze HARDI data, making it difficult to design an experiment that would provide recommendations applicable to all these algorithms. Nevertheless, a number of recent studies have shown that b
-values in the range 2000–3000 s/mm2
provide the best power to resolve crossing fibers (50
). The number of DW directions required has also been the subject of a recent study, suggesting that the minimum number required is at least 28 for low b
~ 1000 s/mm2
), climbing to 45 for intermediate b
~ 3000 s/mm2
), although it is recommended to acquire a greater number of images to boost overall SNR (either as additional DW directions or as multiple repeats of the same DW orientations) (140
). Finally, the optimal distribution of orientations for a given number of DW directions has also been studied extensively (49
)], with a recent study suggesting that force-minimizing strategies [e.g., electro-static repulsion (49
)] and icosahedral schemes perform best in a multifiber setting (141