Estimate response function for multi-shell ex-vivo marmoset data

[Cirong_Liu]

Dear expertise of Mrtrix,

I was a older user of Mrtrix2, and used it for my previous rodent ex-vivo (single-shell) study.

Currently, I am working on a special data, a multi-shell ex-vivo marmoset data (0.15 mm isotropic):
Shell1: 6 b0 + 126 b=2400 (number of average = 1)
Shell2: 6 b0 + 126 b=4800 (number of average = 1)
Shell3: 6 b0 + 126 b=7200 (number of average = 2 for better SNR)
This data was collect for almost 1 week. All DWI have good SNR

With previous successful experience of Mrtrix on my rodents data, the Mrtrix3 become my first option to deal with this current big data.

The issue is about estimate the response function for the multi-shell dataset. The marmoset has a smooth brain with relative thin white matter and thick cortex, more like rodents than humans. Thus, the recommanded msmt_5tt (optimized for humans) is not quite fit with the animal data. After searching the post, I tested the dhollander method, but the results is not even as good as the single-shell with tournier. (see attached figure about dhollander on multi(3)-shell of data, compared with tournier method on the single-shell(7200)

Is there any way to improve the response function for multi-shell data of marmoset data?
In a previous post, the Mrtrix team is working on the infant data, using a 2 tissue approach ( Multi-shell data in neonates)

Is that available yet?

As there is only one marmoset, I think I can do a manual segmentation of the tissue types required, for example, WM, GM and part of CSF. Will that help?

(If possible, I may send a copy of a data for the testing; but the data is huge, 3.6G compressed nii.gz)

---

Thanks!

[jdtournier]

Glad to hear from long-time MRrtix users!

So the issue of getting a clean set of responses for multi-shell analysis in your case is made more complicated for two reasons:

• it’s not human, so doesn’t necessarily match many of the assumptions built into some of the algorithms, but most notably the dhollander approach, which relies on heuristics tuned specifically for in vivo adult human data (@ThijsDhollander might be able to comment further here).

• it’s ex vivo, which changes the contrast properties quite a bit, and again will completely invalidate the assumptions implicit in the dhollander method.

It’s generally accepted that when scanning ex vivo, you need b-values ~4× what you would use in vivo for roughly equivalent contrast. So your b = 7200 s/mm² is equivalent to just under b = 2000 s/mm² in vivo, which isn’t that high. So I don’t understand why the tensor model would fail this badly, either using the full set or just the single-shell variant. Any chance you could post a snapshot of the FA map? It might give use a clue as to what might be going wrong…

Otherwise, to get decent multi-shell responses, you will need to resort to a more manual approach. The simplest is probably to use the manual algorithm, feeding in masks containing voxels that you are confident are representative of ‘pure’ tissue of interest. If you use MRView’s ROI editor and select a few voxels of pure CSF and GM, and save these as CSF_voxels.mif and GM_voxels.mif respectively, you can then estimate those responses with:

dwi2response manual dwi.mif CSF_voxels.mif CSF_response.txt -lmax 0,0,0,0

and similarly for the GM. Note this assumes you’ve concatenated all of your shell data together into one big dataset using mrcat -axis 3.

The WM is a little bit more tricky in that you want to select voxels that contain pure white matter, but also that contain only a single fibre orientation. The simplest is probably to use the same trick as in the old version (and employed in the dwi2response fa algorithm): erode a brain mask, and select the N highest FA voxels within it (note this relies on having a good FA map to begin with…):

maskfilter brain_mask.mif erode -npass 5 - | mrcalc fa.mif - -mult - | mrthreshold - -top 100 WM_voxels.mif

Here I’m eroding the brain mask by 5 voxels, and selecting the top 100 highest FA voxels - you may need to adjust that to fit your data. After that, you can again use the manual algorithm to get the WM response:

dwi2response manual dwi.mif WM_voxels.mif WM_response.txt -lmax 0,4,6,8

In this case, we’re asking for a l_max = 0 for the b = 0 shell, then l_max = 4 for shell 1, l_max = 6 for shell 2, and l_max = 8 for the last shell. Again, these are numbers you can probably afford to play with - higher is probably something we’d consider safer these days, given recent improvements in the response estimation (as provided by amp2response).

Hopefully that’ll allow you to analyse your data - sounds like quite an impressive acquisition. But it would be good to get to the bottom of the issue you’re having with the highest shell…

[Cirong_Liu]

Dear J-Donald,

Thanks a lot for your detailed reply.

For the issue of b=7200, I found the reason, which seems to be caused by rescaling issue (Reco_slope_map and RG of the Bruker). The 6 b=0 images before b=7200 have rescaled by a wrong receiver gain. Previously, I thought the tensor model may be inaccurate, as the assumption of an exponential decay no longer holds for too high b-values. So I excluded b=7200 from my analysis. Now, I fixed the rescaling issue and test the fitting. I also removed it from my question, as this issue is not related from Mrtrix3,

For the main issue “response functions”, I will test out your suggestion. Once it works well, I will post the testing results on the post. I think people working on multi-shell data of rodents and marmosets may face with the same problem.

[jdtournier]

@Cirong_Liu

The 6 b=0 images before b=7200 have rescaled by a wrong receiver gain.

Yes, that would explain it! Glad you got it figure out.

@Cirong_Liu

For the main issue “response functions”, I will test out your suggestion. Once it works well, I will post the testing results on the post. I think people working on multi-shell data of rodents and marmosets may face with the same problem.

Thanks, I’m sure that would be hugely appreciated!

[Cirong_Liu]

Hi, J-Donald,

I have done several tests on the data with your suggestion.

It did a surprisingly good job in segmenting the gray matter and white matter of the marmoset brain. See the thin white matter in the temporal-pole area. Previously, I have tested traditional segmentation tools (fast, ants) Neither of them fully segment-out this thin white matter using MTR (T1) and T2 ex-vivo image. But, the FOD of mrtrix3 did!

---

--------------What I did

_ To generate response functions, here was what I did to select voxels:_

a. White matter - I used two methods. One was as you suggested, I manually chose voxels with highest FA (top 300). Second is using single shell (7200) and applying Tournier method to obtain the white matter voxels (300), and then used the voxels for manual multi-shell dwi2response.

(the blue is manual and the red is Tournier)

b. Gray matter - The marmoset cortex is too thick and contain many voxels. Thus, I manually selected about 1500 voxels with FA lower than 0.1. Note that with the high-resolution ex-vivo data and relative thick cortex of marmoset, we can see the different layers of the cortex that have different FA value and features. This is different from what we have on human invivo data. I am not sure whether this simple method is appropriate or not.

c. CSF - The ventricle of marmoset is not that big. So I manually selected 150 voxels of the CSF from S0 image.

Then I did what you suggested (lmax = 0,4,6,8 for the white matter of different shells; lmax = 0,0,0,0 for gm&csf ). Besides the encouraging results, there are still issues that I am not quite confident:

---

--------------Issues

1. Issue of b=7200 again:

(from left to right: b0 (lmax0), b2400 (lmax4), b4800 (lmax 6), b7200 (lmax8))

Why the radius of circles of b4800 is smaller than b2400, but b7200 is larger than all lower b value? (the response value is listed in the next question)

Is this normal?

As I mentioned previously, once b=7200 involved, the FA of many white matter will >1. I have used the right Receiver Gain to adjust the data, still the same results.

I then consulted expertises in DTI, and identified the reason:
"(Niethammer, et al, 2006, On Diffusion Tensor Estimation) Simple tensor fitting could give rise negative eigenvalues: Because dMRI signal can be really small at b=7200 when G is parallel to a fiber, the fitting can be very sensitive to noise when noise level is in fact above the minimum signal level. The formulation of FA makes it larger than 1 when some EV is negative. ".

Now, for the Mrtrix and the CSD tracking, etc, will it be affected by the issue as the tensor model?

---

2. Manually selected vs. tournier:

The Tournier selected many different voxels than the FA-manual method and resulting different response functions (for 300 white matter voxels).

(left: manual; right:tournier)

The touriner has smaller (radius) but flatter shape than the manual.

Which one is better?

**_Manual:_**
_173670.8311645914 0 0 0 0_
_112471.0979004641 -30904.53437063165 6914.149390457013 0 0_
_84556.50252458362 -30573.53020595047 12082.37401380716 -3937.906991237937 0_
_146785.459380438 -59106.35423326885 26799.40599680597 -11027.75365967198 2882.521653928761_

**_Touriner:_**
_322084.1663155202 0 0 0 0_
_183785.3993442675 -59540.47745553248 13320.75581448077 0 0_
_122420.6147702928 -53843.52800533723 23477.13836021282 -8166.972327980991 0_
_212011.0697191145 -100392.5019609229 50932.57994245004 -23922.30182768965 6450.961159453199_

---

3. About FOD of gray matter and CSF

For tckgen, only the white matter FOD are needed.

What is the usage of FOD of gm and csf?

---

4. Track Orientation Density Imaging (TODI)

I have read the paper about the TODI and also watched the concept video on youtube about the TODI.

The super-resolution image generated will be extremely useful for animal data with thin white matter, including rodents and marmosets. I want follow the concept in the TODI paper and try to implement it on the marmoset data. However, the discussion about implementation I in the paper is general.

Are there any guidelines and suggestions for implementing the new method?

There is a command, tckmap, able to generate old TDI and the new TODI (a lmax can be defined).

Should I use higher lmax for dwi2response and dwi2fod (in the paper, lmax upto 16 was tested for TODI).

Thanks in advance!

[jdtournier]

Nice results! To answer some of your questions:

Yes, you definitely still have scaling issues with that b = 7200 s/mm² shell. If you look at your response coefficients files, the first column corresponds to the mean DW signal within each shell. You can see that the signal goes down as expected, until it hits the last shell, which seems to be roughly twice the expected amplitude. I can’t think of any reasonable explanation for that other than a scaling issue…

Thankfully, MT-CSD is immune to this - since the response files are calibrated directly from the data, it will adjust to such differences in scaling, as long as they happen only between shells. It’s not ideal since the noise will scale along with it, so that last shell might introduce more noise than the others. But other than that, it’ll happily deal with your data.

The tensor model will on the other hand be sensitive to this. I’m not sure how badly that would affect your FA values, but it might be worth not using the last shell for the DTI fit, or better, fix the scaling issue.

Looking at your data, there doesn’t seem to be much point in modelling CSF, since there’s basically none to be seen - in your b=0 image, the ventricles are essentially background, not the expected high amplitude (~free water) signal. This makes sense given that these data were acquired ex vivo… So in your case, I don’t think there’s any point in modelling CSF at all, you can just do a 2 tissue MT-CSD decomposition, using only GM & WM.

As to what is the use of the GM & CSF - well, the main point is to allocate the right signal to the right tissue type. In the regular single-shell CSD, there is no mechanism to prevent the signal from these other tissue types from contaminating the WM ODFs. In MT-CSD, this means that the WM ODFs are much cleaner than they would otherwise be, because the contaminating signal from other tissue types can be allocated to other components in the output. So that’s probably the main motivation here, but there may be lots of other uses for these other components that we haven’t thought of…

One more thought: if you’re using the latest version of MRtrix (3.0_RC1), you can actually use a higher l_max value nowadays than would have been recommended in the past - especially given that you have enough directions for l_max = 14… I’d recommend you increase the value for the WM ODF to at least 10, it might give you slightly cleaner results - although admittedly not massively so.

As to the big difference between the manual and tournier methods, as far as I can tell, it’s mostly a difference in scaling (?). That probably won’t make a great deal of difference to the output, although it will give you slightly smaller ODFs for the tournier method. How big an impact this might have will depend on what you plan to do with these data. But given the scaling issues that are clearly still problematic, I’d be slightly sceptical of the manual results. And even if the scaling wasn’t a problem, we have found that the FA-based method does tend to select different voxels, typically where the b=0 is abnormally low (due to e.g. Gibbs ringing) since that tends to increase FA, while the tournier method is explicitly biased towards high intensity voxels (or rather, high WM ODF amplitude voxels), since these are more likely to represent a full voxel’s worth of WM.

Hope this helps.

[Cirong_Liu]

Hi, J-Donald,

Thanks for your reply. I am much clear now.

The scaling issue is weird since the data was directly converted from 2dseq raw data and the RG and slope were directly read from the raw data (acq and method from the scanner).

The difference between b7200 and other shells (b=0,2400,4800), is that the b7200 have the number of average = 2, and the others are 1. I think this setting may cause the scaling issue. Most possible reason is that when the Bruker acquired the data with 2 averages, it only sum-up the values but didn’t divide the value by 2.

Your answers reminded me of checking the noise.

I calcuated the mean and the std of the background values (noise) for different shells using a mask shown in the red:

One DWI image:
b0 = (mean) 15139.66 +/-  (std)6974.30
b2400 =     14833.24 +/-       6991.54    
b4800 =     14724.05  +/-      6944.20
b7200 =     20905.63 +/-       9887.80

From this table, we can see:

1. The std of noise (also the denominator term for SNR) is not changed across different shells (b0, b2400,b4800)
2. The b7200 indeed has a higher noise scaling (about 1.42x) than the other shells.

Since the b7200 has 2 averages, it should have √2 = 1.42x less noise std than the other shells.
If the scanner did not divide the value by 2, the noise level of b7200 would be 2/√2 =1.42x higher noise than other shells, consistent with the observation.

Based on the above estimation, there are two methods to rescale the b7200.

1. Dividing by 1.42x: the b7200 will have the same noise level as the other shells. However, since the b2700 have two average, this may be not the real noise level. After rescaling, the FA of b7200 is reduced with few voxels > 1, but the overall FA level is still larger than the b2400&b4800.

2. Dividing by 2 : this may be the real noise and signal level for b7200, but in this case, the noise level of b7200 will be different from other shells. After rescaling, The FA of b7200 come to a similar level to b4800!

Which method should I use to rescale the b7200 for Mrtrix(csd)-based tracking?

Should I keep the same noise level across shells, or keep the normal fitting (real noise and signal level of b7200)?

If the later case, should I perform dwidenoise directly on all shells, or perform it seperately on different shells?
(After dividing by 2, the noise level of b7200 will be different from the other shells. However, from my understanding, the PCA will also normalize the data. Thus it should be OK to performed dwidenoise directly on all shells?)

The data will be used for tractography and reconstructing all white matter tracts in marmosets. Meanwhile, I also hope to use it to construct a super-resolution TDI image to help me manually drew or delineated some small white matter ROIs, and also for qualitative visualization.

There will be NO quantitative statistically comparison (for example between two groups). Thus, the real value of the fitting mapping won’t be a big issue for me.

Sorry to bother the issue that less related to Mrtrix.

–

Another question is about dwi2fod, if I want to set the lmax = 12 for example.

dwi2response manual dwi.mif wm.nii.gz wm.txt -lmax 0,8,10,12 
dwi2response manual dwi.mif gm.nii.gz gm.txt -lmax 0,8,10,12 

dwi2fod msmt_csd dwi.mif wm.txt fod_wm.mif gm.txt fod_gm.mif -lmax ???? -mask mask.mif

Should I set lmax of dwi2fod to 12,0 or 12,12, or anything else?

Thanks!

[jdtournier]

OK, that all makes sense. Ideally, you wouldn’t have averaged the two instances of the b7200 shell, and simply concatenated them into one larger shell. This probably doesn’t matter a great deal in your case since these are ex vivo data, but for in vivo data, keeping them separate allows motion correction to perform properly - you can’t undo the motion artefacts after averaging.

The only remaining issue is the √2 noise scaling, which I don’t think will influence things a great deal. But to answer your questions:

@Cirong_Liu

Based on the above estimation, there are two methods to rescale the b7200.

1. Dividing by 1.42x: the b7200 will have the same noise level as the other shells. However, since the b2700 have two average, this may be not the real noise level. After rescaling, the FA of b7200 is reduced with few voxels > 1, but the overall FA level is still larger than the b2400&b4800.
2. Dividing by 2 : this may be the real noise and signal level for b7200, but in this case, the noise level of b7200 will be different from other shells. After rescaling, The FA of b7200 come to a similar level to b4800!

Which method should I use to rescale the b7200 for Mrtrix(csd)-based tracking?

I’d definitely go with option 2 here. Much better to have the correct signal, even if its noise level differs from the other shells.

@Cirong_Liu

Should I keep the same noise level across shells, or keep the normal fitting (real noise and signal level of b7200)?

I’m not sure what you mean here… But I’m assuming it’s a rephrasing of your previous question…?

@Cirong_Liu

If the later case, should I perform dwidenoise directly on all shells, or perform it seperately on different shells?

That’s an interesting problem… dwidenoise expects constant noise level, at least within its default 5×5×5 neighbourhood patch. It will clearly not perform properly on your data if the noise level differs. From this point of view, you might be better off first dividing the b7200 shell by √2 to ensure the same noise level, running dwidenoise, then extracting out the b7200 shell and dividing by √2 again to ensure the correct signal level. You should be able to achieve all this with combinations of mrcat, mrcalc, and dwiextract, e.g. something like (untested):

$ mrcalc b7200.mif 2 -sqrt -div b7200_sqrt2.mif
$ mrcat b0.mif b2400.mif b4800.mif b7200_sqrt2.mif -axis 3 dwi_constant_noise.mif
$ dwidenoise dwi_constant_noise.mif dwi_denoised.mif
$ dwiextract dwi_denoised.mif -shell 7200 b7200_denoised.mif
$ mrcalc b7200_denoised.mif 2 -sqrt -div b7200_denoised_sqrt2.mif
$ dwiextract dwi_denoised.mif -shell 0,2400,4800 other_shells.mif
$ mrcat other_shells.mif b7200_denoised_sqrt2.mif -axis 3 dwi_final.mif

---

@Cirong_Liu

Should I set lmax of dwi2fod to 12,0 or 12,12, or anything else?

I reckon this is something you’ll need to experiment with a bit. The important thing is to get the best response possible, and @rsmith found that it is possible to get non-negligible l_max=10 terms from regular human in vivo data with the new amp2response. Given your data, you might be able to get meaningful values up to l_max=12. But this doesn’t mean you have to use that l_max in the final dwi2fod reconstruction - you’ll have a more accurate response, so the fODFs are likely to be cleaner too (although that will admittedly be a minor difference). For human in vivo data, I tend to think l_max=8 is about the best we can hope for (and is the default), but you may find you can push your analysis to l_max=10 or 12 in your data. I’m not sure how best to evaluate this though, maybe inspect the tractography or TDI output visually? It’s a tricky thing to assess, using l_max=12 will give you sharper fODFs, but that’s likely to come at the expense of less sensitivity to sharp curvature, with bending fibres coming out as two distinct peaks, for example (disclaimer: this is my gut feeling - no hard data). So there will no doubt be difference, but I’m not sure how to tell which is best… Personally, I’d probably stick to l_max=8 just to play it safe…

[Cirong_Liu]

Hi, J-Donald,

Many thanks for your reply and suggestions. I will reprocess the data as you suggested.

Actually, I did a test on 8,12,14 on several well-known tracts, the final tractography is almost similar at least visually (although maybe some quantitative difference).

I will stick to the default lmax = 8 for the later-on analyese. lmax=8 is safter and also faster than higher lmax (for dwi2fod).

Thanks again!

[ThijsDhollander]

It’s late-to-the-party-day, but here we go anyway :

@Cirong_Liu

Currently, I am working on a special data, a multi-shell ex-vivo marmoset data (0.15 mm isotropic):
Shell1: 6 b0 + 126 b=2400 (number of average = 1)
Shell2: 6 b0 + 126 b=4800 (number of average = 1)
Shell3: 6 b0 + 126 b=7200 (number of average = 2 for better SNR)
This data was collect for almost 1 week. All DWI have good SNR

@jdtournier

it’s not human, so doesn’t necessarily match many of the assumptions built into some of the algorithms, but most notably the dhollander approach, which relies on heuristics tuned specifically for in vivo adult human data (@ThijsDhollander might be able to comment further here).

@jdtournier

it’s ex vivo, which changes the contrast properties quite a bit, and again will completely invalidate the assumptions implicit in the dhollander method.

So well, nope, not really a problem after all. The basic strategy and relative tissue properties it relies upon still hold for a marmoset. The challenge with these kinds of data is probably rather the ex-vivo aspect versus the maximum b-value you’ve been able to acquire. This indeed:

@jdtournier

It’s generally accepted that when scanning ex vivo, you need b-values ~4× what you would use in vivo for roughly equivalent contrast. So your b = 7200 s/mm² is equivalent to just under b = 2000 s/mm² in vivo, which isn’t that high.

I’ve heard similar numbers, varying from a factor 3 to 4 indeed. But it also indeed seems that this was at least a bit anticipated here as well. If this lands you with an “in-vivo equivalent” b-value of about 2000, you should be good. What may be needed though, is to tune the -fa parameter of dwi2response dhollander, as I advised here as well: Multi-tissue CSD . A rodent (though in-vivo) is quite a step further away from a human-like brain compared to a marmoset… and even there, no worries given an appropriate tuning of -fa.

@jdtournier

Otherwise, to get decent multi-shell responses, you will need to resort to a more manual approach.

Sometimes not a bad idea indeed, although I’ve found it’s actually surprisingly hard to really select good voxels manually. Even single-fibre ones: you may aim for traditional single-fibre regions, only to select voxels that actually have quite some dispersion (though no “crossing”). The tournier algorithm avoids these consistently (and rightly so).

@Cirong_Liu

The 6 b=0 images before b=7200 have rescaled by a wrong receiver gain. Previously, I thought the tensor model may be inaccurate, as the assumption of an exponential decay no longer holds for too high b-values. So I excluded b=7200 from my analysis. Now, I fixed the rescaling issue and test the fitting. I also removed it from my question, as this issue is not related from Mrtrix3,

Too bad… But this was crucial to dwi2response dhollander working correctly though as well: an initial FA threshold (i.e. that -fa option) is quite an important initial step. If that’s messed up by a faulty FA map, the algorithm may be set up for failure indeed. On the other hand, it’s a pity you now had to exclude this shell, because the highest b-value is definitely the most valuable one for most, if not all, CSD variants.

@Cirong_Liu

It did a surprisingly good job in segmenting the gray matter and white matter of the marmoset brain.

The result looks great though! The reason it works so well compared to traditional segmentation methods is that multi-tissue CSD is actually not a segmentation method (with e.g. spatial priors and based of just 1 or a few intensities): the signal within each voxel is pulled apart into tissue signal contributions independently. So it doesn’t matter how thin any structure is (spatially): if the tissue contribution to the signal in those voxels is there, it should be able to tease it out.

@jdtournier

Looking at your data, there doesn’t seem to be much point in modelling CSF, since there’s basically none to be seen - in your b=0 image, the ventricles are essentially background, not the expected high amplitude (~free water) signal. This makes sense given that these data were acquired ex vivo… So in your case, I don’t think there’s any point in modelling CSF at all, you can just do a 2 tissue MT-CSD decomposition, using only GM & WM.

Makes sense indeed, but may depend on what the ex-vivo brain was suspended in, or even how it was prepared. But I have to agree, judging from the depicted results, it looks ok even if it would just be modelled with WM and GM alone.

[Cirong_Liu]

Hi, Dhollander,

Thanks a lot for your reply! Things become clearer with your comments.

Right now, we are pushing the limits of the MRI scanner and tried to obtain an 80um marmoset diffusion MRI data.

I will test your suggestion on the new data.

Thanks again!

[ThijsDhollander]

@Cirong_Liu

Right now, we are pushing the limits of the MRI scanner and tried to obtain an 80um marmoset diffusion MRI data.

Wow, if that works out well, that would be some very nice data!

[ThijsDhollander]

For future reference, besides the aforementioned rodent data, here’s another example of a (quite high quality) ex-vivo monkey dataset where dwi2response dhollander worked out of the box without issues: Problem with dwi2response msmt_5tt (amp2response) command for ex-vivo monkey data