Dwi2fod Multishell "single-tissue" versus "multi-tissue" CSD

Hello MRTrix team!

My question is similar to this thread, Multishell "single tissue"

I’m working with multi-shell DWI data.

  5 volumes b=0
64 volumes b=1000
64 volumes b=2000
64 volumes b=3000

When calling the dwi2fod function, I am using the msmt_csd algorithm.

I originally used this function:
dwi2fod -nthreads 16 msmt_csd in_wm.mif group_rfunc_wm.txt WM_FOD.mif

Performing the msmt_csd with 3 tissue types is (understandably) taking a much greater amount of processing time per subject. Will the WM_FOD output differ considerably if I perform the FOD estimation using 3 tissue types instead of the original run using only 1 tissue type (WM)?. dwi2fod -nthreads 16 msmt_csd in.mif group_rfunc_wm.txt WM_FOD.mif group_rfunc_gm.txt GM_FOD group_rfunc_csf.txt CSF_FOD.mif

As a follow-up question, will the impact be greater if I am subsequently looking at FOD values in a major WM region (e.g. corpus callosum) versus grey matter (e.g. putamen)?

NOTE: the dwi2response text files were generated using the dhollander algorithm for all 3 tissue types (WM, GM, CSF).

Thanks so much for your time!

Yes, it is expected to be considerably different. I strongly advise against using single-tissue CSD with multi-shell data: your fODFs will actually degrade compared to single-tissue CSD with single-shell data and you will also not see any of the benefits of multi-compartment modelling as shown in Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data - PubMed.

If you have 16 threads at your disposal, the processing time of dwi2fod with three tissue types can’t be that long? Especially if you compare it to steps like motion and eddy current correction, dwi2fod running time should be minor (see MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation - ScienceDirect, for typical runtimes). What is your spatial resolution?

Thanks for your response @bjeurissen !

The native resolution of the diffusion data are 2mm isotropic.
The upsampled data, which are what is being fed into dwi2fod, have a 1 mm isotropic resolution (256 x 256 x 138 dimensions).

I performed the DWI upsampling because it is part of the multi-tissue CSD tutorial ( Fibre density and cross-section - Multi-tissue CSD — MRtrix 3.0 documentation )

The DWI pre-processing runs pretty quickly on 16 threads. I’ve been running a single test subject on the 3 tissue type MSMT-CSD at this resolution, and dwi2fod has been running for approximately 20 hours, and is at 75%. Is it possible that the multi-threading command isn’t being recognized as part of the command syntax?

I am using the most recent version of MRTrix (3.02).


Follow up - I think the problem may relate to the group-average response function.

I tried re-running the same test-subject but using the 3 tissue response functions derived from that same subject for the dwi2fod estimation (command syntax otherwise identical), instead of the group average response functions (derived from average response functions of ~100 other subs). The command ran in about 15 minutes.

Some of the early subjects in the group cohort were processed through a prior version of MRTrix, but all using the dhollander algorithm in dwi2response. Do you think it would be a good idea to rerun the dwi2reseponse for the entire cohort, to make sure they are all being calculated similarly?

Hi Brad,

20 hours is way too long for dwi2fod, even with upsampled data. What would likely need to be happening for this to be the case is for the combination of response functions and DWI data in specific voxels to lead to some form of oscillatory numerical behaviour in the MSMT CSD solver such that it never converges on an FOD solution. It might be interesting for us to be able to see such data to figure out why it’s happening and if anything can be done to improve the robustness of the algorithm.

Given you observed much faster execution using subject-specific response functions, I would suggest comparing the magnitudes and shapes between the two RF sets. E.g. You may have gross differences in intensity scaling between that one subject and the typical image intensities of other subjects.

Some of the early subjects in the group cohort were processed through a prior version of MRtrix, but all using the dhollander algorithm in dwi2response. Do you think it would be a good idea to rerun the dwi2reseponse for the entire cohort, to make sure they are all being calculated similarly?

It would depend on precisely what version of MRtrix3 was used. If it was 3.0.x, then behaviour should be identical; macro updates are only to fix outright bugs and do not change experimental outcomes. The behaviour of specifically the single-fibre WM response function estimation within that algorithm changed between 3.0_RC3 and 3.0.0, so it would be preferable to not use a mixture of those for different subjects.


Thanks for your response @rsmith ! The “early subjects” I mentioned were processed quite a while ago. The version was 3.0_RC2.

Here is the WM response function values from one of those early subjects processed in version 3.0_RC2:

3447.654903262014 0 0 0 0 0
1790.478494717137 -725.5523712037685 163.315124601483 -24.25065540038764 3.130741083344601 -3.09210997343932
1222.44035206819 -706.5014045902841 287.9967951221154 -76.98730773871866 14.9746210225946 -3.593769646578895
976.983091766167 -611.0780390873934 321.6086147846045 -122.5392324291043 35.45435459651641 -8.58570665697982

And here is the same subject when I have re-run dwi2response (dhollander):

# Shells: 0,1000,2000,3000
# command_history: /mrtrix3/bin/amp2response dwi.mif voxels_sfwm.mif safe_vecs.mif response_sfwm.txt -shells '0,1000,2000,3000' -nthreads 16  (version=3.0.2-82-g722d1587)
1363.729113846975 0 0 0 0 0
760.970284411094 -308.0323775782157 68.78931720135283 -9.461863534657233 1.252176481878456 -1.184952131029491
526.0207221262441 -297.655258739154 119.7274140042323 -31.53698799974157 5.55916343247834 -1.047882864659559
422.8809691475108 -255.6856858253959 133.0442433715916 -49.12942009497527 13.81093235209071 -2.971642148784646

These magnitude of the values are quite different, but their shapes using shview are similar.

Old WM values


Recalculated WM values

I imagine that dwi2fod is having trouble resolving the FOD estimation given a group average response function similar to the top set of values using the current version of MRtrix. Looks like re-running the dwi2response and dwi2fod to “update” all subjects might be the best thing.