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

Introduction
Requirements
Simulated data
Correlation
Coherence
Granger causality
Spectral Granger causality
Envelope correlation
Phase locking value
Phase transfer entropy
Method selection and comparison
Scout-level connectivity
Matrix thresholding
On the hard drive
Additional documentation
Scripting
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).

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 $$x$$ to Signal $$y$$ " is identical to connectivity "from Signal $$y$$ to Signal $$x$$ ".

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:

Their anatomical interpretation is limited and ambiguous.
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):

Study priors from the literature, working hypotheses,
Association with a non-neuronal signal (e.g., cortico-muscular coherence),
Signal strength, ROIs as regions with strongest activity in an experimental condition,
Whole-brain cortical parcellation to reduce N.

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:

Select the menu File > Create new protocol > type in "TutorialConnectivity" and select the options:
- Yes, use protocol's default anatomy,
- No, use one channel file per condition.
Right-click on the TutorialConnectivity folder > New subject > Subject01

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:

Signal 1: two oscillatory components. One at 10 Hz (alpha band), and a stronger peak at 25 Hz (beta band).
Signal 2: Same as Signal 1, with strongest peak at 10 Hz.
Signal 3: Same as above, with the two components (10 and 25 Hz) of same magnitude.

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:

In the Process1 tab, leave the file list empty and click on the button [Run]
Select process: Simulate > Simulate AR signals.

Process options:

Subject name: Target subject for the simulated signals. Select Subject01.
Condition name: Target folder for the simulated signals. Set to Simulation.
Number of time samples: Duration of signals, in samples. Set to 12 000.
Sampling frequency: Fs for the simulated signals. Set to 120 Hz.
Interaction specifications: Spectral parameters for the signal components and their interactions in the MVAR model: From, To / Peak frequencies [Hz] / Peak relative magnitudes [0-1]
Set to:
```
1, 1 / 10, 25 / 0.3, 0.5
2, 2 / 10, 25 / 0.7, 0.3
3, 3 / 10, 25 / 0.2, 0.2
1, 3 / 25     / 0.1
```
Display the groundtruth spectral metrics of the MVAR-generated time series: transfer function, cross-spectral power density, magnitude square coherence, directed transfer function (DTF) and partial directed coherence (PDC). The transfer function () characterizes the relationships between signals in the frequency domain. It is a non-symmetric representation, which enables the identification of causal dependencies between signal components. The auto-transfer functions (shown in the graphs along the diagonal below) display the power spectra of each signal. The off-diagonal representations display the interactions between each pair of signals (see Signal 1 >> Signal 3 above). Here, we see the transfer function from signal 1 to signal 3. These transfer functions are our ground truth for connectivity values.
Get coefficients matrix: Shows the coefficients related to the MVAR model. These coefficients can be used in the process Simulate > Simulate AR signals (ARfit) to simulate the same model.

Execution:

Click Run to simulate the signals described by the MVAR model.
In the next sections we will compute different connectivity metrics for these simulated signals. As such, place the simulated data in the Process1 tab, click on [Run] ( ) to open the Pipeline editor, and select the connectivity metric.
This tutorial illustrates only the computation full connectivity graphs (Process 1: NxN). It is however possible to compute the connectivity between 1 signal and the other signals in the same file (Process1: 1xN), or between signals from two different files (Process2: AxB). These options will be illustrated at the end of the tutorial.

Credits:

This process relies on the ARSIM function from the ARFit toolbox:
https://github.com/tapios/arfit
Neumaier A, Schneider T
Estimation of parameters and eigenmodes of multivariate autoregressive models
ACM Transactions on Mathematical Software, 2001
Schneider T, Neumaier A
Algorithm 808: ARfit – A Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models
ACM Transactions on Mathematical Software, 2001

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

Process Connectivity > Correlation NxN:
Time window: Segment of the signal used for the connectivity analysis. Select: All file.
Compute scalar product: If unchecked, the mean of the signals is subtracted before computing the correlation. Uncheck it.
Output options: Select Save individual results.

Result visualization

The results are stored as a N×N connectivity file, icon . Right-click to see its display options:
Display as graph: Plots the connectivity graph using a chord diagram where the color of the edges shows the connectivity metric value. See the connectivity graph tutorial for a detailed explanation of the options of this visualization.
Display as image: Plots the adjacency matrix for the connectivity file.
Display fibers:Additional option available for source connectivity results when a fiber track surface is available, as shown here.

Display options:

Diagonal values: The value of the connectivity metric between a signal and itself is plotted as zero so that it doesn't force scaling the colormap to 1 if the other values are much smaller.
Labels: Click on the figure to see the signal names and connectivity values as the image legend. In order to see the labels corresponding to each column and row: right-click on the figure > Figure > Show labels. If your signal names are very long, try the option Use short labels.
Colormap: By default, the NxN colormap is configured to display the absolute values of the connectivity measures. As in the general case, correlation values can be positive or negative, you need to check this option carefully. If you expect to see negative values, make sure to change the colormap configuration: right-click on the figure > Colormap > Uncheck Absolute values.

Coherence

Coherency or complex coherence, $C_{xy}(f)$ , is a complex-valued metric that measures the linear relationship of two signals in the frequency domain. Its magnitude square coherence (MSC), $|C_{xy}(f)|^2$ , often referred to as coherence, measures the covariance of two signals in the frequency domain. For a pair of signals $$x(t)$$ and $$y(t)$$ , with spectra $$X(f)$$ and $$Y(f)$$ , the MSC is defined as:

$\begin{eqnarray*} C_{xy}(f) &=& \frac{S_{xy}(f)}{\sqrt{ S_{xx}(f)S_{yy}(f)}}\\ |C_{xy}(f)|^2 &=& MSC(f) = \left(\frac{\left |S_{xy}(f) \right |}{\sqrt{ S_{xx}(f)S_{yy}(f) }}\right)^2 = \frac{\left |X(f)Y^*(f) \right |^{2}}{X(f)X^*(f)Y(f)Y^*(f)} \\ S_{xy}(f) &:& \textrm{Cross-spectrum} \\ S_{xx}(f) \quad \textrm{and} \quad S_{yy}(f) &:& \textrm{Auto-spectra or power spectral densities} \\ \end{eqnarray*}$

Two related measures, which alleviate the problem of volume conduction, are imaginary coherence (Nolte et al., 2004), $IC_{xy}(f)$ , and the lagged coherence (Pascual-Maqui, 2007), $LC_{xy}(f)$ , which are defined as:

$\begin{eqnarray*} IC_{xy}(f) &=& \mathrm{Im} \left (C_{xy}(f) \right ) = \frac{\mathrm{Im} \left (S_{xy}(f) \right )}{\sqrt{ S_{xx}(f)S_{yy}(f) }} \\ LC_{xy}(f) &=& \frac{\mathrm{Im} \left (C_{xy}(f) \right )}{\sqrt{ 1 - \left [ \mathrm{Re}\left ( C_{xy}(f) \right ) \right ]^{2} }} = \frac{\mathrm{Im} \left (S_{xy}(f) \right )}{\sqrt{ S_{xx}(f)S_{yy}(f) - \left [ \mathrm{Re}\left ( S_{xy}(f) \right ) \right ]^{2} }} \\ \end{eqnarray*}$

where $\mathrm{Im()}$ and $\mathrm{Re()}$ describe the imaginary and real parts of a complex number.

To calculate coherence values in Brainstorm, select the process.

Process options

Process Connectivity > Coherence NxN:
Time window: Segment of the signal used for the connectivity analysis. Select All file.
Remove evoked response: If checked, removes the average of all the files in input (the "trials") from each file, before computing the connectivity measure. Meaningful only in the context of ERP/ERF analyses. Uncheck it.
Process options: Different measures. Select: Magnitude squared coherence.
Window length: Duration in seconds for the spectrum estimation. Set to: 1s.
Overlap: Percentage of overlap between consecutive windows. Set to: 50%.
Highest frequency: After the computation, removes all the frequencies above this threshold, mostly for visualization purposes. It should be <= Fs/2. Set to: 60 Hz.
Output options: Select: Save individual results.

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:

Display as graph: Plot the connectivity graph at a given frequency point.
Display as image: Plot the connectivity matrix at a given frequency point.
Power spectrum: Plot coherence as a function of frequency for all the possible node pairs.

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.

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.

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 $$x(t)$$ (using a linear autoregressive model) is improved by adding signal $$y(t)$$ (also using a linear autoregressive model). If this is true, signal has a Granger causal effect on the first signal. In other words, independent information of the past of signal improves the prediction of signal $$x(t)$$ obtained with the past of signal 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 $$x(t)$$ and $$y(t)$$ , as well as the addition of a multiple of to .
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 $$w$$ is causing both $$x$$ and $$y$$ , but with a smaller delay for $$y$$ than for $$x$$ , then the GC measure between $$x$$ and $$y$$ would show a non-zero GC for $$y$$ --> $$x$$ , even though $$y$$ is not truly causing $$x$$ (Bressler and Seth, 2011).

Process options

Process Connectivity > Bivariate Granger causality NxN
Time window: Segment of the signal used for the connectivity analysis. Select All file.
Remove evoked response: If checked, removes the average of all the files in input (the "trials") from each file, before computing the connectivity measure. Meaningful only in the context of ERP/ERF analyses. It is recommended by some authors as it meets the zero-mean stationarity requirement (improves stationarity of the system). However, the problem with this approach is that it does not account for trial-to-trial variability. For a discussion see (Wang et al., 2008). Uncheck it.
Maximum Granger model order: The most common criteria used to define the order of the model are the Akaike’s information criterion, the Bayesian-Schwartz’s criterion, and the Hannan-Quinn criterion. Too low orders may lack the necessary details, while too big orders tend to create spurious values of connectivity. While our simulated signals were created with a model of 4, here we used as model order of 6 for a decent connectivity result.
Output options: Select Save individual results.

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.

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

Process: Connectivity > Bivariate Granger causality NxN.

With respect to GC, spectral GC presents two extra parameters:
Maximum frequency resolution: Width of frequency bins in PSD estimation. Set to 1 Hz.
Highest frequency: After the computation, removes all the frequencies above this threshold, mostly for visualization purposes. It should be <= Fs/2. Set to: 60 Hz.

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.

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, $\tilde{x}(t)$ , is a complex signal uniquely associated to a real signal, $$x(t)$$ , that has been useful in signal processing due to its characteristics, more specifically, its module $a_{\tilde{x}}(t)$ , and phase $\phi_{\tilde{x}}(t)$ , correspond to the instantaneous amplitude (or envelope) and instantaneous phase of the associated real signal . The real part of $\tilde{x}(t)$ is its associated real signal $$x(t)$$ , and the imaginary part is the Hilbert transform of the same real signal .

$\begin{eqnarray*} \tilde{x}(t)= x(t) + j\mathcal{H}\left\{ x(t) \right\} = a_{\tilde{x}}(t)e^{j\phi_{\tilde{x}}(t)} \\ \end{eqnarray*}$

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

Process: Connectivity > Envelope Correlation N×N [2020]
Time window: Segment of the signal used for the connectivity analysis. Select: All file.
Remove evoked response: If checked, removes the average of all the files in input (the "trials") from each file, before computing the connectivity measure. Meaningful only in the context of ERP/ERF analyses. Uncheck it.
Time-frequency transformation method: Either Hilbert transform or Morlet wavelets. Each of this methods requires additional parameters that are found in an external panel that opens by clicking on Edit. See the time-frequency tutorial. In this example, the sampling frequency of the signals being 120Hz, make sure you remove the "gamma2" band.
Signal splitting: This process has the capability of splitting the input data into several blocks for performing time-frequency transformation, and then merging them to build a single file. This feature helps to save a huge amount of memory and, at the same time, avoids breaking a long-time recording to short-time signals, which makes inconsistency in dynamic network representation of spontaneous data. The maximum number of blocks which can be specified is 20. Set to 1.
Connectivity measure: This is the connectivity metric that will be used with the envelopes. Select Envelope correlation (orthogonalized).
Time resolution: If Dynamic is selected, connectivity is saved for each window of analysis. If Static is selected, the connectivity results from all the windows are averaged. Select Dynamic.
Estimation window length: Duration in milliseconds to compute the connectivity measure. Set to 5000 ms
Sliding window overlap: Percentage of overlap between consecutive windows to compute the connectivity measure. Set to 50%.
Use the parallel processing toolbox: Enables the use of the parallel processing toolbox in Matlab to accelerate the computational procedure. Uncheck it.
Output configuration: Generally, the above calculation results in a 4-D matrix, where dimensions represent channels (1st and 2nd dimensions), time points (3rd dimension), and frequency (4th dimension). In the case that we analyze event-related data, we have also several files (trials). However, due to the poor signal-to-noise ratio of a single trial, an individual realization of connectivity matrices for each of them is not in our interests. Consequently, we need to average connectivity matrices among all trials of a specific event. The second option of this part performs this averaging.

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

First 5-s window

Last 5-s window

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 $\tilde{x}(t)$ and $\tilde{y}(t)$ , obtained from the TF transformation using the Hilbert transform:

$\begin{eqnarray*} \mathrm{PLV} = \left | E\left [ e^{j\Delta \phi (t)} \right ] \right | \\ \end{eqnarray*}$

with:

$\begin{eqnarray*} \Delta \phi (t) = \phi_{\tilde{x}}(t) - \phi_{\tilde{y}}(t) = arg\left ( \frac{\tilde{x}(t)\tilde{y}^{*}(t)}{\left | \tilde{x}(t) \right |\left | \tilde{y}(t) \right |} \right ) \\ \end{eqnarray*}$

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

Process: Connectivity > Phase locking value NxN
Time window: Segment of the signal used for the connectivity analysis. Select All file.
Frequency bands: Used for the time-frequency transformation with the Hilbert transform method. See the time-frequency tutorial. In this example, the sampling frequency of the signals being 120Hz, make sure you remove the "gamma2" band.
Method: PLV, ciPLV or wPLI (see references).
Keep time information: Computes PLV across trials, thus the result is a PLV time series for each frequency band.
Measure: How the resulting vector is reported. Select Magnitude.
Output options: Select Save individual results.

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.

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 $\phi_{\tilde{x}}(t)$ causes the signal $\phi_{\tilde{y}}(t)$ , the mutual information, between $\phi_{\tilde{y}}(t)$ and the past of $\phi_{\tilde{x}}(t)$ i.e. $\phi_{\tilde{x}}(t')$ is larger than the mutual information of $\phi_{\tilde{y}}(t)$ , the past of $\phi_{\tilde{y}}(t)$ i.e. $\phi_{\tilde{y}}(t')$ and $\phi_{\tilde{x}}(t')$ . This relationship can be seen on the Venn diagram below, where $\mathrm{I()}$ and $\mathrm{H()}$ , indicate mutual information and the individual entropies respectively. Lastly, PTE cannot be negative, and its magnitude does not have a meaningful upper bound.

Process options

Process: Connectivity > Phase Transfer Entropy NxN
Time window: Segment of the signal used for the connectivity analysis. Select All file.
Frequency bands: Used for the time-frequency transformation with the Hilbert transform method. See the time-frequency tutorial. In this example, the sampling frequency of the signals being 120Hz, make sure you remove the "gamma2" band.
Return normalized phase transfer entropy: Divides the each directional PTE value between two signals by the sum of both directional PTE values for those two signals.

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.

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.

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.

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

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.

TF: [Nr x Ntime x Nfreq] matrix containing all the connectivity values. The connectivity matrix R computed in bst_connectivity.m between two sets of signals A and B, is optimized and saved in the TF variable. The size of the R matrix is [Na x Nb x Ntime x Nfreq]. The relation between Na x Nb and Nr depends on the type of optimization in the file. See details below.
RefRowNames: Cell-array of strings {Na x 1}. Labels of the rows of the connectivity matrix R (Y axis in the image display, with the label From). In the case of a AxB process, these labels are the names of the signals extracted from the first list of files (FilesA).
RowNames: Cell-array of strings {Nb x 1}. Labels of the columns of the connectivity matrix R (X axis in the image display, with the label To). In the case of a AxB process, these labels are the names of the signals extracted from the second list of files (FilesB).
Freqs: Double [1 x Nfreq] for a simple list of frequencies; or cell-array {Nfreq x 3} if using frequency bands, where each line represents a band {'band_name', 'frequency definition', 'function'}
Time: [1,Ntime] for time-resolved files; [1,2] for files with no data dimension (e.g. correlation). When no time is available, this variable describes the time segment from which the connectivity measure was computed.
Options.isSymmetric: Boolean indicating whether the matrix R saved in the TF field is symmetric (e.g. NxN correlation or coherence) or not (e.g. Granger causality). If the R matrix is symmetric, it is saved in a compressed format, with only the values from lower triangular matrix.

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

Nolte G, Bai O, Wheaton L, Mari Z, Vorbach S, Hallett M.
Identifying true brain interaction from EEG data using the imaginary part of coherency.
Clinical Neurophysiology. 2004 Oct;115(10):2292–307.
Pascual-Marqui RD.
Coherence and phase synchronization: generalization to pairs of multivariate time series, and removal of zero-lag contributions.
arXiv preprint arXiv:0706.1776. 2007 Jun 12.
Bressler SL, Seth AK.
Wiener–Granger causality: a well established methodology.
Neuroimage. 2011 Sep 15;58(2):323-9.
Dhamala M, Rangarajan G, Ding M.
Estimating Granger causality from Fourier and wavelet transforms of time series data.
Physical review letters. 2008 Jan 10;100(1):018701.
Hipp JF, Hawellek DJ, Corbetta M, Siegel M, Engel AK.
Large-scale cortical correlation structure of spontaneous oscillatory activity.
Nature neuroscience. 2012 Jun;15(6):884-90.
Tass P, Rosenblum MG, Weule J, Kurths J, Pikovsky A, Volkmann J, et al.
Detection of n : m Phase Locking from Noisy Data: Application to Magnetoencephalography.
Phys Rev Lett. 1998 Oct 12;81(15):3291–4.
Bruña R, Maestú F, Pereda E
Phase locking value revisited: teaching new tricks to an old dog
Journal of Neural Engineering, Jun 2018
Vinck M, Oostenveld R, van Wingerden M, Battaglia F, Pennartz CM
An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias
Neuroimage, Apr 2011
Lobier M, Siebenhühner F, Palva S, Palva JM.
Phase transfer entropy: a novel phase-based measure for directed connectivity in networks coupled by oscillatory interactions.
Neuroimage. 2014 Jan 15;85:853-72.
Barzegaran E, Knyazeva MG.
Functional connectivity analysis in EEG source space: The choice of method
Ward LM, editor. PLOS ONE. 2017 Jul 20;12(7):e0181105.
Lai M, Demuru M, Hillebrand A, Fraschini M.
A comparison between scalp- and source-reconstructed EEG networks.
Sci Rep. 2018 Dec;8(1):12269.
Schoffelen J-M, Gross J.
Source connectivity analysis with MEG and EEG.
Hum Brain Mapp. 2009 Jun;30(6):1857–65.
Kriegeskorte N, Simmons WK, Bellgowan PSF, Baker CI.
Circular analysis in systems neuroscience: the dangers of double dipping.
Nat Neurosci. 2009 May;12(5):535–40.
Cohen MX.
Analyzing Neural Time Series Data: Theory and Practice. MIT Press; 2014.
Sadaghiani S, Brookes MJ, Baillet S.
Connectomics of human electrophysiology.
NeuroImage. 2022 Feb;247:118788.

Forum discussions

Granger causality: https://neuroimage.usc.edu/forums/t/12506
Coherence on single trials: https://neuroimage.usc.edu/forums/t/22726
ciPLV: https://neuroimage.usc.edu/forums/t/9751
wPLI: https://neuroimage.usc.edu/forums/t/9434
Coherence and PLV: https://neuroimage.usc.edu/forums/t/33379
Removing ERP response for PLV: https://neuroimage.usc.edu/forums/t/32665
Connectivity of subcortical ROIs: https://neuroimage.usc.edu/forums/t/33479
Non-zero diagonals with unconstrained sources: https://neuroimage.usc.edu/forums/t/33700

Scripting

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

function tutorial_connectivity(reports_dir)
% TUTORIAL_CONNECTIVITY: Script that runs the Brainstorm connectivity tutorial.
%
% INPUTS: 
%    - 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: Raymundo Cassani, 2021-2022
%         Francois Tadel, 2022


%% ===== PARAMETERS =====
% Output folder for reports
if (nargin < 1) || isempty(reports_dir) || ~isfolder(reports_dir)
    reports_dir = [];
end


%% ===== CREATE PROTOCOL  =====
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Create Protocol
ProtocolName = 'TutorialConnectivity';
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
% Start a new report
bst_report('Start');


%% ===== SIMULATE DATA (MVAR MODEL) =====
% Seed for random number generator
rng(111); 
% Process: Simulate AR signals
sFileSim = bst_process('CallProcess', 'process_simulate_ar_spectra', [], [], ...
    'subjectname',  'Subject01', ...
    'condition',    'Simulation', ...
    'samples',      12000, ...
    'srate',        120, ...
    'interactions', ['1, 1 / 10, 25 / 0.3, 0.5' 10 ...
                     '2, 2 / 10, 25 / 0.7, 0.3' 10 ... 
                     '3, 3 / 10, 25 / 0.2, 0.2' 10 ... 
                     '1, 3 / 25 / 0.1']);

% Make sure the display mode is "columns"
bst_set('TSDisplayMode', 'column');
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sFileSim, [], ...
    'type',    'data', ...  % Recordings time series
    'Comment', 'Simulated signals');

                 

%% ===== CORRELATION =====
% Process: Correlation NxN
sFiles = bst_process('CallProcess', 'process_corr1n', sFileSim, [], ...
    'timewindow', [], ...
    'scalarprod', 0, ...
    'outputmode', 1);  % Save individual results (one file per input file)

% Process: Snapshot: Connectivity matrix
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'connectimage', ...  % Connectivity matrix
    'Comment', 'Correlation NxN');


%% ===== COHERENCE =====
% Process: Magnitude-squared coherence:  |C|^2 = |Cxy|^2/(Cxx*Cyy)
sFiles = bst_process('CallProcess', 'process_cohere1n', sFileSim, [], ...
    'timewindow',    [], ...
    'removeevoked',  0, ...
    'cohmeasure',    'mscohere', ...  % Magnitude-squared coherence:  |C|^2 = |Cxy|^2/(Cxx*Cyy)
    'tfmeasure',     'stft', ...  % Fourier transform
    'tfedit',        struct(...
         'Comment',         'Complex', ...
         'TimeBands',       [], ...
         'Freqs',           [], ...
         'StftWinLen',      1, ...
         'StftWinOvr',      50, ...
         'StftFrqMax',      60, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'none', ...
         'Output',          'all', ...
         'SaveKernel',      0), ...
    'timeres',       'none', ...  % None
    'avgwinlength',  1, ...
    'avgwinoverlap', 50, ...
    'outputmode',    'input');  % separately for each file

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'spectrum', ...  % Frequency spectrum
    'Comment', 'MSC NxN');

% Process: Snapshot: Connectivity graph
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'connectgraph', ...  % Connectivity graph
    'Comment', 'MSC NxN');

% Process: Snapshot: Connectivity matrix
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'connectimage', ...  % Connectivity matrix
    'Comment', 'MSC NxN');

% Process: Imaginary coherence:  IC = |imag(C)|
sFiles = bst_process('CallProcess', 'process_cohere1n', sFileSim, [], ...
    'timewindow',    [], ...
    'removeevoked',  0, ...
    'cohmeasure',    'icohere2019', ...  % Imaginary coherence:  IC = |imag(C)|
    'tfmeasure',     'stft', ...  % Fourier transform
    'tfedit',        struct(...
         'Comment',         'Complex', ...
         'TimeBands',       [], ...
         'Freqs',           [], ...
         'StftWinLen',      1, ...
         'StftWinOvr',      50, ...
         'StftFrqMax',      60, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'none', ...
         'Output',          'all', ...
         'SaveKernel',      0), ...
    'timeres',       'none', ...  % None
    'avgwinlength',  1, ...
    'avgwinoverlap', 50, ...
    'outputmode',    'input');  % separately for each file

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'spectrum', ...  % Frequency spectrum
    'Comment', 'Imaginary coherence NxN');

% Process: Lagged coherence / Corrected imaginary coherence:  LC = |imag(C)|/sqrt(1-real(C)^2)
sFiles = bst_process('CallProcess', 'process_cohere1n', sFileSim, [], ...
    'timewindow',    [], ...
    'removeevoked',  0, ...
    'cohmeasure',    'lcohere2019', ...  % Lagged coherence / Corrected imaginary coherence:  LC = |imag(C)|/sqrt(1-real(C)^2)
    'tfmeasure',     'stft', ...  % Fourier transform
    'tfedit',        struct(...
         'Comment',         'Complex', ...
         'TimeBands',       [], ...
         'Freqs',           [], ...
         'StftWinLen',      1, ...
         'StftWinOvr',      50, ...
         'StftFrqMax',      60, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'none', ...
         'Output',          'all', ...
         'SaveKernel',      0), ...
    'timeres',       'none', ...  % None
    'avgwinlength',  1, ...
    'avgwinoverlap', 50, ...
    'outputmode',    'input');  % separately for each file

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'spectrum', ...  % Frequency spectrum
    'Comment', 'Lagged coherence NxN');


%% ===== GRANGER CAUSALITY =====
% Process: Bivariate Granger causality NxN
sFiles = bst_process('CallProcess', 'process_granger1n', sFileSim, [], ...
    'timewindow',   [], ...
    'removeevoked', 0, ...
    'grangerorder', 6, ...
    'outputmode',   1);  % Save individual results (one file per input file)

% Process: Snapshot: Connectivity graph
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'connectgraph', ...  % Connectivity graph
    'Comment', 'Granger causality NxN');

% Process: Snapshot: Connectivity matrix
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'connectimage', ...  % Connectivity matrix
    'Comment', 'Granger causality NxN');


%% ===== SPECTRAL GRANGER CAUSALITY =====
% Process: Bivariate Granger causality (spectral) NxN
sFiles = bst_process('CallProcess', 'process_spgranger1n', sFileSim, [], ...
    'timewindow',   [], ...
    'removeevoked', 0, ...
    'grangerorder', 6, ...
    'maxfreqres',   1, ...
    'maxfreq',      60, ...
    'outputmode',   1);  % Save individual results (one file per input file)

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'spectrum', ...  % Frequency spectrum
    'Comment', 'Spectral Granger causality NxN');


%% ===== ENVELOPE CORRELATION =====
% Process: Envelope Correlation NxN [2023]
sFiles = bst_process('CallProcess', 'process_henv1n', sFileSim, [], ...
    'timewindow',    [], ...
    'removeevoked',  0, ...
    'cohmeasure',    'oenv', ...  % Envelope correlation (orthogonalized)
    'tfmeasure',     'hilbert', ...  % Hilbert transform
    'tfedit',        struct(...
         'Comment',         'Complex', ...
         'TimeBands',       [], ...
         'Freqs',           {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'}}, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'none', ...
         'Output',          'all', ...
         'SaveKernel',      0), ...
    'timeres',       'windowed', ...  % Windowed
    'avgwinlength',  5, ...
    'avgwinoverlap', 50, ...
    'parallel',      0, ...
    'outputmode',    'input');  % separately for each file

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'spectrum', ...  % Frequency spectrum
    'Comment', 'Envelope correlation (Hilbert transform) NxN');

% Process: Envelope Correlation NxN [2023]
sFiles = bst_process('CallProcess', 'process_henv1n', sFileSim, [], ...
    'timewindow',    [], ...
    'removeevoked',  0, ...
    'cohmeasure',    'oenv', ...  % Envelope correlation (orthogonalized)
    'tfmeasure',     'morlet', ...  % Morlet wavelets
    'tfedit',        struct(...
         'Comment',         'Complex,1-60Hz', ...
         'TimeBands',       [], ...
         '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], ...
         'MorletFc',        1, ...
         'MorletFwhmTc',    3, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'none', ...
         'Output',          'all', ...
         'SaveKernel',      0), ...
    'timeres',       'windowed', ...  % Windowed
    'avgwinlength',  5, ...
    'avgwinoverlap', 50, ...
    'parallel',      0, ...
    'outputmode',    'input');  % separately for each file

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'spectrum', ...  % Frequency spectrum
    'Comment', 'Envelope correlation (Morlet wavelet) NxN');


%% ===== PHASE LOCKING VALUE =====
plv_variants = {'plv',   ...  % Phase locking value
                'ciplv', ...  % Lagged phase synchronization / Corrected imaginary PLV
                'wpli'};      % Weighted phase lag index
for ix = 1 : length(plv_variants)
    % Process: Phase locking value
    sFiles = bst_process('CallProcess', 'process_plv1n', sFileSim, [], ...
        'timewindow',    [], ...
        'plvmethod',     plv_variants{ix}, ...
        'plvmeasure',    2, ...  % Magnitude
        'tfmeasure',     'hilbert', ...  % Hilbert transform
        'tfedit',        struct(...
             'Comment',         'Complex', ...
             'TimeBands',       [], ...
             'Freqs',           {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'}}, ...
             'ClusterFuncTime', 'none', ...
             'Measure',         'none', ...
             'Output',          'all', ...
             'SaveKernel',      0), ...
        'timeres',       'none', ...  % None
        'avgwinlength',  1, ...
        'avgwinoverlap', 50, ...
        'outputmode',    'input');  % separately for each file

    % Process: Snapshot: Frequency spectrum
    bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
        'type',    'spectrum', ...  % Frequency spectrum
        'Comment',  ['Phase locking value (' plv_variants{ix} ')  NxN']);
end

%% ===== PHASE TRANSFER ENTROPY =====
% Process: Phase Transfer Entropy NxN
sFiles = bst_process('CallProcess', 'process_pte1n', sFileSim, [], ...
    'timewindow', [], ...
    'freqbands',  {'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'}, ...
    'normalized', 0, ...
    'outputmode', 1);  % Save individual results (one file per input file)

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
    'type',    'spectrum', ...  % Frequency spectrum
    'Comment', 'Phase transfer entropy NxN');


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

TODO

Hossein, Richard

Rewrite bst_henv.m to work with two inputs: https://github.com/brainstorm-tools/brainstorm3/pull/523
Check new process Envelope correlation vs old processes Amplitude envelope correlation
Validate that we can deprecate safely the old processes

Raymundo, Sylvain

Read data: Can we find something more meaningful than what we compute in the example in the dataset? A full connectome in a task vs. baseline maybe?

Francois

Threshold connectivity matrix: show only the top N% measures
Add 1xN and NxN submenus (Sylvain's request)

Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):

Connectivity

Introduction

Definitions

Sensor-level

Source-level

Requirements

Simulated data

Correlation

Process options

Result visualization

Coherence

Process options

Result visualization

Granger causality

Process options

Result visualization

Spectral Granger causality

Process options

Result visualization

Envelope correlation

Process options

Result visualization

Phase locking value

Process options

Result visualization

Phase transfer entropy

Process options

Result visualization

Method selection and comparison

Scout-level connectivity

Matrix thresholding

On the hard drive

File name

File structure

Connectivity matrix encoding

Additional documentation

Related tutorials

Articles

Forum discussions

Scripting

TODO