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 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 crossfrequency 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, phasephase 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, phasephase coupling.However, bicoherence is not a pure phasebased 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 crossfrequency 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 isofrequency 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 crossfrequency 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.

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: y an (t) = A y (t)e iφ y (t) = y R (t) + iy I (t) . (1) 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: 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).
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: 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).
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.: 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: 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 isofrequency 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.: 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).

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

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): 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 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 stopband 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 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 crossfrequency 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(ω) and Y(ω), the DFC can be defined as Shahbazi Avarvand et al. (2018): 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: in which the time series have been divided in K segments.We computed the DFC for ω 1 and ω 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.

Acquisition and Preprocessing
The acquisitions used for the analysis are from healthy subjects acquired at the MEG facility in Naples1 .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.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 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 | 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.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 ≈ 11 Hz and of f 2 ≈ 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 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   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 restingstate 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 5 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 isoand 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 temporoparietal 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.

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

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.

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.

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.

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

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.

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

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

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.

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