Using Radial Search for Singular ROI, as opposed to between 2 ROIs

Hello,

I have a 10M streamline whole-brain tractogram (GMWM-seeded and copped), a set of functional ROIs, and 30 reconstructed bundles from Dipy Recobundles, all in subject space. My goals are to:

  1. Find streamlines that begin/terminate at each fROI, to save as individual .tck or .trk files. I seeded in GMWM since I figured this may cut out some redundant streamlines that don’t reach GM.
  2. Intersect these functional ROI streamline files with Recobundles tracts. This would answer a question such as “What portion of the left arcuate is associated with this reading-activated functional area”?
  3. Make tract profiles plotting FA along these “sub-bundles” and run a group comparison of these values.

Questions on Part 1

I know that for tck2connectome, a radial search is executed to assign streamlines to nodes. Could this be applied to find all streamlines corresponding to a single node, as opposed to between two connectome nodes?

I ask because for something like tckedit, the -include mask seems to be more explicit, meaning no streamlines would likely get assigned to a GM mask because the streamlines are seeded and cropped at GMWM. Maybe it would make sense to create a sphere (x,y,z,radius) surrounding the peak of each functional area, but I would not know what radius would be appropriate, and also fear of overlap due to the regions being close to each other.

Questions on Part 2

Does such functionality exist in MrTrix? I think I have an idea of how to implement this in Dipy, but it is a bit convoluted. Do you have any recommendations for overlapping to streamline files? Part 3 I think can be done on Dipy relatively trivially.

SIFT2 Interpretability for this Application

I can additionally run SIFT2 on the tractograms, but am having trouble interpreting the streamline weights in this context. That is, when I find the streamlines associated with an fROI and later its intersection with segmented tracts, I am oriented to a binary scale; that is, is this streamline passing through the ROI / intersecting with a segmented bundle or not? I also worry about compatibility with DIPY if I do have to use it for tract metrics. I could alternatively remake tractograms with 100M and SIFT to 10M, which would resolve that issue.

Sorry for the loaded questions, but any suggestions would be appreciated!

Best,
Steven

1 Like

Hi Steven,

… and 30 reconstructed bundles from Dipy Recobundles …

If my interpretation is correct, the primary outputs from this process are a set of track files? This is important as the nature of those data dictate what operations are necessary to interact with them.

I seeded in GMWM since I figured this may cut out some redundant streamlines that don’t reach GM.

I’m interpreting from this that you are using ACT, since you stated that your streamlines were cropped at the GM-WM interface. Sometimes people don’t use ACT but just seed based on a GM mask; but that wouldn’t enable cropping at the interface (unless you did that explicitly afterwards using tckedit). If my interpretation is correct then some of my suggestions below may not be accurate.

If ACT is used, then you don’t have “redundant streamlines that don’t reach GM”. Streamlines are forbidden from terminating inside the WM. So seeding at the GM-WM interface does not influence the prevalence of such streamlines.

Intersect these functional ROI streamline files with Recobundles tracts. This would answer a question such as “What portion of the left arcuate is associated with this reading-activated functional area”?

What is the reasoning behind not simply looking at the intersection between the Recobundles tracks and the functional ROIs? It’s not immediately clear to me why a whole-brain tractogram needs to be used as an intermediary between those two data.

Nevertheless:

Do you have any recommendations for overlapping two streamline files?

There isn’t really such a thing as an “intersection” between two streamlines (or two sets of streamlines), since each streamline is infinitessimally thin. The first possibility I would be trying is generating a TDI for each, and then performing some quantification based on comparing those maps.

I know that for tck2connectome , a radial search is executed to assign streamlines to nodes. Could this be applied to find all streamlines corresponding to a single node, as opposed to between two connectome nodes?

I think you’re trying to manipulate the process to get the data you want, whereas you can get the data you want from the existing process.

Once you have a connectome matrix (and the associate streamline-node assignments) via tck2connectome, you can use connectome2tck to extract specific connections of interest. This command is quite flexible; by default it outputs one track file per edge, but that’s not the only possibility. If you specify -files per_node, what you’ll get is one track file per node, where each contains all streamlines that are assigned to an edge for which one of the streamline endpoints is assigned to that node. If you additionally specify -nodes 1 (assuming your blob of interest contains index 1), only a track file for that one node will be produced.

There’s a trick to making this slightly more robust. Say, for instance, that all you have are two separate ROIs, and you’re not interested in sub-division of anything else in the brain. What you want to do is synthesize a parcellation image, where voxels in ROI 1 have the value 1, voxels in ROI 2 have the value 2, and all of the grey matter that is not within either of those two ROIs has the value 3. That way, streamlines that terminate near the functional ROI, but are better characterised as intersecting whatever is adjacent to that ROI, will get assigned to dummy parcel 3 instead of the ROI. This is just mitigation of the imperfections of that radial search mechanism (which I intended to supersede years ago…).

What data manipulations do or do not work here can also depend on the nature of your data; e.g. whether the functional ROIs are big spherical blobs or constrained to the cortical ribbon, and what the distribution of streamlines terminations looks like. The radial search was made the default as it’s a reasonable compromise across a wide range of usage scenarios (the developer doesn’t know a priori the attributes of any data a user may provide it with), but it’s far from perfect.

I ask because for something like tckedit , the -include mask seems to be more explicit, meaning no streamlines would likely get assigned to a GM mask because the streamlines are seeded and cropped at GMWM.

Maybe more “basic” rather than more “explicit”; they’re both explicit, it’s just that one is predicated on an overlap whereas the other is predicated on a maximal distance. Nevertheless, that’s precisely why the radial search was necessary for tck2connectome. But ideally tckedit would have a solution to the same problem.

That is, when I find the streamlines associated with an fROI and later its intersection with segmented tracts, I am oriented to a binary scale; that is, is this streamline passing through the ROI / intersecting with a segmented bundle or not?

Yes, in this framework the attribution of each individual streamline to the pathway of interest is binary; but the streamline weights represent independent information to such. You can still use that information to modulate how much each streamline contributes to whatever metric you are quantifying.

I also worry about compatibility with DIPY if I do have to use it for tract metrics. I could alternatively remake tractograms with 100M and SIFT to 10M, which would resolve that issue.

If resolving the biases that SIFT(2) address are important for your experiment, but you can’t utilise per-streamline weights properly, then yes, you could do that instead. Though you can also look at the comprehensive set of operations within the pipeline and ask whether or not such modulation is actually going to have any effect on your eventual output quantification.

Also, if you instead use dynamic seeding, it’s quite common for reducing streamline count by a factor of 10 to not even be possible for the algorithm, since the magnitude of the bias in the initial tractogram is reduced; so you could do say 50M → 10M instead.

Rob

1 Like

Yup! 30 .trk files

Yes

Good to know, thanks!

There are analyses I am trying to replicate from another paper: https://www.nature.com/articles/s41467-019-11424-1 . One such question they address is “what fraction of fibers eminating from an fROI are from a given tract?” But I agree that for just finding the tract/fROI intersection, the whole-brain tractogram intermediate step may not be necessary.

Perfect

Rest of the answer is a bit long to requote here, but I understand what you’re saying and it certainly gives me plenty to think about and some good next steps to investigate. As always I really appreciate your thorough and thoughtful responses, thanks!

Best,
Steven

Okay, that’s a really useful definition of the task / hypothesis.
Though I would note that compared to your original stated goal:

This would answer a question such as “What portion of the left arcuate is associated with this reading-activated functional area”?

, this is actually the inverse quantity; i.e. “what portion of this reading-activated functional area is associated with the left arcuate”.


For this specific question, being able to use SIFT(2) is theoretically beneficial, as the relative cross-sectional contributions of different streamlines trajectories weigh in on that definition of “fraction”. I.e. If half of the streamlines intersecting the fROI are ascribed to the left arcuate and half are not, but those streamlines in the former group are assigned greater weights by SIFT2 than those in the latter group, then one would ideally want the conclusion of “fraction of “connectivity” involving this fROI belonging to the arcuate” to be greater than 50%.

You still however have the challenge of, for every streamline trajectory deemed to intersect the functional ROI, making a binary decision on whether the streamline does or does not “belong” to one of a set of pre-defined WM pathways (this decision is independent of its SIFT2 weight). There’s no single technical solution for that, and there’s all sorts of confounds to consider; e.g. can one such streamline be assigned to more than one of your bundles of interest? What about if the trajectory corresponds to some known WM pathway that isn’t included in the set of 30 segmented bundles, but nevertheless has “some similarity” with one of those bundles?

As a first pass, making use of existing tools, you could probably do something like:

  • Produce a binary mask for each WM bundle of interest from the output of RecoBundles, using tckmap | mrthreshold

  • Select those streamlines for which one endpoint is ascribed to the fROI.

  • For each of those streamlines:

    • Generate a voxel visitation map using tckmap

    • For each of the 30 bundles, quantify the fraction of the voxels traversed by the streamline that are included in the mask for that bundle (can do with mrcalc | mrstats trickery)

    • Select the bundle for which this overlap is maximal

    • If that overlap exceeds some threshold, assign the streamline to that bundle; if not, assign the streamline to the “everything else” group

    (note how this structure prevents assignment of one streamline to more than one bundle)

  • Take the ratio between the sum of weights of those streamlines attributed both to the fROI and the bundle of interest, and the sum of weights of all streamlines attributed to the fROI.

It ignores orientation information, and will perform poorly in the hypothetical situation where the RecoBundles outputs are spatially conservative estimates and your streamlines traverse adjacent to a bundle of interest (and hence don’t “overlap”), but it’s a reasonable first pass.

Have fun,
Rob

1 Like