Fibre density and cross-section - Multi-tissue CSD

Introduction

This tutorial explains how to perform fixel-based analysis of fibre density and cross-section with fibre orientation distributions (FODs) computing using multi-tissue CSD using single-shell data or multi-shell data. We note that high b-value (>2000s/mm2) data is recommended to aid the interpretation of AFD being related to the intra-axonal space. See the original paper for more details.

All steps in this tutorial have 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 the subject, and all files within identifying the image type. 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

Note

All commands in this tutorial are run from the subjects path up until step 19, where 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. Please post any questions or issues on the MRtrix community forum.

Pre-processsing steps

1. DWI denoising

The effective SNR of diffusion data can be improved considerably by exploiting the redundancy in the data to reduce the effects of thermal noise. This functionality is provided in the command dwidenoise:

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

Note that this denoising step must be performed prior to any other image pre-processing: any form of image interpolation (e.g. re-gridding images following motion correction) will invalidate the statistical properties of the image data that are exploited by dwidenoise, and make the denoising process prone to errors. Therefore this process is applied as the very first step.

2. DWI general pre-processing

The dwipreproc script is provided for performing general pre-processing of diffusion image data - this includes eddy current-induced distortion correction, motion correction, and (possibly) susceptibility-induced distortion correction. Commands for performing this pre-processing are not yet implemented in MRtrix3; the dwipreproc script in its current form is in fact an interface to the relevant commands that are provided as part of the FSL package. Installation of FSL (including eddy) is therefore required to use this script, and citation of the relevant articles is also required (see the dwipreproc help page).

Usage of this script varies depending on the specific nature of the DWI acquisition with respect to EPI phase encoding - full details are available within the DWI distortion correction using dwipreproc page, and the dwipreproc help file.

Here, only a simple example is provided, where a single DWI series is acquired where all volumes have an anterior-posterior (A>>P) phase encoding direction:

foreach * : dwipreproc IN/dwi_denoised.mif IN/dwi_denoised_preproc.mif -rpe_none -pe_dir AP

3. Bias field correction

Bias field correction is important to deal with spatial intensity inhomogeneities. Even though this FBA pipeline will account for these as well (in the later mtnormalise step, which is furthermore crucial to correct for global intensity differences between subjects), performing bias field correction at this stage will allow for more accurate estimation of the tissue response functions as well as the individual subject brain masks.

This can be done in a single step using the dwibiascorrect script in MRtrix. The script uses bias field correction algorthims available in ANTS or FSL. In our experience the N4 algorithm in ANTS gives far superior results. To install N4, install the ANTS package, then perform bias field correction on DW images using:

foreach * : dwibiascorrect -ants IN/dwi_denoised_preproc.mif IN/dwi_denoised_preproc_bias.mif

Fixel-based analysis steps

4. Computing group average tissue response functions

As described here, using the same response function when estimating FOD images for all subjects enables differences in the intra-axonal volume (and therefore DW signal) across subjects to be detected as differences in the FOD amplitude (the AFD). To ensure the response function is representative of your study population, a group average response function can be computed by first estimating a response function per subject, then averaging with the script:

foreach * : dwi2response dhollander IN/dwi_denoised_preproc_bias.mif IN/response_wm.txt IN/response_gm.txt IN/response_csf.txt
average_response */response_wm.txt ../group_average_response_wm.txt
average_response */response_gm.txt ../group_average_response_gm.txt
average_response */response_csf.txt ../group_average_response_csf.txt

5. Upsampling DW images

Upsampling DWI data before computing FODs can increase anatomical contrast and improve downstream spatial normalisation and statistics. We recommend upsampling to a voxel size of 1.25mm (for human brains). If you have data that has smaller voxels than 1.25mm, then we recommend you can skip this step:

foreach * : mrresize IN/dwi_denoised_preproc_bias.mif -vox 1.25 IN/dwi_denoised_preproc_bias_upsampled.mif

6. Compute upsampled brain mask images

Compute a whole brain mask from the upsampled DW images:

foreach * : dwi2mask IN/dwi_denoised_preproc_bias_upsampled.mif IN/dwi_mask_upsampled.mif

7. Fibre Orientation Distribution estimation

When performing analysis of AFD, Constrained Spherical Deconvolution (CSD) should be performed using the group average response functions computed at step 3:

foreach * : dwi2fod msmt_csd IN/dwi_denoised_preproc_bias_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. Intensity normalisation

This step performs global intensity normalisation in the log-domain by scaling all tissue types with a spatially smoothly varying normalisation field:

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 CSD was performed with the same single set of (average) WM, GM and CSF response functions for all subjects, then the resulting output of mtnormalise should make the amplitudes comparable between those subjects as well.

Note that this step is crucial in the FBA pipeline, even if bias field correction was applied during the preprocessing stage, as the latter does not correct for global intensity differences between subjects.

Warning

We recommend you that you check that the normalisation scale (computed during intensity normalisation) is not influenced by the variable of interest in your study. For example if one group contains global (widespread) changes in white matter T2, then this may directly influence the intensity normalisation and therefore bias downstream analysis of apparent fibre density (FD). To check this, you can perform an equivalence test to ensure the overall normalisation scale does not differ between groups. To output these overall normalisation scales for all subjects use mrinfo */wmfod_norm.mif -property lognorm_scale.

9. Generate a study-specific unbiased FOD template

Population template creation is a very time consuming step in a fixel-based analysis. If you have a very large number of subjects in your study, we recommend building the template from a limited subset of 30-40 individuals. Subjects should be chosen to ensure the generated template is representative of your population (e.g. similar number of patients and controls). To build a template, place all FOD images in a single folder. We also highly recommend placing a set of corresponding mask images (with the same prefix as the FOD images) in another folder. Using masks can speed up registration significantly as well as result in a much more accurate template in specific scenarios:

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

Symbolic link all FOD images (and masks) into a single input folder. If you have fewer than 40 subjects in your study, you can 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

Alternatively, if you have more than 40 subjects you can randomly select a subset of the individuals. If your study has multiple groups, then ideally you want to select the same number of subjects from each group to ensure the template is un-biased. 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.25

If you are building a template from your entire study population, run the population_template script use the -warp_dir warps option to output a directory containing all subject warps to the template. Saving the warps here will enable you to skip the next step. Note that the warps used (and therefore output) from the population_template script are 5D images containing both forward and reverse warps (see mrregister for more info). After population template creation is complete, to convert this warp format to a more conventional 4D deformation field format ready for the subsequent steps, run

$ foreach ../template/warps/* : warpconvert -type warpfull2deformation -template ../template/wmfod_template.mif IN PRE/subject2template_warp.mif

10. Register all subject FOD images to the FOD template

Register the FOD image from all subjects to the FOD template image. Note you can skip this step if you built your template from your entire population and saved the warps (see previous step):

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 intersection of all subject masks in template space

Different subjects will have subtly 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 mask intersection. For each subject:

foreach * : mrtransform IN/dwi_mask_upsampled.mif -warp IN/subject2template_warp.mif -interp nearest IN/dwi_mask_in_template_space.mif

Compute the intersection of all warped masks:

mrmath */dwi_mask_in_template_space.mif min ../template/mask_intersection.mif

12. Compute a white matter analysis voxel & fixel mask

Here we first identify all voxels having some white matter by thresholding the DC term (first SH coefficient) of the multi-tissue FOD image:

mrconvert ../template/wmfod_template.mif -coord 3 0 - | mrthreshold - ../template/voxel_mask.mif

Next we segment all fixels from each FOD in the template image (see here for more information about a analysis fixel mask). Note that the fixel image output from this step is stored using the Fixel image (directory) format, which exploits the filesystem to store all fixel data in a directory:

fod2fixel -mask ../template/voxel_mask.mif -fmls_peak_value 0.1 ../template/wmfod_template.mif ../template/fixel_mask

You can visualise the output fixels using the fixel plot tool from mrview, and opening either the index.mif or directions.mif found in ../template/fixel_mask. The automatic thresholding step used above should give you a mask that nicely covers all of white matter, however if not you can always try manually adjusting the threshold with the mrthreshold -abs option.

Note

We recommend having no more than 500,000 fixels in the analysis_fixel_mask (you can check this by mrinfo -size ../template/fixel_mask/directions.mif, and looking at the size of the image along the 1st dimension), otherwise downstream statistical analysis (using fixelcfestats) may run out of RAM). A mask with 500,000 fixels will require a PC with 128GB of RAM for the statistical analysis step. To reduce the number of fixels, try either reducing the number of voxels in the voxel mask by applying a manual threshold using -abs, increasing the -fmls_peak_value, or reducing the extent of upsampling in step 4.

13. Warp FOD images to template space

Note that here we warp FOD images into template space without FOD reorientation. Reorientation will be performed in a separate subsequent step:

foreach * : mrtransform IN/wmfod_norm.mif -warp IN/subject2template_warp.mif -noreorientation IN/fod_in_template_space.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 (see here for details on FOD segmentation). Note that in the following steps we will use a more generic shortened acronym - Fibre Density (FD) instead of AFD, since the following steps can also apply for other measures of fibre density (see the note below). The terminology is also consistent with our recent work:

foreach * : fod2fixel IN/fod_in_template_space.mif -mask ../template/voxel_mask.mif IN/fixel_in_template_space -afd fd.mif

15. Reorient fixel orientations

Here we reorient the direction of all fixels based on the Jacobian matrix (local affine transformation) at each voxel in the warp. Note that in-place fixel reorientation can be performed by specifing the output fixel folder to be the same as the input, and using the -force option:

foreach * : fixelreorient IN/fixel_in_template_space IN/subject2template_warp.mif IN/fixel_in_template_space --force

16. Assign subject fixels to template fixels

In step 8 & 9 we obtained spatial correspondence between subject and template. In step 14 we corrected the fixel orientations to ensure angular correspondence of the segmented peaks of subject and template. Here, for each fixel in the template fixel analysis mask, we identify the corresponding fixel in each voxel of the subject image and assign the FD value of the subject fixel to the corresponding fixel in template space. If no fixel exists in the subject that corresponds to the template fixel then it is assigned a value of zero. See this paper for more information. In the command below, you will note that the output fixel directory is the same for all subjects. This directory now stores data for all subjects at corresponding fixels, ready for input to fixelcfestats in step 20 below:

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

17. Compute fibre cross-section (FC) metric

Apparent fibre density, and other related measures that are influenced by the quantity of restricted water, only permit the investigation of group differences in the number of axons that manifest as a change to within-voxel density. However, depending on the disease type and stage, changes to the number of axons may also manifest as macroscopic differences in brain morphology. This step computes a fixel-based metric related to morphological differences in fibre cross-section, where information is derived entirely from the warps generated during registration (see this paper for more information):

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

Note that the FC files will be used in the next step. However, for group statistical analysis of FC we recommend taking the log(FC) to ensure data are centred about zero and normally distributed. We could place all the log(FC) fixel data files in the same fixel directory as the FC files (as long as they are named differently. However to keep things tidy, 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

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

To account for changes to both within-voxel fibre density and macroscopic atrophy, fibre density and fibre cross-section must be combined (a measure we call fibre density & cross-section, FDC). This enables a more complete picture of group differences in white matter. Note that as discussed in this paper, group differences in FD or FC alone must be interpreted with care in crossing-fibre regions. However group differences in FDC are more directly interpretable. To generate the combined measure we ‘modulate’ the FD by FC:

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

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

Statistical analysis using connectivity-based fixel enhancement exploits connectivity information derived from probabilistic fibre tractography. 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 voxel_mask.mif -mask voxel_mask.mif -select 20000000 tracks_20_million.tck

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

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

Statistical analysis is performed using connectivity-based fixel enhancement, with a separate analysis for FD, 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

Where 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 that for correlation analysis, a column of 1’s will not be automatically included (as per FSL randomise). Note that 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.

Note

We recommend having no more than 500,000 fixels in the analysis fixel mask (you can check this by mrinfo -size ../template/fixel_template/mask.mif, and looking at the size of the image along the 1st dimension), otherwise fixelcfestats will run out of RAM. A mask with 500,000 fixels will require a PC with 128GB of RAM for the statistical analysis step. To reduce RAM requirements, you could reduce the number of fixels by reducing the extent of upsampling, or apply a higher threshold when computing the fixel analysis mask (at the risk of removing WM regions from your analysis); further discussion in the `Commands crashing due to memory requirements`_ section.

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.