fixel2voxel: how to extract primary, secondary, and tertiary voxel-based FD maps

Hi,

I am trying to convert fixel data such as FD, FC and FDC into voxel-based metrics. I believe the fixel2voxel command is best suited for that purpose. The complexity metric is interesting to me as well and this paper:


explains its meaning well.
For the other metrics I am unsure what makes the most sense. I am considering using the sum option for FD, FC (log) and FDC but debating whether it is more appropriate to use the weighted option. Are there some recommendations depending on the type of study? We are planning to look at multiple metrics in a large group of individuals.
I have not found much on the topic in the literature. What I have found in literature, in this paper:

is that they are using the primary, secondary, and tertiary voxel-based FD maps, apparently extracting these with fixel2voxel, but I cannot find how to get these metrics from the Mrtrix documentation.
Is it done using the -none option and -number 3 (to get the 3 largest fiber bundles)? I read in a previous post on this forum that the peaks are actually ordered by peak amplitude and not by fiber density. Has this been changed?

To summarize here are my 2 questions:

  • If/when to use the weighted option of fixel2voxel
  • If FD can be extracted for the 3 main peaks (ordered in terms of FD)

Thank you,

Stefanie

Hi Stefanie,

Okay, there’s a few things to cover here… <stretches>

The complexity metric is interesting to me as well

Free tip: If investigating the complexity metric, I would advocate using the -fmls_no_thresholds option in fod2fixel. Whether or not an FOD lobe other than the largest one does or does not survive thresholding has a drastic effect on this metric. By preserving all positive FOD lobes, no such thresholding effects arise. These manifest as erroneous-looking sudden transitions in the value of complexity between adjacent voxels, and it wouldn’t surprise me if it also falsely increases variance within a cohort.

I am considering using the sum option for FD, FC (log) and FDC but debating whether it is more appropriate to use the weighted option. Are there some recommendations depending on the type of study?

I would say that the recommendations would depend moreso on the physical nature of those metrics rather than any details of the study.

FD is a measure of volume; while it’s unitless due to the way we typically perform spherical deconvolution, it can nevertheless be thought of as having dimensions L3. Now if you want a voxel-wise measure, probably the most sensible is “the total fibre density”, which would be a sum. Note however that rather than going through the process of FOD segmentation, you could just take the l=0 term of the WM FOD spherical harmonic expansion, and that would give you the same measure but without annoying fixel thresholding effects; see also this manuscript on voxel-wise total fibre density. But that’s not to say that this is the only statistic across fixels that can be taken: selecting the maximal FD across fixels in each voxel gives you the maximal FD in each voxel (funnily enough); the mean FD is probably not a good choice, as it will be strongly dependent on the FOD segmentation thresholds as to how many smaller FOD lobes get included and hence drag that average down.

FC is interesting. If you want a voxel-wise measure of morphological change, you could just take the determinant of the Jacobian of the non-linear warp field. This would be ignoring fibre orientation information entirely; which is of lower complexity, but would give a different result to FC in a voxel with a single fixel only. So it’s possible to do a bit better. But here’s where the -weighted option comes into play.

Imagine a voxel with two fixels A and B, with FD values of 0.9 and 0.1 respectively. The non-linear warp field is such that the values of log(FC) for these fixels are -0.5 and 0.5 respectively. Now you want to collapse all of this information into a single scalar value of log(FC) for the voxel. What should that value be? If you use fixel2voxel to take a mean, you’ll get 0.0. If you provide the FD values via the -weighted option, you’ll get -0.4.

Now imagine that fixel C just survives fixel segmentation, with a value of FD of 0.01 and a value of log(FC) of 0.5. With the first method, the voxel mean changes from 0.0 to 0.167; with the second method, it changes from -0.4 to -0.395.

So weighting the averaging process by the fixel fibre densities not only gives greater precedence to larger fixels, but also mitigates the influence of FOD segmentation thresholding.

Now back to your original text. Let’s say you have a voxel with three fixels, each with log(FC) = -0.5. Now you take the sum within the voxel, to produce a value of log(FC) = -1.5. Does that value make any physical sense? FC encodes an expansion / contraction; yet here the result is suggesting that the expansion / contraction of the voxel is three times as severe as that of any individual fixel contained within…

they are using the primary, secondary, and tertiary voxel-based FD maps, apparently extracting these with fixel2voxel, but I cannot find how to get these metrics from the Mrtrix documentation.
Is it done using the -none option and -number 3 (to get the 3 largest fiber bundles)?

Yes, that will give you a 4D image where each of the 3 voxels contains the FD value of the first / second / third tallest FOD peak.

I read in a previous post on this forum that the peaks are actually ordered by peak amplitude and not by fiber density. Has this been changed?

Correct; and no, it has not (yet) been changed. The ordering by peak amplitude is a natural consequence of the way in which the FOD segmentation algorithm works. However given that we advocate use of the FD metric rather than peak amplitude, it probably does make more sense for them to be in the fixel directory output in order of decreasing FD… Any objections to changing this in the next minor update?

That’s not to say that there’s anything stopping you from interacting with the fixel data in your own code, performing whatever reordering or other manipulations that you wish; the format is fully documented. This decision purely influences whether or not certain intuitions or assumptions about the data hold up when performing existing data manipulations.

Cheers
Rob

1 Like

Thank you @rsmith for the very detailed response.

I will definitely implement the -fmls_no_thresholds option in fod2fixel, it indeed makes a lot of sense to compute complexity. I imagine this will also impact the other metrics (sum FC and FD voxel-wise). Any reasons why we would not want to include smaller lobes for these other metrics?

For FD, you recommended using the l=0 term of the WM FOD spherical harmonic expansion. This corresponds to the first volume of the WM FOD image, correct?

For FC, I am still very unsure which measure makes most sense physiologically. When computing FC, we are using the warps (which were used to warp fixels into a common “template fixel space”, shrinking or expanding each fixel in order to align it with the template fixel). So the FC metric is in fact a measure of how much smaller or how much bigger a fixel is compared to the “norm” (is the size of the template fixel pre-determined or is it based on the average fixel size?). Essentially, I would want a metric that is the equivalent (or as close as possible) to summing absolute cross-sections values (in micrometer squared let’s say) of all fixels in a voxel, if that’s at all possible.

Lastly, I did not really understand what the advantage of taking the determinant of the Jacobian of the non-linear warp field would be. Could I just take the value of the subject2template.mif image file at a given voxel as an indication of the total deformation that was applied to the fixels contained in that voxel? I apologize for the perhaps naive question, I am not used to working with warps for anything else than applying to an image in registration and I don’t fully understand what the “determinant of the Jacobian” means.

I greatly appreciate your help,

Stefanie

Hi @rsmith,

I hope you are doing well.

I was wondering if you had a chance to look into this. I would like to get a clearer idea of how to set up my voxel-based “fixel metrics” and ensure I fully understand what they represent before going forward.

Thank you,

Stefanie

Hi Stefanie,

Unfortunately a small amount of time invested in my own work puts me well and truly behind on helping everybody else :scream:

I imagine this will also impact the other metrics (sum FC and FD voxel-wise). Any reasons why we would not want to include smaller lobes for these other metrics?

Probably not, I suppose.

Assuming true non-negativity of the FOD1, with no segmentation threshold, the voxel-wise FD sum should approach the value of the l=0 term of the spherical harmonic expansion (with a sqrt(4pi) scaling factor), whereas any thresholding imposed will cause a deviation from this value, with the magnitude of the deviation depending on the interaction between those thresholds and the particular shape of the FOD. So fixel thresholding would probably just contribute variance. But having said that, if you can derive such from the FOD l=0 term directly, and can avoid fixel segmentation entirely, it becomes a bit of a moot point.

For FC, it would be fine as long as when you compute the voxel-wise measure you take the weighted mean, with the magnitude of each fixel’s contribution being determined by its FD. So smaller fixels that get included due to the absence of segmentation thresholds don’t actually influence the voxel-wise measure very much.

Indeed thinking about it, it should be mathematically possible to derive a voxel-wise aggregate measure of “FC” that also doesn’t rely on fixel segmentation. I’ll likely get the linear algebra wrong linguistically; but conceptually, it’s basically taking the local Jacobian in each voxel, but instead of calculating FC along a finite set of fixel directions and combining such, or just taking the determinant (which is agnostic to fibre direction), instead take some form of inner product (?) between the FOD and the Jacobian to produce a voxel-aggregate FC. I imagine it would be something comparable to what fod2dec does. Hopefully someone understands what I mean…

For FD, you recommended using the l=0 term of the WM FOD spherical harmonic expansion. This corresponds to the first volume of the WM FOD image, correct?

Correct.

For FC, I am still very unsure which measure makes most sense physiologically.

What makes most sense to me for voxel aggregation is a weighted mean, with the weights determined by the fixel FDs. What you’re looking for is “among those fibres in this location, what is the typical / expected expansion / contraction orthogonal to the fibre direction necessary to align to the template?”. The magnitude of expansion / contraction depends on the orientation of the fibre; but if there are more fibres in that voxel in direction A than direction B, then the FC value derived for direction A should have more of an influence on definition of that typical / expected magnitude of expansion / contraction than does the value for direction B.

So the FC metric is in fact a measure of how much smaller or how much bigger a fixel is compared to the “norm ” (is the size of the template fixel pre-determined or is it based on the average fixel size?).

There is no “template fixel norm” here. If you look at the steps necessary for generation of FC data, there is no calculation of FC for the template fixel. FC instead has a sort of explicit normative value of 1.0, which is what you would obtain if you were to transform the template to itself, and therefore there would be no expansion / contraction contained within that transformation.

Your phrasing here also alludes to a “fixel size”, which would reasonably be interpreted as FD. FD does not in any way contribute to the calculation of FC in a standard pipeline, nor does it logically contribute in any way to the interpretation of what FC is. FC is a multiplicative parameter encoding the amount of expansion / contraction of fibres for those fibres in a given direction, which is invariant to the amount of fibres present in that direction. It is only because you are seeking to generate a voxel-aggregate measure that I am proposing that FD should have any influence here.

Essentially, I would want a metric that is the equivalent (or as close as possible) to summing absolute cross-sections values (in micrometer squared let’s say) of all fixels in a voxel, if that’s at all possible.

Okay, I think there’s multiple problems embedded in here that require unpacking.

  1. Terminology (thanks @dave :angry: :grin:). FC does not encode an absolute fibre bundle cross-section, but the change in cross-section that arises in the process of alignment to the template. So it’s both unitless and dimensionless. It says nothing about absolute “values” (which I’ll get to in a moment), only how much expansion or contraction occurs orthogonal to the fibre direction in the process of alignment of those fibres to template space.

  2. Metric selection. If FD gives you something that has dimensions L^3 (i.e. volume) and purports to be proportional to intra-cellular fibre volume, then that is actually closer to “absolute cross-section values” than is FC. Imagine a voxel that is 1 x 1 x 1 mm, containing a single coherently-oriented fibre bundle with a volume fraction of 0.5 (which will require further clarification in point 3). The intra-cellular volume in that voxel is 0.5mm3. The total intra-cellular cross-sectional area within that voxel is the volume divided by the length. If e.g. all fibres are perfectly oriented with the first image axis, the intersections of those fibres with the voxel cube are all 1mm. So the total intra-cellular cross-sectional area is 0.5mm2.

  3. Metric boundedness. If one were using a diffusion model that yielded volume fractions, which are bound between 0 and 1, then it is possible to multiply the volume fraction of any compartment with the volume of a voxel in mm3 to yield an absolute volume of that compartment within the voxel in mm3. However FD does not have this property. We can talk about the dimensions of FD being L^3, since it’s proportional to volume, but it doesn’t have units (one can think of it as being “in units of the response function” if you really want to understand SD properly, but I’ll maybe skip that here). As such, if you were to multiply FD values with voxel volumes in mm3, you will quite regularly observe apparent intra-cellular volumes that are larger than the volume of the image voxel within which it is supposedly being quantified. So we avoided attempting to translate these measures into “absolute values” or ascribe physical units.

Contemplating your question as written, I think the most faithful translation may actually be a voxel-wise FDC aggregate. The purpose of that measure is to centralise data into a common space, and for the volume of each voxel in that space (+ fixel specificity over and above that), quantify the “amount of information-carrying capacity” following such data centralisation. Values of such being larger or smaller in different subjects can be driven by differences in the local packing density of axons per unit volume in that subject prior to warping to the template, or by the absolute cross-section of that bundle being different in that subject, which necessitated expansion or contraction in the cross-sectional plane of that bundle in the process of warping to template space. So FDC gives you something for which “absolute” is a flawed but plausible descriptor. If dealing with voxel-wise data quantification, the distinction between “volume” and “cross-section” is not actually relevant (if you were instead looking for absolute cross-section of macroscopic bundles, that would be a different domain of quantification entirely, so I’m assuming that’s not your intent). But you won’t get to mm2 because of point 3.

Lastly, I did not really understand what the advantage of taking the determinant of the Jacobian of the non-linear warp field would be.

Well, the benefit I suppose would be that you would be immediately calculating a voxel-wise measure from the non-linear warp field; you wouldn’t even need the FOD information, let alone a fixel segmentation of such. But this is not ideal for white matter from a physical perspective Figure 2 of the FBA paper.

Could I just take the value of the subject2template.mif image file at a given voxel as an indication of the total deformation that was applied to the fixels contained in that voxel?

No, the content of that file encodes for each template voxel the location in the subject image from which data should be pulled; the Jacobian, and the determinant of such, are higher-level calculations thereof. Just as the warp2metric command can export FC per fixel, it can also export the Jacobian determinant per voxel.

Rob


1 Even with CSD with a hard non-negativity constraint, that constraint is only applied along a fixed number of directions; it’s possible for the FOD to creep just below zero in directions in between those along which the constraint is applied. This could break the equivalence, as that non-negativity would contribute to the l=0 term, but potentially be explicitly excluded from the voxel-wise FD sum derived from fixel segmentation, which uses a more dense sampling of the FOD (1281 directions) than does CSD (300 directions).

Hi @rsmith ,

I wanted to follow up on this thread. After following the fixel-based analysis on the main Mrtrix wiki, I have absolute effect sizes and p-values as a result of the statistical analysis. If I want to produce 3 maps of the effect size for the primary, secondary, and tertiary fixel, how would I do this? Would fixel2voxel -number 3 abs_effect_size.mif none abs_effect_size_voxel.mif automatically segment the 3 fixels with highest FD? Is the ordering of the peaks independent of the way fixel2voxel picks the fixels to map onto voxels (i.e. how it picks the 3 peaks)

Hi Ana,

In the current public code, the ordering of the fixels in each voxel is by peak amplitude rather than FD. It’s not currently possible to request fixel2voxel to place the fixels in a particular order before extracting them and writing to the voxel grid.

In the eventual 3.1.0 update, this will be changed such that fixels generated by dwi2fod will be ordered by FD (GitHub). If your experiment specifically needs that to be the case, you could try using the code on the dev branch; or, if you want to avoid all of the other changes on dev, you could cherry-pick just this change and re-compile the code.

You do however have to be careful in interpreting such maps in either case. Two fixels in adjacent voxels that are in the same orientation and are part of the same WM bundle are not guaranteed to be in the same image volume, as their respective positions in the order within their respective voxels may differ.

Rob

Hi Rob,

Thanks for the response! I’d like to avoid doing custom code since I’m on a compute cluster and running stuff aside from the preloaded modules tends to be a bit tedious, but thank you for offering a solution!

With respect to ordering by peak amplitude, will this still ensure that the primary/secondary/tertiary voxel maps will correspond to the first/second/third largest fixel in each voxel? If this is the case, I think that’s fine as well.

With the current code (3.0.x), where fixels are ordered by peak amplitude, the voxel maps will correspond to the first/second/third “largest” fixel in each voxel, where “largest” is defined here by tallest peak. So the voxel maps from fixel2voxel will not be in a “random” order, but they also won’t be in order of largest to smallest FD (which is what will happen in 3.1.0).

I’m asking this question, because it very much relates to this topic, and is somewhat a continuation of this thread.

I have done a group analysis looking at fd, log_fc and fdc using the standard pipelines and fixelcfestats. There are a number of areas which show up as significantly different between the groups, with all of the difference showing reduced values in the patient group.

I wanted to explore the data further, in particular to make sure there aren’t any outliers skewing the data in one way or another.

Based on this thread, I generated fd and fdc maps for each subject using fixel2voxel with the “sum” option, and log_fc values with the mean option, weighted by the mean FD map across all subjects.

I then multiplied these by the fwe maps, after those had been thresheld at 0.95 (i.e. p < 0.05) and binarized to obtain mean measures for each subject within these significant areas.

For the fdc and log_fc values, I get the expected result, namely that the mean value for patients is significantly lower than the mean value for controls.

However, for the fd values, I get the opposite result - namely, the value is higher in patients than in controls.

I realized that this may be due to the way MRtrix is doing the statistics and there may be a lower value in the primary component, but higher in the secondary component, and this is in fact the case. Patients have lower primary component, but higher secondary component.

But, now I am unsure how to interpret the result. Why would the fiber density components be going in opposite directions. Is there any biological/physiological way to interpret this? Any thoughts?

Thanks,

Mark

I think it would be better to use the scaled first volume of the FOD (l =0 term of the WM FOD spherical harmonic expansion) instead of the sum of FD, as was explained above. That would decrease errors due to fixel thresholding.

That result might still be valid though. If FD increases slightly but FC is decreased by a lot, I think that could result in lower FDC. I’m not sure though.

Thanks, Stephanie. By the “scaled” first volume, do you mean after applying mtnormalise (which is what we are using for the fixel data)?

Mark

Hi Mark,

I used the first volume of the warped (but not reoriented) FOD (in population template space). Mtnormalise had been applied in a step before registration and warping. I don’t think the reorientation is necessary since we are interested in the total fiber density and the l=0 term of the FOD represents “the orientational average of the FOD” (Calamante et al., 2015)

By scale what I mean is scaling by the “spherical harmonic basis scaling factor: 1.0/sqrt(4*pi)”, there is an example of this in the mrcalc documentation.

Hope this helps!

Hi Mark,

It is not impossible to have a scenario where group differences in FD and FC are of opposite sign; I can certainly think of one case where this has been observed and I think it’s quite likely to be reflective of reality. It’s even plausible, if the magnitude of the effect in one is considerably larger than the other, that you could get a statistically significant effect in FDC in such a scenario.

Your description of your initial analysis did not include which metrics were found to be statistically significant, nor the extent to which different regions were or were not significant for different metrics, and you also didn’t mention whether in those cases where FD and FC had opposite signs whether one was of larger magnitude than the other, so I can’t really infer the likelihood with which such observations may or may not be indicative of some issue. But I would at the very least not rule out the possibility that it may be a reasonable conclusion.

I would suggest avoiding thinking of these as “primary” and “secondary” components; FD and FC are ultimately just two different quantitative metrics. FDC being their product does not deem one in higher regard than the other.

FD being greater in the patient group could plausibly be an artifact of T2 shine-through, which, depending on the spatial distribution, may not be possible to correct based on a smoothly-varying bias field estimate. This is a persistent confound due to the way in which we suggest to intensity normalise, so it’s worth thinking about.

There’s also an alternative way to think about this kind of observation. As shown in the FBA manuscript, interpreting either FD or FC in isolation can produce misleading results; FDC is ultimately the metric that better encapsulates the total information-carrying capacity of the WM. In this manuscript we suggested performing statistical inference on FDC, with non-inferential quantification of the relative contributions of FD vs. FC toward the end result. In circumstances where FD and FC are of opposite sign, the parameter alpha (which is awaiting a rebranding), which encodes the fraction of the total effect in FDC that can be attributed to FD, will be less than 0.0 or greater than 1.0.


While using the SH l=0 term is useful in some circumstances, here the experiment specifically relates to post hoc interrogation of FBA results, so I would suggest primarily sticking with the fixel2voxel approach. It may however be interesting as a separate test to see if the group difference in the l=0 term echoes that of the sum of FD; if it does not, this suggests some difference in the outcomes of either fixel segmentation, or fixel correspondence, which opens a whole other can of worms.

Cheers
Rob