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:

  1. 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.

  2. Select the Gray Mater volume elements for analysis.

  3. 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

    3.3. Create a mask to separate elements inside the ROI

  4. Calculate field in ROI using a weighted average.

    4.1. Get the volumes of the tetrahedra

    4.2. Get the field of interest

    4.3. Calculate an 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', '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(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 = 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 = 'normE'
field = gray_matter.field[field_name][:]

# Calculate the mean
mean_normE = np.average(field[roi], weights=elm_vols[roi])
print('mean ', field_name, ' in M1 ROI: ', mean_normE)

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', '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 = 'normE';
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 checking 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. This is not possible for headreco models ran with the --no-cat option

For this kind of analysis, we have to

  1. Load the subject_overlays/*_central.msh file

    • 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 files have the fields defined in thein nodes

    • SimNIBS also outputs each hemisphere as a FreeSurfer surface file and individual fields as a FreeSurfer curv file.

  2. Define the ROI

    2.1. Load an atlas transformed to subject space

    2.2. Select a region from the atlas

    Note

    Instead of an atlas, it is also possible to define the ROI using the same methods as described in the

  3. Calculate field in ROI using a weighted average.

    3.1. Calculate node areas

    3.2. Get the field of interest

    3.3. Calculate an 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', '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
field_name = 'E_normal'
mean_normE = np.average(gm_surf.field[field_name][roi], weights=node_areas[roi])
print('mean ', field_name, ' in ', region_name, ': ', mean_normE)

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', '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_normal';
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)