Based on this observation, the most obvious approach to resolving crossing fibers is the multi-tensor approach: we assume the voxel contains two fiber populations rather than one, with each population modeled using its own diffusion tensor. In this case, the minimum additional information required is the relative amounts of each fiber population (i.e., their volume fractions) and their respective orientations. The problem is then one of estimating the set of parameters (volume fractions and orientations) that best fit the measured data [4, 30].

It turns out that this is not a simple problem to solve in practice due to the nonlinear formulation of the problem (a sum of exponentials), dictating the use of computationally expensive nonlinear minimization methods . In addition, the problem is generally ill-conditioned, due to the difficulty of differentiating between changes in anisotropy and changes in volume fraction of the constituent fiber populations¹; to get around this, most implementations enforce the condition that each diffusion tensor be at least axially symmetric, and often hold all their diffusivities constant (i.e., assume a fixed mean diffusivity and anisotropy). Nonetheless, a number of methods are based on this general approach; some allow for larger numbers of fiber populations, some include an isotropic “CSF” compartment, and some use advanced Bayesian Monte Carlo Markhov Chain (MCMC) sampling techniques to characterize the uncertainty about the parameters estimates [14, 30–32].

The Combined Hindered and Restricted Model of Diffusion (CHARMED) approach is strongly related to the multi-tensor approach in that it also models each fiber population independently [33, 34]. However, it differs substantially in that it assumes a more biologically plausible model of restricted diffusion in cylinders to represent the DW signal from each fiber populations, rather than a diffusion tensor model. It also includes an extracellular compartment, itself modeled by a diffusion tensor. One limitation is that it requires a more demanding multi-shell (i.e., multiple b-values) HARDI data acquisition, leading to lengthier scan times.

One issue with these methods is that they somehow need to know how many fiber populations to include in the model for each voxel. While a number of approaches have been proposed to do this [2, 14, 31], any errors in this number will inevitably lead to errors in the estimated parameters. Most of these approaches are based on some statistical “goodness of fit” measure, and this will in general lead to underestimation of the number of fiber populations present, particularly in noisy data where any improvements in the model fit will be overwhelmed by the noise. The result is that while the algorithm can be said to be “conservative” in that it only includes additional fiber populations if there is good evidence that these are needed, it will inevitably fail to identify a significant proportion of crossing fiber voxels and therefore model them using the single tensor model , with all the limitations highlighted above.

The multi-tensor approach can be extended by switching from a discrete representation (i.e., a small number of distinct fiber populations) to a continuous representation of the fiber orientation information. One way to understand this is to imagine that rather than trying to estimate the volume fraction and orientation of each fiber population, we are now going to model a large number of fiber populations, each with a fixed orientation. The problem now reduces to figuring out the volume fractions of each of these. If we use a sufficiently dense set of orientations that are uniformly distributed over the sphere, the information can be represented as a distribution of volume fractions over the sphere, as illustrated in Fig. 20.4. This distribution is commonly referred to as the fiber orientation density function (fODF ) or simply the fiber orientation distribution (FOD ).

Essentially, the FOD is a function defined over the sphere, which represents the amount of fibers at any point on the sphere (i.e., aligned with the corresponding orientation). It is often displayed using orientation plots such as those shown in Fig. 20.5, whereby the distance from the origin represents the “density” along the corresponding orientation. Distinct fiber orientations can clearly be seen as distinct peaks in the FOD.

The advantages of using such a distribution to represent the fiber orientation information are threefold. First, it means that the relationship between the measured DW signal and the FOD is linear, allowing the use of much more efficient reconstruction methods based on linear algebra . Second, there is no need to specify the number of directions that might be present in each voxel; the FOD is continuous, and distinct fiber orientations will simply be reflected as distinct peaks in the FOD. Finally, the FOD is suited to describing fiber configurations containing a range of orientations, such as for example bending or fanning fibers; these are clearly not well characterized using a single orientation, as would be used in multi-tensor approaches. In these cases, the corresponding peak in the FOD will simply be broadened accordingly, as expected. Spherical deconvolution approaches are therefore more computationally efficient, while providing a much more general description of the fiber orientation information. These provide estimates of the FOD that are suitable for fiber-tracking , as illustrated for example in Fig. 20.6.

A number of different spherical deconvolution techniques have been proposed [10, 29, 36–39]. The most stable of these impose a non-negativity constraint to explicitly prevent or minimize negative fiber density values, which are clearly non-physical; this has been shown to dramatically improve the robustness to noise of these methods [29, 37]. Some estimate the per-fiber anisotropy on a per-voxel basis [36], other assume fixed diffusivity values [37–39], while others bypass the diffusion tensor model entirely and measure the actual per-fiber DW signal profile from the data themselves [10, 29].

The overwhelming majority of multi-tensor fitting and spherical deconvolution techniques will make the assumption that all fiber bundles can be described by diffusion tensors with the same intrinsic “constant anisotropy” (or at least axial symmetry) to reduce the complexity of the problem. At first sight, this might seem to be a gross over-simplification, particular to those steeped in the more traditional DTI framework, where tensor-derived metrics are often used as surrogate markers of white matter “integrity.” Since these new approaches assume fixed values for these metrics over all white matter bundles, it is now impossible to estimate equivalent measures, let alone detect differences between them.

However, there are very good reasons to believe that this approximation holds in most, if not all cases. The assumption of fixed anisotropy inherently implies that all observed anisotropy differences in the brain are due entirely to crossing fiber effects or other partial volume effects (e.g., with CSF). This is in fact in line with the original DTI literature, where the large variations of anisotropy observed in healthy volunteers were attributed to differences in the coherence of fiber tract directions [24]. The relationship between anisotropy (and more recently radial and axial diffusivities) and white matter “integrity” was suggested based on highly controlled single nerve experiments, where any confounding effects of crossing fibers could explicitly be ruled out [40–44]. While these studies undeniably demonstrated changes in various tensor-derived metrics under different conditions (dysmyelination, axonal injury, changes in axonal diameter, etc.), these changes are dwarfed by the changes induced by the presence of crossing fibers [1, 23]. Moreover, and more importantly, while changes in diffusivities and/or anisotropy can be observed in these experiments, they are not of a magnitude sufficient to alter the results dramatically; as shown in previous studies, while moderate changes in the assumed anisotropy will affect the estimated volume fractions, they have almost no impact on the estimated fiber orientations [10]. In fact, any pathology that would cause severe changes in the anisotropy of the per-fiber DW signal would probably need to involve destruction of the axonal membrane, with a large increase in radial diffusivity and corresponding reduction in the DW signal; this will be interpreted by these models primarily as a dramatic decrease in volume fraction, an interpretation that is in fact fairly accurate given the pathology.

For this reason, this particular assumption is very commonly employed, and with great success. However, the implication is that genuine anisotropy changes are at best extremely difficult to detect from HARDI data (given that they are essentially indistinguishable from partial volume effects), which essentially rules out the possibility of deriving equivalent markers of white matter “integrity”—at least using these approaches. Instead, the current trend is to focus on the estimated partial volume fractions, since these seem to have the most dominant impact on the DW signal, and can be used as estimates of fiber density—a measure that, although different from the more common interpretation of “integrity,” is nonetheless clearly clinically very relevant [45–47].

These methods are all in some way based on the theory of q-space [48], which describes the relationship between the DW signal and the so-called average spin propagator (also called the spin displacement probability density function ). In simple terms, the spin propagator is a function describing the probability P(r|x,Δ) that a water molecule initially at x has moved by displacement r during the diffusion time Δ of the MR measurement. For free diffusion, this is simply a Gaussian distribution with standard deviation equal to the root mean square displacement of the water molecules (as per Chap. 3). For more complex environments, the propagator will have a more complex shape, as illustrated in Fig. 20.7. This makes no assumption about whether diffusion is free or restricted: it simply characterizes the chances of molecules moving along a given direction by a given amount. In general, this will depend on the diffusion time Δ: with free diffusion, the distances moved will increase with diffusion time. However, in restricted geometries the maximum distance moved will clearly be dictated by the size of the restricting compartment. The relationship between the DW signal and the spin propagator is relatively simple and not subject to any particular biological assumptions.² Clearly, the diffusion propagator provides the most complete and accurate description of the diffusion process, and contains all the information that can be extracted from diffusion MRI, which is the reason why a number of methods attempt to estimate it.

According to q-space, the DW signal is related to the spin propagator via a simple Fourier transform. In practice, this means that to estimate the spin propagator along one dimension, a number of DW images need to be acquired, each with a different q-value.³ To get an accurate estimate, a sufficient number of distinct q-values need to be acquired, up to a sufficiently large maximum q-value to ensure near-complete nulling of the DW signal. These two requirements alone are difficult to achieve: for a full three-dimensional characterization of the spin propagator, the number of distinct q-vectors (and hence image volumes) required is of the order of 1000 (approx. 10 per image axis). These requirements also dictate the use of very large gradients and/or long echo times (to achieve the very large maximum q-values needed), leading to a reduction in the SNR. Clearly, the full characterization of the spin propagator is extremely challenging.

Nonetheless, this complete characterization is exactly what Diffusion Spectrum Imaging ( DSI ) aims to achieve [49]. To do this, 515 image volumes are acquired, each with a distinct q-vector with a maximum b-value of 17,000 s/mm². This is clearly a much more demanding acquisition than can realistically be accommodated in routine clinical imaging, and for this reason much of the recent development in this area has been focused on obtaining relevant features of the spin propagator (namely the diffusion orientation density function; dODF) using the much more clinically achievable HARDI acquisition.

Essentially, the dODF is a simplified version of the full spin propagator, which provides only the probability of a water molecule moving along a given direction; any information about how far it moved is discarded. This much more compact version of the spin propagator is obviously still sufficiently informative, especially for the purposes of resolving crossing fibers, motivating the development of these HARDI-based methods.

Q–ball Imaging (QBI ) was the first method proposed to estimate the dODF [50]. It performs a Funk–Radon Transform on the HARDI data to map the DW signal per-voxel onto the corresponding dODF. This has the advantage of being a linear operation, and was subsequently shown to be much more efficiently and robustly implemented using spherical harmonics [51, 52]. More recently, Constant Solid Angle Q–ball Imaging (CSA -QBI ) has been proposed to obtain a more accurate , sharper estimate of the dODF [53].

Other techniques based on q-space that have been proposed include the Diffusion Orientation Transform (DOT ), which provides a contour of the spin propagator evaluated at a particular value of the displacement—in other words, it provides the probability of a water molecule moving along any given direction by a fixed distance [54]. Another early technique is Persistent Angular Structure MRI (PAS-MRI ), which evaluates an ODF from relatively modest data by imposing a maximum entropy constraint [55]. See Fig. 20.8 for an illustration of the results obtained using three different dODF estimation methods.

The theory of q-space does rely on one approximation: that the duration of the DW gradient pulses is negligible—in other words, the DW gradient pulses can be considered to be infinitesimally narrow. More specifically, the distance moved by molecules during the application of each DW gradient pulse should be much smaller than the distance moved between the pulses. The idea is that each DW gradient pulse respectively “tags” and “untags” water molecules based on their current location, so that the only relevant quantity is their actual displacement between these two events. For this approximation to hold in white matter, water molecules should not be allowed to diffuse by any distance approaching the axonal diameter. Even if we assume an acceptable distance to be 1 μm (which is already larger than many axons) and a diffusion coefficient of 10⁻³ mm²/s (which is smaller than would typically be assumed), the DW pulse duration would need to be less than 0.5 ms. To achieve any reasonable b-value with such a short pulse duration would require gradients strengths and rise-times orders of magnitude larger than can be used in the clinic. In practice, the minimum DW pulse duration than can realistically be achieved on a human system is of the order of 10 ms. Clearly, the narrow pulse approximation cannot hold on a clinical system.