Neuronal Oscillations with Non-sinusoidal Morphology Produce Spurious Phase-to-Amplitude Coupling and Directionality

Neuronal oscillations support cognitive processing. Modern views suggest that neuronal oscillations do not only reflect coordinated activity in spatially distributed networks, but also that there is interaction between the oscillations at different frequencies. For example, invasive recordings in animals and humans have found that the amplitude of fast oscillations (>40 Hz) occur non-uniformly within the phase of slower oscillations, forming the so-called cross-frequency coupling (CFC). However, the CFC patterns might be influenced by features in the signal that do not relate to underlying physiological interactions. For example, CFC estimates may be sensitive to spectral correlations due to non-sinusoidal properties of the alpha band wave morphology. To investigate this issue, we performed CFC analysis using experimental and synthetic data. The former consisted in a double-blind magnetoencephalography pharmacological study in which participants received either placebo, 0.5 or 1.5 mg of lorazepam (LZP; GABAergic enhancer) in different experimental sessions. By recording oscillatory brain activity with during rest and working memory (WM), we were able to demonstrate that posterior alpha (8–12 Hz) phase was coupled to beta-low gamma band (20–45 Hz) amplitude envelope during all sessions. Importantly, bicoherence values around the harmonics of the alpha frequency were similar both in magnitude and topographic distribution to the cross-frequency coherence (CFCoh) values observed in the alpha-phase to beta-low gamma coupling. In addition, despite the large CFCoh we found no significant cross-frequency directionality (CFD). Critically, simulations demonstrated that a sizable part of our empirical CFCoh between alpha and beta-low gamma coupling and the lack of CFD could be explained by two-three harmonics aligned in zero phase-lag produced by the physiologically characteristic alpha asymmetry in the amplitude of the peaks relative to the troughs. Furthermore, we showed that periodic signals whose waveform deviate from pure sine waves produce non-zero CFCoh with predictable CFD. Our results reveal the important role of the non-sinusoidal wave morphology on state of the art CFC metrics and we recommend caution with strong physiological interpretations of CFC and suggest basic data quality checks to enhance the mechanistic understanding of CFC.

Introduction

Traditionally brain rhythms have been divided within different frequency bands: delta (1–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), beta (15–30 Hz), and gamma (>30 Hz). These bands have been linked to a variety of normal and pathological behavioral phenotypes (Steriade et al., 1990). Importantly, this classification seems to be well-preserved amongst species despite the large differences in brain size (Buzsáki et al., 2013). Recent studies suggest that different brain rhythms can interact forming the so called cross-frequency coupling (CFC) (Jensen and Colgin, 2007). For example, a seminal work found that broadband gamma power (80–200 Hz) locked to theta rhythm in humans under various tasks (Canolty et al., 2006). In rat hippocampus beta-low gamma power (20–40 Hz) seems to be enhanced is specific phase bins of ongoing theta (Igarashi et al., 2014). Invasive recordings in animals and humans have implicated CFC dynamics to a large variety of functions such as neuronal communication between brain regions (Colgin et al., 2009; von Nicolai et al., 2014) and cortical layers (Spaak et al., 2012), regulation of the firing rate of a particular brain area (Canolty et al., 2012), associative learning (Tort et al., 2009; Igarashi et al., 2014), working memory (WM) information maintenance (Siegel et al., 2009; Axmacher et al., 2010) and neuronal communication of dispersed functional cell assemblies (Canolty et al., 2010, 2012). The overarching idea behind CFC is that independent brain rhythms interact hierarchically: while low frequency oscillations would synchronize large neuronal ensembles, local populations could be synchronized by faster rhythms modulated by the phase of slower oscillations (Buzsaki, 2006). Therefore, it is important to characterize the CFC to understand the dynamics and functions of brain rhythms.

Despite the popularity of the CFC framework, recent studies demonstrated potential caveats and major pitfalls in the interpretation of the results. Kramer et al. (2008) demonstrated that sharp transients inserted in periodic signals produce spurious (i.e., non-zero CFC), which resembles broadband gamma activity. The broadband activity arises because sharp-edged transients in the time domain correspond to high-frequency broadband harmonic components in the frequency domain. This phenomenon has subsequently been reported in rat hippocampus (Tort et al., 2013) (see Figure 8C of the paper). The term spurious CFC refers to the correlations between multiple frequency components (i.e., terms of the Fourier decomposition) caused when the input signal contains non-sinusoidal properties. These non-sinusoidal features can be transient or periodic and result in spectral correlations that will produce a non-zero CFC, even when a true physiological interaction is absent. More recently, Aru et al. (2015) provide an exhaustive review about potential caveats and confounds in the current CFC methodology, such as non-stationarities that may lead to spurious CFCs due to spectral correlations. They describe the difficulty to generate surrogate data that effectively destroys the specific cyclo-stationarities associated to the CFC effect of interest, while retaining the original probability distributions of all unspecific effects (i.e., not associated to the hypothesized CFC). Although the authors do not provide conclusive solutions to deal with the caveats, they propose a number of sanity checks to be performed in order to improve the mechanistic interpretation of the CFC.

Considering the caveats reviewed above (Kramer et al., 2008; Tort et al., 2013; Aru et al., 2015), the phase-to-amplitude coupling may overlap with the spectral harmonics of the alpha frequency. Various studies have shown (Stam et al., 1999; Nikulin et al., 2007, 2010) or observed (Gastaut, 1952; Contreras et al., 1997) that the alpha rhythm does not appear to be a pure sinusoidal wave, with a stronger amplitude of the peaks relative to the troughs. These non-sinusoidal features present in alpha oscillations will be expressed as harmonic components. On that note, harmonic components are, by definition, phase-to-amplitude coupled: the phase of the fundamental is locked to the power of the harmonic(s). This raises the concern whether the CFC patterns recorded in a given sensor are a genuine interaction reflecting a mechanistic physiological process between two independent rhythms, or whether it might be a more trivial consequence of spectral correlations due to the non-sinusoidal morphology of the alpha oscillations.

To investigate the influence of the morphology of neuronal oscillation on CFC metrics, we analyze human magnetoencephalographic (MEG) recordings during resting state and WM. We performed a double-blind randomized cross-over design pharmacological challenge where healthy participants received either placebo or lorazepam (LZP; GABA_A allosteric modulator) of different dosages (0.5 and 1.5 mg) in separate sessions. Pharmacological power modulations during WM have been published elsewhere (Lozano-Soldevilla et al., 2014). In addition, we constructed synthetic data varying systematically the shape of the signal to study its consequences with CFC metrics such as cross-frequency coherence (CFCoh), cross-frequency directionality (CFD), and bicoherence.

Materials and Methods

Experimental Procedures

Participants (n = 25) gave informed consent approved by the local research ethics committee (Commissie Mensgebonden Onderzoek; Arnhem-Nijmegen, number: 2011/199, date: August 5, 2011) and were compensated for their participation. They were engaged in a classic delayed-match to sample visuo-spatial WM task (Vogel and Machizawa, 2004) (Figure S1A) while brain oscillatory activity was acquired with whole-head MEG. Each trial started by presenting a centrally presented visual cue for 1.5 s. The task required participants to encode a sample array composed of colored squares presented in the visual hemifield indicated by a spatial cue. After a 1.5 s delay interval, a probe array was presented and participants had to indicate whether it matched the sample array or not. Participants had to keep their fixation to the center of the screen while covertly attend to the cued hemifield during the whole trial. Before the WM task, 2 min eyes-closed (EC) and 2 min eyes-opened (EO) resting state recordings were acquired (the order was counterbalanced across sessions).

MEG Acquisition and Data Analysis

Specific information about MEG data acquisition can be found in Lozano-Soldevilla et al. (2014) and we describe here the fundamental information. Ongoing brain activity was recorded using a whole-head MEG system with 275 axial sampled at 1200 Hz. To minimize participants head movements while maximizing their comfort, we used an orthopedic neck-collar during the MEG recordings (Push braces, height 10, size 2, model 1.60.2.02; http://www.push.eu). We used the real-time head localizer tool (Stolk et al., 2013) to monitor participant's head position and orientation relative to the MEG sensors during all the recordings.

MEG data was low-pass filtered below 200 Hz and down-sampled to 600 Hz. Subsequently we analyzed our dataset using Fieldtrip (Oostenveld et al., 2011) (http://fieldtriptoolbox.org); an open-source toolbox developed at the Donders Institute for Brain, Cognition, and Behaviour (Nijmegen, The Netherlands) and custom Matlab code (MathWorks, Natick, MA, USA).

Sensor Analysis

Resting-state data

We eliminated line noise by fitting sine and cosine functions between 48–52, 98–102, 148–152 Hz and subsequently subtracting these estimated components. EC and EO recordings were cut in epochs of 120 s duration. Raw data was visually inspected and the pieces containing eye-blink, muscle activity, or sensor artifacts were rejected. Sensors with noisy periods were rejected and linearly interpolated (neighbor average). After partial artifact rejection, data segments shorter or equal to 2 s were excluded leaving the rest of the cleaned data for the main analysis. Remaining eye movements and electrocardiogram (ECG) activity were isolated by independent component analysis (ICA). The independent components (ICs) reflecting potential artifacts were semi-automatically identified by correlating its time-course with the respective vertical EOG and ECG bipolar derivations. Given the stable periodicity of the QRS complex sharp edges, spurious CFC could arise (Kramer et al., 2008). To reduce ECG residual activity after the removal of the identified ICs, we adapted a procedure proposed by Tal and Abeles (2013). Briefly, we band-pass filtered the ECG with a zero-phase two-pass 4th order Butterworth filter (0.5–50 Hz). Then we took the absolute value of the ECG signal and we found the QRS peaks separated by a minimum of 0.66 s. Once all the peaks were identified, we computed peak-triggered averages (–0.3 to 0.4 s) using the time points defined in the ECG bipolar channel of the unfiltered MEG data. The obtained QRS-template was subtracted from the MEG recordings cycle-by-cycle. We observed QRS residuals in most of the participants after ICA nonetheless. Given that it has been empirically demonstrated that the reliability of the Infomax algorithm correlated positively to the amount of data samples used during the decomposition and it correlated negatively to the number of sensors² (Groppe et al., 2009), it is likely that the Infomax did not have enough data points per sensor to perform a more accurate solution (i.e., we only obtained < 4 min recordings with 273 sensors).

Spectral oscillatory power was performed by means of Fast Fourier Transform (FFT) using a fixed window length of 4 s (4–40 Hz in 0.25 Hz steps), 50% overlap, multiplied by a Hanning taper. Prior to spectral analysis, the mean of each epoch was subtracted. Spectra were normalized based on the total power in each sensor within each participant and MEG session before group averaging.

To investigate cross-frequency interactions we chose the CFCoh index (Osipova et al., 2008). CFCoh measures the phase consistency between a low frequency signal and the amplitude envelope of higher frequencies. Let the signal be represented by the time series x₁, x₂, …, x_N where N is the number of time points acquired at specific sampling rate (Fs). First, the time-course of amplitude $y^{\nu }=\left(y_1^{\nu },y_2^{\nu },\dots ,y_N^{\nu }\right)$ was estimated for frequency v by applying a sliding tapered time-window followed by a Fourier transform:

$\begin{array}{ll}{y}_n^v=\frac{1}{Fs}|{{\displaystyle \sum }}_{m=1}^M{h}_m{x}_{n+m-M/2}{e}^{\frac{-i2\pi vm}{Fs}}|& (1)\end{array}$

The function h_m is a Hanning taper M data points-long equaling the length of the sliding time window. The length of the time-window decreases with frequency such that M = w · Fs/v, where w denotes the numbers of cycles per time window (five cycles in the present study). Note that the time-series x can be zero-padded to reduce boundary effects at the beginning and end of the segments. In the present analysis, low frequencies were defined between 5 and 25 Hz and higher frequencies from 5 to 90 Hz in steps of 1 Hz.

Next we proceeded to estimate the Fourier spectra of x and y. Data was segmented into a set of S segments 0.6 s long with 50% overlap. Using a standard FFT algorithm, the complex representations were calculated:

$\begin{array}{ll}{\mathbf{X}}^s=\mathrm{\text{FFT}}\left({\mathbf{h}}^T{\mathbf{x}}^s,n_{FFT}\right)& (2)\end{array}$

and

$\begin{array}{ll}{\mathbf{Y}}^{\upsilon ,s}=\text{ FFT}\left({\mathbf{h}}^T{\mathbf{y}}^{v,s},n_{FFT}\right)& (3)\end{array}$

with x^s representing the temporal evolution of raw signal of segments s, y^{v, s} representing the amplitude envelope at high frequency v and segment s and n_FFT representing the frequency resolution of the x-axis Δf = Fs/n_FFT of the CFCoh comodulogram (n_FFT = 1024). The output vector ${\mathbf{X}}^s=(X_1^s,X_2^s,\dots ,X_{n_{FFT}}^s)$ and ${\mathbf{Y}}^{v,s}=\left(Y_1^{v,s},Y_2^{v,s},\dots ,Y_{n_{FFT}}^{v,s}\right)$ were Fourier transforms centered at frequencies: $f\in \left\{0,\frac{Fs}{n_{FFT}},\dots ,\frac{Fs}{2}\left(1-\frac{1}{n_{FFT}}\right)\right\}$. The cross-spectra ξ of the individual segments was defined as ξ^{v, s} = X^s(Y^{ν, s})* where “*” denoted the complex conjugate.

To quantify the CFC, we used the magnitude squared CFCoh (η²):

$\begin{array}{ll}{\eta }^{\mathbf{2}}(v,f)=\frac{{\left|{\displaystyle {\sum }_{s=1}^S}{\xi }^{v,s}\right|}^2}{{\displaystyle {\sum }_{s=1}^S{\left|{\mathbf{X}}^s\right|}^2}·{\displaystyle {\sum }_{s=1}^S{\left|{\mathbf{Y}}^{v,s}\right|}^2}}& (4)\end{array}$

where frequencies were defined as $f\in \left\{0,\frac{Fs}{n_{FFT}},\dots ,\frac{Fs}{2}\left(1-\frac{1}{n_{FFT}}\right)\right\}$.

To assess the directionality of CFCoh, we employed the new CFD metric (Jiang et al., 2015). CFD is the CFC extension of the phase-slope index (PSI) (Nolte et al., 2008). It rests on the assumption that when the transmission between two signals (sender and receiver) is associated with a fixed delay, the phase difference between the signals will change systematically with frequency as observed in the phase spectrum (Gotman, 1983; Mima et al., 2000). The PSI estimates the slope of the phase difference as a function of frequency in a given band with the direction of interaction attached to its sign. We here applied the PSI to the signal x and the amplitude envelope of the signal y^v at frequency tile (v, f_j):

$\begin{array}{ll}\psi \left(v,f_j\right)=\Im \left({\displaystyle {\sum }_{f_j-\frac{\beta }{2}}^{f_j+\frac{\beta }{2}}{\mathbf{C}}^{\ast }\left(v,f_j\right)}\mathbf{C}\left(v,\left(f_j+\Delta f\right)\right)\right)& (5)\end{array}$

where the complex coherency was defined as:

$\begin{array}{ll}\mathbf{C}\left(v,f_j\right)=\frac{{\sum }_{s=1}^S{\mathbf{X}}^s{({\mathbf{Y}}^{v,s})}^{*}}{\sqrt{{\sum }_{s=1}^S{\left|{\mathbf{X}}^s\right|}^2{\sum }_{s=1}^S{\left|{\mathbf{Y}}^{v,s}\right|}^2}}& (6)\end{array}$

with Δf denoting the frequency resolution and ℑ() representing the imaginary part. We used β to denote the bandwidth for which the phase slope was calculated (β = 4 Hz). $\tilde{\Psi }$ corresponded to the weighted average of the slope and we normalized it by an estimate of its standard deviation through a jackknife procedure (i.e., Ψ = $\tilde{\Psi }$/std($\tilde{\Psi }$)). The authors considered absolute values larger than 2 as significant (Nolte et al., 2008).

In order to assess cross-frequency phase-to-phase coupling we computed the magnitude squared bicoherence index (Kim and Powers, 1979):

$\begin{array}{ll}{b}^2\left(v,f\right)=\frac{{\left|{\sum }_{s=1}^S{X}^s(v)X^s(f)X^{*s}(v+f)\right|}^2}{{\sum }_{s=1}^S{\left|X^s(v)X^s(f)\right|}^2{\sum }_{s=1}^S{\left|X^s(v+f)\right|}^2}& (7)\end{array}$

The bicoherence is close to 0 when the coherence of Fourier frequency triplets (v, f, and v + f) behave randomly across S segments and it is close to 1 when the phase of the numerator (bispectrum) remains constant over segments. As for CFCoh, we epoched the data using fixed window length of 0.6 s, 50% overlapping, multiplied by a Hanning taper.

To further explore the phasic relationship between beta-low gamma power and the alpha rhythm, we took two complementary approaches based on event-related averages. The first one was based on the computation of time-frequency representations (TFR) locked to the alpha peak (Canolty et al., 2006; Osipova et al., 2008). To archive that, the MEG axial gradiometers of each participant and pharmacological session were band-pass filtered within the alpha band (8–12 Hz; zero-phase two-pass 4th order Butterworth filter). Alpha peaks were identified and the unfiltered raw data was cut in 0.6 s epochs centered at the peak. TFRs of power were calculated for each single trial using same parameters as for the CFCoh (i.e., ${\left|y^v\right|}^2$). Induced power estimates were averaged across trials and normalized relative to the mean power per frequency. Complementarily, we also performed event-related averages taking the gamma peak-as a trigger (Bragin et al., 1995).

Resting state source analysis

To localize the LZP-related oscillatory power modulations during EC and EO epochs, we adopted a beamformer approach as defined in Gross et al. (2001). Alpha band source reconstruction was computed by applying a FFT approach (epochs of 0.6 s) centered at 10 Hz using a one Slepian taper resulting in 2 Hz smoothing. The source reconstruction of the beta band activity (14–30 Hz) was centered at 22 Hz and 9 orthogonal Slepian tapers resulting 8 Hz smoothing. The brain volume of each individual was divided into a grid with a 1 cm resolution and. For every grid point a spatial filter was constructed from the lead field and the cross-spectral density (CSD) matrix. The lead fields were calculated from a subject specific realistic single-shell model of the brain (Nolte, 2003). CSD matrices were computed between all MEG pairs for the time period and frequency of interest. Experimental conditions (EC, EO) were combined for the purpose of calculating a common spatial filter, after which the power at each grid point was estimated for each session separately in every participant. We estimated the orientation in which the power of each dipole grid was maximal according by single value decomposition (Schoffelen et al., 2011). Prior to averaging, the source estimates of the individual functional data were normalized onto the MNI brain template (International Consortium for Brain Mapping, Montreal Neurological Institute, Canada) using SPM8 (http://www.fil.ion.ucl.ac.uk/spm).

Working memory (WM) data

Preprocessing, artifact rejection, and sensor spectral analysis procedures were done as in Lozano-Soldevilla et al. (2014). Briefly, trials of 1 s (0.5–1.5 s after memory sample array) were extracted for each pharmacological session. CFCoh, CFD and bicoherence sensor analysis were performed on MEG axial gradiometers with same parameters as those ones used for the resting state data.

Statistical analysis

To statistically assess the drug main effect on oscillatory power, CFCoh, bicoherence, or peak-triggered averages at the sensors space, a cluster-based nonparametric permutation test was applied (Maris and Oostenveld, 2007). This approach provides an effective method for controlling for multiple comparisons.

A repeated-measures ANOVA approach was used to asses drug main effects on selected sensors of interest. P-values were calculated and Greenhouse-Geisser corrections were made to control of sphericity violation when appropriated. Subsequent post-hoc t-tests, in case of main effect and/or interaction, were always corrected by multiple comparisons using the false discovery rate approach (Benjamini and Hochberg, 1995).

Simulations

We constructed synthetic data to demonstrate the role of signal wave morphology on different CFC metrics. All simulations were generated using a sampling rate of 1000 Hz and 180 s duration.

Narrow Gamma Amplitude (60–80 Hz) to Alpha Phase (~10 Hz) CFC

We created a signal containing CFC between the alpha and gamma by adapting the simulation described in example 3 from Kramer et al. (2008):

$\begin{array}{ll}\alpha =cos\left(2\pi ft\right)& (8)\end{array}$

where the α frequency f randomly varied cycle-by-cycle with the form f = 1/(0.1 s ± 0.01 s). In other words, we introduced small temporal variations in its period. The variable t defined the time axis running from 0 to 1 s. To each cycle of the oscillation we add a short 50 ms burst of high frequency (60–80 Hz) white noise that was tapered with a Hanning window. The Hanning window was scaled to a maximum amplitude of 5. The onset of the gamma signal relative to the alpha cycle varied uniformly between 0 and 20°C. Finally, normally distributed random noise with zero mean and standard deviation 0.1 were added.

Narrow-Band Gamma Amplitude (60–80 Hz) to Alpha Phase CFC with Non-sinusoidal Wave Morphology Properties

We generated the same alpha-gamma CFC as in (1) but this time the α waveform was generated with non-sinusoidal wave morphology properties. Concretely, the magnitude of the amplitude of the peaks of the alpha oscillations was higher relative to the troughs:

$\begin{array}{ll}\alpha ={((cos\left(2\pi ft\right)+1)/2)}^b& (9)\end{array}$

This resulted in signals with a “Gaussian-like” wave morphology, where the exponent b controlled the degree of amplitude asymmetry of peaks relative to the troughs. The exponent b was chosen to be 3. The signal α was demeaned and rescaled by a factor of 2. Note that the properties of the non-sinusoidal signal morphology described above differ from the amplitude fluctuation asymmetry properties defined by Mazaheri and Jensen (2008). Amplitude fluctuations presented in this manuscript do not refer to the normalized variance of the peaks and troughs. The alpha morphology model defined here based on the assumption that the magnitude of the amplitude of the peaks (troughs) was higher relative to the magnitude the troughs (peaks). Critically, the variance of the peaks/troughs in our model was highly similar (did not change over time) although the asymmetry in the magnitude of the peaks/troughs was very strong (see Discussion Section).

Beta-Low Gamma Amplitude Harmonics (20–40 Hz) Coupled to Alpha Phase

We simulated a “sawtooth” signal of 10 Hz period with and without frequency fluctuations using the following Equation (Maccarone, 2013):

$\begin{array}{ll}\alpha ={\displaystyle {\sum }_{j=1}^J\frac{1}{j^2}}·cos\left(j2\pi ft+\left(j-1\right)·\phi \right)& (10)\end{array}$

where j indexed the j-th harmonic of frequency f, and the phase term (j − 1) · φ determined the phase of the j-th harmonic. We defined the first two harmonics of the alpha band (20 and 30 Hz, J = 3). The sum of the fundamental and the 1st and 2nd harmonics with the phase term φ = π/2 generated a sawtooth wave whose amplitude raised rapidly decreasing with time in a given cycle (“left-side peak”). This waveform was referred as “negative sawtooth” signal. In contrast, we refered to a “positive sawtooth” (“right-side peak”) a waveform with a constant raising period and a rapidly amplitude decrease (φ = 3π/2). We explored different φ-values [2π, π/4, π/2, 3π/4, π, 5π/4, 3π/2, 7π/4] to investigate systematically a variety of periodic wave morphologies. In comparison to Kramer et al. (2008), who showed spurious CFC due to the sharp edges in a sawtooth signal, we used the $\frac{1}{j^2}$ factor, which results in a more curved sawtooth signal to avoid high broadband frequency activity in the power spectrum (Maccarone, 2013). Importantly, it has been not tested yet whether wave morphologies influences CFD (i.e., non-linear coupling by harmonics of non-sinusoidal oscillations) and under which circumstances this measure provides false positives (Jiang et al., 2015).

Effects of Sharp Edge Artifacts on CFD

Adapted from example 2 in Kramer et al. (2008), we explored how sharp transitions introduced into each alpha cycle would affect not only the CFCoh and bicoherence measures, but critically CFD. We created a 10 Hz sinusoid with a constant period and we added sharp edges that linearly increased the amplitude of the oscillation after 1–5 ms (uniformly distributed over the cycles). The sharp edges began at phase values uniformly distributed between 233 and 265°C (cosine phase).

Results

LZP Modulations in Alpha (8–12 Hz) and Beta (15–30 Hz) Bands during Resting State

Spectral analysis was applied to the MEG axial gradiometers to investigate the power modulations produced by LZP administration during the EO and EC conditions. We found significant differences between drug conditions both during EO (cluster-based nonparametric permutation test, p < 0.001; Figure S2A left) and during EC (p < 0.001; Figure S2A right). Under resting state conditions, LZP reduced alpha power in many sensors, especially in occipital areas. Repeated-measures ANOVA revealed that LZP strongly reduced occipital alpha power during EO [F_{(2, 48)} = 14.48, p < 5 × 10⁻⁵; Placebo = 0.5 mg LZP: t₍₂₄₎ = 1.44; p = 0.16; Placebo > 1.5 mg LZP: t₍₂₄₎ = 4.42, p < 0.001; 0.5 mg LZP > 1.5 mg LZP: t₍₂₄₎ = 3.95, p < 0.001; Figure S2B top]. During EC there was also a robust posterior power reduction by LZP [F_{(2, 48)} = 7.66, p < 0.005; Placebo = 0.5 mg LZP: t₍₂₄₎ = 1.24; p = 0.23; Placebo > 1.5 mg LZP: t₍₂₄₎ = 4.24, p < 0.001; 0.5 mg LZP > 1.5 mg LZP: t₍₂₄₎ = 2.45, p < 0.05; Figure S2B bottom]. To identify the brain sources associated the drug-related power modulations, we pooled the source data over the resting state conditions ((EO + EC)/2) and then considered the relative change: (1.5 mg LZP–Placebo)/Placebo. As expected, the maximum of the LZP alpha power decrease was located over the occipito-parietal cortex (Figure S2C), in line with the LZP-induced power modulations during WM (Lozano-Soldevilla et al., 2014).

While LZP decreased occipital alpha power, fronto-central beta power increased as a function of dosage (Figure S2D) both during EO [F_{(2, 48)} = 12.91, p < 5 × 10⁻⁵; Placebo = 0.5 mg LZP: t₍₂₄₎ = −0.88; p = 0.38; Placebo > 1.5 mg LZP: t₍₂₄₎ = −4.12, p < 0.001; 0.5 mg LZP > 1.5 mg LZP: t₍₂₄₎ = −4.23, p < 0.001; Figure S2E top] and during EC [F_{(2, 48)} = 38.87, p < 1 × 10⁻¹⁰; Placebo > 0.5 mg LZP: t₍₂₄₎ = −3.98; p < 0.001; Placebo > 1.5 mg LZP: t₍₂₄₎ = −7.06, p < 1 × 10⁻⁶; 0.5 mg LZP > 1.5 mg LZP: t₍₂₄₎ = −5.81, p < 1 × 10⁻⁵; Figure S2E bottom]. In line with previous evidence (Jensen et al., 2005; Hall et al., 2010), the LZP-induced beta power-increased was found over sensorimotor cortices (Figure S2F).

In conclusion, LZP produced robust power modulations during resting state conditions: while parieto-occipital alpha power decreased with LZP dosage, sensorimotor beta power increased.

LZP Decreased Alpha-Phase to Beta-Low Gamma Amplitude Coupling during Rest

Given we found a clear alpha power peak in the resting state spectrum, we set out to investigate whether the ongoing alpha phase linked with higher frequency amplitude modulations. Moreover, we also sought whether LZP modulated its amplitude or frequency. First, we calculated the CFCoh metric for all participants and drug sessions during rest. We observed alpha phase (8–12 Hz) to beta-low gamma (20–45 Hz) amplitude coupling over posterior sensors (Figure 1A). Moreover, LZP significantly changed the alpha/beta-low gamma CFC during EC conditions (cluster-based nonparametric permutation test, p < 0.001; Figure 1A right, marked sensors) but not during EO. Repeated-measures ANOVA on sensors of interest revealed that GABAergic enhancement decreased alpha/beta-low gamma CFC under high LZP dosage relative to Placebo and low dosage, with no significant differences between Placebo and low dosage [F_{(2, 48)} = 8.29, p < 0.001; Placebo = 0.5 mg LZP: t₍₂₄₎ = 1.42; p = 0.17; Placebo > 1.5 mg LZP: t₍₂₄₎ = 2.81, p < 0.05; 0.5 mg LZP > 1.5 mg LZP: t₍₂₄₎ = 3.85, p < 0.005; Figures 1B,C]. It is interesting to observe that the beta-low gamma amplitude (Figure 1B) was not detected in the classical spectral analysis (Figure S2B).

FIGURE 1

Figure 1. High LZP dosage decreases alpha phase (8–12 Hz) to beta-low gamma (20–45 Hz) amplitude coupling during EC. (A) Grand mean topographic representation of alpha phase to beta-low gamma amplitude coupling during EO (left) and EC (right) resting state. The color reflect the magnitude of coupling (magnitude squared coherence). The rows reflect the drug levels: Placebo, 0.5 and 1.5 mg LZP. The sensors marked in white displayed significant drug main effect (repeated-measures ANOVA; multiple comparisons controlled by a cluster-based nonparametric permutation test). (B) Grand mean frequency-by-frequency coherence comodulogram of the occipital axial gradiometers that were significantly modulated by LZP [white marked sensors in (A)]. Same color code and scale as in (A). Black rectangles indicate frequency ranges used for the statistical analysis and topographic representations (8–12 Hz; 20–45 Hz). (C) 1.5 mg LZP decreased alpha phase to beta-low gamma amplitude coupling as revealed by repeated-measures ANOVA. **p < 0.005; *p < 0.05. Error bars show the SEM.

It could be that the CFCoh we observe in Figure 1 were produced by harmonics of the alpha oscillations. To test this possibility we computed the magnitude squared bicoherence index for the MEG axial gradiometers. We observed strong bicoherence values within the alpha range during EO peaking over left central and parieto-occipital sensors. The former (left central sensors) were modulated by LZP (cluster-based nonparametric permutation test, p < 0.005; Figures 2A–C left, marked sensors) showing an inverted U-shape [F_{(1.6, 38.9)} = 17.14, p < 5 × 10⁻⁵; Placebo < 0.5 mg LZP: t₍₂₄₎ = −3.02; p < 0.01; Placebo > 1.5 mg LZP: t₍₂₄₎ = 3.73, p < 0.005; 0.5 mg LZP > 1.5 mg LZP: t₍₂₄₎ = 4.88, p < 5 × 10⁻⁴; Figure 2D]. Posterior sensors did not show statistical significant differences as a function of drug dosage. During EC conditions, non-significant differences were found (cluster-based nonparametric permutation test, p = 0.22; Figures 2E–H) but strong bicoherence values were observed over posterior axial gradiometers around ~10 and ~20 Hz. On that note, the CFCoh topographies (Figure 1A) in both resting state conditions remarkably overlapped to both bicoherence (Figure 2) and alpha power topographies (Figure S1A). It is worth to mention that the bicoherence values in the alpha range (8–12 Hz x-axes; 8–12 Hz y-axes; Figure 2, black rectangles) indicated that, over occipital sensors, there were not only consistent phase differences within the ~10 Hz frequency but also between the ~10 and the ~20 Hz frequencies (i.e., ${\left|\sum _{s=1}^S{X}^s(v)X^s(f)X^{*s}(v+f)\right|}^2$), suggesting three harmonics in the power spectrum.

FIGURE 2

Figure 2. Bicoherence modulations by LZP during EO and EC. (A) Grand mean topographic representation of alpha (8–12 Hz) to alpha coupling under Placebo during EO (left) and comodulogram (right). Color code represents magnitude squared bicoherence. Black rectangles in the comodulogram indicate frequency ranges used for the statistical analysis and topographic representations (8–12 Hz; 8–12 Hz). White marked sensors over the left sensorimotor cortex indicate a statistically significant drug main effect (cluster-based nonparametric permutation test, p < 0.005). (B) Topographic representation of bicoherence (left) and comodulogram (right) under 0.5 mg LZP during EO [same conventions as (A)]. (C) Topography representation of bicoherence (left) and comodulogram (right) under 1.5 mg LZP during EO [same conventions as (A)]. (D) LZP modulated sensorimotor (white bold marked sensors) alpha bicoherence in an inverted U-shape as revealed by repeated-measures ANOVA. Error bars show the SEM. (E) Grand mean topographic representation of alpha to alpha coupling under Placebo during EC (left) and comodulogram [right; same conventions as (A)]. The comodulogram displays the average over the occipital white marked sensors from (A)]. (F) Topography and bicoherence comodulogram under 0.5 mg LZP during EC [same conventions as (B)]. (G) Topography and bicoherence comodulogram for 1.5 mg LZP during EC [same conventions as (C)]. (H) LZP did not modulated occipital alpha bicoherence coupling [white bold marked sensors topography in (A)]. ***p < 5 × 10⁻⁴; **p < 0.005; *p < 0.01. Error bars show the SEM.

In summary, parieto-occipital sensors display the strongest phase-to-phase coupling between alpha and its first two harmonics (Figures 2E–G right side).

Figure 1 indicated that there was a CFCoh between alpha and beta-low gamma. Next we asked whether the phase of alpha oscillations drove the power of beta-low gamma oscillations or conversely whether the power of beta-low gamma drove the phase of the alpha oscillations. In addition, a third possibility could appear on where the phase of alpha and the power of beta-low power would interact around zero phase-lag. To achieve that, we used the new CFD metric on MEG axial gradiometers for all participants and sessions during EC. Figure S3A showed the topographic representation of the CFD index (Jiang et al., 2015) over the same frequency ranges as in Figure 1A. We did not observe a clear directionality pattern over posterior sensors. Also, we did not observe significant drug modulations taken our frequency regions of interest (8–12 Hz; 20–45 Hz; cluster-based nonparametric permutation test, p = 1). The CFCoh observed in Figure 1B was not present in the CFD comodulogram in any of the MEG sessions, with very low sigma values ± 0.3 (Figures S3B–D). CFD index did not show the CFCoh coupling pattern found in Figure 2.

The presence of harmonics in the data (Figure 2) could produce spurious CFCoh, raising the question whether the beta-low gamma power band was an independent oscillator or whether it came from the non-sinusoidal properties of the alpha band. To test between these possibilities, complementary CFC metrics were used. To assess the robustness of the beta-low gamma amplitude envelope to alpha phase coupling, we calculated the high frequency TFRs locked to the alpha peaks. If beta-low gamma power is phase locked to the alpha rhythm as the CFCoh index suggested, specific alpha phase segments (i.e., peaks or thoughts) should be associated with power modulations. Furthermore, those modulations should be differently affected by LZP as shown in Figure S2. Contrary to our hypothesis, we did not observe the beta-low gamma power modulations at any alpha cycle and it was not modulated by the drug levels (Figure S4). Despite the clear alpha waves whose central peak (time = 0s) were decreased with LZP (Figures S4D,E,I,J), there were not clear oscillatory beta-low gamma power modulations locked to neither alpha peaks nor troughs. Complementarily, we employed the event-related peak triggered averages around the beta-low gamma frequency range to align signal traces nested to slow oscillations. It is important to note that in case of robust CFC, the high frequency amplitude peak traces (time = 0 s) should be on the top of the specific phase of the ongoing oscillation. Despite the obvious beta-low gamma peaks (Figures S5A,C) no signs of slow oscillations appeared in the average (Figures S5B,D). Given that we have bandpass filtered between the 20 and 45 Hz range, the central averaged peak results from the signal peak alignment whose amplitude is accumulated. The lack of underlying alpha rhythm in the average is less obvious despite the clear spectral peak in the alpha range (Figure S2A), the alpha phase to beta-low gamma CFC (Figure 1B), and the robust alpha peak-triggered average (Figures S4D,I).

Summarizing, event-related oscillatory averages did not allow us to demonstrate clear CFC. This result suggested that the coherence between the alpha band and the beta-low gamma amplitude envelope might be a spurious correlation due to the non-sinusoidal properties of the alpha rhythm.

Occipital Alpha Phase to Beta-Low Gamma Power Coupling during Working Memory Delay Was Not Modulated by LZP

So far, we have focused the analysis on the resting state conditions but it remains to be tested whether similar CFC and drug modulations can be observed under complex cognitive operations. To explore whether the amplitude of high frequency oscillations was phase-locked to the phase of slow rhythms during the WM delay periods, we computed the same CFC measures when participants engaged in a visuo-spatial WM task (Figure S1A). Given previous evidence (Lozano-Soldevilla et al., 2014), the analysis was focused on the WM delay interval (0.5–1.5 s) using the same parameters as done in the resting state analysis.

The strongest coherence appeared in the occipital beta range (16–22 Hz) with very weak modulation in the beta-low gamma band (25–45 Hz; Figures S1B–D). LZP did not modulate the coupling significantly in any frequency range (cluster-based nonparametric permutation test, p > 0.16). Magnitude squared bicoherence showed the strongest coupling in the first two harmonics over occipital sensors (8–12 Hz; 18–22 Hz) but they were not significantly modulated by LZP (cluster-based nonparametric permutation test, p > 0.11; Figure S6). CFD did not show any significant directionality coupling (Figure S7) and LZP did not modulate it.

In conclusion, we observed CFCoh between alpha phase and beta amplitude during the WM delay maintenance period. We found occipital phase-to-phase coherence between alpha and its first two harmonics. CFD did not show convincing evidence of directionality within the frequency range under study. Lastly, the pharmacological GABAergic enhancement did not significantly change any of the CFC metrics used here.

Signal Morphology Produced Spurious CFCoh and CFD

A sizable part of our empirical CFC patters might be explained by harmonics from the alpha band due to its well-known non-sinusoidal wave morphology properties. To test this possibility, we created synthetic data with a variety of wave morphologies to study systematically its influence on CFC metrics. Figure 3A showed CFC between narrow gamma band amplitude (60–80 Hz) and alpha phase (8–12 Hz) as simulated in case 1. In addition, we observed a phase-phase coupling using bicoherence, with no contributions in at the harmonics of the alpha wave. In addition, we observed a robust CFD, consistent with the simulated signal where the gamma bursts were leading the alpha cycles (between 0 and 20°C of the phase of the slow frequency). The observed broadband coherence and directionality extended well beyond the simulated parameters (60–80 Hz), presumably due to spectral leakage. Figure 3B represented the simulated signal from case 2 with a non-sinusoidal morphology of the alpha signal: higher magnitude of the amplitude of the peaks compared to the troughs. The power spectrum revealed three prominent peaks around ~10 Hz (fundamental frequency), ~20 and ~30 Hz (harmonic frequencies). As a consequence of the harmonics, CFCoh and bicoherence not only showed the high gamma to alpha (fundamental component ~10 Hz; x-axes CFCoh comodulogram Figure 3B) coupling but also CFC at the harmonic frequencies. Interestingly we found CFD not only to be robust between high gamma and the ~10 Hz phase but also with the ~20 Hz phase, similar as in the CFCoh. This can be explained because the harmonics were, by definition, phase-coupled to the fundamental frequency and therefore the high gamma coupled to all three components. Strikingly, no CFD was found at the lower frequencies (10–30 Hz), whereas the latter showed strong CFCoh.

FIGURE 3

Figure 3. (A) Simulated alpha sinusoid (~10 Hz) and gamma band (60–80 Hz) CFC. Top: in blue full cosine cycle of a ~10 Hz signal coupled to a gamma 50 ms burst (red). Middle left: time course of the alpha-gamma coupled signal with random noise and the power spectrum in decibels (dB). Middle right: bicoherence of the simulated coupled signal. Color code represents the magnitude squared coherence between a pair of frequencies (x and y axes) and its sum. Bottom left: CFCoh of the simulated signal. Color code represents magnitude squared coherence between y-axes amplitude envelope frequencies and the phase of the raw signal (x-axes). Bottom right: CFD between high frequencies (y-axes) and low frequencies (x-axes). Color code indicates the normalized $\tilde{\Psi }$ weighted average of the CFD. (B) Simulated alpha signal with non-sinusoidal shape (~10 Hz) coupled to the amplitude of the gamma band (60–80 Hz). Top: in blue one asymmetric cycle on a ~10 Hz signal coupled to a gamma 50 ms burst (red). Same figure conventions as in (A). a.u., Arbitrary units.

To further assess the sensitivity of CFC metrics to signal morphology, in simulation case 3 we constructed a sawtooth signal with a fundamental and two harmonics, varying the signal morphology and keeping the frequency constant (i.e., no period fluctuations). Figure 4A showed a curved negative sawtooth with very clean spectral peaks at the alpha frequency and its first two harmonics. We found that CFCoh and bicoherence captured the spurious coupling between frequencies, whereas CFD showed low values not exceeding the significance threshold. Non-surprisingly, adding frequency fluctuations to the same sawtooth simulation yielded significant CFD (Figure 4B).

FIGURE 4

Figure 4. (A) Signal wave morphology and period fluctuations produced spurious CFC and CFD. Top: full cosine cycle of the 10 Hz (blue), 20 Hz (green), and 30 Hz (red), and the resulted signal (black) taking the sum over all the signals (F1 + H1 + H2 = total). Harmonic phase φ was set to π/2 (see simulation case 3 equation for more details). Middle left: time course of the alpha plus harmonics with random noise and the power spectrum in decibels (dB). Middle right: bicoherence of the simulated sawtooth signal. Color code represents the magnitude squared coherence between a pair of frequencies (x- and y-axes) and its sum. Bottom left: CFCoh of the simulated signal. Color code represents magnitude squared coherence between y-axes amplitude envelope frequencies and the phase of the raw signal (x-axes). Bottom right: CFD between high frequencies (y-axes) and low frequencies (x-axes). Color code indicates the normalized $\tilde{\Psi }$ weighted average of the CFD slope. (B) Time course of the simulated alpha sinusoid (F1 = ~10 Hz) and its first two harmonics (H1 = ~20 Hz; H2 = ~30 Hz) with period fluctuations and random noise. Same sawtooth wave as in Figure 2A top was used. Below, the power spectrum in decibels (dB). CFCoh, CFD and bicoherence have the same figure conventions as in Figure 2A. (C) Simulated alpha sinusoid (F1 = ~10 Hz) and its first two harmonics with period fluctuations. Top: full cosine cycle of the ~10 Hz (blue), ~20 Hz (green), and ~30 Hz (red) and the resulted signal (black) taking the sum over all the signals (F1 + H1 + H2 = total). Harmonic phase φ was set to 2π (see simulation case 3) equation for more details). Same figure conventions as in (A). (D) Simulated imperfect sinusoid (see simulation case 4). Top: full cosine cycle of 10 Hz with sharp edges and no period fluctuations. Same figure conventions as in (A).

Significant CFC could be observed in the absence of CFD if the amplitude envelope of the higher frequencies and the analytic phase of the slow frequencies covary in time without phase delays (Jiang et al., 2015). Some type of periodic non-sinusoidal wave morphologies might be unable to detect significant CFD due to the consistent zero phase-lag CFC between the fundamental frequency and its harmonics. To test this possibility, we made a variant of case 3 in which the harmonics were aligned to the fundamental frequency, showing non-sinusoidal shape properties without phase lead, nor lag (Figure 4C). As predicted, we observed spurious CFC between alpha and its harmonics but no significant directionality within the harmonic frequency ranges. To explore whether CFD could be insensitive to spurious CFC due to the non-sinusoidal properties of an alpha-like signal, we systematically modulated the φ-value to simulate a variety of waveforms that are displayed in Figure 5 (see simulation case 3 for details). We found two cases on where CFD between alpha-like non-sinusoidal synthetic signal and its harmonics: the phase φ being π or 2π. In both cases the signal had a non-sinusoidal wave morphology (π, troughs with more amplitude than peaks; 2π, peaks more amplitude than troughs) and near zero phase-lag between the amplitude envelope of the harmonics and the phase of the fundamental frequency. Among other φ values, CFD changed as a function of the sawtooth peak: negative CFD was obtained with a negative sawtooth waveforms and positive CFD was obtained with a positive sawtooth waveforms. Comparing the upper (at 12:00) and lower (at 6:00) panels in Figure 5 one can observe the sign of the directionality being flipped. In conclusion, we identified for the first time the wave morphologies that produced spurious CFCoh and CFD.

FIGURE 5

Figure 5. CFD index changed as a function of wave morphology. Simulated alpha sinusoid (F1 = ~10 Hz) and its first two harmonics (H1 = ~20 Hz; H2 = ~30 Hz) with period fluctuations. Inner circle: full cosine cycles of the ~10 Hz (blue), ~20 Hz (green) and ~30 Hz (red) and the resulted signal (black) taking the sum over all the signals (F1 + H1 + H2 = total). Harmonic phase φ was set to the different phases as indicated in the circle (see simulation case 3 equation for more details) organized in a clock-wise manner. Each time-resolved single cycle was associated to CFCoh, CFD, and bicoherence plots. Figure scales and conventions were the same as in Figure 1.

Sharp Edges in the Signal Produce Spurious CFD

Although CFD should not yield significant directionality when the frequency of a periodic signal does not fluctuate, sharp edge transients could still introduce spurious directionality. To assess the influence of sharp edge artifacts on the CFD index, in simulation case 4 we replicated the example described in Kramer et al. (2008). As seen in Figure 4D, we observed not only strong CFCoh between the fundamental frequency and multiple harmonics, but also robust CFD between multiple frequencies. Therefore, we found that CFD yields multiple false positives when oscillations contain edge artifacts (Jiang et al., 2015). This is true even when the periodicity of the signal did not fluctuate.

Discussion

Here, we analyzed a pharmaco-MEG dataset to investigate how CFC metrics were modulated by LZP during rest and WM maintenance. During rest, we demonstrated that while occipital alpha power decreased with LZP (Link et al., 1991; Schreckenberger et al., 2004; Ahveninen et al., 2007), sensorimotor beta power increased with drug dosage (Baker and Baker, 2003; Jensen et al., 2005; Hall et al., 2010). We found alpha phase to beta-low gamma amplitude coupling over occipital sensors that decreased with 1.5 mg LZP during EC. However, CFD and peak-triggered averages did not provide additional phase-to-amplitude coupling evidence. In addition, we found phase-to-phase coupling between alpha and beta over occipital sensors. The strongest CFCoh we observed during WM conditions was between alpha phase ~10 Hz and beta amplitude ~20 Hz, with robust phase-to-phase alpha-beta coupling without significant drug modulations.

Cross-Frequency Coupling between Alpha and Beta-Low Gamma

Several MEG studies have investigated on the role of phase-to-amplitude coupling under rest and task conditions. In one of the first attempts, Osipova et al. (2008) found that the ongoing gamma (30–70 Hz) power was phase-locked to the alpha band (8–12 Hz) over posterior regions under resting state. This was mainly observed in subjects with strong alpha power. We however found a coupling in lower frequency range over occipital sensors in a bigger sample. A similar gamma frequency range (19–76 Hz) phase-locked to slower alpha (5–13 Hz) was reported by Hawellek et al. (2013). A very recent study found that the phase of spontaneous alpha oscillations in thalamus coupled to ~40 Hz gamma power in parieto-occipital cortex during EC conditions (Roux et al., 2013). ECoC recordings in humans reported the same phenomenon both in lower (Fitzgerald et al., 2013) and higher gamma band frequency during resting state (Bahramisharif et al., 2013). In summary, the absence of high frequency gamma power in our MEG data could be related to due to a lower signal-to-noise ratio for MEG as compared to ECoG data.

Neurophysiological Oscillations with Non-sinusoidal Wave Morphologies Produce Spurious CFC

The results of the simulations presented here raised a challenge to CFC analyses of signals whose waveform is known to be non-sinusoidal, such as the occipital alpha (Stam et al., 1999; Nikulin et al., 2007), hippocampal theta (Belluscio et al., 2012), and the sensorimotor mu-rhythm (Salmelin and Hari, 1994). The latter is in fact quite revealing, given that it was originally called “rhythme en arceau” by Gastaut (1952) due to its arch-like shape.

Our synthetic data showed convincingly that periodic signals deviating from sines and cosines (i.e., the filter kernel) produced harmonics in the power spectrum and these harmonics were, by definition, phase-to-amplitude coupled (Figure 1). Taken this fact into account, we cannot claim with absolute certainty that all CFCoh patterns we observed in our empirical data (Figure 1) were produced by the non-sinusoidal properties of the wave morphology of the alpha rhythm. However, it is undisputable that a sizable part of our CFCoh was explained by harmonics generated by the non-sinusoidal wave morphology properties of the alpha rhythm.

The non-sinusoidal shape properties of the alpha band (i.e., magnitude of the peaks higher relative to the magnitude of the troughs) should not be confused with the alpha asymmetric amplitude fluctuation properties defined by Mazaheri and Jensen (2008). The asymmetric amplitude fluctuations were defined as the normalized ratio of the variance (not the amplitude) of the peaks of an oscillatory activity minus the variance of the troughs (Mazaheri and Jensen, 2008). The authors used a synthetic signal with the form:

$\begin{array}{ll}a=A\left(t\right)*(1+sin\left(2\pi ft\right))& (11)\end{array}$

where A(t) is a slow signal (i.e., 3 Hz) that determined the amplitude modulation of the alpha frequency f. Note that, by construction, this synthetic signal will not produce harmonics because it consisted of a DC offset term and a single sine wave. This type of model might not apply to the empirical data presented here because we observed strong coupling between 10 Hz and its first two harmonics using bicoherence (Figure 2; Figure S6).

Summarizing, there are physiological reasons to believe that some neurophysiological signals have non-sinusoidal properties (i.e., waveforms deviating from sines and cosines). For instance, when assessing the coupling between different frequencies, the potential harmonic contribution to CFC should be taken into account using bicoherence or any other metric sensitive to harmonic.

Non-zero CFC was Produced by a Mismatch between Filter Kernel and Signal Morphology

The results presented here were obtained using only a limited number of CFC metrics that were based on Fourier methods. We expect other methods based on band-pass filtering and Hilbert transform to yield, in principle, similar results (Canolty et al., 2006; Cohen, 2008; Penny et al., 2008; Tort et al., 2010). The reason is because there is an equivalence between wavelet, Fourier and Hilbert methods on amplitude spectral estimation (Bruns, 2004). More importantly, Lepage et al. (2013) demonstrated that the phase variance and bias of the three estimators “are shown, in simulation and in theory, to depend on the extent to which signal present in the data is matched to the phase estimator” (p. 915). While wavelet phase estimator outperformed Fourier and Hilbert transforms when the signal to be analyzed had a broadband pulse-like wave morphology, the Fourier approach yielded more accurate phase estimates in comparison with the other two methods when the signal to be analyzed showed a stationary sinusoidal-like wave morphology (see Lepage et al., 2013 for details). Spectral analysis of periodic signals whose wave morphologies systematically mismatch a given filter kernel (i.e., wavelet, sine/cosine) will produce harmonics, and in consequence, spurious CFC (Lozano-Soldevilla, 2015). This observation has been recently replicated by various independent research groups (Cole et al., 2016; Gerber et al., 2016; Jones, 2016; Scheffer-Teixeira and Tort, 2016).

Following this rationale, if the filter kernel can be adapted to the signal under analysis, the spectral estimation will be more accurate. Future CFC metrics could use filter kernels adapted to non-sinusoidal signals to avoid spurious correlations between frequency components.

Cross-Frequency Directionality Index Produced Spurious Directionality

Present simulations and empirical results showed that spurious CFCoh and CFD may be indeed found between the phase of the fundamental frequency of slow oscillations and the amplitude envelope of its harmonics. Periodic signals with sawtooth-like wave morphologies produce spurious CFD. The CFD proposed by Jiang et al. (2015) is the PSI applied to the cross-frequency domain. PSI has been shown to be resistant to false positives when even when the time series to be assessed contained non-linear interactions (Nolte et al., 2010). The PSI computes directionality by comparing two independent time series, but the CFD was computed here in the same recording (one time series) between different frequency components. The simulations presented here clearly demonstrated that spurious coupling (CFCoh and CFD) between the fundamental frequency and its harmonics can also be narrowband (20–45 Hz). Our synthetic data revealed for the first time that the CFD metric gave false positives not only when the signal contain sharp edge artifacts (Kramer et al., 2008) but also with signal whose wave morphology was sawtooth-like. We want to stress that this is especially problematic if the CFC occurs over the first harmonics as in our case. We recommend future CFD users to check the low frequency signal morphology and other potential sources of non-linear coupling that could produce spurious results (Aru et al., 2015).

One way to reduce this problem would consist in computing CFCoh and CFD between sets of independent recordings (Maris et al., 2011; van der Meij et al., 2012, 2015; Aru et al., 2015). This approach presumably attenuatates contribution of volume conduction that could potentially spread non-sinusoidal periodic signals that could yield spurious correlations. For example, Spaak et al. found stronger CFC between layers in monkey V1 (alpha phase of infragranular layer coupled to gamma power of supragranular layer) in comparison with CFC within the infragranular layer (Spaak et al., 2012). However, a recent study found that the infragranular alpha sources might reflected a mix of local activity of field potentials (LFP) and contaminated activity that was volume-conduced from supragranular layers (Haegens et al., 2015). This adds more complexity not only to the interpretation of LFP and current source density results but also to the usefulness of the CFC in cases on where the cortical connectome unambiguously determined.

Datasets including spikes and LFP could help to dissociate the gamma oscillations from the potential harmonics from non-sinusoidal properties of slower rhythms. In this direction, Colgin et al. (2009) found that both fast and slow gamma oscillations in different circuitry of rat hippocampus locked to spikes at different phases. Specifically, this would mean that both slow and fast gamma are rhythms not only related to spiking activity between regions but also modulated by the phase of the ongoing theta activity. However, this approach is not exempt of problems. It is not easy to establish unambiguously which variables are synchronized using different coherence measures (see Buzsáki and Schomburg, 2015 for an elaborated argumentation). Gamma power can be contaminated by the spike activity, making very challenging to disentangle independent oscillatory gamma estimates from the spike activity (Buzsáki and Schomburg, 2015; Ray and Maunsell, 2015). This is a similar problem described here, but this time the harmonic frequencies are originated from the non-sinusoidal properties of the spike wave morphology. Indeed, an empirical study found that the contribution of these spike-related harmonics can go as low as 10 Hz (Waldert et al., 2013). Although spiking activity can definitely help to understand LFP activity, its interpretation in relation to neuronal oscillations needs to be done carefully.

Conclusion

In conclusion, our findings showed CFC between spontaneous alpha and beta-low gamma over occipital sensors during rest and WM. However, the observed phase-to-amplitude coupling contained harmonic power contributions presumably from the non-sinusoidal properties of the alpha wave morphology. Although alpha harmonics may not fully explain the whole CFC pattern observed here, this seems to be the most parsimonious interpretation. In a broader scope, the present results highlighted the fact that periodic signals with non-sinusoidal wave morphology generates non-zero CFC and spurious CFD (i.e., sawtooth-like morphologies). We suggest a cautionary interpretation of CFC results based on current models that assume that communication between neuronal ensembles are implemented as cross-frequency interactions between gamma and alpha activity (Jensen and Mazaheri, 2010).

Author Contributions

DL designed the study, collected and analyzed the data. DL, Nt, and RO discussed the results and wrote the paper.

Conflict of Interest Statement

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.

Acknowledgments

Results presented in this manuscript have been previously published in chapters 4 and 5 in a dissertation (Lozano-Soldevilla, 2015). This work was supported by the Netherlands Organization for Scientific Research by an Open Competition/MaGW grant n.r. 400-09-491 to Roshan Cools and Ole Jensen. We thank Haiteng Jiang for sharing the code to compute the cross-frequency directionality metric.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fncom.2016.00087

Figure S1. LZP did not modulate alpha-phase to beta-low gamma amplitude coupling during WM delay. (A) Delayed match to sample WM task. Participants had to keep fixations while covertly attend to the cued visual hemifield to encode the target items. The task consisted to compare the cued items, compare them with a probe and decide whether they matched or not. (B) Grand mean topographic representation of alpha-phase to beta (16–22 Hz; left) and to beta-low gamma amplitude coupling (25–45 Hz; right) during WM delay (0.5 to 1.5 s). Color code represents magnitude squared coherence. Note that topographies have different scale. The strongest coherence values appeared over occipito-parietal sensors. Grand mean frequency-by-frequency coherence comodulogram (middle) of the occipital axial gradiometers of interest as in Figure S2A. Black rectangles indicate the x-y frequency selections used in the statistical analysis and topographic representations (8–12 Hz; 16–22 Hz; 25–45 Hz). Comodulograms have the same coherence scale as alpha—beta topography (left). (C,D) are the same as (B) but for 0.5 and 1.5 mg LZP, respectively.

Figure S2. LZP modulations in the alpha (8–12 Hz) and beta (15–30 Hz) bands during rest. (A) Grand mean topographic representation of alpha power during EO (left) and EC (right). Color code represents relative power changes with respect to the mean of the power spectrum (5–40 Hz). The rows reflect the drug levels: Placebo, 0.5 and 1.5 mg LZP. Sensors marked in bold display statistical significant drug main effect (repeated-measures ANOVA) using a cluster-based nonparametric permutation test. (B) Grand average of power spectra [sensors marked white in (A)] estimated for the EO (top) and EC (bottom) and all drug sessions. Inserts show a repeated-measures ANOVA confirming that the occipital alpha power decreased with LZP. Error bars show the SEM. (C) Occipito-parietal sources reflected the alpha power decrease under 1.5 mg LZP dosage relative to placebo (EO and EC averaged). (D) Gran mean topographic representation of beta power during EO (left) and EC (right) during rest [same conventions as in (A)]. (E) Grand average of power spectra [sensors marked white (D)] estimated for the EO (top) and EC (bottom) and all drug sessions. Inserts showed a repeated-measures ANOVA confirming that the sensorimotor beta power increased with LZP dosage; same conventions as in Figure 1B. Error bars show the SEM. (F) The sources in the sensorimotor cortex reflected the beta power increase under 1.5 mg LZP dosage relative to placebo (EO and EC averaged) [same conventions as in (C)]. The peak of the beta source lied out of the brain due to inaccuracies after the individual head model normalization to the MNI space. ^†p < 1 × 10⁻⁶; ^***p < 1 × 10⁻⁵; ^**p < 0.001; ^*p < 0.05.

Figure S3. LZP did not significantly change CFD during EC. (A) Grand mean topographic representation of the CFD (8–12 Hz; 20–45 Hz) for the different drug conditions. White sensors were significantly modulated by LZP in Figure 1. The color code indicates the normalized $\tilde{\Psi }$ weighted average of the CFD slope (Ψ = $\tilde{\Psi }$/std($\tilde{\Psi }$)). (B) Grand mean normalized CFD index averaged over the occipital sensors marked in topographies (A) for each drug session. Note the low sigma values we obtained (± ~0.3; ± 2 was considered significant). (C) CFD weighted average ($\tilde{\Psi }$) and (D) standard deviation of the slope (std($\tilde{\Psi }$)) averaged over the occipital sensors marked in topographies for each drug session.

Figure S4. Time-frequency representations (TFR) of induced power locked to the peaks of the alpha oscillations during rest. (A) Grand mean TFRs from left and right hemisphere MEG axial gradiometers [for white sensors in (E)] time-locked to the alpha peaks (t = 0 s) under placebo session during EO. TFRs were calculated for each single trial and then averaged. Color code represents the relative power normalized to the average power per frequency (y-axes). (B) TFRs of induced power locked to the peak of alpha oscillations for 0.5 mg LZP during EO. Same conventions as in (A). (C) TFRs of induced power locked to the peak of alpha oscillations during 1.5 mg LZP during EO. Same conventions as in (A). (D) Grand mean of event-related fields locked to the alpha peaks. The peaks were found after bandpass filtering of the data (8–12 Hz). MEG axial gradiometers from left and right hemispheres (marked white in topography) during EO for all drug levels. (E) Grand mean topographic representation of alpha peak-triggered average over the period (–0.015 to 0.015 s) around the alpha peak (t = 0 s). Marked sensors indicate statistical significant main effect (repeated-measures ANOVA) between treatments using the cluster-based nonparametric permutation test. (F–J) are the same as (A–E) but for EC respectively.

Figure S5. Beta-low gamma peak-triggered average during rest. (A) Grand mean of event-related fields locked to the beta-low gamma (20–45 Hz) peaks. MEG axial gradiometers from left and right hemispheres [marked white in (B)] during EO for all drug sessions. Vertical dashed lines at ± 0.05 s indicates the lower limit of the bandpass filter used for peak detection (1/0.05 s = 20 Hz). Beyond the ± 0.05 s limit, the underlying alpha oscillations were highly reduced if any. (B) Grand mean topographic representation of beta-low gamma peak-triggered average taken around the peak period (–0.015 to 0.015 s). Marked bold sensors of interest based on Figure S2A. Panels (C,D) are the same as (A,B) but for EC, respectively with same figure conventions.

Figure S6. LZP did not modulate bicoherence values during WM delay. (A) Grand mean topographic representation of alpha–alpha (8–12 Hz; 8–12 Hz; left) and alpha–beta coupling (8–12 Hz; 18–22 Hz; right) during WM delay (0.5–1.5 s). Color code represents magnitude-squared bicoherence. Note that the topographic plots have different scale. The strongest bicoherence values appeared over occipito-parietal sensors. Grand mean frequency-by-frequency bicoherence comodulogram (middle) of occipital axial gradiometers of interest chosen from Figure S2A. Black rectangles indicate the x-y frequency selections used in the statistical analysis and topographic representations (8–12 Hz; 8–12 Hz; 18–22 Hz). Comodulograms have the same scale as alpha–alpha topography (left). Panels (B,C) are the same as (A) but for 0.5 and 1.5 mg LZP, respectively.

Figure S7. LZP did not change CFD during WM delay. (A) Grand mean topographic representation of alpha-beta CFD (16–22 Hz) during WM delay for each drug separately. Color code indicates the normalized $\tilde{\Psi }$ weighted average of the CFD slope [Ψ = $\tilde{\Psi }$/std($\tilde{\Psi }$); color scale same as in (B)]. (B) Grand mean normalized CFD index taken over the occipital sensors of interest chosen from Figure S2A. Note the low sigma values we obtained (~ ±0.6; ±2 is considered significant). (C) CFD weighted average ($\tilde{\Psi }$) and (D) standard deviation of the slope (std($\tilde{\Psi }$)) taken over the same occipital sensors for each drug session respectively. Drug sessions are displayed in different columns.

References

Ahveninen, J., Lin, F. H., Kivisaari, R., Autti, T., Hämäläinen, M., Stufflebeam, S., et al. (2007). MRI-constrained spectral imaging of benzodiazepine modulation of spontaneous neuromagnetic activity in human cortex. Neuroimage 35, 577–582. doi: 10.1016/j.neuroimage.2006.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Aru, J., Aru, J., Priesemann, V., Wibral, M., Lana, L., Pipa, G., et al. (2015). Untangling cross-frequency coupling in neuroscience. Curr. Opin. Neurobiol. 31, 51–61. doi: 10.1016/j.conb.2014.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Axmacher, N., Henseler, M. M., Jensen, O., Weinreich, I., Elger, C. E., and Fell, J. (2010). Cross-frequency coupling supports multi-item working memory in the human hippocampus. Proc. Natl. Acad. Sci. U.S.A. 107, 3228–3233. doi: 10.1073/pnas.0911531107

PubMed Abstract | CrossRef Full Text | Google Scholar

Bahramisharif, A., van Gerven, M. A. J., Aarnoutse, E. J., Mercier, M. R., Schwartz, T. H., Foxe, J. J., et al. (2013). Propagating neocortical gamma bursts are coordinated by traveling alpha waves. J. Neurosci. 33, 18849–18854. doi: 10.1523/jneurosci.2455-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Baker, M. R., and Baker, S. N. (2003). The effect of diazepam on motor cortical oscillations and corticomuscular coherence studied in man. J. Physiol. 546(Pt 3), 931–942. doi: 10.1113/jphysiol.2002.029553

PubMed Abstract | CrossRef Full Text | Google Scholar

Belluscio, M. A., Mizuseki, K., Schmidt, R., Kempter, R., and Buzsáki, G. (2012). Cross-frequency phase–phase coupling between theta and gamma oscillations in the hippocampus. J. Neurosci. 32, 423–435. doi: 10.1523/jneurosci.4122-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B (Stat. Method) 57, 289–300. doi: 10.2307/2346101

CrossRef Full Text | Google Scholar

Bragin, A., Jandó, G., Nádasdy, Z., Hetke, J., Wise, K., and Buzsáki, G. (1995). Gamma (40–100 Hz) oscillation in the hippocampus of the behaving rat. J. Neurosci. 15, 47–60.

PubMed Abstract | Google Scholar

Bruns, A. (2004). Fourier-, Hilbert- and wavelet-based signal analysis: are they really different approaches? J. Neurosci. Methods 137, 321–332. doi: 10.1016/j.jneumeth.2004.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of the Brain. Oxford, UK: Oxford University Press.

Google Scholar

Buzsáki, G., Logothetis, N., and Singer, W. (2013). Scaling brain size, keeping timing: evolutionary preservation of brain rhythms. Neuron 80, 751–764. doi: 10.1016/j.neuron.2013.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Schomburg, E. W. (2015). What does gamma coherence tell us about inter-regional neural communication? Nat. Neurosci. 18, 484–489. doi: 10.1038/nn.3952

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., Ganguly, K., and Carmena, J. M. (2012). Task-dependent changes in cross-level coupling between single neurons and oscillatory activity in multiscale networks. PLoS Comput. Biol. 8:e1002809. doi: 10.1371/journal.pcbi.1002809

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., Ganguly, K., Kennerley, S. W., Cadieu, C. F., Koepsell, K., Wallis, J. D., et al. (2010). Oscillatory phase coupling coordinates anatomically dispersed functional cell assemblies. Proc. Natl. Acad. Sci. U.S.A. 107, 17356–17361. doi: 10.1073/pnas.1008306107

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Cole, S. R., Peterson, E. J., van der Meij, R., de Hemptinne, C., Starr, P. A., and Voytek, B. (2016). Nonsinusoidal oscillations underlie pathological phase-amplitude coupling in the motor cortex in Parkinson's disease. bioRxiv. doi: 10.1101/049304

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Contreras, D., Destexhe, A., Sejnowski, T. J., and Steriade, M. (1997). Spatiotemporal patterns of spindle oscillations in cortex and thalamus. J. Neurosci. 17, 1179–1196.

PubMed Abstract | Google Scholar

Fitzgerald, T. H. B., Valentin, A., Selway, R., and Richardson, M. P. (2013). Cross-frequency coupling within and between the human thalamus and neocortex. Front. Hum. Neurosci. 7:84. doi: 10.3389/fnhum.2013.00084

PubMed Abstract | CrossRef Full Text | Google Scholar

Gastaut, H. (1952). Etude électrocorticographique de la réactivité des rythmes rolandiques. Rev. Neurol. (Paris) 87, 176–182.

PubMed Abstract

Gerber, E. M., Sadeh, B., Ward, A., Knight, R. T., and Deouell, L. Y. (2016). Periodic non-sinusoidal activity can produce cross-frequency coupling in cortical signals in the absence of functional interaction between neural sources. bioRxiv. doi: 10.1101/062190

CrossRef Full Text | Google Scholar

Gotman, J. (1983). Measurement of small time differences between EEG channels: method and application to epileptic seizure propagation. Electroencephalogr. Clin. Neurophysiol. 56, 501–514. doi: 10.1016/0013-4694(83)90235-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Groppe, D. M., Makeig, S., and Kutas, M. (2009). Identifying reliable independent components via split-half comparisons. Neuroimage 45, 1199–1211. doi: 10.1016/j.neuroimage.2008.12.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Gross, J., Kujala, J., Hämäläinen, M., Timmermann, L., Schnitzler, A., and Salmelin, R. (2001). Dynamic imaging of coherent sources: studying neural interactions in the human brain. Proc. Natl. Acad. Sci. U.S.A. 98, 694–699. doi: 10.1073/pnas.98.2.694

PubMed Abstract | CrossRef Full Text | Google Scholar

Haegens, S., Barczak, A., Musacchia, G., Lipton, M. L., Mehta, A. D., Lakatos, P., et al. (2015). Laminar profile and physiology of the α rhythm in primary visual, auditory, and somatosensory regions of neocortex. J. Neurosci. 35, 14341–14352. doi: 10.1523/jneurosci.0600-15.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Hall, S. D., Barnes, G. R., Furlong, P. L., Seri, S., and Hillebrand, A. (2010). Neuronal network pharmacodynamics of GABAergic modulation in the human cortex determined using pharmaco-magnetoencephalography. Hum. Brain Mapp. 31, 581–594. doi: 10.1002/hbm.20889

PubMed Abstract | CrossRef Full Text | Google Scholar

Hawellek, D. J., Schepers, I. M., Roeder, B., Engel, A. K., Siegel, M., and Hipp, J. F. (2013). Altered intrinsic neuronal interactions in the visual cortex of the blind. J. Neurosci. 33, 17072–17080. doi: 10.1523/jneurosci.1625-13.2013

PubMed Abstract | CrossRef Full Text

Igarashi, K. M., Lu, L., Colgin, L. L., Moser, M.-B., and Moser, E. I. (2014). Coordination of entorhinal-hippocampal ensemble activity during associative learning. Nature 510, 143–147. doi: 10.1038/nature13162

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, O., and Mazaheri, A. (2010). Shaping functional architecture by oscillatory alpha activity: gating by inhibition. Front. Hum. Neurosci. 5:12. doi: 10.3389/fnhum.2010.00186

CrossRef Full Text

Jensen, O., Goel, P., Kopell, N., Pohja, M., Hari, R., and Ermentrout, B. (2005). On the human sensorimotor-cortex beta rhythm: sources and modeling. Neuroimage. 26, 347–355. doi: 10.1016/j.neuroimage.2005.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, H., Bahramisharif, A., van Gerven, J. M. A., and Jensen, O. (2015). Measuring directionality between neuronal oscillations of different frequencies. Neuroimage 118, 359–367. doi: 10.1016/j.neuroimage.2015.05.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, S. R. (2016). When brain rhythms aren't ‘rhythmic’: implication for their mechanisms and meaning. Curr. Opin. Neurobiol. 40, 72–80. doi: 10.1016/j.conb.2016.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, Y. C., and Powers, E. J. (1979). Digital bispectral analysis and its applications to nonlinear wave interactions. IEEE Trans. Plasma Sci. 7, 120–131. doi: 10.1109/tps.1979.4317207

CrossRef Full Text | Google Scholar

Kramer, M. A., Tort, A. B. L., and Kopell, N. J. (2008). Sharp edge artifacts and spurious coupling in EEG frequency comodulation measures. J. Neurosci. Methods 170, 352–357. doi: 10.1016/j.jneumeth.2008.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Lepage, K. Q., Kramer, M. A., and Eden, U. T. (2013). Some sampling properties of common phase estimators. Neural Comput. 25, 901–921. doi: 10.1162/NECO_a_00422

PubMed Abstract | CrossRef Full Text | Google Scholar

Link, C. G., Leigh, T. J., and Fell, G. L. (1991). Effects of granisetron and lorazepam, alone and in combination, on the EEG of human volunteers. Br. J. Clin. Pharmacol. 31, 93–97.

PubMed Abstract | Google Scholar

Lozano-Soldevilla, D. (2015). GABAergic Modulations of Gamma and Alpha Oscillations: Consequences for Working Memory Performance. Donders Series 201. Netherlands; Nijmegen: Radboud University Nijmegen. doi: 10.13140/RG.2.1.2714.6322

CrossRef Full Text | Google Scholar

Lozano-Soldevilla, D., ter Huurne, N., Cools, R., and Jensen, O. (2014). GABAergic modulation of visual gamma and alpha oscillations and its consequences for working memory performance. Curr. Biol. 24, 2878–2887. doi: 10.1016/j.cub.2014.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Maccarone, T. J. (2013). The biphase explained: understanding the asymmetries in coupled Fourier components of astronomical time series. MNRAS 435, 3547–3558. doi: 10.1093/mnras/stt1546

CrossRef Full Text | Google Scholar

Maris, E., and Oostenveld, R. (2007). Nonparametric statistical testing of EEG- and MEG-data. J. Neurosci. Methods 164, 177–190. doi: 10.1016/j.jneumeth.2007.03.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Maris, E., van Vugt, M., and Kahana, M. (2011). Spatially distributed patterns of oscillatory coupling between high-frequency amplitudes and low-frequency phases in human iEEG. Neuroimage 54, 836–850. doi: 10.1016/j.neuroimage.2010.09.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazaheri, A., and Jensen, O. (2008). Asymmetric amplitude modulations of brain oscillations generate slow evoked responses. J. Neurosci. 28, 7781–7787. doi: 10.1523/JNEUROSCI.1631-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Mima, T., Steger, J., Schulman, A. E., Gerloff, C., and Hallett, M. (2000). Electroencephalographic measurement of motor cortex control of muscle activity in humans. Clin. Neurophysiol. 111, 326–337. doi: 10.1016/S1388-2457(99)00229-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikulin, V. V., Linkenkaer-Hansen, K., Nolte, G., and Curio, G. (2010). Non-zero mean and asymmetry of neuronal oscillations have different implications for evoked responses. Clin. Neurophysiol. 121, 186–193. doi: 10.1016/j.clinph.2009.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikulin, V. V., Linkenkaer-Hansen, K., Nolte, G., Lemm, S., Müller, K. R., Ilmoniemi, R. J., et al. (2007). A novel mechanism for evoked responses in the human brain. Eur. J. Neurosci. 25, 3146–3154. doi: 10.1111/j.1460-9568.2007.05553.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nolte, G. (2003). The magnetic lead field theorem in the quasi-static approximation and its use for magnetoencephalography forward calculation in realistic volume conductors. Phys. Med. Biol. 48, 3637–3652.

PubMed Abstract | Google Scholar

Nolte, G., Ziehe, A., Krämer, N., Popescu, F., and Müller, K.-R. (2010). Comparison of granger causality and phase slope index. JMLR W & CP 6, 267–276. Available online at: http://www.jmlr.org/proceedings/papers/v6/nolte10a.html

Nolte, G., Ziehe, A., Nikulin, V. V., Schlögl, A., Krämer, N., Brismar, T., et al. (2008). Robustly estimating the flow direction of information in complex physical systems. Phys. Rev. Lett. 100:234101. doi: 10.1103/PhysRevLett.100.234101

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text

Osipova, D., Hermes, D., and Jensen, O. (2008). Gamma power is phase-locked to posterior alpha activity. PLoS ONE 3:e3990. doi: 10.1371/journal.pone.0003990

PubMed Abstract | CrossRef Full Text | Google Scholar

Penny, W. D., Duzel, E., Miller, K. J., and Ojemann, J. G. (2008). Testing for nested oscillation. J. Neurosci. Methods 174, 50–61. doi: 10.1016/j.jneumeth.2008.06.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Ray, S., and Maunsell, J. H. R. (2015). Do gamma oscillations play a role in cerebral cortex? Trends Cogn. Sci. 19, 78–85. doi: 10.1016/j.tics.2014.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, F., Wibral, M., Singer, W., Aru, J., and Uhlhaas, P. J. (2013). The Phase of Thalamic alpha activity modulates cortical gamma-band activity: evidence from resting-state MEG recordings. J. Neurosci. 33, 17827–17835. doi: 10.1523/jneurosci.5778-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmelin, R., and Hari, R. (1994). Spatiotemporal characteristics of sensorimotor neuromagnetic rhythms related to thumb movement. Neuroscience 60, 537–550. doi: 10.1016/0306-4522(94)90263-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheffer-Teixeira, R., and Tort, A. B. L. (2016). Lack of evidence for cross-frequency phase-phase coupling between theta and gamma oscillations in the hippocampus. bioRxiv. doi: 10.1101/045963

CrossRef Full Text | Google Scholar

Schoffelen, J.-M., Poort, J., Oostenveld, R., and Fries, P. (2011). Selective movement preparation is subserved by selective increases in corticomuscular gamma-band coherence. J. Neurosci. 31, 6750–6758. doi: 10.1523/jneurosci.4882-10.2011

PubMed Abstract | CrossRef Full Text

Schreckenberger, M., Lange-Asschenfeld, C., Lochmann, M., Mann, K., Siessmeier, T., Buchholz, H.-G., et al. (2004). The thalamus as the generator and modulator of EEG alpha rhythm: a combined PET/EEG study with lorazepam challenge in humans. Neuroimage 22, 637–644. doi: 10.1016/j.neuroimage.2004.01.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, M., Warden, M. R., and Miller, E. K. (2009). Phase-dependent neuronal coding of objects in short-term memory. Proc. Natl. Acad. Sci. U.S.A. 106, 21341–21346. doi: 10.1073/pnas.0908193106

PubMed Abstract | CrossRef Full Text | Google Scholar

Spaak, E., Bonnefond, M., Maier, A., Leopold, D. A., and Jensen, O. (2012). Layer-specific entrainment of gamma-band neural activity by the alpha rhythm in monkey visual cortex. Curr. Biol. 22, 2313–2318. doi: 10.1016/j.cub.2012.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Stam, C. J., Vliegen, J. H. R., and Nicolai, J. (1999). Investigation of the dynamics underlying periodic complexes in the EEG. Biol. Cybern. 80, 57–69. doi: 10.1007/s004220050504

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., Gloor, P., Llinás, R. R., Lopes da Silva, F. H., and Mesulam, M. M. (1990). Basic mechanisms of cerebral rhythmic activities. Electroencephalogr. Clin. Neurophysiol. 76, 481–508. doi: 10.1016/0013-4694(90)90001-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Stolk, A., Todorovic, A., Schoffelen, J.-M., and Oostenveld, R. (2013). Online and offline tools for head movement compensation in MEG. Neuroimage 68, 39–48. doi: 10.1016/j.neuroimage.2012.11.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Tal, I., and Abeles, M. (2013). Cleaning MEG artifacts using external cues. J. Neurosci. Methods 217, 31–38. doi: 10.1016/j.jneumeth.2013.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Tort, A. B. L., Scheffer-Teixeira, R., Souza, B. C., Draguhn, A., and Brankačk, J. (2013). Theta-associated high-frequency oscillations (110–160Hz) in the hippocampus and neocortex. Prog. Neurobiol. 100, 1–14. doi: 10.1016/j.pneurobio.2012.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Meij, R., Jacobs, J., and Maris, E. (2015). Uncovering phase-coupled oscillatory networks in electrophysiological data. Hum. Brain Mapp. 36, 2655–2680. doi: 10.1002/hbm.22798

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Meij, R., Kahana, M., and Maris, E. (2012). Phase–amplitude coupling in human electrocorticography is spatially distributed and phase diverse. J. Neurosci. 32, 111–123. doi: 10.1523/jneurosci.4816-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogel, E. K., and Machizawa, M. G. (2004). Neural activity predicts individual differences in visual working memory capacity. Nature 428, 748–751. doi: 10.1038/nature02447

PubMed Abstract | CrossRef Full Text | Google Scholar

von Nicolai, C., Engler, G., Sharott, A., Engel, A. K., Moll, C. K., and Siegel, M. (2014). Corticostriatal coordination through coherent phase-amplitude coupling. J. Neurosci. 34, 5938–5948. doi: 10.1523/jneurosci.5007-13.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Waldert, S., Lemon, R. N., and Kraskov, A. (2013). Influence of spiking activity on cortical local field potentials. J. Physiol. 591, 5291–5303. doi: 10.1113/jphysiol.2013.258228

PubMed Abstract | CrossRef Full Text | Google Scholar