-
Notifications
You must be signed in to change notification settings - Fork 2
50. TMS E‐field modeling with Matlab scripting
We use the Matlab toolbox provided by Drs. Sergey N Makarov and Aapo Nummenmaa to calculate the induced electric field in TMS.
The whole script explained below is [here].
Refer to the page using GUI to do TMS e-field modeling.
To model electric fields (e-fields) using realistic brain models and TMS coil models, the following steps will be involved:
-
Prepare "head models", which are realistic 3D models of head/skull/brain geometries and tissue properties, such as conductivity.
-
Prepare a "TMS coil model", which is a realistic 3D model of TMS coil shape with directed current flows.
-
Prepare the relative position between head models and TMS coil model. This process will be adjusted interactively, potentially including guidance from fMRI.
-
Calculate the e-fields.
-
Render calculated electric field strengths, such as the total strength (regardless of x-, y-, z-directional components) or directional component strength (such as component that is oriented tangentially or perpendicular to the cortical surface).
We will model the electric field generated by the MRi-B91 coil (MagVenture).
We will model the electric field generated over a subject coded s002, whose FreeSurfer reconstruction is at /space_lin2/fhlin/eegmri_memory/subjects/s002 at Sunnybrook.
We define environment variables subjects_dir and subject. This example focuses on the e-field distribution in the left hemisphere.
Definitions of head models, including their respective surface geometry (.surf or ?h. files in FreeSurfer; "file_surf") and name (" output_file_surf"). Matlab files named after subject and output_file_surf will be created.
Definition of boundaries, including their respective name ("tissue_name"), conductivity ("tissue_conductivity"), and the name of the enclosing tissue ("tissue_enclosing") will be created as file_bem.
head_surf_stem defines the scalp surface where a TMS coil will be touching upon. It must be an entry among output_file_surf.
TMS coil named as tms_coil_name will be created and prepared.
target_coord is the TMS target location in FreeSurfer surface coordinate system.
flag_nav indicates if a GUI is available for display.
close all; clear all;
subjects_dir='/Users/fhlin/workspace/eegmri_memory/subjects';
%subjects_dir='/space_lin2/fhlin/eegmri_memory/subjects';
subject='s002';
hemi='lh';
file_surf={
{'outer_skin.surf'}, % scalp
{'outer_skull.surf'}, % outer skull
{'inner_skull.surf'}, % inner skull
{'../surf/lh.orig'}, % brain
{'../surf/rh.orig'}, % brain
};
output_file_surf={
'skin',
'skull',
'csf',
'gm_lh',
'gm_rh',
};
head_surf_stem='skin'; %file stem for scalp; must be one of output_file_surf
tissue_name={
'Skin'; %scalp
'Skull'; %outer skull
'CSF'; %inner skull
'GM_LH'; %brain
'GM_RH'; %brain
};
% Electric field calculations in brain stimulation based on finite elements: An optimized processing pipeline for the generation and usage of accurate individual head models Hum. Brain Mapp., 34 (4) (2013), pp. 923-935, 10.1002/hbm.21479
tissue_conductivity=[
0.465; % scalp (below scalp)
0.010; % skull (below outer skull)
1.654; % CSF (below inner skull)
0.275; % brain (below brain)
0.275; % brain (below brain)
];
tissue_enclosing={
'FreeSpace'; % ouside scalp
'Skin'; % outside outer skull
'Skull'; % outside inner skull
'CSF'; % otuside brain
'CSF'; % otuside brain
};
file_bem='tissue_index_bem.txt';
tms_coil_name='MagVenture_MRiB91';
target_coord=[-5.1 56.7 59.4];
flag_nav=0; %open navigation window
We specify locations vertex_coords where e-field will be estimated. These coordinates are in FreeSurfer surface coordinate system.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
setenv('SUBJECTS_DIR',subjects_dir);
switch(lower(hemi))
case 'lh'
[vertex_coords, faces] = read_surf(sprintf('%s/%s/surf/lh.orig',subjects_dir,subject));
case 'rh'
[vertex_coords, faces] = read_surf(sprintf('%s/%s/surf/rh.orig',subjects_dir,subject));
end;
Head models are then prepared. A head model, indexed by head_surf_idx, acting as the boundary for placing a TMS coil will be found.
%prepare BEM
[status, bem_obj]=etc_tms_prepare_bem(subject,file_surf,output_file_surf,tissue_name,tissue_conductivity,tissue_enclosing,file_bem);
[a,head_surf_idx]=ismember(head_surf_stem, output_file_surf);
if(a<eps)
str=sprintf('%s,',output_file_surf{:});
fprintf('Error! Scalp surface [%s] is not among surfaces {%s}\n', head_surf_stem, str(1:end-1));
return;
else
fprintf('Scalp surface [%s] found: <%03d>::<%s>\n', head_surf_stem, head_surf_idx, output_file_surf{head_surf_idx});
end;
[status, bem_t, bem_P, normals, Center, Area, Indicator, name, tissue, cond, enclosingTissueIdx, condin, condout, contrast, tneighbor, RnumberE, ineighborE, EC, file_mesh, file_meshp] =etc_tms_efield_prep_model(file_bem);
Now create the TMS coil model for visualization (coil_P and coil_t) and e-field calculation (strcoil object). Coil position tms_coil_origin and orientation (tms_coil_axis and tms_coil_up) are initialized. The TMS coil transformation matrix (tms_coil_xfm is also initialized. If GUI is enabled (flag_nav=1;), TMS coil location and orientation will be prepared.
%prepare TMS coil (initialization)
[status, strcoil, coil_P, coil_t, coil_tind]=etc_tms_prepare_coil(tms_coil_name);
%%% this is preparing navigation window objects for subsequent rendering, if any
if(flag_nav)
global etc_render_fsbrain
if(isfield(etc_render_fsbrain,'app_tms_nav'))
[status, tms_coil_origin, tms_coil_axis, tms_coil_up, tms_coil_xfm] = etc_tms_init_nav(etc_render_fsbrain.app_tms_nav, strcoil, coil_P, coil_t);
else
[status, tms_coil_origin, tms_coil_axis, tms_coil_up, tms_coil_xfm] = etc_tms_init_nav([], strcoil, coil_P, coil_t,'subject',subject);
end;
else
tms_coil_origin=[0 0 0];
tms_coil_axis=[0 0 -1];
tms_coil_up=[0 1 0];
tms_coil_xfm=eye(4);
end;
Move and orient the TMS to a specified target defined by target_coord. A head surface indexed by head_surf_idx is used to constrain the placement of the coil. In this targeting process, the TMS coil is placed with its axis (tms_coil_axis) passing the target coordnate, the coil plane touching the head surface, and coil orientation with the "up" direction (tms_coil_up) oriented toward the head vertex.
To prepare the target from an MNI coordinate, use the following codes:
Here we assume that a subject has been prepared/shown by etc_render_fsbrain.m. The variable surface_coord is the coordinate in the FreeSurfer surface coordinate system.
MNI_coord= [0 0 0]; %example of the target in MNI coordinate
mni=[MNI_coord(:)' 1]';
click_vertex_vox=inv(etc_render_fsbrain.vol.vox2ras)*inv(etc_render_fsbrain.vol_pre_xfm)*inv(etc_render_fsbrain.talxfm)*mni;
click_vertex_vox=click_vertex_vox(1:3)';
surface_coord=etc_render_fsbrain.vol.tkrvox2ras*[click_vertex_vox(:); 1];
surface_coord=surface_coord(1:3)
The strcoil object is also moved and oriented accordingly.
The TMS coil transformation matrix tms_coil_xfm_moved describes the coil position and orientation after targeting.
If GUI is enabled flag_nav=1;, TMS coil in the GUI is also updated.
% move TMS coil to target
[tms_coil_xfm_moved]=etc_tms_target_xfm_goto(target_coord, bem_obj(head_surf_idx), tms_coil_origin, tms_coil_axis, tms_coil_up, tms_coil_xfm);
%move strcoil
[strcoil, strcoil_xfm] = etc_tms_target_xfm_apply(strcoil, tms_coil_xfm_moved(1:3,4).*1e3, -tms_coil_xfm_moved(1:3,3), tms_coil_xfm_moved(1:3,2), tms_coil_xfm, tms_coil_xfm_moved);
%%% this is updating navigation window objects for subsequent rendering, if any
if(flag_nav)
results = etc_tms_target_xfm_apply_nav(etc_render_fsbrain.app_tms_nav, tms_coil_xfm_moved(1:3,4).*1e3, -tms_coil_xfm_moved(1:3,3), tms_coil_xfm_moved(1:3,2));
end;
tms_coil_center_now=tms_coil_xfm_moved(1:3,4);
tms_coil_axis_now= -tms_coil_xfm_moved(1:3,3);
tms_coil_up_now=tms_coil_xfm_moved(1:3,2);
Calculate the e-field including its total strength, normal component, and tangential component at coords (in millimeter).
%calculate the e-field
coords=vertex_coords./1e3;
[status, efield1]=etc_tms_efield_surf(bem_t, bem_P, normals, Center, Area, Indicator, name, tissue, cond, enclosingTissueIdx, condin, condout, contrast, tneighbor, RnumberE, ineighborE, EC, coords,'tissue_to_plot','GM_LH');
E-field is shown over the cortical surface, if GUI is available.
if(flag_nav)
%the following is for rendering E-field
etc_render_fsbrain.overlay_vertex=efield1.vertices;
etc_render_fsbrain.overlay_value=efield1.Etotal;
etc_render_fsbrain.overlay_smooth=[];
etc_render_fsbrain.overlay_source=1;
if(isempty(etc_render_fsbrain.overlay_vol_stc))
etc_render_fsbrain_overlay_vol_init;
else
etc_render_fsbrain.overlay_vol_stc=etc_render_fsbrain.overlay_value(:);
end;
etc_render_fsbrain.overlay_flag_render=1;
etc_render_fsbrain_handle('update_overlay_vol');
etc_render_fsbrain_handle('redraw');
etc_render_fsbrain_handle('draw_pointer');
end;
To further tune the TMS coil position and orientation, additional adjustment can be added. Here is an example of rotating the TMS coil by 30 degrees. The coil transform matrix reflecting the tuned position/orientation is created tms_coil_xfm_moved_tuned. After tuning, TMS coil, including strcoil is updated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% tuning loop, if any, starts here. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rotation; 30 degrees
[tms_coil_xfm_moved_tuned]=etc_tms_target_xfm_tune(target_coord, bem_obj(head_surf_idx), tms_coil_xfm_moved(1:3,4).*1e3, -tms_coil_xfm_moved(1:3,3), tms_coil_xfm_moved(1:3,2), tms_coil_xfm_moved, 4, 30);
%move strcoil
[strcoil, strcoil_xfm] = etc_tms_target_xfm_apply(strcoil, tms_coil_xfm_moved(1:3,4).*1e3, -tms_coil_xfm_moved(1:3,3), tms_coil_xfm_moved(1:3,2), strcoil_xfm, tms_coil_xfm_moved_tuned);
%%% this is updating navigation window objects for subsequent rendering, if any
if(flag_nav)
results = etc_tms_target_xfm_apply_nav(etc_render_fsbrain.app_tms_nav, tms_coil_xfm_moved_tuned(1:3,4).*1e3, -tms_coil_xfm_moved_tuned(1:3,3), tms_coil_xfm_moved_tuned(1:3,2));
end;
Calculate and display the e-field with the TMS coil at the updated position/orientation.
%calculate the e-field
coords=vertex_coords./1e3;
[status, efield2]=etc_tms_efield_surf(bem_t, bem_P, normals, Center, Area, Indicator, name, tissue, cond, enclosingTissueIdx, condin, condout, contrast, tneighbor, RnumberE, ineighborE, EC, coords,'tissue_to_plot','GM_LH');
if(flag_nav)
%the following is for rendering E-field
etc_render_fsbrain.overlay_vertex=efield2.vertices;
etc_render_fsbrain.overlay_value=efield2.Etotal;
etc_render_fsbrain.overlay_smooth=[];
etc_render_fsbrain.overlay_source=1;
if(isempty(etc_render_fsbrain.overlay_vol_stc))
etc_render_fsbrain_overlay_vol_init;
else
etc_render_fsbrain.overlay_vol_stc=etc_render_fsbrain.overlay_value(:);
end;
etc_render_fsbrain.overlay_flag_render=1;
etc_render_fsbrain_handle('update_overlay_vol');
etc_render_fsbrain_handle('redraw');
etc_render_fsbrain_handle('draw_pointer');
end;
Results are then stored in a file.
save tms_script_test0.mat efield1 efield2