SS3T-CSD b=700, (10 b0 + 30 b700, 2mm voxel), 64 slices, 30 gradient directions, 32 channels, no reverse phase, no topup

[marutimishra]

Hi,
I started with trying the new SS3T-CSD, using the pipeline https://3tissue.github.io/doc/single-subject.html, (except did not do bias field, and will be using a dilated binary brain mask derived from fsl-bbregister-mrivol2vol mapping, mrconvert to .mif, when needed)

Here's the output from my dwi2response dhollander using the voxels option

According to what I have read, red =CSF, blue= Gm and green =WM.
I don't think my data is properly picking up voxels for designated tissue types
Then I used the Tournier and here is the voxels output

It would be great if I could know, whether the data is in the right direction.
Thank you,
Maruti

[marutimishra]

Sorry for the previous multiple comments.

[thijsdhollander]

Hi @marutimishra,

No worries about the previous accidental other posts you created: I was able to delete them, to keep things clean. 👍

I started with trying the new SS3T-CSD, using the pipeline https://3tissue.github.io/doc/single-subject.html, (except did not do bias field, and will be using a dilated binary brain mask derived from fsl-bbregister-mrivol2vol mapping, mrconvert to .mif, when needed)

Yes, that's all excellent. However you obtain your (binary) brain mask is fine for sure. If needed, use a command line as mentioned in the pipeline to regrid your brain mask to the preprocessed (and potentially upsampled, if you choose to do so) dMRI data, e.g. something like

mrgrid YOUR_BINARY_MASK.nii regrid - -template dwi_denoised_unringed_preproc.mif -interp linear -datatype bit | maskfilter - median dwi_mask.mif

If your mask already "lives" on a grid that matches the dMRI data, you're fine of course. Otherwise, adapt the filenames in this example for your scenario at hand.

Here's the output from my dwi2response dhollander using the voxels option

According to what I have read, red =CSF, blue= Gm and green =WM. I don't think my data is properly picking up voxels for designated tissue types

No, it actually looks perfect. 😉 You must've picked up something wrong about the colours, I reckon. It's:

• (single-fibre) WM = blue
• GM = green
• CSF = red

So, as mentioned, your result looks very good!

Then I used the Tournier and here is the voxels output

No, that's not needed. That'll only give you a response for the non-b=0 data, and also only for WM. You can't use that for 3-tissue CSD, because it requires response functions for all 3 tissue types, and also for the b=0 data of each of those tissue types. But also, the mechanism to estimate the WM response in the new dwi2response dhollander in MRtrix3Tissue is improved in accuracy and robustness over the older dwi2response tournier approach; as shown in this recent ISMRM 2019 abstract.

So you should be good with your current dwi2response dhollander outputs. Give it a go and use them with ss3t_csd_beta1. Feel free to post the full command line output from ss3t_csd_beta1 as well as some screenshots from the resulting WM, GM and CSF maps if you want; I can check for you whether they look correct!

Cheers,
Thijs

[marutimishra]

Thank you @thijsdhollander
So I used the regrid option (-voxel 1.5) for both the DWIs data and mask image as mentioned and here is a comparative ss3T_csd output output (right is regrid)

here is the command line output

Then I ran the mtnormalize, and here is the comparative display (right is regrid-normalize)

following is the fod2dec output (regrind-right)

and the default tckgen_0.07

and at tckgen_0.09(left)

Thank you much for the feedback!!
let me know if this is ok for my data.
I propose to conduct tract estimations for FA, MD etc at different ROI's and the number of fibers in a bundle? any recommendations in this direction would be great!
Thank you,
Maruti

[thijsdhollander]

Hi Maruti,

Thanks for providing further outputs. I must say I've very happily surprised how well this has turned out, given a b-value of only b=700, and only 30 gradient directions!! 😲 👍 Sure, I've seen better SS3T-CSD results, but that's for far higher b-values and/or higher angular resolutions of course. It all looks good, given the original data quality; but I'll briefly run through and tick off each step.

So I used the regrid option (-voxel 1.5) for both the DWIs data and mask image as mentioned and here is a comparative ss3T_csd output output (right is regrid)

Very nice result, with all 3 tissue types amazingly well recovered and separated, for such a low b-value. Wow. In any case, what you can also see in that screenshot is the strong bias field across the result that you didn't attempt to correct for in earlier preprocessing steps yet. So it made sense you got your brain mask in another way (because such a bias field would've otherwise really messed up dwi2mask output). In terms of regridding: both outputs (with and without) are excellent. Whether there's much use to upsampling is still a bit under debate; I think it depends on your application. Upsampling is always nice though if you want to generate a higher-resolution crisper looking 3-tissue CSD output and e.g. tissue maps like the ones you show. It sure looks very neat with 1.5 mm voxels!

here is the command line output

All looks good. I asked that generally to see that everything would run well, but additionally also to get an (easy) peek at a few properties of your 3 response functions. In the command line output of SS3T-CSD, it reports their "SDM" ("signal decay metric") values: these should increase from WM to GM to CSF. That CSF has a far higher SDM than WM or GM is easily expected, and the output serves as a simple sanity check there. However, I was particularly curious to see whether, at a b-value as low as b=700, the new dwi2response dhollander was still able to extract a WM response with an SDM lower than that of GM. And apparently, that does still check out comfortably. This is great news in terms of robust applicability of both the new dwi2response dhollander as well as SS3T-CSD via the current ss3t_csd_beta1 implementation, to a wide range of data and data qualities. Your feedback is highly valuable here.

Then I ran the mtnormalize, and here is the comparative display (right is regrid-normalize)

Looks good again. The difference with the maps before mtnormalise is clear and striking (because the bias field was rather strong). The bias field is nicely removed by mtnormalise. While I've come to know mtnormalise as very robust already, it still depends on the quality of the tissue types at the input to allow it to do its job well. This once more confirms that the 3-tissue CSD output at b=700 is of sufficiently high quality to allow for this. 👍

following is the fod2dec output (regrid-right)

All good as with the tissue maps. Still stunned the cortical GM gets filtered out so well, while leaving the anisotropic WM part still in the WM FOD. The FOD image does show the limitations of b=700 in some places, i.e., where there is less WM, and thus less signal, and thus less SNR for the WM FOD, things start to look slightly more noisy. But in any case, it's still a strict upgrade from having only 1 or 2 tissue types: all the cortical GM that gets filtered out is a step forward. So all still good!

and the default tckgen_0.07

Yep, so here the challenges of some noisiness in smaller FODs do start to show a little bit. I think from this screenshot, it's relatively clear you want to increase the -cutoff for tractography a little bit at least. Which brings us to...

and at tckgen_0.09

...which is looking a bit better, I reckon. For subsequent targeted tractography of specific bundles, this should do fine, I think.

Thank you much for the feedback!!

Well, thank you, I should say. This is really useful for me as well in terms of seeing this perform in a real life more challenging (low b-value, low angular resolution) scenario. Very happy to provide some feedback and answers in return. 😉

let me know if this is ok for my data.

So yes, I think this is certainly ok. If you write about this in e.g. a paper, it would always still be wise to mention the low b-value and low angular resolution as "limitations": it is what it is of course. But given these data, you've got something solid to work with here nonetheless.

I propose to conduct tract estimations for FA, MD etc at different ROI's and the number of fibers in a bundle? any recommendations in this direction would be great!

So targeted tractography is your friend here: to extract tracts robustly, you almost always want to have 2 inclusion regions, constraining the tract at both "ends". You can either seed from one end-region and specify the other as an inclusion region; or do the opposite (track from the other, to the first). You can even combine both results to avoid biases (tracking performs sometimes different depending on the "direction" in combination with the shape of certain bundles). Alternatively, you can seed in the whole brain mask and specify both end-regions as 2 separate inclusion regions. That will take longer to run though (and you'll have to allow for a lot of attempted streamlines to reach a given number of accepted streamlines that connect both ROIs). To do this efficiently, you can e.g. first run whole-brain tractography (like you show in the screenshots), but with a lot of streamlines, and then use tckedit to get a particular bundle. Repeated attempts of tckedit can help you fine-tune your ROIs. Once you're certain they're fine-tuned well for a given bundle, you can in principle revert to seeding in the whole brain and including the 2 end ROIs as 2 separate inclusion ROIs. This approach is more efficient in terms of computation time, because tckedit is the convenient way to fine-tune ROIs, and you can then use them "with confidence" for the more computationally expensive "final" tckgen run, if that makes sense.

Once you've got a dense tract (many streamlines) connecting 2 ROIs, use tckmap to generate a track density map of this tract. Next, use mrthreshold with a manually specified -abs (the threshold) value to turn this into a binary tract mask; the -abs value here can be used to clean up the result to your liking. Once you've got a binary tract mask, you can use e.g. mrstats on an FA or MD map, while providing the binary tract mask to -mask, and then e.g. extract the average FA or MD value in the tract. You can also use mrdump in the same way to get a file with all FA or MD values for all voxels in the tract, "dumped" into a text file. This might be useful if you want to process all values of the whole tract in some external statistics or other package. I hope that makes sense! It's easier than it sounds; just a few commands to get used to. But once you master this process, it's incredibly useful to perform a wide range of analyses on your data.

Finally, about "the number of fibers in a bundle", I suppose you're referring to streamline count? That one's not a good idea though. The streamline count is an arbitrary (non-biological) number that is the result of the tractography parameters and settings. You could e.g. extract the value of the first volume of the WM FOD image though (so first get that volume via mrconvert, as you did to create the tissue map). That's apparent fibre density for a high b-value. For lower b-values, I tend to refer to it more as WM-likeness. You could even normalise all 3 tissues to sum voxel-wise to 1; which gives you the 3-tissue signal fractions. See this paper and this abstract / poster for 2 examples of analyses of the 3-tissue signal fractions. We've also found that the test-retest reliability and long-term stability of the 3-tissue signal fractions is very good, as shown in this paper.

Thanks again for providing feedback; and happy to help with these responses for sure!

Cheers,
Thijs

[marutimishra]

Thank you @thijsdhollander
I will apply them and check out. I wanted to understand one more detail

1. As I mentioned previously, I am not able to use topup (or just eddy for that matter) coz I do not have the reverse phase-encoded b0 files (single or few..nothing). But I do have the spin-echo field maps in both the directions. Now the issue is, when I check the EPI factor (104) and echo spacing values (0.58 ms) of these maps, they do not match to those used during the diffusion imaging data collection (0.93 ms and 128). I recall that we can use the spin-echo field maps to calculate the acq.parameter values, but just not clear if we can use them if they mismatch. I also wanted to know how important is it to use them?

Thank you again,
Maruti

[thijsdhollander]

My expertise is a bit more limited on that front, so I'll answer carefully here:

1. Depending on your application, it might not be worth taking the risk to perform a preprocessing step, in particular this kind of EPI distortion correction, if it isn't of much importance for that application. I'd say, if you don't need to get your dMRI data exactly aligned with e.g. a T1w image or something for connectomics, then just don't attempt it if you don't think you've got the proper data. It might just make things worse. If you're just after e.g. targeted tractography or some other kind of direct quantification, it's not all that important to perform this correction.

2. However, having the field maps, it sounds like you should be able to do this. Have you tried giving it a go just like that, and seeing if it improves things wrt distortions? A good feature to check is e.g. the eyeballs; they reveal the distortion or its effective correction quite well typically. However, in any case, I'd suggest to ask your question on the FSL mailing list; they've got the appropriate experts for this topic!