Response function

1

Dear MrTrix experts,

I am using mrtrix on DWI data 60 directions b=1500.
I am experiencing some problems with spherical deconvolution. I have performed both tensor- and SD-based whole brain tractography and looked at some major associative tracts. However it seems like the tensor works much better in reconstructing these tracts, somehow against my expectation.
So I was wondering wether maybe something went wrong in my response function estimation (I used dwi2response). I compared this with the one Donald posted and they don’t look identical. These are the values of my response: 379.9197083 -95.88195801 22.71349716 -4.001269341 0.5601699948.
I attach the response function and the single fibre mask.
I would very much appreciate if you could give me any advise on how to check on this and be sure this is acceptable.

Best

Chiara

2

This definitely doesn’t look right… Can I ask which version of dwi2response you used for this? What does dwi2response -version report? If it gives a sensible version number, then you’re still using the old version, which has since been completely overhauled as part of the big update on March 12. The issues are now explained in detail in the documentation.

On the other hand, if you were running the latest version of MRtrix3, then this is definitely something that we’d need to look into…

3

dwi2response 0.3.12-1027-g14…i guess the old version…
How wrong is this? I mean, how much does this will affect tractography results?
Is there a way I can fix this?

Best

Chiara

4

Well, given how bad your results seem to be, I’d say the old version can be pretty wrong… The proper way to fix this would be to upgrade your installation and try again with the new dwi2response script - but bear in mind that the upgrade isn’t quite as straightforward as a regular update, as explained in the relevant post. If you don’t want to interfere with your current installation, you might find that the best course of action is to install the new version side-by-side, and see what difference it makes.

Brief instructions for a side-by-side install:

$ git clone https://github.com/MRtrix3/mrtrix3.git ~/mrtrix3-updated
$ cd ~/mrtrix3-updated
$ ./configure
$ ./build

and to use executables from the updated install, simply set the PATH temporarily for that shell session:

$ export $PATH=~/mrtrix3-updated/release/bin:~/mrtrix3-updated/scripts:$PATH

If you subsequently decide to fully switch over to the updated version, you can delete your current mrtrix3 folder, and rename the mrtrix3-updated to mrtrix3 - and you’ll also need to change your PATH to match the new location for the executables. Basically, open your ~/.bashrc in a text editor (e.g. gedit ~/.bashrc), and modify the line that sets the PATH for mrtrix3 to read:

export PATH=~/mrtrix3/release/bin:~/mrtrix3/scripts:$PATH

i.e. add release between mrtrix3 and bin…

5

Yes, I suspect this is the sort of GM-like ‘single-fibre voxels’ selection @ThijsDhollander has described to me before.

If you’re absolutely stuck using the software version you have, I would suggest either providing a brain mask to the dwi2response command if you aren’t already, or eroding the mask prior to running dwi2response if you’re already providing one. But the recommended tournier algorithm in the newer dwi2response script does appear to be more robust, particularly in the case of lower b-values, so a full software update would be very much preferable.

6

Yep, it clearly is the now more frequently reported common thing that happens every so often with the Tax algorithm…

Another further suggestion in case you wouldn’t be able to upgrade for some reason: go even a bit more strict, and don’t just give it an (eroded) brain mask, but even an (eroded) FA > 0.5 mask or something. The moment the Tax algorithm has access to a reasonable number of either crossing WM, or GM voxels, this outcome can pop up almost seemingly random. It’s of course not really random, but the algorithm’s very first iteration / initialisation is very unstable, and all the rest simply depends on its outcome. It’ll either go for the true single fibre WM voxels, or specifically for almost anything but the single fibre WM voxels.

7

Thanks for your suggestions,
Yes I am already providing an eroded mask to dwi2response. Also, I have analysed different type of data from different scanners, and the response functions look like this one, though more disc-shaped for high b-values. Is it possible I made some mistake in the preprocessing, or is it more likely that dwi2response is constantly failing?
I followed the old dwi tutorial: /maskfilter mask_mask.nii.gz erode - -npass 6 | ./dwi2response $input response.txt -fslgrad $bvec $bval -mask - -sf sf_mask.nii.

Chiara

8

Is it possible I made some mistake in the preprocessing, or is it more likely that dwi2response is constantly failing?

It can be difficult to assess whether or not the response function estimation has ‘succeeded’ from the response function alone, particularly if you’re analysing data acquired using different protocols. Really the best indicator is the output single-fibre voxel selection mask, and it is a subjective assessment: it is up to you to decide whether or not the voxel selection is appropriate.

I’ve used the pre-FA-threshold trick that Thijs suggested myself in the past and it helped with the data I was dealing with at the time. So you would do something like this:

dwi2tensor $input - -fslgrad $bvec $bval -mask mask_mask.nii.gz | tensor2metric - -fa fa.nii
maskfilter mask_mask.nii.gz erode - -npass 6 | mrcalc fa.nii 0.5 -gt - -mult - | dwi2response $input response.txt -fslgrad $bvec $bval -mask - -sf sf_mask.nii
9

Thanks for your suggestions.
This is the response function I get after applying a binary FA mask thresholded at 0.8 (0.5 was not improving results) on the same data (60 directions, b=1500 s/mm2).
(9 iterations, 302 voxels). (405.3344727 -176.5801544 57.78819275 -11.31396675 1.342096567). I believe this confirms that the issue relates to the error in selecting SF voxels.

10

OK, I guess this makes sense - although it’s not ideal. I’m assuming this is still using the old dwi2response executable? Clearly it’s not doing a great job… which is why we’ve completely overhauled it with the current script. Personally, I’d strongly recommend you update and use the new version, it does seem more stable than the previous approach…

But at least you seem to have confirmed where the problem was. Thanks for updating this thread with your findings, I’m sure it’ll be helpful for others facing similar issues…

11

This is the response function I get after applying a binary FA mask thresholded at 0.8 (0.5 was not improving results) on the same data (60 directions, b=1500 s/mm2).

By the time you’ve applied an FA threshold as high as 0.8, the dwi2response command probably has very few voxels left to process and remove from the single-fibre mask, and becomes relatively redundant; the single-fibre voxel selection becomes dominated by the FA thresholding. So the old dwi2response is really not working full-stop on these data. I hope there isn’t something over and above the inadequacies of the deprecated dwi2response binary going on here… Re-running with the tournier algorithm in the new dwi2response script would be informative in that regard.

12

This is the response function obtained with the new turnier script. (430.9638763 -190.6219556 60.72299934 -11.87131925 1.331388084)
How would you define an ideal response function? And how can I check that?
Why was the previous one non ideal? and Do you consider this ideal instead?
Also values look pretty similar to me. Is there an optimal range for these values?

Robert, what do you mean by something over and above the dwi2response?

thanks
Chiara

13

How would you define an ideal response function? And how can I check that?

If we knew the ideal response function, we wouldn’t need to calculate it every time :stuck_out_tongue:. As mentioned previously, the best sanity check is the single-fibre voxel selection mask image.

Also values look pretty similar to me. Is there an optimal range for these values?

The magnitude of the values will depend on your acquisition. The magnitude of the first term (l=0) is reflective of the typical signal intensities in your image; some scanners provide images with typical intensity in the white matter of ~300, others may be around 1e7. The magnitude of the higher harmonic degree coefficients should decrease monotonically from this, and be negative for every second value. The rate at which these values decrease is also influenced by the b-value; a higher b-value gives greater angular contrast, and hence greater power in the higher harmonic degrees.

Robert, what do you mean by something over and above the dwi2response?

I was surprised by the fact that the dwi2response binary was failing even when provided with an initial mask of FA > 0.5. So I was concerned that there was something else wrong with your data. But I don’t think that’s the case given what the newer script has provided you with.

14

Hi experts

following up on conversation, I am performing tractography on the ODF obtained from the new response function (Tournier).
To recap data: b=1500, 2x2x2 mm3, 60 DW directions. I applied eddy current distortion correction and bvec rotation in FSL.
I would like to perform whole brain CSD deterministic tractography, and dissect some major tracts. However results still look pretty bad to me. In general I have very noisy reconstructions, with short truncated streamlines.
I attach a screenshot of the ODF on coronal view.
Tensorial reconstructions look ok.

I would appreciate any hint on how I can improve tractography or make a quality check to understand what is going wrong.

Best

Chiara

15

Indeed, this is not looking great - although it’s hard to tell with such a small image - maybe you can upload a larger or zoomed-in version? It’s hard to see what’s actually going on at the voxel level…

In terms of what might actually be going wrong, the first thing to check is your gradient directions. Post the output of the following commands:

mrinfo dwi.mif -export_grad_mrtrix enc.b
dirstat enc.b

Note that if your DWI data are stored in a non-MRtrix format (e.g. NIfTI), you’ll need to add the -fslgrad bvecs bvals option to the mrinfo command above.

This will output some summary information about the DW gradient directions used during the acquisition.

16

thanks for your fast reply
Here is the output and the bigger figure

b = 0 [ 10 directions ]

  Bipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 180, range [ 180 - 180 ]
    energy: total = inf, mean = inf, range [ 1.79769e+308 - inf ]

  Unipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 180, range [ 180 - 180 ]
    energy: total = inf, mean = inf, range [ 1.79769e+308 - inf ]

  Spherical Harmonic fit:
    condition numbers for lmax = 2 -> 2: nan


b = 1500 [ 60 directions ]

  Bipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 5.75748, range [ 0.500071 - 11.816 ]
    energy: total = 45056.4, mean = 750.94, range [ 140.738 - 13240.1 ]

  Unipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 26.5455, range [ 25.4903 - 27.6894 ]
    energy: total = 3387.96, mean = 56.4661, range [ 55.776 - 57.4813 ]

  Spherical Harmonic fit:
    condition numbers for lmax = 2 -> 8: 1.00899 1.02505 1.20932 28.8272

FOD.jpeg661×556 446 KB

17

Thanks for the more detailed snapshot - much easier to see how bad things are. It seems to be noisier on the left of the image (subject’s right), is this a consistent observation…?

Regardless, the output of dirstat confirms what I suspected: your DW gradient directions are badly distributed. You can ignore the b=0 output section, it is meaningless for HARDI. For the b=1500 section, this is the interesting section:

  Bipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 5.75748, range [ 0.500071 - 11.816 ]

This states that your directions are on average 5.7° away from their nearest neighbour, and in the worst case, 0.5° apart (!) - this is a very poorly designed gradient scheme…

The reason for the scheme being so poor is probably due to the procedure used to generate it being based on a unipolar electrostatic repulsion model, rather than the required bipolar model (see e.g. Jones 1999 for details). This seems to be what happened here, since the numbers look close to optimal for the unipolar repulsion model:

  Unipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 26.5455, range [ 25.4903 - 27.6894 ]

Unfortunately, that’s not what you want for optimal HARDI… The last line shows that your scheme supports a lmax=6 fit (condition number = 1.2), but completely collapses for lmax = 8 (condition number = 28.8).

For reference, this is what you should get with a properly distributed 60-direction scheme (generated using dirgen):

$ dirgen 60 dir60.txt
dirgen: [done] Optimising directions (power 2, energy: 4001.0599, gradient: 0.00011143311, iteration 361)
$ dirstat dir60.txt 
dir60.txt [ 60 directions ]

  Bipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 18.8232, range [ 18.302 - 19.5847 ]
    energy: total = 8002.12, mean = 133.369, range [ 132.208 - 135.388 ]

  Unipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 19.4942, range [ 18.302 - 30.4278 ]
    energy: total = 3794.16, mean = 63.236, range [ 46.0022 - 84.591 ]

  Spherical Harmonic fit:
    condition numbers for lmax = 2 -> 8: 1.00612 1.01637 1.03319 1.15016

Note that the average nearest-neighbour angle in the bipolar section (the important one here) is ~18°, with a tight range around that value. And the condition numbers stay close to one for all lmax values, which indicates that the fits will be stable up to those values.

As to what you can do about this… Not a great deal unfortunately. The obvious thing to try is to restrict the processing to lmax = 6 - or at the very least, restrict the response function estimation to lmax 6 (which I think it should do by default? How many values do you end up with in your response file? It should be 4 for lmax=6, 5 for lmax=8). I think you’d be safe to run dwi2fod at lmax 8 as long as your response function is estimated up to lmax 6 only - but your best bet is to try it and see…

You will also want to double-check where your data came from, and make sure you fix the acquisition going forwards. In the old days, the MGH works-in-progress package installed on many Siemens scanners suffered from this problem. GE scanners also had such schemes installed on their systems for some time. Philips scanners had other issues if you used their ‘overplus’ option - but I don’t think this is the problem here. And of course, if you’re running a custom sequence or scheme, this could also have introduced the issue…

1 Like
18

Thanks for this useful information. I have a few follow up questions:

  • I do not fully understand what you mean by bipolar or unipolar electrostatic repulsion models. I read the Jones 1999 paper. They do not seem to explicity refer to bipolar or unipolar models, they describe what I understood as the electrostatic repulsion model (for explicit reference I copy below the specific part of the paper I am refering to). Is this description what you are calling the bipolar model? If yes, what would be the definition of the unipolar model?

  • A further literature search on this topic in the literature led me to discussions of the diffusion gradient pulse polarities in the diffusion MRI sequence (eg. same polarities around the spin-echo refocusing pulse, or bipolar pulses during the first half prior to the refocusing pulse). However, I assume this is unrelated to the gradient direction scheme of this discussion.

  • If I have a 60 direction scheme that seems to have been created to satisfy the uniplar model (as you helped me understand), how would you suggest to proceed with “trimming” this to a subset of orientations that brings the subset closer an optimal bipolar model?

Best,
Chiara

PS: From Jones et al., 1999 (page 517, section Gradient Vector Orientations for Estimating the Trace)
“Consider a model in which a line parallel to each gradient vector passes through the center of a sphere, and a
unit electrical charge is placed at both of the points where the line intersects the surface of the sphere. Each gradient vector is represented by a pair of points in this way because a diffusion attenuation measurement with the gradient in a positive direction could equally well have been performed with the gradient in the opposite direction, since the diffusion tensor is symmetric in the absence of charged moieties. The repulsive force between a pair of charges is, according to Coulomb’s law, inversely proportional to the square of the distance between the charges. The algorithm used to arrange the gradient vectors uniformly in 3-dimensional space therefore adjusts the orientations of the sets of orthogonal gradient vectors until the sum of the repulsive forces between every possible pair of charges is minimized.”

19

Yes, that’s exactly what I meant.

Simply that instead of having a pair of points at opposite ends of the line, you would have only one point at one end of the line. What this means example is that for the case of two directions, the unipolar model gives opposite directions as the optimal solution, whereas the bipolar model gives directions at right angles - very different answers…

Yes, that’s another issue altogether…

I’m not sure what you’re asking, and why, but assuming you actually want to pick out a subset of directions that is closer to optimally distributed than your original unipolar repulsion set, the procedure might be to work out for each direction, the bipolar electrostatic repulsion energy with all of the other directions, remove the one with the highest energy, and repeat until the remaining set is ‘good enough’. However, it won’t ever be as good as the equivalent number of properly distributed direction (i.e. via bipolar electrostatic repulsion), and there is no clear point at which the set becomes ‘good enough’, so it’s not ideal.

If you want to split your direction set into two sets that are each near-optimal, there is actually a command that does something like this: dirsplit. It’s a brute-force randomisation approach (not guaranteed to find the absolute optimal solution), but will probably give you something useable.

But this all depends on why you want to do this. If you want to improve your CSD reconstructions, I’d recommend you don’t reduce the direction set, since that will only reduce the SNR. What you can try instead is to reduce the lmax used for the response estimation step to something that your direction set support - so basically add the -lmax 6 option to dwi2response. The CSD process (as implemented in dwi2fod) is itself quite robust to these effects, as long as the response is reasonable.