Fibre density and cross-section - Multi-tissue CSD

Introduction

This tutorial explains how to perform fixel-based analysis of fibre density and cross-section [Raffelt2017] with fibre orientation distributions (FODs) computed using multi-tissue (3-tissue) CSD variants [Jeurissen2014] [Dhollander2016a]. We note that high b-value (>2000s/mm2) data is recommended to aid the interpretation of apparent fibre density (AFD) being related to the intra-axonal space. See [Raffelt2012] for some details about AFD; though note that the interpretation can be altered for multi-tissue (3-tissue) CSD, depending on the context and tissues in the model.

All steps in this tutorial are written as if the commands are being run on a cohort of images, and make extensive use of the for_each script to simplify batch processing. This tutorial also assumes that the imaging dataset is organised with one directory identifying each subject, and all files within identifying the image type (i.e. processing step outcome). For example:

study/subjects/001_control/dwi.mif
study/subjects/002_control/dwi.mif
...
study/subjects/020_control/dwi.mif
study/subjects/021_patient/dwi.mif
...
study/subjects/040_patient/dwi.mif

Note

All commands at the start of this tutorial are run from the subjects path. From the step where tractography is performed on the template onwards, we change directory to the template path.

For all MRtrix scripts and commands, additional information on the command usage and available command-line options can be found by invoking the command with the -help option.

Pre-processsing steps

1. Denoising and unringing

If denoising and/or Gibbs ringing removal are performed as part of the preprocessing, they must be performed prior to any other processing steps: most other processing steps, in particular those that involve interpolation of the data, will invalidate the original properties of the image data that are exploited by dwidenoise and mrdegibbs at this stage, and would render the result prone to errors.

If denoising is included, it’s performed as the first step:

for_each * : dwidenoise IN/dwi.mif IN/dwi_denoised.mif

If Gibbs ringing removal is included, it follows immediately after:

for_each * : mrdegibbs IN/dwi_denoised.mif IN/dwi_denoised_unringed.mif -axes 0,1

Warning

The -axes option to mrdegibbs is used to specify in which plane the slices were acquired. The -axes 0,1 in the example above refers to the x-y plane, which is appropriate for data consisting of a stack of axial slices (assuming a typical human scanner and subject). For typical human data, change this to -axes 0,2 for coronal slices or -axes 1,2 for sagittal slices.

2. Motion and distortion correction

The dwifslpreproc command handles motion and distortion correction for DWI data (including eddy current distortions and optionally susceptibility-induced EPI distortions). Even though the command works seamlessly like any other MRtrix3 command, it is in fact a script that interfaces with the FSL package to perform most of its core functionality and algorithms. For this command to work, FSL (including eddy) needs to be installed. Also remember to cite the relevant articles with respect to the specific algorithms (see the dwifslpreproc help page).

The simplest scenario is to (only) correct for motion and eddy current-induced distortions:

for_each * : dwifslpreproc IN/dwi_denoised_unringed.mif IN/dwi_denoised_unringed_preproc.mif -rpe_none -pe_dir AP

Warning

The -pe_dir option to dwifslpreproc is used to specify the phase encoding direction of the acquisition. The -pe_dir AP in the example above refers to an anterior-posterior phase encoding direction, which is relatively commonly used for acquiring human data. For typical human data, change this to -pe_dir LR for a left-right phase encoding direction or -pe_dir SI for a superior-inferior phase encoding direction.

For different, more advanced, scenarios and acquisitions (e.g. the common case of correcting for susceptibility-induced EPI distortions using a pair of reverse phase-encoded b=0 images), refer to the DWI distortion correction using dwifslpreproc section and the dwifslpreproc documentation.

3. Bias field correction

The multi-tissue FBA pipeline corrects for bias fields (and jointly performs global intensity normalisation) at the later mtnormalise step. The only incentive for running the (less robust and accurate) dwibiascorrect at this stage in the pipeline is to improve brain mask estimation (at the later dwi2mask step, in case severe bias fields are present in the data). However, cases have been reported where running dwibiascorrect at this stage resulted in inferior brain mask estimation later on. This is probably more likely in case bias fields are not as strongly present in the data. Whether dwibiascorrect is run at this stage or not, does not have any significant impact on the performance of mtnormalise later on.

If or when performing DWI bias field correction at this stage, it is achieved by first estimating the bias field from the DWI b=0 data, then applying the field to correct all DW volumes, which is done in a single step using the ants algorithm within the dwibiascorrect script in MRtrix3. The script uses a bias field correction algorithm available in ANTs (the N4 algorithm). Don’t use the fsl algorithm with this script in this fixel-based analysis pipeline. To perform bias field correction on DW images, run:

for_each * : dwibiascorrect ants IN/dwi_denoised_unringed_preproc.mif IN/dwi_denoised_unringed_preproc_unbiased.mif

Fixel-based analysis steps

4. Computing (average) tissue response functions

A robust and fully automated unsupervised method to obtain 3-tissue response functions representing single-fibre white matter, grey matter and CSF from the data itself, is the approach proposed in [Dhollander2016b] with the improvements of [Dhollander2019], which can be run by:

for_each * : dwi2response dhollander IN/dwi_denoised_unringed_preproc_unbiased.mif IN/response_wm.txt IN/response_gm.txt IN/response_csf.txt

It is crucial for fixel-based analysis to only use a single unique set of the (three) response functions to perform (3-tissue) spherical deconvolution of all subjects: as the (3-tissue) spherical deconvolution results will be expressed in function of this set of response functions, they can (in an abstract way) be seen as the units of both the final apparent fibre density metric and the other compartments estimated in the model. One possible way to obtain a unique set of response functions, is to average the response functions obtained from all subjects for each tissue type:

responsemean */response_wm.txt ../group_average_response_wm.txt
responsemean */response_gm.txt ../group_average_response_gm.txt
responsemean */response_csf.txt ../group_average_response_csf.txt

There is however no strict requirement for the final set of response functions to be the average of all subject response functions, for each tissue type (or indeed, it doesn’t even have to be the average per se). In certain very specific cases, it may even be wise to leave out subjects (for this step) where the response functions could not reliably be obtained, or where pathology affected the brain globally.

5. Upsampling DW images

Upsampling DWI data before computing FODs increases anatomical contrast and improves downstream template building, registration, tractography and statistics. We recommend upsampling to an isotropic voxel size of 1.25 mm for human brains (if your original resolution is already higher, you can skip this step):

for_each * : mrgrid IN/dwi_denoised_unringed_preproc_unbiased.mif regrid -vox 1.25 IN/dwi_denoised_unringed_preproc_unbiased_upsampled.mif

6. Compute upsampled brain mask images

Compute a whole brain mask from the upsampled DW images:

for_each * : dwi2mask IN/dwi_denoised_unringed_preproc_unbiased_upsampled.mif IN/dwi_mask_upsampled.mif

Warning

It is absolutely crucial to check at this stage that all individual subject masks include all regions of the brain that are intended to be analysed. Fibre orientation distributions will only be computed within these masks; and at a later step (in template space) the analysis mask will be restricted to the intersection of all masks, so any individual subject mask which excludes a certain region, will result in this region being excluded from the entire analysis (unless a more advanced pipeline is followed; see Mitigating the effects of brain cropping). Masks appearing too generous or otherwise including non-brain regions should generally not cause any concerns at this stage. Hence, if in doubt, it is advised to always err on the side of inclusion (of regions) at this stage.

Note

The earlier dwibiascorrect step is not fundamentally important in the multi-tissue fixel-based analysis pipeline, as the later mtnormalise step performs more robustly (and if dwibiascorrect is included, mtnormalise will later on typically improve the result further). While performing the earlier dwibiascorrect step typically improves dwi2mask performance, cases have been observed where the opposite is true (typically if the data contains only weak bias fields). If required, experiment by either including or excluding dwibiascorrect in the pipeline in function of the best dwi2mask outcome and manually correct the masks if necessary (by adding regions which dwi2mask fails to include).

7. Fibre Orientation Distribution estimation (multi-tissue spherical deconvolution)

When performing fixel-based analysis, multi-tissue constrained spherical deconvolution should be performed using the unique set of (average) tissue response functions obtained before:

for_each * : dwi2fod msmt_csd IN/dwi_denoised_unringed_preproc_unbiased_upsampled.mif ../group_average_response_wm.txt IN/wmfod.mif ../group_average_response_gm.txt IN/gm.mif  ../group_average_response_csf.txt IN/csf.mif -mask IN/dwi_mask_upsampled.mif

8. Joint bias field correction and intensity normalisation

To perform joint bias field correction and global intensity normalisation of the multi-tissue compartment parameters, use mtnormalise:

for_each * : mtnormalise IN/wmfod.mif IN/wmfod_norm.mif IN/gm.mif IN/gm_norm.mif IN/csf.mif IN/csf_norm.mif -mask IN/dwi_mask_upsampled.mif

If multi-tissue CSD was performed with the same single set of (three) tissue response functions for all subjects, then the resulting output of mtnormalise makes the absolute amplitudes comparable between those subjects as well. Note that this step is crucial in the FBA pipeline, even if bias field correction was applied earlier using dwibiascorrect, since dwibiascorrect does not correct for global intensity differences between subjects. The performance of mtnormalise is not significantly impacted by either having run dwibiascorrect before or not. In case prior bias field correction was run in the pipeline, mtnormalise will further correct for residual intensity inhomogeneities.

Warning

mtnormalise results can be sensitive to masks that contain non-brain voxels. The underlying algorithm will attempt to drive the sum of tissue volumes to unity in such voxels - despite not containing brain tissue - which can result in erroneous bias field correction if the number of such voxels is large. For this reason we recommend using conservative (i.e. less spatially extended) masks for the mtnormalise step. Unlike step 6, where inclusion of all brain voxels was encouraged even at the expense of including some non-brain voxels, for bias field estimation exclusion of non-brain voxels is of greater priority than inclusion of all brain voxels.

9. Generate a study-specific unbiased FOD template

Population template creation is one of the most time consuming steps in a fixel-based analysis. If you have a very large number of subjects in your study, you can opt to create the template from a limited subset of 30-40 individuals. Typically, subjects are chosen so the generated template is representative of your population (e.g. similar number of patients and controls, though avoid patients with excessive abnormalities compared to the rest of the population). To build a template, put all FOD images in a single folder and put a set of corresponding mask images (with the same prefix as the FOD images) in another folder (using masks speeds up registration significantly):

mkdir -p ../template/fod_input
mkdir ../template/mask_input

Symbolic link all FOD images (and masks) into a single input folder. To use the entire population to build the template:

for_each * : ln -sr IN/wmfod_norm.mif ../template/fod_input/PRE.mif
for_each * : ln -sr IN/dwi_mask_upsampled.mif ../template/mask_input/PRE.mif

If you opt to create the template from a limited subset of (e.g. 30-40) subjects and your study has multiple groups, then you can aim for a similar number of subjects from each group to make the template more representative of the population as a whole. Assuming the subject directory labels can be used to identify members of each group, you could use:

for_each `ls -d *patient | sort -R | tail -20` : ln -sr IN/wmfod_norm.mif ../template/fod_input/PRE.mif ";" ln -sr IN/dwi_mask_upsampled.mif ../template/mask_input/PRE.mif
for_each `ls -d *control | sort -R | tail -20` : ln -sr IN/wmfod_norm.mif ../template/fod_input/PRE.mif ";" ln -sr IN/dwi_mask_upsampled.mif ../template/mask_input/PRE.mif

Run the template building script as follows:

population_template ../template/fod_input -mask_dir ../template/mask_input ../template/wmfod_template.mif -voxel_size 1.25

The voxel size is typically set to match the voxel size of the input FOD images (which, in this pipeline, would typically be the resolution the preprocessed data was upsampled to, earlier on in the pipeline).

10. Register all subject FOD images to the FOD template

Register the FOD image from each subject to the FOD template:

for_each * : mrregister IN/wmfod_norm.mif -mask1 IN/dwi_mask_upsampled.mif ../template/wmfod_template.mif -nl_warp IN/subject2template_warp.mif IN/template2subject_warp.mif

11. Compute the template mask (intersection of all subject masks in template space)

Different subjects have different brain coverage. To ensure subsequent analysis is performed in voxels that contain data from all subjects, we warp all subject masks into template space and compute the template mask as the intersection of all subject masks in template space. To warp all masks into template space:

for_each * : mrtransform IN/dwi_mask_upsampled.mif -warp IN/subject2template_warp.mif -interp nearest -datatype bit IN/dwi_mask_in_template_space.mif

Compute the template mask as the intersection of all warped masks:

mrmath */dwi_mask_in_template_space.mif min ../template/template_mask.mif -datatype bit

Warning

It is absolutely crucial to check at this stage that the resulting template mask includes all regions of the brain that are intended to be analysed. If this is not the case, the cause will be either an individual subject mask which did not include a certain region, or the template building process or individual subject registrations having gone wrong for one or more subjects. It is advised to go back to these steps, and identify and resolve the cause of the problem before continuing any further.

Note

It is possible at this stage to edit the template mask and remove certain regions that are either not part of the brain (though these are unlikely to survive the intersection step to begin with, as they would have to be present in all subject masks for this to happen), or otherwise not of interest. However, do not add regions to the mask at this stage. If there is a genuine need to do this, go back to the relevant steps which caused the exclusion of these regions to begin with (see also the above warning).

12. Compute a white matter template analysis fixel mask

In this step, we segment fixels from the FOD template. The result is the fixel mask that defines the fixels for which statistical analysis will later on be performed (and hence also which fixels’ statistics can support others via the mechanism of connectivity-based fixel enhancement (CFE) [Raffelt2015]):

fod2fixel -mask ../template/template_mask.mif -fmls_peak_value 0.06 ../template/wmfod_template.mif ../template/fixel_mask

Note

Fixel images, which appear in the pipeline from this step onwards, are stored using the Fixel image (directory) format, which stores all fixel data for a fixel image in a directory (i.e. a folder).

Warning

This step ultimately determines the fixel mask in which statistical analysis will be performed, and hence also which fixels’ statistics can contribute to others via the CFE mechanism; so it may have a substantial impact on the final result. Essentially, it can be detrimental to the result if the threshold value specified via the -fmls_peak_value is too high and hence excludes genuine white matter fixels. This risk is substantially higher in voxels containing crossing fibres (and higher the more fibres are crossing in a single voxel). Even though 0.06 has been observed to be a decent default value for 3-tissue CSD population templates, it is still strongly advised to visualise the output fixel mask using mrview. Do this by opening the index.mif found in ../template/fixel_mask via the fixel plot tool. If, with respect to known or normal anatomy, fixels are missing (especially paying attention to crossing areas), regenerate the mask with a lower value supplied to the -fmls_peak_value option (of course, avoid lowering it too much, as too many false or noisy fixels may be introduced). For an adult human brain template, and using an isotropic template voxel size of 1.25 mm, it is expected to have several hundreds of thousands of fixels in the fixel mask (you can check this by mrinfo -size ../template/fixel_mask/directions.mif, and looking at the size of the image along the first dimension).

13. Warp FOD images to template space

Note that here we warp FOD images into template space without FOD reorientation, as reorientation will be performed in a separate subsequent step (after fixel segmentation):

for_each * : mrtransform IN/wmfod_norm.mif -warp IN/subject2template_warp.mif -reorient_fod no IN/fod_in_template_space_NOT_REORIENTED.mif

14. Segment FOD images to estimate fixels and their apparent fibre density (FD)

Here we segment each FOD lobe to identify the number and orientation of fixels in each voxel. The output also contains the apparent fibre density (AFD) value per fixel (estimated as the FOD lobe integral):

for_each * : fod2fixel -mask ../template/template_mask.mif IN/fod_in_template_space_NOT_REORIENTED.mif IN/fixel_in_template_space_NOT_REORIENTED -afd fd.mif

Note that in the following steps we will use the more generic shortened acronym Fibre Density (FD) to refer to the AFD metric.

15. Reorient fixels

Here we reorient the fixels of all subjects in template space based on the local transformation at each voxel in the warps used previously:

for_each * : fixelreorient IN/fixel_in_template_space_NOT_REORIENTED IN/subject2template_warp.mif IN/fixel_in_template_space

After this step, the fixel_in_template_space_NOT_REORIENTED folders can be safely removed.

16. Assign subject fixels to template fixels

While each subject’s data has already been (spatially) warped to the common template space, and subject fixels have been reoriented accordingly, there is still no specification of which fixels match (across subjects, and between the subject and template fixels). This step establishes exactly that, by matching the fixels of each individual subject to the single common set of template fixels (which then inherently also defines how they match across subjects). This is achieved by, for each fixel in the template fixel mask, identifying the corresponding fixel in the matching voxel of the subject image and assigning the FD value of this corresponding subject fixel to that fixel in template space. If no fixel exists or can be found in a subject that corresponds to a given template fixel then it is assigned a value of zero (as the absence of a subject fixel at this stage is most likely due to a very low, or even zero, FD). This step is performed as follows:

for_each * : fixelcorrespondence IN/fixel_in_template_space/fd.mif ../template/fixel_mask ../template/fd PRE.mif

Note that the output fixel directory ../template/fd is the same for all subjects. This makes sense, since after this operation, there is only a single remaining set of fixels (i.e. the template fixels), with corresponding FD values as obtained from each subject. This resulting directory ../template/fd now stores these data as individual fixel data files: one for each subject, and all with respect to a single set of corresponding template fixels. This way of storing the entire population’s FD data is then ready for input to fixelcfestats later on.

17. Compute the fibre cross-section (FC) metric

The fibre density metric, mapped directly without any modulation to the fixel template space as above, is only sensitive to the original density of intra-axonal space in each voxel. In other words, it ignores the cross-sectional size of the bundle, which is another property that would factor into the bundle’s total intra-axonal space across its full cross-sectional extent, and hence influence its total capacity to carry information. In certain cases, for example, atrophy may impact this cross-sectional size, but not per se the local fibre density metric.

In this step, we compute a fixel-based metric related to morphological differences in fibre cross-section (FC), where information is derived entirely from the warps generated during registration (see [Raffelt2017] for more information):

for_each * : warp2metric IN/subject2template_warp.mif -fc ../template/fixel_mask ../template/fc IN.mif

However, for group statistical analysis of FC we recommend calculating the log(FC) to ensure data are centred around zero and normally distributed. Here, we create a separate fixel directory to store the log(FC) data and copy the fixel index and directions file across:

mkdir ../template/log_fc
cp ../template/fc/index.mif ../template/fc/directions.mif ../template/log_fc
for_each * : mrcalc ../template/fc/IN.mif -log ../template/log_fc/IN.mif

Note

The FC (and hence also the log(FC)) as calculated here, is a relative metric, expressing the local fixel-wise cross-sectional size relative to this study’s population template. While this makes it possible to interpret differences of FC within a single study (because only a single unique template is used in the study), the FC values should not be compared across different studies that each have their own population template. Reporting absolute quantities of FC, or absolute effect sizes of FC, also provides little information; as again, it is only meaningful with respect to the template.

18. Compute a combined measure of fibre density and cross-section (FDC)

The total capacity of a fibre bundle to carry information, is modulated both by the local fibre density at the voxel (fixel) level, as well as its cross-sectional size. Here we compute a combined metric, which factors in the effects of both FD and FC, resulting in a fibre density and cross-section (FDC) metric:

mkdir ../template/fdc
cp ../template/fc/index.mif ../template/fdc
cp ../template/fc/directions.mif ../template/fdc
for_each * : mrcalc ../template/fd/IN.mif ../template/fc/IN.mif -mult ../template/fdc/IN.mif

This is also a nice example of how calculations across multiple fixel data files can be performed. However, note that this is only valid if these both share the same set of original fixels (in this case, the template fixel mask). Because the fixels also have to be stored in exactly the same order for this to work correctly, great care has to be taken the index.mif files (in this case of the ../template/fd and ../template/fc folders) related to all input fixel data files that are used in the mrcalc command are exact copies of each other.

19. Perform whole-brain fibre tractography on the FOD template

Statistical analysis using connectivity-based fixel enhancement (CFE) [Raffelt2015] exploits local connectivity information derived from probabilistic fibre tractography, which acts as a neighbourhood definition for threshold-free enhancement of locally clustered statistic values. To generate a whole-brain tractogram from the FOD template (note the remaining steps from here on are executed from the template directory):

cd ../template
tckgen -angle 22.5 -maxlen 250 -minlen 10 -power 1.0 wmfod_template.mif -seed_image template_mask.mif -mask template_mask.mif -select 20000000 -cutoff 0.06 tracks_20_million.tck

Warning

The appropriate FOD amplitude cutoff for FOD template tractography can vary considerably between different datasets, as well as different versions of MRtrix3 due to historical software bugs. While the value of 0.06 is suggested as a reasonable value for multi-tissue data, it may be beneficial to first generate a smaller number of streamlines (e.g. 100,000) using this value, and visually confirm that the generated streamlines exhibit an appropriate extent of propagation at the ends of white matter pathways, before committing to generation of the dense tractogram.

20. Reduce biases in tractogram densities

Perform SIFT to reduce tractography biases in the whole-brain tractogram:

tcksift tracks_20_million.tck wmfod_template.mif tracks_2_million_sift.tck -term_number 2000000

21. Generate fixel-fixel connectivity matrix

Generation of the fixel-fixel connectivity matrix based on the whole-brain streamlines tractogram is performed as follows:

fixelconnectivity fixel_mask/ tracks_2_million_sift.tck matrix/

The output directory should contain three images: index.mif, fixels.mif and values.mif; these are used to encode the fixel-fixel connectivity that is by its nature sparse.

Warning

Running fixelconnectivity requires quite a lot of memory (for about 500,000 fixels in the template analysis fixel mask and a typical tractogram defining the pairwise connectivity between fixels, 32GB of RAM is a typical memory requirement; you can check the number of fixels in your template analysis fixel mask by mrinfo -size ./fixel_template/directions.mif, and looking at the size of the image along the first dimension; also check for (and avoid) gross false positive connections in the tractogram, e.g. streamlines connecting both hemispheres in large areas where they shouldn’t). Should hardware related limitations force you to reduce the number of fixels, it is not advised to change the threshold to derive the white matter template analysis fixel mask, as it may remove crucial fixels deep in the white matter (e.g. in crossing areas). Rather, consider slightly increasing the template voxel size, or spatially removing regions that are not of interest from the template mask (which was earlier obtained in this pipeline as the intersections of all subject masks in template space).

22. Smooth fixel data using fixel-fixel connectivity

Smoothing of fixel data is performed based on the sparse fixel-fixel connectivity matrix:

fixelfilter fd smooth fd_smooth -matrix matrix/
fixelfilter log_fc smooth log_fc_smooth -matrix matrix/
fixelfilter fdc smooth fdc_smooth -matrix matrix/

By calling the fixelfilter command on each fixel directory in turn, the smoothing filter will be appplied to all fixel data files present in each of those directories; it is not necessary to call the fixelfilter command separately for each individual fixel data file.

23. Perform statistical analysis of FD, FC, and FDC

Statistical analysis using CFE is performed separately for each metric (FD, log(FC), and FDC) as follows:

fixelcfestats fd_smooth/ files.txt design_matrix.txt contrast_matrix.txt matrix/ stats_fd/
fixelcfestats log_fc_smooth/ files.txt design_matrix.txt contrast_matrix.txt matrix/ stats_log_fc/
fixelcfestats fdc_smooth/ files.txt design_matrix.txt contrast_matrix.txt matrix/ stats_fdc/

The input files.txt is a text file containing the filename of each file (i.e. not the full path) to be analysed inside the input fixel directory, each filename on a separate line. The line ordering should correspond to the lines in the file design_matrix.txt.

Note

While previous versions of fixelcfestats performed smoothing of the input fixel data as those data were imported, this is no longer the case. It is expected that the fixel data provided as input to fixelcfestats will already have been appropriately smoothed; e.g. using fixelfilter.

Note

Unlike some other software packages providing a GLM (e.g. FSL randomise), a column of 1’s (corresponding to the “global intercept”: the mean image value when all other design matrix factors are zero) will not be automatically included. If the model for your experiment requires such a factor, it is necessary to include it explicitly in the design matrix.

24. Visualise the results

To view the results load the population FOD template image in mrview, and overlay the fixel images using the vector plot tool. Note that p-value images are saved as (1 - p-value). Therefore to visualise all results at a threshold of p < 0.05, within the mrview fixel plot tool, apply a lower threshold at a value of 0.95.