Virtual channel analysis of epilepsy MEG data

The data for this tutorial can be downloaded from our ftp server

Male, age 9. Right parietal Glioma with parietal extended lesionectomy. Corticography also showed interictal discharges in the frontal lobe, though seizures were of parietal origin. Following the MEG, was operated and is now seizure free and off medication.

MEG data was recorded at Aston Brain Centre (ABC) using both a 275-channel CTF system and using an Elekta 306-channel system. This case report and the data are kindly provided by Professor Stefano Seri. The data has been clinically analysed by the staff of ABC using the software accompanying the MEG system. The FieldTrip analysis demonstrated here is only for educational purposes.

The original MRI that is provided for this patient has been partially processed in the CTF software and is stored in CTF *.mri format. This MRI is not shared for privacy reasons. Nevertheless, here we will show how it was processed in FieldTrip.

mri_orig = ft_read_mri('case1.mri');

The dataset also includes a Polhemus recording of the head surface, which can be used to coregister the MRI to the CTF system.

headshape = ft_read_headshape('case1.pos');
headshape = ft_convert_units(headshape, 'mm');

Check the coregistration of the Polhemus headshape and the anatomical MRI. In the CTF coordinate system the x-axis should be pointing to the nose and the y-axis to the left ear.

ft_determine_coordsys(mri_orig, 'interactive', 'no')
ft_plot_headshape(headshape);

In this case the coregistration is already nearly perfect.

Usually you would start with an anatomical MRI (e.g. stored as a stack of DICOM files) that is not yet coregistered with the MEG. We will improve the coregistration a bit, using the same procedure you could use to coregister from scratch.

The coregistration procedure starts with a coarse manual coregistration, followed by an automatic fine coregistration in which the skin surface from the MRI is fitted to the headshape points.

cfg = [];
cfg.method = 'headshape';
cfg.headshape.interactive = 'yes';
cfg.headshape.icp = 'yes';
cfg.headshape.headshape = headshape;
cfg.coordsys = 'ctf';
mri_realigned = ft_volumerealign(cfg, mri_orig);

The headshape not only covers the scalp, but also the face and nose. Hence the coregistration needs to be done prior to defacing from the anatomical MRI. The translate, rotate and scale parameters specified here were determined experimentally in the graphical user interface of the ft_defacevolume function.

cfg = [];
cfg.translate = [100 0 -60];
cfg.rotate = [-4 30 0];
cfg.scale = [75 150 120];
mri_defaced = ft_defacevolume(cfg, mri_realigned);

For convenience in later plotting, we reslice the MRI so that the axes of the volume are aligned with the axes of the coordinate system.

cfg = [];
mri_resliced = ft_volumereslice(cfg, mri_defaced);

Finally we should do a visual inspection of the realigned, defaced and replaced MRI.

cfg = [];
ft_sourceplot(cfg, mri_resliced)

Note that the patients head is tilted to the right. Apparently the anatomical landmarks at the left and right ear were not clicked symmetrically with the Polhemus. This is not a problem for further processing, as long as we remember that results are expressed in head coordinates relative to the anatomical landmark of this specific recording.

dataset = 'case1.ds';

cfg = [];
cfg.dataset   = dataset;
cfg.hpfilter  = 'yes';
cfg.hpfreq    = 10;
cfg.lpfilter  = 'yes';
cfg.lpfreq    = 70;
cfg.channel   = 'MEG';
data = ft_preprocessing(cfg);

%% visualize the preprocessed data

cfg = [];
cfg.viewmode = 'vertical';
cfg.channel  = 'MEG';
cfg.layout   = 'CTF275.lay';
cfg.event    = ft_read_event(dataset);
ft_databrowser(cfg, data);

%% compute channel-level kurtosis

datak = [];
datak.label    = data.label;
datak.dimord   = 'chan';
datak.kurtosis = kurtosis(data.trial{1}')';

cfg = [];
cfg.comment = 'computed channel-level kurtosis';
datak = ft_annotate(cfg, datak);

cfg = [];
cfg.layout    = 'CTF275.lay';
cfg.parameter = 'kurtosis';
ft_topoplotER(cfg, datak);

% caxis([0 40])

% we will use the defaced MRI, which has been realigned with the CTF system and resliced

mri = ft_read_mri('mri_defaced.mat');

%% segment the brain compartment from the anatomical MRI and make the volume conduction model

cfg = [];
cfg.tissue = 'brain';
seg = ft_volumesegment(cfg, mri);

% save seg seg

cfg = [];
cfg.tissue = 'brain';
brain = ft_prepare_mesh(cfg, seg);

cfg = [];
cfg.method = 'singleshell';
headmodel = ft_prepare_headmodel(cfg, brain);

% save headmodel headmodel

cfg = [];
cfg.grid.resolution = 7;
cfg.grid.unit = 'mm';
cfg.vol = headmodel;
cfg.grad = data.grad; % this being needed here is a silly historical artifact
sourcemodel = ft_prepare_sourcemodel(cfg);

%%
cfg = [];
cfg.xrange = [min(sourcemodel.pos(:,1))-30 max(sourcemodel.pos(:,1))+30];
cfg.yrange = [min(sourcemodel.pos(:,2))-30 max(sourcemodel.pos(:,2))+30];
cfg.zrange = [min(sourcemodel.pos(:,3))-30 max(sourcemodel.pos(:,3))+30];
mri_resliced = ft_volumereslice(cfg, mri_defaced);

% save mri_resliced mri_resliced

figure
ft_plot_vol(headmodel, 'unit', 'mm');
ft_plot_sens(data.grad, 'unit', 'mm', 'coildiameter', 10);
ft_plot_mesh(sourcemodel.pos);
ft_plot_ortho(mri_resliced.anatomy, 'transform', mri_resliced, 'style', 'intersect');

%% compute data covariance for source reconstruction

cfg = [];
cfg.channel = 'MEG';
cfg.covariance = 'yes';
timelock = ft_timelockanalysis(cfg, data);

%% this is not required, but speeds up repeated source reconstructions

cfg = [];
cfg.channel = 'MEG';
cfg.vol  = headmodel;
cfg.grid = sourcemodel;
cfg.normalize = 'yes';
sourcemodel = ft_prepare_leadfield(cfg, timelock);

% save sourcemodel sourcemodel

cfg = [];
cfg.vol  = headmodel;
cfg.grid = sourcemodel;
cfg.method = 'lcmv';
cfg.lcmv.projectmom = 'yes';
sourcep = ft_sourceanalysis(cfg, timelock);

sourcep.kurtosis = nan(size(sourcep.pos,1),1);
sourcep.kurtosisdimord = 'pos';
sel = find(sourcep.inside(:));
for i=1:length(sel)
  disp(i);
  sourcep.kurtosis(sel(i)) = kurtosis(sourcep.avg.mom{sel(i)});
end

cfg = [];
cfg.comment = 'computed source-level kurtosis';
sourcep = ft_annotate(cfg, sourcep);

%% explore the results

cfg = [];
cfg.parameter = 'kurtosis';
sourcepi = ft_sourceinterpolate(cfg, sourcep, mrir);

%%

% cfg = [];
% cfg.funparameter = 'mom';
% ft_sourceplot(cfg, sourcep);

%%

cfg = [];
cfg.funparameter = 'kurtosis';
ft_sourceplot(cfg, sourcep);

%%

cfg = [];
cfg.funparameter = 'kurtosis';
ft_sourceplot(cfg, sourcepi);

%% find peaks and plot timecourse

ispeak = findpeaksn(sourcep.kurtosis);
j = find(ispeak(:));
[m, i] = sort(-sourcep.kurtosis(j));
peaks = j(i);
disp(sourcep.pos(peaks(1:20),:));

% peak 1 is left frontal
% peak 4 is right frontal

s1 = sourcep.avg.mom{peaks(1)};
s4 = sourcep.avg.mom{peaks(4)};
shift = 1.1*(max(s1) - min(s4));

figure;
plot(sourcep.time, s1, 'r');
hold on
plot(sourcep.time, s4 + shift, 'm');
legend({'area 1', 'area 4'})

datas = [];
datas.time = {sourcep.time};
datas.trial = {[s1; s4]};
datas.label = {'area 1', 'area 4'}';

cfg = [];
cfg.comment = 'constructed virtual-channel raw data structure';
datas = ft_annotate(cfg, datas);

%%

cfg = [];
cfg.funparameter = 'kurtosis';
cfg.location = sourcep.pos(peaks(1),:);
ft_sourceplot(cfg, sourcepi);
cfg.location = sourcep.pos(peaks(4),:);
ft_sourceplot(cfg, sourcepi);

%% find maxima in virtual channel time series based on 6 SDs

sd1 = std(s1);
sd4 = std(s4);

tr = abs(s1)>6*sd1 & abs(s4)>6*sd4;
% tr = conv(double(tr), ones(1,60), 'same');
tr = tr>0;
begsample = find(diff([tr 0]==1));
endsample = find(diff([tr 0]==1));

%%

cfg = [];
datac = ft_appenddata(cfg, data, datas);

cfg = [];
cfg.artfctdef.interictal.artifact = [begsample(:) endsample(:)];
cfg.viewmode = 'vertical';
cfg.channel = 'MEG';
cfg.layout = 'CTF275.lay';
% cfg.event = ft_read_event(dataset);
ft_databrowser(cfg, datac);

%% show the provenance of the analysis pipeline

cfg = [];
cfg.filename = 'sourcepi';
cfg.filetype = 'html';
ft_analysispipeline(cfg, sourcepi);

!open sourcepi.html

Female, age 14. Epilepsy. Referral for MEG because EEG did not allow laterlisation or localisation of discharges, though clinically they appeared to come from the left hemisphere. Functional neuroimaging in the form of a PET scan showed a right area of hypo metabolism. Surgical follow-up information about this patient is not available.

MEG data was recorded at Aston Brain Centre (ABC) using both a 275-channel CTF system and using an Elekta 306-channel system. This case report and the data are kindly provided by Professor Stefano Seri. The data has been clinically analysed by the staff of ABC using the software accompanying the MEG system. The FieldTrip analysis demonstrated here is only for educational purposes.