Detection of Cross-Frequency Coupling Between Brain Areas: An Extension of Phase Linearity Measurement

The current paper proposes a method to estimate phase to phase cross-frequency coupling between brain areas, applied to broadband signals, without any a priori hypothesis about the frequency of the synchronized components. N:m synchronization is the only form of cross-frequency synchronization that allows the exchange of information at the time resolution of the faster signal, hence likely to play a fundamental role in large-scale coordination of brain activity. The proposed method, named cross-frequency phase linearity measurement (CF-PLM), builds and expands upon the phase linearity measurement, an iso-frequency connectivity metrics previously published by our group. The main idea lies in using the shape of the interferometric spectrum of the two analyzed signals in order to estimate the strength of cross-frequency coupling. We first provide a theoretical explanation of the metrics. Then, we test the proposed metric on simulated data from coupled oscillators synchronized in iso- and cross-frequency (using both Rössler and Kuramoto oscillator models), and subsequently apply it on real data from brain activity. Results show that the method is useful to estimate n:m synchronization, based solely on the phase of the signals (independently of the amplitude), and no a-priori hypothesis is available about the expected frequencies.

Introduction

Brain areas need to constantly transfer information among themselves to put in place complex behavioral responses to the environment (Bressler, 1995). Functional connectivity is defined as the presence of statistical dependencies between the time-series representing the activity of brain regions (Friston, 1994; Buzsáki and Draguhn, 2004). A variety of mechanisms through which this communication occurs are summarized in Jensen and Colgin (2007), involving only the phase (Tass et al., 1998) or also amplitude (Canolty et al., 2006). Each of these phenomena would underlie a specific neuro-physiological mechanism (for a review, see Engel et al., 2013). In the literature, a wide number of metrics have been proposed to detect each of these kinds of communication (Le Van Quyen and Bragin, 2007; Tort et al., 2010). Furthermore, communication between brain areas can occur either in iso-frequency or in cross-frequency. Cross-frequency coupling (CFC) is the interaction occurring between neuronal populations operating at different frequencies. It has been postulated that this form of synchronization could represent a suitable option to allow large-scale synchronizations across distant areas in the brain (Varela et al., 2001; Canolty and Knight, 2010), yielding the integration of distributed information (Jirsa and Müller, 2013). Moreover, definite (both frequency and spatial) patterns of CFC have been shown to be the neuro-physiological substrate underlying the recruitment of areas needed for the execution of tasks such as specific kinds of learning (Schack and Weiss, 2005; Tort et al., 2009; Kendrick et al., 2011), segregation of interfering inputs (Colgin et al., 2009), perception (Doesburg et al., 2009; Siebenhühner et al., 2016), encoding of reward (Cohen et al., 2009) or sensory processing (Seymour et al., 2017). In human brain activity, two main forms of cross-frequency coupling have been described so far. Firstly, the phase of slow oscillations can modulate the amplitude of faster activity (Vanhatalo et al., 2004; Palva et al., 2005). Furthermore, phase-phase synchronization has also been described, whereby the phases of “n” cycles of a signal are locked to “m” phase cycles of another signal (Tass et al., 1998). This kind of cross-frequency communication, classically defined as n:m synchronization, has been observed previously in human brain data (Nikulin and Brismar, 2006) and is the only mechanism capable of supporting CFC at high temporal resolution (Fell and Axmacher, 2011).

Several metrics have been developed to capture the presence of cross-frequency communication. For instance, phase-amplitude coupling (Tort et al., 2010) can successfully detect the presence of nested-synchronization, while metrics such as bicoherence (Sigl and Chamoun, 1994) can detect cross-frequency, phase-phase coupling. However, bicoherence is not a pure phase-based metrics as its value depends also on the amplitude, preventing an unambiguous interpretation of the involved neuronal mechanisms (Palva et al., 2005). The biphase-locking value, while purely based on the phase, also provides an estimate of the phase-amplitude coupling (Darvas et al., 2009). Metrics such as the phase-locking factor (Palva et al., 2005) detect pure phase to phase locking, but require an accurate a priori hypothesis about the frequencies involved in cross-frequency synchronization. The procedure proposed by Cohen (2008), on the other hand, while not requiring any a priori hypothesis, focuses on phase-amplitude coupling.

Each approach has its own advantages and drawbacks and, when one is dealing with specific task-related data, given that a specific a priori hypothesis is available about the frequencies across which synchronization might be occurring, the application of these metrics is effective (Jirsa and Müller, 2013). However, when dealing with resting-state data, the situation changes because the signals contain several frequency bands interacting with each other (possibly with more than one of the mentioned mechanisms) at once.

Restricting the analysis to phase-to-phase coupling, we have to consider that the bandwidth of the involved signals is so broad and complex to potentially allow the simultaneous occurrence of iso- and cross- frequency synchronizations at once (Varela et al., 2001; Canolty et al., 2006; Jirsa and Müller, 2013; Cabral et al., 2014). Disentangling these complex signals has proven to be elusive when one does not know a priori if, when and where cross-frequency is occurring within the brain. The aim of this article is to develop a reliable estimation of phase-phase cross-frequency communication between the broadband signals of two brain regions, without a priori hypothesis on the frequencies at which such a synchronization might occur. To do this, we build and expand upon the phase linearity measurement (PLM), an iso-frequency phase-based connectivity metrics recently developed by our group (Baselice et al., 2019).

One issue is related to the amount of potential combinations of frequencies and areas that one should test in order to look for CFC throughout the brain and throughout the frequency spectrum. Indeed, an attempt to identify the frequency at which cross-frequency synchronization is present from the data by selecting a number of combinations of possible frequencies has been done (Sauseng et al., 2009), using the level of synchronization across trials in order to statistically estimate where cross-frequency synchronization was present.

To this regard, a new method has been recently proposed, that does not require any a priori hypothesis and can estimate cross-frequency synchronization (Volk et al., 2018). Such an approach estimates from the data the “candidate frequencies” where the CFC might be occurring. However, when performing this procedure, a maximization of the correlation between the signals is performed, hence reintroducing a form of dependency from the amplitude. The issue of the communication between different frequencies has also been addressed using a multiplex network approach (Yu et al., 2017). The idea is that each layer of the multiplex network represents, at a specific frequency, the iso-frequency correlations between brain areas. However, the cross-layers links are not estimated from the data. With the methodology proposed in this article, we aim at providing a data-informed estimate of which brain areas and frequencies are involved in cross-frequency phase-to-phase coupling. The novelty of this work lies in the fact that no a priori information is required about the frequencies and the areas involved in the CFC. On the contrary, our technique allows to start from wide signal spectrum and to detect if cross-frequency is occurring and, if so, to identify which frequency components are involved per each signal. Firstly, we provide a theoretical description of the metric. Secondly, the metric is tested in a number of synthetic analytical models. We first used Rössler oscillators, which capture the non-linearities of the brain. Secondly, in order to simulate the simultaneous presence of iso and cross frequency synchronization, we implemented several Kuramoto oscillators, and introduced a lag between the generated signals. This procedure is known to produce the appearance of synchronization at a lower frequency bandwidth as compared to the original signals (Niebur et al., 1991a). Hence, we tested the ability of the newly proposed methodology, namely cross-frequency phase linearity measurement (CF-PLM), to detect and disentangle both kinds of synchronism. Furthermore, we mixed the previously produced signals linearly, in order to obtain a case in which both iso-frequency and cross-frequency coupling are simultaneously present, and we tested if the newly proposed approach can disentangle such a situation. Finally, we tested the metrics on source-reconstructed MEG data (acquired by the MEG laboratory in Naples), and identified brain areas where cross-frequency is present that are spatially consistent across the tested subjects.

Methods

Definition of the Interferometric Signal

Let us define x(t) and y(t) as the time series related to two brain areas. By applying the Hilbert transform, their analytical expression is obtained:

$\begin{array}{l}{x}_{an}(t)=A_x(t)e^{i{\varphi }_x(t)}=x_R(t)+i{x}_I(t),\\ y_{an}(t)=A_y(t)e^{i{\varphi }_y(t)}=y_R(t)+i{y}_I(t).& (1)\end{array}$

where variables A and ϕ represent the amplitude and the phase, respectively. According to this mathematical description, signals generated by brain areas can be modeled as complex phasors with time-varying amplitude and phase.

According to Baselice et al. (2019), their phase-to-phase connectivity can be measured via a three steps procedure. Firstly, the normalized interferometric component of the two signals z(t) is computed:

$\begin{array}{ll}z(t)=\frac{x_{an}(t)y_{an}^{\prime }(t)}{\left|x_{an}(t)\right|\left|y_{an}(t)\right|}=e^{i\Delta \varphi (t)},& (2)\end{array}$

where the symbol ′ indicates the complex conjugate. Note that the complex interferometric function z(t) has an amplitude equal to 1 [thus it is independent of the amplitudes of the signals A_x(t) and A_y(t)], and a phase term Δϕ(t) = ϕ_x(t) − ϕ_y(t) ∈ [−π, π[, which is the time-varying phase difference between the phases of x_an(t) and y_an(t).

It is known that, due to the Hilbert transform, the analytical expression and therefore the phase of a signal can be computed only in case of a narrow frequency band. Although the considered time series x(t) and y(t) are band limited, their frequency range cannot be considered narrow, and thus an effective estimation of the genuine phase is not guaranteed by the Hilbert transform (Rosenblum et al., 2021). That said, in the following it will only be required the knowledge of the phase difference term Δϕ(t), instead of the genuine phases ϕ_x(t) and ϕ_y(t). It has been shown that, for some tasks such as the quantification of the frequency locking, the Hilbert-based estimation of the phase, i.e., the so-called protophase, could suit well (Kralemann et al., 2008). Therefore, according to Baselice et al. (2019), we will assume the phase term of z(t) as effective for our purpose, i.e., the coupling measurement.

In particular, the behavior of the term Δϕ(t) can be exploited in order to measure phase connectivity between ϕ_x(t) and ϕ_y(t) signals. Hence, the frequency analysis of the function z(t) is carried out. Three different conditions could occur, as reported in Figure 1. In case of no synchrony between the sources, the interferometric phase values appear to be irregularly spread in the [−π, π[range (blue line in Figure 1). In case of phase coupling, the term Δϕ(t) will be characterized by a linear trend. That is, if the two sources have a similar oscillation frequency, the phase of the interferometric signal will be constant or slowly varying in time (red line in Figure 1) while in case of two sources oscillating at different frequencies, a slope will appear (yellow line in Figure 1).

FIGURE 1

Figure 1. Interferometric phase signals in three different conditions: independent sources (blue line), coupled iso-frequency sources (red line) and coupled sources at different frequencies (yellow line).

Once the complex signal z(t) has been computed, the second step for measuring the coupling consists in computing its power spectrum by means of the Fourier transform:

$\begin{array}{ll}\begin{array}{c}{S}_Z(f)=|{\displaystyle {\int }_0^{T}z(t)e^{-i2\pi ft}dt}{|}^2,\end{array}& (3)\end{array}$

where [0, T] is the observation period. In order to have a more reliable evaluation of the PSD function, we implemented the periodogram estimator with a rectangular window and confidence interval of 0.95 for the computation of S_Z(f) (Auger and Flandrin, 1995). The shape of the power spectrum is strongly influenced by the strength of the coupling occurring between the two sources and by their central frequency, and hence it can be exploited to estimate them (Baselice et al., 2019).

Phase Linearity Measurement

In Figure 2, the power spectra occurring in the different scenarios are represented. The blue line does not show any peak, in accordance with the absence of coupling between the sources. This means that the power spectrum of z(t) is almost flat (blu line of Figure 2) if its phase term Δϕ(t) irregularly spreads in the [−π, π] range (blue line of Figure 1).

FIGURE 2

Figure 2. Power spectral densities of the interferometric components (i.e., the power spectrum of the phases of the interferometric signal) in three different conditions: independent sources (blue line), coupled iso-frequency sources (red line), coupled sources at different frequencies (yellow line). The presence of a power peak denotes the coupling between sources, while its position indicates the difference in their resonant frequencies.

The red line shows an evident power peak around f = 0, which is due to a linear behavior of the interferometric phase Δϕ(t), i.e.:

$\begin{array}{ll}\Delta \varphi (t)={\phi }_{xy}+\Delta ft,& (4)\end{array}$

where the term Δf is related to the different central frequencies of the two sources. In the case of iso-frequency coupling (IFC), such a term is relatively small, resulting in a peak centered around f = 0. In this case, the last step for measuring connectivity strength consists in computing the percentage of power within a narrow band [−B, B] around f = 0:

$\begin{array}{ll}\begin{array}{c}PLM=\frac{{\displaystyle {\int }_{-B}^B{S}_Z(f)df}}{{\displaystyle {\int }_{-\infty }^{\infty }{S}_Z(f)df}}.\end{array}& (5)\end{array}$

In Baselice et al. (2019) it has been shown that a B value of 1 Hz is a well balanced trade-off between the discrimination capability and the estimation noise of the algorithm. The PLM approach has shown a good performance in measuring the iso-frequency coupling, i.e., in distinguishing between the case of the blue line and the red line in Figure 2. Nevertheless, it has to be modified in order to make it effective in analyzing the last case, the cross-frequency coupling.

Cross-Frequency PLM

In the CFC condition, a non-minimal frequency difference Δf occurs between the coupled components of the sources, and such difference produces a shift in the interferometric spectrum, as shown by the yellow curve of Figure 2. In this case, the coupling is evident due to the presence of the peak, which is now centered at f = Δf (6 Hz in the reported case) instead of f = 0. This difference makes the PLM (Equation 5) unable to capture the coupling, as the power is no longer concentrated in the [−B, B] band. One should notice that the knowledge of the frequency difference Δf would solve the problem, as the integration could be shifted accordingly into the [Δf −B,Δf +B] frequency range. However, this kind of a priori knowledge is not available at all times, as it is often the case in resting-state, as well as in many task-related settings. This situation can be handled by looking for maxima (i.e., power peaks) in the interferometric power spectrum S_Z(f). Once a local maximum is identified (besides the one centered in 0), its power and position can be easily measured. This is what the proposed methodology implements. In other words, once the PSD function of Equation (3) is computed, the global maximum is identified. By retrieving its position, the difference between the two sources central frequencies Δf is identified. Subsequently, the coupling strength is measured by adapting the upper integral of Equation (5), i.e.:

$\begin{array}{ll}\begin{array}{c}\text{CF-PLM}=\frac{{\displaystyle {\int }_{\Delta f-B}^{\Delta f+B}{S}_Z(f)df}}{{\displaystyle {\int }_{-\infty }^{\infty }{S}_Z(f)df}}.\end{array}& (6)\end{array}$

Volume conduction or field spread is a common source of error in case of coupling metrics for EEG or MEG data. Due to its definition, the CF-PLM is insensitive to this effect. The reason is that the Volume Conduction effect produce a zero-lag copy of the considered signal, with the same frequency. Therefore its contribution to the interferometric PSD will be at f = 0, i.e., out of the [Δf −B, Δf +B] integration band. However, the provided information concerns the frequency difference between the two sources, the central frequencies of the coupled components still have to be determined. Hence, the last step is to identify the oscillation frequencies of the two components involved in the CFC. To this aim, a band-stop Gaussian-shaped frequency filter has been adopted. The stop band is centered at f_H and is 2B large. The central frequency f_H is moved in order to scan the whole frequency range of the acquired signals, e.g., the [0.5, 48] Hz range, as reported in Figure 3 (top). Let us focus on the signal of the first source [i.e., x(t) of Equation (1)]. Once the filter is overlapped to frequency components involved in the coupling and removes them, the peak of the interferometric PSD S_Z(f) disappears, as shown in Figure 3 (center). The filter position f_H will reveal the frequency f_x of x(t) involved in the coupling. The same process is repeated for the second source y(t) for the identification of f_y, according to Figure 3 (bottom). After the frequency scans, the two central frequencies of the components involved in the CFC f_x and f_y are identified, while the amount of coupling is related to the peak energy and is measured via Equation (6).

FIGURE 3

Figure 3. Scheme of the procedure for the identification of frequencies involved in the coupling. When the frequency stop filters are not overlapped to the frequencies involved in the coupling, the peak in the interferometric PSD is present (top). When one filter overlaps with the frequency of the first (center), or the second source (bottom), the peak disappears and there is a reduction in the power.

Results

The proposed methodology has been tested on both synthetic and real datasets. In case of simulated data, two approaches have been adopted for generating the cross-frequency coupled signals, exploiting Rössler attractors and Kuramoto oscillators, respectively. In more detail, the sensibility of the CF-PLM metrics to coupling strength has been analyzed by means of Rössler attractors signals to which a frequency shift has been applied. Furthermore, a modified version of the Kuramoto oscillators implementing signals with different central frequencies has been considered, in order to test the ability of the CF-PLM to identify the two frequencies involved in the coupling. As a third analysis, real data acquisitions have been considered for the final validation of the approach.

Rössler Attractors

Two time series have been generated according to Rosenblum et al. (2001):

$\begin{array}{lll}{ẋ}_{1,2}& =& -2\pi f_{1,2}{y}_{1,2}-z_{1,2}+{\xi }_{1,2}+c(x_{2,1}-x_{1,2})\\ {ẏ}_{1,2}& =& 2\pi f_{1,2}{x}_{1,2}+0.15y_{1,2}\\ {ż}_{1,2}& =& 0.2+z_{1,2}(x_{1,2}-10)& (7)\end{array}$

with coupling strength c varying between 0 and 0.04, which can be considered high coupling according to Rosenblum et al. (2001). A frequency f₁ equal to 10 Hz has been chosen, while the duration and the sampling interval have been set equal to 420 s and 625 Hz, respectively. The two coupled time series x(t) and y(t) have been generated with a central frequency f₁. Subsequently, the cross-frequency has been simulated by applying a frequency shift to the second attractor y(t) via the modulation property of the Fourier transform:

$\begin{array}{ll}{y}_{\Delta f}(t)=y(t)e^{-i2\pi \Delta ft}.& (8)\end{array}$

For this analysis, we considered Δf = 7 Hz. The CF-PLM has been computed between x(t) and y_Δf(t). Several analyses have been conducted aiming at evaluating the sensitivity of the proposed metrics with respect to the coupling strength of the attractors, the Signal to Noise Ratio (SNR), i.e., the ratio between signal and noise powers, and the frequency shift. In Figure 4, the values measured by the CF-PLM as a function of attractor's coupling strength are reported in case of different SNR levels (using white, additive noise). In more detail, a Monte Carlo simulation with 50 iterations has been implemented and the mean values are reported. As expected, the CF-PLM value increases as a function of the coupling strength, for each of the considered noise levels. Moreover, we tested the CF-PLM in case of several frequency shifts, obtaining the curves with the same behavior of those reported in Figure 4. A further analysis has been conducted. In particular, the CF-PLM values have been computed by varying the number of available samples of the two Rössler attractors, again in case of different coupling strength. Results are reported in Figure 5.

FIGURE 4

Figure 4. Mean values of CF-PLM measured in case of two Rössler oscillators varying their coupling strength from 0 (no coupling) to 0.04 (high coupling). Results are reported in case of different SNR values between 0 and 30 dB.

FIGURE 5

Figure 5. CF-PLM values measured in case of two Rössler oscillators varying the number of samples.

Kuramoto Oscillators

Three mutually coupled Kuramoto oscillators, namely s₁(t), s₂(t), and s₃(t), have been generated according to the following model (Cabral et al., 2011):

$\begin{array}{ll}\frac{d{\theta }_n}{dt}=2\pi f_n+k{\displaystyle \sum _{p=1}^3}Csin({\theta }_p(t-\tau )-{\theta }_n(t))& (9)\end{array}$

with n = 1, 2, 3, τ = 0.6 s and C = 1. The central frequencies have been set equal to 10, 10, and 17 Hz, respectively. The first oscillator (s₁(t), f = 10 Hz) has been compared to the second one (s₂(t), f = 10 Hz), to the third one (s₃(t), f = 17 Hz) and to the sum of the last two (s₂(t) + s₃(t), f = 10 Hz and 17 Hz). The three PSDs of the interferometric signals S_Z(f) are reported in Figure 6. The PSD peak in case of the iso-frequency synchronization of two 10 Hz oscillators is clearly visible (the power peak centered at 0 Hz of Figure 6A), as well as in case the cross-frequency synchronization occurring between s₁(t) and s₃(t) (the power peak centered around -7 Hz of Figure 6B). Importantly, the case of multiple components simultaneously synchronized in iso and cross-frequency is correctly handled, with the two power peaks positioned at 0 and -7 Hz in Figure 6C. It has to be underlined that the amplitudes of the considered oscillations plays an important role. In other works, the two peaks shown in Figure 6C have the same height because the oscillators amplitude is similar. In case of different intensities, the power represented by the peaks will be different, and in some circumstances a masking effect could occur.

FIGURE 6

Figure 6. Power spectra of the interferometric signals in case of coupled Kuramoto oscillators: iso-frequency (A), cross-frequency (B), iso and cross-frequency (C).

Since the Kuramoto oscillators are coupled with a time delay between them, the frequency shift depends not only on the natural frequencies of each oscillator but also on the coupling strength between them (Niebur et al., 1991b; Choi et al., 2000). In Figure 7, we show that the position of the interferometric peak is shifted as a function of the coupling strength, thus validating the existence of a cross-frequency interaction between the oscillators. In other words, this shows that the presence of synchrony, at the frequency that is predicted theoretically, is captured by the metric (as opposed to merely capturing n:m phase relationships, whose frequencies are not expected to be dependent from the coupling strength).

FIGURE 7

Figure 7. Position of the interferometric peak while varying the global coupling strength of two Kuramoto oscillators.

Let us now analyze how the frequencies involved in the connectivity process are identified. According to the processing scheme previously described and reported in Figure 3, two stop-band filters are implemented in the frequency domain. The peak power reduction is computed when moving the central frequencies of these filters within the [0, 20] Hz range. Results are reported in Figure 8 for all the considered cases. When the filter removes from the first signal the frequencies involved in the coupling, a power reduction is measured in the PSD peak. As a consequence, a horizontal dark line will appear in the images of Figure 8. Analogously, a vertical line will appear when the corresponding frequency of the second source is removed. The result is a cross-shaped image, with the center identifying the two frequencies involved in the coupling. By looking at Figure 8, it is evident that the maximum power reduction appears at (10, 10) Hz in the case of s₁(t), s₂(t) coupling, at (10, 17) Hz in the case of s₁(t), s₃(t) coupling, while in the simultaneous iso and cross-frequency coupling two couples are correctly identified at (10, 10) and (10, 17) Hz, respectively. All these results are in accordance with what we expected, as the procedure correctly estimates both the connectivity strength and the oscillator frequencies involved in the coupling from the interferometric spectrum.

FIGURE 8

Figure 8. Results of the analysis for the identification of frequencies involved in the coupling in case of different Kuramoto oscillators: iso-frequency (A), cross-frequency (B), and simultaneous iso- and cross- frequency (C). The center of each cross identifies the frequencies of the two oscillators.

In order to have a benchmark, the dual-frequency coherence (DFC) (Mellors et al., 1998), which is a normalized version of the second order bispectrum (Thomson, 1982), and the n:m Phase Synchronization Index (PSI) (Tass et al., 1998; Wacker and Witte, 2011) have been implemented. Given the two acquired signals x(t) and y(t) and their Fourier transform X(ω) and Y(ω), the DFC can be defined as Shahbazi Avarvand et al. (2018):

$\begin{array}{ll}DFC({\omega }_1,{\omega }_2)=\frac{|<X({\omega }_1)Y^{\prime }({\omega }_2)>|}{<|X({\omega }_1){|}^2{>}^{1/2}<|Y({\omega }_2){|}^2{>}^{1/2}},& (10)\end{array}$

in which the symbol < . > stands for the average over trials operation. In our case, we assumed the number of trials to be equal to 1. The n:m PSI is defined as:

$\begin{array}{c}PSI(n,m)=\frac{1}{N}|{\displaystyle \sum _{k=1}^{K}exp[i(n{\varphi }_x^k-m{\varphi }_y^k)]|,}\end{array}$

in which the time series have been divided in K segments. We computed the DFC for ω₁ and ω₂ in the [1, 20] Hz range, and the n:m PSI for n and m values between 1 and 20, obtaining the result reported in Figure 9. Compared to Figure 8B, it is evident that DFC metrics is less effective in determining the frequencies involved in the coupling as two maxima are present at a distance of about 1 Hz, while the n:m PSI is not effective at all.

FIGURE 9

Figure 9. Dual-Frequency Coherence (left) and n:m Phase Synchronization Index (right) results in case of CFC Kuramoto oscillators.

Real Data

Acquisition and Preprocessing

The acquisitions used for the analysis are from healthy subjects acquired at the MEG facility in Naples¹. The detailed procedure used for the processing of the acquisitions has been described in Sorrentino et al. (2018). In brief, subjects were seated in a 163-magnetometers MEG system. The brain activity was recorded twice for 3.5 min, with a small break to minimize the chances of drowsiness. After the anti-aliasing filter, the data were sampled at 1024 Hz, and filtered between 0.5 and 48 Hz with a 4th order Butterworth IIR band-pass filter (Oostenveld et al., 2011). During the acquisitions, the electrocardiogram (ECG) and the electrooculogram (EOG) were also recorded (Gross et al., 2013). Principal component analysis (PCA) was used to reduce the environmental noise (Sadasivan and Narayana, 1996; Oostenveld et al., 2011). Subsequently, noisy channels were removed manually through visual inspection by trained experts. For each subject, supervised independent component analysis (ICA) (Barbati et al., 2004) was performed to eliminate the ECG and, if present, the EOG components from the MEG signals. MEG data were then co-registered to the native MRI of the subjects. We used the volume conduction model proposed by Nolte (2003) and the linearly constrained minimum variance (LCMV) beamformer (Van Veen et al., 1997) to reconstruct the time-series related to the centroids of 90 regions-of-interest (ROIs), derived from the automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002; Gong et al., 2009; Hillebrand et al., 2016). For each source, we projected the time series along the dipole direction explaining most variance by means of singular value decomposition (SVD), obtaining a scalar value per each source.

Connectivity Measurement

The power spectra of the interferometric signal for each couple among the 90 sources have been computed. In the following, we selected two couples of regions, i.e., the couple with the highest CFC peak (among all couples of regions), and a region with an average CFC peak intensity. To show that a high CFC peak is unlikely to appear by chance, we have validated the analysis against surrogates (Lancaster et al., 2018). Hence, the intensity of the highest CFC peak, derived from 10000 random surrogates, obtained by shuffling the phase of the signal in the frequency domain, have been computed. In Figure 10 (top), we show that, in the case of the high observed cross-frequency peak, the peak intensity was above the 99th percentile of the surrogates. On the contrary, in case of the average CFC peak [Figure 10 (bottom)], its intensity was around the 50th percentile, as expected. To further check the validity of our analysis in real subjects, we went on to estimate the anatomical consistency of the CFC peaks per each link across subjects. To do so, we proceeded as follows: for each of the 2 subjects, for each source pair the intensity of the strongest CFC peak has been measured (|Δf| ≥ 2 Hz). Subsequently, we binarized the CFC peak matrix according to a threshold. To avoid a dependency of the result from the choice of the threshold, we have repeated the analysis across a range of thresholds, i.e., from the 1st to the 99th percentile. Then, we have summed across the binary matrices obtaining one matrix per subject which does not depend on any specific threshold. Finally, the subject specific matrices were summed, showing that the presence of a CFC peak is topographically consistent across the two individuals. In the first row of Figure 11, the source couples that are beyond a high number of thresholds, i.e., high percentile values, are reported in progressively intense red. The most evident points of this connectivity matrix are those related to the strongest CFC peaks. We also projected such information on a 3D brain template, as shown in Figure 12. In the second row of Figure 11, the delta-frequencies corresponding to the significant edges are reported.

FIGURE 10

Figure 10. Empirical CDF of the peaks distributions obtained by shuffling the signal phases in the frequency domain (blue lines). Values have been normalized with respect to the highest CFC peak present in the data before shuffling (red line). The analysis refers to two cases: the highest CFC peak found among all couples of regions (top), where the peak intensity is above the 99% of the distribution, and a couple of regions with an average CFC peak intensity (bottom), where the peak is in the middle of the distribution.

FIGURE 11

Figure 11. First row: peak intensity, second row: peak frequency from CF-PLM matrixof the 2 subjects obtained by filtering the interferometric peaks intensity. Rows and columns are the considered sources (between 1 and 96 of the AAL atlas). The threshold values have been selected according to the percentiles of the distribution, spanning between 1 and 100. The most intense red points characterizes the CF-PLM peaks with higher intensity.

FIGURE 12

Figure 12. For one subject, the brain shows the intensity of the cross-frequency coupling in each region. That is, based on the CFC matrix, the weighted-degree has been computed for each region. The higher the value, the more a given region is connected in cross-frequency to other metrics.

Now we are going to focus on the couple of regions with the highest CFC peak. In Figure 13A, the PSD of the interferometric signal related to the right inferior parietal lobule and the orbital part of the right superior frontal gyrus are reported. Results are related to one epoch about 150 seconds long of a single subject. A power peak positioned at 9.5 Hz is clearly visible, which shows the presence of cross-frequency coupling. In the following step the sources frequencies involved in the coupling are determined. As explained earlier, the procedure consists in filtering the two signals and measuring the power reduction of the frequency peak. In Figure 13B, it is evident that the highest reduction is found in case of f₁ ≈ 11 Hz and of f₂ ≈ 1.5 Hz. The result obtained by the DFC, which is reported in Figure 13C, is not effective in identifying the involved frequencies, as the global maximum is hardly distinguishable. This result could be explained as the DFC has not been designed for resting state coupling analysis.

FIGURE 13

Figure 13. Analysis of signals from areas 57 (right inferior parietal lobule) and 42 (orbital part of the right superior frontal gyrus). (A) PSD of the interferometric signal. The peak related to the cross-frequency coupling, located at +9.5 Hz is clearly visible. (B) Frequencies identification via CF-PLM. The frequencies involved are f₁ = 11 Hz (for source 57) and f₂ = 1.5 Hz (for source 42). (C) Frequencies identification via Dual-Frequency Coherence. Although a dark area in case of f₁ = 10 Hz is visible, the two frequencies cannot be identified.

The same couple of regions has been analyzed also in case of the second subject. The PSD of the interferometric signal, which is reported in Figure 14, shows several peaks, although the one located at around 10 Hz is present also in this case. We have to underline that this analysis is done on a single link of one subject, so the amount of available data is very limited and the estimation noise is not minimal. The frequency identification step shows two couplings with the same Δf.

FIGURE 14

Figure 14. Analysis of signals from areas 57 (right inferior parietal lobule) and 42 (orbital part of the right superior frontal gyrus) of the second subject. (A) PSD of the interferometric signal. The peak related to the cross-frequency coupling, located at +10.5 Hz is visible. (B) Frequencies identification via CF-PLM. Two frequency combinations can be identified.

A second pair of brain regions has been considered, namely the left superior frontal gyrus and the left calcarine sulcus. The PSD of the interferometric signal is reported in Figure 15A. Two peaks are evident in this case, one centered in zero, related to the iso-frequency coupling, and another one at -8 Hz, which denotes a cross-frequency connectivity. By focusing on the latter, the identified involved frequencies are reported in Figure 15B, and are around 1 and 9 Hz. Also in this case, results are more convincing than the DFC (Figure 15C), which cannot be exploited for the frequency identification.

FIGURE 15

Figure 15. Analysis of signals from areas 7 (left superior frontal gyrus) and 25 (left calcarine sulcus). (A) PSD of the interferometric signal. Two main peaks are visible, located at 0 Hz (related to the iso-frequency coupling) and -8 Hz (related to the cross-frequency coupling). (B) Frequencies identification via CF-PLM. The frequencies involved are f₁ = 1 Hz (for source 7) and f₂ = 9 Hz (for source 25). (C) Frequencies identification via Dual-Frequency Coherence. A global maximum is not visible, and thus the two involved frequencies cannot be identified.

Discussion

In this article, we present a novel phase-based metrics capturing the occurrence of cross-frequency synchronization in the resting-state brain. The main advance of this work lies in the fact that the proposed procedure detects cross-frequency synchronization reliably, without a priori hypothesis about the frequencies of the synchronized components.

It is important to notice that such a procedure lands itself nicely to study if and where cross-frequency synchronization is occurring in resting-state, when no specific task is taking place, and hence no hypothesis about the frequencies of CFC is available. Furthermore, this procedure only captures phase synchronization, since the amplitude does not affect the estimate (Tass et al., 1998). This is of particular relevance, provided that a number of mechanisms are believed to operate simultaneously in the brain in order to allow communication between neuronal populations operating at different frequencies (Jensen and Colgin, 2007), but n:m synchronization is the only neuronal mechanism by which two neuronal population can influence each other at the temporal accuracy of the fast-operating neuronal population (Palva et al., 2005; Le Van Quyen and Bragin, 2007). Hence, the results provided by our procedure are interpreted in terms of a defined neurophysiological mechanism (i.e., n:m synchronization), and are very noise-resilient while being entirely independent from the amplitude of the signals. The detection of cross-frequency coupling is specifically relevant taking into account the phenomenon of frequency mixing, i.e., the appearance of new frequencies in neuronal circuitry when incoming oscillations are non-linearly integrated (Haufler and Paré, 2019).

Rössler Attractors

Firstly, we simulated synthetic data by using two Rössler attractors, as they retain non-linear properties that are similar to the ones displayed by real M/EEG data (Rosenblum et al., 2000). In order to simulate CFC, we modified one of the two attractors by applying a frequency shift. Similarly to what is shown in our previous work presenting the PLM (Baselice et al., 2019), one can appreciate that the peak in the interferometric spectrum grows monotonically as a function of the strength of synchronization between oscillators (regardless of their frequencies). Furthermore, one can appreciate that the interferometric spectrum peaks are at the frequency corresponding to the shift that had been introduced. Hence, the PSD allows the estimation of both the intensity of the synchronization, as well as the difference between the frequencies of the involved signals. In addition, the resiliency to noise has been tested in this simulation, and a reliable estimate is possible also with realistic amount of noise. By looking at the curves of Figure 5, we can identify 10⁵ as the minimum number of samples for the CF-PLM convergence.

Kuramoto Oscillators

When dealing with real M/EEG signals, the case is more complicated since each signal has a very rich frequency spectrum, where the simultaneous presence of multiple components synchronizing in iso- and cross- frequency occurs. We used a model based on Kuramoto oscillators to explore if the CF-PLM can disentangle these different contributions.

Firstly, we explored the simple synchronization between two oscillators synchronized in iso-frequency (10–10 Hz). As shown by the peak centered in 0 in Figure 6A, the synchronization is correctly captured as expected. Then, we explored the case of cross-frequency synchronization of Kuramoto oscillators. The interferometric spectrum displays one peak in correspondence to the frequency difference of the two oscillators (at -7 Hz in the example in Figure 6B, since the two originating signals are oscillating at 10 and 17 Hz, respectively). In the third case, one oscillator at 10 Hz has been compared to the sum of two more oscillators at 10 and 17 Hz. This simulation is intended to create a single signal where some specific components are synchronized in iso-frequencies, while different ones are synchronized in cross-frequency. As expected, the components in iso-frequency produced a peak at 0 Hz in the interferometric spectrum, while a second peak appears at -7 Hz, capturing the cross-frequency synchronization. Such results show that the proposed methodology can disentangle the cases where multiple components are synchronized simultaneously in iso- and cross-frequency. Similarly to the previous scheme, the CF-PLM produces a noise-resilient estimate of synchronization (despite being based solely on the phase).

Results reported in Figure 8 show that the proposed approach effectively estimates the frequencies involved in the coupling, both iso-frequency (Figure 8A), cross-frequency (Figure 8B) and simultaneous iso and cross-frequency (Figure 8C). More in details, the centers of all the computed crosses are correctly positioned and allow the identification of the frequency components present in the oscillators and involved in the coupling process.

Real Data

The analysis with the real data is intended to show that such a cross-frequency-coupling is happening in resting state, and can be captured by the proposed procedure. Each interferometric peak has been statistically validated against surrogates, making it unlikely that it would occur by chance. Furthermore, if the peaks were occurring by chance, one would expect that no consistent patterns of cross-frequency coupling in different subjects. However, observing the CFC analysis, one can appreciate the cross-frequency-coupling is not a widespread phenomenon happening in any area, but, rather, specific to some areas. More in detail, the first row of Figure 11, clearly shows that the strongest cross-frequency connections do not spread randomly across the matrix, and a texture appears, indicating that a specific CFC network is operating involving specific edges, rather than being randomly spread across the brain, as one would expected for a random phenomenon. Besides this coherent spatial distribution, images in the second row of Figure 11 show that also the frequency components that are correlated are not random, but a pattern emerges. Finally, Figure 12 helps the visualization of the regions where cross frequency coupling occurs consistently. While a systematic description of these patterns goes beyond the scope of the current paper, one should notice that temporo-parietal regions, as well as occipital ones, appear symmetrically involved in cross-frequency communication. Roughly speaking, these regions are involved in perceptive streams processing external stimuli. These findings are in line with previous evidence showing delta-to-alpha cross-synchronization (Isler et al., 2008). Furthermore, delta-to-alpha synchronization was also shown to be influenced by anesthetics (Stankovski et al., 2016), stressing their physiological meaningfulness. Importantly, as the ground truth is ultimately unknown, one should be very cautious at making inference.

Interestingly, in Shahbazi Avarvand et al. (2018) a 10 to 20 Hz coupling has been found in case of EEG data in the sensor space, probably due to alpha band harmonics, but the proposed result did not give this result. This suggests that the two methodologies handle differently the signal harmonics.

However, one good aspect of the proposed procedure lies in the fact that a form of “double-check” is possible. In fact, once the frequency on the peak of the interferometric spectrum is known, and the related component in signal A has been identified, one has already a hypothesis about the frequency of the synchronized component in signal B. Hence, if the procedure is consistent, the filtering of signal B should confirm this hypothesis. In fact, as explained previously, the peak of the interferometric signal appears at a frequency equal to the difference of frequencies between the two components. In conclusion, with the proposed procedure we are able to determine the central frequencies of signals A and B involved in the coupling. We made a number of tests to explore the behavior of the proposed procedure to different lengths of acquisitions and different number of epochs, confirming that this procedure is robust also with little and/or noisy data. One aspect of interest is that our metrics does not need the data to be split into trials, hence taking advantage of the full length of the available data. However, even when only short or limited data segments are available, the new procedure can still retrieve reliable results.

Conclusions

In this manuscript, we propose a new metric that can estimate cross-frequency coupling from broad-band signals with no a priori hypothesis on what the information transfers would be. Since cross-frequency coupling is the only neuronal mechanism that can allow fast communication between neuronal ensemble operating at different frequencies, we believe our metric can help to study the mechanisms of cross-frequency communication in the resting-state, as well as its topography and topology.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Ethics Statement

The studies involving human participants were reviewed and approved by Comitato Etico Campania Centro Prot.n.93C.E./Reg. n.14-17OSS. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

PS, MB, and FB conceived the experiments. RR, MA, and JC conducted the experiments. MA and LG analyzed the results. PS, MA, and FB conceived and developed the methodology. All authors reviewed the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

1. ^Ethics approval by Comitato Etico Campania Centro, Prot.n.93C.E./Reg. n.14-17OSS.

References

Auger, F., and Flandrin, P. (1995). Improving the readability of time-frequency and time-scale representations by the reassignment method. IEEE Trans. Signal Process. 43, 1068–1089.

Google Scholar

Barbati, G., Porcaro, C., Zappasodi, F., Rossini, P. M., and Tecchio, F. (2004). Optimization of an independent component analysis approach for artifact identification and removal in magnetoencephalographic signals. Clin. Neurophysiol. 115, 1220–1232. doi: 10.1016/j.clinph.2003.12.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Baselice, F., Sorriso, A., Rucco, R., and Sorrentino, P. (2019). Phase linearity measurement: a novel index for brain functional connectivity. IEEE Trans. Med. Imag. 38, 873–882. doi: 10.1109/TMI.2018.2873423

PubMed Abstract | CrossRef Full Text | Google Scholar

Bressler, S. L. (1995). Large-scale cortical networks and cognition. Brain Res. Brain Res. Rev. 20, 288–304.

Google Scholar

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Hugues, E., Sporns, O., and Deco, G. (2011). Role of local network oscillations in resting-state functional connectivity. NeuroImage 57, 130–139. doi: 10.1016/j.neuroimage.2011.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Luckhoo, H., Woolrich, M., Joensson, M., Mohseni, H., Baker, A., et al. (2014). Exploring mechanisms of spontaneous functional connectivity in MEG: how delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. NeuroImage 90, 423–435. doi: 10.1016/j.neuroimage.2013.11.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., Edwards, E., Dalal, S. S., Soltani, M., Nagarajan, S. S., Kirsch, H. E., et al. (2006). High gamma power is phase-locked to theta oscillations in human neocortex NIH public access. Science 313, 1626–1628. doi: 10.1126/science.1128115

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., and Knight, R. T. (2010). The functional role of cross-frequency coupling. Trends Cogn. Sci. 14, 506–515. doi: 10.1016/j.tics.2010.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, M. Y., Kim, H. J., Kim, D., and Hong, H. (2000). Synchronization in a system of globally coupled oscillators with time delay. Phys. Rev. E 61, 371–381. doi: 10.1103/PhysRevE.61.371

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, M. X. (2008). Assessing transient cross-frequency coupling in EEG data. J. Neurosci. Methods 168, 494–499. doi: 10.1016/j.jneumeth.2007.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, M. X., Elger, C. E., and Fell, J. (2009). Oscillatory activity and phase-amplitude coupling in the human medial frontal cortex during decision making. J. Cogn. Neurosci. 21, 390–402. doi: 10.1162/jocn.2008.21020

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgin, L. L., Denninger, T., Fyhn, M., Hafting, T., Bonnevie, T., Jensen, O., et al. (2009). Frequency of gamma oscillations routes flow of information in the hippocampus. Nature 462, 353–357. doi: 10.1038/nature08573

PubMed Abstract | CrossRef Full Text | Google Scholar

Darvas, F., Ojemann, J. G., and Sorensen, L. B. (2009). Bi-phase locking-a tool for probing non-linear interaction in the human brain. NeuroImage 46, 123–32. doi: 10.1016/j.neuroimage.2009.01.03

PubMed Abstract | CrossRef Full Text | Google Scholar

Doesburg, S. M., Green, J. J., McDonald, J. J., and Ward, L. M. (2009). Rhythms of consciousness: binocular rivalry reveals large-scale oscillatory network dynamics mediating visual perception. PLoS ONE 4, e6142. doi: 10.1371/journal.pone.0006142

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A. K., Gerloff, C., Hilgetag, C. C., and Nolte, G. (2013). Review intrinsic coupling modes : multiscale interactions in ongoing brain activity. Neuron 80, 867–886. doi: 10.1016/j.neuron.2013.09.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Fell, J., and Axmacher, N. (2011). The role of phase synchronization in memory processes. Nat. Rev. Neurosci. 12, 105–118. doi: 10.1038/nrn2979

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J. (1994). Functional and effective connectivity in neuroimaging: a synthesis. Hum. Brain Map. 2, 56–78. doi: 10.1002/hbm.460020107

CrossRef Full Text | Google Scholar

Gong, G., He, Y., Concha, L., Lebel, C., Gross, D. W., Evans, A. C., et al. (2009). Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cereb. Cortex (New York, N.Y. : 1991) 19, 524–536. doi: 10.1093/cercor/bhn102

PubMed Abstract | CrossRef Full Text | Google Scholar

Gross, J., Baillet, S., Barnes, G. R., Henson, R. N., Hillebrand, A., Jensen, O., et al. (2013). Good practice for conducting and reporting MEG research. NeuroImage 65, 349–363. doi: 10.1016/j.neuroimage.2012.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Haufler, D., and Paré, D. (2019). Detection of multiway gamma coordination reveals how frequency mixing shapes neural dynamics. Neuron 101, 603–614.e6. doi: 10.1016/J.NEURON.2018.12.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillebrand, A., Tewarie, P., van Dellen, E., Yu, M., Carbo, E. W. S., Douw, L., et al. (2016). Direction of information flow in large-scale resting-state networks is frequency-dependent. Proc. Natl. Acad. Sci. U.S.A. 113, 3867–3872. doi: 10.1073/pnas.1515657113

PubMed Abstract | CrossRef Full Text | Google Scholar

Isler, J. R., Grieve, P. G., Czernochowski, D., Stark, R. I., and Friedman, D. (2008). Cross-frequency phase coupling of brain rhythms during the orienting response. Brain Res. 1232, 163–172. doi: 10.1016/j.brainres.2008.07.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, O., and Colgin, L. L. (2007). Cross-frequency coupling between neuronal oscillations. Trends Cogn. Sci. 11, 267–269. doi: 10.1016/j.tics.2007.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V., and Müller, V. (2013). Cross-frequency coupling in real and virtual brain networks. Front. Comput. Neurosci. 7, 78. doi: 10.3389/fncom.2013.00078

PubMed Abstract | CrossRef Full Text | Google Scholar

Kendrick, K. M., Zhan, Y., Fischer, H., Nicol, A. U., Zhang, X., and Feng, J. (2011). Learning alters theta amplitude, theta-gamma coupling and neuronal synchronization in inferotemporal cortex. BMC Neurosci. 12, 55. doi: 10.1186/1471-2202-12-55

PubMed Abstract | CrossRef Full Text | Google Scholar

Kralemann, B., Cimponeriu, L., Rosenblum, M., Pikovsky, A., and Mrowka, R. (2008). Phase dynamics of coupled oscillators reconstructed from data. Phys. Rev. 77, 066205. doi: 10.1103/PhysRevE.77.066205

PubMed Abstract | CrossRef Full Text | Google Scholar

Lancaster, G., Iatsenko, D., Pidde, A., Ticcinelli, V., and Stefanovska, A. (2018). Surrogate data for hypothesis testing of physical systems. Phys. Rep. 748, 1–60. doi: 10.1016/j.physrep.2018.06.001

CrossRef Full Text | Google Scholar

Le Van Quyen, M., and Bragin, A. (2007). Analysis of dynamic brain oscillations: methodological advances. Trends Neurosci. 30, 365–373. doi: 10.1016/J.TINS.2007.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Mellors, R. J., Vernon, F. L., and Thomson, D. J. (1998). Detection of dispersive signals using multitaper dual-frequency coherence. Geophys. J. Int. 135, 146–154. doi: 10.1046/j.1365-246X.1998.00614.x

CrossRef Full Text | Google Scholar

Niebur, E., Schuster, H. G., and Kammen, D. M. (1991a). Collective frequencies and metastability in networks of limit-cycle oscillators with time delay. Phys. Rev. Lett. 67, 2753–2756. doi: 10.1103/PhysRevLett.67.2753

PubMed Abstract | CrossRef Full Text | Google Scholar

Niebur, E., Schuster, H. G., and Kammen, D. M. (1991b). Collective frequencies and metastability in networks of limit-cycle oscillators with time delay. Phys. Rev. Lett. 67, 2753–2756.

PubMed Abstract | Google Scholar

Nikulin, V., and Brismar, T. (2006). Phase synchronization between alpha and beta oscillations in the human electroencephalogram. Neuroscience 137, 647–657. doi: 10.1016/j.neuroscience.2005.10.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Nolte, G. (2003). The magnetic lead field theorem in the quasi-static approximation and its use for magnetoencephalography forward calculation in realistic volume conductors. Phys. Medi. Biol. 48, 3637–3652. doi: 10.1088/0031-9155/48/22/002

PubMed Abstract | CrossRef Full Text | Google Scholar

Oostenveld, R., Fries, P., Maris, E., and Schoffelen, J.-M. (2011). FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput. Intell. Neurosci. 2011, 156869. doi: 10.1155/2011/156869

PubMed Abstract | CrossRef Full Text | Google Scholar

Palva, J. M., Palva, S., and Kaila, K. (2005). Phase synchrony among neuronal oscillations in the human cortex. J. Neurosci. 25, 3962–3972. doi: 10.1523/JNEUROSCI.4250-04.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenblum, M., Pikovsky, A., Kuhn, A. A., and Busch, J. L. (2021). Real time estimation of phase and amplitude with application to neural data. Sci. Rep. 11, 18037. doi: 10.1038/s41598-021-97560-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenblum, M., Pikovsky, A., Kurths, J., Schäfer, C., and Tass, P. (2001). “Chapter 9 phase synchronization: from theory to data analysis,” in Neuro-Informatics and Neural Modelling, Handbook of Biological Physics, vol. 4, eds F. Moss, and S. Gielen (North-Holland), 279–321.

Google Scholar

Rosenblum, M. G., Tass, P. A., and Kurths, J. (2000). Estimation of Synchronization From Noisy Data With Application to Human Brain Activity. (Berlin: Springer), 202–211.

Google Scholar

Sadasivan, P. K., and Narayana Dutt, D. (1996). SVD based technique for noise reduction in electroencephalographic signals. Signal Process. 55, 179–189.

Google Scholar

Sauseng, P., Klimesch, W., Heise, K. F., Gruber, W. R., Holz, E., Karim, A. A., et al. (2009). Brain oscillatory substrates of visual short-term memory capacity. Curr. Biol. 19, 1846–1852. doi: 10.1016/j.cub.2009.08.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Schack, B., and Weiss, S. (2005). Quantification of phase synchronization phenomena and their importance for verbal memory processes. Biol. Cybern. 92, 275–287. doi: 10.1007/s00422-005-0555-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Seymour, R. A., Rippon, G., and Kessler, K. (2017). The detection of phase amplitude coupling during sensory processing. Front. Neurosci. 11, 487. doi: 10.3389/fnins.2017.00487

PubMed Abstract | CrossRef Full Text | Google Scholar

Shahbazi Avarvand, F., Bartz, S., Andreou, C., Samek, W., Leicht, G., Mulert, C., et al. (2018). Localizing bicoherence from eeg and meg. NeuroImage 174, 352–363. doi: 10.1016/j.neuroimage.2018.01.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Siebenhühner, F., Wang, S. H., Palva, J. M., and Palva, S. (2016). Cross-frequency synchronization connects networks of fast and slow oscillations during visual working memory maintenance. eLife 5, e13451. doi: 10.7554/eLife.13451

PubMed Abstract | CrossRef Full Text | Google Scholar

Sigl, J. C., and Chamoun, N. G. (1994). An introduction to bispectral analysis for the electroencephalogram. J. Clin. Monitor. 10, 392–404.

PubMed Abstract | Google Scholar

Sorrentino, P., Rucco, R., Jacini, F., Trojsi, F., Lardone, A., Baselice, F., et al. (2018). Brain functional networks become more connected as amyotrophic lateral sclerosis progresses: a source level magnetoencephalographic study. NeuroImage Clin. 20, 564–571. doi: 10.1016/j.nicl.2018.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Stankovski, T., Petkoski, S., Raeder, J., Smith, A. F., McClintock, P. V., and Stefanovska, A. (2016). Alterations in the coupling functions between cortical and cardio-respiratory oscillations due to anaesthesia with propofol and sevoflurane. Philosoph. Trans. R. Soc. Math. Phys. Eng. Sci. 374, 20150186. doi: 10.1098/rsta.2015.0186

PubMed Abstract | CrossRef Full Text | Google Scholar

Tass, P., Rosenblum, M. G., Weule, J., Kurths, J., Pikovsky, A., Volkmann, J., et al. (1998). Detection of n : m phase locking from noisy data: application to magnetoencephalography. Phys. Rev. Lett. 81, 3291–3294. doi: 10.1103/PhysRevLett.81.3291

CrossRef Full Text | Google Scholar

Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proc. IEEE 70, 1055–1096.

Google Scholar

Tort, A. B. L., Komorowski, R., Eichenbaum, H., and Kopell, N. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. J. Neurophysiol. 104, 1195–1210. doi: 10.1152/jn.00106.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Tort, A. B. L., Komorowski, R. W., Manns, J. R., Kopell, N. J., and Eichenbaum, H. (2009). Theta-gamma coupling increases during the learning of item-context associations. Proc. Natl. Acad. Sci. 106, 20942–20947. doi: 10.1073/pnas.0911331106

PubMed Abstract | CrossRef Full Text | Google Scholar

Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., et al. (2002). Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage 15, 273–289. doi: 10.1006/nimg.2001.0978

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Veen, B. D., Van Drongelen, W., Yuchtman, M., and Suzuki, A. (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44, 867–880.

PubMed Abstract | Google Scholar

Vanhatalo, S., Palva, J. M., Holmes, M. D., Miller, J. W., Voipio, J., and Kaila, K. (2004). Infraslow oscillations modulate excitability and interictal epileptic activity in the human cortex during sleep. Proc. Natl. Acad. Sci. 101, 5053–5057. doi: 10.1073/pnas.0305375101

PubMed Abstract | CrossRef Full Text | Google Scholar

Varela, F., Lachaux, J.-P., Rodriguez, E., and Martinerie, J. (2001). The brainweb: phase synchronization and large-scale integration. Nat. Rev. Neurosci. 2, 229–239. doi: 10.1038/35067550

PubMed Abstract | CrossRef Full Text | Google Scholar

Volk, D., Dubinin, I., Myasnikova, A., Gutkin, B., and Nikulin, V. V. (2018). Generalized cross-frequency decomposition: a method for the extraction of neuronal components coupled at different frequencies. Front. Neuroinf. 12, 72. doi: 10.3389/fninf.2018.00072

PubMed Abstract | CrossRef Full Text | Google Scholar

Wacker, M., and Witte, H. (2011). On the stability of the n:m phase synchronization index. IEEE Trans. Biomed. Eng. 58, 332–338. doi: 10.1109/TBME.2010.2063028

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, M., Engels, M. M. A., Hillebrand, A., van Straaten, E. C. W., Gouw, A. A., Teunissen, C., et al. (2017). Selective impairment of hippocampus and posterior hub areas in Alzheimer's disease: an MEG-based multiplex network study. Brain 140, 1466–1485. doi: 10.1093/brain/awx050

PubMed Abstract | CrossRef Full Text | Google Scholar