Fibre density and cross-section - Single-tissue CSD

Introduction

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_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:

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

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 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

Warning

The -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

Note

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 fa_template_wm_mask.mif.

Warning

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

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. 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

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).

13. 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.10 ../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.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

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.

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/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.

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 -cutoff 0.10 tracks_20_million.tck

Warning

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.10 option in the command line above, in line with the instructions for older versions of MRtrix3.

21. 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

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 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

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

For correlation analysis, a column of 1’s will not be automatically included (as per FSL randomise).

Note

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.

Warning

Running 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). 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.