Analyzing Simulations¶
SimNIBS offers tools to help analyzing your simulation results. Most of the tools have both Python and MATLAB versions. However, their arguments and outputs might be slightly different.
To run the analysis scripts, please follow the same steps needed to run the simulation scripts.
ROI Analysis Using MNI Coordinates¶
The objective of this example is to evaluate the mean electric field in an ROI, defined using MNI coordinates.
This analysis follows the steps:
Load a simulation output.
The default SimNIBS output is a mesh, that is, many nodes scattered in space forming tetrahedra and triangles. Electric fields (and also current density fields) are defined in each tetrahedra and triangle.
Select the gray matter volume elements for analysis:
- Python:
simnibs.Msh.crop_mesh()method
- MATLAB:
mesh_extract_regions function
Define the ROI.
3.1. Use the mni2subject_coords function to transform coordinates from MNI to subject space.
3.2. Get the center of the tetrahedra:
- Python:
simnibs.Msh.elements_baricenters()method
- MATLAB:
mesh_get_tetrahedron_centers function
3.3. Create a mask to separate elements inside the ROI.
Calculate field in ROI using a weighted average.
4.1. Get the volumes of the tetrahedra:
- Python:
simnibs.Msh.elements_volumes_and_areas()method
- MATLAB:
mesh_get_tetrahedron_sizes function
4.2. Get the field of interest:
- Python:
simnibs.Msh.field()attribute
- MATLAB:
get_field_idx function
4.3. Calculate a weighted average:
Python¶
'''
Simple ROI analysis of the electric field from a simulation.
We will calculate the mean electric field in a gray matter ROI defined around M1
'''
import os
import numpy as np
import simnibs
## Load simulation result
# Read the simulation result
head_mesh = simnibs.read_msh(
os.path.join('tdcs_simu', 'ernie_TDCS_1_scalar.msh')
)
# Crop the mesh so we only have gray matter volume elements (tag 2 in the mesh)
gray_matter = head_mesh.crop_mesh(simnibs.ElementTags.GM)
## Define the ROI
# Define M1 from MNI coordinates (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2034289/)
# the first argument is the MNI coordinates
# the second argument is the subject "m2m" folder
ernie_coords = simnibs.mni2subject_coords([-37, -21, 58], 'm2m_ernie')
# we will use a sphere of radius 10 mm
r = 10.
# Electric fields are defined in the center of the elements
# get element centers
elm_centers = gray_matter.elements_baricenters()[:]
# determine the elements in the ROI
roi = np.linalg.norm(elm_centers - ernie_coords, axis=1) < r
# get the element volumes, we will use those for averaging
elm_vols = gray_matter.elements_volumes_and_areas()[:]
## Plot the ROI
gray_matter.add_element_field(roi, 'roi')
gray_matter.view(visible_fields='roi').show()
## Get field and calculate the mean
# get the field of interest
field_name = 'magnE'
field = gray_matter.field[field_name][:]
# Calculate the mean
mean_magnE = np.average(field[roi], weights=elm_vols[roi])
print('mean ', field_name, ' in M1 ROI: ', mean_magnE)
MATLAB¶
% Simple ROI analysis of the electric field from a simulation.
% We will calculate the mean electric field in a gray matter ROI
% The ROI is defined using an MNI coordinates
%% Load Simulation Result
% Read the simulation result
head_mesh = mesh_load_gmsh4(...
fullfile('tdcs_simu', 'ernie_TDCS_1_scalar.msh') ...
);
% Crop the mesh so we only have gray matter volume elements (tag 2 in the mesh)
gray_matter = mesh_extract_regions(head_mesh, 'region_idx', 2);
%% Define the ROI
% Define M1 from MNI coordinates (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2034289/)
% the first argument is the MNI coordinates
% the second argument is the subject "m2m" folder
ernie_coords = mni2subject_coords([-37, -21, 58], fullfile(pwd, 'm2m_ernie'));
% we will use a sphere of radius 10 mm
r = 10;
% Electric fields are defined in the center of the elements
% get element centers
elm_centers = mesh_get_tetrahedron_centers(gray_matter);
% determine the elements in the ROI
roi = sqrt(sum(bsxfun(@minus, elm_centers, ernie_coords).^2, 2)) < r;
% get element volumes, we will use those for averaging
elm_vols = mesh_get_tetrahedron_sizes(gray_matter);
%% Get field and calculate the mean
% Get the field of interest
field_name = 'magnE';
field_idx = get_field_idx(gray_matter, field_name, 'elements');
field = gray_matter.element_data{field_idx}.tetdata;
% Calculate the mean
avg_field_roi = sum(field(roi) .* elm_vols(roi))/sum(elm_vols(roi));
fprintf('mean %s in ROI: %f\n', field_name, avg_field_roi)
ROI Analysis Using Surfaces¶
An alternative analysis method is to use the middle cortical surfaces created by SimNIBS.
Note
This analysis method requires that the simulation has been interpolated to a cortical surface. This can be done either by setting the Interpolate to cortical surface option in the Simulation Options Window. (under Edit -> Simulation Options), or by setting the map_to_surf option in the SESSION structure.
For this kind of analysis, we have to
Load the
subject_overlays/*_central.mshfile.This file contains both left and right hemispheres joined together in that order, as well as all the output fields.
Differently to the original mesh, this file has the fields defined in the nodes
SimNIBS also outputs each hemisphere as a FreeSurfer surface file and individual fields as a FreeSurfer curv file.
Define the ROI.
2.1. Load an atlas transformed to subject space:
- Python:
simnibs.subject_atlas()function
- MATLAB:
subject_atlas function
2.2. Select a region from the atlas.
Calculate field in ROI using a weighted average.
3.1. Calculate node areas:
- Python:
simnibs.Msh.nodes_areas()method
- MATLAB:
mesh_get_node_areas function
3.2. Get the field of interest:
- Python:
simnibs.Msh.field()attribute
- MATLAB:
get_field_idx function
3.3. Calculate a weighted average.
Python¶
'''
ROI analysis of the electric field from a simulation using an atlas.
We will calculate the mean electric field in a gray matter ROI defined using an atlas
'''
import os
import numpy as np
import simnibs
## Input ##
# Read the simulation result mapped to the gray matter surface
gm_surf = simnibs.read_msh(
os.path.join('tdcs_simu', 'subject_overlays',
'ernie_TDCS_1_scalar_central.msh')
)
# Load the atlas and define the brain region of interest
atlas = simnibs.subject_atlas('HCP_MMP1', 'm2m_ernie')
region_name = 'lh.4'
roi = atlas[region_name]
# plot the roi
gm_surf.add_node_field(roi, 'ROI')
gm_surf.view(visible_fields='ROI').show()
# calculate the node areas, we will use those later for averaging
node_areas = gm_surf.nodes_areas()
# finally, calculate the mean of the field strength
field_name = 'E_magn'
mean_magnE = np.average(gm_surf.field[field_name][roi], weights=node_areas[roi])
print('mean ', field_name, ' in ', region_name, ': ', mean_magnE)
MATLAB¶
% ROI analysis of the electric field from a simulation using an atlas.
% We will calculate the mean electric field in a gray matter
% ROI defined using an atlas
% Read the simulation results interpolated to the middle gray matter surface
surf = mesh_load_gmsh4(...
fullfile('tdcs_simu', 'subject_overlays', ...
'ernie_TDCS_1_scalar_central.msh')...
);
% Load the atlas and define the brain region of interest
[labels, snames] = subject_atlas(surf, fullfile(pwd, 'm2m_ernie'), 'HCP_MMP1');
region_name = 'lh.4';
roi_idx=find(strcmpi(snames, region_name));
node_idx = labels.node_data{end}.data==roi_idx;
% Plot the ROI
surf.node_data{end+1}.data = int8(node_idx);
surf.node_data{end}.name = region_name;
mesh_show_surface(surf, 'field_idx', region_name)
% calculate the node areas, we will use those later for averaging
nodes_areas = mesh_get_node_areas(surf);
% Get the field of interest
field_name = 'E_magn';
field_idx = get_field_idx(surf, field_name, 'node');
field = surf.node_data{field_idx}.data;
% Calculate the mean
avg_field_roi = sum(field(node_idx).*nodes_areas(node_idx))/sum(nodes_areas(node_idx));
fprintf('Average %s in %s: %f\n', field_name, region_name, avg_field_roi)
