MEG source reconstruction and sensor topography don't agree

[yaqing]

Dear List,

I’ve been puzzled by a problem with source reconstruction for a while and hope someone here can share their intuition.

My experiment involves comprehension of short spoken sentences. As a sanity check for my data, I’m plotting evoked responses at sentence onset. For the sensor-space data, every participant seems to have an activation in bilateral temporal sensors between 140 to 200ms, even if their butterfly plots don’t look very neat, like the two examples below. For the majority of participants, source localization of this evoked response reveals peak activation in the auditory areas as expected (Example 1). However, for a few participants (Example 2), the reconstructed source activations 1) do not show the onset response in temporal areas, 2) are often sporadic, and importantly 3) have very high amplitude compared to those with a clear auditory onset.

I have tried going back and forth several times with preprocessing, especially ICA, but it doesn’t seem to improve.

For reference, here’s a brief outline of my pipeline (MEGIN system, data filtered online between 0.1 to 330Hz during acquisition, Fs=1kHz):

A. Sensor-space preprocessing
–filter continuous HPI
–maxfilter and project head_device_t to a common head
–remove bad data chunks by visual inspection
–filter sensor space data between 0.1 to 40Hz
–ICA (picard, n_component=0.98) using data filtered between 1 to 40 Hz to remove EKG and EOG components, then apply the ICA result to the filtered signal above
–epoch at sentence onset between -0.3 to 1.2s, baseline correction by subtracting the mean value before 0s, and decimate by a factor of 5
–average across sentences to create evoked response

B. Main steps in source reconstruction
–compute noise covariance from 3 minutes of empty room recording (ER preprocessed using the same pipeline as experimental data, including ICA)
–create forward solution then inverse operator using loose=‘auto’, depth=0.8, spacing=oct6
–apply the inverse operator to the evoked response using the dSPM method with SNR=3, lambda2=1/SNR^2
–morph to fsaverage

Has anyone had a similar problem before, or have an idea why I get results like Example 2? For the moment I exclude such participants in my source-space group analysis because their super-high activation level will bias the results. But I hope I can find a solution to include them in all analysis since their sensor-space data don’t look particularly wrong.

Thanks,
Yaqing

Example 1, not super clean sensor-space response but reasonably good source reconstruction:

Example 2, better-looking sensor-space onset response but source reconstruction don’t seem to reflect it:

[richard (Richard Höchenberger)]

Hello,

I wonder if in the 2nd example, the regularization didn’t work.

You should inspect the noise covariance matrices and whitened data by:

• plotting them
• plotting the whitened evoked data

Please feel free to share the figures you get so we can take a look.

Good luck,
Richard

[yaqing]

Hi Richard,

Thanks a lot for your suggestion. Indeed, the regularization in Example 2 doesn’t seem to work. The covariance matrices look similar but the whitened signal in Example 2 looks way off.

Is there a good way to fix this?

Example 1, where source reconstruction was good:

Example 2, where it didn’t work:

[richard (Richard Höchenberger)]

@larsoner @agramfort
In the above example, it appears that noise covariance rank estimation did work, but whitening of the evoked signal did not. Any suggestions what could be the cause of this? Would trying out another regularization method potentially help? The data appears to have been Maxwell-filtered.

@yaqing Could you try calling compute_covariance() with method="auto"?

Richard

[yaqing]

Thanks, if I understand correctly, compute_covariance() takes the resting state (here my baseline) for computing the noise covariance. I followed the tutorial you linked above and the whitened data now looks good for Example 2.

I read that using baseline for the noise covariance suppresses activities resembling baseline in the source estimate. Does this indicate that some activities (oscillations?) in the baseline are causing the problem? Is it generally recommended to use resting state rather than empty room for noise covariance if I want to look at evoked data?

To complicate the question even more, the response I’m interested in is not really the onset but rather words within sentences (participants with problematic onset source reconstruction also have “bad” source reconstruction with these words). The baseline I’m using to calculate the evoked response to these words are response immediately before their onset, which are not technically resting states.

Example 1:

Example 2:

[mscheltienne (Mathieu Scheltienne)]

I think you should be able to estimate the noise covariance from both an empty-room recording or from the pre-stimulus periods. Is that the only point you changed between both processing, to get the 2 different whitened data plots?

[richard (Richard Höchenberger)]

@yaqing

Thanks, if I understand correctly, compute_covariance() takes the resting state (here my baseline) for computing the noise covariance. I followed the tutorial you linked above and the whitened data now looks good for Example 2.

The pre-stimulus period is not resting-state. It’s just that: the pre-stimulus period.

You can calculate the covariance of these segments and use them as “noise covariance” to decorrelate and scale the data (whitening).

@yaqing

I read that using baseline for the noise covariance suppresses activities resembling baseline in the source estimate.

Yes; this is often exactly what is desired.

@yaqing

Does this indicate that some activities (oscillations?) in the baseline are causing the problem?

I actually don’t know how you’re calculating the noise covariance now (please share your code), so I cannot answer this.

@yaqing

Is it generally recommended to use resting state rather than empty room for noise covariance if I want to look at evoked data?

I cannot give you a definitive answer here, but usually, when looking at evoked responses, using a pre-stimulus baseline period to estimate the noise covariance will yield better whitening results compared to using empty-room data.

@yaqing

To complicate the question even more, the response I’m interested in is not really the onset but rather words within sentences (participants with problematic onset source reconstruction also have “bad” source reconstruction with these words). The baseline I’m using to calculate the evoked response to these words are response immediately before their onset, which are not technically resting states.

Again, you’re not measuring resting-state activity in your paradigm anyway.

What you could do is take a time period before the participants get active in each trial. But I’m not sure I fully understand your experimental paradigm, so …

Best wishes,
Richard

[yaqing]

Sorry for the confusion. Yes, it should be prestimulus period not resting state.

My experiment paradigm is, briefly, participants listen to simple sentences (1 sentence=1 trial) one after another, separated by ~2s inter-trial interval with no acoustic input. Throughout the experiment they look at a fixation screen. What I’m showing above is evoked response from -0.3 to 1.2s relative to the sentence onset, using -0.3-0s as the baseline. In this case, baseline period is the silent ITI without acoustic input. To get such evoked response, I do:

epochs_g1 = mne.Epochs(raw1, events1, tmin=-0.3, tmax=1.2, event_id=event_IDs, baseline=(None, 0), picks="meg", preload=True, decim=5)
evoked_g1 = epochs_g1[cond1].average()

(because each trial has multiple triggers which I will explain below, cond1 is to choose sentence onsets)

In the original post, I computed noise covariance from 3 minutes of empty room recording before the participant arrived. After preprocessing the ER recording in the same way as my data (Maxfilter, bandpass filter, ICA), I do:

noise_cov = mne.compute_raw_covariance(er_ICAed, tmin=0, tmax=None)
noise_cov.plot(er_ICAed.info, proj=True)
evoked_g1.plot_white(noise_cov, time_unit="s")

Here er_ICAed is the preprocessed ER. This gave me the problematic whitened data as well as source estimation. Following your suggestion, I tried computing noise covariance from the prestimulus baseline:

cov_bl = mne.compute_covariance(epochs_g1[cond1], tmax=0)
cov_bl.plot(epochs_g1[cond1].info, proj=True)
evoked_g1.plot_white(cov_bl, time_unit = 's')

And now the whitened data looks normal. So @mscheltienne if I understand your question correctly, yes that’s the only thing I changed to get the two different whitened data. Now my question 1 is, why doesn’t it work with noise covariance computed from the ER, and is there a way to tune the regularization parameter to make it work?

A following question concerns a different type of epoch. In each of the sentence, I have two critical words, say the 4th and the 8th words. I want to compare evoked response to these critical words across conditions. There are two ways to baseline-correct the critical word epochs that both seem reasonable:

1. use the prestimulus period as the baseline. This way I get to find effects before word onsets, e.g. predictions.
2. use the period before word onsets (i.e. part of response to the 3rd or the 7th word) as baseline. This is probably gives a more “classic” evoked response such as N400.

I have not yet decided which baseline correction is better suited for my analysis. My question 2 is, is one method of noise covariance better than the other, depending on the baseline I choose for the evoked response?

[richard (Richard Höchenberger)]

Hello, intuitively, I would use the pre-stimulus period as baseline. If you use the period just before the critical words, you may remove relevant brain activity. But you can try both approaches and check what the results look like

I’ve never used empty-room recordings as a source of noise when analyzing evoked paradigms, so I cannot really help with your question regarding this. I wonder if @sophie might have an opinion here?

Best wishes,
Richard

[mscheltienne (Mathieu Scheltienne)]

As we already discussed, intuitively, I would also use the pre-stimulus period as baseline.

You do have other recordings where it worked with noise covariance computed from the ER? If so, could you check the whitened data as well? If it looks correct in your other recording, it would hint to a bug within MNE that should be fixed.

Mathieu

[yaqing]

Yes, the whitened data using the ER for noise covariance looked normal on some participants but not on others. The two examples are shown above, after Richard’s first response. The noise covariance matrices look similar between the two, and the sensor-space signal looks (to my bare eyes) even better in the Ex.2 where source reconstruction didn’t work. Both the abnormally high amplitude in reconstruction and the oscillatory pattern in the whitened data seem to hint towards numerically unstable solution…

[richard (Richard Höchenberger)]

@yaqing

Both the abnormally high amplitude in reconstruction and the oscillatory pattern in the whitened data seem to hint towards numerically unstable solution…

@larsoner might (not?) like this!

[mscheltienne (Mathieu Scheltienne)]

I’ve requested the problematic files from @yaqing I’ll try to reproduce locally first and will follow-up on an gh issue. Thanks @richard for the feedback!

[richard (Richard Höchenberger)]

@mscheltienne Great! I suggest you reach out to Eric or Alex in case you get stuck.

Best of luck!
Richard

[yaqing]

A quick update: I looked into another participant whose source reconstruction failed. However, this participant’s whitened data doesn’t look bad to me… I will share both participants’ data with @mscheltienne

[richard (Richard Höchenberger)]

@yaqing When creating the inverse operator, do you pass rank="info"? If not, could you try to see if this helps? Also pass this when creating the noise covariance matrix.

[sophie (Sophie Herbst)]

@yaqing my first choice would also be to compute the noise covariance from the pre-stimulus time window, as these data are closer to your data of interest compared to empty rooms. Resting state recordings would also work.
However, I would only use silent periods, not ones with spoken words for computing the noise cov, even if those words are not the ones of interest. Your 2sec window before the sentence onset seems well suited for that.

The regularization problems sounds like what I encountered previously, resulting from a mismatch in rank between the data on which I computed the noise cov and the data of interest. To solve this, provide the info of the data on which you computed the noise cov to the inverse operator function and pass rank=“info”.

[yaqing]

Thank you for the suggestion, also @richard

Passing rank=“info” seems to help to some extent but not completely resolving the issue with the empty room covariance, at least for the example participant 2 above.

The source reconstruction now has a much smaller magnitude than before (~200 compared to >10e3) and a hint of auditory response, but still very different from those “normal” participants whose magnitudes are below 100. The two figures below show the result from this participant using the new inverse operator. There are two panels for two conditions (for simplicity, the examples above are all from condition 1. The auditory stimuli are identical in these two conditions, and the raw data went through identical preprocessing). There are still oscillatory patterns in the reconstructed source activity.

The whitened data still look wrong, although the magnitude also reduced compared to using the old empty room covariance:

It seems to me that using the prestimulus period for noise covariance is the way to go. I will try it out next.

[richard (Richard Höchenberger)]

@yaqing Also please share the code you’re using to create the noise cov (as you did above) and inverse operator. Otherwise we’re left to do guesswork here

[yaqing]

Yes, to compute the noise covariance from empty room with rank info, I now do:

noise_cov = mne.compute_raw_covariance(er_ICAed, tmin=0, tmax=None, rank='info')

Then, to compute the inverse operator also with rank info, briefly:

info = mne.io.read_info(raw_fname)
src = mne.setup_source_space(subject, spacing="oct6", add_dist="patch", subjects_dir=subjects_dir)
fwd = mne.make_forward_solution(raw_fname, trans=trans, src=src, bem=bem, meg=True, eeg=False, mindist=1.0, n_jobs=-1, verbose=True)
inverse_operator = mne.minimum_norm.make_inverse_operator(info, fwd, noise_cov, loose='auto', depth=0.8, rank='info')

Here raw_fname points to the preprocessed raw signal, that is, Maxfilter+band pass filter+ICA.