Fibre density and cross-section - Multi-tissue CSD¶
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 foreach 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
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
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:
foreach * : dwidenoise IN/dwi.mif IN/dwi_denoised.mif
If Gibbs ringing removal is included, it follows immediately after:
foreach * : mrdegibbs IN/dwi_denoised.mif IN/dwi_denoised_unringed.mif -axes 0,1
-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:
foreach * : dwifslpreproc IN/dwi_denoised_unringed.mif IN/dwi_denoised_unringed_preproc.mif -rpe_none -pe_dir AP
-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:
foreach * : 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:
foreach * : 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.3 mm for human brains (if your original resolution is already higher, you can skip this step):
foreach * : mrresize IN/dwi_denoised_unringed_preproc_unbiased.mif -vox 1.3 IN/dwi_denoised_unringed_preproc_unbiased_upsampled.mif
6. Compute upsampled brain mask images¶
Compute a whole brain mask from the upsampled DW images:
foreach * : dwi2mask IN/dwi_denoised_unringed_preproc_unbiased_upsampled.mif IN/dwi_mask_upsampled.mif
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. 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.
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:
foreach * : 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 (in the log-domain), use mtnormalise:
foreach * : 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.
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:
foreach * : ln -sr IN/wmfod_norm.mif ../template/fod_input/PRE.mif foreach * : 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:
foreach `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 foreach `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.3
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:
foreach * : 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:
foreach * : 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
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.
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
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).
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.3 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):
foreach * : mrtransform IN/wmfod_norm.mif -warp IN/subject2template_warp.mif -noreorientation 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):
foreach * : 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:
foreach * : 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:
foreach * : 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):
foreach * : 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 foreach * : mrcalc ../template/fc/IN.mif -log ../template/log_fc/IN.mif
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 foreach * : 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/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
The command line above assumes you’re working with MRtrix3 RC3 or above. An important bug in the tractography code was fixed in that version of the software. If you are not able to update your installation, and are still working with an older version of MRtrix3, you should remove the -cutoff 0.06 option in the command line above, in line with the instructions for older versions of MRtrix3.
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. 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 files.txt design_matrix.txt contrast_matrix.txt tracks_2_million_sift.tck stats_fd fixelcfestats log_fc files.txt design_matrix.txt contrast_matrix.txt tracks_2_million_sift.tck stats_log_fc fixelcfestats fdc files.txt design_matrix.txt contrast_matrix.txt tracks_2_million_sift.tck stats_fdc
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
For correlation analysis, a column of 1’s will not be automatically included (as per FSL randomise).
fixelcfestats currently only accepts a single contrast. However if the opposite (negative) contrast is also required (i.e. a two-tailed test), then use the
-neg option. Several output files will generated all starting with the supplied prefix.
fixelcfestats requires 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, 128GB 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 (and avoid) for 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. 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 p-values < 0.05, threshold the fixels using the fixel plot tool at 0.95.