Testing BEM created EEG lead fields
Introduction
The following section compares lead fields calculated by means of BEM implementation and theoretical model, by correlating their values for different densities of the head model vertices. It is expected that solutions for a single shell EEG realistic head model converge to the theoretical one with an increasing number of points in the mesh.
For the simplest case, the BEM and the theoretical solutions for EEG lead field can be compared directly, because theoretical EEG forward solution for a spherical shell is known. Let's calculate the electric field on top of a sphere, generated by one dipole source with different positions along the z-axis:
% Create a sphere of radius '100' conductivity '1' and center '[0,0,0] (mm)' within a volume structure vol = []; vol.r = 100; vol.c = 1; vol.o = [0 0 0]; % Create a set of 42 electrodes on the outer surface [pnt, tri] = icosahedron42; sens.pnt = vol.r * pnt; sens.label = {}; nsens = size(sens.pnt,1); for ii = 1:nsens sens.label{ii} = sprintf('vertex%03d', ii); end % Define positions of dipoles (along z axis) zp = linspace(0,100,50)'; pos = [zeros(size(zp)) zeros(size(zp)) zp]; % Define the corresponding spatial grid cfg = []; cfg.inwardshift = 5; cfg.grid.pos = pos; gridp = ft_prepare_sourcemodel(cfg, vol, sens);
BEM model, single sphere:
Building the geometrical head model with BEM
The following code generates different BEM (Dipoli) system matrixes, describing head properties which are depending on the increasing number of triangles used for each mesh (from 50 - the lowest - to 2000 - the highest -)
This implements Oostendorp and van Oosterom, 1989, eq. 2.3, 2.4
vertic = [50 80 100 200 2000]; for ll=1:length(vertic) volbem = vol; cfg = []; cfg.method = 'dipoli'; cfg.conductivity = 1; cfg.numvertices = vertic(ll); cfg.isolatedsource = false; [volbem] = ft_prepare_bemmodel(cfg, volbem); vols{ll} = volbem; end
Simulation and results
The following steps generate the lead fields with 1) the BEM model integrated in the forward equations, 2) the theoretical forward solution for a dipole in spherical conductor. Sources' positions are defined in the variable gridp, head properties are defined in vols{ll} and sensors' positions in sens.
% calculate BEM lead fields: for ll=1:length(vertic) cfg = []; cfg.grid = gridp; cfg.vol = vols{ll}; cfg.elec = sens; gridbem = ft_prepare_leadfield(cfg); grids{ll} = gridbem; end % calculate theoretical singleshell forward solution cfg = []; cfg.grid = gridp; cfg.vol = vol; cfg.elec = sens; gridSph = ft_prepare_leadfield(cfg);
Result:
The pinky arrow describes the correlation curves of meshes with increasing number of triangles. The last mesh (2000 vertices) has a flat correlation curve at value y=1.