[FieldTrip] Head model for EEG source localization - Possible bug in BEMCP?

Debora Desideri deb.desideri at gmail.com
Sat Mar 12 15:06:38 CET 2016


Dear all,
I am working on EEG source localization of a very simple fist clenching
task (right hand, 10 seconds of contraction alternated to 10 seconds of
relaxation to be used as contrast condition) using DICS beamformer
(frequency 10 Hz). I have the individial MRI and EEG electrode positions of
my participants, so I wrote a script to create the individual forward model
following the tutorials on the fieldtrip webpage and some useful
discussions on the mailing list. Below I report the full code used for the
forward model and the code used for the source analysis. Attached, you
find the  results of the source localization and
​some ​
figures of different steps in the computation of the forward model.

According to my task, I expect to identify a desynchronization on the left
motor cortex, but what I am getting are very weird results(see attachment
​s)​
.  I tryed to use the standard MRI in template\headmodel\standard_mri, and
I get similar results.

After a full day of self-debugging, I managed to get the result I expect
(attached) using the method "concentricspheres" instead of "bemcp" in the
computation of the head model, both with the MRI of my participants and
with the standard MRI.

I assume there must be something wrong with "bemcp", also becasue when I
use the standard bem headmodel (template\headmodel\standard_bem) created
with the method "dipoli", everything looks fine again.


Did anybody experienced a similar issue before? Is there anything wrong I
am doing in the pipeline?

Thanks a lot in advace for your support!
Best,
Debora


Fieldtrip version: 20150828
Matlab version: R2015a (64 bit)
OS: Windows 8 (64 bit)



%%%% --------------------- Code for  Forward Model ---------------- %%%


%% Read  MRI data (DICOM)
mri = ft_read_mri(experiment.mri); % experiment.mri contains the path to
individual MRIs

%% Convert anatomical MRI to MNI coordinate system (= RAS)

cfg = [];
cfg.method = 'interactive';
cfg.coordsys = 'spm'; % MNI coordinate system (=RAS)
cfg.viewmode = 'ortho';
mri_spm = ft_volumerealign(cfg,mri);

%% Re-slice the MRI volume so that all the 3 dimensions have the same
number of voxels

cfg = [];
cfg.resolution = 1; % 1 mm thick slices
cfg.dim = [256 256 256];

mri_spm = ft_volumereslice(cfg,mri_spm);

%% Visualize MRI data to check that everything is ok

cfg = [];
cfg.method = 'ortho';
ft_sourceplot(cfg,mri_spm)

%% Segment the anatomical data

cfg = [];
cfg.output = {'brain' 'skull' 'scalp'};
cfg.brainsmooth = 2;
cfg.brainthreshold = 0.25;
cfg.scalpthreshold = 0.025;
mri_seg    = ft_volumesegment(cfg, mri_spm);

%% Check segmentation brain
cfg = [];
cfg.method = 'ortho';
cfg.funparameter = 'brain';
ft_sourceplot(cfg,mri_seg)

%% Check segmentation skull
cfg = [];
cfg.method = 'ortho';
cfg.funparameter = 'skull';
ft_sourceplot(cfg,mri_seg)

%% Check segmentation scalp
cfg = [];
cfg.method = 'ortho';
cfg.funparameter = 'scalp';
ft_sourceplot(cfg,mri_seg)


%% Prepare mesh: creates surfaces at the boarders of the different
tissue-types
cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = {'brain' 'skull' 'scalp'};
cfg.numvertices = [6000 4000 2000];
bnd = ft_prepare_mesh(cfg, mri_seg);

%% Visualize surfaces
figure
ft_plot_mesh(bnd(1),'facecolor','none') % brain
figure
ft_plot_mesh(bnd(2),'facecolor','none') % skull
figure
ft_plot_mesh(bnd(3),'facecolor','none') % scalp

%% Visualize full mesh
figure
ft_plot_mesh(bnd(3), 'facecolor',[0.2 0.2 0.2], 'facealpha', 0.3,
'edgecolor', [1 1 1], 'edgealpha', 0.05);
hold on;
ft_plot_mesh(bnd(2),'edgecolor','none','facealpha',0.4);
hold on;
ft_plot_mesh(bnd(1),'edgecolor','none','facecolor',[0.4 0.6 0.4]);


%% Prepare head model (volume conduction model)
cfg = [];
cfg.method = 'bemcp';
headmodel = ft_prepare_headmodel(cfg,bnd);


%% % --------------------- Align the EEG electrodes!!!-------------------%%%

%% (1) Read electrode position from Localizer file

elec = ft_read_sens(experiment.electrodes); % experiment.electrodes
contains the path to the individual electrode positions



%% (2) Check if alignment is necessary

figure
%plot scalp
ft_plot_mesh(headmodel.bnd(3),'edgecolor','none','facealpha',0.8,'facecolor',[0.6
0.6 0.8])
hold on
% plot electrodes
ft_plot_sens(elec,'style','sk')


%% (3) If the electrodes are not aligned to the head surface, interactively
align them

cfg = [];
cfg.method = 'interactive';
cfg.elec = elec;

cfg.headshape = headmodel.bnd(3); % scalp
elec_aligned = ft_electroderealign(cfg);


%%  Load hdr and order elec_aligned accordign to original order of signal
acquisition
%%% This is necessary because the eeg elctrod eposition shave been acquired
%%% in a different order with respect to theorder of the signal recording

load([experiment.path '\Data'],'hdr') % hdr contains some information about
the original order of the channels and so on...

 [~, original_order] = ismember(hdr.elec.label,elec_aligned.label);


elec_aligned.chanpos = elec_aligned.chanpos(original_order,:);
elec_aligned.elecpos = elec_aligned.elecpos(original_order,:);
elec_aligned.label = elec_aligned.label(original_order);

%% Check if the alignment is ok

figure
%plot scalp
ft_plot_mesh(headmodel.bnd(3),'edgecolor','none','facealpha',0.8,'facecolor',[0.6
0.6 0.8])
hold on
% plot electrodes
ft_plot_sens(elec_aligned,'style','sk')



%% -------------------------------- Source model -------------------------%%

%% Make the individual subject's grid


cfg = [];
cfg.grid.warpmni = 'yes';
cfg.resolution = 6; % resolution of the template grid in mm
cfg.grid.nonlinear = 'yes'; % use non linear normalization
cfg.mri = mri_spm; % use subject's mri in mni coordinates
cfg.grid.unit = 'mm';
sourcemodel = ft_prepare_sourcemodel(cfg);


%% Plot the single subject'e head model and grid positions

figure
hold on
ft_plot_vol(headmodel,'edgecolor','none','facealpha',0.3);
ft_plot_mesh(sourcemodel.pos(sourcemodel.inside,:));

%% Double check that source model is aligned to the segmented volume and to
the electrodes
figure
hold on
plot3(sourcemodel.pos(sourcemodel.inside,1),sourcemodel.pos(sourcemodel.inside,2),sourcemodel.pos(sourcemodel.inside,3),'.');
ft_plot_vol(headmodel,'facecolor','brain','edgecolor','none','facealpha',0.3);
ft_plot_sens(elec_aligned,'label','label','style','y*','Markersize',10);


%% ----------------------  FORWARD SOLUTION ------------------------ %%


%% Prepare leadfield

cfg = [];
cfg.elec = elec_aligned;
cfg.grid = sourcemodel;
cfg.headmodel = headmodel;
leadfield = ft_prepare_leadfield(cfg);


%% Test sourc elocalization



 load([experiment.path '\DATA'], 'preexperiment_freqsource')

 % preexperiment_freqsource contains the output of ft_freqanalysis

        cfg                    = [];
        cfg.elec               = elec_aligned;
        cfg.method             = 'dics';
        cfg.frequency          = 10;
        cfg.grid               = leadfield;
        cfg.headmodel          = headmodel;
        cfg.dics.keepfilter    = 'yes';
        cfg.dics.lambda       = '10%';
        src                    = ft_sourceanalysis(cfg,
preexperiment_freqsource);
        filter                 = src.avg.filter;


        markers = [224 192]; % markers are some number in the trialinfo,
224 is related to fist clenching, 192 to relaxation
        for mm = 1:length(markers)

           cfg = [];
           cfg.trials = find(preexperiment_freqsource.trialinfo ==
markers(mm));
           data = ft_selectdata(cfg,preexperiment_freqsource);

            %%% DICS
            cfg                    = [];
            cfg.elec               = elec_aligned;
            cfg.method             = 'dics';
            cfg.frequency          = 10;
            cfg.grid               = leadfield;
            cfg.headmodel          = headmodel;
            cfg.dics.projectnoise  = 'yes';
            cfg.dics.lambda       = '10%';
            cfg.grid.filter        = filter;
            Source{mm}            = ft_sourceanalysis(cfg, data);

        end


 %% Interpolate results to anatomy

    source = Source{1};
    source1 = Source{1}; %clenching
    source2 =  Source{2}; % release
    source.avg.cont = (source2.avg.pow./source1.avg.pow)-1;

    cfg = [];
    cfg.parameter         = 'cont';
    source               = ft_sourceinterpolate(cfg, source, mri_spm);



 %% Plot surface

    cfg = [];
    cfg.method            = 'surface';
    cfg.funparameter      = 'cont';
    cfg.funcolormap = 'jet';
    ft_sourceplot(cfg, source)


--------------------------------------------------------------------

Debora Desideri, PhD Student
BNP Lab - Neurology Department
University Hospital Tübingen
Eberhard Karls University Tübingen
Hoppe-Seyler Str. 3
D-72076 Tübingen
Tel: +49 7071/29 80478
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 5_CONCENTRICSPHERES_Results_Source.PNG
Type: image/png
Size: 126346 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0010.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 0_IndividualMRI_spm.PNG
Type: image/png
Size: 80674 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0011.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 1_Segmentaiton_BrainSkullScalp.PNG
Type: image/png
Size: 41612 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0012.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 2a_Output ft_plot_mesh_BrainSkullScalp.PNG
Type: image/png
Size: 66519 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0013.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 2b_Output ft_plot_mesh_allSurfaces.PNG
Type: image/png
Size: 40960 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0014.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 3_BEMCP_SourceModel.PNG
Type: image/png
Size: 19959 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0015.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 3_CONCENTRICSPHERES_SourceModel.PNG
Type: image/png
Size: 19714 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0016.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 4_BEMCP_Alignment_HeadModel_SourceModel_Electrodes.PNG
Type: image/png
Size: 58428 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0017.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 4_CONCENTRICSPHERES_Alignment_HeadModel_SourceModel_Electrodes.PNG
Type: image/png
Size: 67775 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0018.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 5_BEMCP_Results_Source.PNG
Type: image/png
Size: 143105 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160312/a9845405/attachment-0019.png>


More information about the fieldtrip mailing list