Brain-fingerprinting

Authors: Jason da Silva Castanheira, Raymundo Cassani

This tutorial introduces the concept of neurophysiological brain-fingerprinting and its implementation in Brainstorm. To follow this tutorial, please first complete the OMEGA tutorial as we will derive the brain-fingerprints of its five participants.

Contents

Introduction
Preparing the data for brain-fingerprinting
The brain-fingerprinting method
Differentiation accuracy
Differentiability
Quantifying the most salient features of the brain-fingerprint
Applications of the brain-fingerprint
Additional documentation
Scripting

Introduction

Brain-fingerprinting is a method to assess inter-individual differences in brain activity. It was first proposed by Finn and colleagues in 2015. This tutorial will derive spectral brain-fingerprints using neural power spectra at computed within each parcel of the Destrieux atlas. Once the individuals' brain-fingerprints have been obtained, differentiation accuracy and differentiability are computed for the cohort. Finally, we analyze relative contribution of each of the cortical parcels (i.e., brain regions) and frequency bins to participant differentiation

This demonstration uses the pre-processed data from five participants obtained from the outputs of the OMEGA tutorial. Please run the entire tutorial before using the tutorial_brain_fingerprint.m script (shown at the end of this page).

Preparing the data for brain-fingerprinting

The brain-fingerprint of an individual can be derived from functional connectomes or neural power spectrum, as detailed in da Silva Castanheira et al., 2021. Note that this is not an exhaustive list of all possible neurophysiological features that can be used to define brain-fingerprints.

From the estimated surface source maps, we compute the power spectrum density (i.e., PSD) using the Welch method at every vertex of the cortical surface, and average the resulting spectra within each region defined by the Destrieux atlas (see Scout page for more details). Note, this is done to downsample the spatial features of the brain-fingerprint. If we did not downsample the PSD to a cortical atlas, the resulting feature space would be too large to work with (e.g., 15,000 vertices × 300 frequency bins). The choice of atlas and neurophysiological features (e.g., functional connectomes vs power spectra) to define the brain-fingerprint is dependent on the user's specific hypotheses. See the section of the script entitled Compute PSD for all ROIs of an atlas

Brain-fingerprinting requires at least two data segments for every individual. This can be two separate recordings (i.e., between-session fingerprinting), or two segments within the same recording session (i.e., within-session fingerprinting). For this demonstration, we will split the MEG recordings of the OMEGA tutorial into two parts (data segments). As such, for each participant, there will be two brain-fingerprints (b-fp1 and b-fp2), each consisting of a vector of nParcellations×nFrequencies elements. See the image below for a summary of the pipeline:

AlternativeText

The brain-fingerprinting method

The brain-fingerprinting method is based on the correlational similarity of participants across data segments (see image above). For each participant in a given cohort, we compute the Pearson correlation coefficient between the brain-fingerprint of the first segment and second segment of all participants in the cohort, including the probe participant. This yields a participant similarity matrix, where off-diagonal elements represent the similarity of a participant to all other participants and diagonal elements represent the similarity of a participant's brain-fingerprint across data segments (i.e., similarity to their own brain-fingerprint). See the image below for an example participant similarity matrix.

Participant differentiation consists of a lookup procedure along the rows (or columns) of the participant similarity matrix. A participant is said to be correctly differentiated if a participant's brain-fingerprint is more similar to themselves (self-similarity, also referred to as Iself) than all other participants in the cohort (other-similarity, also referred to as Iothers). In other words, if the brain-fingerprint of participant i is most similar to their brain-fingerprint taken at the second time point, this participant is said to be correctly differentiated.

Differentiation accuracy

Differentiation accuracy represents the percent ratio of the number of participants correctly differentiated across the cohort (i.e., differentiation accuracy). Brain-fingerprinting is typically repeated for all possible pairs of data segments. In the case of this tutorial, there are only two data segments, the first and second half of the recording. We therefore obtain two differentiation accuracies along the rows and columns of the participant similarity matrix respectively. In the case where three data segments are available, we can compute six differentiation accuracies: the accuracy for the columns and rows of data segment 1 vs segment 2, and for each possible pair of recordings (i.e., data segment 2 vs 3 and data segment 1 vs 3).

Differentiability

Differentiability is a measure which describes how easily a given participant is to differentiate from a cohort of individuals. This measure consists of scaling (z-scoring) the self-similarity of brain-fingerprints against the other-similarity (off-diagonal elements). A participant with a high differentiability score will have a higher self-similarity relative to their similarity to others in the cohort (other-similarity) and therefore are easily differentiated (see the last row in the example participant similarity matrix below). In contrast, a participant whose brain-fingerprint is most similar to others in the cohort will be more challenging to differentiate (see the second row of the participant similarity matrix below). AlternativeText

Note that other measures to quantify the relative easy of differentiating a participant exist: Idiff represents the difference between Iself and Iothers. See Amico & Goñi 2018 and Sareen at el., 2021.

The participant similarity matrix is generally symmetric (although it does not necessarily have to be). This implies that the differentiability across rows and columns are generally strongly correlated, as such, the participant differentiability can be computed as the mean of these two per participant

Quantifying the most salient features of the brain-fingerprint

Beyond differentiation accuracy and differentiability, we can quantify the relative contribution of cortical parcels (i.e., brain regions) and frequency bins to participant differentiation.

To do so, we calculated intraclass correlations (ICC). ICC quantifies the ratio of within-participant variance and between-participant variance, with participants as their own raters across data segments (see Amico et al., 2018 & da Silva Castanheira et al., 2021 for details). Features that contribute the most to participant differentiation should be highly consistent within individuals, but show great inter-individual variance. Larger values of ICC indicate that a neurophysiological feature (e.g., region of interest or frequency band) contributes more to participant differentiation than smaller values.

See the subsection entitled IntraClass Correlations (ICC) under the section SIMILARITY AND DIFFERENTIABILITY in the script.

In this tutorial we plot the average ICC value per cortical parcel within each frequency band of interest as defined in the script. See image below:

% Bands for ICC cortex plot
BandNames = {'theta', 'alpha', 'beta', 'gamma', 'highgamma'};
BandLowerFreqs = [ 4,  8, 13, 30,  50]; % Hz, inclusive
BandUpperFreqs = [ 8, 13, 30, 50, 150]; % Hz, exclusive

In the Brainstorm database, the ICC values are saved: as (i) Matrix files (which can be exported to tabular files) and, as (ii) cortical maps (which allow to visualize their spatial distribution). Both set of files are found in the Group analysis folder. The image below presents the database (left), and ICC values for theta band, as cortical maps (center), and as matrix (right).

AlternativeText

Applications of the brain-fingerprint

Beyond biometric differentiation, brain-fingerprints advance the neurobiological origins of individual traits and inform clinical. Below we summarize some of the most recent applications of brain-fingerprints.

Decoding individual traits: Previous work demonstrates that beyond biometric identification, brain-fingerprints can predict individual differences in intelligence test performance (e.g., Finn et al., 2015; Amico et al., 2018) and demographics (da Silva Castanheira et al., 2021).
Clinical application: Recent work has explored the clinical utility of bran-fingerprints. This work has shown across multiple neurological diseases that the brain-fingerprints of patients show increased variability (lower self-similairty/ Iself ) in brain-fingerprints of clinical populations (e.g., Kaufmann et al., 2017, 2018; Sorrentino et al., 2021; Troisi Lopez et al., 2023; da Silva Castanheira et al., 2024). A potential explanation for the de-individualization of brain activity in observed in Parkinson’s disease may be an inherent instability in the arrhythmic brain activity of patients over short periods. In contrast, rhythmic brain activity appears particularly characteristic of patients with PD (da Silva Castanheira et al., 2024). Future work is exploring the application of brain-fingerprints to the design of personalized medicine.

Additional documentation

Please cite previous work on brain-fingerprinting:

da Silva Castanheira, J., Orozco Perez, H.D., Misic, B. et al. Brief segments of neurophysiological activity enable individual differentiation. Nat Commun 12, 5713 (2021). https://doi.org/10.1038/s41467-021-25895-8
Finn, E., Shen, X., Scheinost, D. et al. Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity. Nat Neurosci 18, 1664–1671 (2015). https://doi.org/10.1038/nn.4135

Tutorials

Tutorial: Source estimation
Tutorial: Scouts
Tutorial: Power spectrum
Tutorial: Scripting

Scripting

The following script from the Brainstorm distribution reproduces the analysis presented in this tutorial page: brainstorm3/toolbox/script/tutorial_brain_fingerprint.m

function tutorial_brain_fingerprint(ProtocolNameOmega, reports_dir)
% TUTORIAL_BRAIN_FINGERPRINT: Script that reproduces the results of the online tutorial "Brain-fingerprint".
%
% CORRESPONDING ONLINE TUTORIAL:
%     https://neuroimage.usc.edu/brainstorm/Tutorials/BrainFingerprint
%
% INPUTS:
%    - ProtocolNameOmega : Name of the protocol created with all the data imported (TutorialOmega)
%    - reports_dir       : Directory where to save the execution report (instead of displaying it)
%
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
% 
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Jason da Silva Castanheira, Raymundo Cassani, 2024


%% ===== CHECK PROTOCOL =====
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Output folder for reports
if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
    reports_dir = [];
end
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin < 1) || isempty(ProtocolNameOmega)
    ProtocolNameOmega = 'TutorialOmega';
end


%% ===== BRAIN-FINGERPRINT PARAMETERS =====
% Subjects
SubjectNames = {'sub-0002', 'sub-0003', 'sub-0004', 'sub-0006', 'sub-0007'};
nSubjects = length(SubjectNames);
% Frequency range
LowerFreq =   4; % Hz, inclusive
UpperFreq = 150; % Hz, exclusive
% Cortical parcellation (atlas)
    % Note we downsample the cortical surface using an atlas for computation efficiency
AtlasName = 'Destrieux';
AtlasScouts = {'G_Ins_lg_and_S_cent_ins L', 'G_Ins_lg_and_S_cent_ins R', 'G_and_S_cingul-Ant L', 'G_and_S_cingul-Ant R', 'G_and_S_cingul-Mid-Ant L', 'G_and_S_cingul-Mid-Ant R', 'G_and_S_cingul-Mid-Post L', 'G_and_S_cingul-Mid-Post R', 'G_and_S_frontomargin L', 'G_and_S_frontomargin R', 'G_and_S_occipital_inf L', 'G_and_S_occipital_inf R', 'G_and_S_paracentral L', 'G_and_S_paracentral R', 'G_and_S_subcentral L', 'G_and_S_subcentral R', 'G_and_S_transv_frontopol L', 'G_and_S_transv_frontopol R', 'G_cingul-Post-dorsal L', 'G_cingul-Post-dorsal R', 'G_cingul-Post-ventral L', 'G_cingul-Post-ventral R', 'G_cuneus L', 'G_cuneus R', 'G_front_inf-Opercular L', 'G_front_inf-Opercular R', 'G_front_inf-Orbital L', 'G_front_inf-Orbital R', 'G_front_inf-Triangul L', 'G_front_inf-Triangul R', 'G_front_middle L', 'G_front_middle R', 'G_front_sup L', 'G_front_sup R', 'G_insular_short L', 'G_insular_short R', 'G_oc-temp_lat-fusifor L', 'G_oc-temp_lat-fusifor R', 'G_oc-temp_med-Lingual L', 'G_oc-temp_med-Lingual R', 'G_oc-temp_med-Parahip L', 'G_oc-temp_med-Parahip R', 'G_occipital_middle L', 'G_occipital_middle R', 'G_occipital_sup L', 'G_occipital_sup R', 'G_orbital L', 'G_orbital R', 'G_pariet_inf-Angular L', 'G_pariet_inf-Angular R', 'G_pariet_inf-Supramar L', 'G_pariet_inf-Supramar R', 'G_parietal_sup L', 'G_parietal_sup R', 'G_postcentral L', 'G_postcentral R', 'G_precentral L', 'G_precentral R', 'G_precuneus L', 'G_precuneus R', 'G_rectus L', 'G_rectus R', 'G_subcallosal L', 'G_subcallosal R', 'G_temp_sup-G_T_transv L', 'G_temp_sup-G_T_transv R', 'G_temp_sup-Lateral L', 'G_temp_sup-Lateral R', 'G_temp_sup-Plan_polar L', 'G_temp_sup-Plan_polar R', 'G_temp_sup-Plan_tempo L', 'G_temp_sup-Plan_tempo R', 'G_temporal_inf L', 'G_temporal_inf R', 'G_temporal_middle L', 'G_temporal_middle R', 'Lat_Fis-ant-Horizont L', 'Lat_Fis-ant-Horizont R', 'Lat_Fis-ant-Vertical L', 'Lat_Fis-ant-Vertical R', 'Lat_Fis-post L', 'Lat_Fis-post R', 'Pole_occipital L', 'Pole_occipital R', 'Pole_temporal L', 'Pole_temporal R', 'S_calcarine L', 'S_calcarine R', 'S_central L', 'S_central R', 'S_cingul-Marginalis L', 'S_cingul-Marginalis R', 'S_circular_insula_ant L', 'S_circular_insula_ant R', 'S_circular_insula_inf L', 'S_circular_insula_inf R', 'S_circular_insula_sup L', 'S_circular_insula_sup R', 'S_collat_transv_ant L', 'S_collat_transv_ant R', 'S_collat_transv_post L', 'S_collat_transv_post R', 'S_front_inf L', 'S_front_inf R', 'S_front_middle L', 'S_front_middle R', 'S_front_sup L', 'S_front_sup R', 'S_interm_prim-Jensen L', 'S_interm_prim-Jensen R', 'S_intrapariet_and_P_trans L', 'S_intrapariet_and_P_trans R', 'S_oc-temp_lat L', 'S_oc-temp_lat R', 'S_oc-temp_med_and_Lingual L', 'S_oc-temp_med_and_Lingual R', 'S_oc_middle_and_Lunatus L', 'S_oc_middle_and_Lunatus R', 'S_oc_sup_and_transversal L', 'S_oc_sup_and_transversal R', 'S_occipital_ant L', 'S_occipital_ant R', 'S_orbital-H_Shaped L', 'S_orbital-H_Shaped R', 'S_orbital_lateral L', 'S_orbital_lateral R', 'S_orbital_med-olfact L', 'S_orbital_med-olfact R', 'S_parieto_occipital L', 'S_parieto_occipital R', 'S_pericallosal L', 'S_pericallosal R', 'S_postcentral L', 'S_postcentral R', 'S_precentral-inf-part L', 'S_precentral-inf-part R', 'S_precentral-sup-part L', 'S_precentral-sup-part R', 'S_suborbital L', 'S_suborbital R', 'S_subparietal L', 'S_subparietal R', 'S_temporal_inf L', 'S_temporal_inf R', 'S_temporal_sup L', 'S_temporal_sup R', 'S_temporal_transverse L', 'S_temporal_transverse R'};
% Bands for ICC cortex plot
BandNames = {'theta', 'alpha', 'beta', 'gamma', 'highgamma'};
BandLowerFreqs = [ 4,  8, 13, 30,  50]; % Hz, inclusive
BandUpperFreqs = [ 8, 13, 30, 50, 150]; % Hz, exclusive

%% ===== VERIFY REQUIRED PROTOCOL =====
% Check Protocol that it exists
iProtocolOmega = bst_get('Protocol', ProtocolNameOmega);
if isempty(iProtocolOmega)
    error(['Unknown protocol: ' ProtocolNameOmega]);
end
% Select input protocol
gui_brainstorm('SetCurrentProtocol', iProtocolOmega);

% Verify the Subjects
ProtocolSubjects = bst_get('ProtocolSubjects');
ProtocolSubjectNames = {ProtocolSubjects.Subject.Name};
if ~all(ismember(SubjectNames, ProtocolSubjectNames))
    error(['All requested subjects must be present in the ' ProtocolNameOmega ' protocol.']);
end


%% ===== FIND FILES =====
bst_report('Start');

% Get raw and source files for each Subject
sRawDataFiles = [];
sSourcesFiles = [];

for iSubject = 1 : nSubjects
    % Process: Select data files in: sub-000*/*
    sFiles = bst_process('CallProcess', 'process_select_files_data', [], [], ...
        'subjectname',   SubjectNames{iSubject}, ...
        'condition',     '', ...
        'tag',           '', ...
        'includebad',    0, ...
        'includeintra',  0, ...
        'includecommon', 0);
    sRawDataFiles = [sRawDataFiles, sFiles];

    % Process: Select results files in: sub-000*/*
    sFiles = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   SubjectNames{iSubject}, ...
        'condition',     '', ...
        'tag',           '', ...
        'includebad',    0, ...
        'includeintra',  0, ...
        'includecommon', 0);
    sSourcesFiles = [sSourcesFiles, sFiles];
end


%% Compute PSD for all ROIs of an atlas
% For each subject, their recordings are split in two parts (data segments)
% these two data segments will be used to create PSDs
% These two PSDs are the features which will define the brain-fingerprint

nSegments = 2; % Training and Validation
timeIni   = zeros(nSubjects, nSegments);
timeFin   = zeros(nSubjects, nSegments);
sPsdFiles = repmat(db_template('processfile'), nSubjects, nSegments);
for iSubject = 1 : nSubjects
    sData = load(file_fullpath(sRawDataFiles(iSubject).FileName), 'Time');
    halfTime = diff(sData.Time) / 2;
    % First segment
    timeIni(iSubject, 1) = sData.Time(1) + 30;
    timeFin(iSubject, 1) = halfTime      - 30;
    % Second segment
    timeIni(iSubject, 2) = halfTime      + 30;
    timeFin(iSubject, 2) = sData.Time(2) - 30;
    % Compute PSD
    for iSegment = 1 : size(timeIni, 2)
        % Process: Power spectrum density (Welch)
        sPsdFiles(iSubject, iSegment) = bst_process('CallProcess', 'process_psd', sSourcesFiles(iSubject).FileName, [], ...
            'timewindow',  [timeIni(iSubject, iSegment), timeFin(iSubject, iSegment)], ...
            'win_length',  2, ...
            'win_overlap', 50, ...
            'units',       'physical', ...  % Physical: U2/Hz
            'clusters',    {AtlasName, AtlasScouts}, ...
            'scoutfunc',   1, ...  % Mean
            'win_std',     0, ...
            'edit',        struct(...
                 'Comment',         'Scouts,Power', ...
                 'TimeBands',       [], ...
                 'Freqs',           [], ...
                 'ClusterFuncTime', 'before', ...
                 'Measure',         'power', ...
                 'Output',          'all', ...
                 'SaveKernel',      0));
    end
end


%% ===== VECTORIZE PSD DATA FOR FINGERPRINTING =====
% Find size of requested PSD data
sPsdMat = in_bst_timefreq(sPsdFiles(1,1).FileName, 0, 'RowNames', 'Freqs');
Freqs = sPsdMat.Freqs;
ixLowerFreq = find(Freqs >= LowerFreq, 1, 'first');
ixUpperFreq = find(Freqs <  UpperFreq, 1, 'last');
Freqs = sPsdMat.Freqs(ixLowerFreq:ixUpperFreq);
nFreqs  = (ixUpperFreq - ixLowerFreq + 1);
ScoutNames = sPsdMat.RowNames;
nScouts = length(ScoutNames);

% Subject vectors
trainingVectors   = zeros(iSubject, nScouts * nFreqs);
validationVectors = zeros(iSubject, nScouts * nFreqs);

% Vectorize Scout PSDs
for iSubject = 1 : nSubjects
    for iSegment = 1 : nSegments
        sPsdMat = in_bst_timefreq(sPsdFiles(iSubject, iSegment).FileName, 0, 'TF');
        if iSegment == 1
            trainingVectors(iSubject,:) = reshape(squeeze(sPsdMat.TF(:, :, ixLowerFreq:ixUpperFreq)), 1, []);
        elseif iSegment == 2
            validationVectors(iSubject,:) = reshape(squeeze(sPsdMat.TF(:, :, ixLowerFreq:ixUpperFreq)), 1, []);
        end
    end
end


%% ==== SIMILARITY AND DIFFERENTIABILITY =====
% Subject similarity matrix
% It is generally symmetric although it does not necessarily have to be!
SubjectCorrMatrix= corr(trainingVectors', validationVectors');
% Differentiability
diff_1= (diag(SubjectCorrMatrix)-mean(SubjectCorrMatrix,1) )/ std(SubjectCorrMatrix,0,1);  % along columns
diff_2= (diag(SubjectCorrMatrix)-mean(SubjectCorrMatrix,2)')/ std(SubjectCorrMatrix,0,2)'; % along rows
% Differentiability across rows and columns are generally strongly correlated
% Take the mean for a summary statistic per participant
Differentiability = (diff_1 + diff_2)/2; % Mean differentiability derived from rows and columns

% IntrClass Correlations (ICC)
zTraining  = (trainingVectors - mean(trainingVectors, 2)) ./ std(trainingVectors, 0, 2);
zValidation= (validationVectors - mean(validationVectors, 2)) ./ std(validationVectors, 0, 2);
df_b = nSubjects-1;             % degrees of freedom for s2
df_w = nSubjects*(nSegments-1); % degrees of freedom for r
nFeatures = nScouts * nFreqs;
icc = zeros(1, nFeatures);
for iFeature = 1 : nFeatures
    x = [zTraining(:, iFeature), zValidation(:, iFeature)];
    x_w_mean = mean(x, 2);
    x_g_mean = mean(x(:));
    ss_t = sum((x - x_g_mean).^2, 'all');
    ss_w = sum((x - x_w_mean).^2, 'all');
    ss_b = ss_t - ss_w;
    ms_b = ss_b / df_b;
    ms_w = ss_w / df_w;
    icc(iFeature) = (ms_b - ms_w) ./ (ms_b + ((nSegments-1).*ms_w));
end
iccFreqsScouts = reshape(icc, [nFreqs, nScouts]);
% Compute Scout ICC for frequency bands
nBands = length(BandNames);
iccBands = zeros(nScouts, nBands);
for iBand = 1 : nBands
    ixLowerFreq = find(Freqs >= BandLowerFreqs(iBand), 1, 'first');
    ixUpperFreq = find(Freqs <  BandUpperFreqs(iBand), 1, 'last');
    iccBands(: ,iBand) = mean(iccFreqsScouts(ixLowerFreq:ixUpperFreq, :), 1)';
end

%% ===== SAVE OUTCOME IN BRAINSTORM DATABASE =====
% Retrieve 'Group_analysis' subject
sNormSubj = bst_get('NormalizedSubject');
% Study to save files
[sOutputStudy, iOutputStudy] = bst_get('StudyWithCondition', bst_fullfile(sNormSubj.Name, bst_get('DirAnalysisIntra')));

% Similarity matrix
sSimilarityMat = db_template('timefreqmat');
% Reshape: [nA x nB x nTime x nFreq] => [nA*nB x nTime x nFreq]
sSimilarityMat.TF = reshape(SubjectCorrMatrix, [], 1, 1);
sSimilarityMat.Comment      = 'Similarity matrix';
sSimilarityMat.DataType     = 'matrix';
sSimilarityMat.Time     = [0, 1];
sSimilarityMat.RefRowNames = cellfun(@(x) ['Train ', x], SubjectNames, 'UniformOutput', false);
sSimilarityMat.RowNames    = cellfun(@(x) ['Validation ', x], SubjectNames, 'UniformOutput', false);
% Output filename
SimilarityFile = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'timefreq_connectn_corr');
% Save file
bst_save(SimilarityFile, sSimilarityMat, 'v6');
% Add file to database structure
db_add_data(iOutputStudy, SimilarityFile, sSimilarityMat);

% Differentiability matrix
sDiffMat = db_template('matrixmat');
sDiffMat.Value   = Differentiability;
sDiffMat.Comment = 'Differentiability';
sDiffMat.Time     = [0, 1];
sDiffMat.Description = SubjectNames;
% Output filename
DiffFile = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'matrix');
% Save file
bst_save(DiffFile, sDiffMat, 'v6');
% Add file to database structure
db_add_data(iOutputStudy, DiffFile, sDiffMat);

% Save band ICC values (scout level)
for iBand = 1 : nBands
    sIccBandMat = db_template('matrixmat');
    sIccBandMat.Value   = iccBands(:, iBand);
    sIccBandMat.Comment = ['ICC_' BandNames{iBand}];
    sIccBandMat.Time    = [0, 1];
    sIccBandMat.Description = ScoutNames;
    % Output filename
    IccBandFile = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'matrix');
    % Save file
    bst_save(IccBandFile, sIccBandMat, 'v6');
    % Add file to database structure
    db_add_data(iOutputStudy, IccBandFile, sIccBandMat);
end

% Cortex surface display for band ICC values
% Verify destination surface
sSubjectGroup = bst_get('Subject', sOutputStudy.BrainStormSubject);
sSurfFile = sSubjectGroup.Surface(sSubjectGroup.iCortex);
sSurfMat = in_tess_bst(sSurfFile.FileName);
iAtlas = find(ismember({sSurfMat.Atlas.Name}, AtlasName));
if isempty(iAtlas)
    error('Default surface in "%s" does not contain "%s" atlas', bst_get('NormalizedSubjectName'), AtlasName);
end
sAtlas = sSurfMat.Atlas(iAtlas);
if ~all(ismember({sAtlas.Scouts.Label}, ScoutNames))
    error('Atlas "%s" in "%s" does not contain the same scouts as ICC', AtlasName, bst_get('NormalizedSubjectName'));
end
% Create one full results file per band ICC values
IccBandFiles = {};
for iBand = 1 : nBands
    sIccBandMat = db_template('resultsmat');
    sIccBandMat.SurfaceFile = sSurfFile.FileName;
    sIccBandMat.ImageGridAmp = nan(size(sSurfMat.Vertices, 1), 1);
    sIccBandMat.Comment = ['ICC_' BandNames{iBand}];
    sIccBandMat.Time    = [0, 1];
    % Assign ICC values to vertices in Scout
    for ix = 1 : nScouts
        iScout = find(ismember({sAtlas.Scouts.Label}, ScoutNames(ix)));
        sIccBandMat.ImageGridAmp(sAtlas.Scouts(iScout).Vertices) = deal(iccBands(ix, iBand));
    end
    % Output filename
    IccBandFile = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'results');
    IccBandFiles = [IccBandFiles, IccBandFile];
    % Save file
    bst_save(IccBandFile, sIccBandMat, 'v6');
    % Add file to database structure
    db_add_data(iOutputStudy, IccBandFile, sIccBandMat);
end

% Reload database
db_reload_studies(iOutputStudy)


%% ===== SNAPSHOTS =====
% Process: Snapshot: Similarity matrix
bst_process('CallProcess', 'process_snapshot', SimilarityFile, [], ...
    'type',    'connectimage', ...  % Connectivity matrix
    'Comment', 'Similarity matrix');

% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', DiffFile, [], ...
    'type',           'data', ...  % Recordings time series
    'time',           0, ...
    'Comment',        'Differentiability');

for iBand = 1 : nBands
    % Process: Snapshot: Sources (one time)
    bst_process('CallProcess', 'process_snapshot', IccBandFiles{iBand}, [], ...
        'type',           'sources', ...  % Sources (one time)
        'orient',         1, ...  % left
        'time',           0, ...
        'contact_time',   [0, 0.1], ...
        'contact_nimage', 12, ...
        'threshold',      0, ...
        'surfsmooth',     30, ...
        'mni',            [0, 0, 0], ...
        'Comment',        ['ICC_' BandNames{iBand}]);
end

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, reports_dir);
else
    bst_report('Open', ReportFile);
end

Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):