- Email: [email protected]
Crossing fibres in tract-based spatial statistics
NeuroImage 49 (2010) 249–256
Contents lists available at ScienceDirect
NeuroImage j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g
Crossing fibres in tract-based spatial statistics Saad Jbabdi a,⁎, Timothy E.J. Behrens a,b, Stephen M. Smith a a b
Oxford Centre for Functional Magnetic Resonance Imaging of the Brain, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK Department of Experimental Psychology, University of Oxford, South Parks Road, Oxford, UK
a r t i c l e
i n f o
Article history: Received 18 May 2009 Revised 6 August 2009 Accepted 17 August 2009 Available online 25 August 2009
a b s t r a c t Voxelwise analysis of white matter properties typically relies on scalar measurements derived, for example, from a tensor model fit to diffusion MRI data. These are spatially matched across subjects prior to statistical modelling. In this paper, we show why and how this can be improved through the use of directionally dependent measurements. In the case where different orientations relate to different fibre populations (e.g., in the presence of crossing fibres), distinguishing and matching those populations of fibres across subjects are important prior to any statistical modelling. It allows one to compare measurements that are related to the same fibres across subjects. We show how this framework applies to the parameters of a crossing fibre model and discuss its implications for voxelwise analysis of the white matter. © 2009 Elsevier Inc. All rights reserved.
Introduction Diffusion MRI is widely used to probe white matter variations across individuals. The tensor model (Basser et al., 1994) is still very popular for in vivo human studies, as it provides semi-quantitative scalar measures such as fractional anisotropy (FA) and mean apparent diffusivity (MD) that have been related to white matter microstructural “integrity” (e.g., Song et al., 2003). Voxelwise analysis of diffusion MRI aims at modelling or quantifying changes in white matter across individuals. This relies on the precise matching of anatomical locations across subjects, which is necessary in order for the experimenter to be confident that any result from the analysis is due to genuine changes in white matter microstructure rather than variations in brain structures' shape, size or position. For example, tract-based spatial statistics (TBSS; Smith et al., 2006) attempts to achieve this by restricting the statistical comparisons to the centres of white matter tracts after non-linear registration of different subjects into a common space. TBSS uses FA measurements to realign subjects and extract the centres of white matter tracts. The tensor model is based on the assumption that water selfdiffusion has a Gaussian diffusion profile. It does not assume any particular arrangement of the tissue microstructure, nor does it explicitly relate micro-structural features to the diffusion profile. This means that changes (e.g., across space or across subjects) in tensorderived scalar quantities such as FA and MD cannot always be related to micro-architecture in a straightforward, unambiguous and quantitative manner (Beaulieu, 2002). One claim that one can make from
⁎ Corresponding author. Fax: 44 1865 222717. E-mail address: [email protected] (S. Jbabdi). 1053-8119/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2009.08.039
a tensor model, with relative confidence, is that the principal diffusion direction is aligned with the main orientation of axonal fibres, at least when all the axons within a voxel have the same orientation. More complex models of the diffusion MR signal go beyond the Gaussian assumption of the tensor model, and some attempt to relate the signal to aspects of the underlying white matter microarchitecture. For example, one may consider the distinction between intra-cellular and extra-cellular diffusion and assume the former to be restricted by regular geometries such as cylinders (axons) and the latter hindered as a result of axonal packing (Assaf and Basser, 2005). Another extension of the tensor model may simply consist of considering the presence of distinct populations of fibre bundles with distinct orientations within a single imaging voxel. Using such crossing fibre models, one can infer on the orientations from the diffusion measurements (Tournier et al., 2004; Hosey et al., 2005; Parker and Alexander, 2005; Behrens et al., 2007). One may also consider that the contribution of each fibre population to the diffusion MR signal is related to the amount of space that is occupied by each fibre population. In the context of voxelwise white matter analysis, these partial volume fractions (PVFs) may be more interpretable than FA in locations where white matter bundles supporting distinct brain functions contribute to the signal within one voxel. We argue that voxelwise studies of the white matter could benefit from the use of such models. Including information about the relative amounts of specialised fibres within a voxel could enhance the interpretability of a finding. One would hope to be able to associate a local change in white matter with a particular fibre population. Our interpretation of any local change in white matter would therefore become more tract specific. In fact, interpreting FA changes in locations where functionally distinct fibre bundles are present is sometimes more challenging than simply knowing which fibres are
250
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
inducing the observed changes. The directionality of the changes, i.e., whether, for example, FA correlates positively or negatively with a variable of interest can also be ambiguous. In one of the first studies that reported a strong and tract-specific correlation between FA and a behavioural measure, Tuch et al. (2005) have noticed that FA along the optic radiations was negatively correlated with individual skills in a visuo-motor task. The effect was strong, but the sign of the correlations was counter-intuitive. The authors acknowledged the fact that this may be due to crossing fibres (for an illustration of this issue, see Fig. 1). This crossing fibre issue was also noted by Pierpaoli et al. (2001) and more recently by Wheeler-Kingshott and Cercignani (2009), where it was suggested that changes in tensor-derived scalar measurements (such as FA and diffusivity) should only be interpreted as being related to axonal integrity in locations where fibre bundles are aligned coherently. In this paper, we will focus on an extension of the TBSS framework to accounting for crossing fibre models. We will use as an example the crossing fibre model outlined in (Behrens et al., 2007). Extending TBSS to crossing fibre models requires some careful considerations with regard to the data, which we present in detail in the Methods section. The methodology is extremely straightforward and integrates naturally within the TBSS framework. Although we primarily discuss the extension of TBSS to crossing fibre models, the same methodology could potentially be applied within other frameworks for white matter analysis (Yushkevich et al., 2008; Eckstein et al., in press; Goodlett et al., 2009; O'Donnell et al., 2009).
Scalar reassignments based on directional data
Methods
This first step reassigns labels (i.e., swapping f1 and f2) within subjects by considering relationships across voxels. The purpose of this stage is to have smooth spatial orientation maps, such that adjacent voxels are likely to associate f1 with the same fibre population. We do not smooth the data but rather swap the labels (“1” and “2”) so that the associated vector maps are as smooth as possible. We use “front evolution” (Sethian, 2002) in order to perform this spatial regularisation. Starting from a region (a group of voxels) for which no reassignment is already done, the front evolution proceeds by repeating the following steps:
Throughout this paper, we will use the crossing fibre model described by Behrens et al. (2007). In this context, the diffusion signal is modelled as a weighted sum of signals accounting for the contribution of various compartments in each voxel: an infinitely anisotropic component for each fibre orientation and a single isotropic component. The weights of the signals from the anisotropic compartments will be denoted f1, f2, etc. The method described below may be generalized to any other crossing fibre model.
The idea behind using “directional data” within TBSS (or any other voxelwise method) is to be able to associate a voxelwise scalar measurement with an orientation. In the case of the crossing fibre model used here, in each voxel one may have two scalars f1 and f2, indicating the contribution of compartments 1 and 2 to the overall signal. In this context, f1 and f2 are associated with two different orientations x1 and x2, potentially indicating two distinct fibre populations within the voxel. It is important, when comparing f1 and f2 across subjects, that these scalars are associated with the same respective fibre populations in all subjects. It may happen that, for example, for a proportion of the subjects, f1 is associated with cortico-spinal tracts (e.g., x1 running along the z-direction) and for other subjects f1 accounts for associative pathways (e.g., x1 running along the y-direction). Assigning the labels 1 and 2 to the different fibres on the basis of the relative size of f1 and f2 (as happens by default in FSL) will not necessarily guarantee consistent labelling across subjects, particularly when the two tracts are of comparable “strength”. This can be corrected using label reassignments, i.e., swapping f1 and f2 for some voxels and some subjects such that these labels become consistent across subjects. This method may be applied to any number of directional measurements, as long as it is meaningful to match their orientations across subjects. We perform these reassignments in two steps. Intra-subject reassignment (spatial regularisation of fibre labelling)
1) For each voxel in the vicinity of the front, calculate the permutation of labels (1, 2, etc.) that aligns best (in terms of angle) their associated orientations with those of their neighbouring, previously considered, voxels near the front. 2) Find the voxel at the vicinity of the front, which aligns best (after permutation) with the neighbouring voxels from the front, and include that voxel in the front (which means that this voxel's associated measurements will not be reassigned any further). Starting a single front from a group of voxels and repeating the above two steps until no other voxel can be included in the front, we get a reassignment of all brain voxels that produces a smooth associated vector field. Fig. 2 shows the result of this procedure on a simulated vector field. Inter-subject reassignment (cross-subject regularisation of fibre labelling) Once scalars have been reassigned across space for each subject, we still need to reassign them so that the orientations are consistent across subjects. We do this in two steps:
Fig. 1. Illustration of a potential effect of crossing fibres on FA changes. The top row shows an FA increase due to an increase in the amount of vertical fibres. The increase in anisotropy is due to the fact that the vertical fibre population becomes the majority within the fibre crossing area. In the bottom row, an increase in this fibre population leads to a decrease in FA because the population that was a minority (in an anisotropic environment) became equal in proportion to the horizontal fibres, leading to a more isotropic environment at the scale of a voxel.
1) For each voxel, we calculate the modal (most common) orientation associated with each of the classes (or fibre populations). 2) For each subject (and each voxel), we permute the class labels and choose the permutation that gives the best alignment with the cross-subject modes. Fig. 2 shows the effect of this final reassignment on a simulated crossing fibre data set and compares the results of the two-stage
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
251
Fig. 2. An illustration of the two-stage process of label reassignment in simulated crossing fibre data with random noise on the orientations (five subjects simulated, one of them shown). The left panel shows the effect of permuting orientations 1 and 2 on the basis of intra-subject cross-voxel reassignment (the front evolution approach), followed by the cross-subject reassignment. The starting region for the front evolution is shown in the middle top (shaded area) and corresponds to the voxels with no crossing fibres. Note that after the first step, the vector maps are smooth but still imperfect, especially at the interfaces. These imperfections, however, are not consistent across subjects; hence, the second step achieves the desired result, i.e., grouping together coherent orientations. (In this simulation, only the central squares contain crossing fibres, so the remaining voxels do not show any coherence of the secondary orientation.) On the right panel, only the inter-subject reassignment is performed. This shows clearly the importance of the spatial regularisation (via the front evolution). Note that these simulations exaggerate this issue (compared with typical real data), for clarity.
reassignment to the case where only the second stage is performed (i.e., no intra-subject reassignment). Fig. 3 shows an example on real data, where we focus on the reassignment of skeleton voxels along the arcuate fasciculus. Practicalities Prior to the reassignment of labels of fibre populations across subjects and space, all subjects' diffusion-derived measurements need to be realigned onto a common space. The first step of TBSS consists of registering the FA images of each subject to an FA template in MNI152 standard space. This non-linear registration is carried out using FNIRT, which is part of FSL4.1 (www.fmrib.ox.ac. uk/fsl/fnirt). FNIRT generates a warp field for each subject that can be used to transform all the scalar measures from the crossing fibre model (f1, f2, etc.) along with their corresponding orientation images. Note that the transformation of the vector fields needs to be done with care. In particular, and as noted by Alexander et al. (2001), the warp fields contain an implicit local rotation that needs to be estimated and applied locally to each vector after registration. This is done using vecreg, which is also part of FSL4.1 (www.fmrib. ox.ac.uk/fsl/fdt). Finally, TBSS proceeds with projecting FA values onto the group-mean-FA white matter skeleton. The same projection is applied to the scalar and vector values in the crossing fibre model for each subject. Although the reassignment process is not restricted in validity to the white matter skeleton, we restrict it in our analysis to the TBSS white matter skeleton. Furthermore, we restrict it to the subset of the skeleton for which the data support more than one fibre population.1 1 The front evolution is performed on the whole skeleton, but only voxels supporting more than one fibre are reassigned. For the inter-subject reassignment, we reassign voxels where the group supports crossing fibres, regardless of whether there are crossing fibres or not at the individual level.
Determining whether or not the data in any given voxel support more than one fibre is achieved using the automatic relevance determination (ARD) framework detailed by Behrens et al. (2007). The ARD allows us to fit a crossing fibre model to the data (in a Bayesian framework), where the volume fraction of the secondary fibre population (f2) is shrunk to zero when the data do not support the presence of a secondary fibre population, thereby avoiding over-fitting. In locations where the ARD selects a single fibre model, the average orientation of the secondary population is random because changing this orientation does not affect the predicted signal. Therefore, we do not reassign these orientations because they can, by chance, realign better than the primary ones, even if they do not correspond to a true underlying white matter tract. Finally, the front evolution outlined above requires the specification of an initial region (or set of regions) where the front starts. We choose voxels where the ARD supports a single fibre population and where f1 is higher than a threshold value of 0.1. The reasoning behind the choice of not reassigning voxels that support a single orientation during the front evolution is as follows. If we were to reassign a voxel A with a single fibre orientation, there are two possibilities: (1) its neighbour B has a single fibre orientation (f2 is small), and (2) its neighbour B has two fibre orientations. In case 1, we do not want any reassignment for voxel A because f2 is related to an arbitrary orientation for both voxels A and B. In case 2, we do want to reassign x2 for voxel B if it is oriented along x1 in A; in such case, it is likely that these two orientations relate to the same tract. Reorienting x1 in A instead is less robust because x1 is in general more likely to be spatially smooth than x2. Using this approach, there might still be some issues at the interfaces between tracts that do not have the same orientations (as noticeable in Fig. 2, left), but this is generally resolved when performing the inter-subject reassignment.
252
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
Fig. 3. Example of the effect of the reassignment on real data in the arcuate fasciculus in 25 subjects. In panel a, we highlight subjects in which the orientation associated with f1 is not associated with the arcuate fasciculus in a subset of the skeleton (where it appears red rather than green), and in panel b, we show that the reassignment corrects for that. (c) Location of the arcuate region within an axial slice (top) and the modal orientation across subjects (bottom). Red/green/blue indicate left–right, rostral–caudal and inferior– superior, respectively.
Application We applied this approach to a group of 65 healthy adults (29 ± 8 years, 33 females). The data were acquired using a standard diffusion protocol as described by Tomassini et al. (2007). Here we use the PVFs in a conventional general linear model analysis, with age and gender as covariates. Inference is carried out using cluster thresholding, with a cluster forming threshold of t N 2.4. The null distribution of the cluster-size statistic was built up over 5000 permutations of the subjects, with the maximum cluster size recorded at each iteration. The 95th percentile was then used as the cluster-size threshold, i.e., clusters were thresholded at a level of p b 0.05, fully corrected for multiple comparisons across space. Fig. 4 shows the results of the correlations against age. Most of the voxels that show a significant effect are located within the frontal
white matter. All these clusters show a significant linear decrease of the scalar measurements with age. Two points are worth noting. First, while FA shows a widespread pattern of significant clusters, the PVF results seem to be more restricted, suggesting that FA is more sensitive but perhaps less specific to an effect of age. Second, the resulting significant clusters for f1 and f2 are spatially segregated and seem to be associated with a specific group of fibre tracts that are running anteriorly towards the frontal lobe. This suggests that the PVFs may provide more specific measurements for this type of analysis, relating age effects to actual tracts rather than to voxel locations (but see the Interpretation of the PVFs section). Fig. 5 illustrates the effect of the reassignments of f1/f2 on the regression analysis. Most of the voxels on the skeleton have a reassigned value in at least one subject, but some locations are more likely to be reassigned than others. In order to understand why this is,
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
253
Fig. 4. (a) Voxels of the skeleton showing a significant negative correlation between f1 (red), f2 (blue), FA (green) and age. (b) The expanded ROI view shows the two distinct fibre populations present in the f2 cluster (averaged across subjects). The negative correlation with age is associated here with the blue orientations.
Fig. 5. Effect of the label reassignment on the outcome of the regression analysis. Top left: spatial map showing the percentage of subjects for which there was a reassignment of f1/f2 within each voxel of the skeleton. It is clear from this figure that most of the reassignments occur in locations of crossing fibres and towards the periphery of the cortex. Top right: histogram summary of the top-left figure showing the distribution of swapped voxels against the percentage of subjects that were swapped. For example, the histogram shows that almost 14% of the total number of voxels on the skeleton were not swapped in any of the subjects, whereas roughly 2% were swapped in 20% of subjects. Bottom left: difference of the t-value between unswapped and swapped data. The t-values for both swapped and unswapped data sets were calculated in all voxels showing a significant negative correlation between f1 and age, then the difference between these two t-statistics was calculated for each voxel. These differences in t-values have then been sorted to visualise the asymmetry between negative and positive values (i.e., positive values indicate that the t-statistic is higher for the swapped than for the unswapped data). It is clear from this plot that the effect of the swap was to increase, on average, the effect size. The bar plots on the right show the average positive and negative differences in t-statistics (between swapped and unswapped data) for the various contrasts used in this paper. Again, we see that on average, the effects are larger after reassignment, although to a lesser extent than for the age versus f1 contrast. Note: all the t-statistics used in this figure are taken from skeleton voxels that are reported as significant (at a cluster level) for both reassigned and non-reassigned data, so that no bias is introduced in the choice of those voxels.
254
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
process. The effect size after reassignments, for f1 versus age, is most often higher than it was prior to reassignment. This suggests that the reassignment is acting sensibly, as we generally expect a zero or negative correlation between f1 and age, and so any improvement in the analysis methodology should increase the strength of this (negative) correlation (Salat et al., 2005). We also observed an average increase in effect size for the other contrasts, although to a lesser extent (Fig. 5). Finally, we have used the same general linear model to test for gender differences in the white matter skeleton (keeping the age regressor as a confound). This analysis showed a significant cluster for the contrast male N female for f1 and FA (at a similar location), but not f2 (Fig. 6, left). However, the contrast female N male only showed a significant cluster for f2 (right); this illustrates the decreased sensitivity of FA in crossing fibre regions, when the effect is concerned with secondary fibre populations. The relation between FA and f1, f2 is discussed further below. Discussion Interpretation of the partial volume fractions
Fig. 6. Gender effect on FA, f1 and f2. Left: clusters where f1 is significantly higher in males than females (the same clusters are significant for FA, but not f2; the FA cluster is not shown here because it is similar to f1). Right: clusters where f2 is significantly higher for females than males (no significant clusters were found for either f1 or FA) for this contrast. Bottom scatter plots: black = male, grey = female.
it is worth mentioning that when fitting a crossing fibre model such as the one used here, the choice of which fibre population represents f1 or f2 is arbitrary. The default in FDT (part of FSL4.1), which we used to fit the model here, is to sort fibre populations such that f1 is always higher than f2. This means that in voxels where two fibre populations are equally likely to co-exist (i.e., have equal contributions to the signal), it is more likely that these populations would not be matched across subjects since there will be more ambiguity as to which has a higher PVF. Fig. 5 confirms that this happens more often in crossing fibre areas. The figure also shows the potential increase in sensitivity of the regression analysis by means of the reassignment
In this paper, we illustrate the use of crossing fibres in TBSS using the partial volume model presented by Hosey et al. (2005) and Behrens et al. (2007). The diffusion signal is modelled as a weighted sum of signals accounting for the contribution of various compartments in each voxel: an infinitely anisotropic component for each fibre orientation and a single isotropic component. The weights, or partial volume fractions, are inferred from the data and represent the contribution of each compartment to the signal. In the discussion of Behrens et al. (2007), a clear point is made about the interpretation of the PVFs; unlike the more biologically motivated models of restricted diffusion such as CHARMED (Assaf and Basser, 2005) where the PVFs directly relate to compartment volumes within the tissues, the model used here does not make this assumption explicitly. While changes in the PVFs of the CHARMED model account for changes in the intraaxonal compartments and hence can be interpreted as changes in the relative density of each fibre population, such interpretation is not possible with the model used here. However, we would still argue that the biological interpretation, at least in crossing fibre areas, can be more clearly made (with our crossing fibre model) than with the single tensor model. Clearly, voxel-based analyses of white matter such as TBSS would greatly benefit from biologically motivated models. Such models, however, require more data than conventionally acquired in a clinical
Fig. 7. Relationship between PVFs and tensor-based scalars. The mean and standard deviation of the cross-subject correlation between PVFs and tensor-based values is calculated on the whole skeleton and sorted according to the angle (estimated) between the two crossing fibres. At high angles, the crossing fibre effect is more striking (and, presumably, all quantities more “cleanly” estimated). We can see, for example, that fractional anisotropy (FA) is much more correlated with f1 than f2. The mode (MO; see text) is correlated with both f1 and f2, but the correlations have opposite signs.
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
setting (e.g., multiple b-values), and hence, with current scanner technology, are generally only practical in studies involving a limited number of willing healthy volunteers. The simple crossing fibre model presented here benefits from its simplicity in that it requires less data. Below, we present a simple analysis that motivates the use of the simple partial volume model for white matter analysis in the presence of crossing fibres. At a local (voxelwise) level, we report here a closer look at the relationship between tensor-derived quantities and partial volume fractions from our crossing fibre model. Fig. 7 shows subject-wise correlations between PVFs and FA. We also report correlations with a tensor-derived quantity termed mode of anisotropy (MO), which was introduced in the context of DTI by Ennis and Kindlmann (2006). The mode distinguishes cases where anisotropy is linear (mode N 0, indicating a single fibre population) from cases where anisotropy is planar (mode b 0, potential crossing fibres). The correlations between tensor quantities and PVFs were calculated for all voxels of the white matter skeleton (across all 65 subjects) and then sorted according to the angle estimated between the two fibre orientations. Results are shown in Fig. 7. The first thing to note is that while FA relates strongly to f1 and increasingly so at a higher crossing fibre angle, it does not, on the other hand, relate strongly to f2 (or when it does, the correlation is negative). This is because FA does not distinguish well between cases of prolate anisotropy (single fibre, FA ∼ 1, MO ∼ 1) and oblate (planar) anisotropy (putative fibre crossing, FA ∼ 1/sqrt(2) ∼ 0.7, MO ∼ − 1). This is important as it shows that FA is relatively insensitive to changes in the white matter relating to secondary fibre populations (f2). By contrast, the mode of the diffusion tensor is explicitly designed to differentiate these two forms of anisotropy. This can be seen in Fig. 7 (right). The mode correlates positively with f1 and negatively with f2, particularly when the angle between the fibres is close to perpendicular (when the tensor is closest to oblate). Examining these two graphs (Fig. 7), one might wonder whether the complicated computations required to fit multi-fibre models are justified, given that the shape of the diffusion tensor carries related information. However, it is important to consider the profound differences between the two approaches. In a tensor representation, diffusivities are computed along angles that are constrained to be perpendicular. In cases when fibres cross at 90°, it is indeed possible that particular combinations of tensor eigenvalues might be used to assess the relative contribution of two fibre populations. However, in cases where fibre populations may cross at any angle, this angle confounds inference about individual sub-populations. For example, the exact same tensor might be recovered from two unequal populations crossing at 90° as would be recovered from two equal populations crossing at 45°. By contrast, multi-fibre models attempt to recover the individual contribution of each fibre, independent of the crossing angle. Hence, f1 and f2 clearly inform us about distinct fibre populations within a voxel (subject to the data quality supporting this complexity of model), and changes in f1/f2 may be interpreted as changes that are specific to those fibre populations. However, as noted above, we must bear in mind they do not relate explicitly to specific aspects of the microstructure. Finally, it is also clear that parallel (major eigenvalue) and transverse (mean of second and third eigenvalues) diffusivities are difficult to interpret in regions of crossing fibres for similar reasons (for more explicit investigations of these measures, see WheelerKingshott and Cercignani, 2009). What can we gain from directional measurements in TBSS? One may doubt the necessity of accounting for crossing fibres in voxelwise white matter analyses. FA is generally very sensitive (as shown in Fig. 4), and the potential problems with FA only occur in crossing fibre areas. However, one should keep in mind that most of
255
the white matter voxels support more than one fibre population with a distinct tract orientation, given data are readily available from modern MRI scanners and the crossing fibre model that we used here. Furthermore, while we have used reasonably high-resolution data in this manuscript ( 2-mm isotropic), it is very common in clinical routine to acquire data that have poorer resolution (e.g., 2.5–3 mm isotropic). This, we believe, will make the crossing fibre issue even more important. Partial volume effects will be greater, with different fibre bundles contributing to the data within single voxels (e.g., the corpus callosum can mix with the cingulum bundle within one voxel). In addition, as we have stated in the Introduction and shown on an example in the results, accounting for distinct fibre populations within a voxel hopefully allows the investigator to associate white matter changes with distinct white matter tracts, potentially subserving different functions, i.e., connecting functionally distinct cortical or sub-cortical regions. Finally, while FA is clearly confounded and hard to interpret at crossing fibre regions, other tensor parameters (such as MO) may be as sensitive to relative changes in fibre populations as a crossing fibre model. Fig. 7 illustrates this point as it shows a strong correlation of MO positively with f1 and negatively with f2. Hence, although the diffusion tensor-derived quantities may not be directly interpretable in terms of the underlying crossing fibres, they might still provide good statistics for voxelwise analyses. The practical value of using crossing fibre models, besides interpretability, still needs to be investigated and will probably be study specific and/or tract specific. Improvements There are clear potential further improvements to the present method and a few open questions that we would like to raise. Firstly, we have shown, as an example application, the use of PVFs from a simple crossing fibre model. This method, however, is by no means restricted to this particular model. Despite the promising results, which suggest a potential relevance of the PVFs in terms of the underlying micro-architecture, this framework would probably meet its best use if applied together with a more biologically motivated crossing fibre model, such as CHARMED. Contributions from intraaxonal water could then be interpreted in terms of axonal density. However, one should note that modelling the biophysics of water diffusion in tissues is far from resolved and still controversial (Le Bihan, 2007). Alternative models to intra/extra-cellular diffusion also exist and have had supporting evidence from empirical data. Notwithstanding the problem of finding the optimal crossing fibre model for a given data set, we believe that voxelwise white matter analyses will, in general, benefit from crossing fibre modelling both in terms of sensitivity and interpretability. On a more technical side, there are important potential refinements to the current methodology that are worth mentioning. One obvious example is the use of orientation information not only for reassigning associated scalar values across subjects but also in the alignment stage. For example, Park et al. (2003) have shown the value of including full diffusion tensor information in driving the non-linear registration. This might be of relevance here since if the FA-based initial alignment associates two different tracts across subjects, the second reassignment step will not resolve this issue. Another potential improvement concerns the skeleton projection stage. In TBSS, the projection is based on FA images. For each subject, the highest FA value is searched perpendicular to the skeleton, and this value is projected onto the skeleton, thereby assumed to correspond to an identical location across subjects. FA decreases due to crossing fibres might bias such an approach to projecting locations of single fibre population, rendering the modelling of crossing fibres less efficient in this framework. Alternatives may be driving the projection using f1 + f2 or the intra-cellular volume fraction in the case of biophysical models.
256
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
Acknowledgments This work was supported by the Medical Research Council. The authors would like to thank the many researchers in the FMRIB Centre who have contributed to the diffusion database that has been used in this article. References Alexander, D.C., Pierpaoli, C., Basser, P.J., Gee, J.C., 2001. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans. Med. Imaging 20, 1131–1139. Assaf, Y., Basser, P.J., 2005. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. NeuroImage 27, 48–58. Basser, P., Mattiello, J., Bihan, D.L., 1994. Estimation of the effective self-diffusion tensor from the NMR spin echo. J. Magn. Reson. 103, 247–254. Beaulieu, C., 2002. The basis of anisotropic water diffusion in the nervous system—a technical review. NMR. Biomed. 15, 435–455. Behrens, T.E., Berg, H.J., Jbabdi, S., Rushworth, M.F., Woolrich, M.W., 2007. Probabilistic diffusion tractography with multiple fibre orientations: what can we gain? NeuroImage 34, 144–155. Eckstein, I., Shattuck, D.W., Stein, J.L., McMahon, K.L., de Zubicaray, G., Wright, M.J., Thompson, P.M., Toga, A.W., 2009. Active fibers: matching deformable tract templates to diffusion tensor images. NeuroImage 47 (Suppl. 2), T82–T89. Ennis, D.B., Kindlmann, G., 2006. Orthogonal tensor invariants and the analysis of diffusion tensor magnetic resonance images. Magn. Reson. Med. 55, 136–146. Goodlett, C.B., Fletcher, P.T., Gilmore, J.H., Gerig, G., 2009. Group analysis of DTI fiber tract statistics with application to neurodevelopment. NeuroImage 45, S133–S142. Hosey, T., Williams, G., Ansorge, R., 2005. Inference of multiple fiber orientations in high angular resolution diffusion imaging. Magn. Reson. Med. 54, 1480–1489. Le Bihan, D., 2007. The ‘wet mind’: water and functional neuroimaging. Phys. Med. Biol. 52, R57–90. O'Donnell, L.J., Westin, C.F., Golby, A.J., 2009. Tract-based morphometry for white matter group analysis. NeuroImage 45, 832–844.
Park, H.J., Kubicki, M., Shenton, M.E., Guimond, A., McCarley, R.W., Maier, S.E., Kikinis, R., Jolesz, F.A., Westin, C.F., 2003. Spatial normalization of diffusion tensor MRI using multiple channels. NeuroImage 20, 1995–2009. Parker, G., Alexander, D., 2005. Probabilistic anatomical connectivity derived from the microscopic persistent angular structure of cerebral tissue. Philos. Trans. R. Soc. B 360, 893–902. Pierpaoli, C., Barnett, A., Pajevic, S., Chen, R., Penix, L.R., Virta, A., Basser, P., 2001. Water diffusion changes in Wallerian degeneration and their dependence on white matter architecture. NeuroImage 13, 1174–1185. Salat, D.H., Tuch, D.S., Hevelone, N.D., Fischl, B., Corkin, S., Rosas, H.D., Dale, A.M., 2005. Age-related changes in prefrontal white matter measured by diffusion tensor imaging. Ann. N. Y. Acad. Sci. 1064, 37–49. Sethian, J.A., 2002. Level Set Methods and Fast Marching Methods. Cambridge University Press, Cambridge. Smith, S.M., Jenkinson, M., Johansen-Berg, H., Rueckert, D., Nichols, T.E., Mackay, C.E., Watkins, K.E., Ciccarelli, O., Cader, M.Z., Matthews, P.M., Behrens, T.E.J., 2006. Tractbased spatial statistics: voxelwise analysis of multi-subject diffusion data. NeuroImage 31, 1487–1505. Song, S.K., Sun, S.W., Ju, W.K., Lin, S.J., Cross, A.H., Neufeld, A.H., 2003. Diffusion tensor imaging detects and differentiates axon and myelin degeneration in mouse optic nerve after retinal ischemia. NeuroImage 20, 1714–1722. Tomassini, V., Jbabdi, S., Klein, J.C., Behrens, T.E.J., Pozzilli, C., Matthews, P.M., Rushworth, M.F.S., Johansen-Berg, H., 2007. Diffusion-weighted imaging tractography-based parcellation of the human lateral premotor cortex identifies dorsal and ventral subregions with anatomical and functional specializations. J. Neurosci. 27, 10259–10269. Tournier, J.D., Calamante, F., Gadian, D.G., Connelly, A., 2004. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage 23, 1176–1185. Tuch, D.S., Salat, D.H., Wisco, J.J., Zaleta, A.K., Hevelone, N.D., Rosas, H.D., 2005. Choice reaction time performance correlates with diffusion anisotropy in white matter pathways supporting visuospatial attention. Proc. Natl. Acad. Sci. U. S. A. 102, 12212–12217. Wheeler-Kingshott, C.A., Cercignani, M., 2009. About “axial” and “radial” diffusivities. Magn. Reson. Med. 61, 1255–1260. Yushkevich, P.A., Zhang, H., Simon, T.J., Gee, J.C., 2008. Structure-specific statistical mapping of white matter tracts. NeuroImage 41, 448–461.
Contents lists available at ScienceDirect
NeuroImage j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g
Crossing fibres in tract-based spatial statistics Saad Jbabdi a,⁎, Timothy E.J. Behrens a,b, Stephen M. Smith a a b
Oxford Centre for Functional Magnetic Resonance Imaging of the Brain, John Radcliffe Hospital, University of Oxford, Oxford OX3 9DU, UK Department of Experimental Psychology, University of Oxford, South Parks Road, Oxford, UK
a r t i c l e
i n f o
Article history: Received 18 May 2009 Revised 6 August 2009 Accepted 17 August 2009 Available online 25 August 2009
a b s t r a c t Voxelwise analysis of white matter properties typically relies on scalar measurements derived, for example, from a tensor model fit to diffusion MRI data. These are spatially matched across subjects prior to statistical modelling. In this paper, we show why and how this can be improved through the use of directionally dependent measurements. In the case where different orientations relate to different fibre populations (e.g., in the presence of crossing fibres), distinguishing and matching those populations of fibres across subjects are important prior to any statistical modelling. It allows one to compare measurements that are related to the same fibres across subjects. We show how this framework applies to the parameters of a crossing fibre model and discuss its implications for voxelwise analysis of the white matter. © 2009 Elsevier Inc. All rights reserved.
Introduction Diffusion MRI is widely used to probe white matter variations across individuals. The tensor model (Basser et al., 1994) is still very popular for in vivo human studies, as it provides semi-quantitative scalar measures such as fractional anisotropy (FA) and mean apparent diffusivity (MD) that have been related to white matter microstructural “integrity” (e.g., Song et al., 2003). Voxelwise analysis of diffusion MRI aims at modelling or quantifying changes in white matter across individuals. This relies on the precise matching of anatomical locations across subjects, which is necessary in order for the experimenter to be confident that any result from the analysis is due to genuine changes in white matter microstructure rather than variations in brain structures' shape, size or position. For example, tract-based spatial statistics (TBSS; Smith et al., 2006) attempts to achieve this by restricting the statistical comparisons to the centres of white matter tracts after non-linear registration of different subjects into a common space. TBSS uses FA measurements to realign subjects and extract the centres of white matter tracts. The tensor model is based on the assumption that water selfdiffusion has a Gaussian diffusion profile. It does not assume any particular arrangement of the tissue microstructure, nor does it explicitly relate micro-structural features to the diffusion profile. This means that changes (e.g., across space or across subjects) in tensorderived scalar quantities such as FA and MD cannot always be related to micro-architecture in a straightforward, unambiguous and quantitative manner (Beaulieu, 2002). One claim that one can make from
⁎ Corresponding author. Fax: 44 1865 222717. E-mail address: [email protected] (S. Jbabdi). 1053-8119/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2009.08.039
a tensor model, with relative confidence, is that the principal diffusion direction is aligned with the main orientation of axonal fibres, at least when all the axons within a voxel have the same orientation. More complex models of the diffusion MR signal go beyond the Gaussian assumption of the tensor model, and some attempt to relate the signal to aspects of the underlying white matter microarchitecture. For example, one may consider the distinction between intra-cellular and extra-cellular diffusion and assume the former to be restricted by regular geometries such as cylinders (axons) and the latter hindered as a result of axonal packing (Assaf and Basser, 2005). Another extension of the tensor model may simply consist of considering the presence of distinct populations of fibre bundles with distinct orientations within a single imaging voxel. Using such crossing fibre models, one can infer on the orientations from the diffusion measurements (Tournier et al., 2004; Hosey et al., 2005; Parker and Alexander, 2005; Behrens et al., 2007). One may also consider that the contribution of each fibre population to the diffusion MR signal is related to the amount of space that is occupied by each fibre population. In the context of voxelwise white matter analysis, these partial volume fractions (PVFs) may be more interpretable than FA in locations where white matter bundles supporting distinct brain functions contribute to the signal within one voxel. We argue that voxelwise studies of the white matter could benefit from the use of such models. Including information about the relative amounts of specialised fibres within a voxel could enhance the interpretability of a finding. One would hope to be able to associate a local change in white matter with a particular fibre population. Our interpretation of any local change in white matter would therefore become more tract specific. In fact, interpreting FA changes in locations where functionally distinct fibre bundles are present is sometimes more challenging than simply knowing which fibres are
250
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
inducing the observed changes. The directionality of the changes, i.e., whether, for example, FA correlates positively or negatively with a variable of interest can also be ambiguous. In one of the first studies that reported a strong and tract-specific correlation between FA and a behavioural measure, Tuch et al. (2005) have noticed that FA along the optic radiations was negatively correlated with individual skills in a visuo-motor task. The effect was strong, but the sign of the correlations was counter-intuitive. The authors acknowledged the fact that this may be due to crossing fibres (for an illustration of this issue, see Fig. 1). This crossing fibre issue was also noted by Pierpaoli et al. (2001) and more recently by Wheeler-Kingshott and Cercignani (2009), where it was suggested that changes in tensor-derived scalar measurements (such as FA and diffusivity) should only be interpreted as being related to axonal integrity in locations where fibre bundles are aligned coherently. In this paper, we will focus on an extension of the TBSS framework to accounting for crossing fibre models. We will use as an example the crossing fibre model outlined in (Behrens et al., 2007). Extending TBSS to crossing fibre models requires some careful considerations with regard to the data, which we present in detail in the Methods section. The methodology is extremely straightforward and integrates naturally within the TBSS framework. Although we primarily discuss the extension of TBSS to crossing fibre models, the same methodology could potentially be applied within other frameworks for white matter analysis (Yushkevich et al., 2008; Eckstein et al., in press; Goodlett et al., 2009; O'Donnell et al., 2009).
Scalar reassignments based on directional data
Methods
This first step reassigns labels (i.e., swapping f1 and f2) within subjects by considering relationships across voxels. The purpose of this stage is to have smooth spatial orientation maps, such that adjacent voxels are likely to associate f1 with the same fibre population. We do not smooth the data but rather swap the labels (“1” and “2”) so that the associated vector maps are as smooth as possible. We use “front evolution” (Sethian, 2002) in order to perform this spatial regularisation. Starting from a region (a group of voxels) for which no reassignment is already done, the front evolution proceeds by repeating the following steps:
Throughout this paper, we will use the crossing fibre model described by Behrens et al. (2007). In this context, the diffusion signal is modelled as a weighted sum of signals accounting for the contribution of various compartments in each voxel: an infinitely anisotropic component for each fibre orientation and a single isotropic component. The weights of the signals from the anisotropic compartments will be denoted f1, f2, etc. The method described below may be generalized to any other crossing fibre model.
The idea behind using “directional data” within TBSS (or any other voxelwise method) is to be able to associate a voxelwise scalar measurement with an orientation. In the case of the crossing fibre model used here, in each voxel one may have two scalars f1 and f2, indicating the contribution of compartments 1 and 2 to the overall signal. In this context, f1 and f2 are associated with two different orientations x1 and x2, potentially indicating two distinct fibre populations within the voxel. It is important, when comparing f1 and f2 across subjects, that these scalars are associated with the same respective fibre populations in all subjects. It may happen that, for example, for a proportion of the subjects, f1 is associated with cortico-spinal tracts (e.g., x1 running along the z-direction) and for other subjects f1 accounts for associative pathways (e.g., x1 running along the y-direction). Assigning the labels 1 and 2 to the different fibres on the basis of the relative size of f1 and f2 (as happens by default in FSL) will not necessarily guarantee consistent labelling across subjects, particularly when the two tracts are of comparable “strength”. This can be corrected using label reassignments, i.e., swapping f1 and f2 for some voxels and some subjects such that these labels become consistent across subjects. This method may be applied to any number of directional measurements, as long as it is meaningful to match their orientations across subjects. We perform these reassignments in two steps. Intra-subject reassignment (spatial regularisation of fibre labelling)
1) For each voxel in the vicinity of the front, calculate the permutation of labels (1, 2, etc.) that aligns best (in terms of angle) their associated orientations with those of their neighbouring, previously considered, voxels near the front. 2) Find the voxel at the vicinity of the front, which aligns best (after permutation) with the neighbouring voxels from the front, and include that voxel in the front (which means that this voxel's associated measurements will not be reassigned any further). Starting a single front from a group of voxels and repeating the above two steps until no other voxel can be included in the front, we get a reassignment of all brain voxels that produces a smooth associated vector field. Fig. 2 shows the result of this procedure on a simulated vector field. Inter-subject reassignment (cross-subject regularisation of fibre labelling) Once scalars have been reassigned across space for each subject, we still need to reassign them so that the orientations are consistent across subjects. We do this in two steps:
Fig. 1. Illustration of a potential effect of crossing fibres on FA changes. The top row shows an FA increase due to an increase in the amount of vertical fibres. The increase in anisotropy is due to the fact that the vertical fibre population becomes the majority within the fibre crossing area. In the bottom row, an increase in this fibre population leads to a decrease in FA because the population that was a minority (in an anisotropic environment) became equal in proportion to the horizontal fibres, leading to a more isotropic environment at the scale of a voxel.
1) For each voxel, we calculate the modal (most common) orientation associated with each of the classes (or fibre populations). 2) For each subject (and each voxel), we permute the class labels and choose the permutation that gives the best alignment with the cross-subject modes. Fig. 2 shows the effect of this final reassignment on a simulated crossing fibre data set and compares the results of the two-stage
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
251
Fig. 2. An illustration of the two-stage process of label reassignment in simulated crossing fibre data with random noise on the orientations (five subjects simulated, one of them shown). The left panel shows the effect of permuting orientations 1 and 2 on the basis of intra-subject cross-voxel reassignment (the front evolution approach), followed by the cross-subject reassignment. The starting region for the front evolution is shown in the middle top (shaded area) and corresponds to the voxels with no crossing fibres. Note that after the first step, the vector maps are smooth but still imperfect, especially at the interfaces. These imperfections, however, are not consistent across subjects; hence, the second step achieves the desired result, i.e., grouping together coherent orientations. (In this simulation, only the central squares contain crossing fibres, so the remaining voxels do not show any coherence of the secondary orientation.) On the right panel, only the inter-subject reassignment is performed. This shows clearly the importance of the spatial regularisation (via the front evolution). Note that these simulations exaggerate this issue (compared with typical real data), for clarity.
reassignment to the case where only the second stage is performed (i.e., no intra-subject reassignment). Fig. 3 shows an example on real data, where we focus on the reassignment of skeleton voxels along the arcuate fasciculus. Practicalities Prior to the reassignment of labels of fibre populations across subjects and space, all subjects' diffusion-derived measurements need to be realigned onto a common space. The first step of TBSS consists of registering the FA images of each subject to an FA template in MNI152 standard space. This non-linear registration is carried out using FNIRT, which is part of FSL4.1 (www.fmrib.ox.ac. uk/fsl/fnirt). FNIRT generates a warp field for each subject that can be used to transform all the scalar measures from the crossing fibre model (f1, f2, etc.) along with their corresponding orientation images. Note that the transformation of the vector fields needs to be done with care. In particular, and as noted by Alexander et al. (2001), the warp fields contain an implicit local rotation that needs to be estimated and applied locally to each vector after registration. This is done using vecreg, which is also part of FSL4.1 (www.fmrib. ox.ac.uk/fsl/fdt). Finally, TBSS proceeds with projecting FA values onto the group-mean-FA white matter skeleton. The same projection is applied to the scalar and vector values in the crossing fibre model for each subject. Although the reassignment process is not restricted in validity to the white matter skeleton, we restrict it in our analysis to the TBSS white matter skeleton. Furthermore, we restrict it to the subset of the skeleton for which the data support more than one fibre population.1 1 The front evolution is performed on the whole skeleton, but only voxels supporting more than one fibre are reassigned. For the inter-subject reassignment, we reassign voxels where the group supports crossing fibres, regardless of whether there are crossing fibres or not at the individual level.
Determining whether or not the data in any given voxel support more than one fibre is achieved using the automatic relevance determination (ARD) framework detailed by Behrens et al. (2007). The ARD allows us to fit a crossing fibre model to the data (in a Bayesian framework), where the volume fraction of the secondary fibre population (f2) is shrunk to zero when the data do not support the presence of a secondary fibre population, thereby avoiding over-fitting. In locations where the ARD selects a single fibre model, the average orientation of the secondary population is random because changing this orientation does not affect the predicted signal. Therefore, we do not reassign these orientations because they can, by chance, realign better than the primary ones, even if they do not correspond to a true underlying white matter tract. Finally, the front evolution outlined above requires the specification of an initial region (or set of regions) where the front starts. We choose voxels where the ARD supports a single fibre population and where f1 is higher than a threshold value of 0.1. The reasoning behind the choice of not reassigning voxels that support a single orientation during the front evolution is as follows. If we were to reassign a voxel A with a single fibre orientation, there are two possibilities: (1) its neighbour B has a single fibre orientation (f2 is small), and (2) its neighbour B has two fibre orientations. In case 1, we do not want any reassignment for voxel A because f2 is related to an arbitrary orientation for both voxels A and B. In case 2, we do want to reassign x2 for voxel B if it is oriented along x1 in A; in such case, it is likely that these two orientations relate to the same tract. Reorienting x1 in A instead is less robust because x1 is in general more likely to be spatially smooth than x2. Using this approach, there might still be some issues at the interfaces between tracts that do not have the same orientations (as noticeable in Fig. 2, left), but this is generally resolved when performing the inter-subject reassignment.
252
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
Fig. 3. Example of the effect of the reassignment on real data in the arcuate fasciculus in 25 subjects. In panel a, we highlight subjects in which the orientation associated with f1 is not associated with the arcuate fasciculus in a subset of the skeleton (where it appears red rather than green), and in panel b, we show that the reassignment corrects for that. (c) Location of the arcuate region within an axial slice (top) and the modal orientation across subjects (bottom). Red/green/blue indicate left–right, rostral–caudal and inferior– superior, respectively.
Application We applied this approach to a group of 65 healthy adults (29 ± 8 years, 33 females). The data were acquired using a standard diffusion protocol as described by Tomassini et al. (2007). Here we use the PVFs in a conventional general linear model analysis, with age and gender as covariates. Inference is carried out using cluster thresholding, with a cluster forming threshold of t N 2.4. The null distribution of the cluster-size statistic was built up over 5000 permutations of the subjects, with the maximum cluster size recorded at each iteration. The 95th percentile was then used as the cluster-size threshold, i.e., clusters were thresholded at a level of p b 0.05, fully corrected for multiple comparisons across space. Fig. 4 shows the results of the correlations against age. Most of the voxels that show a significant effect are located within the frontal
white matter. All these clusters show a significant linear decrease of the scalar measurements with age. Two points are worth noting. First, while FA shows a widespread pattern of significant clusters, the PVF results seem to be more restricted, suggesting that FA is more sensitive but perhaps less specific to an effect of age. Second, the resulting significant clusters for f1 and f2 are spatially segregated and seem to be associated with a specific group of fibre tracts that are running anteriorly towards the frontal lobe. This suggests that the PVFs may provide more specific measurements for this type of analysis, relating age effects to actual tracts rather than to voxel locations (but see the Interpretation of the PVFs section). Fig. 5 illustrates the effect of the reassignments of f1/f2 on the regression analysis. Most of the voxels on the skeleton have a reassigned value in at least one subject, but some locations are more likely to be reassigned than others. In order to understand why this is,
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
253
Fig. 4. (a) Voxels of the skeleton showing a significant negative correlation between f1 (red), f2 (blue), FA (green) and age. (b) The expanded ROI view shows the two distinct fibre populations present in the f2 cluster (averaged across subjects). The negative correlation with age is associated here with the blue orientations.
Fig. 5. Effect of the label reassignment on the outcome of the regression analysis. Top left: spatial map showing the percentage of subjects for which there was a reassignment of f1/f2 within each voxel of the skeleton. It is clear from this figure that most of the reassignments occur in locations of crossing fibres and towards the periphery of the cortex. Top right: histogram summary of the top-left figure showing the distribution of swapped voxels against the percentage of subjects that were swapped. For example, the histogram shows that almost 14% of the total number of voxels on the skeleton were not swapped in any of the subjects, whereas roughly 2% were swapped in 20% of subjects. Bottom left: difference of the t-value between unswapped and swapped data. The t-values for both swapped and unswapped data sets were calculated in all voxels showing a significant negative correlation between f1 and age, then the difference between these two t-statistics was calculated for each voxel. These differences in t-values have then been sorted to visualise the asymmetry between negative and positive values (i.e., positive values indicate that the t-statistic is higher for the swapped than for the unswapped data). It is clear from this plot that the effect of the swap was to increase, on average, the effect size. The bar plots on the right show the average positive and negative differences in t-statistics (between swapped and unswapped data) for the various contrasts used in this paper. Again, we see that on average, the effects are larger after reassignment, although to a lesser extent than for the age versus f1 contrast. Note: all the t-statistics used in this figure are taken from skeleton voxels that are reported as significant (at a cluster level) for both reassigned and non-reassigned data, so that no bias is introduced in the choice of those voxels.
254
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
process. The effect size after reassignments, for f1 versus age, is most often higher than it was prior to reassignment. This suggests that the reassignment is acting sensibly, as we generally expect a zero or negative correlation between f1 and age, and so any improvement in the analysis methodology should increase the strength of this (negative) correlation (Salat et al., 2005). We also observed an average increase in effect size for the other contrasts, although to a lesser extent (Fig. 5). Finally, we have used the same general linear model to test for gender differences in the white matter skeleton (keeping the age regressor as a confound). This analysis showed a significant cluster for the contrast male N female for f1 and FA (at a similar location), but not f2 (Fig. 6, left). However, the contrast female N male only showed a significant cluster for f2 (right); this illustrates the decreased sensitivity of FA in crossing fibre regions, when the effect is concerned with secondary fibre populations. The relation between FA and f1, f2 is discussed further below. Discussion Interpretation of the partial volume fractions
Fig. 6. Gender effect on FA, f1 and f2. Left: clusters where f1 is significantly higher in males than females (the same clusters are significant for FA, but not f2; the FA cluster is not shown here because it is similar to f1). Right: clusters where f2 is significantly higher for females than males (no significant clusters were found for either f1 or FA) for this contrast. Bottom scatter plots: black = male, grey = female.
it is worth mentioning that when fitting a crossing fibre model such as the one used here, the choice of which fibre population represents f1 or f2 is arbitrary. The default in FDT (part of FSL4.1), which we used to fit the model here, is to sort fibre populations such that f1 is always higher than f2. This means that in voxels where two fibre populations are equally likely to co-exist (i.e., have equal contributions to the signal), it is more likely that these populations would not be matched across subjects since there will be more ambiguity as to which has a higher PVF. Fig. 5 confirms that this happens more often in crossing fibre areas. The figure also shows the potential increase in sensitivity of the regression analysis by means of the reassignment
In this paper, we illustrate the use of crossing fibres in TBSS using the partial volume model presented by Hosey et al. (2005) and Behrens et al. (2007). The diffusion signal is modelled as a weighted sum of signals accounting for the contribution of various compartments in each voxel: an infinitely anisotropic component for each fibre orientation and a single isotropic component. The weights, or partial volume fractions, are inferred from the data and represent the contribution of each compartment to the signal. In the discussion of Behrens et al. (2007), a clear point is made about the interpretation of the PVFs; unlike the more biologically motivated models of restricted diffusion such as CHARMED (Assaf and Basser, 2005) where the PVFs directly relate to compartment volumes within the tissues, the model used here does not make this assumption explicitly. While changes in the PVFs of the CHARMED model account for changes in the intraaxonal compartments and hence can be interpreted as changes in the relative density of each fibre population, such interpretation is not possible with the model used here. However, we would still argue that the biological interpretation, at least in crossing fibre areas, can be more clearly made (with our crossing fibre model) than with the single tensor model. Clearly, voxel-based analyses of white matter such as TBSS would greatly benefit from biologically motivated models. Such models, however, require more data than conventionally acquired in a clinical
Fig. 7. Relationship between PVFs and tensor-based scalars. The mean and standard deviation of the cross-subject correlation between PVFs and tensor-based values is calculated on the whole skeleton and sorted according to the angle (estimated) between the two crossing fibres. At high angles, the crossing fibre effect is more striking (and, presumably, all quantities more “cleanly” estimated). We can see, for example, that fractional anisotropy (FA) is much more correlated with f1 than f2. The mode (MO; see text) is correlated with both f1 and f2, but the correlations have opposite signs.
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
setting (e.g., multiple b-values), and hence, with current scanner technology, are generally only practical in studies involving a limited number of willing healthy volunteers. The simple crossing fibre model presented here benefits from its simplicity in that it requires less data. Below, we present a simple analysis that motivates the use of the simple partial volume model for white matter analysis in the presence of crossing fibres. At a local (voxelwise) level, we report here a closer look at the relationship between tensor-derived quantities and partial volume fractions from our crossing fibre model. Fig. 7 shows subject-wise correlations between PVFs and FA. We also report correlations with a tensor-derived quantity termed mode of anisotropy (MO), which was introduced in the context of DTI by Ennis and Kindlmann (2006). The mode distinguishes cases where anisotropy is linear (mode N 0, indicating a single fibre population) from cases where anisotropy is planar (mode b 0, potential crossing fibres). The correlations between tensor quantities and PVFs were calculated for all voxels of the white matter skeleton (across all 65 subjects) and then sorted according to the angle estimated between the two fibre orientations. Results are shown in Fig. 7. The first thing to note is that while FA relates strongly to f1 and increasingly so at a higher crossing fibre angle, it does not, on the other hand, relate strongly to f2 (or when it does, the correlation is negative). This is because FA does not distinguish well between cases of prolate anisotropy (single fibre, FA ∼ 1, MO ∼ 1) and oblate (planar) anisotropy (putative fibre crossing, FA ∼ 1/sqrt(2) ∼ 0.7, MO ∼ − 1). This is important as it shows that FA is relatively insensitive to changes in the white matter relating to secondary fibre populations (f2). By contrast, the mode of the diffusion tensor is explicitly designed to differentiate these two forms of anisotropy. This can be seen in Fig. 7 (right). The mode correlates positively with f1 and negatively with f2, particularly when the angle between the fibres is close to perpendicular (when the tensor is closest to oblate). Examining these two graphs (Fig. 7), one might wonder whether the complicated computations required to fit multi-fibre models are justified, given that the shape of the diffusion tensor carries related information. However, it is important to consider the profound differences between the two approaches. In a tensor representation, diffusivities are computed along angles that are constrained to be perpendicular. In cases when fibres cross at 90°, it is indeed possible that particular combinations of tensor eigenvalues might be used to assess the relative contribution of two fibre populations. However, in cases where fibre populations may cross at any angle, this angle confounds inference about individual sub-populations. For example, the exact same tensor might be recovered from two unequal populations crossing at 90° as would be recovered from two equal populations crossing at 45°. By contrast, multi-fibre models attempt to recover the individual contribution of each fibre, independent of the crossing angle. Hence, f1 and f2 clearly inform us about distinct fibre populations within a voxel (subject to the data quality supporting this complexity of model), and changes in f1/f2 may be interpreted as changes that are specific to those fibre populations. However, as noted above, we must bear in mind they do not relate explicitly to specific aspects of the microstructure. Finally, it is also clear that parallel (major eigenvalue) and transverse (mean of second and third eigenvalues) diffusivities are difficult to interpret in regions of crossing fibres for similar reasons (for more explicit investigations of these measures, see WheelerKingshott and Cercignani, 2009). What can we gain from directional measurements in TBSS? One may doubt the necessity of accounting for crossing fibres in voxelwise white matter analyses. FA is generally very sensitive (as shown in Fig. 4), and the potential problems with FA only occur in crossing fibre areas. However, one should keep in mind that most of
255
the white matter voxels support more than one fibre population with a distinct tract orientation, given data are readily available from modern MRI scanners and the crossing fibre model that we used here. Furthermore, while we have used reasonably high-resolution data in this manuscript ( 2-mm isotropic), it is very common in clinical routine to acquire data that have poorer resolution (e.g., 2.5–3 mm isotropic). This, we believe, will make the crossing fibre issue even more important. Partial volume effects will be greater, with different fibre bundles contributing to the data within single voxels (e.g., the corpus callosum can mix with the cingulum bundle within one voxel). In addition, as we have stated in the Introduction and shown on an example in the results, accounting for distinct fibre populations within a voxel hopefully allows the investigator to associate white matter changes with distinct white matter tracts, potentially subserving different functions, i.e., connecting functionally distinct cortical or sub-cortical regions. Finally, while FA is clearly confounded and hard to interpret at crossing fibre regions, other tensor parameters (such as MO) may be as sensitive to relative changes in fibre populations as a crossing fibre model. Fig. 7 illustrates this point as it shows a strong correlation of MO positively with f1 and negatively with f2. Hence, although the diffusion tensor-derived quantities may not be directly interpretable in terms of the underlying crossing fibres, they might still provide good statistics for voxelwise analyses. The practical value of using crossing fibre models, besides interpretability, still needs to be investigated and will probably be study specific and/or tract specific. Improvements There are clear potential further improvements to the present method and a few open questions that we would like to raise. Firstly, we have shown, as an example application, the use of PVFs from a simple crossing fibre model. This method, however, is by no means restricted to this particular model. Despite the promising results, which suggest a potential relevance of the PVFs in terms of the underlying micro-architecture, this framework would probably meet its best use if applied together with a more biologically motivated crossing fibre model, such as CHARMED. Contributions from intraaxonal water could then be interpreted in terms of axonal density. However, one should note that modelling the biophysics of water diffusion in tissues is far from resolved and still controversial (Le Bihan, 2007). Alternative models to intra/extra-cellular diffusion also exist and have had supporting evidence from empirical data. Notwithstanding the problem of finding the optimal crossing fibre model for a given data set, we believe that voxelwise white matter analyses will, in general, benefit from crossing fibre modelling both in terms of sensitivity and interpretability. On a more technical side, there are important potential refinements to the current methodology that are worth mentioning. One obvious example is the use of orientation information not only for reassigning associated scalar values across subjects but also in the alignment stage. For example, Park et al. (2003) have shown the value of including full diffusion tensor information in driving the non-linear registration. This might be of relevance here since if the FA-based initial alignment associates two different tracts across subjects, the second reassignment step will not resolve this issue. Another potential improvement concerns the skeleton projection stage. In TBSS, the projection is based on FA images. For each subject, the highest FA value is searched perpendicular to the skeleton, and this value is projected onto the skeleton, thereby assumed to correspond to an identical location across subjects. FA decreases due to crossing fibres might bias such an approach to projecting locations of single fibre population, rendering the modelling of crossing fibres less efficient in this framework. Alternatives may be driving the projection using f1 + f2 or the intra-cellular volume fraction in the case of biophysical models.
256
S. Jbabdi et al. / NeuroImage 49 (2010) 249–256
Acknowledgments This work was supported by the Medical Research Council. The authors would like to thank the many researchers in the FMRIB Centre who have contributed to the diffusion database that has been used in this article. References Alexander, D.C., Pierpaoli, C., Basser, P.J., Gee, J.C., 2001. Spatial transformations of diffusion tensor magnetic resonance images. IEEE Trans. Med. Imaging 20, 1131–1139. Assaf, Y., Basser, P.J., 2005. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. NeuroImage 27, 48–58. Basser, P., Mattiello, J., Bihan, D.L., 1994. Estimation of the effective self-diffusion tensor from the NMR spin echo. J. Magn. Reson. 103, 247–254. Beaulieu, C., 2002. The basis of anisotropic water diffusion in the nervous system—a technical review. NMR. Biomed. 15, 435–455. Behrens, T.E., Berg, H.J., Jbabdi, S., Rushworth, M.F., Woolrich, M.W., 2007. Probabilistic diffusion tractography with multiple fibre orientations: what can we gain? NeuroImage 34, 144–155. Eckstein, I., Shattuck, D.W., Stein, J.L., McMahon, K.L., de Zubicaray, G., Wright, M.J., Thompson, P.M., Toga, A.W., 2009. Active fibers: matching deformable tract templates to diffusion tensor images. NeuroImage 47 (Suppl. 2), T82–T89. Ennis, D.B., Kindlmann, G., 2006. Orthogonal tensor invariants and the analysis of diffusion tensor magnetic resonance images. Magn. Reson. Med. 55, 136–146. Goodlett, C.B., Fletcher, P.T., Gilmore, J.H., Gerig, G., 2009. Group analysis of DTI fiber tract statistics with application to neurodevelopment. NeuroImage 45, S133–S142. Hosey, T., Williams, G., Ansorge, R., 2005. Inference of multiple fiber orientations in high angular resolution diffusion imaging. Magn. Reson. Med. 54, 1480–1489. Le Bihan, D., 2007. The ‘wet mind’: water and functional neuroimaging. Phys. Med. Biol. 52, R57–90. O'Donnell, L.J., Westin, C.F., Golby, A.J., 2009. Tract-based morphometry for white matter group analysis. NeuroImage 45, 832–844.
Park, H.J., Kubicki, M., Shenton, M.E., Guimond, A., McCarley, R.W., Maier, S.E., Kikinis, R., Jolesz, F.A., Westin, C.F., 2003. Spatial normalization of diffusion tensor MRI using multiple channels. NeuroImage 20, 1995–2009. Parker, G., Alexander, D., 2005. Probabilistic anatomical connectivity derived from the microscopic persistent angular structure of cerebral tissue. Philos. Trans. R. Soc. B 360, 893–902. Pierpaoli, C., Barnett, A., Pajevic, S., Chen, R., Penix, L.R., Virta, A., Basser, P., 2001. Water diffusion changes in Wallerian degeneration and their dependence on white matter architecture. NeuroImage 13, 1174–1185. Salat, D.H., Tuch, D.S., Hevelone, N.D., Fischl, B., Corkin, S., Rosas, H.D., Dale, A.M., 2005. Age-related changes in prefrontal white matter measured by diffusion tensor imaging. Ann. N. Y. Acad. Sci. 1064, 37–49. Sethian, J.A., 2002. Level Set Methods and Fast Marching Methods. Cambridge University Press, Cambridge. Smith, S.M., Jenkinson, M., Johansen-Berg, H., Rueckert, D., Nichols, T.E., Mackay, C.E., Watkins, K.E., Ciccarelli, O., Cader, M.Z., Matthews, P.M., Behrens, T.E.J., 2006. Tractbased spatial statistics: voxelwise analysis of multi-subject diffusion data. NeuroImage 31, 1487–1505. Song, S.K., Sun, S.W., Ju, W.K., Lin, S.J., Cross, A.H., Neufeld, A.H., 2003. Diffusion tensor imaging detects and differentiates axon and myelin degeneration in mouse optic nerve after retinal ischemia. NeuroImage 20, 1714–1722. Tomassini, V., Jbabdi, S., Klein, J.C., Behrens, T.E.J., Pozzilli, C., Matthews, P.M., Rushworth, M.F.S., Johansen-Berg, H., 2007. Diffusion-weighted imaging tractography-based parcellation of the human lateral premotor cortex identifies dorsal and ventral subregions with anatomical and functional specializations. J. Neurosci. 27, 10259–10269. Tournier, J.D., Calamante, F., Gadian, D.G., Connelly, A., 2004. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage 23, 1176–1185. Tuch, D.S., Salat, D.H., Wisco, J.J., Zaleta, A.K., Hevelone, N.D., Rosas, H.D., 2005. Choice reaction time performance correlates with diffusion anisotropy in white matter pathways supporting visuospatial attention. Proc. Natl. Acad. Sci. U. S. A. 102, 12212–12217. Wheeler-Kingshott, C.A., Cercignani, M., 2009. About “axial” and “radial” diffusivities. Magn. Reson. Med. 61, 1255–1260. Yushkevich, P.A., Zhang, H., Simon, T.J., Gee, J.C., 2008. Structure-specific statistical mapping of white matter tracts. NeuroImage 41, 448–461.