Fibre density and cross-section - Single-tissue CSD¶
This tutorial explains how to perform fixel-based analysis of fibre density and cross-section [Raffelt2017] using single-tissue spherical deconvolution. 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 more details about single-tissue AFD.
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_patient/dwi.mif study/subjects/001_patient/wmfod.mif study/subjects/002_control/dwi.mif study/subjects/002_control/wmfod.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 dwipreproc 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 dwipreproc help page).
The simplest scenario is to (only) correct for motion and eddy current-induced distortions:
foreach * : dwipreproc IN/dwi_denoised_unringed.mif IN/dwi_denoised_unringed_preproc.mif -rpe_none -pe_dir AP
-pe_dir option to dwipreproc 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 dwipreproc section and the dwipreproc documentation.
3. Estimate a temporary brain mask¶
Compute a brain mask:
foreach * : dwi2mask IN/dwi_denoised_preproc.mif IN/dwi_temp_mask.mif
AFD-specific pre-processsing steps¶
To enable robust quantitative comparisons of AFD across subjects three additional steps are required. Note these can potentially be skipped if analysing certain other DWI fixel-based measures related to fibre density (for example CHARMED).
4. Bias field correction¶
Because we recommend a global intensity normalisation, bias field correction is required as a pre-processing step to eliminate low frequency intensity inhomogeneities across the image. DWI bias field correction is perfomed by first estimating the bias field from the DWI b=0 data, then applying the field to correct all DW volumes. This can be done in a single step using the dwibiascorrect script with the
-ants option in MRtrix. The script uses a bias field correction algorithm available in ANTS (the N4 algorithm). Don’t use the
-fsl option with this script in this fixel-based analysis pipeline. To perform bias field correction on DW images, run:
foreach * : dwibiascorrect -ants IN/IN/dwi_denoised_unringed_preproc.mif IN/IN/dwi_denoised_unringed_preproc_unbiased.mif
5. Global intensity normalisation across subjects¶
As outlined here, a global intensity normalisation is required for AFD analysis. For a single-tissue pipeline, a possible approach is to use the dwiintensitynorm script. The script performs normalisation on all subjects within a study (using a group-wise registration), and therefore the input and output arguments are directories containing all study images. First create directories to store all the input and output images. From the subjects directory:
mkdir -p ../dwiintensitynorm/dwi_input mkdir ../dwiintensitynorm/mask_input
You could copy all files into this directory, however symbolic linking them will save space:
foreach * : ln -sr IN/dwi_denoised_unringed_preproc_unbiased.mif ../dwiintensitynorm/dwi_input/IN.mif foreach * : ln -sr IN/dwi_temp_mask.mif ../dwiintensitynorm/mask_input/IN.mif
Perform intensity normalisation:
dwiintensitynorm ../dwiintensitynorm/dwi_input/ ../dwiintensitynorm/mask_input/ ../dwiintensitynorm/dwi_output/ ../dwiintensitynorm/fa_template.mif ../dwiintensitynorm/fa_template_wm_mask.mif
Link the output files back to the subject directories:
foreach ../dwiintensitynorm/dwi_output/* : ln -sr IN PRE/dwi_denoised_unringed_preproc_unbiased_normalised.mif
The dwiintensitynorm script also outputs the study-specific FA template and white matter mask. It is recommended that you check that the white matter mask is appropriate (i.e. does not contain CSF or voxels external to the brain. It needs to be a rough WM mask). If you feel the white matter mask needs to be larger or smaller you can re-run
dwiintensitynorm with a different
-fa_threshold option. Note that if your input brain masks include CSF then this can cause spurious high FA values outside the brain which are then included in the template white matter mask.
Keeping the FA template image and white matter mask is also handy if additional subjects are added to the study at a later date. New subjects can be intensity normalised in a single step by piping the following commands together. Run from the subjects directory:
dwi2tensor new_subject/dwi_denoised_unringed_preproc_unbiased.mif -mask new_subject/dwi_temp_mask.mif - | tensor2metric - -fa - | mrregister -force ../dwiintensitynorm/fa_template.mif - -mask2 new_subject/dwi_temp_mask.mif -nl_scale 0.5,0.75,1.0 -nl_niter 5,5,15 -nl_warp - /tmp/dummy_file.mif | mrtransform ../dwiintensitynorm/fa_template_wm_mask.mif -template new_subject/dwi_denoised_unringed_preproc_unbiased.mif -warp - - | dwinormalise new_subject/dwi_denoised_unringed_preproc_unbiased.mif - ../dwiintensitynorm/dwi_output/new_subject.mif
The above command may also be useful if you wish to alter the mask and re-apply the intensity normalisation to all subjects in the study. For example, you may wish to edit the mask using the ROI tool in
mrview to remove white matter regions that you hypothesise are affected by the disease (e.g. removing the corticospinal tract in a study of motor neurone disease due to T2 hyperintensity). You also may wish to redefine the mask entirely, for example in an elderly population (with enlarged ventricles) it may be possible, or even preferable, to normalise using the median b=0 CSF. This could be performed by manually masking partial-volume-free CSF voxels, then running the above command with the CSF mask instead of the
We strongly recommend you that you check the scale factors applied during intensity normalisation are not influenced by the variable of interest in your study. For example if one group contains global changes in white matter T2 then this may directly influence the intensity normalisation and therefore bias downstream AFD analysis. To check this, you can perform an equivalence test to ensure mean scale factors are the same between groups. To output the scale factor applied for all subjects use
mrinfo ../dwiintensitynorm/dwi_output/* -property dwi_norm_scale_factor.
Fixel-based analysis steps¶
6. Computing an (average) white matter response function¶
A robust and fully automated (unsupervised) method to obtain single-shell response functions representing single-fibre white matter from your data, is the approach proposed in [Tournier2013], which can be run by:
foreach * : dwi2response tournier IN/dwi_denoised_unringed_preproc_unbiased_normalised.mif IN/response.txt
It is crucial for fixel-based analysis to only use a single unique response function to perform spherical deconvolution of all subjects: as all resulting fibre orientation distributions will be expressed in function of it, it can (in an abstract way) be seen as the unit of the final apparent fibre density metric. A possible way to obtain a unique response function, is to average the response functions obtained from all subjects:
average_response */response.txt ../group_average_response.txt
There is however no strict requirement for the (one) final response function to be the average of all subject response functions. In certain very specific cases, it may even be wise to leave out subjects (for this step) where a response function could not reliably be obtained, or where pathology affected the brain globally.
7. Upsampling DW images¶
Upsampling DWI data before computing FODs can increase anatomical contrast and improve 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_normalised.mif -vox 1.3 IN/dwi_denoised_unringed_preproc_unbiased_normalised_upsampled.mif
8. Compute upsampled brain mask images¶
Compute a whole brain mask from the upsampled DW images:
foreach * : dwi2mask IN/dwi_denoised_unringed_preproc_unbiased_normalised_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. Manually correct the masks if necessary.
9. Fibre Orientation Distribution estimation (spherical deconvolution)¶
When performing fixel-based analysis, constrained spherical deconvolution (CSD) should be performed using the unique (average) white matter response function obtained before. Note that
dwi2fod csd can be used, however here we use
dwi2fod msmt_csd (even with single shell data) to benefit from the hard non-negativity constraint, which has been observed to lead to more robust outcomes:
foreach * : dwiextract IN/dwi_denoised_unringed_preproc_unbiased_normalised_upsampled.mif - \| dwi2fod msmt_csd - ../group_average_response.txt IN/wmfod.mif -mask IN/dwi_mask_upsampled.mif
10. 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.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.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.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).
11. Register all subject FOD images to the FOD template¶
Register the FOD image from each subject to the FOD template:
foreach * : mrregister IN/wmfod.mif -mask1 IN/dwi_mask_upsampled.mif ../template/wmfod_template.mif -nl_warp IN/subject2template_warp.mif IN/template2subject_warp.mif
12. 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).
13. Compute a white matter template analysis fixel mask¶
Here we perform a 2-step threshold to identify template white matter fixels to be included in the analysis. Fixels in the template fixel analysis mask are also used to identify the best fixel correspondence across all subjects (i.e. match fixels across subjects within a voxel).
Compute a template AFD peaks fixel image:
fod2fixel ../template/wmfod_template.mif -mask ../template/template_mask.mif ../template/fixel_temp -peak peaks.mif
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).
Next view the peaks file using the fixel plot tool in mrview and identify an appropriate threshold that removes peaks from grey matter, yet does not introduce any large gaps in your white matter (e.g. 0.25; this should be tuned manually to be suitable for your data).
Threshold the peaks fixel image:
mrthreshold ../template/fixel_temp/peaks.mif -abs 0.25 ../template/fixel_temp/mask.mif
Generate an analysis voxel mask from the fixel mask. The median filter in this step should remove spurious voxels outside the brain, and fill in the holes in deep white matter where you have small peaks due to 3-fibre crossings:
fixel2voxel ../template/fixel_temp/mask.mif max - | mrfilter - median ../template/voxel_mask.mif rm -rf ../template/fixel_temp
Recompute the fixel mask using the analysis voxel mask. Using the mask allows us to use a lower AFD threshold than possible in the steps above, to ensure we have included fixels with low AFD inside white matter (e.g. regions with fibre crossings):
fod2fixel -mask ../template/voxel_mask.mif -fmls_peak_value 0.1 ../template/wmfod_template.mif ../template/fixel_mask
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.1 has been observed to be a reasonable value for single-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. At the same time, try to avoid introducing too many noisy fixels, especially in gray matter areas. For single-tissue CSD, this balance is however not always easy to find. 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).
14. 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.mif -warp IN/subject2template_warp.mif -noreorientation IN/fod_in_template_space_NOT_REORIENTED.mif
15. 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.
16. 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.
17. 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.
18. 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.
19. 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.
20. 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 tracks_20_million.tck
21. Reduce biases in tractogram densities¶
Perform SIFT to reduce tractography biases in the whole-brain tractogram:
tcksift tracks_20_million.tck fod_template.mif tracks_2_million_sift.tck -term_number 2000000
22. 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 input_tracks_2_million_sift.tck stats_fd fixelcfestats log_fc files.txt design_matrix.txt contrast_matrix.txt input_tracks_2_million_sift.tck stats_log_fc fixelcfestats fdc files.txt design_matrix.txt contrast_matrix.txt input_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, 128GB of RAM is required; you can check the number of fixels in your template analysis fixel mask by
mrinfo -size ../template/fixel_template/mask.mif, and looking at the size of the image along the first dimension). 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). Refer to this section for more information: Commands crashing due to memory requirements.
23. 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.