Hello everyone,
I am working on inverse problem in EEG. As seen in a paper ( Sun, R., Sohrabpour et al. Deep neural networks constrained by neural mass models improve electrophysiological source imaging of spatiotemporal brain dynamics. https://doi.org/10.1073/pnas.2201128119 ), I want to “map” EEG data from a subject on a template montage to solve the inverse problem on fsaverage head model using spherical spline interpolation.
I am using the sample dataset, more specifically the visual right stimuli data. To check the process I perform the inverse problem resolution on the subject specific head model and morph the solution on fsaverage. I then use spherical spline interpolation to interpolate the EEG data of the subject to standard_1020 electrode montage and solve the inverse problem on fsaverage head model.
When I visualize the results (see image), they are not as expected 
. It looks like the activity is moved upward in the brain and I do not understand why.
My questions are :
- Is there something I am missing in the process or an explanation for those results ?
- Is it a common technique to interpolate subject specific data to a template head model ?
Any information on this topic would be amazing, thank you very much !
Sarah
Results :
Below is the code I am running (a bit long sorry) :
- MNE version: 1.5.1
- operating system : Ubuntu 22.4.3 LTS
import os
import mne
import numpy as np
from mne.datasets import sample
from mne.channels.interpolation import _make_interpolation_matrix
#### load the subject data ####
data_path = sample.data_path()
fname_raw = data_path / "MEG" / "sample" / "sample_audvis_raw.fif"
raw = mne.io.read_raw_fif(fname_raw)
raw.pick_types(meg=False, eeg=True, stim=False, exclude=["EEG 053"]).load_data() # exclude bad channel
event_file = data_path / "MEG" / "sample" / "sample_audvis_raw-eve.fif"
events = mne.read_events(event_file)
epochs = mne.Epochs(raw, events, baseline=(None, 0.0), preload=True)
epochs.set_eeg_reference(projection=True)
reject_criteria = dict(eeg=100e-6)
epochs.drop_bad(reject=reject_criteria)
## we work on the right visual stimuli
vis_right = epochs[f"{4}"].average()
_, time_p100 = vis_right.get_peak(mode = 'pos',time_as_index=True)
#### load template data ####
standard_montage_kind = "standard_1020"
standard_montage = mne.channels.make_standard_montage(standard_montage_kind)
# remove duplicate electrodes (modified combinatorial nomenclature)
exclude_mdn = ["T3", "T4", "T5", "T6"]
ids_duplicate = [ np.where([ch == e for ch in standard_montage.ch_names])[0][0]
for e in exclude_mdn ]
ch_names = list(np.delete(standard_montage.ch_names, ids_duplicate))
ssi_info = mne.create_info(
ch_names, raw.info["sfreq"], ch_types="eeg", verbose=0
)
ssi_info.set_montage(standard_montage)
### electrode loactions ###
epos_subj = np.array([ch["loc"][:3] for ch in raw.info["chs"]])
epos_templ = np.array([ch["loc"][:3] for ch in ssi_info["chs"]])
### INTERPOLATION ###
interp_mat = _make_interpolation_matrix(epos_s, epos_t)
epoch_data_ssi = np.matmul(interp_mat, epochs._data[...,:,:])
epochs_ssi = mne.EpochsArray(epoch_data_ssi, ssi_info, events = epochs.events, tmin = epochs.tmin, event_id = epochs.event_id, baseline=epochs.baseline)
vis_right_ssi = epochs_ssi[f"{4}"].average()
### inverse problem ###
# on subject
subject = "sample"
subjects_dir = data_path / "subjects"
# build forward model
if not os.path.exists( f"./fsaverage_ori_sample_ico3-{standard_montage}-fwd.fif"):
trans = f"{data_path}/MEG/sample/sample_audvis_raw-trans.fif"
conductivity = (0.3, 0.006, 0.3)
sampling = "ico3"
model = mne.make_bem_model(
subject=subject, ico=4, conductivity=conductivity, subjects_dir=subjects_dir
)
bem = mne.make_bem_solution(model)
src = mne.setup_source_space(
subject,
spacing=sampling,
surface="white",
subjects_dir=subjects_dir,
add_dist=False,
n_jobs=-1,
verbose=0,
)
fwd = mne.make_forward_solution(
raw.info,
trans=trans,
src=src,
bem=bem,
eeg=True,
mindist=5.0,
n_jobs=-1,
verbose=0,
)
mne.write_forward_solution(
f"./fsaverage_ori_sample_ico3-{standard_montage}-fwd.fif",
fwd,
)
else:
fwd = mne.read_forward_solution(
f"./fsaverage_ori_sample_ico3-{standard_montage}-fwd.fif"
)
fwd = mne.convert_forward_solution(
fwd, surf_ori=True, force_fixed=True, use_cps=True, verbose=0
)
# for the template head model
subject = "fsaverage"
subjects_dir = os.path.dirname(mne.datasets.fetch_fsaverage(verbose=0))
if not os.path.exists( f"./fsaverage_{standard_montage_kind}_ico3-fwd.fif"):
conductivity = (0.3, 0.006, 0.3)
sampling = "ico3"
model = mne.make_bem_model(
subject="fsaverage", ico=4, conductivity=conductivity, subjects_dir=subjects_dir
)
bem = mne.make_bem_solution(model)
src = mne.setup_source_space(
"fsaverage",
spacing=sampling,
surface="white",
subjects_dir=subjects_dir,
add_dist=False,
n_jobs=-1,
verbose=0,
)
fwd_templ = mne.make_forward_solution(
ssi_info,
trans="fsaverage",
src=src,
bem=bem,
eeg=True,
mindist=5.0,
n_jobs=-1,
verbose=0,
)
mne.write_forward_solution(
f"./fsaverage_{standard_montage_kind}_ico3-fwd.fif",
fwd_templ, overwrite=True
)
else:
fwd_templ = mne.read_forward_solution(
f"./fsaverage_{standard_montage_kind}_ico3-fwd.fif"
)
fwd_templ = mne.convert_forward_solution(
fwd_templ, surf_ori=True, force_fixed=True, use_cps=True, verbose=0
)
### inverse solution ###
method = "sLORETA"
snr = 3.0
lambda2 = 1.0 / snr**2
cov = mne.compute_covariance(epochs, tmax=0.0)
inv = mne.minimum_norm.make_inverse_operator(
vis_right.info, fwd, cov, loose=0, depth=0, fixed=True, verbose=True
)
stc = mne.minimum_norm.apply_inverse(
vis_right,
inv,
lambda2,
method=method,
pick_ori=None,
return_residual=False,
verbose=True,
)
vis_right_ssi.set_eeg_reference(projection=True)
epochs_ssi.set_eeg_reference(projection=True)
cov_ssi = mne.compute_covariance(epochs_ssi, tmax=0.0)
inv_ssi = mne.minimum_norm.make_inverse_operator(
vis_right_ssi.info, fwd_templ, cov_ssi, loose=0, depth=0, fixed=True, verbose=False
)
stc_ssi = mne.minimum_norm.apply_inverse(
vis_right_ssi,
inv_ssi,
lambda2,
method=method,
pick_ori=None,
return_residual=False,
verbose=False,
)
