I’m trying to follow the code for performing fixel-based analysis (FBA).
This is for DWI data of b=1000, 32-directional, and no reverse phase-encoding readout.
I note that within the related papers to both FBA and AFD, for each subject the spherical harmonic images are to be generated for lmax=4, for enabling co-registration to a population template.
Prior to this step I generated the response functions with the dhollander option and lmax values set to 0,4.
However, when running dwi2fod msmt_csd -lmax 0,4 -shell 0,1000 with the response functions from WM, GM, CSF etc - I receive the following error: “Number of lmaxes specified does not match number of tissues”
If I let the algorithm also go to default and do not specify the lmax, it rather runs smoothly (I’ve passed the -info and -debug options).
I’ve checked everything ranging from the embedded gradient direction files, and the response functions for each tissue appear to be in order. This makes me believe it may be a bug… For an example subject:
If you don’t have multi-shell data you only need a single WM response function file. Note that we will be adding a new method to perform multi-shell CSD using single shell data in the very near future (@ThijsDhollander :nudge: :nudge:). However until then you will be limted to only performing single tissue CSD.
For single tissue CSD you should use the tournier method in dwi2response, since it will only give you the b=1000 shell (and not the b=0).
Note that the FBA instructions still recommend using the dwi2fod msmt_csd algorithm even for single shell data, since it uses a hard non-negativity constraint. It’s a bit of a hack (feeding single shell data into the msmt_csd algorithm), and at some point we will probably make the csd algorithm option use a hard constraint too. I know this is a little confusing. I hope to improve the FBA documentation in the near future.
mrregister, which is called inside the population_template script, will automatically use the lower degree coefficients (up to lmax 4). So just compute FODs at lmax 8. Therefore you should not need to specify either the lmax or the shells.
In dwi2response (particularly in the dhollander algorithm), the -lmax option indicates the maximum spherical harmonic degree per b-value shell, for the white matter response function. (The GM and CSF responses are assumed to be isotropic for all b-value shells in this algorithm)
In dwi2fod, the -lmax option specifies the maximum spherical harmonic degree of the estimated FOD for each tissue.
This is why, when you specify ‘-lmax 0,4’, dwi2response proceeds quite happily (since there are two b-value ‘shells’ in your data when you include the b=0 images), but dwi2fod fails (since you’re providing three tissue response functions, and therefore requesting three output FOD images, but have only specified two values for lmax).
As Dave said, there’s a more fundamental problem with trying to estimate three tissues with only two b-value shells using the stock-standard MSMT CSD algorithm; I just wanted to clarify why the reported error appears where it does.
I should of made clearer that I’m unable to perform the co-registrations between the diffusion and T1-images (no reverse b0’s were acquired). Hence, my enthusiasm for testing out dhollander…
I guess this actually means that nevertheless I cannot benefit from the constraint of dwi2fod msmt_csd , because I won’t have accurate response fxns of the different tissue types??
Meaning I should (for this dataset) simply be using dwi2response tournier and dwi2fod csd until the new algorithms come in.
I have a doubt about the acquisition and the method.
With this b-value is it possible to perform this kind of analysis? If I understand correctly, one of the conditions of the paper is b values around 3000. Is that correct?
Good question. Personally, I don’t see any reason why you can’t do a fixel-based analysis on regular b=1000 data. What it does mean though is that the simple interpretation of your fibre density values as being proportional to the intra-axonal volume won’t be as valid. The extra-axonal signal won’t be attenuated to negligible levels, so it will also contribute to the overall signal, and hence to the apparent fibre density. Note that some people would argue that we need b>6000 to fully attenuate the extra-axonal compartment, in which case this confound will also be relevant for b=3000, although to a much lower extent.
On the upside, that signal only corresponds to ~20℅ by volume in healthy white matter, and it will still be more strongly attenuated than the intra-axonal signal, so its contribution to the overall DW signal should be relatively minor.
So in my opinion, the downside of lower b-values (other than a minor loss of angular resolution) is that any differences observed in a fixel-based analysis are slightly harder to interpret: the difference may be related to changes in the extra-axonal space, rather than to changes in the intra-axonal volume… This is problematic because we’d expect intra-axonal volume to correlate very strongly with fibre density, so introducing a confound with the extra-axonal signal dilutes the ‘purity’ of that interpretation. But all that means is that you have to be a bit careful not to over-interpret changes you might observe in your analysis…
I should of made clearer that I’m unable to perform the co-registrations between the diffusion and T1-images (no reverse b0’s were acquired). Hence, my enthusiasm for testing out dhollander…
No stress, that was perfectly clear. It just wipes dwi2response msmt_5tt off the list of options.
I guess this actually means that nevertheless I cannot benefit from the constraint of dwi2fod msmt_csd, because I won’t have accurate response fxns of the different tissue types?
You can still benefit from the hard non-negativity constraint that’s intrinsic to the msmt algorithm, as opposed to the soft non-negativity constraint in the csd algorithm; this is what @Dave was alluding to. What you can’t benefit from at present is the multi-tissue capability of the msmt algorithm, since you need at least 3 b-values to disentangle the three primary brain ‘tissue’ types. (Technically you could use two response functions to estimate two tissues during deconvolution; but we haven’t properly demonstrated how & why this could / should be used, particularly since this concept will be likely superseded by the 2S3T algorithm anyways.)
Meaning I should (for this dataset) simply be using dwi2response tournier and dwi2fod csd until the new algorithms come in.
You can definitely still use dwi2fod msmt_csd as mentioned already. You could technically use dwi2response dhollander, and then simply extract the b=1000 shell of the WM response function to use in deconvolution: I’d expect this to behave more or less like deriving an approximate WM mask using FA>0.2 and using this as the input mask to dwi2response tournier (since that’s close to what’s going on internally anyways), @ThijsDhollander might have additional details on the pragmatic performance of WM response function estimation between the two. One possible advantage of sticking to dwi2response tournier for now is that it’s published and citable - might be weird trying to justify in a manuscript why an unpublished multi-shell response function estimation technique was used for single-shell deconvolution.
dwi2fod: [DEBUG] Shells: b = { 0(1) 1000.0108224759756(32) }
dwi2fod: [DEBUG] loading matrix file "rsp.txt"...
dwi2fod: [DEBUG] found 1x4 matrix in file "rsp.txt"
dwi2fod: [ERROR] number of rows in response function must match number of b-value shells```
Obviously I would generate a 2 x 4 matrix function to circumvent this, but I am unable to with the ```tournier``` algorithm.
I have even passed through the -shell 1000 and -lmax 6 options
Cheers,
Alistair
Yes, we know this is a bit confusing currently. The problem is that dwi2fod tries to use the b=0 image here as well, due to the msmt_csd algorithm, but you only offered a response that has no b=0 part in it, due to the tournier algorithm. What you are trying to do here is essentially single-shell single-tissue CSD, but just using the msmt_csd algorithm’s code. That is ok, there’s nothing wrong with that. But, the msmt_csd algorithm then needs to only get the single-shell part of your data, and not the b=0 image(s) as well. To resolve this for your example here, replace your step 2 by:
Note that if you do this, you best don’t describe this in a paper as if having done MSMT-CSD; because as mentioned, you’ve only been doing single-shell single-tissue CSD, but while exploiting the benefit of the “hard constraint” in the msmt_csd code.
Again, I/we fully agree this is confusing at the moment; this will eventually become much more clear when the need for a single-shell single-tissue “exploitation” of the msmt_csd algorithm will disappear.
Also, I will post something else below to give you another good option for your scenario: doing MSMT-CSD with 2 tissue types (WM and CSF) on your data. Bear with me for a moment, as I gather some text from elsewhere to explain this…
There we are with a second reply. I’ve actually gotten a very similar question (as compared to your original question) from another user; but they accidentally sent it to me in a private message; so I’ll be doing an effort to quote the relevant bits for you. This user was also in a scenario of a low b-value and single-shell data (with the presence of a b=0 image).
Here goes part of one of my answers on parts of their question:
They then went on to ask me, among other questions, these things:
To which I then gave this list of answers. Keep in mind that they also found themselves in a scenario with a low b-value (55 directions at b=1500 in their case).
They then went on to try these 2 scenarios (as quoted from the user):
Note both are valid options: the first one is using tournier for the WM response and csd for single-shell single-tissue CSD, the second one is using dhollander for WM GM CSF responses and msmt_csd for multi-tissue CSD with only WM and CSF (i.e. ditching the GM response at this stage).
The user then went on to show me screenshots of both results. Both looked good, but the second one definitely improved over the first one in 2 ways:
Benefiting of the hard-constraint in the (deep) WM.
Benefit of the multi-tissue aspect, notably the inclusion of CSF in the model, as a reasonable free-water removal mechanism. Essentially makes all false positive “WM FODs” in the CSF disappear, and improves the quality of the WM FODs when partial voluming with extra-cellular free water happens.
As mentioned, the extra benefit of such free-water elimination is actually more relevant for lower b-values. I’ll ask “the user” (who is @stellamarissanchez89) if they want to post those (or other) screen shots of their outcome here as well. It may indeed help a few other users too.
To clarify: the first one is indeed tournier, but followed by csd (single-shell single-tissue), whereas the second one is indeed dhollander, but followed by msmt_csd using only the WM and CSF responses (i.e. ditching the GM one).
Interesting idea to include the CSF response in dwi2fod msmt_csd to remove false positive WM FODs even in case of single shell (+b0) data. Would you advise to use this approach as well in case of extensive white matter lesions? Or would there be a chance that no FOD remains in the lesion voxels?
I’m convinced it will never hurt to do so, and improve your results slightly anyway, so I would certainly advise it. However, don’t expect too much of it though: given a realistic CSF response, it will never introduce drastic changes to the FODs in white matter hyperintense lesions, since the CSF response at the b≠0 shell should only show very little signal. Or to put in in another way: it’ll only fully remove FODs of the size/kind you’d see in the actual CSF (in the ventricles) when doing simple single-shell single-tissue CSD.
I’ve got an ISMRM abstract accepted that relates very closely to this topic (the one I also mentioned in this and following posts). As mentioned over there, I’ll make the text available, e.g. on ResearchGate, ASAP after the proceedings are available for all attendees.
Another side note: you may still see your FODs “visually” (when viewing in mrview) more drastically changing in size when switching from single-shell single-tissue CSD to MSMT-CSD (e.g. with WM and CSF), but those changes will be visually prevalent over all (healthy) WM. This has nothing to do with the multi-tissue aspect, but rather with the hard non-negativity constraint of the MSMT-CSD implementation (whereas the original (SSST-)CSD implementation has a “soft” constraint/regularisation). The fact that this makes your FODs look smaller, is a bit misleading: they actually also become slightly wider, so it’s just the peak amplitude that becomes smaller. The lobes’ integrals will generally be very similar in magnitude.
Finally, as mentioned a few times in this thread above as well, in the future you’ll be able to use single-shell 3-tissue (SS3T) CSD as an option too, which can give you the full 3 tissue types from just single-shell + b=0 data. This will work better with higher b-values (e.g. 2500-3000 and up), and depends more crucially on very decent SNR for lower b-values. However, “in the future” is still a little bit of a wait for now… as I’m sadly occupied with way to many different things at the moment.
Just to be clear: I fully agree with @bjeurissen here, but it’s probably not very likely to happen. If it would happen, the “lesion” you’re looking at essentially has the properties of a (100%) CSF-like fluid. In that case, there’s indeed no evidence of any WM left any more (in the measured data), so the WM FOD should optimally be zero indeed. However, in that case, you would probably not be calling it a “white matter hyperintensity” (WMH), since it would appear pitch black on a FLAIR. A “fluid filled cavity” or something would probably better describe it then.
(you can download a PDF version of the whole thing from there, or log in to ResearchGate for a more pleasant browsing experience)
In context of your question, @tbilliet: for typical white matter lesions (hyperintensities on a FLAIR) of the kind that you’d encounter in multiple sclerosis or in an elderly population (whether healthy or suffering from mild cognitive impairment or Alzheimer’s disease, etc…), you’d typically see reduced WM AFD (FOD amplitude) combined with both CSF-like and GM-like isotropic signal contributions. The amounts of these may vary depending on the exact lesion content, but it’s quite interesting to see that the ways in which they seem to do so make quite good sense with respect to what may have effectively caused them.
If you were to work with only WM and CSF responses, you’d still see a clear CSF-like compartment appearing, but it would only clean up the high b-value signal a very tiny bit; the remainder (of GM-like signal) would still mostly end up in the WM FOD. Given what has been observed micro-structurally in many of those lesions, gliosis is quite a sensible explanation for such signals. This would mean the apparent fibre density from single-tissue CSD (or 2-tissue WM-CSF CSD) may show only a moderate (to close to none) effect: intra-axonal space reducing, yet being taken up by a lot of other mess (e.g. astrocytes) that still result in a decent chunk of DWI signal, which ends up partially compensating the AFD decrease that you would expect. Introducing the GM response resolves a lot, if not very close to all, of this. The remaining FODs’ shapes make good sense with respect to “what was there before the damage”, and their sizes have decreased according to the damage done.
