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

Pierpaolo Sorrentino, Michele Ambrosanio, Rosaria Rucco, Joana Cabral, Leonardo Gollo, Michael Breakspear, Fabio Baselice

► To cite this version:

Pierpaolo Sorrentino, Michele Ambrosanio, Rosaria Rucco, Joana Cabral, Leonardo Gollo, et al. Detection of Cross-Frequency Coupling Between Brain Areas: An Extension of Phase Linearity Measurement. Frontiers in Neuroscience, 2022, 16, 10.3389/fnins.2022.846623. [https://doi.org/10.3389/fnins.2022.846623](https://doi.org/10.3389/fnins.2022.846623). [CC-BY](https://creativecommons.org/licenses/by/4.0/)

HAL Id: hal-03786867

https://hal.science/hal-03786867

Submitted on 23 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

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

Pierpaolo Sorrentino<sup>1,2</sup>, Michele Ambrosanio<sup>3</sup>, Rosaria Rucco<sup>2</sup>, Joana Cabral<sup>4,5</sup>, Leonardo L. Gollo<sup>6,7</sup>, Michael Breakspear<sup>7,8</sup> and Fabio Baselice<sup>3*</sup>

**OPEN ACCESS**

**Edited by:**
Pedro Antonio Valdes-Sosa,
University of Electronic Science and
Technology of China, China

**Reviewed by:**
Laura Marzetti,
University of Studies G. d'Annunzio
Chieti and Pescara, Italy
Guido Nolte,
University Medical Center
Hamburg-Eppendorf, Germany

*Correspondence:
Fabio Baselice
fabio.baselice@uniparthenope.it

**Specialty section:**
This article was submitted to
Brain Imaging Methods,
a section of the journal
Frontiers in Neuroscience

Received: 31 December 2021

Accepted: 21 February 2022

Published: 25 April 2022

**Citation:**
Sorrentino P, Ambrosanio M, Rucco R, Cabral J, Gollo LL, Breakspear M and Baselice F (2022) Detection of Cross-Frequency Coupling Between Brain Areas: An Extension of Phase Linearity Measurement. *Front. Neurosci.* 16:846623.
doi: 10.3389/fnins.2022.846623

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.

**Keywords:** cross frequency coupling, brain network, brain functional connectivity, phase coupling, phase linearity measurement, PLM

## 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

<footnote n="1"><sup>1</sup> Systems Neuroscience Institute, Marseille, France, <sup>2</sup> Hermitage Capodimonte Hospital, Naples, Italy, <sup>3</sup> Engineering Department, University of Naples Parthenope, Naples, Italy, <sup>4</sup> Life and Health Sciences Research Institute (ICVS), University of Minho, Braga, Portugal, <sup>5</sup> Department of Psychiatry, University of Oxford, Oxford, United Kingdom, <sup>6</sup> Turner Institute for Brain and Mental Health, Monash University, Melbourne, VIC, Australia, <sup>7</sup> QIMR Berghofer Medical Research Institute, Brisbane, QLD, Australia, <sup>8</sup> Hunter Medical Research Institute, Newcastle, NSW, Australia</footnote>
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 neurophysiological 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{aligned}
x_{an}(t) &= A_x(t)e^{i\phi_x(t)} = x_R(t) + ix_I(t), \\
y_{an}(t) &= A_y(t)e^{i\phi_y(t)} = y_R(t) + iy_I(t).
\end{aligned} 
\quad (1) $$

where variables $A$ and $\phi$ 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:

$$ z(t) = \frac{x_{an}(t)y'_{an}(t)}{|x_{an}(t)|y_{an}(t)|} = e^{i\Delta\phi(t)}, \quad (2) $$

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 $\Delta\phi(t) = \phi_x(t) - \phi_y(t) \in [-\pi, \pi[$, 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 $\Delta\phi(t)$, instead of the genuine phases $\phi_x(t)$ and $\phi_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 $\Delta\phi(t)$ can be exploited in order to measure phase connectivity between $\phi_x(t)$ and $\phi_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 $[-\pi, \pi[$ range (blue line in Figure 1). In case of phase coupling, the term $\Delta\phi(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).

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:

$$ S_Z(f) = \left| \int_0^T z(t)e^{-i2\pi ft} dt \right|^2, \quad (3) $$

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 (blue line of Figure 2) if its phase term $\Delta\phi(t)$ irregularly spreads in the $[-\pi, \pi[$ range (blue line of Figure 1).

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

$$ \Delta\phi(t) = \varphi_{xy} + \Delta ft, \quad (4) $$

<figure><img src="image_3.png" /><figcaption>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).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

where the term $\Delta 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$:

$$PLM = \frac{\int_{-B}^{B} S_Z(f) df}{\int_{-\infty}^{\infty} S_Z(f) df} \qquad (5)$$

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 $\Delta 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 = \Delta 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 $\Delta f$ would solve the problem, as the integration could be shifted accordingly into the $[\Delta f - B, \Delta 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

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

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 $\Delta f$ is identified. Subsequently, the coupling strength is measured by adapting the upper integral of Equation (5), i.e.:

$$CF-PLM = \frac{\int_{\Delta f - B}^{\Delta f + B} S_Z(f) df}{\int_{-\infty}^{\infty} S_Z(f) df} \qquad (6)$$

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 $[\Delta f - B, \Delta 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).

## 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{aligned}
\dot{x}_{1,2} &= -2\pi f_{1,2}y_{1,2} - z_{1,2} + \xi_{1,2} + c(x_{2,1} - x_{1,2}) \\
\dot{y}_{1,2} &= 2\pi f_{1,2}x_{1,2} + 0.15y_{1,2} \\
\dot{z}_{1,2} &= 0.2 + z_{1,2}(x_{1,2} - 10)
\end{aligned}
\quad (7) $$

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_1$ 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_1$. 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:

$$ y_{\Delta f}(t) = y(t)e^{-i2\pi \Delta ft}. \quad (8) $$

For this analysis, we considered $\Delta f = 7$ Hz. The CF-PLM has been computed between $x(t)$ and $y_{\Delta 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.

### Kuramoto Oscillators

Three mutually coupled Kuramoto oscillators, namely $s_1(t)$, $s_2(t)$, and $s_3(t)$, have been generated according to the following model (Cabral et al., 2011):

$$ \frac{d\theta_n}{dt} = 2\pi f_n + k \sum_{p=1}^{3} C\sin(\theta_p(t - \tau) - \theta_n(t)) \quad (9) $$

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>
with $n = 1, 2, 3$, $\tau = 0.6$ s and $C = 1$. The central frequencies have been set equal to 10, 10, and 17 Hz, respectively. The first oscillator ($s_1(t), f = 10$ Hz) has been compared to the second one ($s_2(t), f = 10$ Hz), to the third one ($s_3(t), f = 17$ Hz) and to the sum of the last two ($s_2(t) + s_3(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_1(t)$ and $s_3(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.

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

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

<figure><img src="image_7.png" /><figcaption>**FIGURE 5** | CF-PLM values measured in case of two Rössler oscillators varying the number of samples.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**FIGURE 7** | Position of the interferometric peak while varying the global coupling strength of two Kuramoto oscillators.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**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).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

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_1(t)$, $s_2(t)$ coupling, at (10, 17) Hz in the case of $s_1(t)$, $s_3(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.

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(\omega)$ and $Y(\omega)$, the DFC can be defined as Shahbazi Avarvand et al. (2018):

$$DFC(\omega_1, \omega_2) = \frac{|\langle X(\omega_1)Y'(\omega_2) \rangle|}{\langle |X(\omega_1)|^2 \rangle^{1/2} \langle |Y(\omega_2)|^2 \rangle^{1/2}}, \quad (10)$$

in which the symbol $\langle \cdot \rangle$ 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:

$$PSI(n, m) = \frac{1}{N} \left| \sum_{k=1}^{K} \exp[i(n\phi_x^k - m\phi_y^k)] \right|, \quad (11)$$

in which the time series have been divided in K segments. We computed the DFC for $\omega_1$ and $\omega_2$ 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.

## Real Data

### Acquisition and Preprocessing

The acquisitions used for the analysis are from healthy subjects acquired at the MEG facility in Naples<fnref n="1" />. 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

<footnote n="1">¹Ethics approval by Comitato Etico Campania Centro, Prot.n.93C.E./Reg. n.14-17OSS.</footnote>
<figure><img src="image_11.png" /><figcaption>**FIGURE 9** | Dual-Frequency Coherence (left) and n:m Phase Synchronization Index (right) results in case of CFC Kuramoto oscillators.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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 ($|\Delta f| \ge 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.

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
<figure><img src="image_13.png" /><figcaption>**FIGURE 11** | First row: peak intensity, second row: peak frequency from CF-PLM matrix of 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 characterize the CF-PLM peaks with higher intensity.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

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_1 \approx 11$ Hz and of $f_2 \approx 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.

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 $\Delta f$.

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
A

<figure><img src="image_15.png" /></figure>

B

<figure><img src="image_16.png" /></figure>

C

<figure><img src="image_17.png" /></figure>

**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_1 = 11$ Hz (for source 57) and $f_2 = 1.5$ Hz (for source 42). (C) Frequencies identification via Dual-Frequency Coherence. Although a dark area in case of $f_1 = 10$ Hz is visible, the two frequencies cannot be identified.

A

<figure><img src="image_18.png" /></figure>

B

<figure><img src="image_19.png" /></figure>

**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

<figure><img src="image_20.png" /></figure>

B

<figure><img src="image_21.png" /></figure>

C

<figure><img src="image_22.png" /></figure>

**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 = 1$ Hz (for source 7) and $f_2 = 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.
more convincing than the DFC (Figure 15C), which cannot be exploited for the frequency identification.

## 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 expect 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 temporal-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.

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

* 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

* 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

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

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

* 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

* 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

* 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

* 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

* 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

* 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

* 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

* 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
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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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

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

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

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.

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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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

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

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

Copyright © 2022 Sorrentino, Ambrosanio, Rucco, Cabral, Gollo, Breakspear and Baselice. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.