Testing BEM created EEG lead fields

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:

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

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.