Compare two protocols for TEDANA

[roeysc]

Dear community members,

I have recently started working with Multi-Band Multi-Echo (MBME) data in a pipeline including fMRIprep and TEDANA (see details of scan protocols below). To test the pipeline, I collected fixation (“resting state”) data with two protocols.

I used fMRIprep 23.2.0, which outputs the single-echos after head motion correction + suceptibility distortion correctoin (no slice-timing correction due to short TR), and subsequently ran TEDANA (23.0.1) with the default tedpca=‘ica’ option.

As I have no experience in analyzing ME data, I would love to get your insight in order to choose the best protocol. My concerns are the following (I hope it’s okay they’re all in a single post):

1. I was worried that the high accelaration of MB4+GRAPPA2 might lead to bad motion artifacts. However, at least the carpet plots don’t show much difference compared with the MB4 protocol:

MB4 (before and after denoising)

MB4 + GRAPPA2 (before and after denoising)

I tend towards the MB4+ GRAPPA2 protocol, since its highest echo time is lower (42 ms < 63 ms). The Adaptive Mask for both protocols is 100% at 3, meaning all voxels within the brain mask had good signal for all echos. Is this reasonable?

2. It seems that using curve-fit has a strong impact on the denoising. Compare with the last image above:

MB4 + GRAPPA2 + curve-fit (after denoising)

Is this as surprising as I feel it is?

3. Both protocols led to ~27 ICs, with ~20% variance explained in the final accepted ICs. Is this much lower than expected for fixation data?

However, if I use tedpca=kundu, about 50% of the variance is retained:

These include many ICs with a broad spectrum and very diffuse whole-brain maps, so I lean against using it here. Do you agree?

To summarize, I tend to go with MB4 + GRAPPA2 using TEDANA with curve-fit.

I would love to hear it you think I’m missing something here.
In case it’s helpful, I uploaded a Drive folder with tedana_report.html file and the figures folder for each protocol and pipeline.

Many thanks for your time!
Best,
Roey

Scan protocol
Data collected on 3T Siemens MAGNETOM Skyra scanner with 32 channel headcoils, with fixation for 7:10 minutes.

• Protocol 1 used MB accelaration factor of 4 with 3 echo times (14, 39, 63 ms). TR=1000 ms (466 volumes).
• Protocol 2 used MB accelaration factor of 4 + GRAPPA factor of 2 with 3 echo times (12, 27, 42 ms). TR=760 ms (613 volumes).

[dowdlelt]

Lots of good questions here - My initial response is that your denoised data doesn’t look like much happened (until you added curve fit), but I favor faster sampling in general, even though MB4 x GRAPPA 2 is considered by many to be a pretty high accel rate. What is your voxel size? partial fourier? coverage (whole brain?). And even with short TR, I think it is valuable to apply that. It is likely to help the PCA-ICA step in fact!

1. I was worried that the high accelaration of MB4+GRAPPA2 might lead to bad motion artifacts. However, at least the carpet plots don’t show much difference compared with the MB4 protocol:

You will certainly be more sensitive to motion with this, I recommend using a FLEET reference scan (assuming you are using SIEMENS with CMRR or MGH sequence for this) - it is selectable in the same SIEMENS panel where you select the GRAPPA factor. It will be less sensitive to motion during the reference scan that occurs at the beginning. I choose FLEET over GRE because it is a) shorter and b) sounds like the rest of the scan so there is less potential for surprise movement when the tone of the scanner changes.

That said - the faster sampling due to MB4 helps. I’d be more worried about the g-factor related increase in thermal noise. Your data with GRAPPA 2 will have more noise towards the center of the brain especially, so keep that in mind if that is a region you are very keen on investigating.

2. It seems that using curve-fit has a strong impact on the denoising. Compare with the last image above:

This is very surprising to me, so much so that I would want to be certain that something else didn’t change. In my hands, particularly at 3T, curve fit has had a minimal impact, and the only differences that could really be observed existing in ventral temporal and orbitofrontal cortex where T2* values are extremely low (“dropout”).

So, yes, that is odd. I’ve certainly never seen it have a big impact on denoising, since it only changes S0 and T2* voxel-wise estimates. Not sure what to make of that.

3. Both protocols led to ~27 ICs, with ~20% variance explained in the final accepted ICs. Is this much lower than expected for fixation data?

I feel that given the number of TRs, I would want something more like 130 ~200 ICs (though I feel that kundu may have gone too far), but its been a minute since I’ve been directly engaged with this type of data. 27 BOLD comps might be ok, but 27 total feels much too low to me. Maybe there is a need to try a different selection metric (mdl versus aic) - @handwerkerd can probably chime in here with something more informative. There is also an upcoming pull request adding robust ica that may be ideal for this sort of problem, but that is still a work in progress.

Overall, I agree with your conclusion, I think your selected protocol is ok - but I do have several questions. Is the main analysis going to be resting state? Have you looked at seed based correlations in your original/denoised data? Does it look like how you would expect? What areas of the brain are you most interested in investigating? Is this the expected length of your final runs and will there just be one per subject?

I’ll try and take a closer look at your uploaded folders. It may also be valuable, if possible, to upload the data that was input into tedana, as well as its full output (denoised data, T2* map, etc). I know that will be a bit of bulk.

[roeysc]

Dear @dowdlelt (Logan!),

Thank you so much for this detailed reply. Now I really want to get to the bottom of this.

What is your voxel size? partial fourier? coverage (whole brain?)

MB4 protocol: 3.2 isotropic. Phase partial fourier of 6/8. Whole brain.
MB4_GRP2 protocol: 3.2 isotropic. Phase partial fourier of 7/8. Whole brain.

I recommend using a FLEET reference scan

We’ll try this!

27 total feels much too low to me. Maybe there is a need to try a different selection metric (mdl versus aic)

Thanks for the suggestion. Unfortunarely this gave rise to even less components (only 16, 4 accepted. Similar results for loglin and curvefit):

Actually the MB4 sequence (no GRAPPA) with aic and curvefit gave 39 PCs (20 accepted):

I’m surprised that the data with less volumes (466<613) gave more PCs.

Is the main analysis going to be resting state? Have you looked at seed based correlations in your original/denoised data? Does it look like how you would expect? What areas of the brain are you most interested in investigating? Is this the expected length of your final runs and will there just be one per subject?

The main analysis will be based on two scans:

• 2-3 scans of resting state (7:10 minutes each). This will be used for subject-specific parcellation of resting state networks.
• 1 scan of task (15:00 minutes). This is a continuous task (no block design), where we will be interested in activations around specific events, as well as changes in connectivity between different parts of the task.

So we are not focusing on specific brain regions in this analysis.
I haven’t yet looked at seed-based correlations, but that’s definitely down the line.

It may also be valuable, if possible, to upload the data that was input into tedana, as well as its full output (denoised data, T2* map, etc).

I’m happy to share this with you directly!

Many thanks for all you invaluable help,
Roey

[dowdlelt]

One explanation for the less volumes, more PCs could be that partial fourier + noise. By which I mean, more partial fourier (6/8ths versus 7/8ths) will lead to a little more smoothness in the data (a small, but noticeable effect) and GRAPPA 2 means the longer timeseries will have more noise in general. I’d still favor GRAPPA 2 - you’ll have less distortion, and less dropout in addition to the faster sampling rate. Also, seems like you may be more focused on cortex, broadly speaking, so higher noise levels in center of brain may not be largest concern.

regarding aic I had things backwards in my head - thats already the more liberal option, so it makes sense that mdl was “worse”. aic should return the most. Seems like this high accel, large voxel data may be misbehaving - not your fault, but still an issue.

I don’t think switching to FLEET will solve your problems, but it should improve TSNR in general and is a good choice.

Feel free to email me (check my profile to fill this in right, its my firstname.lastname@gmail.com, all lowercase) if you want to set up data sharing.

In general, as a next step, you could consider trying something like specifying the number of component directly, and we’ve had good luck with selecting ~1/3rd of the timeseries length as a guess. Maybe worth a shot, to see if denoising is better - I’m mostly concerned because its clear you will have some noise over such a long scan, but the effects you see are so small, tedana is doing very little. It may also require using something like RICA from the wonderful @e.urunuela to do manual component selection - annoying, but powerful.

I tend to recommend kundu - but as you’ve shown in your case it almost gives you too much to deal with, so for now I think no to that for now.

[roeysc]

Thanks so much for your help, @dowdlelt.
Could you please tell me how I might set the number of components directly? I could not find documentation of such an option on TEDANA’s manual.
Thanks!
Roey

[dowdlelt]

Of course, should have mentioned it directly - its the --tedpca option, detailed here, relevant bit:

–tedpca
Method with which to select components in TEDPCA. PCA decomposition with the mdl, kic and aic options is based on a Moving Average (stationary Gaussian) process and are ordered from most to least aggressive. ‘kundu’ or ‘kundu-stabilize’ are selection methods that were distributed with MEICA. Users may also provide a float from 0 to 1, in which case components will be selected based on the cumulative variance explained or an integer greater than 1 in which case the specificed number of components will be selected.

You could also specify variance explained with a float. Unfortunately, guidance kind of breaks down at this point, it is hard to recommend any specific settings - we often think about 1/3rd of the timepoints works well, but its all so dataset dependent!

I’ve seen the shared data, I’ll try and take a look.

[roeysc]

Thanks for pointing me to the --tedpca option! I missed the integer part there.
Unfortunately this didn’t help in terms of getting cleaner data.

I did upload the TEDANA reports for the additional data (in the same Drive folder). As you predicted, using FLEET doesn’t seem to make a strong difference (but may have helped with tSNR as you said). There it seems like there is a strong global signal that is not braing cleaned well by TEDANA.

MB4_GRP2 (curvefit+aic)
optimally combined + denoised:

I’m alittle worried about this second dataset. Perhaps these large flctuations are respiratory, and hard for TEDANA to handle? See these fMRIprep QC plots:

These large fluctuations might stem from the subject closing their eyes every once in a while, since overlaying the BOLD signal from some voxels in the eye shows clar correspondence between the GS (red) and the eye signal (blue), shown in another scan of that subject:

Thanks again for all your help!
Roey

[roeysc]

@dowdlelt

large voxel data may be misbehaving - not your fault, but still an issue.

Returning to Kundu et al., (2013), I’m surprised that a voxel size of 3.2 mm would be an issue, as they used voxels of 3.75 mm.

[dowdlelt]

I wasn’t very clear there - Less about the specific voxel size, but maybe the high level of acceleration. The PCA/ICA steps depend a lot on the noise levels of the data - so you have high SNR, but also a fair amount of noise from acceleration plus extra cross voxel correlations due to that same acceleration.

[roeysc]

I see. Thanks!
I will try scanning a full Fourier protocol to reduce the voxels correlations.

Do you think despiking the data after HMC and SDC but before tedana is also worth a try ?

[dowdlelt]

I suppose I wasn’t clear enough - I wasn’t intending that to be a recommendation, full fourier is going to take much more time (in the TE dimension), and I’m not sure you could spare it. I think 7/8ths is fine. It may be worth adding an additional echo onto the faster one - which would sacrifice a bit more time, but may give you better signal to work with. 3 echoes is necessary, but 4, 5, 6 etc can be helpful depending on the situation.

I’m still planning on taking a closer look at what you have shared - just wanted to make this note quickly.

[roeysc]

Hi NeuroStarsers, especially @dowdlelt!

I would greatly appreciate your insights regarding the latest data I’ve collected. Here I chose to reduce the MB factor from 4 to 3, collecting full Fourier data (at the expanse of a longer TR of 1093 ms).

The tedana QC plots (full reports are here on Drive):

You can see that a larger percentage of the variance is epxlained by accepted components, although the number of components is still surprisingly low (35 components for ~400 volumes).

The thing I’m worried about is the large amount of global signal artifacts, here in the optimally combined data:

GSR really gets rid of most of it:

And it’s a little hard to tell whether TE-dependent denoising makes a large difference:

I would love to hear what you think about it.
Also, since I don’t have much experience in tedana, I’m not sure if the following signal for the rejected components looks reasonable; I’m worried that the repeated patterns hereare indicative of some underlying problem I’m missing:

As always, your help is greatly appreciated!
Enjoy the weekend

[dowdlelt]

I’m still chewing on this - sorry for all the back and forth, and the bit of confusion, but the data is surprising me, in that its not behaving as I would expect, and nor is tedana.

Is there a reason to prefer fast sampling? I’m a fan personally, but I tend to think in 7T terms these days where things are a bit different. As I think I said before, maybe it is worth dragging out those echo times - with more echoes (to be clear, even if you need 7/8ths partial fourier - its not evil, technically).

I’m thinking in particular of the Lynch paper - https://www.sciencedirect.com/science/article/pii/S2211124720315291?via%3Dihub#sec4

From their methods:

Multi-echo, multi-band resting-state fMRI scans were collected using a T2∗-weighted echo-planar sequence covering the full brain (TR: 1355 ms; TE1: 13.40 ms, TE2: 31.11 ms, TE3: 48.82 ms, TE4: 66.53 ms, and TE5: 84.24 ms; FOV: 216 mm; flip angle: 68°; 2.4mm isotropic; 72 slices; AP phase encoding direction; in-plane acceleration factor: 2; and multi-band acceleration factor: 6) with 640 volumes acquired per scan for a total acquisition time of 14 min and 27 s.

I’m not sure that you would need to match exactly those parameters (perhaps stick with MB4, no need to go to 6!) - your relaxed voxel size may allow for less acceleration for example. I’m mostly recommending this because I don’t feel like the combined data looks quite like I would expect. I think those longer echoes could matter a lot. Some part of me recalls that the last echo should be approximately 2 x T2*, ie about ~66 or so? This could explain why your “slower” protocol had somewhat better behavior.

I think if you take your “faster” protocol, and add a few more echoes, out towards 70ms, you may see much better performance. @handwerkerd tends to be a bit smarter about this, so I’m curious what he thinks. I’m just getting the feeling from your data that the faster sampling is only hurtning you - and with 5 echoes (or so?), and a TR of 1~1.3 seconds you’d see very nice signal with that voxel size.

[roeysc]

Thank you again for these insights, @dowdlelt!

As for the sampling, we do need a TR or around 1s for the purpose of our experiment.
I wasn’t aware of this 2 x T2* rule-of-thumb. Thanks for bringing it to my attention.

I downloaded some online ME datasets to get a better feeling of what I might expect in our own data. This included the Multi-Echo Cambridge dataset with
TE: 12, 28, 44, 60 ms
TR = 2.47
GRAPPA 3 (I think not MB at all)

I ran fmriprep (this time with slice-timing correction due to the slow TR) and ran tedana with either all 4 echos or only the first 3 (which is more similar to our own data). I found the results quite perplexing. For example, for one subject (04570):

Tedana, curvefit, AIC, no GSR, 3 echos

with the first excluded component being (53% variance explained):

Tedana, curvefit, AIC, no GSR, 4 echos

with the 2nd accepted component (27% variance explained) being quite similar to the excluded one when we use only 3 echos:

Indeed, there are many more ICs compared with our data (~75) with only 139 volumes per scan.

I’m confused by this extreme difference in tedana’s behavior when taking 3 or 4 echos.

I will also add that for some subjects (like 11310), the carpet plot for the rejected components does have a similar repetitive texture to it, like we see in our data:

Thanks for staying with me here

[dowdlelt]

A quick note because I found the manuscript I was thinking of - Theoretical optimization of multi-echo fMRI data acquisition

If my quick read is right, then the target is ~1.5 x T2* of the tissue, so ~60 is perhaps even on the low end. In other words, the more echoes, the better, and the longer the better (up to a point).

I think the analysis you’ve done there in instructive for this - that last echo is adding a lot of useful information about the BOLD contrast, and therefor helping the fit. I’ll try and take a closer look shortly.

[roeysc]

Thanks, Logan!
I also uploaded the reports (and figures folders) as a fourth dataset here.

Thank you for the reference! I will look into it.

[roeysc]

with the 2nd accepted component (27% variance explained) being quite similar to the excluded one when we use only 3 echos:

So you think that this large component, with this specific timecourse is in fact BOLD related and should be rtained? My intuition was that this is indeed a bad component that should be discarded, but perhaps I’m wrong. Looking at the denoised data, the 3-echo version actually seems cleaner to me (number of echos appears in the top left label):

and similarly for another subject:

[dowdlelt]

I read that completely backwards. That is a bit goofy and it is hard to pin down any specific behavior - but I agree, noise. I will say that if the noise sources like that are inaccurately retained, then polynomials in the model (or detrending/filtering) would remove the bulk of that variance - but the benefit of higher BOLD contrast due to the 4th echo being retained would remain.

These sorts of issues are also a reason tedana devs are working on some details that would allow enable more information to be passed into the model. For example, one of those artifact looks like motion to me (perhaps), and if the component had a high correlation with motion parameters, it could then be rejected in a more principled way. Alas, nothing ready quite yet.

I still think that a MB4/GRAPPA 2 4~5 echo version of your data would be closer to ideal, but I need to look over the analyses you’ve done thus far to little closer to make sure I"m not saying something too early in the morning again.

[roeysc]

For example, one of those artifact looks like motion to me (perhaps), and if the component had a high correlation with motion parameters, it could then be rejected in a more principled way. Alas, nothing ready quite yet.

Thanks, @dowdlelt. I did see a post from @handwerkerd regarding the potential integration of both approaches.

I will look into adding another echo or two to our data (hopefully without GRAPPA, as I’m worried about subject motion during the reference scan, but we’ll see about that).

[dowdlelt]

I’m less worried about motion with FLEET - its quick and somewhat robust against motion - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4580494/ but of course would still be worse than no GRAPPA in a motion case (probably - you are already using multiband which also requires a slow reference scan - the SBRef)

GRAPPA=2 would mean nearly 2x as many echoes, which is a lot of signal to leave out, imo. But I’m used to pretty low motion subjects so maybe my advice doesn’t apply here.