When I computed the forward solution of my fitted dipole and compared the predicted evoked with the measured evoked, I saw a huge difference in topography. What could be the reason for it? Is it faulty head model or something else?
the left plot (labelled āmeasured fieldā) and the right plot (labelled ādifferenceā) look identical to me. Are you sure youāre actually plotting the difference correctly? Looking at the left and middle plots, Iām not sure Iād call that a āhugeā difference, assuming the colormaps are on the same scale.
Hi @ drammock,
I used the code in the MNE tutorial. As youāve mentioned, I am not sure why the measured field and the one labelled difference look nearly the same.
best_idx = np.argmax(dip.gof)
best_time = dip.times[best_idx]
print('Highest GOF %0.1f%% at t=%0.1f ms with confidence volume %0.1f cm^3'
% (dip.gof[best_idx], best_time * 1000,
dip.conf['vol'][best_idx] * 100 ** 3))
# remember to create a subplot for the colorbar
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=[10., 3.4],
gridspec_kw=dict(width_ratios=[1, 1, 1, 0.1],
top=0.85))
vmin, vmax = -400, 400 # make sure each plot has same colour range
# first plot the topography at the time of the best fitting (single) dipole
plot_params = dict(times=best_time, ch_type='mag', outlines='skirt',
colorbar=False, time_unit='s')
evoked.plot_topomap(time_format='Measured field', axes=axes[0], **plot_params)
# compare this to the predicted field
pred_evoked.plot_topomap(time_format='Predicted field', axes=axes[1],
**plot_params)
# Subtract predicted from measured data (apply equal weights)
diff = combine_evoked([evoked, pred_evoked], weights=[1, -1])
plot_params['colorbar'] = True
diff.plot_topomap(time_format='Difference', axes=axes[2:], **plot_params)
fig.suptitle('Comparison of measured and predicted fields '
'at {:.0f} ms'.format(best_time * 1000.), fontsize=16)
fig.tight_layout()
This is what my BEM looks like. It is not perfect, so I was wondering if it could be the reason for the poor Goodness of fit (only 5%)
@agramfort , I am sharing the files with you here. If possible, would you please take a look at what I could have possibly done wrong? The GOF for my fitted dipole is too low
import matplotlib
import pathlib
import mne
import mne_bids
matplotlib.use('Qt5Agg')
import numpy as np
cd E:\
fname='phantom_rebuilt_MEG_4232022_L6.fif'
raw=mne.io.read_raw_fif(fname)
raw.load_data()
raw.filter(l_freq=1,h_freq=70)
picks = mne.pick_types(raw.info, meg= True, eog=False,
stim=False, exclude='bads')
raw.notch_filter(np.arange(60, 241, 60), picks=picks, filter_length='auto',
phase='zero')
event=mne.read_events('E:\L6M.eve')
for i in range(1286):
event[i][1]=0
for i in range(1286):
event[i][2]=1 # making all spikes as event 1
for i in range(1285):
if event[i+1][0]- event[i][0]< 10:
event[i+1][0]=event[i][0]
event[i+1][2]=event[i][2] # if two spikes are really close to each other, making them same event
event_dict = {'1':1}
tmin, tmax = -.2, 0.3
epochs = mne.Epochs(raw, events=event, event_id=event_dict,tmin=-0.2, tmax=0.3,event_repeated='drop', baseline=None)
epochs.plot()
epochs.save('revised.fif',overwrite=True)
subjects_dir = 'E:\myMRI-20220510T004802Z-001\myMRI'
mne.viz.plot_bem(subject='E:\myMRI-20220510T004802Z-001\myMRI\_my_subject', subjects_dir=subjects_dir,
orientation='coronal')
trans='C.fif'
subject = 'E:\myMRI-20220510T004802Z-001\myMRI\_my_subject'
src = mne.setup_source_space(subject='E:\myMRI-20220510T004802Z-001\myMRI\_my_subject',
spacing='oct6',
subjects_dir=subjects_dir,
add_dist=False)
conductivity = (0.3, 0.006, 0.3)
model = mne.make_bem_model(subject='E:\myMRI-20220510T004802Z-001\myMRI\_my_subject', ico=4,
conductivity=conductivity,
subjects_dir=subjects_dir)
bem_sol = mne.make_bem_solution(model)
bem_sol
bem_fname = 'E:\myMRI-20220510T004802Z-001\myMRI\sampleBEM.fif'
mne.bem.write_bem_solution(bem_fname, bem_sol,overwrite=True)
cov=mne.compute_covariance(epochs,tmin=-0.2,tmax=-0.05)
evoked=epochs.average()
evoked.plot()
evoked_full = evoked.copy()
evoked.crop(0.00, 0.01)
dip = mne.fit_dipole(evoked, cov, bem=bem_sol, trans=trans)[0]
dip.plot_locations(trans=trans, subject=subject, subjects_dir=subjects_dir, mode='orthoview',coord_frame='mri',idx='gof')
Here is the google drive link to all the files
https://drive.google.com/drive/folders/1rZL6Zj-nZfzUlLNJKF2ThfTp0oY_93Et?usp=sharing
@Bamford - Not entirely sure of the discrepancy between the diff and original. I looked at your data and noticed a few items that might help. Your spikes (Iām assuming epilepsy) are defined with a maximum at time t=0. Therefore your baseline shouldnāt be -0.2-0. If you change it to (-0.2,-0.1), you will get better results. Also at the end of your data you have a lot of noise - see image below. This needs to be removed from the data or it will dominate the signal. Also - I am not sure why, but the EEG doesnāt seem to be baseline correcting properly - I donāt know if this is a reference issue (??) - but the data near the end is floating off of baseline.
epo = mne.read_epoch('C.fif', preload=True)
epo.apply_baseline(baseline=(-0.2, -0.1))
epo.plot()
Produces this image - these are the last trials and the baseline does not appear to be set (??).
If I just take the first 100 epochs - you get a pretty strong activity vs baseline. The data also appears to be more dipolar (I work with MEG - but I think this is appropriate for an EEG dipole). Fitting this data will give you a better dipole fit / GOF.
epo2 = epo[0:100]
epo2.filter(1.0, None)
epo2.plot()
ave=epo2.average()
One other major component that is going to affect your data: your 3 layer BEM does not resolve the inner and outer-skull boundary well. In EEG this is pretty important since the skull is ~40X the resistance of the other layers. With a typical T1w MRI, it is not uncommon to have this issue.
And the final problem with this type of analysis is that the underlying neurophysiology may not be represented by a single dipole. With these large discharges, it is possible that there is a large synchronous network firing together. With that amount of tissue active - you might not be able to represent all of that activity with a single dipole (typically these results will send the dipole deeper). Its possible that distributed modelling (MNE, Loretaā¦) would be preferred.
Those are some initial thoughts - but I canāt comment too much on the dipole simulation.
-Jeff
3 Likes
@Bamford - I just looked over what you had posted again. It appears that you are doing phantom MEG (āphantom_rebuilt_MEGā) and in your picks, you have meg=True. But in your data it is listed as EEG
Out[74]:
<Info | 13 non-empty values
bads: 21 items (E67, E68, E73, E74, E92, E83, E82, E91, E94, E95, ...)
ch_names: E1, E2, E3, E4, E5, E6, E7, E8, E9, E10, E11, E12, E13, E14, ...
chs: 256 EEG, 1 misc
custom_ref_applied: True
device_info: 1 item (dict)
dig: 261 items (3 Cardinal, 257 EEG, 1 Extra)
Not sure if this is the expected data. And could result in the poor dipole fit if the sensor type is different from the actual data.
-Jeff
3 Likes
hI @jstout211 Thank you for taking your time out to go through the files. I seem to have attached a wrong epoch file ( the one with EEG here). I have updated the google drive with the new epoch file [ newmeg.fif].
I was able to get a GOF of 97% with the way suggested by you for EEG file, however the GOF for MEG files is still low. Could you give it a check?
Thanks
@Bamford - I looked at your MEG data. The sensor quality is a little noisy and inconsistent between channels. I marked about 10 as bad before making the epochs. I donāt know how much this really contributes to your issue - but just wanted to point it out.
I averaged all of the epochs and got the same topography map as you posted. It looks like your sensor data go into a steady state rythmn where the filter rolloff bleeds into the next trial (you can see the persistent activity between the two trials below). I thought that using the initial stimuli may get around this problem, since it wonāt be affected by prior trials. I set the baseline as 5 seconds prior to the initial event and everything looks relatively flat until the dipole stimuli. Unfortunately you still get the same topography.
The topography maps (at least the mag data) does not appear to show great phantom dipoles. It should appear to look like a butterfly pattern with a positive and negative lobe (similar to your dipole map), but the topography still looks a little noisy / non-dipolar for having such a high signal to noise level. I think your gradiometer data looks more appropriate, but there is some distant noise there as well.
Another comment on your data - you have an MRI that you are using to fit the dipole - but the MEG data does not have any neurophysiological activity. I would assume you are better off using a spherical head model that represents your phantom head. I dontā really know how to do this in MNE.
I also am not sure what is causing the topography to be distorted, which would lead to a poor diple fit. There are many factors - possibly channel tuning (??). The magnetometers are going to be more sensitive to artifact compared with the gradiometers. Just as a quick test, I would select just the gradiometers (raw.pick_types(meg=āgradā)) and redo the forward model and the dipole fit using the MRI. See if this helps. It could be that the mags are just driven outside of their linearity range.
Did you perform SSS / tSSS on the data? Its possible that this might be able to help. Are you running MEG and EEG concurrently on this phantom? Or are the EEG leads still attached to the phantom, but not to the electronics - they can act like antennae and distort the fields.
Sorry I couldnāt help any more with this - I am kind of stumped as to why you are getting this kind of topography.
-Jeff
2 Likes