adaptive mask uses different thresholds in t2smap.py and tedana.py, which causes an inconsistency problem in the OFC areas

[zswang-brainimaging]

Summary

In the older version of t2samp, "mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=3)" causes signal loss problem in the prefrontal areas, especially in the OFC, which is very obvious and problematic. In the latest version of t2smap.py, "mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=3)" has been changed to "mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=1)", meaning that 3 good echo data are no longer required. Yes, using threshold =1, the issue of signal loss in the prefrontal areas are resolved, where optimally combined image desc-optcom_bold.nii looks much better now. But in the latest version of tedana.py, "mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=3)" still uses threshold =3 (3 good echo data are required), which causes the same signal loss problem in the prefrontal areas, especially in the OFC. Try to modify threshold =1, but running tedana.py failed, with some errors, one of which is: "IndexError: index 0 is out of bounds for axis 1 with size 0".

Additional Detail

if modifying " mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=3)" in tedana.py into
" mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=1)", the running procedure will fail, but if modifying " mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=3)" to " mask, masksum = utils.make_adaptive_mask(catd, mask=mask, getsum=True, threshold=2)" , the running procedure can be completed. The result image dn_ts_OC or ts_OC looks better than ones generated when using threshold =3, but still not as good as the image "desc-optcom_bold.nii.gz" generated using t2smap.py when settting threshold=1 in the adaptive mask.

Next Steps

[emdupre]

Hi @zswang-brainimaging --

Thanks for opening this ! Indeed, the recent release explicitly changes the behavior in the t2smap as compared to the tedana workflow. This was a conscious decision based on discussion in #617 that the two workflows have two slightly different goals.

In t2smap, we're just generating the optimal combination with no intention of downstream (multi-echo informed) denoising. In this case, we use any echo with reasonable signal (as determined by threshold=1) to construct our adaptive mask. By contrast, in tedana the goal is explicitly to denoise, so we need to have signal in multiple echos available (threshold=3) to calculate our echo-related denoising metrics for each estimated component. As you've seen you can reduce the number of echos necessary (with e.g. threshold=2), but this reduces the accuracy of the echo-related denoising metrics.
EDIT: Note that you can't reduce tedana to threshold=1 as we cannot calculate the echo-related denoising metrics with only a single echo !

All that to say, the behavior you're seeing is as expected, and the inconsistencies between the workflows come from their divergent goals !

Does this help to explain the "inconsistency" ? Let us know.

[zswang-brainimaging]

Hi @emdupre

   Thank you so much for detailed reply!
    So now understood better your strategy/algorithm to deal with denoising the ME data.  But at the meantime, this kind of strategy/algorithm (requires at least 3 good echo data) will reduce the advantages of using multil-echo data although  noise can be reduced.  Based on our observations,  prefrontal (e.g. OFC) signals are strongest in echo 1,  then disappear very quickly after echo 1 (almost nothing even in echo 2).  So how to deal with this problem if 3 good echo data are required ?  One of big reasons to use ME fMRI is to get OFC signals even if they are weak or contain noise, but now they disappear completely  in tedana algorithm.  So we may need to find a better way to deal with this.   For those inherently weak signal regions such as OFC, can we use a different way to deal with?  For example, can we use the optimally combined data to replace echo2  and echo3 data within these regions only,  to avoid echo loss before running the denosing procedures ?  

 And how "good" echo and "bad" echo are judged ?  Why we need to do this kind of judgment before running ICA?  ICA is central here to apply for removal of noise, right ?  Usually in an ICA-based denoising algorithm,  PCA is applied to all the raw data  (they may be pre-masked to reduce memory requirement using a skull-stripping mask for example) to reduce the dimension.  Then ICA is applied to the PCA-compressed data.  Then good components and bad components are identified.  After removing bad components,  images are reconstructed from the good components.  If we have 5 echo data, for example, why we do not use all the 5 echo data for PCA and ICA?   After ICA, we get 5 noise-removed echo data, then we can apply t2smap to the 5 echo data to optimally combine them, right  ? 

 My understanding is that before doing optimal combinations of the multiple echo data, S0-related noise/artifacts should be removed.  Applying ICA  to here is already a good idea:   Remove S0-related noise from all the echo data, then run t2smap to combine all the noise-removed-by-ICA  echo data.  This should be straightforward.  But in the tedana package, many complicated algorithms are applied before running ICA, even good echo and bad echo, and good components and bad components are judged before running ICA.   Hard to understand these procedures.  

[zswang-brainimaging]

Thank you so much for detailed reply!
So now understood better your strategy/algorithm to deal with denoising the ME data.
But at the meantime, this kind of strategy/algorithm (requires at least 3 good echo data) will reduce the advantages of using
multil-echo data although noise can be reduced. Based on our observations, prefrontal (e.g. OFC) signals are strongest in echo
1, then disappear very quickly after echo 1 (almost nothing even in echo 2). So how to deal with this problem if 3 good echo
data are required ? One of big reasons to use ME fMRI is to get OFC signals even if they are weak or contain noise, but now
they disappear completely in tedana algorithm. So we may need to find a better way to deal with this. For those inherently
weak signal regions such as OFC, can we use a different way to deal with? For example, can we use the optimally combined data
to replace echo2 and echo3 data within these regions only, to avoid echo loss before running the denosing procedures ?

 And how "good" echo and "bad" echo are judged ?  Why we need to do this kind of judgment before running ICA?  ICA is 

central here to apply for removal of noise, right ? Usually in an ICA-based denoising algorithm, PCA is applied to all the raw
data (they may be pre-masked to reduce memory requirement using a skull-stripping mask for example) to reduce the
dimension. Then ICA is applied to the PCA-compressed data. Then good components and bad components are identified.
After removing bad components, images are reconstructed from the good components. If we have 5 echo data, for example,
why we do not use all the 5 echo data for PCA and ICA? After ICA, we get 5 noise-removed echo data, then we can apply
t2smap to the 5 echo data to optimally combine them, right ?

 My understanding is that before doing optimal combinations of the multiple echo data, S0-related noise/artifacts should be

removed. Applying ICA to here is already a good idea: Remove S0-related noise from all the echo data, then run t2smap to
combine all the noise-removed-by-ICA echo data. This should be straightforward. But in the tedana package, many
complicated algorithms are applied before running ICA, even good echo and bad echo, and good components and bad
components are judged before running ICA. Hard to understand these procedures.

[emdupre]

I'll try to unpack this a bit here to help in the discussion !

And how "good" echo and "bad" echo are judged ? Why we need to do this kind of judgment before running ICA? ICA is central here to apply for removal of noise, right ? Usually in an ICA-based denoising algorithm, PCA is applied to all the raw
data (they may be pre-masked to reduce memory requirement using a skull-stripping mask for example) to reduce the
dimension. Then ICA is applied to the PCA-compressed data.

This is what's happening in tedana. It might be most helpful to look at this figure from the documentation. It's a little out of date (and your opening this issue is a great reminder that we should update it, thank you !), but we can follow the general idea.

First, we take the input data and fit our decay model to calculate the optimal combination. This is done in both the t2smap and tedana workflows, but with different thresholds as you've identified. t2smap stops here, but tedana keeps going.

The optimally combined data is then passed to PCA. In the original Kundu methods (also known as ME-ICA, see Kundu et al 2017, NIMG, for example), this "TEDPCA" included a decision tree that integrated echo-specific information to judge whether PCA-derived components were more likely signal or noise. Although we retain the option to call this PCA-specific decision tree (see the docstring for tedana) by default we do something more standard: an automatic variance selection criteria based on minimum-description length (adapted from GIFT). This is just to remove low variance components and is commonly done for most neuroimaging ICAs (see MELODIC, for another example).

The data after PCA are then passed to ICA and decomposed into independent components. The independent components are then evaluated as BOLD or non-BOLD based on several criteria, including projection back onto the original echos. This is why we need a minimum number of good echos available -- otherwise it's unclear how to evaluate the components.

Then good components and bad components are identified. After removing bad components, images are reconstructed from the good components.

This is exactly what's happening, but the key difference from getting five denoised time courses for five echo data is that the optimal combination is what went into our PCA + ICA decomposition. So we get back the reconstructed optimal combination. If we wanted the individual echos, we would need to reverse the optimal combination procedure, but this still wouldn't solve the issue you've identified, since we needed several good echos to get reasonable denoising (at least based on our current criteria !).

There are multiple paths forward, here. The one we're thinking most about is to develop a more flexible decision tree that will allow "less rigorous" denoising (compared to our current pipeline, but maybe not to your particular use case) and therefore allow you to maintain more signal. The challenge there is that we (1) need to build that flexibility into the decision tree (something @tsalo and @handwerkerd are working on !) and (2) need to have a means to allow users to transparently report exactly what decision criteria they chose to operate under.

Hopefully that clarifies a bit what is happening, even if it doesn't immediately address your use case.

[zswang-brainimaging]

@emdupre :

     Again, thank you so much for taking your invaluable time to answer my questions patiently and tirelessly .  

I have learnt a lot from your detailed explanations. Really Appreciate ! I still have some questions that I want to
discuss with you. I try to simplify them this time:

1. As you said, PCA+ ICA were applied to the optimally combined data (data_oc), not all the original echo data instead.
   This idea seems different from the paper "Differentiating BOLD and non-BOLD signals in fMRI time series using multi-
   echo EPI" (Prantik Kundu et al, NeuroImage 2012, 1759-1770). In that paper, Equation 7 defined that ME-ICA was applied
   to all the echo data by concatenating all the Ne echo date along space. After ME-ICA, components were identified from this
   type of outcome. Not sure why tedana was not implemented this way.

2. In tedana.py, Can we still use threshold=1 to generate adaptive mask and optimally combined data (data_oc) as we do in
   t2smap.py? Since we will run PCA+ICA based on data_oc, this data need to contain as much information as possible,
   right? Then, for those brain regions (like OFC) with weak signals (they should be considered weak, not necessarily bad,
   right?), if their corresponding masksum=1 (meaning that the signals exist only in echo 1 most likely), then can we just
   use one echo data (like echo 1) for the component evaluation purpose?

3. Again, based on our observations, OFC and other prefrontal signals only exist in Echo 1. They are weak, but not necessarily
   bad. At least for those regions. the evaluations of the components cannot be much relied on the number of good echo
   data. The data in these regions may need to be treated differently.

[emdupre]

1. As you said, PCA+ ICA were applied to the optimally combined data (data_oc), not all the original echo data instead.
   This idea seems different from the paper "Differentiating BOLD and non-BOLD signals in fMRI time series using multi-
   echo EPI" (Prantik Kundu et al, NeuroImage 2012, 1759-1770). In that paper, Equation 7 defined that ME-ICA was applied
   to all the echo data by concatenating all the Ne echo date along space. After ME-ICA, components were identified from this
   type of outcome. Not sure why tedana was not implemented this way.

This was not consistent across ME-ICA implementations, and in later papers (see eg Kundu et al., 2017, NIMG) the decomposition on optimally combined data was preferred, see eg the section on Making ME-fMRI practical in that paper which states:

tedana.py first does ME-PCA to isolate signals for ICA, by default on optimally combined data. It then finds ICA components using FastICA, implemented with Numpy and Scipy routines.

If you look in the original ME-ICA codebase, you can see that this was controlled with an option called --sourceTEs, where -1 was the optimal combination (the default) and 0 was the z-concatenated echos (the option you're referring to). tedana in fact used to include this option of z-concatenated echos, but this was removed based on discussion here. This is one of the tedana project differences mentioned here.

2. In tedana.py, Can we still use threshold=1 to generate adaptive mask and optimally combined data (data_oc) as we do in
   t2smap.py?

The workflow itself will fail, as it requires signal in at least two echos to complete the current decision tree.

Then, for those brain regions (like OFC) with weak signals (they should be considered weak, not necessarily bad,
right?),

Absolutely ! I don't mean to say that they're bad, but I also can't say that they're good -- I just don't have that kind of information available, at least based on our current metrics.

if their corresponding masksum=1 (meaning that the signals exist only in echo 1 most likely), then can we just
use one echo data (like echo 1) for the component evaluation purpose?

Our current decision tree (see here, for example) needs more than one component. You could just use that data, regardless of any decision tree criteria, but then it's hard to know what that data would mean. It would be denoised in those voxels with sufficient signal, but it would not be denoised in those regions with insufficient signal, like OFC. Again, though, this is something that could be addressed in a more flexible decision tree, where you could compare echo-specific and more general (e.g. AROMA-like) component evaluation.

3. Again, based on our observations, OFC and other prefrontal signals only exist in Echo 1. They are weak, but not necessarily
   bad. At least for those regions. the evaluations of the components cannot be much relied on the number of good echo
   data. The data in these regions may need to be treated differently.

I agree with you that I don't want to pass judgment on those OFC signals ! And I'd also agree that they likely need to be treated differently. I think the immediate question is : what would that "different treatment" really look like, and how could we achieve it in a principled way ?

---

Personally, I'm pinning a lot of my hopes on updating the decision tree, here, or maybe the decomposition itself. Again, this doesn't help your immediate issues, but for that I see at least a few paths forward:

1. Use the optimally combined data and denoise using principled component evaluation criteria that are not echo-specific, such as those implemented by ICA-AROMA. You can create optimally combined data and run AROMA directly in the latest fMRIPrep version.
2. Run the original ME-ICA implementation with --sourceTEs set to 0. This would best capture the behavior documented in the 2012 paper.

Of course, we'd love your continued feedback as we update the decision tree ! But this is an ongoing process, and I know that your data processing needs likely have some time pressure.

[emdupre]

I just wanted to explicitly link this discussion to a previous one (#600 (comment)), since @tsalo also provides a nice walk-through there of the current approach in the tedana workflow !

[zswang-brainimaging]

@emdupre :

   Again, thank you so much for your last detailed explanations.  Very helpful for understanding the algorithm here.  

After we run more tests, it has become clear that t2smap now is working well to keep the OFC signals when using
threshold=1 in the adaptive mask function. But tedana with ME-ICA still has the same problem of OFC signal loss
when setting threshold=3 in the adaptive mask.

So can the algorithm be slightly modified to use threshold =1 for those prefrontal areas with weak signals ? (meaning
that components are checked for quality just using echo 1 data since other echo have do signals in those regions ?).

Thank you!

[zswang-brainimaging]

"since other echo have do signals in those regions ?" : a typo, actually, I meant "since other echo have no signals in those regions ?"

[zswang-brainimaging]

Dear all: For your following note: “ If you look in the original ME-ICA codebase, you can see that this was controlled with an option called --sourceTEs, where -1 was the optimal combination (the default) and 0 was the z-concatenated echos (the option you're referring to). tedana in fact used to include this option of z-concatenated echos, but this was removed based on discussion here. This is one of the tedana project differences mentioned here.” Could you please email me this original code with sourceTEs option ? Thank you so much! David

…

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.