Brain rhythms - as recorded in the local field potential (LFP) or scalp electroencephalogram (EEG) - are believed to play a critical role in coordinating brain networks. By modulating neural excitability, these rhythmic fluctuations provide an effective means to control the timing of neuronal firing (Engel et al., 2001; Buzsáki and Draguhn, 2004). Oscillatory rhythms have been categorized into different frequency bands (e.g., theta [4–10 Hz], gamma [30–80 Hz]) and associated with many functions: the theta band with memory, plasticity, and navigation (Engel et al., 2001); the gamma band with local coupling and competition (Kopell et al., 2000; Börgers et al., 2008). In addition, gamma and high-gamma (80–200 Hz) activity have been identified as surrogate markers of neuronal firing (Rasch et al., 2008; Mukamel et al., 2005; Fries et al., 2001; Pesaran et al., 2002; Whittingstall and Logothetis, 2009; Ray and Maunsell, 2011), observable in the EEG and LFP.

In general, lower frequency rhythms engage larger brain areas and modulate spatially localized fast activity (Bragin et al., 1995; Chrobak and Buzsáki, 1998; von Stein and Sarnthein, 2000; Lakatos et al., 2005; Lakatos et al., 2008). For example, the phase of low frequency rhythms has been shown to modulate and coordinate neural spiking (Vinck et al., 2010; Hyafil et al., 2015b; Fries et al., 2007) via local circuit mechanisms that provide discrete windows of increased excitability. This interaction, in which fast activity is coupled to slower rhythms, is a common type of cross-frequency coupling (CFC). This particular type of CFC has been shown to carry behaviorally relevant information (e.g., related to position [Jensen and Lisman, 2000; Agarwal et al., 2014], memory [Siegel et al., 2009], decision making and coordination [Dean et al., 2012; Pesaran et al., 2008; Wong et al., 2016; Hawellek et al., 2016]). More generally, CFC has been observed in many brain areas (Bragin et al., 1995; Chrobak and Buzsáki, 1998; Csicsvari et al., 2003; Tort et al., 2008; Mormann et al., 2005; Canolty et al., 2006), and linked to specific circuit and dynamical mechanisms (Hyafil et al., 2015b). The degree of CFC in those areas has been linked to working memory, neuronal computation, communication, learning and emotion (Tort et al., 2009; Jensen et al., 2016; Canolty and Knight, 2010; Dejean et al., 2016; Karalis et al., 2016; Likhtik et al., 2014; Jones and Wilson, 2005; Lisman, 2005; Sirota et al., 2008), and clinical disorders (Gordon, 2016; Widge et al., 2017; Voytek and Knight, 2015; Başar et al., 2016; Mathalon and Sohal, 2015), including epilepsy (Weiss et al., 2015). Although the cellular mechanisms giving rise to some neural rhythms are relatively well understood (e.g. gamma [Whittington et al., 2000; Whittington et al., 2011; Mann and Mody, 2010]), the neuronal substrate of CFC itself remains obscure.

Analysis of CFC focuses on relationships between the amplitude, phase, and frequency of two rhythms from different frequency bands. The notion of CFC, therefore, subsumes more specific types of coupling, including: phase-phase coupling (PPC), phase-amplitude coupling (PAC), and amplitude-amplitude coupling (AAC) (Hyafil et al., 2015b). PAC has been observed in rodent striatum and hippocampus (Tort et al., 2008) and human cortex (Canolty et al., 2006), AAC has been observed between the alpha and gamma rhythms in dorsal and ventral cortices (Popov et al., 2018), and between theta and gamma rhythms during spatial navigation (Shirvalkar et al., 2010), and both PAC and AAC have been observed between alpha and gamma rhythms (Osipova et al., 2008). Many quantitative measures exist to characterize different types of CFC, including: mean vector length or modulation index (Canolty et al., 2006; Tort et al., 2010), phase-locking value (Mormann et al., 2005; Lachaux et al., 1999; Vanhatalo et al., 2004), envelope-to-signal correlation (Bruns and Eckhorn, 2004), analysis of amplitude spectra (Cohen, 2008), coherence between amplitude and signal (Colgin et al., 2009), coherence between the time course of power and signal (Osipova et al., 2008), and eigendecomposition of multichannel covariance matrices (Cohen, 2017). Overall, these different measures have been developed from different principles and made suitable for different purposes, as shown in comparative studies (Tort et al., 2010; Cohen, 2008; Penny et al., 2008; Onslow et al., 2011).

Despite the richness of this methodological toolbox, it has limitations. For example, because each method focuses on one type of CFC, the choice of method restricts the type of CFC detectable in data. Applying a method to detect PAC in data with both PAC and AAC may: (i) falsely report no PAC in the data, or (ii) miss the presence of significant AAC in the same data. Changes in the low frequency power can also affect measures of PAC; increases in low frequency power can increase the signal to noise ratio of phase and amplitude variables, increasing the measure of PAC, even when the phase-amplitude coupling remains constant (Aru et al., 2015; van Wijk et al., 2015; Jensen et al., 2016). Furthermore, many experimental or clinical factors (e.g., stimulation parameters, age or sex of subject) can impact CFC in ways that are difficult to characterize with existing methods (Cole and Voytek, 2017). These observations suggest that an accurate measure of PAC would control for confounding variables, including the power of low frequency oscillations.

To that end, we propose here a generalized linear model (GLM) framework to assess CFC between the high-frequency amplitude and, simultaneously, the low frequency phase and amplitude. This formal statistical inference framework builds upon previous work (Kramer and Eden, 2013; Penny et al., 2008; Voytek et al., 2013; van Wijk et al., 2015) to address the limitations of existing CFC measures. In what follows, we show that this framework successfully detects CFC in simulated signals. We compare this method to the modulation index, and show that in signals with CFC dependent on the low-frequency amplitude, the proposed method more accurately detects PAC than the modulation index. We apply this framework to in vivo recordings from human and rodent cortex to show examples of PAC and AAC detected in real data, and how to incorporate new covariates directly into the model framework.

We first examine the performance of the CFC measure through simulation examples. In doing so, we show that the statistics ${\text{𝐑}}_{\text{PAC}}$ and ${\text{𝐑}}_{\text{AAC}}$ accurately detect different types of cross-frequency coupling, increase with the intensity of coupling, and detect weak PAC coupled to the low frequency amplitude. We show that the proposed method is less sensitive to changes in low frequency power, and outperforms an existing PAC measure that lacks dependence on the low frequency amplitude. We conclude with example applications to human and rodent in vivo recordings, and show how to extend the modeling framework to include a new covariate.

The absence of CFC produces no significant detections of coupling

We first consider simulated signals without CFC. To create these signals, we follow the procedure in Materials and methods: Synthetic Time Series with PAC with the modulation intensity set to zero ($I_{\text{PAC}}=0$). In the resulting signals, $A_{\text{high}}$ is approximately constant and does not depend on ${\varphi }_{\text{low}}$ or $A_{\text{low}}$ (Figure 5A). We estimate the ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ model, the $A_{\text{low}}$ model, and the ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ model from these data; we show example fits of the model surfaces in Figure 5B. We observe that the models exhibit small modulations in the estimated high frequency amplitude envelope as a function of the low frequency phase and amplitude.

Figure 5

Download asset Open asset

To assess the distribution of significant R values in the case of no cross-frequency coupling, we simulate 1000 instances of the pink noise signals (each of 20 s) and apply the R measures to each instance, plotting significant R values in Figure 5C. We find that for all 1000 instances, $p_{\text{PAC}}$ and $p_{\text{AAC}}$ are less than 0.05 in only 0.6% and 0.2% of the simulations, respectively, indicating no significant evidence of PAC or AAC, as expected.

We also applied these simulated signals to assess the performance of two standard model comparison procedures for GLMs. Simulating 1000 instances of pink noise signals (each of 20 s) with no induced PAC or AAC, we performed a chi-squared test for nested models (Kramer and Eden, 2016) between models $A_{\text{low}}$ and $A_{\text{low}},{\varphi }_{\text{low}}$, and detected significant PAC (p < 0.05) in 59.7% of simulations. Similarly, performing a chi-squared test for nested models between models ${\varphi }_{\text{low}}$ and $A_{\text{low}},{\varphi }_{\text{low}}$, we detected significant AAC (p < 0.05) in 41.5% of simulations. Using an AIC-based model comparison, we found a decrease in AIC from the $A_{\text{low}}$ model to the $A_{\text{low}},{\varphi }_{\text{low}}$ model (consistent with significant PAC) in 98.6% of simulations, and a decrease in AIC from the ${\varphi }_{\text{low}}$ model to the $A_{\text{low}},{\varphi }_{\text{low}}$ model (consistent with significant AAC) in 87.2% of simulations. By contrast, we rarely detect significant PAC (<0.6% of simulations) or AAC (<0.2% of simulations) in the pink noise signals using the two statistics ${\text{𝐑}}_{\text{PAC}}$ and ${\text{𝐑}}_{\text{AAC}}$ implemented here. We conclude that, in this modeling regime, two deviance-based model comparison procedures for GLMs are less robust measures of significant PAC and AAC.

${\text{𝐑}}_{\text{PAC}}$ and modulation index are both sensitive to weak modulations

To investigate the ability of the proposed method and the modulation index to detect weak coupling between the low frequency phase and high frequency amplitude, we perform the following simulations. For each intensity value $I_{\text{PAC}}$ between 0 and 0.5 (in steps of 0.025), we simulate 1000 signals (see Materials and methods) and compute ${\text{𝐑}}_{\text{PAC}}$ and a measure of PAC in common use: the modulation index MI (Tort et al., 2010) (Figure 6). We find that both MI and ${\text{𝐑}}_{\text{PAC}}$, while small, increase with $I_{\text{PAC}}$; in this way, both measures are sensitive to small values of $I_{\text{PAC}}$. However, we note that ${\text{𝐑}}_{\text{PAC}}$ is not significant for very small intensity values ($I_{\text{PAC}}\le 0.3$), while MI is significant at these small intensities. Significant ${\text{𝐑}}_{\text{PAC}}$ appears when the MI exceeds 0.7 × 10^-3, a value below the range of MI values detected in many existing studies (Tort et al., 2008; Zhong et al., 2017; Jackson et al., 2019; Axmacher et al., 2010; Tort et al., 2018). We conclude that, while the modulation index may be more sensitive than ${\text{𝐑}}_{\text{PAC}}$ to very weak phase-amplitude coupling, ${\text{𝐑}}_{\text{PAC}}$ can detect phase-amplitude coupling at MI values consistent with those observed in the literature.

Figure 6

Download asset Open asset

The proposed method is less affected by fluctuations in low-frequency amplitude and AAC

Increases in low frequency power can increase measures of phase-amplitude coupling, although the underlying PAC remains unchanged (Aru et al., 2015; Cole and Voytek, 2017). Characterizing the impact of this confounding effect is important both to understand measure performance and to produce accurate interpretations of analyzed data. To examine this phenomenon, we perform the following simulation. First, we simulate a signal $V$ with fixed PAC (intensity $I_{\text{PAC}}=1$, see Materials and methods). Second, we filter $V$ into its low and high frequency components $V_{\text{low}}$ and $V_{\text{high}}$, respectively. Then, we create a new signal $V^{*}$ as follows:

(8) ${\displaystyle {\displaystyle V^{\ast }=2V_{\text{low}}+V_{\text{high}}+V_{\text{noise}},}}$

where $V_{\text{noise}}$ is a pink noise term (see Materials and methods). We note that we only alter the low frequency component of $V$ and do not alter the PAC. To analyze the PAC in this new signal we compute ${\text{𝐑}}_{\text{PAC}}$ and MI.

We show in Figure 7 population results (1000 realizations each of the simulated signals $V$ and $V^{*}$) for the R and MI values. We observe that increases in the amplitude of $V_{\text{low}}$ produce increases in MI and ${\text{𝐑}}_{\text{PAC}}$. However, this increase is more dramatic for MI than for ${\text{𝐑}}_{\text{PAC}}$; we note that the distributions of ${\text{𝐑}}_{\text{PAC}}$ almost completely overlap (Figure 7A), while the distribution of MI shifts to larger values when the amplitude of ${\displaystyle V_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ increases (Figure 7B). We conclude that the statistic ${\text{𝐑}}_{\text{PAC}}$ — which includes the low frequency amplitude as a predictor in the GLM — is more robust to increases in low frequency power than a method that only includes the low frequency phase.

Figure 7

Download asset Open asset

We also investigate the effect of increases in amplitude-amplitude coupling (AAC) on the two measures of PAC. As before, we simulate a signal $V$ with fixed PAC (intensity $I_{\text{PAC}}=1$) and no AAC (intensity $I_{\text{AAC}}=0$). We then simulate a second signal $V^{*}$ with the same fixed PAC as $V$, and with additional AAC (intensity $I_{\text{AAC}}=10$). We simulate 1000 realizations of $V$ and $V^{*}$ and compute the corresponding ${\text{𝐑}}_{\text{PAC}}$ and MI values. We observe that the increase in AAC produces a small increase in the distribution of ${\text{𝐑}}_{\text{PAC}}$ values (Figure 7C), but a large increase in the distribution of MI values (Figure 7D). We conclude that the statistic ${\text{𝐑}}_{\text{PAC}}$ is more robust to increases in AAC than MI.

These simulations show that at a fixed, non-zero PAC, the modulation index increases with increased ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ and AAC. We now consider the scenario of increased ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ and AAC in the absence of PAC. To do so, we simulate 1000 signals of 200 s duration, with no PAC (intensity ${\displaystyle I_{\mathrm{P}\mathrm{A}\mathrm{C}}=0}$). For each signal, at time 100 s (i.e., the midpoint of the simulation) we increase the low frequency amplitude by a factor of 10 (consistent with observations from an experiment in rodent cortex, as described below), and include AAC between the low and high frequency signals (intensity ${\displaystyle I_{\mathrm{A}\mathrm{A}\mathrm{C}}=0}$ for ${\displaystyle t<100\mathrm{s}}$ and intensity ${\displaystyle I_{\mathrm{A}\mathrm{A}\mathrm{C}}=2}$ for ${\displaystyle t\ge 100\mathrm{s}}$). We find that, in the absence of PAC, ${\text{𝐑}}_{\text{PAC}}$ detects significant PAC (p<0.05) in 0.4% of the simulated signals, while MI detects significant PAC in 34.3% of simulated signals. We conclude that in the presence of increased low frequency amplitude and amplitude-amplitude coupling, MI may detect PAC where none exists, while ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$, which accounts for fluctuations in low frequency amplitude, does not.

Sparse PAC is detected when coupled to the low frequency amplitude

While the modulation index has been successfully applied in many contexts (Canolty and Knight, 2010; Hyafil et al., 2015b), instances may exist where this measure is not optimal. For example, because the modulation index was not designed to account for the low frequency amplitude, it may fail to detect PAC when ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ depends not only on ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, but also on ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$. For example, since the modulation index considers the distribution of $A_{\text{high}}$ at all observed values of ${\varphi }_{\text{low}}$, it may fail to detect coupling events that occur sparsely at only a subset of appropriate ${\displaystyle {\varphi }_{\text{low}}}$ occurrences. ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$, on the other hand, may detect these sparse events if these events are coupled to ${\displaystyle A_{\text{low}}}$, as ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$ accounts for fluctuations in low frequency amplitude. To illustrate this, we consider a simulation scenario in which PAC occurs sparsely in time.

We create a signal ${\displaystyle V}$ with PAC, and corresponding modulation signal M with intensity value ${\displaystyle I_{\mathrm{P}\mathrm{A}\mathrm{C}}=1.0}$ (see Materials and methods, Figure 8A–B). We then modify this signal to reduce the number of PAC events in a way that depends on ${\displaystyle A_{\text{low}}}$. To do so, we preserve PAC at the peaks of ${\displaystyle V_{\text{low}}}$ (i.e., when ${\displaystyle {\varphi }_{\text{low}}=0}$), but now only when these peaks are large, more specifically in the top 5% of peak values.

Figure 8

Download asset Open asset

We define a threshold value ${\displaystyle T}$ to be the 95^th quantile of the peak ${\displaystyle V_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ values, and modify the modulation signal M as follows. When M exceeds 1 (i.e., when ${\displaystyle {\varphi }_{\text{low}}=0}$) and the low frequency amplitude exceeds ${\displaystyle T}$ (i.e., ${\displaystyle A_{\text{low}}\ge T}$), we make no change to M. Alternatively, when M exceeds one and the low frequency amplitude lies below ${\displaystyle T}$ (i.e., ${\displaystyle A_{\text{low}}<T)}$, we decrease M to 1 (Figure 8C). In this way, we create a modified modulation signal ${\displaystyle {\mathbf{\text{M}}}_1}$ such that in the resulting signal ${\displaystyle V_1}$, when ${\displaystyle {\varphi }_{\text{low}}=0}$ and ${\displaystyle A_{\text{low}}}$ is large enough, ${\displaystyle A_{\text{high}}}$ is increased; and when ${\displaystyle {\varphi }_{\text{low}}=0}$ and ${\displaystyle A_{\text{low}}}$ is not large enough, there is no change to ${\displaystyle A_{\text{high}}}$. This signal ${\displaystyle V_1}$ hence has fewer phase-amplitude coupling events than the number of times ${\displaystyle {\varphi }_{\text{low}}=0}$.

We generate 1000 realizations of the simulated signals ${\displaystyle V_1}$, and compute ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ and MI. We find that while MI detects significant PAC in only 37% of simulations, ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$ detects significant PAC in 72% of simulations. In this case, although the PAC occurs infrequently, these occurrences are coupled to ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, and ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$, which accounts for changes in ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, successfully detects these events much more frequently. We conclude that when the PAC is dependent on ${\displaystyle A_{\text{low}}}$, ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ more accurately detects these sparse coupling events.

The CFC model detects simultaneous PAC and AAC missed in an existing method

To further illustrate the utility of the proposed method, we consider another scenario in which ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ impacts the occurrence of PAC. More specifically, we consider a case in which ${\displaystyle A_{\text{high}}}$ increases at a fixed low frequency phase for high values of ${\displaystyle A_{\text{low}}}$, and ${\displaystyle A_{\text{high}}}$ decreases at the same phase for small values of ${\displaystyle A_{\text{low}}}$. In this case, we expect that the modulation index may fail to detect the coupling because the distribution of ${\displaystyle A_{\text{high}}}$ over ${\displaystyle {\varphi }_{\text{low}}}$ would appear uniform when averaged over all values of ${\displaystyle A_{\text{low}}}$; the dependence of ${\displaystyle A_{\text{high}}}$ on ${\displaystyle {\varphi }_{\text{low}}}$ would only become apparent after accounting for ${\displaystyle A_{\text{low}}}$.

To implement this scenario, we consider the modulation signal M (see Materials and methods) with an intensity value ${\displaystyle I_{\mathrm{P}\mathrm{A}\mathrm{C}}=1}$. We consider all peaks of ${\displaystyle A_{\text{low}}}$ and set the threshold ${\displaystyle T}$ to be the 50^th quantile (Figure 9A). We then modify the modulation signal M as follows. When M exceeds 1 (i.e., when ${\displaystyle {\varphi }_{\text{low}}=0}$) and the low frequency amplitude exceeds ${\displaystyle T}$ (i.e., ${\displaystyle A_{\text{low}}\ge T}$), we make no change to M. Alternatively, when M exceeds one and the low frequency amplitude lies below ${\displaystyle T}$ (i.e. ${\displaystyle A_{\text{low}}<T)}$, we decrease M to 0 (Figure 9B). In this way, we create a modified modulation signal M such that when ${\displaystyle {\varphi }_{\text{low}}=0}$ and ${\displaystyle A_{\text{low}}}$ is large enough, ${\displaystyle A_{\text{high}}}$ is increased; and when ${\displaystyle {\varphi }_{\text{low}}=0}$ and ${\displaystyle A_{\text{low}}}$ is small enough, ${\displaystyle A_{\text{high}}}$ is decreased (Figure 9C).

Figure 9

Download asset Open asset

Using this method, we simulate 1000 realizations of this signal, and calculate MI and ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ for each signal (Figure 9D). We find that ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$ detects significant PAC in nearly all (96%) of the simulations, while MI detects significant PAC in only 58% of the simulations. We conclude that, in this simulation, ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$ more accurately detects PAC coupled to low frequency amplitude.

A simple stochastic spiking neural model illustrates the utility of the proposed method

In the previous simulations, we created synthetic data without a biophysically principled generative model. Here we consider an alternative simulation strategy with a more direct connection to neural dynamics. While many biophysically motivated models of cross-frequency coupling exist (Sase et al., 2017; Chehelcheraghi et al., 2017; Sotero, 2016; Hyafil et al., 2015a; Lepage and Vijayan, 2015; Onslow et al., 2014; Fontolan et al., 2013; Malerba and Kopell, 2013; Jirsa and Müller, 2013; Spaak et al., 2012; Wulff et al., 2009; Tort et al., 2007), we consider here a relatively simple stochastic spiking neuron model (Aljadeff et al., 2016). In this stochastic model, we generate a spike train (${\displaystyle V_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$) in which an externally imposed signal ${\displaystyle V_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ modulates the probability of spiking as a function of ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ and ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$. We note that high frequency activity is thought to represent the aggregate spiking activity of local neural populations (Ray and Maunsell, 2011; Buzsáki and Wang, 2012; Ray et al., 2008a; Jia and Kohn, 2011); while here we simulate the activity of a single neuron, the spike train still produces temporally focal events of high frequency activity. In this framework, we allow the target phase (${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast }}$) modulating ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ to change as a function of ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$: when ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ is large, the probability of spiking is highest near ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}=\pm \pi }$, and when ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ is small, the probability of spiking is highest near ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}=0.}$ More precisely, we define ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast }}$ as

${\displaystyle {\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast }=\pi (1+A_{\mathrm{l}\mathrm{o}\mathrm{w}})}}$

where ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ is a sinusoid oscillating between 1 and 2 with period 0.1 Hz. We define the spiking probability, ${\displaystyle \lambda }$, as

${\displaystyle {\displaystyle \lambda ={\lambda }_0\exp [-\left(1+\frac{s({\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}-{\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast }{)}^2}{2{\sigma }^2}\right)],}}$

where ${\displaystyle \sigma =0.01}$, ${\displaystyle s(\varphi )}$ is a triangle wave, and we choose ${\displaystyle {\lambda }_0}$ so that the maximum value of ${\displaystyle \lambda }$ is 2. We note that the spiking probability ${\displaystyle \lambda }$ is zero except near times when the phase of the low frequency signal ${\displaystyle ({\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}})}$ is near ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast }}$. We then define ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ as:

${\displaystyle {\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}=S+n,}}$

where ${\displaystyle S}$ is the binary sequence generated by the stochastic spiking neuron model, and ${\displaystyle n}$ is Gaussian noise with mean zero and standard deviation 0.1. In this scenario, the distribution of ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ over ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ appears uniform when averaged over all values of ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$. We therefore expect the modulation index to remain small, despite the presence of PAC with maximal phase dependent on ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$. However, we expect that ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$, which accounts for fluctuations in low frequency amplitude, will detect this PAC. We show an example signal from this simulation in Figure 10A. As expected, we find that ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$ detects PAC (${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}=0.172}$, ${\displaystyle p=0.02}$); we note that the (${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$) surface exhibits a single peak near ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}=0}$ at small values of ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, and at ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}=\pm \pi }$ at large value of ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ (Figure 10B). The (${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$) surface deviates significantly from the ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ surface, resulting in a large ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$ value. However, the non-uniform shape of the (${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$) surface is lost when we fail to account for ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$. In this scenario, the distribution of ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ over ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ appears uniform, resulting in a low MI value (Figure 10C).

Figure 10

Download asset Open asset

Application to in vivo human seizure data

To evaluate the performance of the proposed method on in vivo data, we first consider an example recording from human cortex during a seizure (see Materials and methods: Human subject data). Visual inspection of the LFP data (Figure 11A) reveals the emergence of large amplitude voltage fluctuations during the approximately 80 s seizure. We compute the spectrogram over the entire seizure, using windows of width 0.8 s with 0.002 s overlap, and identify a distinct 10 s interval of increased power in the 4–7 Hz band (Figure 11B). We analyze this section of the voltage trace ${\displaystyle V}$, filtering into ${\displaystyle V_{\text{high}}}$ (100–140 Hz) and ${\displaystyle V_{\text{low}}}$ (4–7 Hz), and extracting ${\displaystyle A_{\text{high}}}$, ${\displaystyle A_{\text{low}}}$, and ${\displaystyle {\varphi }_{\text{low}}}$ as in Methods (Figure 11C). Visual inspection reveals the occurrence of large amplitude, low frequency oscillations and small amplitude, high frequency oscillations.

Figure 11

Download asset Open asset

We find during this interval significant phase-amplitude coupling computed using ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ (${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}=1.55}$, ${\displaystyle p_{\text{PAC}}=0.005}$, Figure 12), and using the modulation index (${\displaystyle \mathbf{\text{MI}}=0.03}$, ${\displaystyle p_{\text{MI}}=5.0\times {10}^{-4}}$). To examine the phase-amplitude coupling in more detail, we isolate a 2 s segment (Figure 11D) and display the signal ${\displaystyle V}$, the high frequency signal ${\displaystyle V_{\text{high}}}$, the low frequency phase ${\displaystyle {\varphi }_{\text{low}}}$, and the low frequency amplitude ${\displaystyle A_{\text{low}}}$. We observe that when ${\displaystyle {\varphi }_{\text{low}}}$ is near ${\displaystyle \pi }$, the amplitude of ${\displaystyle V_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ tends to increase, consistent with the presence of PAC and a significant value of ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}}$ and MI.

Figure 12

Download asset Open asset

We also find significant amplitude-amplitude coupling computed using ${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}}$ (${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}=0.85}$, ${\displaystyle p_{\text{AAC}}=0.005}$). Comparing ${\displaystyle A_{\text{high}}}$ and ${\displaystyle A_{\text{low}}}$ over the 10 s interval (each smoothed using a 1 s moving average filter and normalized), we observe that both ${\displaystyle A_{\text{high}}}$ and ${\displaystyle A_{\text{low}}}$ steadily increase over the duration of the interval (Figure 11E).

Using the flexibility of GLMs to improve detection of phase-amplitude coupling in vivo

One advantage of the proposed framework is its flexibility: covariates are easily added to the generalized linear model and tested for significance. For example, we could include covariates for trial, sex, and stimulus parameters and explore their effects on PAC, AAC, or both.

Here, we illustrate this flexibility through continued analysis of the rodent data. We select a single electrode recording from these data, and hypothesize that the condition, either pre-stimulation or post-stimulation, affects the coupling. To incorporate this new covariate into the framework, we consider the concatenated voltage recordings from the pre-stimulation condition ${\displaystyle V_{\mathrm{p}\mathrm{r}\mathrm{e}}}$ and the post-stimulation condition ${\displaystyle V_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}}$:

${\displaystyle {\displaystyle V=[V_{\mathrm{p}\mathrm{r}\mathrm{e}},V_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}].}}$

From ${\displaystyle V}$, we obtain the corresponding high frequency signal ${\displaystyle V_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ and low frequency signal ${\displaystyle V_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, and subsequently the high frequency amplitude ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$, low frequency phase ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, and low frequency amplitude ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$. We use these data to generate two new models:

(9) ${\displaystyle {\displaystyle A_{\text{high}}|{\varphi }_{\text{low}},A_{\text{low}},P\sim \text{Gamma}[\mu ,\nu ],}}$

${\displaystyle {\displaystyle \log \mu =\sum _{k=1}^n{\beta }_k{f}_k({\varphi }_{\text{low}})+{\beta }_{n+1}{A}_{\text{low}}+{\beta }_{n+2}{A}_{\text{low}}\sin ({\varphi }_{\text{low}})+{\beta }_{n+3}{A}_{\text{low}}\cos ({\varphi }_{\text{low}})+P(\sum _{j=1}^n{\beta }_{n+3+j}{f}_j({\varphi }_{\text{low}})+{\beta }_{2n+4}{A}_{\text{low}})}}$

(10) ${\displaystyle {\displaystyle A_{\text{high}}|{\varphi }_{\text{low}},A_{\text{low}},P\sim \text{Gamma}[\mu ,\nu ]}}$

${\displaystyle {\displaystyle \log \mu =\sum _{k=1}^n{\beta }_k{f}_k({\varphi }_{\text{ low}})+{\beta }_{n+1}{A}_{\text{low}}+{\beta }_{n+2}{A}_{\text{low}}\sin ({\varphi }_{\text{low }})+{\beta }_{n+3}{A}_{\text{low}}\cos ({\varphi }_{\text{low}})+P({\beta }_{n+4}{A}_{\text{ low}}),}}$

where ${\displaystyle P}$ is an indicator function specifying whether the signal is in the pre-stimulation (${\displaystyle P=0}$) or post-stimulation (${\displaystyle P=1}$) condition. The effect of the indicator function is to include the effect of stimulus condition on the high frequency amplitude. The models in Equations 9 and 10 now include the effect of low frequency amplitude, low frequency phase, and condition on high frequency amplitude. To determine whether the condition has an effect on PAC, we test whether the term ${\displaystyle P(\sum _{j=1}^n{\beta }_{n+3+j}{f}_j({\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}))}$ in Equation 9 is significant, that is whether there is a significant difference between the models in Equations 9 and 10. If the difference between the two models is very small, we gain no improvement in modeling ${\displaystyle A_{\text{high}}}$ by including the interaction between ${\displaystyle P}$ and ${\displaystyle {\varphi }_{\text{low}}}$. In that case, the impact of ${\displaystyle {\varphi }_{\text{low}}}$ on ${\displaystyle A_{\text{high}}}$ can be modeled without considering stimulus condition ${\displaystyle P}$, that is the impact of stimulus condition on PAC is negligible.

To measure the difference between the models in Equations 9 and 10, we construct a surface ${\displaystyle S_{P{\varphi }_{\text{low}}}}$ from the model in Equation 9, and a surface ${\displaystyle S_P}$ from the model in Equation 10 in the (${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$, P) space, assessing the models at ${\displaystyle \mathit{\text{P}}=1}$. We compute ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC},\text{condition}}}$, which measures the impact of stimulus condition on PAC, as:

(11) ${\displaystyle {\displaystyle {\mathbf{\text{R}}}_{\text{PAC},\text{ condition}}=max[\mathrm{a}\mathrm{b}\mathrm{s}[1-{\mathrm{S}}_{\mathrm{P}}/{\mathrm{S}}_{\mathrm{P}{\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}]].}}$

We find for the example rodent data an ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C},\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}}$ value of 0.3608, with a p-value of 0.0005. Hence, we find evidence for a significant effect of stimulus on PAC.

To further explore this assessment of stimulus condition on PAC, we simulate 1000 instances of a 40 s signal divided into two conditions: no PAC for the first 20 s (${\displaystyle I_{\text{PAC}}=0}$) and non-zero PAC for the final 20 s ${\displaystyle (I_{\text{PAC}}=1)}$. We design this simulation to mimic an increase in PAC from pre-stimulation to post-stimulation (Figure 13A). Using the models in Equations 9 and 10, and computing ${\displaystyle {\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C},\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}}$, we find ${\displaystyle p<0.05}$ for 100% of simulated signals. We also simulate 1000 instances of a 40 s signal with no PAC (${\displaystyle I_{\text{PAC}}=0}$) for the entire 40 s, that is PAC does not change from pre-stimulation to post-stimulation (Figure 13B), and find in this case ${\displaystyle p<0.05}$ for only 4.6% of simulations. Finally, we simulate 1000 instances of a 40 s signal with fixed PAC ${\displaystyle (I_{\text{PAC}}=1)}$, and with a doubling of the low frequency amplitude occuring at 20 s (i.e., pre-stimulation the low frequency amplitude is 1, and post-stimulation the low frequency amplitude is 2). We find ${\displaystyle p<0.05}$ for only 3.6% of simulations. We conclude that this method effectively determines whether stimulation condition significantly changes PAC.

Figure 13

Download asset Open asset

This example illustrates the flexibility of the statistical modeling framework. Extending this framework is straightforward, and new extensions allow a common principled approach to test the impact of new predictors. Here we considered an indicator function that divides the data into two states (pre- and post-stimulation). We note that the models are easily extended to account for multiple discrete predictors such as gender and participation in a drug trial, or for continuous predictors such as age and time since stimulus.

In this paper, we proposed a new method for measuring cross-frequency coupling that accounts for both phase-amplitude coupling and amplitude-amplitude coupling, along with a principled statistical modeling framework to assess the significance of this coupling. We have shown that this method effectively detects CFC, both as PAC and AAC, and is more sensitive to weak PAC obscured by or coupled to low-frequency amplitude fluctuations. Compared to an existing method, the modulation index (Tort et al., 2010), the newly proposed method more accurately detects scenarios in which PAC is coupled to the low-frequency amplitude. Finally, we applied this method to in vivo data to illustrate examples of PAC and AAC in real systems, and show how to extend the modeling framework to include a new covariate.

One of the most important features of the new method is an increased ability to detect weak PAC coupled to AAC. For example, when sparse PAC events occur only when the low frequency amplitude (${\displaystyle A_{\text{low}}}$) is large, the proposed method detects this coupling while another method not accounting for ${\displaystyle A_{\text{low}}}$ misses it. While PAC often occurs in neural data, and has been associated with numerous neurological functions (Canolty and Knight, 2010; Hyafil et al., 2015b), the simultaneous occurrence of PAC and AAC is less well studied (Osipova et al., 2008). Here, we showed examples of simultaneous PAC and AAC recorded from human cortex during a seizure, and we note that this phenomena has been simulated in other works (Mazzoni et al., 2010).

While the exact mechanisms that support CFC are not well understood (Hyafil et al., 2015b), the general mechanisms of low and high frequency rhythms have been proposed. Low frequency rhythms are associated with the aggregate activity of large neural populations and modulations of neuronal excitability (Engel et al., 2001; Varela et al., 2001; Buzsáki and Draguhn, 2004), while high frequency rhythms provided a surrogate measure of neuronal spiking (Rasch et al., 2008; Mukamel et al., 2005; Fries et al., 2001; Pesaran et al., 2002; Whittingstall and Logothetis, 2009; Ray and Maunsell, 2011; Ray et al., 2008b). These two observations provide a physical interpretation for PAC: when a low frequency rhythm modulates the excitability of a neural population, we expect spiking to occur (i.e., an increase in ${\displaystyle A_{\text{high}}}$) at a particular phase of the low frequency rhythm (${\displaystyle {\varphi }_{\text{low}}}$) when excitation is maximal. These notions also provide a physical interpretation for AAC: increases in ${\displaystyle A_{\text{low}}}$ produce larger modulations in neural excitability, and therefore increased intervals of neuronal spiking (i.e., increases in ${\displaystyle A_{\text{high}}}$). Alternatively, decreases in ${\displaystyle A_{\text{low}}}$ reduce excitability and neuronal spiking (i.e., decreases in ${\displaystyle A_{\text{high}}}$).

The function of concurrent PAC and AAC, both for healthy brain function and during a seizure as illustrated here, is not well understood. As PAC occurs normally in healthy brain signals, for example during working memory, neuronal computation, communication, learning and emotion (Tort et al., 2009; Jensen et al., 2016; Canolty and Knight, 2010; Dejean et al., 2016; Karalis et al., 2016; Likhtik et al., 2014; Jones and Wilson, 2005; Lisman, 2005; Sirota et al., 2008), these preliminary results may suggest a pathological aspect of strong AAC occurring concurrently with PAC.

Proposed functions of PAC include multi-item encoding, long-distance communication, and sensory parsing (Hyafil et al., 2015b). Each of these functions takes advantage of the low frequency phase, encoding different objects or pieces of information in distinct phase intervals of ${\displaystyle {\varphi }_{\text{low}}}$. PAC can be interpreted as a type of focused attention; ${\displaystyle A_{\text{high}}}$ modulation occurring only in a particular interval of ${\displaystyle {\varphi }_{\text{low}}}$ organizes neural activity - and presumably information - into discrete packets of time. Similarly, a proposed function of AAC is to encode the number of represented items, or the amount of information encoded in the modulated signal (Hyafil et al., 2015b). A pathological increase in AAC may support the transmission of more information than is needed, overloading the communication of relevant information with irrelevant noise. The attention-based function of PAC, that is having reduced high frequency amplitude at phases not containing the targeted information, may be lost if the amplitude of the high frequency oscillation is increased across wide intervals of low frequency phase.

Like all measures of CFC, the proposed method possesses specific limitations. We discuss five limitations here. First, the choice of spline basis to represent the low frequency phase may be inaccurate, for example if the PAC changes rapidly with ${\displaystyle {\varphi }_{\text{low}}}$. Second, the value of ${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}}$ depends on the range of ${\displaystyle A_{\text{low}}}$ observed. This is due to the linear relationship between ${\displaystyle A_{\text{low}}}$ and ${\displaystyle A_{\text{high}}}$ in the ${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$ model, which causes the maximum distance between the surfaces ${\displaystyle S_{A_{\mathrm{l}\mathrm{o}\mathrm{w}}}}$ and ${\displaystyle S_{A_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}}$ to occur at the largest or smallest value of ${\displaystyle A_{\text{low}}}$. To mitigate the impact of extreme ${\displaystyle A_{\text{low}}}$ values on ${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}}$, we evaluate the surfaces ${\displaystyle S_{A_{\mathrm{l}\mathrm{o}\mathrm{w}}}}$ and ${\displaystyle S_{A_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}}$ over the 5^th to 95^th quantiles of ${\displaystyle A_{\text{low}}}$. We note that an alternative metric of AAC could instead evaluate the slope of the ${\displaystyle S_{A_{\mathrm{l}\mathrm{o}\mathrm{w}}}}$ surface; to maintain consistency of the PAC and AAC measures, we chose not to implement this alternative measure here. Third, the frequency bands for ${\displaystyle V_{\text{high}}}$ and ${\displaystyle V_{\text{low}}}$ must be established before R values are calculated. Hence, if the wrong frequency bands are chosen, coupling may be missed. It is possible, though computationally expensive, to scan over all reasonable frequency bands for both ${\displaystyle V_{\text{high}}}$ and ${\displaystyle V_{\text{low}}}$, calculating R values for each frequency band pair. Fourth, we note that the proposed modeling framework assumes the data contain approximately sinusoidal signals, which have been appropriately isolated for analysis. In general, CFC measures are sensitive to non-sinusoidal signals, which may confound interpretation of cross-frequency analyses (Cole and Voytek, 2017; Kramer et al., 2008; Aru et al., 2015). While the modeling framework proposed here does not directly account for the confounds introduced by non-sinusoidal signals, the inclusion of additional predictors (e.g. detections of sharp changes in the unfiltered data) in the model may help mitigate these effects. Fifth, we simulate time series with known PAC and AAC, and then test whether the proposed analysis framework detects this coupling. The simulated relationships between ${\displaystyle A_{\text{high}}}$ and (${\displaystyle {\varphi }_{\text{low}}}$,${\displaystyle A_{\text{low}}}$) may result in time series with simpler structure than those observed in vivo. For example, a latent signal may drive both ${\displaystyle A_{\text{high}}}$ and ${\displaystyle {\varphi }_{\text{low}}}$, and in this way establish nonlinear relationships between the two observables ${\displaystyle A_{\text{high}}}$ and ${\displaystyle {\varphi }_{\text{low}}}$. We note that, if this were the case, the latent signal could also be incorporated in the statistical modeling framework (Yousefi et al., 2019).

We chose the statistics ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ and ${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}}$ for two reasons. First, we found that two common methods of model comparison for GLMs provide less robust measures of significance than ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ and ${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}}$. While the statistics ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ and ${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}}$ are less powerful than standard model comparison tests, the large amount of data typically assessed in CFC analysis may compensate for this loss. We showed that the statistics ${\displaystyle {\mathbf{\text{R}}}_{\text{PAC}}}$ and ${\displaystyle {\mathbf{\text{R}}}_{\text{AAC}}}$ performed well in simulations, and we note that these statistics are directly interpretable. While many model comparison methods exist - and another method may provide specific advantages - we found that the framework implemented here is sufficiently powerful, interpretable, and robust for real-world neural data analysis.

The proposed method can easily be extended by inclusion of additional predictors in the GLM. Polynomial ${\displaystyle A_{\text{low}}}$ predictors, rather than the current linear ${\displaystyle A_{\text{low}}}$ predictors, may better capture the relationship between ${\displaystyle A_{\text{low}}}$ and ${\displaystyle A_{\text{high}}}$. One could also include different types of covariates, for example classes of drugs administered to a patient, or time since an administered stimulus during an experiment. To capture more complex relationships between the predictors (${\displaystyle A_{\mathrm{l}\mathrm{o}\mathrm{w}}}$, ${\displaystyle {\varphi }_{\mathrm{l}\mathrm{o}\mathrm{w}}}$) and ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$, the GLM could be replaced by a more general form of Generalized Additive Model (GAM). Choosing GAMs would remove the restriction that the conditional mean ${\displaystyle A_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}}$ must be linear in each of the model parameters (which would allow us to estimate knot locations directly from the data, for example), at the cost of greater computational time to estimate these parameters. The code developed to implement the method is flexible and modular, which facilitates modifications and extensions motivated by the particular data analysis scenario. This modular code, available at https://github.com/Eden-Kramer-Lab/GLM-CFC, also allows the user to change latent assumptions, such as choice of frequency bands and filtering method. The code is freely available for reuse and further development.

Rhythms, and particularly the interactions of different frequency rhythms, are an important component for a complete understanding of neural activity. While the mechanisms and functions of some rhythms are well understood, how and why rhythms interact remains uncertain. A first step in addressing these uncertainties is the application of appropriate data analysis tools. Here we provide a new tool to measure coupling between different brain rhythms: the method utilizes a statistical modeling framework that is flexible and captures subtle differences in cross-frequency coupling. We hope that this method will better enable practicing neuroscientists to measure and relate brain rhythms, and ultimately better understand brain function and interactions.

In vivo human data available at https://github.com/Eden-Kramer-Lab/GLM-CFC (copy archived at https://github.com/elifesciences-publications/GLM-CFC). In vivo rat data available at https://github.com/tne-lab/cl-example-data (copy archived at https://github.com/elifesciences-publications/cl-example-data).

1. 1. Agarwal G
   2. Stevenson IH
   3. Berényi A
   4. Mizuseki K
   5. Buzsáki G
   6. Sommer FT

   (2014) Spatially distributed local fields in the Hippocampus encode rat position

   Science 344:626–630.

   https://doi.org/10.1126/science.1250444

   • PubMed
   • Google Scholar
2. 1. Aljadeff J
   2. Lansdell BJ
   3. Fairhall AL
   4. Kleinfeld D

   (2016) Analysis of neuronal spike trains, deconstructed

   Neuron 91:221–259.

   https://doi.org/10.1016/j.neuron.2016.05.039

   • PubMed
   • Google Scholar
3. 1. Aru J
   2. Aru J
   3. Priesemann V
   4. Wibral M
   5. Lana L
   6. Pipa G
   7. Singer W
   8. Vicente R

   (2015) Untangling cross-frequency coupling in neuroscience

   Current Opinion in Neurobiology 31:51–61.

   https://doi.org/10.1016/j.conb.2014.08.002

   • PubMed
   • Google Scholar
4. 1. Axmacher N
   2. Henseler MM
   3. Jensen O
   4. Weinreich I
   5. Elger CE
   6. Fell J

   (2010) Cross-frequency coupling supports multi-item working memory in the human Hippocampus

   PNAS 107:3228–3233.

   https://doi.org/10.1073/pnas.0911531107

   • PubMed
   • Google Scholar
5. 1. Başar E
   2. Schmiedt-Fehr C
   3. Mathes B
   4. Femir B
   5. Emek-Savaş DD
   6. Tülay E
   7. Tan D
   8. Düzgün A
   9. Güntekin B
   10. Özerdem A
   11. Yener G
   12. Başar-Eroğlu C

   (2016) What does the broken brain say to the neuroscientist? oscillations and connectivity in schizophrenia, Alzheimer’s disease, and bipolar disorder

   International Journal of Psychophysiology 103:135–148.

   https://doi.org/10.1016/j.ijpsycho.2015.02.004

   • PubMed
   • Google Scholar
6. 1. Blackwood E
   2. Lo M
   3. Widge SA

   (2018) Continuous phase estimation for phase-locked neural stimulation using an autoregressive model for signal prediction

   Conference of the IEEE Engineering in Medicine and Biology Society 2018:4736–4739.

   https://doi.org/10.1109/EMBC.2018.8513232

   • Google Scholar
7. 1. Börgers C
   2. Epstein S
   3. Kopell NJ

   (2008) Gamma oscillations mediate stimulus competition and attentional selection in a cortical network model

   PNAS 105:18023–18028.

   https://doi.org/10.1073/pnas.0809511105

   • PubMed
   • Google Scholar
8. 1. Bragin A
   2. Jandó G
   3. Nádasdy Z
   4. Hetke J
   5. Wise K
   6. Buzsáki G

   (1995) Gamma (40-100 hz) oscillation in the Hippocampus of the behaving rat

   The Journal of Neuroscience 15:47–60.

   https://doi.org/10.1523/JNEUROSCI.15-01-00047.1995

   • PubMed
   • Google Scholar
9. 1. Bruns A
   2. Eckhorn R

   (2004) Task-related coupling from high- to low-frequency signals among visual cortical Areas in human subdural recordings

   International Journal of Psychophysiology 51:97–116.

   https://doi.org/10.1016/j.ijpsycho.2003.07.001

   • PubMed
   • Google Scholar
10. 1. Buzsáki G
    2. Draguhn A

    (2004) Neuronal oscillations in cortical networks

    Science 304:1926–1929.

    https://doi.org/10.1126/science.1099745

    • PubMed
    • Google Scholar
11. 1. Buzsáki G
    2. Wang XJ

    (2012) Mechanisms of gamma oscillations

    Annual Review of Neuroscience 35:203–225.

    https://doi.org/10.1146/annurev-neuro-062111-150444

    • PubMed
    • Google Scholar
12. 1. Canolty RT
    2. Edwards E
    3. Dalal SS
    4. Soltani M
    5. Nagarajan SS
    6. Kirsch HE
    7. Berger MS
    8. Barbaro NM
    9. Knight RT

    (2006) High gamma power is phase-locked to theta oscillations in human neocortex

    Science 313:1626–1628.

    https://doi.org/10.1126/science.1128115

    • PubMed
    • Google Scholar
13. 1. Canolty RT
    2. Knight RT

    (2010) The functional role of cross-frequency coupling

    Trends in Cognitive Sciences 14:506–515.

    https://doi.org/10.1016/j.tics.2010.09.001

    • PubMed
    • Google Scholar
14. 1. Chehelcheraghi M
    2. van Leeuwen C
    3. Steur E
    4. Nakatani C

    (2017) A neural mass model of cross frequency coupling

    PLOS ONE 12:e0173776.

    https://doi.org/10.1371/journal.pone.0173776

    • PubMed
    • Google Scholar
15. 1. Chrobak JJ
    2. Buzsáki G

    (1998) Gamma oscillations in the entorhinal cortex of the freely behaving rat

    The Journal of Neuroscience 18:388–398.

    https://doi.org/10.1523/JNEUROSCI.18-01-00388.1998

    • PubMed
    • Google Scholar
16. 1. Cohen MX

    (2008) Assessing transient cross-frequency coupling in EEG data

    Journal of Neuroscience Methods 168:494–499.

    https://doi.org/10.1016/j.jneumeth.2007.10.012

    • PubMed
    • Google Scholar
17. 1. Cohen MX

    (2017) Multivariate cross-frequency coupling via generalized eigendecomposition

    eLife 6:e21792.

    https://doi.org/10.7554/eLife.21792

    • PubMed
    • Google Scholar
18. 1. Cole SR
    2. Voytek B

    (2017) Brain oscillations and the importance of waveform shape

    Trends in Cognitive Sciences 21:137–149.

    https://doi.org/10.1016/j.tics.2016.12.008

    • PubMed
    • Google Scholar
19. 1. Colgin LL
    2. Denninger T
    3. Fyhn M
    4. Hafting T
    5. Bonnevie T
    6. Jensen O
    7. Moser MB
    8. Moser EI

    (2009) Frequency of gamma oscillations routes flow of information in the Hippocampus

    Nature 462:353–357.

    https://doi.org/10.1038/nature08573

    • PubMed
    • Google Scholar
20. 1. Csicsvari J
    2. Jamieson B
    3. Wise KD
    4. Buzsáki G

    (2003) Mechanisms of gamma oscillations in the Hippocampus of the behaving rat

    Neuron 37:311–322.

    https://doi.org/10.1016/S0896-6273(02)01169-8

    • PubMed
    • Google Scholar
21. 1. Dean HL
    2. Hagan MA
    3. Pesaran B

    (2012) Only coherent spiking in posterior parietal cortex coordinates looking and reaching

    Neuron 73:829–841.

    https://doi.org/10.1016/j.neuron.2011.12.035

    • PubMed
    • Google Scholar
22. 1. Dejean C
    2. Courtin J
    3. Karalis N
    4. Chaudun F
    5. Wurtz H
    6. Bienvenu TC
    7. Herry C

    (2016) Prefrontal neuronal assemblies temporally control fear behaviour

    Nature 535:420–424.

    https://doi.org/10.1038/nature18630

    • PubMed
    • Google Scholar
23. 1. Engel AK
    2. Fries P
    3. Singer W

    (2001) Dynamic predictions: oscillations and synchrony in top-down processing

    Nature Reviews Neuroscience 2:704–716.

    https://doi.org/10.1038/35094565

    • PubMed
    • Google Scholar
24. 1. Fontolan L
    2. Krupa M
    3. Hyafil A
    4. Gutkin B

    (2013) Analytical insights on theta-gamma coupled neural oscillators

    The Journal of Mathematical Neuroscience 3:16.

    https://doi.org/10.1186/2190-8567-3-16

    • PubMed
    • Google Scholar
25. 1. Fries P
    2. Reynolds JH
    3. Rorie AE
    4. Desimone R

    (2001) Modulation of oscillatory neuronal synchronization by selective visual attention

    Science 291:1560–1563.

    https://doi.org/10.1126/science.1055465

    • PubMed
    • Google Scholar
26. 1. Fries P
    2. Nikolić D
    3. Singer W

    (2007) The gamma cycle

    Trends in Neurosciences 30:309–316.

    https://doi.org/10.1016/j.tins.2007.05.005

    • PubMed
    • Google Scholar
27. 1. Gordon JA

    (2016) On being a circuit psychiatrist

    Nature Neuroscience 19:1385–1386.

    https://doi.org/10.1038/nn.4419

    • PubMed
    • Google Scholar
28. 1. Hawellek DJ
    2. Wong YT
    3. Pesaran B

    (2016) Temporal coding of reward-guided choice in the posterior parietal cortex

    PNAS 113:13492–13497.

    https://doi.org/10.1073/pnas.1606479113

    • PubMed
    • Google Scholar
29. 1. Hyafil A
    2. Fontolan L
    3. Kabdebon C
    4. Gutkin B
    5. Giraud AL

    (2015a) Speech encoding by coupled cortical theta and gamma oscillations

    eLife 4:e06213.

    https://doi.org/10.7554/eLife.06213

    • PubMed
    • Google Scholar
30. 1. Hyafil A
    2. Giraud AL
    3. Fontolan L
    4. Gutkin B

    (2015b) Neural Cross-Frequency coupling: connecting architectures, mechanisms, and functions

    Trends in Neurosciences 38:725–740.

    https://doi.org/10.1016/j.tins.2015.09.001

    • PubMed
    • Google Scholar
31. 1. Jackson N
    2. Cole SR
    3. Voytek B
    4. Swann NC

    (2019) Characteristics of waveform shape in Parkinson's Disease Detected with Scalp Electroencephalography

    Eneuro 6:ENEURO.0151-19.2019.

    https://doi.org/10.1523/ENEURO.0151-19.2019

    • PubMed
    • Google Scholar
32. 1. Jensen O
    2. Spaak E
    3. Park H

    (2016) Discriminating valid from spurious indices of Phase-Amplitude coupling

    Eneuro 3:ENEURO.0334-16.2016.

    https://doi.org/10.1523/ENEURO.0334-16.2016

    • PubMed
    • Google Scholar
33. 1. Jensen O
    2. Lisman JE

    (2000) Position reconstruction from an ensemble of hippocampal place cells: contribution of theta phase coding

    Journal of Neurophysiology 83:2602–2609.

    https://doi.org/10.1152/jn.2000.83.5.2602

    • PubMed
    • Google Scholar
34. 1. Jia X
    2. Kohn A

    (2011) Gamma rhythms in the brain

    PLOS Biology 9:e1001045.

    https://doi.org/10.1371/journal.pbio.1001045

    • PubMed
    • Google Scholar
35. 1. Jirsa V
    2. Müller V

    (2013) Cross-frequency coupling in real and virtual brain networks

    Frontiers in Computational Neuroscience 7:78.

    https://doi.org/10.3389/fncom.2013.00078

    • PubMed
    • Google Scholar
36. 1. Jones MW
    2. Wilson MA

    (2005) Theta rhythms coordinate hippocampal–prefrontal interactions in a spatial memory task

    PLOS Biology 11:e402.

    https://doi.org/10.1371/journal.pbio.0030402

    • Google Scholar
37. 1. Karalis N
    2. Dejean C
    3. Chaudun F
    4. Khoder S
    5. Rozeske RR
    6. Wurtz H
    7. Bagur S
    8. Benchenane K
    9. Sirota A
    10. Courtin J
    11. Herry C

    (2016) 4-Hz oscillations synchronize prefrontal-amygdala circuits during fear behavior

    Nature Neuroscience 19:605–612.

    https://doi.org/10.1038/nn.4251

    • PubMed
    • Google Scholar
38. 1. Kopell N
    2. Ermentrout GB
    3. Whittington MA
    4. Traub RD

    (2000) Gamma rhythms and beta rhythms have different synchronization properties

    PNAS 97:1867–1872.

    https://doi.org/10.1073/pnas.97.4.1867

    • PubMed
    • Google Scholar
39. 1. Kramer MA
    2. Tort AB
    3. Kopell NJ

    (2008) Sharp edge artifacts and spurious coupling in EEG frequency comodulation measures

    Journal of Neuroscience Methods 170:352–357.

    https://doi.org/10.1016/j.jneumeth.2008.01.020

    • PubMed
    • Google Scholar
40. 1. Kramer MA
    2. Eden UT

    (2013) Assessment of cross-frequency coupling with confidence using generalized linear models

    Journal of Neuroscience Methods 220:64–74.

    https://doi.org/10.1016/j.jneumeth.2013.08.006

    • PubMed
    • Google Scholar
41. Book

    1. Kramer MA
    2. Eden UT

    (2016)

    Case Studies in Neural Data Analysis: A Guide for the Practicing Neuroscientist

    The MIT Press.

    • Google Scholar
42. 1. Lachaux JP
    2. Rodriguez E
    3. Martinerie J
    4. Varela FJ

    (1999) Measuring phase synchrony in brain signals

    Human Brain Mapping 8:194–208.

    https://doi.org/10.1002/(SICI)1097-0193(1999)8:4<194::AID-HBM4>3.0.CO;2-C

    • PubMed
    • Google Scholar
43. 1. Lakatos P
    2. Shah AS
    3. Knuth KH
    4. Ulbert I
    5. Karmos G
    6. Schroeder CE

    (2005) An oscillatory hierarchy controlling neuronal excitability and stimulus processing in the auditory cortex

    Journal of Neurophysiology 94:1904–1911.

    https://doi.org/10.1152/jn.00263.2005

    • PubMed
    • Google Scholar
44. 1. Lakatos P
    2. Karmos G
    3. Mehta AD
    4. Ulbert I
    5. Schroeder CE

    (2008) Entrainment of neuronal oscillations as a mechanism of attentional selection

    Science 320:110–113.

    https://doi.org/10.1126/science.1154735

    • PubMed
    • Google Scholar
45. 1. Lepage KQ
    2. Vijayan S

    (2015) A Time-Series model of phase amplitude cross frequency coupling and comparison of spectral characteristics with neural data

    BioMed Research International 2015:1–8.

    https://doi.org/10.1155/2015/140837

    • Google Scholar
46. 1. Likhtik E
    2. Stujenske JM
    3. Topiwala MA
    4. Harris AZ
    5. Gordon JA

    (2014) Prefrontal entrainment of amygdala activity signals safety in learned fear and innate anxiety

    Nature Neuroscience 17:106–113.

    https://doi.org/10.1038/nn.3582

    • PubMed
    • Google Scholar
47. 1. Lisman J

    (2005) The theta/gamma discrete phase code occuring during the hippocampal phase precession may be a more general brain coding scheme

    Hippocampus 15:913–922.

    https://doi.org/10.1002/hipo.20121

    • PubMed
    • Google Scholar
48. 1. Malerba P
    2. Kopell N

    (2013) Phase resetting reduces theta-gamma rhythmic interaction to a one-dimensional map

    Journal of Mathematical Biology 66:1361–1386.

    https://doi.org/10.1007/s00285-012-0534-9

    • PubMed
    • Google Scholar
49. 1. Mann EO
    2. Mody I

    (2010) Control of hippocampal gamma oscillation frequency by tonic inhibition and excitation of interneurons

    Nature Neuroscience 13:205–212.

    https://doi.org/10.1038/nn.2464

    • PubMed
    • Google Scholar
50. 1. Mathalon DH
    2. Sohal VS

    (2015) Neural oscillations and synchrony in brain dysfunction and neuropsychiatric disorders: it's about time

    JAMA Psychiatry 72:840.

    https://doi.org/10.1001/jamapsychiatry.2015.0483

    • PubMed
    • Google Scholar
51. 1. Mazzoni A
    2. Whittingstall K
    3. Brunel N
    4. Logothetis NK
    5. Panzeri S

    (2010) Understanding the relationships between spike rate and Delta/gamma frequency bands of LFPs and EEGs using a local cortical network model

    NeuroImage 52:956–972.

    https://doi.org/10.1016/j.neuroimage.2009.12.040

    • PubMed
    • Google Scholar
52. 1. Mormann F
    2. Fell J
    3. Axmacher N
    4. Weber B
    5. Lehnertz K
    6. Elger CE
    7. Fernández G

    (2005) Phase/amplitude reset and theta–gamma interaction in the human medial temporal lobe during a continuous word recognition memory task

    Hippocampus 15:890–900.

    https://doi.org/10.1002/hipo.20117

    • PubMed
    • Google Scholar
53. 1. Mukamel R
    2. Gelbard H
    3. Arieli A
    4. Hasson U
    5. Fried I
    6. Malach R

    (2005) Coupling between neuronal firing, field potentials, and FMRI in human auditory cortex

    Science 309:951–954.

    https://doi.org/10.1126/science.1110913

    • PubMed
    • Google Scholar
54. Software

    1. Nadalin J
    2. Kramer M

    (2019) A GLM based approach to measure CFC, version 0649e2f

    GitHub.

    https://github.com/Eden-Kramer-Lab/GLM-CFC
55. 1. Onslow AC
    2. Bogacz R
    3. Jones MW

    (2011) Quantifying phase-amplitude coupling in neuronal network oscillations

    Progress in Biophysics and Molecular Biology 105:49–57.

    https://doi.org/10.1016/j.pbiomolbio.2010.09.007

    • PubMed
    • Google Scholar
56. 1. Onslow AC
    2. Jones MW
    3. Bogacz R

    (2014) A canonical circuit for generating phase-amplitude coupling

    PLOS ONE 9:e102591.

    https://doi.org/10.1371/journal.pone.0102591

    • PubMed
    • Google Scholar
57. 1. Osipova D
    2. Hermes D
    3. Jensen O

    (2008) Gamma power is phase-locked to posterior alpha activity

    PLOS ONE 3:e3990.

    https://doi.org/10.1371/journal.pone.0003990

    • PubMed
    • Google Scholar
58. 1. Penny WD
    2. Duzel E
    3. Miller KJ
    4. Ojemann JG

    (2008) Testing for nested oscillation

    Journal of Neuroscience Methods 174:50–61.

    https://doi.org/10.1016/j.jneumeth.2008.06.035

    • PubMed
    • Google Scholar
59. 1. Pesaran B
    2. Pezaris JS
    3. Sahani M
    4. Mitra PP
    5. Andersen RA

    (2002) Temporal structure in neuronal activity during working memory in macaque parietal cortex

    Nature Neuroscience 5:805–811.

    https://doi.org/10.1038/nn890

    • PubMed
    • Google Scholar
60. 1. Pesaran B
    2. Nelson MJ
    3. Andersen RA

    (2008) Free choice activates a decision circuit between frontal and parietal cortex

    Nature 453:406–409.

    https://doi.org/10.1038/nature06849

    • PubMed
    • Google Scholar
61. 1. Popov T
    2. Jensen O
    3. Schoffelen JM

    (2018) Dorsal and ventral cortices are coupled by cross-frequency interactions during working memory

    NeuroImage 178:277–286.

    https://doi.org/10.1016/j.neuroimage.2018.05.054

    • PubMed
    • Google Scholar
62. 1. Rasch MJ
    2. Gretton A
    3. Murayama Y
    4. Maass W
    5. Logothetis NK

    (2008) Inferring spike trains from local field potentials

    Journal of Neurophysiology 99:1461–1476.

    https://doi.org/10.1152/jn.00919.2007

    • PubMed
    • Google Scholar
63. 1. Ray S
    2. Crone NE
    3. Niebur E
    4. Franaszczuk PJ
    5. Hsiao SS

    (2008a) Neural correlates of high-gamma oscillations (60-200 hz) in macaque local field potentials and their potential implications in electrocorticography

    Journal of Neuroscience 28:11526–11536.

    https://doi.org/10.1523/JNEUROSCI.2848-08.2008

    • PubMed
    • Google Scholar
64. 1. Ray S
    2. Hsiao SS
    3. Crone NE
    4. Franaszczuk PJ
    5. Niebur E

    (2008b) Effect of stimulus intensity on the spike-local field potential relationship in the secondary somatosensory cortex

    Journal of Neuroscience 28:7334–7343.

    https://doi.org/10.1523/JNEUROSCI.1588-08.2008

    • PubMed
    • Google Scholar
65. 1. Ray S
    2. Maunsell JH

    (2011) Different origins of gamma rhythm and high-gamma activity in macaque visual cortex

    PLOS Biology 9:e1000610.

    https://doi.org/10.1371/journal.pbio.1000610

    • PubMed
    • Google Scholar
66. 1. Sase T
    2. Katori Y
    3. Komuro M
    4. Aihara K

    (2017) Bifurcation analysis on Phase-Amplitude Cross-Frequency coupling in neural networks with dynamic synapses

    Frontiers in Computational Neuroscience 11:18.

    https://doi.org/10.3389/fncom.2017.00018

    • PubMed
    • Google Scholar
67. 1. Scheffer-Teixeira R
    2. Belchior H
    3. Leão RN
    4. Ribeiro S
    5. Tort AB

    (2013) On high-frequency field oscillations (>100 hz) and the spectral leakage of spiking activity

    Journal of Neuroscience 33:1535–1539.

    https://doi.org/10.1523/JNEUROSCI.4217-12.2013

    • PubMed
    • Google Scholar
68. 1. Shirvalkar PR
    2. Rapp PR
    3. Shapiro ML

    (2010) Bidirectional changes to hippocampal theta-gamma comodulation predict memory for recent spatial episodes

    PNAS 107:7054–7059.

    https://doi.org/10.1073/pnas.0911184107

    • PubMed
    • Google Scholar
69. 1. Siegel M
    2. Warden MR
    3. Miller EK

    (2009) Phase-dependent neuronal coding of objects in short-term memory

    PNAS 106:21341–21346.

    https://doi.org/10.1073/pnas.0908193106

    • PubMed
    • Google Scholar
70. 1. Sirota A
    2. Montgomery S
    3. Fujisawa S
    4. Isomura Y
    5. Zugaro M
    6. Buzsáki G

    (2008) Entrainment of neocortical neurons and gamma oscillations by the hippocampal theta rhythm

    Neuron 60:683–697.

    https://doi.org/10.1016/j.neuron.2008.09.014

    • PubMed
    • Google Scholar
71. 1. Sotero RC

    (2016) Topology, Cross-Frequency, and Same-Frequency band interactions shape the generation of Phase-Amplitude coupling in a neural mass model of a cortical column

    PLOS Computational Biology 12:e1005180.

    https://doi.org/10.1371/journal.pcbi.1005180

    • PubMed
    • Google Scholar
72. 1. Spaak E
    2. Bonnefond M
    3. Maier A
    4. Leopold DA
    5. Jensen O

    (2012) Layer-specific entrainment of γ-band neural activity by the α rhythm in monkey visual cortex

    Current Biology 22:2313–2318.

    https://doi.org/10.1016/j.cub.2012.10.020

    • PubMed
    • Google Scholar
73. 1. Theiler J
    2. Eubank S
    3. Longtin A
    4. Galdrikian B
    5. Doyne Farmer J

    (1992) Testing for nonlinearity in time series: the method of surrogate data

    Physica D: Nonlinear Phenomena 58:77–94.

    https://doi.org/10.1016/0167-2789(92)90102-S

    • Google Scholar
74. 1. Tort AB
    2. Rotstein HG
    3. Dugladze T
    4. Gloveli T
    5. Kopell NJ

    (2007) On the formation of gamma-coherent cell assemblies by Oriens lacunosum-moleculare interneurons in the Hippocampus

    PNAS 104:13490–13495.

    https://doi.org/10.1073/pnas.0705708104

    • PubMed
    • Google Scholar
75. 1. Tort AB
    2. Kramer MA
    3. Thorn C
    4. Gibson DJ
    5. Kubota Y
    6. Graybiel AM
    7. Kopell NJ

    (2008) Dynamic cross-frequency couplings of local field potential oscillations in rat striatum and Hippocampus during performance of a T-maze task

    PNAS 105:20517–20522.

    https://doi.org/10.1073/pnas.0810524105

    • PubMed
    • Google Scholar
76. 1. Tort AB
    2. Komorowski RW
    3. Manns JR
    4. Kopell NJ
    5. Eichenbaum H

    (2009) Theta-gamma coupling increases during the learning of item-context associations

    PNAS 106:20942–20947.

    https://doi.org/10.1073/pnas.0911331106

    • PubMed
    • Google Scholar
77. 1. Tort AB
    2. Komorowski R
    3. Eichenbaum H
    4. Kopell N

    (2010) Measuring phase-amplitude coupling between neuronal oscillations of different frequencies

    Journal of Neurophysiology 104:1195–1210.

    https://doi.org/10.1152/jn.00106.2010

    • PubMed
    • Google Scholar
78. 1. Tort ABL
    2. Brankačk J
    3. Draguhn A

    (2018) Respiration-Entrained brain rhythms are global but often overlooked

    Trends in Neurosciences 41:186–197.

    https://doi.org/10.1016/j.tins.2018.01.007

    • PubMed
    • Google Scholar
79. 1. van Wijk BC
    2. Jha A
    3. Penny W
    4. Litvak V

    (2015) Parametric estimation of cross-frequency coupling

    Journal of Neuroscience Methods 243:94–102.

    https://doi.org/10.1016/j.jneumeth.2015.01.032

    • PubMed
    • Google Scholar
80. 1. Vanhatalo S
    2. Palva JM
    3. Holmes MD
    4. Miller JW
    5. Voipio J
    6. Kaila K

    (2004) Infraslow oscillations modulate excitability and interictal epileptic activity in the human cortex during sleep

    PNAS 101:5053–5057.

    https://doi.org/10.1073/pnas.0305375101

    • PubMed
    • Google Scholar
81. 1. Varela F
    2. Lachaux JP
    3. Rodriguez E
    4. Martinerie J

    (2001) The brainweb: phase synchronization and large-scale integration

    Nature Reviews Neuroscience 2:229–239.

    https://doi.org/10.1038/35067550

    • PubMed
    • Google Scholar
82. 1. Vinck M
    2. Lima B
    3. Womelsdorf T
    4. Oostenveld R
    5. Singer W
    6. Neuenschwander S
    7. Fries P

    (2010) Gamma-phase shifting in awake monkey visual cortex

    Journal of Neuroscience 30:1250–1257.

    https://doi.org/10.1523/JNEUROSCI.1623-09.2010

    • PubMed
    • Google Scholar
83. 1. von Stein A
    2. Sarnthein J

    (2000) Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronization

    International Journal of Psychophysiology 38:301–313.

    https://doi.org/10.1016/S0167-8760(00)00172-0

    • PubMed
    • Google Scholar
84. 1. Voytek B
    2. D'Esposito M
    3. Crone N
    4. Knight RT

    (2013) A method for event-related phase/amplitude coupling

    NeuroImage 64:416–424.

    https://doi.org/10.1016/j.neuroimage.2012.09.023

    • PubMed
    • Google Scholar
85. 1. Voytek B
    2. Knight RT

    (2015) Dynamic network communication as a unifying neural basis for cognition, development, aging, and disease

    Biological Psychiatry 77:1089–1097.

    https://doi.org/10.1016/j.biopsych.2015.04.016

    • PubMed
    • Google Scholar
86. 1. Wagner FB
    2. Eskandar EN
    3. Cosgrove GR
    4. Madsen JR
    5. Blum AS
    6. Potter NS
    7. Hochberg LR
    8. Cash SS
    9. Truccolo W

    (2015) Microscale spatiotemporal dynamics during neocortical propagation of human focal seizures

    NeuroImage 122:114–130.

    https://doi.org/10.1016/j.neuroimage.2015.08.019

    • PubMed
    • Google Scholar
87. 1. Weiss SA
    2. Lemesiou A
    3. Connors R
    4. Banks GP
    5. McKhann GM
    6. Goodman RR
    7. Zhao B
    8. Filippi CG
    9. Nowell M
    10. Rodionov R
    11. Diehl B
    12. McEvoy AW
    13. Walker MC
    14. Trevelyan AJ
    15. Bateman LM
    16. Emerson RG
    17. Schevon CA

    (2015) Seizure localization using ictal phase-locked high gamma

    Neurology 84:2320–2328.

    https://doi.org/10.1212/WNL.0000000000001656

    • Google Scholar
88. 1. Whittingstall K
    2. Logothetis NK

    (2009) Frequency-band coupling in surface EEG reflects spiking activity in monkey visual cortex

    Neuron 64:281–289.

    https://doi.org/10.1016/j.neuron.2009.08.016

    • PubMed
    • Google Scholar
89. 1. Whittington MA
    2. Traub RD
    3. Kopell N
    4. Ermentrout B
    5. Buhl EH

    (2000) Inhibition-based rhythms: experimental and mathematical observations on network dynamics

    International Journal of Psychophysiology 38:315–336.

    https://doi.org/10.1016/S0167-8760(00)00173-2

    • PubMed
    • Google Scholar
90. 1. Whittington MA
    2. Cunningham MO
    3. LeBeau FE
    4. Racca C
    5. Traub RD

    (2011) Multiple origins of the cortical γ rhythm

    Developmental Neurobiology 71:92–106.

    https://doi.org/10.1002/dneu.20814

    • PubMed
    • Google Scholar
91. 1. Widge AS
    2. Ellard KK
    3. Paulk AC
    4. Basu I
    5. Yousefi A
    6. Zorowitz S
    7. Gilmour A
    8. Afzal A
    9. Deckersbach T
    10. Cash SS
    11. Kramer MA
    12. Eden UT
    13. Dougherty DD
    14. Eskandar EN

    (2017) Treating refractory mental illness with closed-loop brain stimulation: progress towards a patient-specific transdiagnostic approach

    Experimental Neurology 287:461–472.

    https://doi.org/10.1016/j.expneurol.2016.07.021

    • PubMed
    • Google Scholar
92. 1. Wong YT
    2. Fabiszak MM
    3. Novikov Y
    4. Daw ND
    5. Pesaran B

    (2016) Coherent neuronal ensembles are rapidly recruited when making a look-reach decision

    Nature Neuroscience 19:327–334.

    https://doi.org/10.1038/nn.4210

    • PubMed
    • Google Scholar
93. 1. Wulff P
    2. Ponomarenko AA
    3. Bartos M
    4. Korotkova TM
    5. Fuchs EC
    6. Bähner F
    7. Both M
    8. Tort AB
    9. Kopell NJ
    10. Wisden W
    11. Monyer H

    (2009) Hippocampal theta rhythm and its coupling with gamma oscillations require fast inhibition onto parvalbumin-positive interneurons

    PNAS 106:3561–3566.

    https://doi.org/10.1073/pnas.0813176106

    • PubMed
    • Google Scholar
94. 1. Yousefi A
    2. Basu I
    3. Paulk AC
    4. Peled N
    5. Eskandar EN
    6. Dougherty DD
    7. Cash SS
    8. Widge AS
    9. Eden UT

    (2019) Decoding hidden cognitive states from behavior and physiology using a bayesian approach

    Neural Computation 31:1751–1788.

    https://doi.org/10.1162/neco_a_01196

    • PubMed
    • Google Scholar
95. 1. Zhong W
    2. Ciatipis M
    3. Wolfenstetter T
    4. Jessberger J
    5. Müller C
    6. Ponsel S
    7. Yanovsky Y
    8. Brankačk J
    9. Tort ABL
    10. Draguhn A

    (2017) Selective entrainment of gamma subbands by different slow network oscillations

    PNAS 114:4519–4524.

    https://doi.org/10.1073/pnas.1617249114

    • PubMed
    • Google Scholar