Recommended tractography options for low b-value?

[kjamison]

I’ve been trying out various combinations of tckgen and tcksift options on our dataset and so far I haven’t found a satisfactory combination. I was wondering I might be missing something or if there was some kind of general guidance available for this situation.

The dataset: 1.8x1.8x2.5mm, 64 directions in a single b=800 shell (+ 1 b0), and includes a standard EPI fieldmap as well as a 1mm T1 anatomical image. We performed basic topup and eddy correction, though I have noticed some residual disagreement between diffusion and T1 around the frontal pole for many subjects. We have run freesurfer on the T1, and performed manual editing of these freesurfer reconstructions. Many of our subjects have white matter abnormalities, so we would prefer to use these manually edited segmentations where possible. T1 was registered to diffusion after freesurfer.

We’ve been reading about the reduced false positive rates using CSD+deterministic tractography+filtering, and based on this post on single shell data, our current rough-draft pipeline is:

dwi2response dhollander DWI.mif RF_wm_dhollander.txt RF_gm_dhollander.txt RF_csf_dhollander.txt -voxels RF_voxels_dhollander.mif
dwi2fod msmt_csd DWI.mif RF_wm_dhollander.txt RF_wm_dhollander.mif RF_csf_dhollander.txt RF_csf_dhollander.mif
tckgen RF_wm_dhollander.mif CSD_sdstream_1M.tck -algorithm SD_STREAM -seed_dynamic RF_wm_dhollander.mif -select 1M
tcksift CSD_sdstream_1M.tck RF_wm_dhollander.mif CSD_sdstream_1M_sift500K.tck -term_number 500k

I realize 1M->500K is not a very rigorous filter, and we could use larger values, but I’m hesitant to commit significant compute resources (we have >100 subjects) until I see outputs that suggest I’m on the right path. The current outputs are dominated by short fibers within gray matter . I would like to restrict these in a meaningful way without being overly prescriptive. I have played with running 5ttgen freesurfer (to make use of our manual edits) and using the -act options in tckgen and tcksift. Just running tckgen with -act has taken tractography from 40 minutes to what looks like 6 or more hours per subject, not including tcksift.

As I said, if this is the “correct” approach, I’m willing to go for it, but I was hoping for some expert input before doing so, especially since there are many additional algorithm or parameter choices that I haven’t even considered that might be having more of an impact in this case.

[mblesac]

Hi,

I can think in a couple of things, that maybe are useful for your pipeline:

I have noticed some residual disagreement between diffusion and T1 around the frontal pole for many subjects.

Have you tried to align the b0 and the T1w using non-linear registration, this should improve your alignment.

One parameter you could add to tckgen is -cutoff 0.1 for me it worked in the past removing a lot of noisy tracts, specially when using the 2 tissues response function approach.

Another approach that maybe you are interested to try is the recently published SS3T-CSD algorithm.

What I don’t really understand, is why by adding -act to your command, the time increases this much. Have you check the alignment of your 5TT file with your B0? is it good? The only think It comes to my mind is that act is rejecting a lot of tracts due to a wrong termination (such as inside WM).

We’ve been reading about the reduced false positive rates using CSD+deterministic tractography+filtering

I usually prefer the other approach, to create a very dense probabilistic connectome (10M for example) and then apply SIFT2. If my understanding is correct the FP should receive a weight very close to 0 and very small compared to real path. I hope this helps.

Best regards,

Manuel

[kjamison]

Thanks for the insight!

@mblesac

What I don’t really understand, is why by adding -act to your command, the time increases this much. Have you check the alignment of your 5TT file with your B0? is it good? The only think It comes to my mind is that act is rejecting a lot of tracts due to a wrong termination (such as inside WM).

Yeah I am confused about this as well. It actually ended up taking 15 hours, compared to 40 minutes for the exact same tckgen but without “-act”. Registration seems to be fine. Here’s the GM segment overlaid on diffusion FA in mrview:

The output also seems a little TOO anatomically constrained, so I’m a little hesitant to continue on that route:

@mblesac

I usually prefer the other approach, to create a very dense probabilistic connectome (10M for example) and then apply SIFT2. If my understanding is correct the FP should receive a weight very close to 0 and very small compared to real path.

I may explore this route. I did try tcksift2 but the output ended up being all 1s. I just used the default parameters

tcksift2 CSD_sdstream_5M.tck RF_wm_dhollander.mif CSD_sdstream_5M_sift2.txt

Is there any specific option that is required?

[mblesac]

Hi,

My guessing would be that this is due to the combination of the deterministics algorithm and the dinamyc seeding mechanism. The seeding mechanism keeps starting tracts in the empty spaces that are rejected because some of the constrains, act adds more constrains making it even more difficult. Could you try the same command changing the seeding strategy or the tracking algorithm? and adding act? to check the time. I don’t know if I’m correct, probably @rsmith could give you better advice.

Also from SIFT2 , from my poimt of view the output is not very surprissing due the small number of streamlines and the lack of the mask option (or act). Could you try with more tracts and the mask option? Again @rsmith would have a better opinion on this.

I hope this helps,

Best regards,

Manuel

[kjamison]

Regarding SIFT2 what do you mean small number of streamlines? I had a few mistakes in what I posted, which I’ve updated now. I’m starting with 5M tracks, which seems like enough to do SOME sorting. Do I need to set a different regularization parameter? The doc says default -reg_tikhonov is 0, and default -reg_tv is 0.1, so I assume it’s using TV, but I have no intuition for how to set that (edit: or if some of the other tcksift2 arguments are necessary).

[mblesac]

Hi,

Sorry, I don’t know the details about that. I would add the -mask flag and see the results. 5M should be enough to see different values.

Best regards,

Manuel

[ThijsDhollander]

@mblesac

Another approach that maybe you are interested to try is the recently published SS3T-CSD algorithm .

Thanks for the shout-out Manuel.

@kjamison it’s definitely worthwhile trying, but it’s always going to be slightly challenging to perform any FOD-based tractography with the very low b-value, and the single b=0 image. That said, the 64 gradient directions may make up for that a bit. I’ve seen this before, but it’s an odd acquisition: at such a low b-value, the angular contrast is so low (and thus the signal so “smooth”) that 64 directions is a bit overkill.

@kjamison

Yeah, the sdstream result there definitely doesn’t look good. One of the main things that stands out to me is the almost complete lack of streamlines in the splenium (!). Combined with your massive running time with ACT there, it does look like a lot of streamlines are crashing and burning in the middle of the WM.

I think you’d best stick to probabilistic tractography by default for your data, and look into tuning (all) the other choices and parameters. Looking at your first result:

@kjamison

.

I agree with Manuel there, you’d need to look into the -cutoff parameter for sure.

@mblesac

One parameter you could add to tckgen is -cutoff 0.1

That might be slightly drastic though, I’m guessing. This might date back to before a bug with SH amplitudes was fixed. I’d suggest starting maybe with a cutoff value of 0.07 or something, and adjusting up or down based on what you see. Also, make sure to “first” get these parameters right, by running without ACT while testing maybe. This has the benefit of running faster (so you get faster feedback from your results to adjust parameters), but also that it’s a bit more direct or simple. It’s harder to understand what goes “wrong” when all mechanisms combine; it’s sometimes easier to diagnose stuff by focusing on parts of the challenge instead, and gradually introduce more complexity when all other boxes are ticked.

Even though running SS3T-CSD on such low b-value data is extra challenging sometimes, it actually brings greater benefits for these low b-value data in particular for tractography. The problem you’ve got in finding a reasonable cutoff is due to the contrast / intensity of non-WM tissues at that b-value. If you can even filter out some of that, it helps massively with separating WM from the rest via a certain cutoff value.

Well, see how you go. Happy to help out or provide input if you’ve got any further questions.

Cheers,
Thijs

[kjamison]

Thanks!

Here is the cross-section with 1M fibers at cutoff 0.10:

Here is 0.15;

0.15 seems drastic, but 0.10 seems reasonable to me.

I’ll check out SS3T-CSD. I ran it but haven’t tried tractography yet. One note: I ran it with -nthreads 0 to avoid eating up cluster resources (and get an apples-to-apples timing comparison), and it does not seem to have respected that constraint. CPU usage during various subprocesses was definitely using up all available cores. Is there some other way to restrict CPU usage?

Additional note on SIFT2: I found that running tcksift2 on the 5M sdstream tracks gave me a vector of all 1s. This was true when I provided an explicit -proc_mask as well. It seems to give a range of values when I additionally provide -fd_thresh 0.1, which I did on a whim thinking it was equivalent to the FA cutoff I used for tckgen, but now that I think about it, I’m unclear on what exactly it did. Is this a reasonable parameter choice? Is there some other approach to get reasonable sift2 values in this scenario?

[ThijsDhollander]

@kjamison

0.15 seems drastic, but 0.10 seems reasonable to me.

Wow, definitely higher values than I expected above indeed. This is no doubt due to the low b-value in addition to absence of a GM compartment. Your 0.15 results doesn’t even look all “that” drastic to me. But it’s a tough choice either way with those amplitudes you’re facing and trying to “separate” between true and false ones.

@kjamison

I’ll check out SS3T-CSD. I ran it but haven’t tried tractography yet.

Nice; as mentioned, it’ll still be far from perfect due to the data quality, but I’m keen to learn how it goes.

@kjamison

One note: I ran it with -nthreads 0 to avoid eating up cluster resources (and get an apples-to-apples timing comparison), and it does not seem to have respected that constraint. CPU usage during various subprocesses was definitely using up all available cores.

Hmm, good to know! I don’t have immediate time to look into it, but it’s well noted. I’m not entirely surprised though; but I’ll need to look into it at some point.

@kjamison

Is there some other way to restrict CPU usage?

Not 100% sure it’ll fix things, but see if you can set the NumberOfThreads configuration option in a configuration file: List of MRtrix3 configuration file options — MRtrix 3.0 documentation . You can store it in a system wide location or create a user specific one, as detailed here: Configuration file — MRtrix 3.0 documentation .

[ThijsDhollander]

@ThijsDhollander

kjamison:

One note: I ran it with -nthreads 0 to avoid eating up cluster resources (and get an apples-to-apples timing comparison), and it does not seem to have respected that constraint. CPU usage during various subprocesses was definitely using up all available cores.

Hmm, good to know! I don’t have immediate time to look into it, but it’s well noted. I’m not entirely surprised though; but I’ll need to look into it at some point.

@kjamison: sorry, I was wrong there, this should work perfectly out of the box. I had misread your actual option there: the problem is with -nthreads 0. If you set that to zero, it effectively aims to use as many threads as cores (so it’ll use most of all cores at most times). You should set it to any non-zero number if you want to limit eating up too much resources at once. Timing should be roughly 3x MSMT-CSD, but might be slightly slower depending on how fast disk access is on your system at hand; especially for larger files (or e.g. larger spatial resolutions / smaller voxels).

[kjamison]

OK I’ll try that. For the record, the documentation says “(set to 0 to disable multi-threading).” where “disable multi-threading” sounds like single core to me, so that’s a bit confusing.

[ThijsDhollander]

@kjamison

For the record, the documentation says “(set to 0 to disable multi-threading).” where “disable multi-threading” sounds like single core to me, so that’s a bit confusing.

You’re absolutely right! To be clear: it’s the documentation that is wrong or at least confusing there. This is an MRtrix3-wide thing in fact; for both scripts as well as binary commands. I’ve created an issue to get this resolved eventually.

[kjamison]

I’ve been using “-nthreads 0” for all of the other mrtrix tools and it has restricted to a single core. So the documentation seems to reflect the behavior for most of the tools.

[ThijsDhollander]

@kjamison

most of the tools

And again, you’re right. Ok, did another double-check, I was wrong above when I said

@ThijsDhollander

for both scripts as well as binary commands

It seems to be only the case for scripts. I’ve confirmed the difference by testing a few binaries (e.g. fod2fixel and dwi2fod, etc…) that indeed run single-threaded with -nthreads 0 and a few scripts (e.g. dwi2response tournier that spends some time at some point using fod2fixel) that indeed run multi-threaded with -nthreads 0.

Finally, do note that there’s also a good number of commands (even binary ones) that don’t ever use multi-threading; regardless of specifying -nthreads ... or not.

So in general, if you want to be in full control without ambiguity, use -nthreads with a non-zero number. If you don’t specify -nthreads at all, it’ll use as many threads as cores on your system; i.e. “maximum” (more or less) multi-threading; but only for commands that have multi-threading functionality at some point in their operation.

[jdtournier]

@kjamison

I’ll check out SS3T-CSD. I ran it but haven’t tried tractography yet. One note: I ran it with -nthreads 0 to avoid eating up cluster resources (and get an apples-to-apples timing comparison), and it does not seem to have respected that constraint. CPU usage during various subprocesses was definitely using up all available cores. Is there some other way to restrict CPU usage?

Just to clarify things here:

• -nthreads 0 should explicitly disable multi-threading, in C++ and Python commands.

• I’ve had a quick look through the code, this should all be handled properly. The Python scripting backend will invoke all commands with -nthreads 0 if provided, even if this is not reported on the terminal output (you can see the full command actually invoked if you run with -debug). I’ve checked this with dwi2response, it behaves as expected.

• Just for the record, running with -threads 1 is not the same as -nthreads 0. In the first case, this sets the number of worker threads to 1, but there may very well be other threads launched to interact with that (e.g. in a pipeline infrastructure, we might still have source and sink threads launched to respectively provide data to and consume the output of the worker thread that does the heavy lifting). With -nthreads 0, everything should run sequentially on the main thread.

So as far as I can tell from the code and my limited testing, everything is OK. @kjamison, can you confirm whether this happens only with SS3T-CSD, or do we need to investigate anything within MRtrix3 itself?

[jdtournier]

On further testing, I note that this option is indeed ignored on our development branch – which SS3T-CSD is based on. We’ll fix that prior to release, thanks for pointing it out.

[ThijsDhollander]

Good to hear this conversation was helpful. Thanks @kjamison for spotting and flagging this.

[kjamison]

Happy to help! While I have your attention, can I maybe ask about this issue again?

@kjamison

I found that running tcksift2 on the 5M sdstream tracks gave me a vector of all 1s. This was true when I provided an explicit -proc_mask as well. It seems to give a range of values when I additionally provide -fd_thresh 0.1 , which I did on a whim thinking it was equivalent to the FA cutoff I used for tckgen , but now that I think about it, I’m unclear on what exactly it did. Is this a correct interpretation? Is there some other approach to get actual sift2 values in this scenario?

[ThijsDhollander]

No, that option isn’t an FA threshold. The one in tckgen for FOD-based algorithms actually also isn’t an FA threshold. The one in tckgen is an FOD amplitude threshold, so lower amplitudes can be disregarded. The one in tcksift2 appears to be an FOD lobe integral threshold, as the option seems to describe at least. Not sure why you’d get all 1s without such a threshold though.

[rsmith]

Just running tckgen with -act has taken tractography from 40 minutes to what looks like 6 or more hours per subject, not including tcksift.

This is most likely a shift in the fraction of generated streamlines that are deemed satisfactory by all criteria, such that a greater number of streamlines need to be generated in order to produce the desired number. Running tckgen with -info might provide some insight.

My guessing would be that this is due to the combination of the deterministics algorithm and the dinamyc seeding mechanism.

That would potentially have an effect; not sure to what extent I experimented with that combination. There should be at least some protection against seeding probabilities becoming degenerate; but I wouldn’t know for sure unless I poked around at the innards when running against such data. If you wanted, you could try uncommenting this line and recompiling, and see if the extra outputs give any insight. But generally you might be better off not combining dynamic seeding with a deterministic algorithm.

If my understanding is correct the FP should receive a weight very close to 0 and very small compared to real path.

Not really. SIFT2 doesn’t have a tailored mechanism by which to classify streamlines as false positives. All it can do is modulate the relative densities of streamlines pathways in order to equalise them with the FODs.

I did try tcksift2 but the output ended up being all 1s.

Never seen that before…

It seems to give a range of values when I additionally provide -fd_thresh 0.1

What would cause this to have a drastic influence (as would the slightly unusual tracking thresholds required) is if your response function is not calibrated to your data, resulting in FODs that are much larger / smaller than expected. However from the tracking you’re requiring a higher threshold, which suggests that the FODs are larger than expected, yet providing this option should be making things more stringent with respect to the FOD segmentation. Running tcksift2 with the -output_debug option might provide some insight into what’s going on.