-seed_dynamic and -seed_random_per_voxel

1

I am running some analyses on patients with stroke. I have a couple of questions, but before that let me tell you the context of the analyses.

I use a lesion filling algorithm and obtain a healthy brain on a MPRAGE, which has almost identical folding in the healthy (non-lesioned) parts of the brain as my subject, but also newly imputed folding patterns on the lesioned area. I obtained the 5TT image from this “imputed” MPRAGE. Then I run tractography with the -act -seed_dynamic options. It takes about 2 days to obtain 15M streamlines, after which I run tcksift to bring it down to 2M. Note, the underlying DWI is the original lesioned one, but the 5TT driving the seeds is from the imputed MPRAGE.

Here are my questions:

  1. I am not sure sure I understand the difference between -act -seed_dynamic and -act -seed_random_per_voxel (with gmwm mask), or even -act -seed_gmwmi. Doesn’t the seed always start at the gm/wm interface once the user chooses -act? Would it be ok to seed explicitly a fixed number of streamlines per voxel in the gm/wm interface, or, as I suspect, -seed_dynamic is always the best option?
  2. Since -seed_dynamic already applied the SIFT concepts, and for this perhaps takes longer time, is it ok if the user is satisfied with an an initial 15M streamlines before applying tcksift, instead of a larger number, say 30M, obtained with other seed mechanisms which reach those numbers more quickly with more artifacts?
  3. This is a lazy question. I’ve noticed that AFD can be used for connectome analyses instead of streamline counts. What command can obtain the afd connectome from a tck file?

Thank a lot for this and other helps, and thanks again for this piece of cake software.

2

Hi Dorian,

Then I run tractography with the -act -seed_dynamic options. It takes about 2 days to obtain 15M streamlines

That’s a bit long… Are you using exotic tracking parameters at all? Are you getting a very low rate of success (# of seeds v.s. # of streamlines)? What’s your spatial resolution?

Doesn’t the seed always start at the gm/wm interface once the user chooses -act?

No. Seeding from the GM/WM interface is a capability that is enabled (in a more robust manner at least) when using ACT, but it is not enforced that all streamlines must commence from the GM/WM interface if using ACT. Seeds right in the middle of the white matter are perfectly acceptable. The determination of streamline seeds is most often independent of whether or not ACT is used; only once a seed has been determined is it validated against the 5TT image (if present). GM/WM interface seeding is the exception to this statement, since it explicitly uses the 5TT image as part of the seeding mechanism.

I am not sure sure I understand the difference between -act -seed_dynamic and -act -seed_random_per_voxel (with gmwm mask), or even -act -seed_gmwmi.

  • -seed_dynamic is very different. Details are in the appendix of the SIFT2 paper. Seeds are determined based on the FOD image, (almost) independently of the 5TT image; they are only checked against the 5TT image to make sure the seed point is valid anatomically.

  • If you were to determine a GM/WM interface mask, and use this as input to -seed_random_per_voxel, this would perform similarly to -seed_gmwmi, but not exactly the same:

    • -seed_random_per_voxel explicitly gives you a fixed number of seeds in each voxel of the mask (most seeding mechanisms in MRtrix3 give a total number of seeds, and they are distributed randomly throughout the whole mask, i.e. not necessarily the same number of seeds in each voxel).

    • -seed_gmwmi explicitly moves each seed point onto the best estimate of the precise location of the GM/WM interface, based on tissue partial volume fractions. Any “candidate” seed point within the input seed mask that cannot be converged toward this interface is discarded. Using -seed_random_per_voxel with -act, this mechanism would not be used; instead, any seed point within the grey matter (about half of them if your input mask contains GM/WM interface voxels) would be discarded, since it would be impossible for ACT to propagate the streamline.

Since -seed_dynamic already applied the SIFT concepts, and for this perhaps takes longer time, is it ok if the user is satisfied with an an initial 15M streamlines before applying tcksift, instead of a larger number, say 30M, obtained with other seed mechanisms which reach those numbers more quickly with more artifacts?

Sure. Compared to WM or GM/WM interface seeding, dynamic seeding will give you:

  • A better fit to the data for the same proportion of streamlines removed.

  • A comparable fit to the data for a reduced proportion of streamlines removed.

It’s all a balancing act though.

This is a lazy question. I’ve noticed that AFD can be used for connectome analyses instead of streamline counts. What command can obtain the afd connectome from a tck file?

Man, I really gotta finish that paper…

In simplistic terms, you can thing of a track-count-weighted connectome following SIFT as an “AFD-weighted connectome”. The terminology of that statement is not as precise as I like all things to be; but it at least captures the intent of the entire SIFT model / method(s).

3

Thanks for the detailed response.

Here are the details.

Preprocessed DWI:

************************************************
Image:               "DWI_denois_preproc_nobias.mif"
************************************************
  Dimensions:        96 x 96 x 55 x 80
  Voxel size:        2.5 x 2.5 x 2.5 x 1
  Data strides:      [ -2 3 4 1 ]
  Format:            MRtrix
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1           0           0      -110.1
                               -0           1           0      -87.53
                               -0           0           1      -52.51
  command_history:   /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrconvert "dwi_post_eddy.nii.gz" "result.mif" "-stride" "-2,3,4,1" "-fslgrad" "dwi_post_eddy.eddy_rotated_bvecs" "bvals" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
                     /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrconvert "result.mif" "/data/jag//mrtrix/DWI_denois_preproc.mif" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
                     /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrconvert "/data/jag//mrtrix/DWI_denois_preproc.mif" "/data/jag//mrtrix/dwibiascorrect-tmp-1OUN2F/in.mif" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
                     /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrcalc "in.mif" "bias.nii" "-div" "result.mif" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
                     /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrconvert "result.mif" "/data/jag//mrtrix/DWI_denois_preproc_nobias.mif" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
  comments:          FSL5.0
  dw_scheme:         0,0,0,0
  [80 entries]       0,0,0,0
                     ...
                     -0.9950762097,-0.07915032335,-0.05965369467,1100
                     0.4215337843,0.5028266271,-0.7546354429,1100
  mrtrix_version:    3.0_RC1-83-gc16a3d0e

FOD file:

************************************************
Image:               "WMfod.mif"
************************************************
  Dimensions:        96 x 96 x 55 x 45
  Voxel size:        2.5 x 2.5 x 2.5 x 1
  Data strides:      [ -2 3 4 1 ]
  Format:            MRtrix
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1           0           0      -110.1
                               -0           1           0      -87.53
                               -0           0           1      -52.51
  command_history:   /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrconvert "dwi_post_eddy.nii.gz" "result.mif" "-stride" "-2,3,4,1" "-fslgrad" "dwi_postdy_rotated_bvecs" "bvals" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
  [6 entries]        /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrconvert "result.mif" "/data/jag/mrtrix/DW_preproc.mif" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
                     ...
                     /share/apps/MRtrix3/2017-05-30/mrtrix3/bin/mrconvert "result.mif" "/data/jag/mrtrix/DW_preproc_nobias.mif" "-nthreads" "4" "-quiet"  (version=3.0_RC1-83-gc16a3d0e)
                     dwi2fod "msmt_csd" "DWI_denois_preproc_nobias.mif" "resp_wm.txt" "WMfod.mif" "resp_gm.txt" "gmfod.mif" "resp_csf.txt" "csffod.mif" "dwimask.nii.gz" "-quiet" "-nthreads" "4"  (version=3.0_RC1-83-gc16a3d0e)
  comments:          FSL5.0
  mrtrix_version:    3.0_RC1-83-gc16a3d0e
  prior_dw_scheme:   0,0,0,0
  [80 entries]       0,0,0,0
                     ...
                     -0.9950762097,-0.07915032335,-0.05965369467,1100
                     0.4215337843,0.5028266271,-0.7546354429,1100

Track file before sift.

***********************************
  Tracks file: "FullACTtracts.tck"
    act:                  5tt_on_DWI_1mm.nii.gz
    backtrack:            1
    count:                9523144
    crop_at_gmwmi:        1
    downsample_factor:    3
    fod_power:            0.25
    init_threshold:       0.100000001
    lmax:                 8
    max_angle:            60
    max_dist:             300
    max_num_seeds:        15000000000
    max_num_tracks:       15000000
    max_seed_attempts:    1000
    max_trials:           1000
    method:               iFOD2
    min_dist:             10
    mrtrix_version:       3.0_RC1-83-gc16a3d0e
    output_step_size:     0.5
    rk4:                  0
    samples_per_step:     4
    seed_dynamic:         WMfod.mif
    sh_precomputed:       1
    source:               WMfod.mif
    step_size:            0.5
    stop_on_all_include:  0
    threshold:            0.100000001
    timestamp:            1496684523.2961580753
    total_count:          17893779
    unidirectional:       0
    ROI:                  seed WMfod.mif

As you can see, there is about twice seeded as selected streamlines. Note that this is a brain with stroke (~20% of the brain missing). If you see anything out of place, please let me know so I can fix things before starting the analyses on all subjects.

This is an example of final tractography after sift, the dark part is where the lesion is located from the manual segmentation.

image.jpg715×841 90.2 KB

4

As you can see, there is about twice seeded as selected streamlines.

That’s about normal for whole-brain tracking with ACT with back-tracking enabled.

The only bits that caught my eye are these:

max_angle:            60
step_size:            0.5

This will slow down the tracking a bit; being a second-order method, iFOD2 doesn’t require as small a step size as some other algorithms. The default is half a voxel, or 1.25mm in your case. A large angle (particularly in combination with a small step size) may also make the tracking more prone to looping around in the white matter, hence extending the typical streamline length and correspondingly the execution time for a fixed number of streamlines. These doesn’t fully explain why your tracking is taking so long though (unless you’re running it single-threaded, or on a very old machine?).

I don’t see anything glaringly wrong with your result.

5

Thank you @rsmith

How about running iFOD2 with 1mm step size and 45 degree angle? Would this be more balanced or did you have other values in mind?

Dorian

6

The default is step size = 0.5 x voxel_size and 45 degrees, though this was tested primarily on lower-resolution data when we came up with the defaults. 1mm and 45 degrees would probably work OK. But there’s no right answer.

One thing I forgot about with regards to execution time is the iFOD2 calibrator mechanism. I’m guessing that with an angle of 60 degrees, a greater number of calibration paths will need to be sampled at each step to set up the rejection sampling, which will correspondingly increase the computation time (you can see this information by using the -info option with ``tckgen`).