Slice-timing-correction induced time shift holds potential for error in subsequent fMRI analysis

[theoschaefer]

Hi everyone,

with this post, we would like to point to a potentially serious issue introduced by fMRIPrep’s slice-timing correction (STC) which affects further processing of fMRI data and might therefore be relevant for a large number of users.

STC in fMRIPrep is referenced to the middle slice which leads to a time-shift in the onsets by 0.5 TR. Importantly, users might not be explicitly aware of this time shift and might therefore not consider it for subsequent analyses, especially since defaults of other software packages are sometimes misaligned (e.g. nilearn’s function “nilearn.glm.first_level.FirstLevelModel'' has a default microtime onset of 0). This may inadvertently lead to model misspecification for a potentially large number of users. Currently, it is not clear how many users are affected by this issue but an informal survey of fMRIPrep users at our institute makes us think this may be quite prevalent. The impact is currently also not very clear but correcting for the effect (i.e. considering the time shift in subsequent analyses) strongly improved the results of all users involved at our institute, so we think the impact may be high. The upside of this effect is that accurately taking the time shift into account should only improve model specification and therefore lead to stronger results.

While the reference to the middle slice is currently documented at two locations in fMRIPrep (here and here), it does not appear in the usage notes or in the output documentation.

For this reason, it might first be necessary to explicitly inform users in the documentation about the time shift at many different locations and - importantly - explain the consequence for statistical analyses after fMRIPrep is completed. This could include an explicit warning box along the lines of “Caution! Slice timing correction shifts the data by half a TR. This means that you need to adjust your onset files, the microtime onset or the parameters in your first-level model function”. The information could be included in “Processing Pipeline Details”, “Outputs of fMRIprep”, “Citing fMRIPrep” and in the methods template.

Here are some additional ideas on how to deal with the time-shift issue:

1. More information in the outputs, e.g. by emphasizing the new reference time in the output or by providing new time frames in the output (e.g. [1., 3., 5., ...] instead of [0., 2., 4., …] for TR = 2)

2. Adjust the STC reference to the first slice to make fMRIPrep more modular. This has to be accompanied by additional documentation to avoid similar errors in the other direction.

3. Write out new time-shifted onset files to keep everything downwards compatible and still use the optimal processing pipeline. This also has to be accompanied by additional documentation of the output files. This may, however, currently not be possible with the current version of BIDS.

4. Include an argument in the fmriprep function which specifies the preferred reference slice.

Looking forward to hearing from you!

@martinhebart
@nitschalex
@oliver-contier
@theoschaefer

[effigies]

Hi Theo, thanks for this.

1. More information in the outputs, e.g. by emphasizing the new reference time in the output or by providing new time frames in the output (e.g. [1., 3., 5., ...] instead of [0., 2., 4., …] for TR = 2)

This seems doable. There's a toffset field in NIfTI. I'm not sure if this is respected by most tools. My initial thought had been to use slice_code to indicate how slice timing was performed, but toffset seems like the better fix.

CIFTI images do have explicit time series axes, so we should update the SeriesStart field.

GIFTI time series do not have any representation of an offset. We will need a BIDS metadata field for this.

In all cases, we need to ensure that downstream tools handle the internal and BIDS metadata correctly.

2. Adjust the STC reference to the first slice to make fMRIPrep more modular. This has to be accompanied by additional documentation to avoid similar errors in the other direction.

4. Include an argument in the fmriprep function which specifies the preferred reference slice.

4 seems to cover 2. I think that's reasonable.

3. Write out new time-shifted onset files to keep everything downwards compatible and still use the optimal processing pipeline. This also has to be accompanied by additional documentation of the output files. This may, however, currently not be possible with the current version of BIDS.

Given that fMRIPrep doesn't use or modify events.tsv, this would be a significant change in where downstream tools are expecting to find events files. I could also imagine situations where somebody preprocesses their fMRI data and then separately updates their events files; it would be burdensome to re-compute.

That said, a small tool that takes the original dataset and the fMRIPrep outputs and produces corrected events files would be sensible.

[neurolabusc]

2. Adjust the STC reference to the first slice to make fMRIPrep more modular. This has to be accompanied by additional documentation to avoid similar errors in the other direction.

Using the first slice makes the timing intuitive. However, there is a rationale for using the middle slice in the 3D volume in particular, for single-band sequential (rather than inter-leaved) acquisitions: the slice where no temporal shift is required will require no interpolation. Typically, regions of interest are near the center of our field of view, so it makes sense to manipulate this central most interesting data the least. This is somewhat related to my recent comments on slice time correction.

I do agree that the NIfTI toffset should record any change.

[effigies]

Opened discussion for BIDS metadata on bids-standard/bids-specification#836.

[theoschaefer]

Thanks for your answers and for cross-posting the issue!

We could start to write documentation snippets as we think that explicitly informing users in the documentation about the time-shift consequences might be one of the first steps. We would then push our suggestions to the repository.

Adding the information to the Nifti toffset field would certainly make the metadata more informative. How many software packages consider the toffset field automatically? Or would this still require users to actively retrieve this information out of the metadata?

We agree that changing the original events.tsv might lead to additional difficulties for downstream tools or custom scripts. Perhaps the tool you propose, Chris, could either save corrected events.tsv separately in the fMRIPrep derivatives folder or output new frame times for the acquired volumes (e.g. 1,3,5 instead of 0,2,4 s). Also for the frame times, the question would then be if this is used by software packages (e.g. the nilearn function make_first_level_design_matrix requires frame times as an explicit input) - or at least generating these outputs might make the user aware of the time shift and perhaps help them to double check their timings.

[effigies]

We could start to write documentation snippets as we think that explicitly informing users in the documentation about the time-shift consequences might be one of the first steps. We would then push our suggestions to the repository.

Sounds good to me. Please feel free to ask questions if you need help getting going on that.

Adding the information to the Nifti toffset field would certainly make the metadata more informative. How many software packages consider the toffset field automatically? Or would this still require users to actively retrieve this information out of the metadata?

We would need to look into this. We should definitely encourage tools to use toffset by default and to allow users to provide an offset. I know nilearn does not use toffset.

We agree that changing the original events.tsv might lead to additional difficulties for downstream tools or custom scripts. Perhaps the tool you propose, Chris, could either save corrected events.tsv separately in the fMRIPrep derivatives folder or output new frame times for the acquired volumes (e.g. 1,3,5 instead of 0,2,4 s).

I would personally keep them in a new derivatives dataset, to keep the provenance of the fMRIPrep derivatives clean. But clean organization is in the eye of the beholder, and this will almost certainly be a decision that will be made by the user.

Also for the frame times, the question would then be if this is used by software packages (e.g. the nilearn function make_first_level_design_matrix requires frame times as an explicit input) - or at least generating these outputs might make the user aware of the time shift and perhaps help them to double check their timings.

Yeah, user awareness seems essential. While a BIDS-aware tool should respect the BIDS metadata, it's going to take time to percolate through, and not all tools are going to be BIDS-aware.

[HippocampusGirl]

Would it be possible to also add the number of non-steady state volumes that are skipped to the toffset? Or is it meant only as the offset within a volume?

[effigies]

From the spec: The toffset field can be used to indicate a nonzero start point for the time axis. That is, time point #m is at t=toffset+m*pixdim[4] for m=0..dim[4]-1.

This seems clear to me that toffset should give the timing of the first volume, regardless of whether it's steady state or not.

If the goal is to indicate that the first steady state volume is t0, and events.tsv are defined relative to that, then it would make sense for the raw toffset to be set to -n_steadystate * pixdim[4]. And we should make sure that we do respect incoming toffset and StartTime when making these adjustments.

[kellyhennigan]

Thanks for taking this issue seriously as it does impact timing and subsequent modeling choices. Would it be possible to add which slice (or timepoint) the slice timing correction aligns as a user option? So that users can choose for example to align their data to the first or last slice, rather than the middle slice?

[effigies]

AFNI's 3dTshift allows a lot of flexibility with -tzero or -slice:

  -tzero zzz    = align each slice to time offset 'zzz';
                  the value of 'zzz' must be between the
                  minimum and maximum slice temporal offsets.
            N.B.: The default alignment time is the average
                  of the 'tpattern' values (either from the
                  dataset header or from the -tpattern option)

  -slice nnn    = align each slice to the time offset of slice
                  number 'nnn' - only one of the -tzero and
                  -slice options can be used.

I suspect we don't want to provide this much flexibility, but what options should we allow? I'm a little inclined to only allow the first slice or the middle of acquisition. Are there other use cases we should seriously consider?

To make a concrete proposal, I'd say let's create a new --stc-target flag with possible values first and mid. (Not really attached to any of these, so happy for others to propose alternatives.)

[effigies]

In the absence of feedback, I'm going to go with --slice-time-ref {first,mid}, defaulting to mid. The alternative would be to make it a float from 0-1, following nilearn's:

    slice_time_ref : float, optional
        This parameter indicates the time of the reference slice used in the
        slice timing preprocessing step of the experimental runs. It is
        expressed as a percentage of the t_r (time repetition), so it can have
        values between 0. and 1. Default=0.

[martinhebart]

Sorry, I’ve been a bit more passive on this thread.

I really like both suggestions. I think I’d either keep it as no free parameter to choose from or provide full flexibility. Thus, from your suggestions I’d probably slightly prefer the alternative that’s similar to nilearn’s option (with a default of 0.5?)

[effigies]

That works for me. I was slightly leaning toward only providing the two options, because I can't think of a reason I'd ever want to use a different value. But even a little weight on the 0-1 scale is enough to tip me to the other side.

[martinhebart]

I think other reference time points can become relevant if people have a region of interest they care about most - then they could minimize the interpolation error by choosing the slice as a reference time point that covers that ROI. However, in practice almost nobody does this, so it’s a special and rare use case, so perhaps your choice of two parameters (or keeping it a fixed parameter?) is good for now. My guess is it depends whether you ever want to extend it to more than two cases? If yes, I’d opt for full flexibility. If no, then two parameters.

[martinhebart]

Having thought about it a bit more, I’m with your original suggestion of the two options {first, middle}. If people really need more flexibility, one could one day still additionally add the float option.

[effigies]

I'd rather keep the flag options the same as much as possible once they're added, so with a convincing (if rare) use case for other options, I'm inclined to go with the float. I will add documentation that lets people know they almost always want 0 or 0.5.

Edit: One thing we can do is have first or middle be aliases for 0 and 0.5.

[martinhebart]

Yes that’s what I meant (with the alias). It wouldn’t be a change, only an addition. But either ways are fine, flag {mid,first} is probably less confusing, whereas 0-1 gives full flexibility.

edit: your suggestion of clear documentation in terms of what people would almost always want is also great!

[HippocampusGirl]

Is there any chance that 0.5 could refer to different slices on different systems depending on how the float is parsed/slice number is rounded?

[martinhebart]

My understanding is that it would be aligned to the middle of the TR, so no slice reference is chosen but a time reference. But you’re right: it’s important to communicate that it’s not automatically the middle slice if that’s someone cares about.

[kellyhennigan]

I vote for having the full float 0-1 option (with a default of .5). I really appreciate having the option to align to the first slice (or really time 0) as this is the convention we've been using in our lab. Thank you!

…

On Thu, Aug 26, 2021 at 12:07 PM martinhebart ***@***.***> wrote: My understanding is that it would be aligned to the middle of the TR, so no slice reference is chosen but a time reference. But you’re right: it’s important to communicate that it’s not automatically the middle slice if that’s someone cares about. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2477 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMT5NNRP6XKYWW4C6XSRULT62GFNANCNFSM5A62JUPQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[effigies]

This is worth thinking through. What AFNI 3dTshift does by default is take the tpattern (BIDS SliceTiming array) and take the average:

       TS_tzero = 0.0 ;
       for( ii=0 ; ii < nzz ; ii++ ) TS_tzero += TS_tpat[ii] ;
       TS_tzero /= nzz ;

If we can assume slices are evenly spaced (and further that the first slice is at 0), then:

tzero = min(SliceTiming) + 0.5 * (max(SliceTiming) - min(SliceTiming)) =  0.5 * max(SliceTiming)

If the main assumption is violated, then there's no linear transformation that preserves the average at 0.5. To preserve backwards-compatibility in the LTS series (20.2.x), we'll need to check 0.5 and allow AFNI to find the average. Going forward, I suspect we want to explicitly set tzero to be the midpoint of the SliceTiming range.

Thoughts?

[martinhebart]

It’s probably even more accurate to use the midpoint of the SliceTiming range when dealing with sparse-sampling designs (where TA is often widely different from TR - used in auditory studies), where the middle of the TR is not even the time point one would want to interpolate to.

[neurolabusc]

As an aside, the long TRs associated with sparse MRI may not be suitable for slice time correction. The temporal interpolation may often be at a longer frequency than the HRF you are trying to capture.

[effigies]

All: I would appreciate a review of #2520.

[martinhebart]

I left some comments!

[effigies]

I believe this has been substantially addressed, and I will close it now. We're always open to additional PRs for improved documentation, metadata, or corrected behavior.

Thanks again for bringing driving this discussion!