Bispectral pairwise interacting source analysis for identifying systems of cross-frequency interacting brain sources from electroencephalographic or magnetoencephalographic signals

Brain cognitive functions arise through the coordinated activity of several brain regions, which actually form complex dynamical systems operating at multiple frequencies. These systems often consist of interacting subsystems, whose characterization is of importance for a complete understanding of the brain interaction processes. To address this issue, we present a technique, namely the bispectral pairwise interacting source analysis (biPISA), for analyzing systems of cross-frequency interacting brain sources when multichannel electroencephalographic (EEG) or magnetoencephalographic (MEG) data are available. Specifically, the biPISA makes it possible to identify one or many subsystems of cross-frequency interacting sources by decomposing the antisymmetric components of the cross-bispectra between EEG or MEG signals, based on the assumption that interactions are pairwise. Thanks to the properties of the antisymmetric components of the cross-bispectra, biPISA is also robust to spurious interactions arising from mixing artifacts, i.e., volume conduction or field spread, which always affect EEG or MEG functional connectivity estimates. This method is an extension of the pairwise interacting source analysis (PISA), which was originally introduced for investigating interactions at the same frequency, to the study of cross-frequency interactions. The effectiveness of this approach is demonstrated in simulations for up to three interacting source pairs and for real MEG recordings of spontaneous brain activity. Simulations show that the performances of biPISA in estimating the phase difference between the interacting sources are affected by the increasing level of noise rather than by the number of the interacting subsystems. The analysis of real MEG data reveals an interaction between two pairs of sources of central mu and beta rhythms, localizing in the proximity of the left and right central sulci.

Figure 1

The two-step computation of the antisymmetric bispectral tensor between channel signals, ${\mathbf{B}}_{\left[\mathbf{i}\right|\mathbf{j}\left|\mathbf{k}\right]}$, from the antisymmetric bispectral tensor between source signals, ${\mathcal{B}}_{\left[\mathit{m}\right|\mathit{n}\left|\mathit{p}\right]}$.

Figure 2

Simulated data generation. (a) Two independent oscillators at 6 and 10 Hz (top signals) were generated by band-pass filtering white Gaussian noise around the respective frequencies, with 1 Hz bandwidth. A 16-Hz oscillator (bottom signal) was generated through a multiplicative interaction, i.e., a time point by time point multiplication, between the oscillators at 6 and 10 Hz, followed by band-pass filtering around 16 Hz with 1 Hz bandwidth. Short segments of the oscillator time courses and the respective power spectral densities (PSDs) are shown in panels (b) and (c). The oscillators generated in this way resulted in quadratic phase coupling, and they were used for the construction of pairwise interacting source time courses according to two different scenarios of interaction. In the first scenario [panel (d)], the time course of source 1 was obtained by summing the time courses of all the three oscillators, and the time course of source 2 was set to a time delayed copy of the time course of source 1, with a time delay $\tau$. Up to three interacting source pairs were generated, with $\tau$ being 5 ms for the first pair, 10 ms for the second pair, and 15 ms for the third pair. As a result, both sources contained the three components at frequencies 6, 10, and 16 Hz. In the second scenario of interaction [panel (e)], the time course of source 1 was set to the time course of the oscillator at 6 Hz, whereas the time course of source 2 was set to the sum of the time courses of the 10- and 16-Hz oscillators. In order to account for a time delay between the two sources, in this scenario of interaction the 16-Hz oscillator was generated by multiplying the oscillator at 10 Hz by a time delayed copy of the oscillator at 6 Hz, with the time delay $\tau$ being as in the previous case.

Figure 3

The plots show, for all combinations of pair numbers, i.e., $Q$, and SNRs, the Frobenius norms of the 30 normalized matrices (ordered in abscissa), i.e., whose entries have been normalized as in Eq. (14), returned by the SVD-based factorization described in the paper. Data are from one representative simulation repetition, corresponding to scenario I of interaction between sources.

Figure 4

Scenario I of interaction. Histograms of $1/(1-r)$, with $r$ being the smallest of the $M$ canonical correlation coefficients between the true and estimated $M$-dimensional subspaces spanned by interacting source topographies, for all combinations of pair numbers, i.e., $Q$, and SNRs. $m$ denotes the mean value. Data from 1000 simulation repetitions.

Figure 5

Scenario II of interaction. Histograms of $1/(1-r)$, with $r$ being the smallest of the $M$ canonical correlation coefficients between the true and estimated $M$-dimensional subspaces spanned by interacting source topographies, for all combinations of pair numbers, i.e., $Q$, and SNRs. $m$ denotes the mean value. Data from 1000 simulation repetitions.

Figure 6

Contrast $C_q$ between the magnitudes of the coefficients ${\alpha }_q$ and ${\beta }_q$ obtained in the two simulated scenarios of interaction. The black dots denote the mean value; the error bars denote the standard deviations.

Figure 7

Scenario I of interaction. Histograms of the estimated phase differences between interacting sources, for all combinations of pair numbers, i.e., $Q$, and SNRs. The vertical dashed line denotes the true values for phase differences, i.e., $\mathrm{\Delta }{\varphi }_1=0.314,\mathrm{\Delta }{\varphi }_2=0.628$, and $\mathrm{\Delta }{\varphi }_3=0.942$. $m$ and $s$ denote the mean value and the standard deviation, respectively. Data from 1000 simulation repetitions.

Figure 8

(Top) The magnitude of the antisymmetric component of bicoherence, $\left|b_{\left[i\right|j\left|k\right]}(f_1,f_2)\right|$, is shown as function of frequencies $f_1$ and $f_2$. (Bottom) A detailed view of $\left|b_{\left[i\right|j\left|k\right]}(f_1,f_2)\right|$ for the diagonal axis $f_1=f_2$.

Figure 9

Frobenius norm of the normalized matrices (ordered in abscissa) returned by the SVD-based factorization of tensor slices for the antisymmetric component of the cross-bispectrum at $(f_1,f_2)=(11\text{Hz},11\text{Hz})$.

Figure 10

Topographies (first and third lines) and the corresponding reconstruction of cortical activity (second and fourth lines) for pairwise interacting sources of brain mu (11 Hz) and beta (22 Hz) rhythms. The maps have been scaled between $-1$ and 1 for topographies and between 0 and 1 for cortical activity. The maps use arbitrary units.