Fmriprep: Susceptibility distortion correction (SDC) using Spin Echo Fieldmaps

[sameera2004]

Hi,

I pre-processed resting state data using fMRIprep 20.1.1. This data set contain 20 subjects, 2 sessions for each subject, and 4 runs for each session. I used SpinEcho field maps to do the susceptibility distortion correction (SDC). The SDC with SpinEcho’s worked for majority of the data. However, it worsen the SDC’s in certain subjects. Below is the command that I used to run fMRIprep;

singularity run --cleanenv --bind /mnt/fMRIprep /mnt/tools/fMRIprep_sing/fmriprep-20.1.1.simg /mnt/fMRIprep/scR21 /mnt/fMRIprep/scR21_outputs -w /mnt/fMRIprep/scR21_outputs/derivatives/scratch --fs-license-file $PWD/license.txt --write-graph --cifti-output 91k participant --participant-label 50015 50016 50020 50024 50027 50033 50042 50137 50144 --ignore slicetiming --output-spaces T1w MNI152NLin2009cAsym:res-2 fsaverage fsLR
I am wondering if anybody else has experienced similar issues. Please let me know.

Thanks in advance
Best Regards
Sameera

[oesteban]

Hi @sameera2004, I assume that the above is before SDC and the below shows the image after SDC? If that is the case, you might have a wrong PhaseEncodingDirection parameter in your sidecar JSON. Seems the correction was executed in the opposed direction.

[jxvansne]

Hi @oesteban, Sameera works in my lab. I believe the two images she posted were an AP and a PA phase encode direction run (we alternate directions from one run to the next, HCP style) after SDC for a single subject.

When I first saw this back in June, I also thought the PE direction had been somehow specified wrong. We went back and flipped the specified direction in the json files (i.e., we changed j- to j and j to j-), and this produced even worse results than you see here. I could post the before-and-after SDC results for both directions if you need to see it, but it might take a few days to dig up / redo since we did this back and June, and I'd want to be certain I was posting the correct results.

Ultimately this issue resulted in us using the HCP pipeline rather than fMRIprep for a forthcoming manuscript. I'm in the process of trying to make a final pipeline choice for multiband fMRI in my lab... a side-by-side review of HCP and fMRIprep results showed that on the whole, HCP was doing a better job across more subjects: even though HCP completely failed one subject in a dataset of ~20 subjects (a subject that fMRIprep had no issue with at all), there would overall been more data loss due to unacceptable SDC results with fMRIprep, which was entirely due to SDC results like what you see in the above posts from Sameera.

Let us know what we can post here (.json files, the fMRIprep .html outputs showing SDC, anything else) to help troubleshoot this issue.

Thanks,
-Jared

[effigies]

Would you be able to share a problematic subject?

[oesteban]

a side-by-side review of HCP and fMRIprep results showed that on the whole, HCP was doing a better job across more subjects

I'm curious about how you run this QC step. It might be interesting for fMRIPrep to provide some visualization based on your experience. Could you share a bit more about how you did this?

@effigies
Would you be able to share a problematic subject?

This would be really important. Is this possible for you?

[jxvansne]

a side-by-side review of HCP and fMRIprep results showed that on the whole, HCP was doing a better job across more subjects

I'm curious about how you run this QC step. It might be interesting for fMRIPrep to provide some visualization based on your experience. Could you share a bit more about how you did this?

We have a MATLAB OOP post-processing pipeline that can now take in either HCP or fMRIprep pre-processed results and do additional analysis. The first major output of that pipeline are the QC plots that Sameera cropped the bottom right quarter of above. So we use those three orthogonal views of the mean EPI image (average over time for a single run) to quickly evaluate SDC / coregistration / normalization. So for 8 runs of data for 20 subjects I want through and did a side-by-side comparison of HCP and fMRIprep. I'm mostly looking at the mid-saggital slices for the length and shape of the corpus callosum, as well as the shape and size of the OFC dropout region and anterior mPFC.

@effigies
Would you be able to share a problematic subject?

This would be really important. Is this possible for you?

We're working on this, which is why we hadn't responded yet. I think we should be able to post defaced data in the coming weeks. We'll be looking to provide T1w, T2w, AP/PA spin echo fieldmaps, B0 fieldmaps, and the EPI data itself (defaced). Am I missing anything else?

[effigies]

We're working on this, which is why we hadn't responded yet. I think we should be able to post defaced data in the coming weeks. We'll be looking to provide T1w, T2w, AP/PA spin echo fieldmaps, B0 fieldmaps, and the EPI data itself (defaced). Am I missing anything else?

Any pertinent JSON files. The simplest thing to ensure consistent behavior with what you've seen is the entire dataset (including all dataset-level JSON and TSV files) and just drop all but one or two subjects.

[sameera2004]

Hi @effigies and @oesteban ,

I am de-identifying problematic subject data. Yes we will include JSON and TSV files too, after removing identifiable personal health information from them. Please let me know where do you want me to upload those data.
Thanks
Sameera

[effigies]

The simplest way is to upload to OpenNeuro. You can keep the dataset private and share it with me. My login is my GitHub username @ gmail.com.

[sameera2004]

Hi @effigies,

I uploaded to OpenNeuro and shared the dataset with you. Please let me know if you have any issues.
Thanks
Sameera

[effigies]

Hi @sameera2004, I didn't find an invitation. Could you share the link to your dataset?

[sameera2004]

Hi @effigies,

Here's the link to data set https://openneuro.org/datasets/ds003130. I just sent another invitation.
Thanks
Sameera

[effigies]

Thanks. I can see the dataset and will have a look.

[sameera2004]

Hi @effigies and @oesteban,

Any updates regarding this SDC issue?

Thanks
Sameera

[effigies]

Sorry, these last couple weeks have been pretty full. I'll be looking into this this week.

[effigies]

Hi @sameera2004, we've had a look at this now. We want to verify the intended phase-encoding directions of each of these runs:

import nibabel as nb
import bids
layout = bids.BIDSLayout('/data/bids/ds003130-download/')
imgs = layout.get(extension='.nii.gz', datatype=['fmap', 'func'])
dirs = {"A": {"j": "PA", "j-": "AP"}}  # Truncated for simplicity
for bfile in imgs:
     img = bfile.get_image() 
     ornt = ''.join(nb.aff2axcodes(img.affine)) 
     md = bfile.get_metadata() 
     peaxis = md['PhaseEncodingDirection'] 
     pedir = dirs[ornt[1]][peaxis] 
     print(f"{bfile.filename:<46} {ornt} {peaxis:<3} {pedir}")

Output:

sub-50015_ses-1_dir-AP_run-1_epi.nii.gz        LAS j-  AP
sub-50015_ses-1_dir-PA_run-1_epi.nii.gz        LAS j   PA
sub-50015_ses-1_task-RSFC_run-1_bold.nii.gz    LAS j-  AP
sub-50015_ses-1_task-RSFC_run-1_sbref.nii.gz   LAS j-  AP
sub-50015_ses-1_task-RSFC_run-2_bold.nii.gz    LAS j   PA
sub-50015_ses-1_task-RSFC_run-2_sbref.nii.gz   LAS j   PA
sub-50015_ses-1_task-RSFC_run-3_bold.nii.gz    LAS j-  AP
sub-50015_ses-1_task-RSFC_run-3_sbref.nii.gz   LAS j-  AP
sub-50015_ses-1_task-RSFC_run-4_bold.nii.gz    LAS j   PA
sub-50015_ses-1_task-RSFC_run-4_sbref.nii.gz   LAS j   PA
sub-50015_ses-2_dir-AP_run-1_epi.nii.gz        LAS j-  AP
sub-50015_ses-2_dir-PA_run-1_epi.nii.gz        LAS j   PA
sub-50015_ses-2_task-RSFC_run-1_bold.nii.gz    LAS j-  AP
sub-50015_ses-2_task-RSFC_run-1_sbref.nii.gz   LAS j-  AP
sub-50015_ses-2_task-RSFC_run-2_bold.nii.gz    LAS j   PA
sub-50015_ses-2_task-RSFC_run-2_sbref.nii.gz   LAS j   PA
sub-50015_ses-2_task-RSFC_run-3_bold.nii.gz    LAS j-  AP
sub-50015_ses-2_task-RSFC_run-3_sbref.nii.gz   LAS j-  AP
sub-50015_ses-2_task-RSFC_run-4_bold.nii.gz    LAS j   PA
sub-50015_ses-2_task-RSFC_run-4_sbref.nii.gz   LAS j   PA

cc @mgxd @oesteban

[jxvansne]

That looks right, although Sameera can confirm for sure. We initially also thought that the PE direction may have been messed up, and actually re-processed a subject or two in fMRIprep with reverse PE from what we did initially. That made things far worse, so we concluded that we had the PE direction correct in the first place (which would anyway be consistent with the many other subjects we have whose SDC worked fine)

[sameera2004]

Hi @effigies,

Yes these phase-encoding directions are correct.
Thanks
Sameera

[oesteban]

Okay, I had an in-depth into the reports the day before yesterday and I'm going to write this comment side-by-side with the report again. This is because I think fMRIPrep is doing very well in some of the runs and pretty bad in others, which is surprising. I have to say that the reports I'm looking at only show the RSFC task (btw, I think this should be rest to be BIDS compliant).

In the beginning, I thought that there is some problem with the PhaseEncodingDirection of some runs, and I reported that impression back to the fMRIPrep team. However, when coming back to the data I realized that I was wrong and @sameera2004 and @jxvansne are indeed reporting an issue with fMRIPrep.

Before digging in, I want to comment on a @jxvansne's previous message because it is very related to my current judgment on the issue:

I believe the two images she posted were an AP and a PA phase encode direction run (we alternate directions from one run to the next, HCP style) after SDC for a single subject.

AFAIK, the HCP encodings are LR and RL. Given the LR symmetry of most of the brains, it makes sense to acquire EPI in alternate encoding polarities (and then even combine them for the case of dMRI). However, I'd be wary of doing this when the encoding is along the P/A axis. Because of the dropout, the correction can be particularly bad in one of the polarities and I believe this is happening. I will refer to this below.

Now, let's dive in (reports are referred to https://oesteban.github.io/fmriprep-issue-2210/):

• sub-50015_ses-1_task-RSFC_run-1_desc-sdc_bold.svg - This run looks great to me, within what can possibly be done using PEPOLAR methods.
• sub-50015_ses-1_task-RSFC_run-2_desc-sdc_bold.svg - Here I think is happening what I mentioned before. Overall, the low-frequency component of the fieldmap seems to be correctly estimated, and the outline of the brain and the brainstem seem to fit better the anatomical segmentation contours. However, if you look at the genu of the corpus callosum for x=1 and maybe that frontal extreme of the ventricle in x=11, you see a high-frequency component that is not correctly estimated and is pulling the subcortical structures towards the frontal pole.
• sub-50015_ses-1_task-RSFC_run-3_desc-sdc_bold.svg - Like ses-1/run-1, this one seems pretty decent to me.
• sub-50015_ses-1_task-RSFC_run-4_desc-sdc_bold.svg - Again, we have another example of this problem. If we look at the genu of the CC (x=1) you see how it is smeared towards the front. However, if you look at the very bright ventricles in z=-24, they are fitting in the contours incredibly well. Again, low-frequency components are correctly estimated, and high-frequency components not so great.
• sub-50015_ses-2_task-RSFC_run-1_desc-sdc_bold.svg - Good one.
• sub-50015_ses-2_task-RSFC_run-2_desc-sdc_bold.svg - Again, the low-frequency component is okay (could be better, TBH) and the high-frequency is very bad.
• sub-50015_ses-2_task-RSFC_run-3_desc-sdc_bold.svg - Good.
• sub-50015_ses-2_task-RSFC_run-4_desc-sdc_bold.svg - As ses-2/run-2, but even less good for the low-frequency component (see the ventricles in z=-24, one is better fit but the other is possibly worse).

So I believe the fundamental difference/problem is that HCP Pipelines is using FSL TOPUP in this estimation and fMRIPrep uses AFNI's 3dQwarp. This is not to say that TOPUP is better than 3dQwarp, this is to say that, to get visually plausible results with this PE strategy, TOPUP is better - i.e., a visually plausible result does not mean a better result per se in this case, because the dropout will really make your fMRI very dubious in the affected regions anyways. That said, we will need to push for a TOPUP implementation of SDC within fMRIPrep because it is also likely that TOPUP is not just visually better in this case.

This also leads me to say one of the principles for fMRIPrep is not to replace HCP Pipelines [or insert any other very specific pipeline here]. If you are acquiring HCP-like data, then you'll probably be better off using HCP Pipelines.

[sameera2004]

Hi @oesteban and @effigies

Thank you very much for testing our data and reporting it. Yes, we got same results when we use fMRIprep to preprocess sub-50015 data. We saw these issues on the same runs, sub-50015_ses-1_task-RSFC_run-2_desc-sdc_bold.svg and sub-50015_ses-1_task-RSFC_run-4_desc-sdc_bold.svg, sub-50015_ses-2_task-RSFC_run-2_desc-sdc_bold.svg and sub-50015_ses-2_task-RSFC_run-4_desc-sdc_bold.svg. No issues with other runs. We experienced similar issues on some other subjects/runs.

Are you planning to use FSL TOPUP implementation of SDC within fMRIprep in near future?

Best,
Sameera

[effigies]

There's ongoing work on an implementation in https://github.com/nipreps/sdcflows/pull/106.

[jxvansne]

Hi, is work on this still ongoing? I realize @effigies referenced sdcflows #106 and #117 above, but both of those are closed and I'm not sure where I should be following to be notified if TOPUP gets an fMRIprep implementation (or any other workaround for the issue we're experiencing)?

Thanks again,
-Jared

[oesteban]

Yes, SDCFlows is going through a pretty deep overhaul and that development version includes TOPUP. However, it will also require quite an effort to propagate the new API into fMRIPrep. For simplicity, we are testing the new solutions in https://github.com/nipreps/dmriprep first.