Commentary: Evaluation of Phase-Amplitude Coupling in Resting State Magnetoencephalographic Signals: Effect of Surrogates and Evaluation Approach

Using Human Connectome Project (HCP) data, Gohel et al. (2016) (GEA) recently reported that “resting-state MEG signals failed to exhibit ubiquitous phase-amplitude coupling (PAC) phenomenon, contrary to what has been suggested” by Florin and Baillet (2015) (FB). GEA argued that the original PAC findings by FB were driven by false positives resulting from the use of inappropriate methods. In this commentary, we first correct GEA’s mischaracterization of the approach actually used by FB. We then investigated the PAC computations in Gohel et al. (2016) (GEA) and demonstrate that with FB’s approach, it is actually possible to detect PAC in the Human Connectome Project (HCP) resting-state data. Finally, when making the data processing as similar as possible to GEA we still found significant PAC across large portions of the brain.


INTRODUCTION
Using Human Connectome Project (HCP) data, Gohel et al. (2016) (GEA) recently reported that "resting-state MEG signals failed to exhibit ubiquitous phase-amplitude coupling (PAC) phenomenon, contrary to what has been suggested" by Florin and Baillet (2015) (FB). GEA argued that the original PAC findings by FB were driven by false positives resulting from the use of inappropriate methods. In this commentary, we first correct GEA's mischaracterization of the approach actually used by FB. We then investigated the PAC computations in Gohel et al. (2016) (GEA) and demonstrate that with FB's approach, it is actually possible to detect PAC in the Human Connectome Project (HCP) resting-state data. Finally, when making the data processing as similar as possible to GEA we still found significant PAC across large portions of the brain. Florin and Baillet (2015) GEA claimed that FB's significance threshold was "generated after combining all surrogate PAC scores from all frequency pairs" and argued that such pooling "of PAC scores across the frequency pairs is not appropriate, " because it would bias the finding of significant couplings toward low frequencies. This claim is incorrect: The PAC detection threshold in FB was actually computed separately for each frequency pair (page 28 in FB: "the value at the 95th quantile of PAC scores for every low/high frequency pair was determined as the significance threshold at p < 0.05 for actual recordings").

Artifact Removal
From GEA's statement that at least 10 independent components were removed from the data, it seems that they used their own preprocessing pipeline, not the distributed preprocessed version of the HCP data, because this latter had fewer independent components removed. After direct communication with GEA, and in an attempt to reproduce their published method, we used the preprocessed HCP data set.

Source Modeling
GEA used a single-sphere head model, while FB used multiplesphere forward modeling. Even though the latter is best-practice in MEG analysis (Gross et al., 2013), we tried to reproduce for the purpose of this commentary the approach of GEA by performing analyses with a single-sphere head model. However, their arbitrary selection of 58 sources where they computed PAC measures was not described with sufficient details to enable the identical reproduction of their analyses. In addition, such a small sample of sources can be expected to have biased GEA's results.

Frequency Selection
In the FB original analysis we used all frequency pairs (f ϕ , f a ), with f ϕ ∈ [2, 48] Hz and f a ∈ [80, 150] Hz for the calculation of PAC. The binning width for f ϕ was 0.5 Hz from 2 to 12 Hz and 2 Hz from 12 to 48 Hz. The bins for f a were logarithmically spaced. Following Canolty et al. (2006), we used chirplets to extract the phase and amplitude with a chirp-factor of 0, hence yielding wavelet-like decompositions.
GEA chose three low frequencies (3, 6, 10 Hz) and as high frequencies the sum over two frequency ranges, 53:10:93 and 83:10:143 Hz. To extract phase and amplitude they used Morlet wavelets with a mother wavelet of full width at half maximum (FWHM) of 1s and a center frequency of 1 Hz. This choice leads to 46 Hz wide FWHM filter at 53 Hz. This results in considerable smearing of the high-frequency binning.

Statistical Detection of PAC
GEA used different statistical methods to determine statistical significance of PAC values. For the purpose of the present rebuttal, we only used one approach from GEA, which was also used in FB: We generated pink noise time series of the same lengths as the original data and calculated PAC on these surrogate data. We generated 500 pink noise time series as in GEA. For the determination of the significance threshold, GEA used a more conservative 99% quantile, while FB used 95%.

METHOD COMPARISON WITH HCP DATA
We used the preprocessed data from four HCP subjects (100307,108323,113922,877168). All subjects were included by GEA and we only used the one run described in the supplementary material of GEA. For preprocessing, as many ICA components were removed as described in the icaclass folder of the preprocessed data. The further analyses were conducted with Brainstorm (Tadel et al., 2011). For source reconstruction we used Fieldtrip's (Oostenveld et al., 2011) implementation to derive a single sphere head model and estimate source activations with a weighted minimum norm estimator. This is presumably the same approach as GEA used. We then computed PAC according to the method proposed by Özkurt et al. (2011) and used in FB and GEA at each of the 8,004 individual source locations.
As could be expected, using a higher significance threshold lead to fewer significant findings than in FB (Table 1). Comparing Chirplets to Morlets, fewer vertices with significant PAC were found with Morlet wavelets. Importantly the binning of low frequencies and summing of the gamma amplitudes as in GEA led to almost no PAC detection. Overall sparse sampling of the low frequencies had a greater influence on the detection of PAC than the sparse sampling in the gamma range.

DISCUSSION
We investigated the claim from GEA that the results of FB were only due to false positives. We did so by using four of the HCP datasets used by GEA, staying as close as possible to their approach.

NO SUMMING IN GAMMA FREQUENCY RANGE; LF ACCORDING TO GEA
Sign. sources 68.6 ± 7.7 21.4 ± 9.4 13.9 ± 9.4 5.4 ± 8.4 Coupling δ/θ − γ 55.5 ± 9.9 15.4 ± 9.2 13.2 ± 9.6 5.4 ± 8.4 Coupling δ − γ 37.3 ± 6.8 9.8 ± 6.7 9.8 ± 11.6 5.0 ± 8.6 LF FROM 3 TO 10 Hz WITH 0. methodology. We believe that not restricting PAC analysis to few frequencies pairs provides a more encompassing picture of potential oscillatory coupling in the brain. We used pink (1/f ) Gaussian noise for the determination of the significance threshold as in FB. GEA demonstrated that using pink Gaussian noise is less conservative than the non-colored noise approach. While this might be the case, pink noise is the preferred option, because neural signals exhibit a 1/f frequency spectrum and non-parametric statistical testing requires that one should target the original data properties of non inferential interest with surrogate data. Of note for future research is that the non-sinusoidal shape of brain oscillations might affect the detection of cross-frequency coupling depending on the method used to extract the phase and amplitude coupling data (Yeh et al., 2016;Cole and Voytek, 2017). In principle, chirplets and wavelets can be constructed to account for the expected non-sinusoidal shape of brain oscillations (Mallat, 2009).
In summary, we can only restate that ubiquitous PAC in resting state MEG data is present across the brain. This was also the case for data from the HCP repository if all potential low-frequency high-frequency pairs are tested exhaustively.

ETHICS STATMENT
We used publicly available, open access data from the human connectome project.