Movement compensation of CTF data by Maxwell filtering *without* external interference suppression?

[jeythekey]

I would like to do movement compensation of CTF raw data (several runs of length 15-20 Minutes), in order to concatenate epochs within them in sensor level for decoding analyses. I still want to be able to apply the CTF gradient compensation, since it appears to me to give superior noise suppression over pure tSSS in my dataset. I already extracted the continuous head position with mne.chpi._calculate_head_pos_ctf (#4088) by @bloyl.

However, from the old discussion in #2277 (@staulu, @larsoner, @wronk), I understand that

1. it is not possible to perform CTF compensation after SSP has been applied (by the way, also, when only applying it to the data channels?), and
2. it is not yet (?) clear, whether SSP can be applied to the gradient compensated data.

So, my questions is, if setting st_duration=None in maxwell_filter will prevent maxwell_filter from doing external interference suppression and only do the movement compensation? And whether that procedure might actually be OK (for leadfield computation etc.)? I tried to understand the code (thereby encountering some other questions/confusions), but unfortunately cannot answer this question with confidence myself.

A subsequent questions would be, if it is better to first apply movement compensation and then the gradient compensation or the other way round. (My intuition would opt for the first way.)

Thanks in advance for any thoughts/help/advice!!

[larsoner]

So, my questions is, if setting st_duration=None in maxwell_filter will prevent maxwell_filter from doing external interference suppression and only do the movement compensation?

No, it will just not do tSSS. SSS is still applied. The SSS decomposition of the data is used for movement compensation, so they can't be easily separated.

it is not yet (?) clear, whether SSP can be applied to the gradient compensated data.

Yes you should be able to do SSP after gradient compensation. (You cannot currently do SSS after gradient compensation, though. In principle it is possible but would take some thinking and testing to get it working.)

if it is better to first apply movement compensation and then the gradient compensation or the other way round.

Currently your only choice is movement compensation then gradient compensation. SSS/movement compensation after gradient compensation has not been implemented or tested.

[jeythekey]

Thanks, @larsoner! That was what I actually expected. Please allow a query: You say:

Currently your only choice is movement compensation then gradient compensation.

However, back in 2015, @staulu stated (highlights mine):

Your comment made me realize that it is actually not feasible to perform the reference sensor based compensation if SSS has already been applied to the whole CTF array (actual sensors + reference sensors) because all the reconstructed signals will be of internal origin and the assumptions that justify using external Taylor's expansion for the extrapolation of reference signals don't hold true anymore.

--> So, is this only true, if I include the ref_meg sensors in the SSS and am I fine otherwise?

He then concluded similarly to you:

[...] or 2) first perform the 3rd order gradiometrization and then possibly do SSS depending on the details of the compensation. I think this topic requires some extensive thinking and planning, maybe in a future project :)

I actually got this technically running by just commenting out the check-line (L. 281). But of course, that doesn't say anything ...

--> What concretely would need to be taken into account to come to a conclusion? Is there already anything related underway?

[larsoner]

Hmm. Actually it would be a bit weird to not include the reference sensors in the expansion, but then do include them in the gradient comp. If SSS removed all external interference from the non-ref sensors, then you did gradient compensation, you would probably end up putting it back.

The other option, SSS on all channels then compensation, I don't know about (they Taylor expansion assumption, etc.).

These are the sorts of things that would need to be discussed. But in principle the compensation-then-SSS could be implemented and tested regardless. Basically it needs to be applied to the data, then also applied to the SSS model before it's (pseude)inverted in the code. Nobody is working on this at the moment.

[jeythekey]

Thanks again for your answers. So, to summarize, as of the current state SSS and gradient comp cannot really be combined, right?

1. (t)SSS --> gradient compensation: Problematic due to Taylor assumptions probably not met anymore.
2. gradient compensation --> (t)SSS: Needs adaptions in the core of the maxwell code (applying the gradient comp also to the SSS model).

Option 2 is theoretically sound and just needs to be implemented, for option 1 further theoretical insight is needed.

Is this correct?

[larsoner]

That's my understanding

[staulu]

I basically agree. I have a couple of comments:

1. In the case of (t)SSS --> gradient compensation, one assumes that the data contain signals of external origin that can be reconstructed from the reference sensor measurement. However, if the actual MEG channels have been processed by (t)SSS, then extrapolating and subtracting external interference based on reference sensors (gradient compensation) will most probably increase interference in the (t)SSS-processed data.
2. In principle, if the gradient compensation is applied to purely subtract homogeneous fields or low-order gradients of the external field from the measurement, then one shouldn't necessarily need to modify the SSS model because such a subtraction, if done properly, shouldn't fundamentally alter the spatial pattern of the brain signal. But whether the SSS model needs to be modified or not depends on the details of the gradient compensation. Anyway, this should be a doable option.

[larsoner]

@jeythekey what you did by commenting out (a bit more than) this check is basically what (2) is without modifying the SSS basis. Can you try it out on some of your CTF datasets to see how well it works?

[jeythekey]

Thanks, @larsoner and @staulu for your helpful comments.

I now tested all the variants. I uploaded a notebook here (cleaned pdf) (the ipynb-file is here).

To summarize:

• tSSS alone (without ref_meg!) performs similarly to gradient compensation alone, but seems to suppress more higher frequency noise (at least for most channels, see below ...)
• Additional gradient compensation before tSSS removes slow drifts (0th order vs. >1)
• No difference between 1st order and 2nd order gradient comp (but there is definitely a difference without SSS.)
• Subtle difference between 1st order and 3rd order gradient comp (3rd order a little less slow drifts).

--> @staulu: Do you think, 3rd order gradient (the default) would be ok? Can this be judged from raw data figures? What tests would be needed to judge? Do we need a simulation with source reconstruction? I also pasted some info from the CTF manual to the notebook, if that should help.

---

• tSSS --> gradient compensation is, as expected, not usable.
• However, also only tSSS on the whole array (as suggested by before), including ref_meg, introduces a lot of noise to the data channels.

---

Now, the problem ...:

• Applying any destination other than None (even using raw.info['dev_head_t']['trans'][:3, 3]) during Maxwell filtering makes some - mainly occipital - sensors extremely noisy > 25 Hz (broad band). head_pos alone, with destination=None does not have this effect.

--> WHAT IS GOING ON HERE?

... And maybe a solution?

• Applying OTP before Maxwell filtering alleviates the broad-band-noise problem a lot, and completely until 45 Hz. Cool! :) Still, the question above remains.

--> Now, one additional question concerning the last step to @larsoner: In the help it is stated, that st_correlation of maxwell_filter should be increased. Do you have any experience on what would be a suitable value? Or what is the criterion to select this parameter?

[larsoner]

Applying any destination other than None (even using raw.info['dev_head_t']['trans'][:3, 3]) during Maxwell filtering makes some - mainly occipital - sensors extremely noisy > 25 Hz (broad band). head_pos alone, with destination=None does not have this effect.

Using destination=None should be the same as using destination=fname_raw (since we apparently don't support directly doing destination=raw.info['dev_head_t], though we probably should), can you check? If it isn't, then we have a bug.

If it is, then the question becomes why destination=raw.info['dev_head_t']['trans'][:3, 3] is so different. The only difference there should be the rotation matrix raw.info['dev_head_t']['trans'][:3, :3] -- is this very different from np.eye(3)? You can use something like the following to get how many degrees off it is:

np.rad2deg(mne.transforms._angle_between_quats(rot_to_quat(np.eye(3)), rot_to_quat(raw.info['dev_head_t']['trans'][:3, :3])))

But the fact that OTP seems to help a lot might mean you have a channel or two with some bad artifacts -- have you browsed the raw data to see if there are some you could mark bad (if not using OTP)?

Do you have any experience on what would be a suitable value? Or what is the criterion to select this parameter?

We are currently working on what it should be, but in general if 0.98 works well without OTP, you probably need something over 0.99 like 0.995 or 0.999. But it depends on the data.

[jeythekey]

Thanks a lot for your helpful comments!!

1. destination=raw.info['dev_head_t] is actually supported and is the same as destination=None. (So one could simply update the doc string). [But I don't see why that could be connected to the problem, since providing raw.info['dev_head_t']['trans'][:3, 3] is just a np.ndarray, not any more a Transform class, or am I mistaken?]

2. Anyway, the rotation would be 6°; I don't know if that is a lot?

3. I marked all suspicious channels by hand as 'bads', and let maxwell reconstruct them, but also excluding them before did not change anything. But based on the success of OTP I would also suspect that I missed some bads. It's only strange they do not occur at all in the spectrum of the non-maxwell-data, isn't it? (Unfortunately, I did not get autoreject to work, due to some circular ref_meg checkings, but I didn't have the time yet to look into that in detail ...)

4. I also tested this with the st_threshold of 0.999 and the noise on the same channels got a little bit worse again between 20 and 40 Hz, but still nothing compared to the noise without OTP.

So, what is your feeling, should I just go for it and pre-process all my data with OTP + movement correction? :) Or would you rather recommend to do additional tests before, like a simulation with a dummy? Or trying around with different st_thresholds? (But I don't have a clue how to judge what is a good value ...)

Again, thank you a lot.

PS: pasting from maxwell.py to illustrate the first point:

def _check_destination(destination, info, head_frame):
    """Triage our reconstruction trans."""
    if destination is None:
        return info['dev_head_t']
    # ...
    elif isinstance(destination, Transform):
        recon_trans = destination

[larsoner]

Ahh right yes the docstring is just out of date.

[But I don't see why that could be connected to the problem, since providing raw.info['dev_head_t']['trans'][:3, 3] is just a np.ndarray, not any more a Transform class, or am I mistaken?]

The only difference between giving raw.info['dev_head_t']['trans'][:3, 3] (which is problematic) and raw.info['dev_head_t'] is whether or not the rotation is present. I wanted you to try raw.info['dev_head_t'] just to make sure it indeed is okay, and equivalent to destination=None, which it thankfully is. That means that there is something about this 6 deg rotation that causes the noise increase. That's not a big rotation, so I'm surprised it makes a big difference.

How do things work without tSSS (st_duration=None) / with plain SSS? Same problem depending on whether or not destination has the rotation or not?

[staulu]

Interesting discussion going on here. Regarding the question about the order of gradient compensation and subsequent tSSS, I didn't see any significant differences in the raw data figures between, e.g., the 2nd and 3rd order. Using the default order (3rd) should be fine. A relevant test is simply comparing the spectra across different processing variants, but this is already being done based on the notebook. This investigation may not need a simulation but it would be good to use a well-defined brain response to verify that the signal of interest stays intact. Some phantom data could also be very informative, if available.