# Group Analysis¶

SimNIBS supports group analysis of electric fields in both FsAverage and MNI spaces.

In this example, we will simulate the average electric field for a simple tDCS montage across five subjects using the FsAverage space.

## Example Dataset¶

This example dataset is composed of a subgroup of a cohort available at OpenfMRI. The data was processed in SimNIBS 2.1 using headreco. For more information, please see the OSF repository and Saturnino et al., 2018.

In this example, we have five subjects, named sub01, sub09, sub10, sub12, and sub15.

## Objective¶

Suppose we apply tDCS on the five subjects using an anode over C3 and a cathode over AF4. We want to obtain the mean and standard deviation of the normal component of the electric field (that is, the incoming or outgoing electric field component) in the middle gray matter surface across all subjects.

Norm (top) and Norma (bottom) of the electric field, from Antonenko et al. 2019

## Set up and Run Simulations¶

There are several ways to set-up and run Simulations in SimNIBS

### GUI¶

Set-up the simulation for each subject in the Graphical User Interface. In this case, remember to tick the Transform to FsAverage option in the Simulation Options Window. (under Edit -> Simulation Options)

### Python¶

Write a Python script. In this case, remember to set map_to_fsavg to True in the SESSION structure. See Scripting Simulations for more information.

'''
This example runs tDCS simulations with a bipolar montage for five subjects
The dataset with the five head models is avaliable at https://osf.io/ah5eu/
please look at the "group_average" for how to do a simple analysis of the group data
'''
import os
from simnibs import sim_struct, run_simnibs

# Set the subjects
subjects = ['sub01', 'sub09', 'sub10', 'sub12', 'sub15']

# Set a TDCSLIST structure with the simulation set-up
tdcslist = sim_struct.TDCSLIST()
tdcslist.currents = [0.001, -0.001]

anode.channelnr = 1
anode.centre = 'C3'
anode.pos_ydir = 'C1'
anode.shape = 'rect'
anode.dimensions = [50, 50]
anode.thickness = 4

cathode.channelnr = 2
cathode.centre = 'AF4'
cathode.pos_ydir = 'F6'
cathode.shape = 'rect'
cathode.dimensions = [50, 70]
cathode.thickness = 4

# Run the simulation in each subject
for sub in subjects:
# ALWAYS create a new SESSION when changing subjects
s = sim_struct.SESSION()
s.map_to_fsavg = True
s.map_to_MNI = True
s.fields = 'eEjJ'
s.fnamehead = os.path.join(sub, sub + '.msh')
s.pathfem = os.path.join(sub, 'bipolar')
# Don't open in gmsh
s.open_in_gmsh = False
# Add the tdcslist we defined above
# Run the sumulation
run_simnibs(s)


### MATLAB¶

Write a MATLAB script. In this case, remember to set map_to_fsavg to True in the SESSION structure. See Scripting Simulations for more information.

%  This example runs tDCS simulations with a bipolar montage for five subjects
%  The dataset with the five head models is avaliable at https://osf.io/ah5eu/
%  please look at the "group_average" for how to do a simple analysis of the group data

% Set the subjects
subjects = {'sub01', 'sub09', 'sub10', 'sub12', 'sub15'};

% Start a SESSION
S = sim_struct('SESSION');
S.map_to_fsavg = true;
S.map_to_mni = true;
S.fields = 'eEjJ';

% Set a TDCSLIST with the simulation set-up
S.poslist{1} = sim_struct('TDCSLIST');
S.poslist{1}.currents = [0.001, -0.001];

S.poslist{1}.electrode(1).channelnr = 1;
S.poslist{1}.electrode(1).centre = 'C3';
S.poslist{1}.electrode(1).pos_ydir = 'C1';
S.poslist{1}.electrode(1).shape = 'rect';
S.poslist{1}.electrode(1).dimensions = [50, 50];
S.poslist{1}.electrode(1).thickness = 4;

S.poslist{1}.electrode(2).channelnr = 2;
S.poslist{1}.electrode(2).centre = 'AF4';
S.poslist{1}.electrode(2).pos_ydir = 'F6';
S.poslist{1}.electrode(2).shape = 'rect';
S.poslist{1}.electrode(2).dimensions = [50, 70];
S.poslist{1}.electrode(2).thickness = 4;

% Run the simulation in each subject
for i = 1:length(subjects)
sub = subjects{i};
S.pathfem = fullfile(sub, 'bipolar'); % Output directory
run_simnibs(S);
end


## Calculate Mean¶

When the simulations are over, we need to collect their results to calculate an average. In SimNIBS, we can do it either in Python or MATLAB. Please notice that, while for setting up simulations Python and MATLAB share a similar interface, in post-processing the interfaces can be very different.

### Python¶

'''
This example wil go throgh simulations and calculate
the average and the standard deviation of the normal component
of the electric field in FsAverage space

It is a follow-up to the "run_simulations" example
'''
import os
import numpy as np
import simnibs

subjects = ['sub01', 'sub09', 'sub10', 'sub12', 'sub15']
results_folder = os.path.join('bipolar', 'fsavg_overlays')
fsavg_msh_name = '_TDCS_1_scalar_fsavg.msh'
field_name = 'E_normal'

fields = []
for sub in subjects:
# read mesh with results transformed to fsaverage space
os.path.join(sub, results_folder, sub + fsavg_msh_name)
)
# save the field in each subject
fields.append(results_fsavg.field[field_name].value)

## Calculate and plot averages
# Calculate
fields = np.vstack(fields)
avg_field = np.mean(fields, axis=0)
std_field = np.std(fields, axis=0)

# Plot
results_fsavg.nodedata = [] # cleanup fields

# show surface with the fields
results_fsavg.view(visible_fields='E_normal_avg').show()

## Calculate average in an ROI defined using an atlas
# load atlas and define a region
atlas = simnibs.get_atlas('HCP_MMP1')
region_name = 'lh.4'
roi = atlas[region_name]
# visualize region
results_fsavg.view(visible_fields=region_name).show()

# calculate mean field using a weighted mean
node_areas = results_fsavg.nodes_areas()
avg_field_roi = np.average(avg_field[roi], weights=node_areas[roi])
print(f'Average {field_name} in {region_name}: ', avg_field_roi)


### MATLAB¶

% This example wil go throgh simulations and calculate
% the average and the standard deviation of the normal component
% of the electric field in FsAverage space
%
% It is a follow-up to the "run_simulations" example

subjects = {'sub01', 'sub09', 'sub10', 'sub12', 'sub15'};
results_folder = fullfile('bipolar', 'fsavg_overlays');
fsavg_msh_name = '_TDCS_1_scalar_fsavg.msh';
field_name = 'E_normal';

fields = {};
for i = 1:length(subjects)
sub = subjects{i};
% load mesh with results transformed to fsaverage space
m = mesh_load_gmsh4(fullfile(pwd, sub, results_folder, [sub fsavg_msh_name]));
% Save the field of each subject
fields{i} = m.node_data{get_field_idx(m, field_name, 'node')}.data;
end
%% Calculate and plot averages
% Calculate
fields = cell2mat(fields);
avg_field = mean(fields, 2);
std_field = std(fields, 0, 2);
% Plot
m.node_data = {}; %cleanup fields
m.node_data{1}.data = avg_field; % add average field
m.node_data{1}.name = [field_name '_avg'];
m.node_data{2}.data = std_field; % add std field
m.node_data{2}.name = [field_name '_std'];

% show surfaces with fields
mesh_show_surface(m, 'field_idx', [field_name '_avg'])
mesh_show_surface(m, 'field_idx', [field_name '_std'])

%% Calculate average in an ROI defined using an atlas
% load atlas and define a region
region_name = 'lh.4';
roi_idx=find(strcmpi(snames, region_name));
node_idx = m.node_data{end}.data==roi_idx;
% visualize region
m.node_data{end+1}.data = int8(node_idx);
m.node_data{end}.name = region_name;
mesh_show_surface(m, 'field_idx', region_name)

% calculate a weighted mean using the node areas
nodes_areas = mesh_get_node_areas(m);
avg_field_roi = sum(avg_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)