Anatomically constrained tractography using ODF from BEDPOSTX

Dear experts,

I was wondering whether it is possible to use Anatomically constrained tractography framework (ACT) with ODF created by BEDPOSTX. Is it sensible? How to do that?

I have rather low-orientation-resolution DWI data: 30 directions (half-sphere direction coding), 2 averages, b=900. The dwi2response automatically lowers lmax to 4, estimated coefficients are 593.64, -63.89, 6.82. The response function does not look very good:

I thought whether BEDPOSTX could be more suitable with these data than CSD. But I would like to use ACT and other mrtrix tools for connectivity analysis. What do you think?

My first question is which version of dwi2response are you using? We recently replaced the previous executable with a more flexible script (full details here) with lots more options, and the default parameters should work better than that…

However, the fact that the algorithm reverted to a lower lmax does suggest that your sampling scheme is far from optimal. What scanner/sequence were these data acquired with?

Just out of interest, have you actually tried using for dwi2fod…? It might not behave too horribly - although you’d really need to fix that response function to get something more sensible.

The second issue is whether BEDPOSTX would be expected to produce better results. I’m definitely the wrong person to ask… I don’t think it would do much better than CSD (assuming your response function is sensible), but others might disagree. Converting the output of BEDPOSTX into an FOD format suitable for use in tckgen is not impossible, but it would be a bit of work - and it’s not something I feel we need to implement any time soon…

The more direct use of bedpostx output in tckgen would be:

  • Obtain the fibre directions - dyads<i> files from the looks of things
  • Manually rotate these directions from image space into scanner space (well, I’m assuming these are provided in image space)
  • Concatenate the images to form a single 4D volume, with 3N volumes (where N is the maximum number of fibres in a voxel)
  • Use tckgen -algorithm fact. This expects the image intensities to represent X-Y-Z components of a fibre direction in each triplet of volumes.

However this would only be a deterministic algorithm; it would only use the mean of the samples calculated rather than the full distribution.

But as Donald said, it’d be better to figure out why dwi2response is giving that result. I would look very closely at your pre-processing. In cases where the gradient table is incorrectly oriented, the resulting response function will often look like this, so that’s where I’d look first; there’s been a couple of bugs in various iterations of pre-processing scripts that could potentially have caused this kind of corruption. The best way to have a look is to run tckgen -algorithm tensor_det on your data before and after pre-processing, and see if the tracks go in the wrong directions in the latter.

…this should actually not have any effect on the response function estimation, I reckon.

Nonetheless, definitely use the latest version; I have a very strong feeling this may be an outcome of the previous dwi2response command rather than the current script. Especially on low b-value and only average angular resulution data, I’ve seen it behave like this myself.

Dear all,
thank you for all your suggestions.

My DWI was acquired as follows:
30 directions, 2 averages for each direction, b=900, voxel size 2mm3. Siemens Trio 3T, VB17. It is standard Siemens DWI sequence. The diffusion sensitizing directions are distributed on half-sphere, principially like on the right panel at the fsl-eddy page:

http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/EDDY

Since I cannot use eddy for data preprocessing (“half-sphere” diffusion encoding), I preprocessed the data by FLIRT and epi unwarp (using fieldmap) by script obtained from here: http://www.bic.mni.mcgill.ca/~thayashi/dti.html
I did not rotate b-vecs and not averaged the acquisitions, only stacked them together.

In my previous post I was indeed using older version of dwi2response.
After updating mrtrix, I tried FA, Tournier and Tax algorithms for dwi2response. The worst results I obtained by Tax algorithm. The response function in Tax looks much like in my previous post. The best shape I got with Tournier algorithm, which is slightly better than FA algorithm.

There are my response functions:
Tax, coef 589.16 -83.61 12.07

FA, coef 549.73 -164.22 28.18

Tournier, coef 716.43 -244.36 50.91

I also tried to determine response function from raw uncorrected data and by Tournier algorithm I get almost the same results:
706.65 -256.54 52.69

Both Tournier response functions (for uncorrected and corrected data) look optically the same. It seems that there isn’t problem with my preprocessing, however, the preprocessing makes much of improvement in this case.

Here is a comparison of FOD generated by:

  1. old dwi2response
  2. new dwi2response - Tournier method

    There is clearly much better angular resolution with new dwi2response with Tournier method. I suppose that the tractography results could be reasonable now.

I would like to ask:
1.
I was wondering, why optically the amplitudes of FOD lobes were higher in old dwi2response in comparison to new dwi2response (to show approximately similar amplitudes of FOD lobes, my screenshots were generated for old dwi2response with scale 2, but for new dwi2response I had to use scale 3). Could you comment on? Or is it expected behaviour?

Is here any objective evaluation criterion of response function quality? I suppose the more higher value of lmax, the better. But what in case of 2 responses with identical lmax? The higher value of SH lmax coefficient, the better?

Does the “half-sphere” diffusion scheme make any problems/consequences for response function determination / FOD determination by constrained spherical deconvolution?

From your CSD papers I understood that your acquisition schemes are “full-sphere” acquisitions. But as the diffusion process is symmetric, the “half-sphere” acquisitions could be theoretically of benefit for having better angular resolution (not speaking now about the disadvantage that “half-sphere” data could not be corrected by fsl-eddy). Could you comment on?

Does automatic selection of lmax for my “half-sphere” data works OK?
According to your paper
NMR Biomed. 2013; 26: 1775–1786
you are starting to able to estimate lmax=6 with 28 DWI directions.

In case of more averages, do you recommend to average data, or prefer to stack them together? What is better option for response function / FOD determination? Does it matter?

Concerning my original query about BEDPOSTX, thank you for the replies and suggestions. Indeed I would need to use probabilistic anatomically-constrained tractography (ACT), therefore deterministic FACT algorithm is not much of interest for me. I would therefore stick to CSD for FOD determination. But indeed I think it could be useful to have the option to use BEDPOSTX for ACT, to be able to compare both FOD determination methods in terms of performance and usability in different raw data acquisition setups.

Since I cannot use eddy for data preprocessing (“half-sphere” diffusion encoding)

Although this may not permit complete correction of eddy current distortions in eddy, it wouldn’t surprise me if it still performs better than a linear registration approach, as it will generate appropriate target volumes for realigning each diffusion direction independently. I’d suggest at least giving it a test.

In my previous post I was indeed using older version of dwi2response.

Given the new results, I think your gradient table is OK; I’d suggest that it is indeed the nature of the Tax algorithm combined with your low b-value that’s causing the highly isotropic responses.

This approach initially treats every voxel as a single-fibre voxel; crossing fibre voxels are iteratively removed from the mask as the response function is progressively sharpened and multiple FOD peaks are identified. However, the lower angular contrast in low b-value data makes these crossings harder to identify. Consequently, many genuine crossing-fibre voxels will not be removed from the single-fibre mask. When the signal profiles within this mask are then averaged to produce the response function, the profile is not as sharp as it should be: maximal sharpness arises from selection of pure single-fibre voxels and alignment of their fibre orientations, so the presence of signal from crossing fibres inevitably broadens this shape.

  1. I was wondering, why optically the amplitudes of FOD lobes were higher in old dwi2response in comparison to new dwi2response

Yes, this is expected given the difference in mechanisms between the tax and tournier algorithms. In the tax method, all voxels that are not removed from the single-fibre mask contribute to the response function estimation; so the overall response function amplitude will be an average of that in all single-fibre voxels. In the tournier method, a cost function is defined to quantify the ‘single-fibre-ness’ of each voxel, and only the top 300 voxels are averaged to produce the response function. Since the cost function includes the amplitude of the primary fibre peak, the algorithm preferentially uses voxels where the single-fibre FOD amplitude is largest (typically the splenium). Therefore the response function from this method will have a larger overall amplitude. Since the amplitude of the FOD reflects the magnitude of the diffusion signal relative to the response function, and the response function resulting from this method is greater, the resulting FOD amplitudes will be smaller.

2.Is here any objective evaluation criterion of response function quality? I suppose the more higher value of lmax, the better. But what in case of 2 responses with identical lmax? The higher value of SH lmax coefficient, the better?

I’m not sure if there’s a black-and-white answer to any of these; @jdtournier might pitch in at some point and has more experience than I do. But as an example: if you have a lower b-value, the signal power in the higher harmonic degrees is smaller; when deconvolution is then applied (at least in the linear case), there’s a reciprocal relationship between the power in the response function and the power in the resulting FOD; so there’s a risk of introducing large noisy coefficients into the higher harmonic degrees of the FOD compared to if you had truncated the response function at a lower harmonic degree. But the presence of the non-negativity constraint complicates this relationship.

In general, sharper is considered better, as long as that sharpness is derived from the data itself. For instance, it’s possible to define the thinnest possible ‘pancake’ response function and use that; but the resulting FODs tend to be very noisy.

3.Does the “half-sphere” diffusion scheme make any problems/consequences for response function determination / FOD determination by constrained spherical deconvolution?

4.Does automatic selection of lmax for my “half-sphere” data works OK?

No and yes. All of the mathematical transformations between amplitudes on the sphere and spherical harmonics in MRtrix3 assume antipodal symmetry, and so will work with direction schemes defined either on the sphere or the half-sphere.

However, to be clear:

From your CSD papers I understood that your acquisition schemes are “full-sphere” acquisitions. But as the diffusion process is symmetric, the “half-sphere” acquisitions could be theoretically of benefit for having better angular resolution (not speaking now about the disadvantage that “half-sphere” data could not be corrected by fsl-eddy). Could you comment on?

A full-sphere acquisition does not imply acquisition of volumes where the diffusion directions are perfectly opposite to one another, and in fact this is usually explicitly avoided. A full-sphere acquisition is simply one where the directions of diffusion sensitisation are not constrained to the half-sphere. Typical derivation of such a scheme would involve optimising a set of directions on the half-sphere (assuming e.g. a unipolar electrostatic repulsion model; e.g. dirgen), and subsequently reversing the polarity of around half of those directions so that they are spread relatively homogeneously over the full sphere (e.g. dirflip). So there is no loss in angular resolution in using a full-scheme acquisition; the benefit is better modelling of eddy currents, and in fact long-term eddy currents may be reduced.

5.In case of more averages, do you recommend to average data, or prefer to stack them together?

Though it probably doesn’t make much difference in practise currently, conceptually (and for better compatibility with future developments) stacking volumes with equal diffusion sensitisation directions is preferable. Averaging magnitude-transformed volumes alters the noise distribution from Rician to non-stationary non-central chi, and so would not perform optimally in any technique that assumes a Rician distribution. However it’s even more preferable to acquire more unique directions.

In the past, having multiple image volumes with equivalent diffusion sensitisation directions could cause issues as the automatic selection of lmax was based on the number of diffusion volumes only; however we now have a more robust heuristic in place that should deal with such data more appropriately.

6… But indeed I think it could be useful to have the option to use BEDPOSTX for ACT, to be able to compare both FOD determination methods in terms of performance and usability in different raw data acquisition setups.

I’ve added an entry to the GitHub issue tracker regarding this feature. It’s at a relatively low state of priority right now given the other issues we need to address in order to properly release the software, but I am intending on implementing a few different tracking algorithms at some point for my own evaluation purposes, so there’s a chance I might be able to implement this during my own research, and later add to MRtrix3.

Cheers
Rob

@rsmith has covered everything already, but just thought I’d comment on a couple of aspects.

The model in dirgen is a bipolar electrostatic repulsion model (by default). A unipolar model gives terrible direction schemes, since it totally ignores antipodal symmetry.

@rsmith is correct that the Tournier approach does select the larger amplitude voxels preferentially (more likely to consist of pure WM, and higher SNR), but I’m not sure that’s sufficient to explain the magnitude of the effect you’re showing. The first term in the responses you show only differ by ~30% (these terms correspond to the mean amplitude of the response function). I’m not sure I completely understand why the difference would appear to be so great, but maybe the fact that the FODs come out so much broader also makes them look larger…?

Thank you for the suggestion. From the eddy documentation I initially supposed that eddy is not suitable at all for the “half-sphere” data. However, thoroughful study of recent eddy paper


showed me that the main potential problem is that there would remain some significant residual eddy current effect due to in average (from eddy current perspective) non-symmetrical diffusion sensitization gradients in “half-sphere” acquisitions, which can be faced by using second level eddy current modelling. Therefore, I tried eddy on my 2x30 directions “half-sphere” data and it indeed worked very well. I found that:

  • The movement correction worked perfectly, much better than eddy_correct
  • There is almost no difference in data between using or not using second level modelling

From this I conjecture that the eddy current effect in these DWI data is very small.
This lead me to conclude that eddy can be very well used for preprocessing of half-sphere data with 30 directions. As this 30 direction “half-sphere” data have been acquired by standard Trio 3T DWI VB17 protocol, I hope that this info could be useful for anyone who acquired data this way and now is looking for optimal data preprocessing.

1 Like

The model in dirgen is a bipolar electrostatic repulsion model (by default).

Derp. I knew what I was talking about; just got the words mixed up…

From this I conjecture that the eddy current effect in these DWI data is very small.

Are you using a twice-refocusing diffusion sequence protocol? In such cases the residual eddy current distortions can be sub-voxel at low resolutions, in which case fitting non-linear eddy current distortions would be difficult even with a full-sphere acquisition.