@CameronTEllis, are you in touch with anyone for a review?

Given that they are small changes, are you able to review them? Or should I ask someone more familiar with the code?

Please ask someone more familiar.

@lcnature Would it be possible for you to find time to review this please?

@lcnature, we would really appreciate your help reviewing this PR.

The recommendation to use nbstripout is not automatically enforced. While we accepted non-stripped notebooks, we very much prefer not to pollute Git with non-code data, such as images. The long-term solution to using notebooks is something like sphinx-gallery (see PR #270).

@CameronTEllis @mihaic I think the code looks pretty good to me. I checked the error message. It is from my test code of the gen_design function in utils. I think it arises because in the past, the down-sampling after convolution in fmrisim was a direct down-sampling, while it is an averaging in a time window in the new code. I included that original test code based on the old fmrisim implementation. I was not to sure exactly which approach match the MRI physics better so I asked mrihelp yesterday. Mark told me that the acquisition time of a slice (ignoring multi-band) is the TR divided by the number of slices. So for any single slice, it seems that the old approach of direct down-sampling makes sense. But for slice-timing corrected image, interpolation is involved to infer the BOLD signal at one particular time for all slices. To me, this would be a bit close to the average signal within that TR, but not exactly. So I am still not sure whether direct down-sampling or averaging within a TR is better for your code -- both make some senses. Maybe it is up to you to choose one of them. If you stick with averaging, then I think deleting Lines 124-128 (

@CameronTEllis, when implementing @lcnature's suggestion, please document your choice in a docstring.

@lcnature,

This is an interesting nuance. As I see it, the contrast is between saying the signal within a TR is either represented by the activity at single time point or it reflects the average of the activity within a TR. Practically, these are small differences because the HRF is slow relative to the length of a TR. However, as you point out there are substantive theoretical differences.

In an EPI we are indeed measuring voxel intensity in a slice by slice manner and this acquisition time is very brief. So in some sense the signal that a voxel represents should come from a single timepoint of the HRF. I believe the first approach is correct, but inaccurate if applied by downsampling as suggested.

If we downsampled by taking every nth timepoint, we would (re-)introduce a substantial bias (remember this is what we previously had as the default). For instance, if you have a signal function (estimate of the HRF convolved with the event sequence) that is 100,000 elements, and you want to down sample to 100 TRs, each of 1s duration, then you might take the 0th element, 100th element, 200th ... etc to represent your expected activity for each TR. This would effectively mean you are treating this and every voxel like it as something like the first acquired slice. You could pick a different nth element (the 50th or the 99th) but either way you are making an assumption that all signal is being acquired at the same time. This is also a problem because it can effectively shift the timing of events: sampling the 0th point is a TR different that sampling the 99th time point, so your peak response would come later in the 0th version than the 99th version.

As mentioned, we are in fact sampling data at multiple times throughout a TR. This means that taking the average of the signal function over that time window is approximately the same as if the voxels containing signal come from a distribution of voxels uniformly distributed across the slices. However, this is also often wrong, especially when we are doing whole brain acquisitions, since there is often non brain in the first and last slices. If we are simulating interleaved or multiband slice acquisition sequences this would be less of a problem but still present.

The most optimal way of designing the simulator would be to sample at a different time within a TR depending on the slice of the voxel being modelled. In other words, if the voxel is in the slice acquired first then the signal function should be sampled early in the TR, if the slice is last it should be sampled late in the TR. I could add this option to the code (provide a slice time of the voxel for each signal that is being simulated). This would be good additional functionality that would provide a lot of flexibility; although as mentioned before it would make little practical difference compared to averaging (I would default to take the middle time point).

If you agree with this logic I will go ahead and make this change.

Thanks,
Cameron

Hi @CameronTEllis
Yes, I agree that for most of the usage, the two would not make much difference. The only case that they may differ a lot is when people use super long TR in some auditory task in order to make the auditory stimuli and the scanner noise far apart in time.

And yes, if the simulator capture the acquisition time difference between slices, it would be much closer to the reality. But of course this can be left for future extension if you want to do it.

But I think the worry about choosing to sample at 0th or 99th interval is less of concern. This is because if the simulator assume the same acquisition time for all slices, the user should make sure the simulated data are used to represent slice-timing corrected data. And in practice when doing slice-timing correction, you always specify which slice all other slices are aligned to. As long as the same timing choice is made for both the simulator and the real data, choosing the 0th or 99th slice are both correct.
The down-sampling here is different from the issue with down-sampling in the previous PR. Because last time the down-sampling was done to the event-timing series, in which sparse binary events may be entirely missed by the down-sampling. But here we are down-sampling at the simulated fMRI signal, which is already a smooth signal. (and that is why down-sampling and averaging are probably very close for most of the usage cases)

Maybe another (minor) argument for using direct down-sampling here is that respiratory or heart beat noise are sampled at discrete time instead of averaged within a TR in your code, if I understand correctly. Maybe it makes slightly more sense to use the same approach for both signal and noise.

Other people behind mrihelp just stopped by my office and some suggestion they gave include checking what other packages like FSL choose to do... They are happy to answer any questions you have, btw.

@lcnature,

Thank you for the follow up. I think your suggestion to simulate the data 'as if' it is slice time corrected is a good one. FSL does slice time correction by aligning to the middle time point in a TR so if I downsample at the middle time point then this will resolve this. This is a much simpler solution than providing a slice time option.

You are right about the difference in downsampling the event time course rather than the convolved response, my mistake.

Fortunately this will be a very minor change since the mean and middle time point are equal when the trend within a TR is linear (which is not always the case but is often close).

I will implement this now and re-commit the code. I will make sure to note this in the docstring

Thanks,
Cameron

@mihaic The PR looks fine for me now! Thanks for your effort @CameronTEllis !

@lcnature, you have the right to squash and merge the commits to master yourself. First you have to update your review to formally approve the PR.