Using synthetic field map from synbold-disco with fmriprep

[bryjack0890]

I’m working with an older clinical fMRI resting-state dataset that does not contain field maps. I would like to correct for distortions but do not like the syn-sdc option in fmriprep (produced wonky images) and have thus opted to use synbold-disco GitHub - MASILab/SynBOLD-DisCo: Synthetic BOLD images for distortion correction of fMRI without additional calibration scans to create synthetic field maps .I want to further preprocess the data using fmriprep and am unsure how best to do that.

The tool says it produces the undistorted image using topup by setting the synthetic image’s bandwidth to infinite. I however do not know how to do this (or if it’s possible) for fmriprep. Does anyone have a recommendation about how I might be able to marry these techniques?

The tool also creates undistorted 4D BOLD images. However, using fmriprep presumably requires me to ignore the undistorted 4D BOLD images that are outputted by synbold-disco and instead use the “BOLD_s_3D.nii.gz: Synthetic BOLD native space” synthetic image output as a field map. This is mostly a guess however so I’m also curious about whether it would be incorrect for me to just feed the undistorted image into fmriprep. At that point, the 4D BOLD would be motion corrected but other than worrying about the motion parameters during denoising, which would be taken care of after fmriprep, I don’t see an issue. Am I over/under-thinking this? Has anyone tried to do something similar before?

Thanks in advance. -Bryan

[jsein]

Hi,

You raise an interesting question here: how to tweak fmriprep to make it use SynBOLD-DISCO without knowing it?
From the top of my head, I see a couple of options that could be interesting to test to see if that could work:

• feed FMRIPREP with the undistorded, motion corrected BOLD_u.nii.gz (4D timeserie corrected by SynBOLD-DISCO) and tell fmriprep to not run SDC (susceptibility distortion correction) on this data. That would work, but with the disadvantage of interpolating twice (= more smoothing) the data (once with SynBOLD-DISCO and once with fmriprep). Also motion correction will be done twice in that case: once with SynBOLD-DISCO and once with FMRIPREP, so you could not use the motion parameters calculated by fmriprep but instead need to use the motion parameters from Syn-BOLD-DISCO for your GLM and motion outliers. Or you could try to not run motion correction with SynB0-DISCO but I guess that would affect the accuracy of SDC with SynB0-DISCO.

• Feed FMRIPREP with the native BOLD image and give in the fmap/ folder: BOLD_s_3D.nii.gz and BOLD_d_3D_smooth.nii.gz images as being the blip-up and blip-down images to be used with the PEPOLAR method. The PEPOLAR method is calling topup, as does SynBOLD-DISCO, but I am not sure if SDCFlows (used by FMRIPREP for SDC) would understand how to fill correctly the acquisition parameters file for topup in that case.

• Feed FMRIPREP with the native BOLD image and give in the fmap/ folder: the fieldmap in Hz and undistorded magnitude image generated by topup run by SYNB0-DISCO. In that case FMRIPREP would use the Direct B0 mapping method for SDC. That could work! To be tested.

EDIT: for the last method, SynBOLD-DISCO does not output a fieldmap from topup , one would need to run topup again with the output of SynBOLD-DISCO with the option --fout and giving in acqparams.txt the correct value for the total readout time of the bold serie. To be tested, I am not sure it will work.

[bryjack0890]

The first two options were what came to mind for me and didn’t seem reasonable given 1) the additional smoothing/mc and 2) not knowing how fmriprep’s topup would perform by adding them to the fmap folder.

#3 however is very interesting and didn’t cross my mind. I’ll test that out over the weekend and will report back.

[jsein]

I did make a try with option 3 above, for a dataset where I also add acquired two SE-EPI images with opposite phase encoding direction for SDC. Thus I can compare the FMRIPREP pipeline with Synbold-DISCO to the original (gold standard method) with the two SE-EPI images for SDC.

I ran one subject, as a proof of principle and I tend to prefer the images preprocessed with Synbold-DISCO! This last statement needs much more exploration, with quantitative evaluation on more subjects, but it looks interesting.

Here is what I did:

• First I re-ran topup on the output of SynDISCO folder, with the option -m myfield and feeding it with a corrected ‘acqparams.txt’ with the correct total readout time value in the first line, 4th column.
• Then I remove the 2 SE-EPI images from the BIDS fmap folder for my subject, and added the file BOLD_s_3D.nii.gz (renamed as sub-XX_magnitude.nii.gz) and myfield.nii.gz (renamed as sub-XX_fieldmap.nii.gz), and I kept one of the _epi.json file that I renamed sub-XX_fieldmap.json, adding into it the field: Units: Hz, and keeping the field IntendedFor that was already there.
• I then launched FMRIPREP as usual, starting from a clean temporary folder to be sure to get a correct comparison with the FMRIPREP execution using the pair of SE-EPI images for SDC.

=> it worked! FMRIPREP did use the new fieldmap calculated from Synbold-DISCO outputs, doing SDC with the “Direct Fieldmap” method.

Here is a visual comparison between the output of FMRIPREP with 2SE-EPI images for SDC, or with the fieldmap calculated by Synbold-DISCO outputs.

Here is the fmriprep (v23.0.0) command I used each time:

singularity run -B /scratch/jsein/BIDS:/work,$HOME/.templateflow:/opt/templateflow --cleanenv /scratch/jsein/my_images/fmriprep-23.0.0.simg \
		 --fs-license-file /work/freesurfer/license.txt /work/$study /work/$study/derivatives/fmriprep  \
		 participant --participant-label $sub \
		 -w /work/temp_data_${study}\
		 --mem-mb 50000 --omp-nthreads 10 --nthreads 12  \
		 --fd-spike-threshold 0.5 --dvars-spike-threshold 2.0 --bold2t1w-dof 9  \
		 --output-spaces MNI152NLin2009cAsym T1w --ignore slicetiming sbref  \
		 --fs-subjects-dir /work/$study/derivatives/fmriprep/sourcedata/freesurfer --debug compcor --track-carbon --country-code FRA

Fieldmap estimation:

• Method with Synbold-DISCO (“Direct Fieldmap”):
  

  

• Method with two SE-EPI images (PEPOLAR):
  

  

At first sight, the fieldmap estimated by the PEPOLAR method with the 2 SE-EPI images makes more sense.

Interestingly, here is the output of the bold-to-T1w registration:

• Method with Synbold-DISCO (“Direct Fieldmap”):
  

  

• Method with two SE-EPI images (PEPOLAR):
  

  

I think it is already a good news: for your dataset with no fieldmap, it looks like you may be able to use the Synbold-DISCO method and use the OUTPUTS for fmriprep!

I am currently running the method 2 (giving FMRIPREP the synthetized BOLD image from Synbold-DISCO with a total readout time of 0 for usage with the PEPOLAR method), I will keep you updated on how it went.

[schillkg]

Hi @jsein and @bryjack0890 - very interesting discussion - We are happy to see some critical evaluation of synbold-disco. We want to continually evaluate and improve this algorithm!

It seems this algorithm may be most useful to the field (or better-integrated into fmriprep) if we included the field map as output using --fout? If so, we can recompile the container with this as an output.

Our expertise is in diffusion MRI, with less experience with distortion correction (and correction of other artifacts) for fMRI. Happy to also make any other suggested changes to enable better distortion correction and validation.

Thank you,
Kurt

PS - we would always advocate for the use of a field map or reverse PE pair if those are available before using our algorithm. We envision it’s use in the case where these are unavailable and you want to use physics-based correction (i.e. topup) instead of a simple nonlinear registration.

[jsein]

Hi @schillkg , welcome to Neurostars and many thanks for sharing those nice tools: synbold-disco and synb0-disco to the community!

To continue the discussion, here are the results with method 2: using FMRIPREP PEPOLAR method to run directly topup on the two images provided by synbold-disco.

First off, one thing to know is that a "TotalReadoutTime" of 0 is not accepted through FMRIPREP for the _epi images to be used for SDC in the fmap/ folder. So I gave a "TotalReadoutTime" of 0.000001 seconds to the synthetic image BOLD_s_3D.nii.gz) that I renamed sub-SUB_dir-AP_epi.json and I gave it "PhaseEncodingDirection": "j-" since my distorted bold images were with "PhaseEncodingDirection": "j"

FMRIPREP ran with no error with those images but unfortunately the corrected bold images are not great.

Here is an overview of the output images:

Estimated fieldmap:

SDC corrected bold images:

Registration of bold image to T1w image:

I tracked down a bit what went wary, and here is what I saw:

After topup run bold-disco style:
Corrected bold image:

Fieldmap image estimated by topup:

After topup run by FMRIPREP: (same image scaling for display as above)
Corrected bold image:

Estimated fieldmap:

=> the corrected image seems similarly corrected but the estimated fieldmap is drastically different! This is something to track down further…

For reference: topup command ran in synbold-disco style:

singularity exec -e -B /scratch/jsein/BIDS/$study \
/scratch/jsein/my_images/fmriprep-23.0.0.simg \
topup --imain=/scratch/jsein/BIDS/$study/derivatives/synbold/sub-${sub}/BOLD_all.nii.gz \
--datain=/scratch/jsein/BIDS/$study/derivatives/synbold/sub-${sub}/acqparamsJS.txt \
--config=b02b0.cnf --subsamp=1,1,1,1,1,1,1,1,1 --miter=10,10,10,10,10,20,20,30,30 \
--lambda=0.00033,0.000067,0.0000067,0.000001,0.00000033,0.000000033,0.0000000033,0.000000000033,0.00000000000067 --scale=0 \
--out=/scratch/jsein/BIDS/$study/derivatives/synbold/sub-${sub}/my_topup_results \
--fout=/scratch/jsein/BIDS/$study/derivatives/synbold/sub-${sub}/my_field \
--iout=/scratch/jsein/BIDS/$study/derivatives/synbold/sub-${sub}/my_unwarped_images \
--jacout=jac --rbmout=xfm --dfout=warpfield -v

Topup command run by FMRIPREP:

topup --config=/opt/conda/lib/python3.9/site-packages/sdcflows/data/flirtsch/b02b0.cnf \
--datain=/work/temp_data_NEMO_DISCOPE/fmriprep_23_0_wf/single_subject_pilote1_wf/fmap_preproc_wf/wf_B0map/topup/sub-pilote1_dir-PA_epi_regrid001_merged_sliced_volreg_encfile.txt\
--imain=/work/temp_data_NEMO_DISCOPE/fmriprep_23_0_wf/single_subject_pilote1_wf/fmap_preproc_wf/wf_B0map/setwise_avg/sub-pilote1_dir-PA_epi_regrid001_merged_sliced_volreg.nii.gz\
--out=sub-pilote1_dir-PA_epi_regrid001_merged_sliced_volreg_base\
--iout=sub-pilote1_dir-PA_epi_regrid001_merged_sliced_volreg_corrected.nii.gz\
--fout=sub-pilote1_dir-PA_epi_regrid001_merged_sliced_volreg_field.nii.gz\
--jacout=jac\
--logout=sub-pilote1_dir-PA_epi_regrid001_merged_sliced_volreg_topup.log\
 --rbmout=xfm --dfout=warpfield