Hell MNE forum!
With a student, Anand, we are comparing different source reconstruction methods for simulated dipoles in auditory areas. While it works with e/sLoreta, we see are systematic displacement with MNE, and some displacement with dSPM. Can you help us figure out what the problem is? We already tried playing with the parameters of the inverse. I am pasting Anand’s code below.
Our code:
# %%
#From https://mne.tools/stable/auto_examples/simulation/simulate_raw_data.html#sphx-glr-auto-examples-simulation-simulate-raw-data-py
import matplotlib.pyplot as plt
import numpy as np
import mne
from mne import Epochs, compute_covariance, find_events, make_ad_hoc_cov
from mne.datasets import sample
from mne.simulation import (
add_ecg,
add_eog,
add_noise,
simulate_raw,
simulate_sparse_stc,
)
from mne.minimum_norm import apply_inverse, make_inverse_operator
from mne.beamformer import apply_lcmv, make_lcmv
from mne.datasets import fetch_fsaverage, sample
# ### Loading data for experimental setup from mne sample
#
# %%
data_path = sample.data_path()
meg_path = data_path / "MEG" / "sample"
raw_fname = meg_path / "sample_audvis_raw.fif"
fwd_fname = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
subjects_dir = data_path / "subjects"
raw = mne.io.read_raw_fif(raw_fname)
noise_cov=mne.read_cov(meg_path / "sample_audvis-cov.fif")
subject = "sample"
annot = "aparc.a2009s"
sfreq =raw.info["sfreq"]
tstep = 1 /sfreq
fwd = mne.read_forward_solution(fwd_fname)
src = fwd["src"]
# ## parameter setup
# %%
epoch_duration = 1 # duration of each epoch
times = raw.times[: int(raw.info["sfreq"] * epoch_duration)]
n_dipoles = 19 # number of dipoles to create
extent = 0.0 #exten=0 means a single dipole, try to create a sequentially larger sources later
rng = np.random.RandomState(0) # random state (make reproducible)
# ### Selecting the labels for vertices of dipoles
#
# %%
label_primary_aud_cortex = mne.read_labels_from_annot(
subject, annot, subjects_dir=subjects_dir, regexp="G_temp_sup-G_T_transv-lh", verbose=False )# can also use G_occipital_middle-lh , G_temp_sup-G_T_transv-lh
label_primary_aud_cortex = label_primary_aud_cortex*n_dipoles
vertices_in_aud_cortex =np.intersect1d(src[0]["vertno"], label_primary_aud_cortex[0].vertices) #
indices_in_label = np.where(np.isin(label_primary_aud_cortex[0].vertices, vertices_in_aud_cortex))[0]
print(f"Number of dipoles is {'lower/equal' if n_dipoles <= vertices_in_aud_cortex.shape[0] else 'higher'} than the total vertices in the label border")
label_dipole = []
saving_vertices=[]
for one_dipole in np.random.choice(indices_in_label,n_dipoles, replace=False): # Creating a label object for
# One dipole source # each vertices, that can be passed down to simulator
saving_vertices .append( label_primary_aud_cortex[0].vertices[one_dipole])
label_dipole.append(mne.label.select_sources(
subject,
label_primary_aud_cortex[0],
location=one_dipole,
extent=extent,
subjects_dir=subjects_dir,
random_state=rng,
))
# ## Defining wavelet function
#
# %%
## function to make simple N100 like wavelet with damped sine
def simple_evoked(times):
base_amplitude = 50e-9
min_amplitude = 4e-8
frequency = 7
peak_time = 0.1
std_dev = 0.003 * abs(np.random.randn(1))
amplitude = max(min_amplitude, base_amplitude * abs(np.random.randn(1)))
phase_shift = -np.pi / 2 - 2 * np.pi * frequency * peak_time
sine_wave = np.sin(2 * np.pi * frequency * times + phase_shift)
gaussian_envelope = np.exp(-((times - peak_time ) ** 2) / std_dev)
signal = amplitude * sine_wave * gaussian_envelope
return signal
# ## Defining events and signal associated with each events , here one source will only have one kind of wave repeated at each trial
# %%
signal_matrix = np.zeros((n_dipoles,times.shape[0]))
trials_no = 50
n_events = trials_no * n_dipoles
events = np.zeros((n_events, 3), int)
events[:, 0] = 100+400 * np.arange(n_events)
for i in range(n_dipoles):
signal_matrix[int(i),:]=simple_evoked(times)
events[trials_no *(i):trials_no *(i+1), 2] = i+1
np.random.shuffle(events[:, 2]) # nonsequential event timings
# ## Adding data to sourcesimulator object
# %%
source_simulator_dipole = mne.simulation.SourceSimulator(src, tstep=tstep)
for i in range(n_dipoles):
source_simulator_dipole.add_data(label_dipole[i], signal_matrix[i,:], events[events[:,2]==i+1])
# %%
raw_sim = mne.simulation.simulate_raw(raw.info, source_simulator_dipole, forward=fwd)
raw_dipole = raw_sim.pick(picks=["meg", "stim"])
# cov = mne.make_ad_hoc_cov(raw_sim.info)
cov = noise_cov # I am loading noise cov data from sample provided by mne
mne.simulation.add_noise(
raw_dipole, cov, random_state=rng) #iir_filter=[0.2, -0.2, 0.04]
# ## Saving evoked and epochs in a dict
# %%
epochs_dict = dict()
evoked_dict = dict()
for i in range(n_dipoles):
epochs_dict[f'epoch_{i+1}']= mne.Epochs(raw_dipole, events, i+1, tmin=-0.2, tmax=0.6, baseline=(None,0.0))
evoked_dict[f'epoch_{i+1}'] = epochs_dict[f'epoch_{i+1}'].average().pick("mag").filter(l_freq=None, h_freq=40.0, picks=None, filter_length='auto',verbose=False)
# # Inverse modelling
# %%
snr = 3.0
inv_methods_list = ['sLORETA',"MNE","eLORETA","dSPM"]
lambda2 = 1.0 / snr**2
inv_meth_dict = dict()
for inv_method in inv_methods_list: #Probably use another way to store data later if there are 100s of sources, rather than dict of dict
stc_est_dipole_dict=dict() # Also change i and j in all loops
peak_dict = dict()
for i,j in evoked_dict.items():
inverse_operator = make_inverse_operator(j.info, fwd, cov, loose="auto", depth=0.8, fixed=True)
# inverse_operator = make_inverse_operator(j.info, fwd, cov, loose="auto", depth=None, fixed="auto")
stc_est_dipole_dict[i] = apply_inverse(j, inverse_operator, lambda2, inv_method, pick_ori=None, method_params={"max_iter":100})
peak_dict[i] = stc_est_dipole_dict[i].get_peak(hemi="lh")[0]
inv_meth_dict[f"{inv_method}"] = {"peak_coord":peak_dict,"stc_est_dict":stc_est_dipole_dict}
# ## Visualisation of peak coordinates wrt simulation vertices
# %%
for inv_method in inv_methods_list:
peak_dict = inv_meth_dict[inv_method]["peak_coord"]
brain_kwargs = dict(alpha=0.2, background="white", cortex="low_contrast",title=f"{inv_method}Vertices in Aud Cortex Present in SOurce Space")
correct=(np.array(list(peak_dict.values())) ==np.array(saving_vertices))
brain = mne.viz.Brain("sample", subjects_dir=subjects_dir, **brain_kwargs)
brain.add_foci(np.array(saving_vertices),hemi="lh",coords_as_verts=True,color="blue",scale_factor=0.2,alpha=0.3)
brain.add_foci(
np.array(list(peak_dict.values())),
coords_as_verts=True,
hemi="lh",
color="red",
scale_factor=0.4,
alpha=0.1,
)
brain.add_label(
label_primary_aud_cortex[0], hemi="lh", color="green", subdir="/home/hp/mne_data/MNE-sample-data/subjects/fsaverage/label/lh.aparc.a2009s",borders=True)
brain.save_image(f"images_to_mne_forum/{inv_method}_sct_peak_original_coord_peak.png")
#brain.close()
# %%
brain.save_image(f"images_to_mne_forum/zoomed_{inv_method}_sct_peak_original_coord_peak.png")