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 measures are designed to probe how brain regions (or nodes) interact as a network. A distinction is made between structural (fiber pathways), functional (non-directed statistical associations) and effective (causal interactions, or "directed functional connectivity") connectivity between regions. Here we explain how to compute various connectivity metrics for non-directed and directed functional connectivity analyses with Brainstorm, both with simulated (ground-truth) and empirical data.

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

Contents

Introduction
Pre-requisites
Simulated data
Correlation
Coherence
Granger causality
Spectral Granger causality
Envelope correlation
Phase-Locking Value
Phase-Transfer Entropy
Method selection and comparison
Scout-level connectivity
Thresholding of connectivity estimates
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 a connectome) can be represented as a connectivity graph (left panel below), with each brain region as a node (x, y, z,...), and the connectivity measures shown above each edge of the graph. The connectome can also be represented by a connectivity array, a.k.a. an adjacency matrix (right panel below).

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

Directed vs. non-directed: The direction of interactions between signals (as a measure of statistical causation) is estimated with directed metrics. Non-directed metrics produce symmetrical connectivity graphs/matrices: the connectivity measure "from Signal $$x$$ to Signal $$y$$ " is identical to the connectivity measure "from Signal $$y$$ to Signal $$x$$ ".

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

Full (NxN) vs. seeded (1xN) connectivity: In a full connectivity analysis, the connectivity metric is computed for all the possible pairs of nodes between N time series (noted N×N here). Alternatively, seeded connectivity (noted 1×N) is measured 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, or Morlet wavelets.

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

Sensor-level analyses

Their anatomical interpretation is limited and ambiguous.
Sensor data is affected by field spread and volume conduction of brain activity across the scalp/sensor array. Hence, activity from one single brain area is picked up at multiple, possibly distant sensor locations. Connectivity measures may show strong interactions between sensor time series, which may be wrongly interpreted as inter-regional brain connections.

Source-level analyses

Connectivity measures between source time series are neuroanatomically interpretable and can be derived across participants.

We recommended that the outcomes of sensor and source connectivity analyses are mutually compatible; see 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.6 GB. 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) may be required for practical reasons, depending on computing resources.

Regions of interest

Grouping individual brain sources into ROIs, defined on the cortical surface or in the brain volume helps contain the computational resources required for connectivity analyses (see above). ROIs can be defined based on neuroscience priors and the purpose of the study (Schhoffen and Gross, 2009). You may consider:

Study priors from the literature in your field, specific working hypotheses in terms of neuro-anatomical brain structures expected to be involved.
Association with a non-neuronal signal (e.g., cortico-muscular coherence).
Signal strength, by defining ROIs from brain regions with strongest activity in current data.
Whole-brain cortical parcellation to reduce N.

The tutorial Corticomuscular coherence explains the computation of connectivity measures between one sensor (EMG) and 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 below.

Studies have reported similarities and differences between electrophysiological connectomes (from MEG/EEG) and fMRI connectomes. See how the electrophysiological connectomes provide unique insights on how functional communication is implemented in the brain in Sadaghiani et al., 2022).

On whole-brain connectivity analyses and the issue of circular analysis, please see (Kriegeskorte et al., 2009).

Pre-requisites

Here we focus on the specifics on connectivity analyses with Brainstorm. Please make sure you are familiar with Brainstorm and go through all introduction tutorials first, before following the present tutorial.

We first use simulated, synthetic 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 in the 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 synthetic time series simulated with known ground truth interactions. We will use a multivariate autoregressive (MVAR) model that consists of the following three signals:

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

The 25-Hz component of Signal 3 will be driven in part by the 25-Hz component of Signal 1 (denoted Signal 1>>Signal 3).

Let's now generate the 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 where to store the simulated data in your database: Set to Simulation.
Number of time samples for each simulated time series: Set to 12,000.
Sampling frequency (Fs): 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 generate the time series with the MVAR model.
In the next sections, we will compute different connectivity metrics from these synthetic time series. Drag and drop the simulated data in the Process1 tab, click on [Run] ( ) to open the Pipeline editor, and select the connectivity metric.
This tutorial illustrates the computation of full connectivity graphs (Process 1: NxN). It is also possible to obtain connectivity measures between one of the time series of a given data file and the other time series in the same file (use Process1: 1xN), or between time series from two different data files (use Process2: AxB), as shown further below.

References used:

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 relatively simple, non-directed connectivity metric of association between two time series. It is relatively limited without further preprocessing of the input time series: correlation is sensitive to volume conduction and is not frequency specific.

Process options

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

Result visualization

The results are stored as a N×N connectivity file, with a new icon . Right-click for display options:
Display as graph: Plots the connectivity graph using a chord diagram where the color of the edges shows the value of the connectivity metric. See the connectivity graph tutorial for a detailed explanation of the options of this visualization.
Display as image: Plots the connectivity measures in an adjacency matrix.
Display fibers: Displays source connectivity edges as virtual white fiber tracts, as shown here.

Display options:

Show labels: Click on the figure to display the names of the times series and the connectivity values in the figure legend. To display the labels corresponding to each column and row: right-click on the figure > Figure > Show labels. If the time series names are too long, use the option Use short labels.
Hide self-connectivity: In a square NxN connectivity matrix, the diagonal values correpond to the values of the connectivity metric between a time series and itself. These self-connectivity values are typically much higher than the other values of the matrix. When checked, these values are hidden so that to display the off-diagonal values more clearly.
Colormap: By default, the NxN colormap is configured to display the absolute values of the connectivity measures. Correlation values may be positive or negative. To display the 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. Magnitude square coherence (MSC), $|C_{xy}(f)|^2$ , or coherence, measures the covariance of two signals in the frequency domain. For a pair of time series $$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 address the problem of volume conduction, are imaginary coherence (Nolte et al., 2004), $IC_{xy}(f)$ , and 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()}$ are the imaginary and real parts of a complex number, respectively.

Process options

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

Result visualization

Coherence is frequency specific: a connectivity graph and a connectivity matrix are produced for each frequency bin of the spectrum. Right-click on the coherence result file to show display options:

Display as graph: Plots a connectivity graph for each frequency bin.
Display as image: Plot a connectivity matrix for each frequency bin.
Power spectrum: Plot coherence as a function of frequency, between each time series.

Open the 3 types of visualization: they are linked by clicking on the spectral representation of the coherence, which will change the current frequency bin and update the connectivity measures displayed in the graph and matrix. The value of the current frequency bin can also be changed in the Time panel.

Other variants of coherence can be derived following the above procedure. The figures below display the coherence spectra obtained from the imaginary coherence (left) and the lagged coherence (right) measures.

Granger causality

Granger causality (GC) is a measure of directed functional connectivity based on the Wiener-Granger causality framework. GC measure linear dependencies between time series, and tests whether the prediction of the future of signal $$x(t)$$ (approximated by a linear autoregressive model) is improved by considering signal $$y(t)$$ (also approximated by a linear autoregressive model). If there is such improvement, one concludes that signal has a Granger causal effect on the other signal. In sum, some independent information from the past of signal improves the prediction of the future of signal $$x(t)$$ , with respect to only observing its past . GC takes nonnegative values; its is zero when no Granger causality can be attributed. The term independent means that GC is invariant with the scale of the input signals and $$y(t)$$ .
See Granger causality - mathematical background for a more complete background.

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: Time segment of the time series used for the connectivity analysis. Select All file.
Remove evoked response: See same option for Coherence above. This option is recommended by some authors as it satisfies the zero-mean stationarity of the GC model, but does not account for trial-to-trial variability; see (Wang et al., 2008). Leave it Unchecked.
Maximum Granger model order: The most common criteria used to define the order of the autoregressive models used by GC to approximate the input time series are the Akaike’s information criterion, the Bayesian-Schwartz’s criterion, and the Hannan-Quinn criterion. Low model orders will yield coarse approximations of the input time series. High model orders may yield spurious GC connectivity estimates. The simulated time series were created with an AR model of order 4; please use a maximum Granger model order of 6.
Output options: Select Save individual results.

Results visualization

In the connectivity graph (left panel below) the directionality is shown with an arrowhead at the center for the arc connecting two nodes. The connectivity matrix (right panel below) is not symmetric because GC values are direction specific. The upper right elements of the matrix indicate there is a GC influence from signal 1 to signal 3.

Spectral Granger causality

The Spectral Granger causality is a measure of directed functional connectivity that was developed to indicate frequency specific influences between time series (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: computes GC values under the specificied frequency. Should be set <= Fs/2. Use: 60 Hz.

Result visualization

As with coherence, spectral GC can be plotted as a function of frequency. The display below shows a peak around 25 Hz that corresponds to the expected causal interaction of Signal 1 on Signal 3.

Envelope correlation

In the time-frequency tutorial, we have introduced Morlet wavelets and the Hilbert transform as methods to decompose time series in the time-frequency (TF) domain.

The outcome of the Hilbert transform is an analytic signal, $\tilde{x}(t)$ , which is a complex time series uniquely associated to the original data time series, $$x(t)$$ , which module $a_{\tilde{x}}(t)$ , and phase $\phi_{\tilde{x}}(t)$ , correspond to the instantaneous amplitude (or envelope) and instantaneous phase of the original time series , respectively. The real part of $\tilde{x}(t)$ is the original time series , and the imaginary part is the Hilbert transform of that same time series $$x(t)$$ .

$\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 estimations of instantaneous amplitude and phase. While it is technically feasible to derive these elements from broadband time series, they would not be interpretable (see Cohen, 2014).

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

An optional step consists in orthogonalizing the pairs of envelope time series used to compute the amplitude envelope correlation, by first removing the real part of their coherence to reduce volume conduction (at the sensor level) or cross-talk effects (at the source level) (Hipp et al., 2012).

Process options

Process: Connectivity > Envelope Correlation N×N [2020]
Time window: Time segment of the signal used for the connectivity analysis. Select: All file.
Remove evoked response: See Coherence above. Leave it unchecked.
Time-frequency transformation method: Hilbert transform or Morlet wavelets. Each method requires additional parameters that will be displayed in a separate panel after clicking on Edit. See the time-frequency tutorial. In the present example, the sampling frequency of the signals being 120Hz, make sure you remove the "gamma2" band from the targetted frequency bands.
Signal splitting decomposes the input data into time blocks before performing time-frequency transformations, and merging the latter back into a single output file. This feature saves on memory, especially when processing long recordings or trials. The maximum number of blocks is 20. Set to 1.
Connectivity measure is the connectivity metric applied to the resulting signal envelopes. Select Envelope correlation (orthogonalized).
Time resolution: If Dynamic, connectivity is saved for each time window. If Static, the connectivity measures are averaged across all windows. Select Dynamic.
Estimation window length: Duration in milliseconds over which the connectivity measure is computed. Set to 5000 ms
Sliding window overlap: Percentage of overlap between consecutive time windows to compute the connectivity measure. Set to 50%.
Use the parallel processing toolbox: If available, will call Matlab's parallel processing to accelerate computations. Leave it Unchecked.
Output configuration: The process outcome is a 4-D array, with channels (1st and 2nd dimensions), time points (3rd dimension), and frequency (4th dimension). With event-related data, the outcomes from each trial would be stored in separate connectivity files. Here we will average these connectivity outcomes across all trials.

Results visualization

As with coherence and spectral Granger causality, amplitude envelope correlations (AEC) can be displayed as a function of frequency, and as a function of time if the Time resolution option was set to Dynamic. Below are the AEC results obtained with the Hilbert transform (left panels) and with Morlet wavelets (right panels) for the first 5-s time window (top panels) and the 5-s last window (bottom panels).

First 5-s window

Last 5-s window

Phase-Locking Value

An alternative class of connectivity metrics considers the relative instantaneous phase between two time series as a marker of connectivity. Phase-locking or phase synchronization are two such measures (Tass et al., 1998). In principle, phase-based measures are robust against fluctuations in signal amplitude, although low signal-to-noise conditions remain challenging for these measures (Lachaux et al., 1999; Mormann et al., 2000).

The phase-locking value (PLV) is a popular metric defined as the length of the average vector of many unit vectors whose phase angle corresponds to the phase difference between two time series (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 0. If the phases of the two signals are strongly coupled, the length of the average vector will be close to 1. For event-related studies, we would expect phase differences across trials to be uniformly distributed, unless the phase of the brain signals is locked to the event onset, which would then indicate a form of phase locking between the two time series. Considering a pair of narrow-band analytic signals $\tilde{x}(t)$ and $\tilde{y}(t)$ , obtained from the Hilbert transform (see above):

$\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*}$

Several variants of PLV have been contritributed. Brainstorm features
ciPLV (Bruña 2018) and wPLI (Vinck 2011) variants contributed to Brainstorm by Daniele Marinazzo.

PLV values tend to be overestimated when derived from few (<50) time samples. As a consequence, when comparing PLV values between conditions, please make sure that they are derived from about the same number of samples in each condition.

Process options

Process: Connectivity > Phase locking value NxN
Time window: Time 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, please remove the "gamma2" band.
Method: PLV, ciPLV or wPLI (see references).
Keep time information: Computes PLV across trials instead of across time: the outcome is a PLV time series for each frequency band.
Measure specifies how the resulting phase-locking vector is reported. Select Magnitude.
Output options: Select Save individual results.

Results visualization

PLV is frequency resolved, and here 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 is indeed our observation, with a peak at 22 Hz (center of beta band) in the PLV spectrum.

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 timeseries influences the ability to predict the future of the target timeseries. 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 represented on the Venn diagram below, where $\mathrm{I()}$ and $\mathrm{H()}$ , indicate mutual information and the individual entropies, respectively. PTE values are positive and unbounded.

Process options

Process: Connectivity > Phase Transfer Entropy NxN
Time window: Time 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 120 Hz, make sure the "gamma2" band is removed from the option.
Return normalized phase transfer entropy: Divides each directional PTE value between two signals by the sum of both directional PTE values for those two signals.

Results visualization

Here we computed PTE for the delta, theta, alpha, beta and gamma bands. From the synthetic data used here, we expect to observe a higher PTE value in the beta band, From Signal 1 To Signal 3. This is indeed shown with a peak at 22 Hz (center of beta band) in the PTE spectrum.

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 pairs of time series. When computing whole-brain connectomes (i.e., NxN connectivity matrices between all brain sources), additional technical considerations need to be taken into account: the number of time series is typically too large to be computationally tractable on most workstations (see introduction above); if using unconstrained source maps, with 3 time series at each brain location, this requires an extra step to derive connectivity estimates.

The tutorial Corticomuscular coherence explains the computation of connectivity measures between one sensor and brain source maps, with sensor x sources and sensor x scouts scenarios, both with constrained and unconstrained source orientation models.

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

Scout function: In this configuration, one connectivity result (i.e., coherence spectrum) is computed for each pair of scouts, for each elementary source orientation. It is necessary to specify two parameters to indicate how the data is aggregated in each scout: The scout function (for instance, the mean), and the within-scout aggregation procedure (before or after the computation of the connectivity measure).

Unconstrained maps / maximum: The graphs below show two scouts (Scout-1 and Scout-2), with 3 orientations each (x,y,z). Coherence is computed between each pair of elementary sources in all three orientations (1x-2x, 1x-2y, 1x-2z, 1y-2x, ..., 1z-2z), yielding 9 coherence spectra. From these 9 values at each frequency bin, only the maximum value is kept to summarize the connectivity estimate between the two scout regions" The connectvity outcome measure is therefore one coherence spectrum per pair of scouts. The choice of the maximum statistic is empirical: it is not rotation invariant (if the specifications of the NAS/LPA/RPA fiducials change, the x,y,z axes also change, which may influence conenctivity measures quantitatively). Nevertheless, our empirical tests showed this option produced smooth and robust (to noise) connectivity estimates.

Before: The scout function (e.g., mean) is applied on each direction on the elementary source time series in the scout, which results in one time series per elementary direction, per scout. These scouts time series are then used to compute coherence measures, and the resulting coherence spectra are finally aggregated across the x, y and z dimensions, yielding one single coherence spectrum for each pair of scouts.

After: With this option, coherence is first computed between each pair of elementary sources (number of sources in each scout x 3 orientations) for two scouts, yielding 9 coherence spectra. Then, these coherence spectra are aggregated across dimensions to obtain one single coherence spectrum for each pair of elementary sources within the scouts. Finally, the scout function (e.g., mean) is applied to the coherence spectra for each elementary source within the scouts to obtain one single coherence spectrum for each pair of scouts. This is a considerably more greedy option that the "before" option above, as it computes coherence measures between all 45000x45000 elementary source signals. This requires more computing resources; see the introduction for an example.

The result is a connectivity file () that contains Scouts x Scouts coherence spectra.

The connectivity graph or the adjacency matrix are displayed for each frequency bin of the coherence spectra. For more details, see the connectivity graph tutorial.

Thresholding of connectivity estimates

Interactive tools

Brainstorm provides tools to apply empirical thresholds to display the resulting connectivity estimates. Interactions can be filtered interactively based on their magnitudes, the distance between the nodes of the graph or their category. See the tutorial: Connectivity graphs.

Threshold by percentile

The process Test > Threshold by percentile allows the selection of the top n% connectivity estimates in a single input file.

All connectivity estimates (left panel below). The top-5% connectivity matrix (right panel below) only keeps the 5 percentile of all the connectivity estimates.

Statistical significance

Inference based threshold can be derived to assess the statistical significance of the connectivity estimates, for instance using non-parametric paired permutation tests.

Connectivity outcomes can be computed for different experimental conditions, or between an active state and a baseline. Non-parametric paired permutation tests are available in Brainstorm to compare between two groups of repeated measures. More information in the Statistics tutorial.

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 the first connectivity file computed > File > View file contents.

The data structure is the same as for Brainstorm time-frequency files. Only the fields specific to connectivity analyses are documented below:

TF: [Nr x Ntime x Nfreq] array containing all 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 (i.e., not compressed): for each time and frequency, the list of values from the first dimension of the TF variable is 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
Coherence/Lagged coherence vs. Envelope correlation/Lagged coherence: https://neuroimage.usc.edu/forums/t/lagged-coherence/39440

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

Granger Causality implementation: https://github.com/brainstorm-tools/brainstorm3/pull/433

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?

Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):

Connectivity

Introduction

Definitions

Sensor-level analyses

Source-level analyses

Pre-requisites

Simulated data

Correlation

Process options

Result visualization

Coherence

Process options

Result visualization

Granger causality

Process options

Results visualization

Spectral Granger causality

Process options

Result visualization

Envelope correlation

Process options

Results visualization

Phase-Locking Value

Process options

Results visualization

Phase-Transfer Entropy

Process options

Results visualization

Method selection and comparison

Scout-level connectivity

Thresholding of connectivity estimates

Interactive tools

Threshold by percentile

Statistical significance

On the hard drive

File name

File structure

Connectivity matrix encoding

Additional documentation

Related tutorials

Articles

Forum discussions

Scripting

TODO