FDT Tractography Practical

In this practical you are going to learn how to use the outputs of BEDPOSTX, FSL's tool for estimating diffusion along multiple fibre populations. First, we explore the outputs of bedpostx, and run some tractography analyses on them (using the bedpostx model outputs to estimate structural connectivity from a seed to the rest of the brain). We will explore the outputs of bedpostx where a single fibre population has been modelled, and then look at how our tractography results change when we estimate three fibre populations, and use different fibre estimation models. We will also look at how to constrain the tractography search using anatomical masks.

Contents:

Bedpostx output
Exploring the results of bedpostx, estimating a single fibre population per voxel.
Tractography
Generating connectivity distributions from single voxels
Connectivity-based segmentation
Using connectivity to classify voxels, and running probabilistic tractography from seed masks
Crossing-Fibre Tractography
Generating connectivity distributions using multi-fibre directions
Anatomical priors in tractography
Using exclusion/waypoint/termination masks to inform tractography
Bedpostx model specification (optional)
Choosing the right model to improve data fitting
Using surfaces to constrain tractography (optional)
Using surfaces to avoid spurious connections

Bedpostx output

cd ~/fsl_course_data/fdt2/subj1_model2_1fibre.bedpostX

This is what a bedpostx directory looks like if you use only a single fibre orientation per voxel (i.e. bedpostx -n 1). You will see a three-fibres bedpostx directory later.

As it takes several hours to run (unless you run it on a GPU), bedpostx has already been run for you on the data in subj1. It is possible to run bedpostx from either the FDT GUI or via the command line. Before you run it you need to have a directory setup that contains data files with standard names, as you can see in the FDT User Guide on the FSL Wiki. In addition to the necessary input files, the current directory also contains all the output files essential for tractography (as output by bedpostx).

For now, ignore all the files beginning with MASK. They will be used later.

Open mean_f1samples in FSLeyes. This is the mean of the distribution on the anisotropy parameter in the partial volume model - it should look very similar to the FA map you saw in the dtifit practical.

Add dyads1, then open the overlay display panel () and change the Overlay data type to 3-direction vector image (Line). This is the mean of the distribution of fibre orientations - it should look very similar to the dti_V1 from the dtifit practical.

Now add dyads1_dispersion. You need to adjust the display range - set the minimum to 0.005 and maximum to 0.03. Also change the colour map to Red-Yellow. The dyads1_dispersion image encodes for each voxel the uncertainty (i.e. spread of the distribution) around the mean fibre orientation dyads1. As you may expect, highly anisotropic voxels in white matter have very low dispersion, while grey matter and CSF voxels have orders of magnitude higher dispersion. Notice that some voxels in white matter also have high dispersion. An example of such voxels can be found in the cerebral white matter (centrum semiovale). Set Y=60 and look at the coronal view.

Why do some voxels in white matter also have high dispersion?
Not really. That could be the case if our b-values were really high and the SNR of our data really low. In that case, though, we would observe a general increase in fibre dispersion, not a white matter-specific one.
Correct! High dispersion values in white matter indicate complex fibre arrangement (i.e. crossing fibres).
Wrong! It is actually the other way round.

In addition to the mean fibre orientation distribution, bedpostx returns samples from this distribution that are utilised in probabilistic tractography. These are contained in merged_th1samples and merged_ph1samples. Fibre orientations are described in spherical coordinates* for each orientation sample. Viewing the merged_ files with FSLeyes is not particularly informative. In general, the dyads are used for visualisation.

*Spherical polar coordinates


Tractography

Simple tracking - connectivity from a single seed point.

Voxel in Internal Capsule

In FSLeyes, open the mean_f1samples map again. Find the co-ordinates of a voxel in the internal capsule (the white matter tract shown in the picture; running between the caudate/thalamus and the putamen/pallidum).

Bring up Fdt - it should default to PROBTRACKX with single voxel as the selected seed mode.

Select subj1_model2_1fibre.bedpostX as the Bedpostx directory and enter your selected voxel as the seed.

Choose an output name, and press Go.

Back in FSLeyes, overlay the connectivity distribution on top of the mean_f1samples map - note the pathway's name ends with the coordinates of the seed voxel (basename_X_Y_Z.nii.gz). Adjust the min and max display thresholds to look at the relative probabilities of paths. Change the colour map to Red-Yellow. Note the narrowness of the spatial distribution of the tractography results - this is because there is very low uncertainty in the Internal Capsule.

digraph G { rankdir=LR node [shape=box]; 1 [label="basename for samples files"]; 2 [label="binary BET mask file in diffusion space"]; 3 [label="seed voxel coordinate (text/ascii file)"]; 4 [label="reference vol to define seed space in simple mode"]; 5 [label="simple mode for single seed voxel"]; 6 [label="PROBTRACKX"]; 7 [label="output file"]; 8 [label="directory to put the final volumes in"]; 1 -> 6; 2 -> 6; 3 -> 6; 4 -> 6; 5 -> 6; 6 -> 7; 6 -> 8; }

Tracking from a voxel in a non-diffusion space.

Seed points can also be specified in a space different from the native diffusion one (e.g., high resolution T1-weighted anatomical space).

Open subj1_model2_1fibre.bedpostX/struct_brain in FSLeyes and find an interesting seed location.

In the Fdt GUI, remaining in single voxel seed mode, select seed space is not diffusion. This should bring up two new fields, allowing you to specify the seed space:

Run this tractography and load the output (basename_X_Y_Z) into FSLeyes on top of the structural image.

digraph G { rankdir=LR node [shape=box]; 1 [label="basename for samples files"]; 2 [label="binary BET mask file in diffusion space"]; 3 [label="seed voxel coordinate (text/ascii file)"]; 4 [label="reference vol to define seed space in simple mode"]; 5 [label="simple mode for single seed voxel"]; 6 [label="transformation matrix from seed space to DTI space (FLIRT matrix)", color=red, fontcolor=red]; 7 [label="PROBTRACKX"]; 8 [label="output file"]; 9 [label="directory to put the final volumes in"]; 1 -> 7; 2 -> 7; 3 -> 7; 4 -> 7; 5 -> 7; 6 -> 7; 7 -> 8; 7 -> 9; }

Tracking from voxels in MNI152 standard space

You can also specify the seed points in standard space mm (e.g. FMRI activation blobs output from FEAT). Either linear or nonlinear transformations to standard space can be utilised. In the former case (linear), use standard2diff.mat as a transform and standard as a reference image. In the latter (nonlinear), use standard2diff_warp and diff2standard_warp as transforms and standard as a reference image. Try the nonlinear option (and don't forget to tick the nonlinear checkbox).

Don't try too many tracts now because you will see a big difference when we run tractography with three fibres per voxel.

digraph G { rankdir=LR node [shape=box]; 1 [label="basename for samples files"]; 2 [label="binary BET mask file in diffusion space"]; 3 [label="seed voxel coordinate (text/ascii file)"]; 4 [label="reference vol to define seed space in simple mode"]; 5 [label="simple mode for single seed voxel"]; 6 [label="transformation matrix from seed space to DTI space (FLIRT matrix)", color=red, fontcolor=red]; 7 [label="transformation warp from DTI space to seed space (FNIRT warpfield)", color=red, fontcolor=red]; 8 [label="PROBTRACKX"]; 9 [label="output file"]; 10 [label="directory to put the final volumes in"]; 1 -> 8; 2 -> 8; 3 -> 8; 4 -> 8; 5 -> 8; 6 -> 8; 7 -> 8; 8 -> 9; 8 -> 10; }

Connectivity-based segmentation

cd ~/fsl_course_data/fdt2/subj1_model2_1fibre.bedpostX

A connectivity-based seed classification analysis has been run for you and is saved in a sub-directory called THAL2CTX_right within your subj1.bedpostX directory. This analysis was run via the FDT GUI, as described in the wiki.

This analysis was carried out in standard space. cd back into subj1_model2_1fibre.bedpostX. Open the standard brain in FSLeyes and overlay the seed mask, MASK_average_thal_right - you should see that the seed is in the thalamus.

In the subj1_model2_1fibre.bedpostX directory, there is a collection of files called MASK_average_*. If you overlay these files with a standard brain, you will see the target masks used in this analysis. Overlay some of these masks in your open FSLeyes window - they are different cortical zones (you will need to change the colour maps to be able to view these as separate masks).

The objective of this type of analysis is to ask the following question: For each location in the seed mask, what is the relative probability of connection to each of my target masks? The output is therefore a single image for each target mask in which only voxels within the seed mask contain data. At each voxel within the seed mask the voxel value is the number of samples which reached this target mask from this seed voxel. These outputs can be found in THAL2CTX_right in the subj1_model2_1fibre.bedpostX directory.

Now overlay in FSLeyes, in turn, each of the seeds_to... files in THAL2CTX_right. Compare the spatial distributions and probabilities of connection from different thalamic voxels to different cortical target zones.

Run find_the_biggest on the outputs of seeds to targets to generate a hard segmentation of the thalamus and overlay this segmentation onto the standard brain. Something like:

find_the_biggest THAL2CTX_right/seeds_to_* biggest_segmentation

The different numbers in the output image (biggest_segmentation) each relate to a different target mask. You may want to change the colour map to something like Random.


Crossing fibre tractography

cd ~/fsl_course_data/fdt2/subj1_model2_3fibres.bedpostX

This is a three-fibre bedpostx analysis output directory.

Just like before, mean_f1samples tells you the estimated proportion of the diffusion signal that can be accounted for by the first fibre orientation. Now, mean_f2samples and mean_f3samples tell you the same for the second and third orientations. dyads1, dyads2 and dyads3 are the mean estimates for the orientation vectors. Load these six images into FSLeyes.

First turn off dyads1, dyads2 and dyads3 and change mean_f2samples and mean_f3samples to use different colour maps. Threshold both of them at 0.05 (that is, set the minimum value in the display range to 0.05 - we will not actually threshold/change the data, just the display).

You should see that in much of the white matter, there are two fibre populations, but that in non-white matter tissue, the second and third fibres have very small volume fraction (e.g., <10e-6). This is due to the ARD (automatic relevance determination, part of the Bayesian analysis of the crossing fibres). Notice that in some white matter regions where a single coherent orientation is expected (e.g., corpus callosum), the second and third volume fractions are also very small.

Now turn back on dyads1, dyads2 and dyads3 and display them as line vectors ( -> Overlay data type -> 3-direction vector image (Line)). Change the X, Y, and Z Colour to red for dyads1, to green for dyads2 and to blue for dyads3. You will note that there are green and blue lines everywhere (even where there is no evidence for a second fibre in mean_f2samples), but you should notice the green and blue lines are well ordered in the regions where 3 fibres are supported (where mean_f2samples and mean_f3samples survive the threshold), but random elsewhere.

In the bedpostx folder, there are two images called dyads2_thr0.05 and dyads3_thr0.05. These contain average orientations only in those voxels where 2 or 3 fibres are supported (i.e., their anisotropic fibre fractions are > 0.05). You can change the threshold by using the maskdyads command in your shell.

Reopen FSLeyes with this command:

fsleyes mean_f1samples.nii.gz \
  dyads1.nii.gz         -ot linevector -xc 1 0 0 -yc 1 0 0 -zc 1 0 0 -lw 2 \
  dyads2_thr0.05.nii.gz -ot linevector -xc 0 1 0 -yc 0 1 0 -zc 0 1 0 -lw 2 \
  dyads3_thr0.05.nii.gz -ot linevector -xc 0 0 1 -yc 0 0 1 -zc 0 0 1 -lw 2

You should now be able to see some beautiful crossing fibre architecture!

Now try to visualise these crossing fibres in standard space. For that you will have to transform the dyads, applying the registration transformation in xfms/diff2standard_warp, using the program vecreg. Type:

vecreg -i dyads1 -o dyads1_to_standard \
  -r standard -w xfms/diff2standard_warp

Do the same for dyads2_thr0.05 and dyads3_thr0.05. Open FSLeyes with a standard brain and the three registered dyads.

Try some of the same tractography you have already done in the previous section but, this time, specify subj1_model2_3fibres.bedpostX as your subject directory. You should find that you can see many pathways that were previously invisible. For example, see if you can find the medial portion of the SLF in the 3-fibre data (see example in lecture).

Try tracking from cortical grey matter - this is likely to produce variable results due to high uncertainty on fibre direction in the cortex (try moving into the white matter adjacent to the grey matter area of interest).


Anatomical priors in tractography

cd ~/fsl_course_data/fdt2/subj1_model2_3fibres.bedpostX

Here we will make use of target masks to inform tractography. Using a seed mask in the lateral geniculate nucleus of the thalamus (LGN), we will construct the pathway connecting LGN to the ipsilateral primary visual cortex (V1) that passes via Meyer's loop.

Voxel in Internal Capsule

Run tractography from the LGN mask (file called SEED_LGN_LEFT.nii.gz). You will need to change the seed mode from single voxel to single mask. Note that this mask is in standard space (i.e. use the transformation as in the example above). You may want to reduce the number of samples to 500 to speed up tractography (under the Options tab).

Use a waypoint mask in V1 (TARGET_V1_LEFT.nii.gz) to isolate only those tracts that reach V1. Use this V1 mask also as a termination mask to avoid tracts that reach V1 and then continue to other parts of the brain.

This pathway is called the optic radiation. Try to visualise it, on top of the standard image, in FSLeyes.

There are still some artefactual tracts in the results. How can you get rid of those?
Wrong!
Wrong! A stop mask will still retain the spurious fibres.
Correct! Draw an exclusion mask in FSLeyes and include it in the tractography. (or use the file EXCLUSION.nii.gz)

Bedpostx model specification (optional)

cd ~/fsl_course_data/fdt2

In this section we will see how the choice of the bedpostx model can affect the fitting results.

We will compare model 1 (ball & sticks, single diffusion coefficient) to model 2 (ball and sticks, continuous gamma distribution of diffusion coefficients).

Model 2 should be better suited to deal with multi-shell data, as the apparent diffusion coefficient can change with the b-value. Using model 2 tends to reduce overfitting (i.e., spurious fibre orientations).

Open FSLeyes by typing the following command in your terminal:

fsleyes \
  subj1_model2_3fibres.bedpostX/str2diff.nii.gz -or 0 250 \
  subj1_model1_3fibres.bedpostX/mean_f2samples.nii.gz     \
    -cm blue-lightblue -or 0.05 0.2                       \
  subj1_model2_3fibres.bedpostX/mean_f2samples.nii.gz     \
    -cm red-yellow -or 0.05 0.2                           \
  subj1_model1_3fibres.bedpostX/dyads1.nii.gz             \
    -ot linevector -xc 0 0 0 -yc 0 0 0 -zc 0 0 0 -lw 2    \
  subj1_model1_3fibres.bedpostX/dyads2_thr0.05.nii.gz     \
    -ot linevector -xc 1 1 1 -yc 1 1 1 -zc 1 1 1 -lw 2 &

This loads the fitted orientations using model 1 (single diffusion coefficient), as well as the volume fractions of the secondary fibre orientation for both models 1 and 2.

First, deselect the dyads files, and examine the volume fractions (mean_f2samples). Model 1 should be displayed in Blue and model 2 in Red-Yellow.

The first thing that is immediately visible is that many more voxels survive a 0.05 thresholding on the volume fractions of the secondary fibre when using model 1. This is particularly true at interfaces between grey/white/CSF tissues, where partial volume is likely to occur, and it indicates that the model is overfitting the data.

Now examine the dyads in the voxels where a secondary fibre is estimated when using model 1 but not model 2 (i.e., the blue voxels). You can observe that in most of them, especially in those closer to tissue boundaries (grey/white/CSF interface) the second dyad is perpendicular to the first one, but does not show spatial continuity when compared to neighbouring dyads.

Partial volume effects can cause the signal to decay non mono-exponentially when increasing the b-value. Since model 1 has a single diffusion coefficient, it can only capture this feature of the data by adding a perpendicular second dyad, which fits the data better but does not correspond to a genuine (anatomical) fibre orientation.


Using surfaces to constrain tractography (optional)

cd ~/fsl_course_data/fdt2/surfaces

All previous examples use volumetric seeds and constraint masks. Probtrackx2 also supports surface files (e.g., generated using FreeSurfer or Caret)

In this example use of surfaces, we consider the case where FreeSurfer has been used to generate a cortical surface, and show how such surface can be used in probtrackx2.

In the current folder you will find the high-resolution T1-weighted volume produced by the FreeSurfer recon-all command (orig.nii.gz), the necessary transformation matrices to go from FreeSurfer anatomical space to diffusion space and the white/grey matter boundary surfaces for the left (lh.cortex.gii) and right (rh.cortex.gii) hemispheres. For details on how to generate these files, see the relevant documentation.

Here, we will see an example where, in order to avoid spurious connections that connect adjacent gyri, we impose the convoluted WM/GM boundary surface as a termination mask. Tracking will be performed again in anatomical space as defined by the orig.nii.gz volume. We will first run an example without any surface constraints from a seed close to the cortex.

The command:

echo "97 96 125" > svox.txt

creates a text file called svox.txt that contains the string 97 96 125. This could also be created using your preferred "plain" text editor.

In your terminal type:

echo "97 96 125" > svox.txt
probtrackx2 \
  --samples=../subj1_model2_3fibres.bedpostX/merged \
  --mask=../subj1_model2_3fibres.bedpostX/nodif_brain_mask.nii.gz \
  -x svox.txt \
  --xfm=freesurfer2fa.mat --seedref=orig.nii.gz \
  --loopcheck --simple --forcedir \
  --opd -V 1 --dir=Tracts --out=no_surf

The next commands will include as termination masks the wm/gm boundaries of both hemispheres:

echo lh.cortex.gii > stop.txt
echo rh.cortex.gii >> stop.txt
probtrackx2 \
  --samples=../subj1_model2_3fibres.bedpostX/merged \
  --mask=../subj1_model2_3fibres.bedpostX/nodif_brain_mask.nii.gz \
  -x svox.txt \
  --xfm=freesurfer2fa.mat --seedref=orig.nii.gz \
  --loopcheck --simple --forcedir --opd -V 1 \
  --dir=Tracts --out=surf --meshspace=freesurfer --stop=stop.txt

After the results have been computed, open FSLeyes with the following command:

fsleyes -ds world \
  orig.nii.gz rh.cort.nii.gz -cm green \
  lh.cort.nii.gz -cm copper \
  Tracts/no_surf_97_96_125 -cm blue-lightblue -or 10 1000 \
  Tracts/surf_97_96_125 -cm red-yellow -or 10 1000

This loads the tractography results overlaid on top of the high-resolution anatomical volume together with the surface boundaries. The results obtained when using the surfaces as termination masks should be displayed in red-yellow and those obtained without surface contraints should be in blue-lightblue.

Go to slice Y=91. Observe how, when surfaces are not used as termination masks, the tracts in blue wrongly connect two adjacent gyri by passing through a sulcus.


The End.