Lead field calculation help

Hello,

I’m new to BrainStorm and am trying to use it for something (hopefully) straightforward. I just want to calculate a lead field for an anatomical MRI, in the coordinate system defined by the anatomical image if possible. I have scalp/cortex tesselations from BrainSuite. I don’t want to import/analyze my EEG data in BrainStorm, create a database, etc. (my data is not in a convenient format so I have my own analysis software). I have EEG electrode coordinates (xyz, text file), but they are measured on a styrofoam head so I don’t want to warp the MRI to match them but would rather warp the electrode locations to match the MRI. BrainStorm doesn’t need to do that for me, unless it can (the tutorial seems to indicate it can only take MRI into electrode space and not vice versa). It seems if I only want the lead field, I shouldn’t need to input any EEG data of any kind?

Because of my needs, I’m guessing I want to use the real “nuts and bolts” numerical functions of BrainStorm, avoiding all the GUI and data-management stuff. (Hopefully I could just call them from my own MATLAB code as library functions). I don’t know where to start on that and would need help understanding the output format, since I want to use the “answer” (lead field matrix and list of solution points) in further calculations. I am very new to EEG but not MATLAB, so writing software to nudge my data into the appropriate form and vice versa isn’t a problem, though I would like to avoid filling up large BrainStorm-specific structs as much as possible (since I don’t plan to use 90% of BrainStorm’s functionality).

Any help/guidance you can provide with this would be greatly appreciated. Thank you!

Kevin

Hi Kevin,

Sure you may use BST’s basic scripts to do your thing. Regarding computation of EEG lead-fields, the very basic script is eeg_sph.m which computes forward models using a spherical geometry for head tissues. See the script header for help. When it says ‘Channel structure’, it means you need to arrange your sensor information in the specific BST structure format which is described when you type help_data_channel at the Matlab prompt when BST is in your path.

I think it’s a good way to start. Let me know how things are moving.

Hello Sylvain,

Thank you for pointing me towards eeg_sph. Getting my electrode coords. into a channel structure shouldn’t be a problem. However, the Param structure needs some anatomical info (radii and such) that I need to obtain from my BrainSuite parameterizations. Is there a function that accepts those tesselations/surfaces and returns those radii?

Thanks for your continued help,
Kevin

Hi Kevin,

for a simple approach, you may fit a sphere to your electrode locations to begin with. Hence, you’ll just need to replicate the same values for Radii and Centre in the Param structure.

We use the following to adjust a sphere to a set of points (vertices or electrode locations) this is an exerpt from bst_headmodeler which is the core routine for computing the various foward models in BrainStorm:


 % Run parameters fit ----------------------------------------------------
mass = mean(Vertices); % center of mass of the scalp vertex locations   
R0 = mean(norlig(Vertices - ones(nscalp,1)*mass)); % Average distance between the center of mass and the scalp points
vec0 = [mass,R0];
                
[minn,brp] = fminsearch('dist_sph',vec0,[],Vertices);
OPTIONS.HeadCenter = minn(1:end-1); 
if isempty(OPTIONS.Radii)
  %OPTIONS.Radii = minn(end);
  OPTIONS.Radii = minn(end)*[1/1.14 1/1.08 1]; % 3-shell model
end

I am attaching the dis_sph.m function.

Hope this helps.

Sylvain.

Hello Sylvain,

Thank you for your continued help and the code snippets. I think I may have what I need to get started (at least with a single sphere model), since I also found in another post a bit of code to make the BrainStorm channel structure. I will let you know how I’m doing.

Best,
Kevin

Hello Sylvain,

Now I’m ready to run eeg_sph.m but I’m getting an error. I’ve loaded my channel info into a 1 X 129 struct with the appropriate members (electrode 129, Cz, is the reference). The Param struct is also 1x129 - all 129 elements are the same, containing Radii, Conductivity, Center, and EEGType members, set accordingly I hope (I used a modification of your code snippet to fit the sphere to the electrode locations).

I’m trying to call eeg_sph with these structures and a 3X100 matrix of random numbers for L (scaled to be random locations inside the head; I just want to try to get some output), I get a “Matrix dimensions must agree” error here:

Rq = L’;
center = [Param.Center]’;
Re = Re - center;

Rq = Rq - repmat(center(1,:),size(Rq,1),1);

I’ve gone through line-by-line and it seems like there must be an error here, or a MATLAB inconsistency - I’m using R2007a. The problem is L’ is 100 X 3, and center is 387 X 1. The Re subtraction is fine, but the Rq line has inconsistent dimensions; if eeg_sph is trying to subtract the position of the center from the positions of all my sources, shouldn’t it be

Rq = Rq - repmat(center(1:3)’,size(L’,1),1)?

center(1,:slight_smile: only gets the x position of the center, not (x,y,z), so even if the line as written had the right dimensions, it seems it wouldn’t have the desired behavior. What am I missing here?

Thanks,
Kevin

Hi Again Sylvain,

I’m not sure if you saw my last post, but if I make the change I indicated to eeg_sph.m I can obtain a gain matrix of the expected dimensions. This must mean some of my input data is transposed or has unexpected dimensions. It seems the problem is with Params(i).Center? Mine is a row vector; should it be a column vector or?

I’ve almost got it (in fact maybe I do have it, given my modification of eeg_sph.m) and I really want to thank you for your continued help and clear answers.

Kevin

No problem. I think indeed you need to transpose the center coordinates.

Let me know whether you are facing any other issue.

Hi Sylvain,

It turned out almost everything in the Param struct needed to be a column vector - including Param.Center and Param.Loc. After that, eeg_sph runs fine. I now have another issue - for my purposes I would actually like a scalar lead field rather than a vector one. Is there any (sensible) way to obtain this via BrainStorm?

Thanks again,
Kevin

What do you mean by ‘scalar’ here? The lead fields computed by BrainStorm are vectors as they link all sources from the model to the surface sensors. Maybe you need just one source and one sensor ? This should not pose any problem in principle.

Hi Sylvain,

Nice to hear from you again! You can ignore my previous question (I was sort of hoping I could get an ELECTRA-type ‘scalar’ lead field, relating LFP to scalp potential, out of BrainStorm). I’ve got another one. I’d eventually like to be able to do what I’m doing with eeg_sph for a realistic head model. I can obtain cortical surface/scalp/etc. parameterizations from BrainSuite (already figured that out), and again I’d like to sidestep the BrainStorm databases and such and use the low level functions. Ideally I would be feeding in just sensor locations and parameterized surfaces and using the low-level BEM/FEM solvers in BrainStorm to do the numerics and give me the lead field. Is there a place you can point me to get me started on this?

As always, thanks for your help!
Kevin

Hi Kevin,

I am not sure if my comments will help you, but it seems you want to calculate a leadfield without any GUI. As I am doing it, I will post some code that might be usefull for you although I am using MEG data. Here is the code:

%---------------------------------------------------------%
% save scalp in brainstorm tess format %
%---------------------------------------------------------%

Vertices{1} = sfpreg2./100; %% take head shape points
sfpface=convhulln(sfpreg2’); %% create faces
Faces{1} = sfpface;
Comment{1} = ‘scalp’;
command = sprintf(‘save %s Vertices Faces Comment’,[outpath,vp{v},’_brst_scalp.mat’]); % save it for later use in bst_headmodeler
eval(command);

%------------------------------------------------------------%
% compute leadfield matrix and normalize it (depth weighting)%
%------------------------------------------------------------%
options = bst_headmodeler; % create structure
options.Method = {‘meg_os’}; % use meg overlapping sphere
options.ChannelType = ‘MEG’;
options.ImageGridBlockSize=7205; % put one point more than our cortical mesh to get mesh data
options.ChannelLoc=sensregreg1(1:3,:)./100; % my MEG sensor positions in m
options.ChannelOrient=sensregreg1(4:6,:)./100; % my MEG sensor orientations in m
options.SourceModel = -1;
options.VolumeSourceGrid = 0;
options.GridLoc=brain.vert’./100; % here I put the vertex points of my surface mesh
options.GridOrient = brain.normal’; % here normals
options.ApplyGridOrient=0; % 3 orientation per dipole
%options.ApplyGridOrient=1; % 1 orientation per dipole
options.SourceLoc = options.GridLoc;
options.SourceOrient = options.GridOrient;
options.Scalp.FileName = [outpath,vp{v},’_brst_scalp.mat’];
options.Scalp.iGrid = 1;
options.Verbose = 1;

[G,options] = bst_headmodeler(options); % G should contain your leadfield

Hope that helps,

best,

Stephan