Hello!
Context
I’m studying an MEG dataset with an n-back task. Without delving too much into the experimental paradigm, here are the preprocessing steps (source article: 10.1038/s42003-023-05448-z ) :
- Maxwell filtering (already done on the data shared with me)
- Notch filter (50 Hz)
- Band pass filter (1-45 Hz)
- ICA denoising of eye and cardiac artefacts
- Downsampling to 250 Hz (from 1kHz)
- Source reconstruction with LCMV beamforming (forward model: single shell with a projection to 5mm grid)
- Atlasing
- Symmetric source orthogonalization
Their preprocessing used MATLAB libraries, but I aimed to reproduce the steps in Python with MNE.
The issue
The lab kindly shared with me both processed data and unprocessed data. Here is the PSD of the channels in their preprocessed data for a single participant:
Remark the dip at 50 Hz, which is of course due to the notch filter.
When applying all steps up to (and including) downsampling in my own pipeline:
raw_icaed = mne.io.read_raw("/path/to/subject/raw_nBack_tsss_mc_annot_filtered_icaed_annotated.fif")
raw_icaed.compute_psd(method='welch',fmax=120).plot()
The notch filter is preserved, as expected.
But when looking at the PSD of the source-reconstructed data (without any atlasing or symmetric source orthogonalization, since they should not affect frequency content I believe):
info_ = mne.create_info(['Reg_' + str(i) for i in range(1,stc.data.shape[0]+1)], 250, ch_types='bio', verbose=None)
raw_stc = mne.io.RawArray(stc.data, info_)
raw_stc.compute_psd(picks='bio',method='welch').plot(picks='bio')
This is precisely my issue: I expected the frequency content to be (somewhat) preserved, but it is clear this is not true. This discrepancy is likely explaining why I cannot replicate results.
My question: how can I preserve frequencies and tend to the results from the article? (ie, what did I do wrong?)
Details of beamforming
The beamforming used in the article uses Bayesian PCA to estimate covariance matrices.
Since such a beamforming is not implemented in MNE Python, I used an empty room recording, kindly provided by the lab.
raw_icaed = mne.io.read_raw("/path/to/subject/raw_nBack_tsss_mc_annot_filtered_icaed_annotated.fif")
raw_er = mne.io.read_raw("/path/to/subject/er_sss.fif")
noise_cov_raw = mne.compute_raw_covariance(raw_er,method='ledoit_wolf', rank='info')
noise_cov_raw.plot(raw_icaed.info)
data_cov_raw = mne.compute_raw_covariance(raw_icaed,reject_by_annotation=True, method='ledoit_wolf', rank={'meg':67})
data_cov_raw.plot(raw_icaed.info)
Here is how the forward model is setup:
info = raw_icaed.info
info = mne.pick_info(info, mne.pick_types(info, meg=True, eeg=False, exclude=[]))
vol_src = mne.setup_volume_source_space(subject_name,
subjects_dir=subjects_mri_dir,
surface= Path(subjects_mri_dir, subject_name, "bem", "inner_skull.surf"))
model = mne.make_bem_model(subject=subject_name, subjects_dir=subjects_mri_dir, ico=4, conductivity=(0.33,))
bem_sol = mne.make_bem_solution(model)
fwd = mne.make_forward_solution(info, fname_trans, vol_src, bem_sol)
and here is the beamforming:
# Compute source rec result
filters = mne.beamformer.make_lcmv(raw_icaed.info, forward=fwd, data_cov=data_cov_raw, noise_cov=noise_cov_raw, rank={'meg':72},reg=0.05, reduce_rank=True)
stc = mne.beamformer.apply_lcmv_raw(raw_icaed, filters)
EDIT:
This behaviour is due to the aggregation of (X,Y,Z) components during beamforming.
The beamforming solution in vector space preserves frequencies:
However the default behaviour of the beamformer is to aggregate by Euclidean norm, which then leads to frequency loss:
Using max-power, by contrast, preserves frequency profiles:
Might be worth mentioning somewhere, if relevant in general cases!
System informations
- MNE version: 1.8.0
- operating system: Ubuntu 22.04.5 LTS
Have a nice day!