Connectivity

Authors: Hossein Shahabi, Raymundo Cassani, Takfarinas Medani, François Tadel, Sylvain Baillet

Brain functions (e.g., in cognition, behavior and perception) stem from the coordinated activity of multiple regions. Brain connectivity investigates how these different regions (or nodes) interact as a network. Depending on which connectivity characteristic is studied, a distinction is made between structural (fiber pathways), functional (non-directed statistical associations) and effective (causal interactions) connectivity between regions. Effective connectivity is often referred as directed functional connectivity. In this tutorial we will see how to compute different connectivity metrics for non-directed and directed functional analyses using Brainstorm, first with simulated data and later with real data.

We encourage the interested reader to learn more about the specific aspects of electrophysiology for studying human connectomics.

Contents

  1. Introduction
  2. Requirements
  3. Simulated data
  4. Correlation
  5. Coherence
  6. Granger causality
  7. Spectral Granger causality
  8. Envelope correlation
  9. Phase locking value
  10. Phase transfer entropy
  11. Method selection and comparison
  12. Scout-level connectivity
  13. Matrix thresholding
  14. On the hard drive
  15. Additional documentation
  16. Scripting
  17. TODO

Introduction

Definitions

Connectivity analyses are commonly performed by computing a bivariate measure between pairs of regional time series of interest. The outcome (a.k.a connectome) can be presented as a connectivity graph (left image), where each region is represented as a node (x, y, z,...), and the values of the connectivity metric shown next to the edge linking between two nodes. The connectome can also be represented by a connectivity matrix, a.k.a. adjacency matrix (right image).

cnx_graph_matrix.png

Sensors or sources: The signals used for connectivity analyses can be from sensor data (EEG/MEG signals) or from source time seriess (voxels or scouts).

Directed vs. non-directed: The direction of the interaction between signals (as statistical causation) can be measured with directed metrics. Non-directed metrics produce symmetrical connectivity graphs/matrices as connectivity "from Signal 

latex error! exitcode was 2 (signal 0), transscript follows:
to Signal 
latex error! exitcode was 2 (signal 0), transscript follows:
" is identical to connectivity "from Signal 
latex error! exitcode was 2 (signal 0), transscript follows:
to Signal 
latex error! exitcode was 2 (signal 0), transscript follows:
".

Experimental condition: Depending on the neuroscience question, connectivity analyses can be performed on resting-state (spontaneous) or task (e.g., trials) data.

Full (NxN) vs. seeded (1xN) connectivity: In a full connectivity analysis, the connectivity metric is computed for all the possible node pairs between N time series (noted N×N here). Alternatively, seeded connectivity (noted 1×N) is performed between one time series of interest (a seed, e.g., one brain region or a behavioral marker) and N other regions/time series.

Time-frequency transformations: Some connectivity metrics rely on a time-frequency representation of the signals. These latter are obtained with approaches such as the short-time Fourier transform, Hilbert transform, and Morlet wavelets.

Just as in other areas of electrophysiology studies, connectivity analyses need to be guided by mechanistic hypotheses concerning the expected effects.

Sensor-level

Sensor connectivity analyses present two important limitations:

  1. Their anatomical interpretation is limited and ambiguous.

  2. Sensor data is severely corrupted by field spread and volume conduction. Hence, activity from one single brain area is detected at multiple, often distant, sensor locations, which may be wrongly interpreted as network connections.

Source-level

Source connectivity analyses are neuroanatomically interpretable and can be derived across participants, following spatial normalization and registration.

It is recommended to verify that the outcomes of sensor and source connectivity analyses are compatible with one another (Lai et al., 2018).

Full brain connectomes

Whole-brain connectivity analyses at the typical resolution of cortical surfaces in Brainstorm involve thousands of source locations, making the N×N derivations impractical. For instance, 15000 cortical vertices would yield a connectivity matrix of 15000x15000x8 bytes = 1.6Gb. If unconstrained cortical sources are used and coherence is computed across 50 frequency bins, the memory allocation increases to 45000x45000x50x8 = 754 Gb per participant/condition/trial/etc. Reducing the N (via e.g., cortical parcellations, ROIs) is therefore essential.

Regions of interest

One solution for reducing the complexity of the connectivity analysis is to group the sources by ROIs, defined on the cortical surface or in the volume. ROIs can be defined based on study priors and other considerations such as the source estimation method, experimental task, and data available (Schhoffen and Gross, 2009):

The tutorial Corticomuscular coherence explains the computation of connectivity measures between one sensor and the minimum norm source maps: sensor x sources and sensor x scouts. The computation of ROI-based connectomes (scouts x scouts) is described in the section Scout-level connectivity of this tutorial page.

Diverse studies have shown some overlap between connectomes derived from electrophysiological signals (MEG/EEG) and the ones derived from fMRI, which is reasonably expected as both are the result of the undergoing biological system. However, due to its nature, the electrophysiological connectomes provide unique insights on how functional communication is implemented in the brain (Sadaghiani et al., 2022).

Whole-brain connectivity estimates may be exposed to the issue of circular analysis (Kriegeskorte et al., 2009).

Requirements

Here we skip most of the interface details and focus on the specifics on connectivity analyses with Brainstorm. So please make sure you are familiar with Brainstorm and go through all introduction tutorials first.

We first use simulated data to emphasize the theoretical aspects of each connectivity metric with respect to groundtruth outcomes. Real, empirical MEG data are featured later in the tutorial (same auditory oddball dataset as other tutorial sections).

Let's start by creating a new protocol in the Brainstorm database:

Simulated data

To compare connectivity metrics, let's use simulated time series with known ground truth interactions using a multivariate autoregressive (MVAR) model. The model we'll use consists of the following three signals:

We will simulate the fact that the component of Signal 3 at 25 Hz is driven in part by that of Signal 1 (denoted Signal 1>>Signal 3).

Let's now generate those three time series:

Process options:

Execution:

Credits:



Correlation

Correlation is a non-directed connectivity metric that can be used to show similarity, dependence or association among two random variables or signals. While this metric has been widely used in electrophysiology, it should not be considered the best technique to evaluate connectivity. Due to its nature, correlation fails to alleviate the problem of volume conduction and cannot explain the association in different frequency bands. However, it still can provide valuable information in case we deal with a few narrow-banded signals.

Process options

Result visualization

Display options:



Coherence

Coherency or complex coherence, 

latex error! exitcode was 2 (signal 0), transscript follows:
, is a complex-valued metric that measures the linear relationship of two signals in the frequency domain. Its magnitude square coherence (MSC), 
latex error! exitcode was 2 (signal 0), transscript follows:
, often referred to as coherence, measures the covariance of two signals in the frequency domain. For a pair of signals 
latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
, with spectra 
latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
, the MSC is defined as:

Two related measures, which alleviate the problem of volume conduction, are imaginary coherence (Nolte et al., 2004), 

latex error! exitcode was 2 (signal 0), transscript follows:
, and the lagged coherence (Pascual-Maqui, 2007), 
latex error! exitcode was 2 (signal 0), transscript follows:
, which are defined as:

where 

latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
describe the imaginary and real parts of a complex number.

To calculate coherence values in Brainstorm, select the process.

Process options

Result visualization

Coherence is a function of frequency, as such, for each frequency point there is a connectivity graph and a connectivity matrix. Right-click on the coherence result file to see its display options:

Open the 3 representations. These representations are linked such as by clicking on the spectral representation of the coherence, we change the frequency that is displayed in the connectivity graph and matrix. This frequency can be also changed in the Time panel.

res_cohere1n_a.png

res_cohere1n_a2.png

res_cohere1n_b.png

res_cohere1n_c.png

In the same way, we can compute the other types of coherence. The figure below presents the spectra for the imaginary coherence (left) and the lagged coherence (right). Both, imaginary and lagged coherence aim to address the volume conduction problem, although they present small differences.

res_cohere1n_d.png

res_cohere1n_e.png



Granger causality

Granger causality (GC) is a method of directed functional connectivity, which is base on the Wiener-Granger causality methodology. GC is a measure of linear dependence, which tests whether the prediction of signal 

latex error! exitcode was 2 (signal 0), transscript follows:
(using a linear autoregressive model) is improved by adding signal 
latex error! exitcode was 2 (signal 0), transscript follows:
(also using a linear autoregressive model). If this is true, signal 
latex error! exitcode was 2 (signal 0), transscript follows:
has a Granger causal effect on the first signal. In other words, independent information of the past of signal 
latex error! exitcode was 2 (signal 0), transscript follows:
improves the prediction of signal 
latex error! exitcode was 2 (signal 0), transscript follows:
obtained with the past of signal 
latex error! exitcode was 2 (signal 0), transscript follows:
alone. GC is nonnegative, and zero when there is no Granger causality. As only the past of the signals is considered, the GC metric is directional. The term independent is emphasized because it creates some interesting properties for GC, such as, that it's invariant under rescaling of 
latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
, as well as the addition of a multiple of 
latex error! exitcode was 2 (signal 0), transscript follows:
to 
latex error! exitcode was 2 (signal 0), transscript follows:
.
See Granger causality - mathematical background for a complete formulation of the method.

Despite the name, Granger causality indicates directionality but not true causality.
For example, if a variable 

latex error! exitcode was 2 (signal 0), transscript follows:
is causing both 
latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
, but with a smaller delay for 
latex error! exitcode was 2 (signal 0), transscript follows:
than for 
latex error! exitcode was 2 (signal 0), transscript follows:
, then the GC measure between 
latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
would show a non-zero GC for 
latex error! exitcode was 2 (signal 0), transscript follows:
--> 
latex error! exitcode was 2 (signal 0), transscript follows:
, even though 
latex error! exitcode was 2 (signal 0), transscript follows:
is not truly causing 
latex error! exitcode was 2 (signal 0), transscript follows:
(Bressler and Seth, 2011).

Process options

Result visualization

In the connectivity graph (left) the directionality is shown with an arrow head at the center for the arc connecting nodes. As GC metric is not symmetric, the connectivity matrix (right) is not symmetric. The upper right element of this matrix shows there is a signal flow from signal 1 to signal 3.

res_granger1n_a.png

res_granger1n_b.png



Spectral Granger causality

GC lacks of resolution in the frequency domain, as such, the spectral Granger causality was developed (Dhamala et al., 2008).

Process options

Result visualization

As with coherence, spectral GC can be plotted as a function of frequency. The plot below clearly shows a peak around 25 Hz for the interaction from signal 1 to signal 3, as expected.

res_spgranger1n.png

Envelope correlation

In the time-frequency tutorial the Morlet wavelets and Hilbert transform were introduced as methods to decompose signals in the time-frequency (TF) domain. The result of this TF transformation can be seen as a set of narrowband complex signals, which are analytic signals.

The analytic signal, 

latex error! exitcode was 2 (signal 0), transscript follows:
, is a complex signal uniquely associated to a real signal, 
latex error! exitcode was 2 (signal 0), transscript follows:
, that has been useful in signal processing due to its characteristics, more specifically, its module 
latex error! exitcode was 2 (signal 0), transscript follows:
, and phase 
latex error! exitcode was 2 (signal 0), transscript follows:
, correspond to the instantaneous amplitude (or envelope) and instantaneous phase of the associated real signal 
latex error! exitcode was 2 (signal 0), transscript follows:
. The real part of 
latex error! exitcode was 2 (signal 0), transscript follows:
is its associated real signal 
latex error! exitcode was 2 (signal 0), transscript follows:
, and the imaginary part is the Hilbert transform of the same real signal 
latex error! exitcode was 2 (signal 0), transscript follows:
.

The analytic signal of oscillatory or narrowband signals provide meaningful and interpretable results for the instantaneous amplitude and phase. While it could be computed for broadband signals, the instantaneous parameters would be difficult to interpret, as they would be the contributions of several oscillatory signals (Cohen, 2014).

The instantaneous amplitude (or envelope) of these band analytic signals can be used to carry out pairwise connectivity analysis with metrics such as correlation and coherence (including lagged coherence).

In computing the envelope correlation, an optional step is to orthogonalize the envelopes by removing their real part of coherence before the correlation (Hipp et al., 2012). This orthogonalization process alleviates the effect of volume conduction in MEG/EEG signals.

Process options

Result visualization

Similar to the results from coherence and spectral Granger causality, the envelope correlation can be plotted as a function of frequency, and as a function of time if the Time resolution option is set to Dynamic. Below, the results obtained with the Hilbert transform (left) and with Morlet wavelet (right) for the first 5-s window (top) and the 5-s last window (bottom).

res_henv1n_h.png

First 5-s window

res_henv1n_w.png

res_henv1n_h2.png

Last 5-s window

res_henv1n_w2.png

Phase locking value

An alternative class of connectivity metrics considers only the relative instantaneous phase between the two signals, i.e., phase-locking or synchronization (Tass et al., 1998). Phase locking is a fundamental concept in dynamical systems that has been used in control systems (the phase-locked loop) and in the analysis of nonlinear, chaotic and non-stationary systems. Since the brain is a nonlinear dynamical system, phase locking is an appropriate approach to quantifying connectivity. A more pragmatic argument for its use in studies of LFPs, EEG, and MEG is that it is robust to fluctuations in amplitude that may contain less information about interactions than does the relative phase (Lachaux et al., 1999; Mormann et al., 2000).

The most commonly used phase connectivity metric is the phase-locking value (PLV), which is defined as the length of the average vector of many unit vectors whose phase angle corresponds to the phase difference between two signals (Tass et al., 1998). If the distribution of the phase difference between the two signals is uniform, the length of such an average vector will be zero. Conversely, if the phases of the two signals are strongly coupled, the length of the average vector will approach unity. For event-related studies, we would expect the phase difference across trials to be uniform distributed unless the phase is locked to the stimulus. In that case, we may have nonuniform marginals which could in principle lead to indications of phase locking between two signals. Considering a pair of narrow-band analytic signals 

latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
, obtained from the TF transformation using the Hilbert transform: 

latex error! exitcode was 2 (signal 0), transscript follows:

with: 

latex error! exitcode was 2 (signal 0), transscript follows:

plv.png

Other measures have been added rencently thanks to the contribution of Daniele Marinazzo:
ciPLV (Bruña 2018) and wPLI (Vinck 2011).

PLV value tend to be overestimated when there are few (<50) phase samples for its computation. As a consequence, when comparing PLV values between conditions, make sure that these have the approximately the same number of samples.

Process options

Result visualization

PLV is frequency resolved, and it was computed for the delta, theta, alpha, beta and gamma bands. With the simulated data, we expect a higher PLV value in the beta band (15 to 29 Hz), between signal 1 and signal 3. This result is seen as a peak at 22 Hz (center of beta band) shown in PLV as a function of frequency.

res_plv1n.png

Phase transfer entropy

Phase transfer entropy (PTE) is a directed connectivity metric that quantifies the transfer entropy (TE) between two instantaneous phase time series (Lobier et al., 2014). Similar to GC, TE estimates whether including the past of both source and target time-series influences the ability to predict the future of the target time-series. In PTE, if a phase signal 

latex error! exitcode was 2 (signal 0), transscript follows:
causes the signal 
latex error! exitcode was 2 (signal 0), transscript follows:
, the mutual information, between 
latex error! exitcode was 2 (signal 0), transscript follows:
and the past of 
latex error! exitcode was 2 (signal 0), transscript follows:
i.e. 
latex error! exitcode was 2 (signal 0), transscript follows:
is larger than the mutual information of 
latex error! exitcode was 2 (signal 0), transscript follows:
, the past of 
latex error! exitcode was 2 (signal 0), transscript follows:
i.e. 
latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
. This relationship can be seen on the Venn diagram below, where 
latex error! exitcode was 2 (signal 0), transscript follows:
and 
latex error! exitcode was 2 (signal 0), transscript follows:
, indicate mutual information and the individual entropies respectively. Lastly, PTE cannot be negative, and its magnitude does not have a meaningful upper bound.

pte.png

Process options

Result visualization

PTE was computed for the delta, theta, alpha, beta and gamma bands. With the simulated data, we expect a higher PTE value in the beta band, From signal 1 To signal 3, as PTE is directed metric. This is confirmed with a peak at 22 Hz (center of beta band) shown in PTE frequency representation.

res_pte1n.png .

Method selection and comparison

The following table list the available connectivity metrics in Brainstorm and their description.

Metric

Directionality

Domain

1×N

N×N

Time resolved

Process

Info

Correlation

Non-directed

Time

✅

✅

✅

bst_corrn.m

Link

Coherence

Non-directed

Frequency

✅

✅

✅

bst_cohn.m

Link

Granger causality

Directed

Time

✅

✅

❌

bst_granger.m

Link

Spectral Granger causality

Directed

Frequency

✅

✅

❌

bst_granger_spectral.m

Link

Envelope Correlation (2020)

Non-directed

T-F

✅

✅

✅

bst_henv.m

Link

Phase locking value

Non-directed

Phase

✅

✅

❌

bst_connectivity.m

Link

Phase transfer entropy

Directed

Phase

❌

✅

❌

PhaseTE_MF.m

Link

Scout-level connectivity

The sections above explain the computation of various connectivity measures between a few pairs of signals. When computing whole-brain connectomes (i.e. NxN connectivity matrices between all the brain sources), there are additional technical questions to take into account: the number of signals is typically too large for keeping the full resolution of the minimum norm source maps, as explained in the introduction of this page, and the unconstrained source maps with 3 signals at each location require an extra step of simplification.

The tutorial Corticomuscular coherence explains the computation of connectivity measures between one sensor and the source maps: sensor x sources and sensor x scouts, both in the constrained and unconstrained cases.

This section explains the dimension reduction using ROIs (named "scouts" in Brainstorm), using the coherence NxN (scouts x scouts) as an example, in the case of unconstrained source maps (three orthogonoal orientations at each location).

Scout function: In this configuration, one connectivity result (i.e. coherence spectrum) is computed for each pair or scouts in each orientation. It is necessary to provide two parameters that define how the data is aggregated per scout: The scout function (mean is often used), and when the within-scout aggregation takes place (before or after the coherence computation).

Unconstrained maps / maximum: The graphs below show two scouts (Scout-1 and Scout-2), with 3 orientations each (x,y,z). The coherence is computed for each pair of orientations (1x-2x, 1x-2y, 1x-2z, 1y-2x, ..., 1z-2z), leading to 9 coherence spectra. From these 9 values at each frequency bin, only the maximum value is kept to represent the connectivity between the two scouts, for a final output of one coherence spectrum per pair of scouts. The choice of the maximum statistic is empirical: it lacks rotational invariance (if you change the position of the NAS/LPA/RPA fiducials, you get different x,y,z axes and therefore different results), but among the solutions we tested, it is the one that led to the smoothest and most reproducible spatial maps.

Before: The scout function is applied for each direction on the vertices' source time series that make up a scout; resulting in one time series per direction per scout. Then, the scouts time series are used to compute coherence, and the coherence spectra are aggregated across dimensions, to obtain one coherence spectrum per scout. https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_nxn_coh_sct_bef.png

After: Coherence is computed between each pair of dipoles (number of vertices x 3 orientations). Then, the scout function is applied on the coherence spectra for each direction of the vertices within a scout. Finally these spectra are aggregated across dimensions to obtain a coherence spectrum per scout. This option computes the coherence between 45000x45000 source signals, instead of a handful of times with the "before" option. The computation is therefore much longer and demanding in terms of RAM memory. See the introduction for an example.https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_nxn_coh_sct_aft.png


The result is a NxN connectivity file (), It contains (Scouts x Scouts) coherence spectra.

diagram_nxn_coh_sct_end_small.gif

Such a visualization is not practical, thus the connectivity graph or the adjacent matrix are displayed for each frequency bin in the coherence spectra. For more details, see the connectivity graph tutorial.

Matrix thresholding

The connectivity graphs obtained with the methods described above indicate non-zero connectivity values for most connections. Most of these values are not significant and should be excluded from any further graph analysis or report.

Even though the Brainstorm interface provides tools to apply a fixed threshold to the connectivity matrices, based on an arbitrary value of the connectivity metric or a percentage of the strongest connections, we recommend assessing the significance of a connection using permutation tests across multiple participants.

In the context of a group analysis: compute the same connectivity measure for each subject, for two different experimental conditions, or between an active state and a baseline. Then run a non-parametric paired permutation test to compare the two conditions across subjects. More information in the Statistics tutorial.

At the moment, we do not provide any solution for within-subject statistical thresholding of the connectivity matrices, i.e. testing using multiple trials within the same subject.

On the hard drive

File name

The connectivity file structure is an extension of the time-frequency structure. The file names start with timefreq_, followed by connect1 (1xN) or connectn (NxN or AxB), the connectivity method and a time stamp. Example: timefreq_connectn_corr_220120_1350.mat.

File structure

Right click on of the first connectivity file computed here > File > View file contents.

The data structure is the same as for the time-frequency files. Only the fields that some specificity related with the connectivity analysis are documented here.

Connectivity matrix encoding

Let's consider the structure TfMat, loaded from a connectivity file, for example by right-clicking on the file > File > Export to Matlab. The connectivity matrix R can be obtained with function GetConnectMatrix. The size of R is [Na x Nb x Ntime x Nfreq]. 

R = bst_memory('GetConnectMatrix', TfMat);

If the matrix is not symmetrical (ie. not compressed): for each time and frequency, the list of values from the first dimension of the TF variable are reshaped into a 2D matrix: [Na x Nb].

If the matrix is symmetrical and compressed: only the lower triangular matrix is saved in the TF variable. The full matrix is first reconstructed with function process_compress_sym>Expand, then reshaped into [Na x Nb].

Saving the connectivity matrix R back in the TF variable is possible by reshaping to [Nx1] (R(:)) and then compressing again the matrix: 

TfMat.TF = process_compress_sym('Compress', R(:));

Additional documentation

Articles

Forum discussions

Scripting

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

1 function tutorial_connectivity(reports_dir) 2 % TUTORIAL_CONNECTIVITY: Script that runs the Brainstorm connectivity tutorial. 3 % 4 % INPUTS: 5 % - reports_dir : Directory where to save the execution report (instead of displaying it) 6 7 % @============================================================================= 8 % This function is part of the Brainstorm software: 9 % https://neuroimage.usc.edu/brainstorm 10 % 11 % Copyright (c) University of Southern California & McGill University 12 % This software is distributed under the terms of the GNU General Public License 13 % as published by the Free Software Foundation. Further details on the GPLv3 14 % license can be found at http://www.gnu.org/copyleft/gpl.html. 15 % 16 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 17 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 18 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 19 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 20 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 21 % 22 % For more information type "brainstorm license" at command prompt. 23 % =============================================================================@ 24 % 25 % Author: Raymundo Cassani, 2021-2022 26 % Francois Tadel, 2022 27 28 29 %% ===== PARAMETERS ===== 30 % Output folder for reports 31 if (nargin < 1) || isempty(reports_dir) || ~isfolder(reports_dir) 32 reports_dir = []; 33 end 34 35 36 %% ===== CREATE PROTOCOL ===== 37 % Start brainstorm without the GUI 38 if ~brainstorm('status') 39 brainstorm nogui 40 end 41 % Create Protocol 42 ProtocolName = 'TutorialConnectivity'; 43 % Delete existing protocol 44 gui_brainstorm('DeleteProtocol', ProtocolName); 45 % Create new protocol 46 gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); 47 % Start a new report 48 bst_report('Start'); 49 50 51 %% ===== SIMULATE DATA (MVAR MODEL) ===== 52 % Seed for random number generator 53 rng(111); 54 % Process: Simulate AR signals 55 sFileSim = bst_process('CallProcess', 'process_simulate_ar_spectra', [], [], ... 56 'subjectname', 'Subject01', ... 57 'condition', 'Simulation', ... 58 'samples', 12000, ... 59 'srate', 120, ... 60 'interactions', ['1, 1 / 10, 25 / 0.3, 0.5' 10 ... 61 '2, 2 / 10, 25 / 0.7, 0.3' 10 ... 62 '3, 3 / 10, 25 / 0.2, 0.2' 10 ... 63 '1, 3 / 25 / 0.1']); 64 65 % Make sure the display mode is "columns" 66 bst_set('TSDisplayMode', 'column'); 67 % Process: Snapshot: Recordings time series 68 bst_process('CallProcess', 'process_snapshot', sFileSim, [], ... 69 'type', 'data', ... % Recordings time series 70 'Comment', 'Simulated signals'); 71 72 73 74 %% ===== CORRELATION ===== 75 % Process: Correlation NxN 76 sFiles = bst_process('CallProcess', 'process_corr1n', sFileSim, [], ... 77 'timewindow', [], ... 78 'scalarprod', 0, ... 79 'outputmode', 1); % Save individual results (one file per input file) 80 81 % Process: Snapshot: Connectivity matrix 82 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 83 'type', 'connectimage', ... % Connectivity matrix 84 'Comment', 'Correlation NxN'); 85 86 87 %% ===== COHERENCE ===== 88 % Process: Magnitude-squared coherence: |C|^2 = |Cxy|^2/(Cxx*Cyy) 89 sFiles = bst_process('CallProcess', 'process_cohere1n', sFileSim, [], ... 90 'timewindow', [], ... 91 'removeevoked', 0, ... 92 'cohmeasure', 'mscohere', ... % Magnitude-squared coherence: |C|^2 = |Cxy|^2/(Cxx*Cyy) 93 'tfmeasure', 'stft', ... % Fourier transform 94 'tfedit', struct(... 95 'Comment', 'Complex', ... 96 'TimeBands', [], ... 97 'Freqs', [], ... 98 'StftWinLen', 1, ... 99 'StftWinOvr', 50, ... 100 'StftFrqMax', 60, ... 101 'ClusterFuncTime', 'none', ... 102 'Measure', 'none', ... 103 'Output', 'all', ... 104 'SaveKernel', 0), ... 105 'timeres', 'none', ... % None 106 'avgwinlength', 1, ... 107 'avgwinoverlap', 50, ... 108 'outputmode', 'input'); % separately for each file 109 110 % Process: Snapshot: Frequency spectrum 111 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 112 'type', 'spectrum', ... % Frequency spectrum 113 'Comment', 'MSC NxN'); 114 115 % Process: Snapshot: Connectivity graph 116 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 117 'type', 'connectgraph', ... % Connectivity graph 118 'Comment', 'MSC NxN'); 119 120 % Process: Snapshot: Connectivity matrix 121 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 122 'type', 'connectimage', ... % Connectivity matrix 123 'Comment', 'MSC NxN'); 124 125 % Process: Imaginary coherence: IC = |imag(C)| 126 sFiles = bst_process('CallProcess', 'process_cohere1n', sFileSim, [], ... 127 'timewindow', [], ... 128 'removeevoked', 0, ... 129 'cohmeasure', 'icohere2019', ... % Imaginary coherence: IC = |imag(C)| 130 'tfmeasure', 'stft', ... % Fourier transform 131 'tfedit', struct(... 132 'Comment', 'Complex', ... 133 'TimeBands', [], ... 134 'Freqs', [], ... 135 'StftWinLen', 1, ... 136 'StftWinOvr', 50, ... 137 'StftFrqMax', 60, ... 138 'ClusterFuncTime', 'none', ... 139 'Measure', 'none', ... 140 'Output', 'all', ... 141 'SaveKernel', 0), ... 142 'timeres', 'none', ... % None 143 'avgwinlength', 1, ... 144 'avgwinoverlap', 50, ... 145 'outputmode', 'input'); % separately for each file 146 147 % Process: Snapshot: Frequency spectrum 148 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 149 'type', 'spectrum', ... % Frequency spectrum 150 'Comment', 'Imaginary coherence NxN'); 151 152 % Process: Lagged coherence / Corrected imaginary coherence: LC = |imag(C)|/sqrt(1-real(C)^2) 153 sFiles = bst_process('CallProcess', 'process_cohere1n', sFileSim, [], ... 154 'timewindow', [], ... 155 'removeevoked', 0, ... 156 'cohmeasure', 'lcohere2019', ... % Lagged coherence / Corrected imaginary coherence: LC = |imag(C)|/sqrt(1-real(C)^2) 157 'tfmeasure', 'stft', ... % Fourier transform 158 'tfedit', struct(... 159 'Comment', 'Complex', ... 160 'TimeBands', [], ... 161 'Freqs', [], ... 162 'StftWinLen', 1, ... 163 'StftWinOvr', 50, ... 164 'StftFrqMax', 60, ... 165 'ClusterFuncTime', 'none', ... 166 'Measure', 'none', ... 167 'Output', 'all', ... 168 'SaveKernel', 0), ... 169 'timeres', 'none', ... % None 170 'avgwinlength', 1, ... 171 'avgwinoverlap', 50, ... 172 'outputmode', 'input'); % separately for each file 173 174 % Process: Snapshot: Frequency spectrum 175 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 176 'type', 'spectrum', ... % Frequency spectrum 177 'Comment', 'Lagged coherence NxN'); 178 179 180 %% ===== GRANGER CAUSALITY ===== 181 % Process: Bivariate Granger causality NxN 182 sFiles = bst_process('CallProcess', 'process_granger1n', sFileSim, [], ... 183 'timewindow', [], ... 184 'removeevoked', 0, ... 185 'grangermethod', 'mvgc', ... 186 'grangerorder', 6, ... 187 'outputmode', 1); % Save individual results (one file per input file) 188 189 % Process: Snapshot: Connectivity graph 190 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 191 'type', 'connectgraph', ... % Connectivity graph 192 'Comment', 'Granger causality NxN'); 193 194 % Process: Snapshot: Connectivity matrix 195 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 196 'type', 'connectimage', ... % Connectivity matrix 197 'Comment', 'Granger causality NxN'); 198 199 200 %% ===== SPECTRAL GRANGER CAUSALITY ===== 201 % Process: Bivariate Granger causality (spectral) NxN 202 sFiles = bst_process('CallProcess', 'process_spgranger1n', sFileSim, [], ... 203 'timewindow', [], ... 204 'removeevoked', 0, ... 205 'grangermethod', 'mvgc', ... 206 'grangerorder', 6, ... 207 'maxfreqres', 1, ... 208 'maxfreq', 60, ... 209 'outputmode', 1); % Save individual results (one file per input file) 210 211 % Process: Snapshot: Frequency spectrum 212 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 213 'type', 'spectrum', ... % Frequency spectrum 214 'Comment', 'Spectral Granger causality NxN'); 215 216 217 %% ===== ENVELOPE CORRELATION ===== 218 % Process: Envelope Correlation NxN [2023] 219 sFiles = bst_process('CallProcess', 'process_henv1n', sFileSim, [], ... 220 'timewindow', [], ... 221 'removeevoked', 0, ... 222 'cohmeasure', 'oenv', ... % Envelope correlation (orthogonalized) 223 'tfmeasure', 'hilbert', ... % Hilbert transform 224 'tfedit', struct(... 225 'Comment', 'Complex', ... 226 'TimeBands', [], ... 227 'Freqs', {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'}}, ... 228 'ClusterFuncTime', 'none', ... 229 'Measure', 'none', ... 230 'Output', 'all', ... 231 'SaveKernel', 0), ... 232 'timeres', 'windowed', ... % Windowed 233 'avgwinlength', 5, ... 234 'avgwinoverlap', 50, ... 235 'parallel', 0, ... 236 'outputmode', 'input'); % separately for each file 237 238 % Process: Snapshot: Frequency spectrum 239 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 240 'type', 'spectrum', ... % Frequency spectrum 241 'Comment', 'Envelope correlation (Hilbert transform) NxN'); 242 243 % Process: Envelope Correlation NxN [2023] 244 sFiles = bst_process('CallProcess', 'process_henv1n', sFileSim, [], ... 245 'timewindow', [], ... 246 'removeevoked', 0, ... 247 'cohmeasure', 'oenv', ... % Envelope correlation (orthogonalized) 248 'tfmeasure', 'morlet', ... % Morlet wavelets 249 'tfedit', struct(... 250 'Comment', 'Complex,1-60Hz', ... 251 'TimeBands', [], ... 252 'Freqs', [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], ... 253 'MorletFc', 1, ... 254 'MorletFwhmTc', 3, ... 255 'ClusterFuncTime', 'none', ... 256 'Measure', 'none', ... 257 'Output', 'all', ... 258 'SaveKernel', 0), ... 259 'timeres', 'windowed', ... % Windowed 260 'avgwinlength', 5, ... 261 'avgwinoverlap', 50, ... 262 'parallel', 0, ... 263 'outputmode', 'input'); % separately for each file 264 265 % Process: Snapshot: Frequency spectrum 266 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 267 'type', 'spectrum', ... % Frequency spectrum 268 'Comment', 'Envelope correlation (Morlet wavelet) NxN'); 269 270 271 %% ===== PHASE LOCKING VALUE ===== 272 plv_variants = {'plv', ... % Phase locking value 273 'ciplv', ... % Lagged phase synchronization / Corrected imaginary PLV 274 'wpli'}; % Weighted phase lag index 275 for ix = 1 : length(plv_variants) 276 % Process: Phase locking value 277 sFiles = bst_process('CallProcess', 'process_plv1n', sFileSim, [], ... 278 'timewindow', [], ... 279 'plvmethod', plv_variants{ix}, ... 280 'plvmeasure', 2, ... % Magnitude 281 'tfmeasure', 'hilbert', ... % Hilbert transform 282 'tfedit', struct(... 283 'Comment', 'Complex', ... 284 'TimeBands', [], ... 285 'Freqs', {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'}}, ... 286 'ClusterFuncTime', 'none', ... 287 'Measure', 'none', ... 288 'Output', 'all', ... 289 'SaveKernel', 0), ... 290 'timeres', 'none', ... % None 291 'avgwinlength', 1, ... 292 'avgwinoverlap', 50, ... 293 'outputmode', 'input'); % separately for each file 294 295 % Process: Snapshot: Frequency spectrum 296 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 297 'type', 'spectrum', ... % Frequency spectrum 298 'Comment', ['Phase locking value (' plv_variants{ix} ') NxN']); 299 end 300 301 %% ===== PHASE TRANSFER ENTROPY ===== 302 % Process: Phase Transfer Entropy NxN 303 sFiles = bst_process('CallProcess', 'process_pte1n', sFileSim, [], ... 304 'timewindow', [], ... 305 'freqbands', {'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'}, ... 306 'normalized', 0, ... 307 'outputmode', 1); % Save individual results (one file per input file) 308 309 % Process: Snapshot: Frequency spectrum 310 bst_process('CallProcess', 'process_snapshot', sFiles, [], ... 311 'type', 'spectrum', ... % Frequency spectrum 312 'Comment', 'Phase transfer entropy NxN'); 313 314 315 %% ===== SAVE REPORT ===== 316 % Save and display report 317 ReportFile = bst_report('Save', []); 318 if ~isempty(reports_dir) && ~isempty(ReportFile) 319 bst_report('Export', ReportFile, reports_dir); 320 else 321 bst_report('Open', ReportFile); 322 end




Feedback on the documentation (typos, unclear sections, missing information)
For questions, bug reports, and feature requests, please use the Brainstorm Forum.
Email address (if you expect an answer):


TODO

Hossein, Richard

Raymundo, Sylvain

Francois

Tutorials/Connectivity (last edited 2022-05-24 12:04:42 by FrancoisTadel)