A statistical framework to assess crossfrequency coupling while accounting for confounding analysis effects
Abstract
Cross frequency coupling (CFC) is emerging as a fundamental feature of brain activity, correlated with brain function and dysfunction. Many different types of CFC have been identified through application of numerous data analysis methods, each developed to characterize a specific CFC type. Choosing an inappropriate method weakens statistical power and introduces opportunities for confounding effects. To address this, we propose a statistical modeling framework to estimate high frequency amplitude as a function of both the low frequency amplitude and low frequency phase; the result is a measure of phaseamplitude coupling that accounts for changes in the low frequency amplitude. We show in simulations that the proposed method successfully detects CFC between the low frequency phase or amplitude and the high frequency amplitude, and outperforms an existing method in biologicallymotivated examples. Applying the method to in vivo data, we illustrate examples of CFC during a seizure and in response to electrical stimuli.
https://doi.org/10.7554/eLife.44287.001Introduction
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 highgamma (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 crossfrequency 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: phasephase coupling (PPC), phaseamplitude coupling (PAC), and amplitudeamplitude 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), phaselocking value (Mormann et al., 2005; Lachaux et al., 1999; Vanhatalo et al., 2004), envelopetosignal 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 phaseamplitude 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 highfrequency 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 lowfrequency 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.
Materials and methods
Estimation of the phase and amplitude envelope
View detailed protocolTo study CFC we estimate three quantities: the phase of the low frequency signal, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$; the amplitude envelope of the high frequency signal, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$; and the amplitude envelope of the low frequency signal, $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$. To do so, we first bandpass filter the data into low frequency (4–7 Hz) and high frequency (100–140 Hz) signals, $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, respectively, using a leastsquares linearphase FIR filter of order 375 for the high frequency signal, and order 50 for the low frequency signal. Here we choose specific high and low frequency ranges of interest, motivated by previous in vivo observations (Canolty et al., 2006; Tort et al., 2008; SchefferTeixeira et al., 2013). However, we note that this method is flexible and not dependent on this choice. We select a wide high frequency band consistent with recommendations from the literature (Aru et al., 2015) and the mechanistic explanation that extracellular spikes produce this broadband high frequency activity (SchefferTeixeira et al., 2013). We use the Hilbert transform to compute the analytic signals of $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, and from these compute the phase and amplitude of the low frequency signal $({A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and ${\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}})$ and the amplitude of the high frequency signal $({A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}})$.
Modeling framework to assess CFC
Generalized linear models (GLMs) provide a principled framework to assess CFC (Penny et al., 2008; Kramer and Eden, 2013; van Wijk et al., 2015). Here, we present three models to analyze different types of CFC. The fundamental logic behind this approach is to model the distribution of $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ as a function of different predictors. In existing measures of PAC, the distribution of $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ versus $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is assessed using a variety of different metrics (e.g., Tort et al., 2010). Here, we estimate statistical models to fit $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ as a function of $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, and their combinations. If these models fit the data sufficiently well, then we estimate distances between the modeled surfaces to measure the impact of each predictor.
The $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model
Request a detailed protocolThe $\varphi}_{low$ model relates $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, the response variable, to a linear combination of $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, the predictor variable, expressed in a spline basis:
where the conditional distribution of $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ given $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is modeled as a Gamma random variable with mean parameter $\mu $ and shape parameter $\nu $, and ${\beta}_{k}$ are undetermined coefficients, which we refer to collectively as $\beta}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$. We choose this distribution as it guarantees real, positive amplitude values; we note that this distribution provides an acceptable fit to the example human data analyzed here (Figure 1). The functions $\{{f}_{1},\mathrm{\cdots},{f}_{n}\}$ correspond to spline basis functions, with $n$ control points equally spaced between 0 and $2\pi $, used to approximate $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$. We note that the spline functions sum to 1, and therefore we omit a constant offset term. We use a tension parameter of 0.5, which controls the smoothness of the splines. We note that, because the link function of the conditional mean of the response variable $({A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}})$ varies linearly with the model coefficients ${\beta}_{k}$ the model is a GLM, though the spline basis functions situate the model in the larger class of Generalized Additive Models (GAMs). Here we fix $n=10$, which is a reasonable choice for smooth PAC with one or two broad peaks (Kramer and Eden, 2013). To support this choice, we apply an AICbased selection procedure to 1000 simulated instances of signals of duration 20 s with phaseamplitude coupling and amplitudeamplitude coupling (see Materials and methods: Synthetic Time Series with PAC and Synthetic Time Series with AAC, below, for simulation details). For each simulation, we fit the model in Equation 1 to these data for 27 different values of $n$ from $n=4$ to $n=30$. For each simulated signal, we record the value of $n$ that minimizes the AIC, defined as
where $\mathrm{\Delta}$ is the deviance from the model in Equation 1. The values of $n$ that minimize the AIC tend to lie between $n=7$ and $n=12$ (Figure 2). These simulations support the choice of $n=10$ as a sufficient number of splines.
For a more detailed discussion and simulation examples of the PAC model, see Kramer and Eden (2013). We note that the choices of distribution and link function differ from those in Penny et al. (2008) and van Wijk et al. (2015), where the normal distribution and identity link are used instead.
The $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model
Request a detailed protocolThe $A}_{low$ model relates the high frequency amplitude to the low frequency amplitude:
where the conditional distribution of $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ given $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is modeled as a Gamma random variable with mean parameter $\mu $ and shape parameter $\nu $. The predictor consists of a single variable and a constant, and the length of the coefficient vector $\beta}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}}=\{{\beta}_{1},{\beta}_{2}\$ is 2.
The $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model
Request a detailed protocolThe $A}_{low},{\varphi}_{low$ model extends the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model in Equation 1 by including three additional predictors in the GLM: $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, the low frequency amplitude; and interaction terms between the low frequency amplitude and the low frequency phase: ${A}_{\mathrm{l}\mathrm{o}\mathrm{w}}\mathrm{sin}({\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}})$, and ${A}_{\mathrm{l}\mathrm{o}\mathrm{w}}\mathrm{cos}({\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}})$. These new terms allow assessment of phaseamplitude coupling while accounting for linear amplitudeamplitude dependence and more complicated phasedependent relationships on the low frequency amplitude without introducing many more parameters. Compared to the original $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model in Equation 1, including these new terms increases the number of variables to $n+3$, and the length of the coefficient vector $\beta}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to $n+3$. These changes result in the following model:
Here, the conditional distribution of $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ given $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is modeled as a Gamma random variable with mean parameter $\mu $ and shape parameter $\nu $, and ${\beta}_{k}$ are undetermined coefficients. We note that we only consider two interaction terms, rather than the spline basis function of phase, to limit the number of parameters in the model.
The statistics $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$
Request a detailed protocolWe compute two measures of CFC, $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$ which use the three models defined in the previous section. We evaluate each model in the threedimensional space ($\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) and calculate the statistics $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$. We use the MATLAB $(RRID:SC{R}_{0}01622)$ function fitglm to estimate the models; we note that this procedure estimates the dispersion directly for the gamma distribution. In what follows, we first discuss the three model surfaces estimated from the data, and then how we use these surfaces to compute the statistics $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$.
To create the surface $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$, which fits the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model in the threedimensional ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) space, we first compute estimates of the parameters $\beta}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ in Equation 3. We then estimate $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ by fixing $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ at one of 640 evenly spaced values between the 5th and 95th quantiles of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ observed; we choose these quantiles to avoid extremely small or large values of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$. Finally, at the fixed $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, we compute the high frequency amplitude values from the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model over 100 evenly spaced values of $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ between $\pi $ and $\pi $. This results in a twodimensional curve $C}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ in the twodimensional ($\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, ${A}_{\text{high}}$) space with fixed ${A}_{\text{low}}$. We repeat this procedure for all 640 values of ${A}_{\text{low}}$ to create a surface $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ in the threedimensional space ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) (Figure 3C). To create the surface $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$, which fits the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model in the threedimensional ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) space, we estimate the coefficient vector $\beta}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ for the model in Equation 2. We then estimate the high frequency amplitude over 640 evenly spaced values between the 5th and 95th quantiles of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ observed, again to avoid extremely small or large values of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$. This creates a mean response function which appears as a curve $C}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ in the twodimensional ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) space. We extend this twodimensional curve to a threedimensional surface $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ by extending $C}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ along the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ dimension (Figure 3A).
To create the surface $S}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$, which fits the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model in the threedimensional ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) space, we first estimate the coefficients $\beta}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ for the model in Equation 1. From this, we then compute estimates for the high frequency amplitude using the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model with 100 evenly spaced values of $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ between $\pi $ and $\pi $. This results in the mean response function of the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model. We extend this curve $C}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ in the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ dimension to create a surface $S}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ in the threedimensional ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) space. The surface $S}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ has the same structure as the curve $C}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ in the ($\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) space, and remains constant along the dimension $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ (Figure 3B).
The statistic $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ measures the effect of low frequency phase on high frequency amplitude, while accounting for fluctuations in the low frequency amplitude. To compute this statistic, we note that the model in Equation 3 measures the combined effect of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ on $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, while the model in Equation 2 measures only the effect of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ on $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$. Hence, to isolate the effect of $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ on $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, while accounting for $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, we compare the difference in fits between the models in Equations 2 and 3. We fit the mean response functions of the models in Equations 2 and 3, and calculate $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ as the maximum absolute fractional difference between the resulting surfaces $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ (Figure 3D):
That is we measure the largest distance between the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ models. We expect fluctuations in $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ not present in $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to be the result of $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, that is PAC. In the absence of PAC, we expect the surfaces $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to be very close, resulting in a small value of $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$. However, in the presence of PAC, we expect $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to deviate from $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$, resulting in a large value of $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$. We note that this measure, unlike R_{2} metrics for linear regression, is not meant to measure the goodnessoffit of these models to the data, but rather the differences in fits between the two models. We also note that $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ is an unbounded measure, as it equals the maximum absolute fractional difference between distributions, which may exceed 1.
To compute the statistic $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$, which measures the effect of low frequency amplitude on high frequency amplitude while accounting for fluctuations in the low frequency phase, we compare the difference in fits of the model in Equation 3 from the model in Equation 1. We note that the model in Equation 3 predicts $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ as a function of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, while the model in Equation 1 predicts $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ as a function of $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ only. Therefore we expect a difference in fits between the models in Equations 1 and 3 results from the effects of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ on $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$. We fit the mean response functions of the models in Equations 1 and 3 in the threedimensional ($\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) space, and calculate $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$ as the maximum absolute fractional difference between the resulting surfaces $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and $S}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ (Figure 3E):
That is we measure the distance between the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ models. We expect fluctuations in $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ not present in $S}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to be the result of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, that is AAC. In the absence of AAC, we expect the surfaces $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and $S}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to be very close, resulting in a small value for $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$. Alternatively, in the presence of AAC, we expect $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to deviate from $S}_{{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$, resulting in a large value of $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$.
Estimating 95% confidence intervals for $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$
Request a detailed protocolWe compute 95% confidence intervals for $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$ via a parametric bootstrap method (Kramer and Eden, 2013). Given a vector of estimated coefficients ${\beta}_{\text{x}}$ for $x=\{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}};\text{}{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}};\mathrm{o}\mathrm{r}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{A}}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}\}$, we use its estimated covariance and estimated mean to generate 10,000 normally distributed coefficient sample vectors ${\beta}_{\text{x}}^{j}$, $j\in \{0,\mathrm{\dots},10000\}$. For each ${\beta}_{\text{x}}^{j}$, we then compute the high frequency amplitude values from the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, or $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, ${S}_{\text{x}}^{j}$. Finally, we compute the statistics $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}^{j$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}}^{j$ for each $j$ as,
The 95% confidence intervals for the statistics are the values of $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}^{j$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}}^{j$ at the 0.025 and 0.975 quantiles (Kramer and Eden, 2013).
Assessing significance of AAC and PAC with bootstrap pvalues
Request a detailed protocolTo assess whether evidence exists for significant PAC or AAC, we implement a bootstrap procedure to compute pvalues as follows. Given two signals $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, and the resulting estimated statistics $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$ we apply the Amplitude Adjusted Fourier Transform (AAFT) algorithm (Theiler et al., 1992) on $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ to generate a surrogate signal $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}^{i$. In the AAFT algorithm, we first reorder the values of ${V}_{\text{high}}$ by creating a random Gaussian signal $W$ and ordering the values of $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ to match $W$. For example, if the highest value of $W$ occurs at index $j$, then the highest value of $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ will be reordered to occur at index $j$. Next, we apply the Fourier Transform (FT) to the reordered $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ and randomize the phase of the frequency domain signal. This signal is then inverse Fourier transformed and rescaled to have the same amplitude distribution as the original signal $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$. In this way, the algorithm produces a permutation $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}^{i$ of $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ such that the power spectrum and amplitude distribution of the original signal are preserved.
We create 1000 such surrogate signals $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}^{i$, and calculate $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}^{i$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}}^{i$ between $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and each $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}^{i$. We define the pvalues ${p}_{\text{PAC}}$ and ${p}_{\text{AAC}}$ as the proportion of values in ${\{{\text{\mathbf{R}}}_{\text{PAC}}^{i}\}}_{i=1}^{1000}$ and ${\{{\text{\mathbf{R}}}_{\text{AAC}}^{i}\}}_{i=1}^{1000}$ greater than the estimated statistics $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$, respectively. If the proportion is zero, we set $p=0.0005$.
We calculate pvalues for the modulation index in the same way. The modulation index calculates the distribution of high frequency amplitudes versus low frequency phases and measures the distance from this distribution to a uniform distribution of amplitudes. Given the signals $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, and the resulting modulation index MI between them, we calculate the modulation index between $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and 1000 surrogate permutations of $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ using the AAFT algorithm. We set ${p}_{\text{MI}}$ to be the proportion of these resulting values greater than the MI value estimated from the original signals.
Synthetic time series with PAC
Request a detailed protocolWe construct synthetic time series to examine the performance of the proposed method as follows. First, we simulate 20 s of pink noise data such that the power spectrum scales as $1/f$. We then filter these data into low (4–7 Hz) and high (100–140 Hz) frequency bands, as described in Materials and methods: Estimation of the phase and amplitude envelope, creating signals ${V}_{\text{low}}$ and ${V}_{\text{high}}$. Next, we couple the amplitude of the high frequency signal to the phase of the low frequency signal. To do so, we first locate the peaks of ${V}_{\text{low}}$ and determine the times ${t}_{k},k=\{1,2,3,\mathrm{\dots},K\}$, of the $K$ relative extrema. We note that these times correspond approximately to ${\varphi}_{\text{low}}=0$. We then create a smooth modulation signal M which consists of a 42 ms Hanning window of height $1+{I}_{\text{PAC}}$ centered at each ${t}_{k}$, and a value of 1 at all other times (Figure 4A). The intensity parameter ${I}_{\text{PAC}}$ in the modulation signal corresponds to the strength of PAC. ${I}_{\text{PAC}}=0.0$ corresponds to no PAC, while ${I}_{\text{PAC}}=1.0$ results in a 100% increase in the high frequency amplitude at each ${t}_{k}$, creating strong PAC. We create a new signal ${V}_{\text{high}}^{\prime}$ with the same phase as ${V}_{\text{high}}$, but with amplitude dependent on the phase of ${V}_{\text{low}}$ by setting,
We create the final voltage trace $V$ as
where ${V}_{\text{pink}}$ is a new instance of pink noise multiplied by a small constant $c=0.01$. In the signal $V$, brief increases of the high frequency activity occur at a specific phase (0 radians) of the low frequency signal (Figure 4B).
Synthetic time series with AAC
Request a detailed protocolTo generate synthetic time series with dependence on the low frequency amplitude, we follow the procedure in the preceding section to generate ${V}_{\text{low}}$, ${V}_{\text{high}}$, and ${A}_{\text{low}}$. We then induce amplitudeamplitude coupling between the low and high frequency components by creating a new signal ${V}_{\text{high}}^{*}$ such that
where $I}_{\mathrm{A}\mathrm{A}\mathrm{C}$ is the intensity parameter corresponding to the strength of amplitudeamplitude coupling. We define the final voltage trace $V$ as
where ${V}_{\text{pink}}$ is a new instance of pink noise multiplied by a small constant $c=0.01$ (Figure 4C).
Human subject data
Request a detailed protocolA patient (male, age 32 years) with medically intractable focal epilepsy underwent clinically indicated intracranial cortical recordings for epilepsy monitoring. In addition to clinical electrode implantation, the patient was also implanted with a 10 × 10 (4 mm ×4 mm) NeuroPort microelectrode array (MEA; Blackrock Microsystems, Utah) in a neocortical area expected to be resected with high probability in the temporal gyrus. The MEA consists of 96 platinumtipped silicon probes, with a length of 1.5 mm, corresponding to neocortical layer III as confirmed by histology after resection. Signals from the MEA were acquired continuously at 30 kHz per channel. Seizure onset times were determined by an experienced encephalographer (S.S.C.) through inspection of the macroelectrode recordings, referral to the clinical report, and clinical manifestations recorded on video. For a detailed clinical summary, see patient P2 of Wagner et al. (2015). For these data, we analyze the 100–140 Hz and 4–7 Hz frequency bands to illustrate the proposed method; a more rigorous study of CFC in these data may require a more principled choice of high frequency band. All patients were enrolled after informed consent and consent to publish was obtained, and approval was granted by local Institutional Review Boards at Massachusetts General Hospital and Brigham Women’s Hospitals (Partners Human Research Committee), and at Boston University according to National Institutes of Health guidelines.
Code availability
Request a detailed protocolThe code to perform this analysis is available for reuse and further development at https://github.com/EdenKramerLab/GLMCFC (Nadalin and Kramer, 2019; copy archived at https://github.com/elifesciencespublications/GLMCFC).
Results
We first examine the performance of the CFC measure through simulation examples. In doing so, we show that the statistics ${\text{\mathbf{R}}}_{\text{PAC}}$ and ${\text{\mathbf{R}}}_{\text{AAC}}$ accurately detect different types of crossfrequency 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 $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, the ${A}_{\text{low}}$ model, and the $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.
To assess the distribution of significant R values in the case of no crossfrequency 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 chisquared 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 chisquared 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 AICbased 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{\mathbf{R}}}_{\text{PAC}}$ and ${\text{\mathbf{R}}}_{\text{AAC}}$ implemented here. We conclude that, in this modeling regime, two deviancebased model comparison procedures for GLMs are less robust measures of significant PAC and AAC.
The proposed method accurately detects PAC
We next consider signals that possess phaseamplitude coupling, but lack amplitudeamplitude coupling. To do so, we simulate a 20 s signal with ${A}_{\text{high}}$ modulated by ${\varphi}_{\text{low}}$ (Figure 5D); more specifically, ${A}_{\text{high}}$ increases when ${\varphi}_{\text{low}}$ is near 0 radians (see Materials and methods, ${I}_{\mathrm{P}\mathrm{A}\mathrm{C}}=1$). We then estimate the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, and the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model from these data; we show example fits in Figure 5E. We find that in the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model $A}_{\text{high}$ is higher when $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is close to 0 radians, and the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model follows this trend. We note that, because the data do not depend on the low frequency amplitude (${A}_{\text{low}})$, the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ models have very similar shapes in the (${\varphi}_{\text{low}}$, ${A}_{\text{low}}$, ${A}_{\text{high}}$) space, and the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model is nearly flat.
Simulating 1000 instances of these 20 s signals with induced phaseamplitude coupling, we find ${p}_{\text{AAC}}<0.05$ for only 0.6% of the simulations, while ${p}_{\text{PAC}}<0.05$ for 96.5% of the simulations. We find that the significant values of ${\text{\mathbf{R}}}_{\text{PAC}}$ lie well above 0 (Figure 5F), and that as the intensity of the simulated phaseamplitude coupling increases, so does the statistic ${\text{\mathbf{R}}}_{\text{PAC}}$ (Figure 5G). We conclude that the proposed method accurately detects the presence of phaseamplitude coupling in these simulated data.
The proposed method accurately detects AAC
We next consider signals with amplitudeamplitude coupling, but without phaseamplitude coupling. We simulate a 20 s signal such that ${A}_{\text{high}}$ is modulated by ${A}_{\text{low}}$ (see Materials and methods, ${I}_{\mathrm{A}\mathrm{A}\mathrm{C}}=1$); when $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is large, so is $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ (Figure 5H). We then estimate the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, and the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model (example fits in Figure 5I). We find that the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model increases along the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ axis, and that the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model closely follows this trend, while the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model remains mostly flat, as expected.
Simulating 1000 instances of these signals we find that ${p}_{\text{AAC}}<0.05$ for 97.9% of simulations, while ${p}_{\text{PAC}}<0.05$ for 0.3% of simulations. The significant values of ${\text{\mathbf{R}}}_{\text{AAC}}$ lie above 0 (Figure 5J), and increases in the intensity of AAC produce increases in ${\text{\mathbf{R}}}_{\text{AAC}}$ (Figure 5K). We conclude that the proposed method accurately detects the presence of amplitudeamplitude coupling.
The proposed method accurately detects the simultaneous occurrence of PAC and AAC
We now consider signals that possess both phaseamplitude coupling and amplitudeamplitude coupling. To do so, we simulate time series data with both AAC and PAC (Figure 5L). In this case, ${A}_{\text{high}}$ increases when ${\varphi}_{\text{low}}$ is near 0 radians and when ${A}_{\text{low}}$ is large (see Materials and methods, ${I}_{\mathrm{P}\mathrm{A}\mathrm{C}}=1$ and ${I}_{\mathrm{A}\mathrm{A}\mathrm{C}}=1$). We then estimate the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, the ${A}_{\text{low}}$ model, and the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model from the data and visualize the results (Figure 5M). We find that the $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model increases near ${\varphi}_{\text{low}}=0$, and that the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model increases linearly with ${A}_{\text{low}}$. The $A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model exhibits both of these behaviors, increasing at ${\varphi}_{\text{low}}=0$ and as ${A}_{\text{low}}$ increases.
Simulating 1000 instances of signals with both AAC and PAC present, we find that ${p}_{\text{AAC}}<0.05$ in 96.7% of simulations and ${p}_{\text{PAC}}<0.05$ in 98.1% of simulations. The distributions of significant ${\text{\mathbf{R}}}_{\text{PAC}}$ and ${\text{\mathbf{R}}}_{\text{AAC}}$ values lie above 0, consistent with the presence of both PAC and AAC (Figure 5N), and as the intensity of PAC and AAC increases, so do the values of ${\text{\mathbf{R}}}_{\text{PAC}}$ and ${\text{\mathbf{R}}}_{\text{AAC}}$ (Figure 5O). We conclude that the model successfully detects the concurrent presence of PAC and AAC.
${\text{\mathbf{R}}}_{\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{\mathbf{R}}}_{\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{\mathbf{R}}}_{\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{\mathbf{R}}}_{\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{\mathbf{R}}}_{\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{\mathbf{R}}}_{\text{PAC}}$ to very weak phaseamplitude coupling, ${\text{\mathbf{R}}}_{\text{PAC}}$ can detect phaseamplitude coupling at MI values consistent with those observed in the literature.
The proposed method is less affected by fluctuations in lowfrequency amplitude and AAC
Increases in low frequency power can increase measures of phaseamplitude 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:
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{\mathbf{R}}}_{\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{\mathbf{R}}}_{\text{PAC}}$. However, this increase is more dramatic for MI than for ${\text{\mathbf{R}}}_{\text{PAC}}$; we note that the distributions of ${\text{\mathbf{R}}}_{\text{PAC}}$ almost completely overlap (Figure 7A), while the distribution of MI shifts to larger values when the amplitude of $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ increases (Figure 7B). We conclude that the statistic ${\text{\mathbf{R}}}_{\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.
We also investigate the effect of increases in amplitudeamplitude 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{\mathbf{R}}}_{\text{PAC}}$ and MI values. We observe that the increase in AAC produces a small increase in the distribution of ${\text{\mathbf{R}}}_{\text{PAC}}$ values (Figure 7C), but a large increase in the distribution of MI values (Figure 7D). We conclude that the statistic ${\text{\mathbf{R}}}_{\text{PAC}}$ is more robust to increases in AAC than MI.
These simulations show that at a fixed, nonzero PAC, the modulation index increases with increased $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and AAC. We now consider the scenario of increased $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 ${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 ${I}_{\mathrm{A}\mathrm{A}\mathrm{C}}=0$ for $t<100\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$ and intensity ${I}_{\mathrm{A}\mathrm{A}\mathrm{C}}=2$ for $t\ge 100\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}$). We find that, in the absence of PAC, ${\text{\mathbf{R}}}_{\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 amplitudeamplitude coupling, MI may detect PAC where none exists, while $\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 $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ depends not only on $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, but also on $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 $\varphi}_{\text{low}$ occurrences. $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$, on the other hand, may detect these sparse events if these events are coupled to $A}_{\text{low}$, as $\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 $V$ with PAC, and corresponding modulation signal M with intensity value ${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 $A}_{\text{low}$. To do so, we preserve PAC at the peaks of $V}_{\text{low}$ (i.e., when ${\varphi}_{\text{low}}=0$), but now only when these peaks are large, more specifically in the top 5% of peak values.
We define a threshold value $T$ to be the 95^{th} quantile of the peak $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ values, and modify the modulation signal M as follows. When M exceeds 1 (i.e., when ${\varphi}_{\text{low}}=0$) and the low frequency amplitude exceeds $T$ (i.e., ${A}_{\text{low}}\ge T$), we make no change to M. Alternatively, when M exceeds one and the low frequency amplitude lies below $T$ (i.e., ${A}_{\text{low}}<T)$, we decrease M to 1 (Figure 8C). In this way, we create a modified modulation signal $\mathbf{\text{M}}}_{1$ such that in the resulting signal $V}_{1$, when ${\varphi}_{\text{low}}=0$ and $A}_{\text{low}$ is large enough, $A}_{\text{high}$ is increased; and when ${\varphi}_{\text{low}}=0$ and $A}_{\text{low}$ is not large enough, there is no change to $A}_{\text{high}$. This signal $V}_{1$ hence has fewer phaseamplitude coupling events than the number of times ${\varphi}_{\text{low}}=0$.
We generate 1000 realizations of the simulated signals $V}_{1$, and compute $\mathbf{\text{R}}}_{\text{PAC}$ and MI. We find that while MI detects significant PAC in only 37% of simulations, $\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 $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, and $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$, which accounts for changes in $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, successfully detects these events much more frequently. We conclude that when the PAC is dependent on $A}_{\text{low}$, $\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 $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ impacts the occurrence of PAC. More specifically, we consider a case in which $A}_{\text{high}$ increases at a fixed low frequency phase for high values of $A}_{\text{low}$, and $A}_{\text{high}$ decreases at the same phase for small values of $A}_{\text{low}$. In this case, we expect that the modulation index may fail to detect the coupling because the distribution of $A}_{\text{high}$ over $\varphi}_{\text{low}$ would appear uniform when averaged over all values of $A}_{\text{low}$; the dependence of $A}_{\text{high}$ on $\varphi}_{\text{low}$ would only become apparent after accounting for $A}_{\text{low}$.
To implement this scenario, we consider the modulation signal M (see Materials and methods) with an intensity value ${I}_{\mathrm{P}\mathrm{A}\mathrm{C}}=1$. We consider all peaks of $A}_{\text{low}$ and set the threshold $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 ${\varphi}_{\text{low}}=0$) and the low frequency amplitude exceeds $T$ (i.e., ${A}_{\text{low}}\ge T$), we make no change to M. Alternatively, when M exceeds one and the low frequency amplitude lies below $T$ (i.e. ${A}_{\text{low}}<T)$, we decrease M to 0 (Figure 9B). In this way, we create a modified modulation signal M such that when ${\varphi}_{\text{low}}=0$ and $A}_{\text{low}$ is large enough, $A}_{\text{high}$ is increased; and when ${\varphi}_{\text{low}}=0$ and $A}_{\text{low}$ is small enough, $A}_{\text{high}$ is decreased (Figure 9C).
Using this method, we simulate 1000 realizations of this signal, and calculate MI and $\mathbf{\text{R}}}_{\text{PAC}$ for each signal (Figure 9D). We find that $\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, $\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 crossfrequency 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 ($V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$) in which an externally imposed signal $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$ modulates the probability of spiking as a function of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $\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 ($\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast$) modulating $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ to change as a function of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$: when $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is large, the probability of spiking is highest near ${\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}=\pm \pi$, and when $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is small, the probability of spiking is highest near ${\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}=0.$ More precisely, we define $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast$ as
where $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, $\lambda$, as
where $\sigma =0.01$, $s(\varphi )$ is a triangle wave, and we choose $\lambda}_{0$ so that the maximum value of $\lambda$ is 2. We note that the spiking probability $\lambda$ is zero except near times when the phase of the low frequency signal $({\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}})$ is near $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\ast$. We then define $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ as:
where $S$ is the binary sequence generated by the stochastic spiking neuron model, and $n$ is Gaussian noise with mean zero and standard deviation 0.1. In this scenario, the distribution of $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ over $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ appears uniform when averaged over all values of $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 $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$. However, we expect that $\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 $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ detects PAC (${\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}}=0.172$, $p=0.02$); we note that the ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$) surface exhibits a single peak near ${\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}=0$ at small values of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, and at ${\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}=\pm \pi$ at large value of $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ (Figure 10B). The ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$) surface deviates significantly from the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ surface, resulting in a large $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ value. However, the nonuniform shape of the ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$) surface is lost when we fail to account for $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$. In this scenario, the distribution of $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ over $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$ appears uniform, resulting in a low MI value (Figure 10C).
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 $V$, filtering into $V}_{\text{high}$ (100–140 Hz) and $V}_{\text{low}$ (4–7 Hz), and extracting $A}_{\text{high}$, $A}_{\text{low}$, and $\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.
We find during this interval significant phaseamplitude coupling computed using $\mathbf{\text{R}}}_{\text{PAC}$ (${\mathbf{\text{R}}}_{\text{PAC}}=1.55$, ${p}_{\text{PAC}}=0.005$, Figure 12), and using the modulation index ($\mathbf{\text{MI}}=0.03$, $p}_{\text{MI}}=5.0\times {10}^{4$). To examine the phaseamplitude coupling in more detail, we isolate a 2 s segment (Figure 11D) and display the signal $V$, the high frequency signal $V}_{\text{high}$, the low frequency phase $\varphi}_{\text{low}$, and the low frequency amplitude $A}_{\text{low}$. We observe that when $\varphi}_{\text{low}$ is near $\pi$, the amplitude of $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ tends to increase, consistent with the presence of PAC and a significant value of $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and MI.
We also find significant amplitudeamplitude coupling computed using $\mathbf{\text{R}}}_{\text{AAC}$ (${\mathbf{\text{R}}}_{\text{AAC}}=0.85$, ${p}_{\text{AAC}}=0.005$). Comparing $A}_{\text{high}$ and $A}_{\text{low}$ over the 10 s interval (each smoothed using a 1 s moving average filter and normalized), we observe that both $A}_{\text{high}$ and $A}_{\text{low}$ steadily increase over the duration of the interval (Figure 11E).
Application to in vivo rodent data
As a second example to illustrate the performance of the new method, we consider LFP recordings from from the infralimbic cortex (IL) and basolateral amygdala (BLA) of an outbred LongEvans rat before and after the delivery of an experimental electrical stimulation intervention described in Blackwood et al. (2018). Eight microwires in each region, referenced as bipolar pairs, sampled the LFP at 30 kHz, and electrical stimulation was delivered to change interregional coupling (see Blackwood et al., 2018 for a detailed description of the experiment). Here we examine how crossfrequency coupling between low frequency (5–8 Hz) IL signals and high frequency (70–110 Hz) BLA signals changes from the prestimulation to the poststimulation condition. To do so, we filter the data $V$ into low and high frequency signals (see Materials and methods), and compute the MI, $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$ between each possible BLAIL pairing, sixteen in total.
We find three separate BLAIL pairings where $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C}$ reports no significant PAC pre or poststimulation, but MI reports significant coupling poststimulation. Investigating further, we note that in all three cases, the amplitude of the low frequency IL signal increases from pre to poststimulation, and $\mathbf{\text{R}}}_{\mathrm{A}\mathrm{A}\mathrm{C}$, the measure of amplitudeamplitude coupling, increases from pre to poststimulation. These observations are consistent with the simulations in Results: The proposed method is less affected by fluctuations in lowfrequency amplitude and AAC, in which we showed that increases in the low frequency amplitude and AAC produced increases in MI, although the PAC remained fixed. We therefore propose that, consistent with these simulation results, the increase in MI observed in these data may result from changes in the low frequency amplitude and AAC, not in PAC.
Using the flexibility of GLMs to improve detection of phaseamplitude 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 prestimulation or poststimulation, affects the coupling. To incorporate this new covariate into the framework, we consider the concatenated voltage recordings from the prestimulation condition $V}_{\mathrm{p}\mathrm{r}\mathrm{e}$ and the poststimulation condition $V}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}$:
From $V$, we obtain the corresponding high frequency signal $V}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ and low frequency signal $V}_{\mathrm{l}\mathrm{o}\mathrm{w}$, and subsequently the high frequency amplitude $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, low frequency phase $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, and low frequency amplitude $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$. We use these data to generate two new models:
where $P$ is an indicator function specifying whether the signal is in the prestimulation ($P=0$) or poststimulation ($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 $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 $A}_{\text{high}$ by including the interaction between $P$ and $\varphi}_{\text{low}$. In that case, the impact of $\varphi}_{\text{low}$ on $A}_{\text{high}$ can be modeled without considering stimulus condition $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 $S}_{P{\varphi}_{\text{low}}$ from the model in Equation 9, and a surface $S}_{P$ from the model in Equation 10 in the ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $A}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, P) space, assessing the models at $\mathit{\text{P}}=1$. We compute $\mathbf{\text{R}}}_{\text{PAC},\phantom{\rule{thinmathspace}{0ex}}\text{condition}$, which measures the impact of stimulus condition on PAC, as:
We find for the example rodent data an $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C},\text{}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$ value of 0.3608, with a pvalue 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 (${I}_{\text{PAC}}=0$) and nonzero PAC for the final 20 s $({I}_{\text{PAC}}=1)$. We design this simulation to mimic an increase in PAC from prestimulation to poststimulation (Figure 13A). Using the models in Equations 9 and 10, and computing $\mathbf{\text{R}}}_{\mathrm{P}\mathrm{A}\mathrm{C},\text{}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$, we find $p<0.05$ for 100% of simulated signals. We also simulate 1000 instances of a 40 s signal with no PAC (${I}_{\text{PAC}}=0$) for the entire 40 s, that is PAC does not change from prestimulation to poststimulation (Figure 13B), and find in this case $p<0.05$ for only 4.6% of simulations. Finally, we simulate 1000 instances of a 40 s signal with fixed PAC $({I}_{\text{PAC}}=1)$, and with a doubling of the low frequency amplitude occuring at 20 s (i.e., prestimulation the low frequency amplitude is 1, and poststimulation the low frequency amplitude is 2). We find $p<0.05$ for only 3.6% of simulations. We conclude that this method effectively determines whether stimulation condition significantly changes PAC.
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 poststimulation). 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.
Discussion
In this paper, we proposed a new method for measuring crossfrequency coupling that accounts for both phaseamplitude coupling and amplitudeamplitude 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 lowfrequency 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 lowfrequency 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 ($A}_{\text{low}$) is large, the proposed method detects this coupling while another method not accounting for $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 $A}_{\text{high}$) at a particular phase of the low frequency rhythm ($\varphi}_{\text{low}$) when excitation is maximal. These notions also provide a physical interpretation for AAC: increases in $A}_{\text{low}$ produce larger modulations in neural excitability, and therefore increased intervals of neuronal spiking (i.e., increases in $A}_{\text{high}$). Alternatively, decreases in $A}_{\text{low}$ reduce excitability and neuronal spiking (i.e., decreases in $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 multiitem encoding, longdistance 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 $\varphi}_{\text{low}$. PAC can be interpreted as a type of focused attention; $A}_{\text{high}$ modulation occurring only in a particular interval of $\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 attentionbased 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 $\varphi}_{\text{low}$. Second, the value of $\mathbf{\text{R}}}_{\text{AAC}$ depends on the range of $A}_{\text{low}$ observed. This is due to the linear relationship between $A}_{\text{low}$ and $A}_{\text{high}$ in the $A}_{\mathrm{l}\mathrm{o}\mathrm{w}$ model, which causes the maximum distance between the surfaces $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ to occur at the largest or smallest value of $A}_{\text{low}$. To mitigate the impact of extreme $A}_{\text{low}$ values on $\mathbf{\text{R}}}_{\text{AAC}$, we evaluate the surfaces $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and $S}_{{A}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ over the 5^{th} to 95^{th} quantiles of $A}_{\text{low}$. We note that an alternative metric of AAC could instead evaluate the slope of the $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 $V}_{\text{high}$ and $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 $V}_{\text{high}$ and $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 nonsinusoidal signals, which may confound interpretation of crossfrequency 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 nonsinusoidal 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 $A}_{\text{high}$ and ($\varphi}_{\text{low}$,$A}_{\text{low}$) may result in time series with simpler structure than those observed in vivo. For example, a latent signal may drive both $A}_{\text{high}$ and $\varphi}_{\text{low}$, and in this way establish nonlinear relationships between the two observables $A}_{\text{high}$ and $\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 $\mathbf{\text{R}}}_{\text{PAC}$ and $\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 $\mathbf{\text{R}}}_{\text{PAC}$ and $\mathbf{\text{R}}}_{\text{AAC}$. While the statistics $\mathbf{\text{R}}}_{\text{PAC}$ and $\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 $\mathbf{\text{R}}}_{\text{PAC}$ and $\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 realworld neural data analysis.
The proposed method can easily be extended by inclusion of additional predictors in the GLM. Polynomial $A}_{\text{low}$ predictors, rather than the current linear $A}_{\text{low}$ predictors, may better capture the relationship between $A}_{\text{low}$ and $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 ($A}_{\mathrm{l}\mathrm{o}\mathrm{w}$, $\varphi}_{\mathrm{l}\mathrm{o}\mathrm{w}$) and $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 $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/EdenKramerLab/GLMCFC, 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 crossfrequency 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.
Data availability
In vivo human data available at https://github.com/EdenKramerLab/GLMCFC (copy archived at https://github.com/elifesciencespublications/GLMCFC). In vivo rat data available at https://github.com/tnelab/clexampledata (copy archived at https://github.com/elifesciencespublications/clexampledata).
References

Untangling crossfrequency coupling in neuroscienceCurrent Opinion in Neurobiology 31:51–61.https://doi.org/10.1016/j.conb.2014.08.002

What does the broken brain say to the neuroscientist? oscillations and connectivity in schizophrenia, Alzheimer’s disease, and bipolar disorderInternational Journal of Psychophysiology 103:135–148.https://doi.org/10.1016/j.ijpsycho.2015.02.004

Continuous phase estimation for phaselocked neural stimulation using an autoregressive model for signal predictionConference of the IEEE Engineering in Medicine and Biology Society 2018:4736–4739.https://doi.org/10.1109/EMBC.2018.8513232

Gamma (40100 hz) oscillation in the Hippocampus of the behaving ratThe Journal of Neuroscience 15:47–60.https://doi.org/10.1523/JNEUROSCI.150100047.1995

Taskrelated coupling from high to lowfrequency signals among visual cortical Areas in human subdural recordingsInternational Journal of Psychophysiology 51:97–116.https://doi.org/10.1016/j.ijpsycho.2003.07.001

Neuronal oscillations in cortical networksScience 304:1926–1929.https://doi.org/10.1126/science.1099745

Mechanisms of gamma oscillationsAnnual Review of Neuroscience 35:203–225.https://doi.org/10.1146/annurevneuro062111150444

The functional role of crossfrequency couplingTrends in Cognitive Sciences 14:506–515.https://doi.org/10.1016/j.tics.2010.09.001

Gamma oscillations in the entorhinal cortex of the freely behaving ratThe Journal of Neuroscience 18:388–398.https://doi.org/10.1523/JNEUROSCI.180100388.1998

Assessing transient crossfrequency coupling in EEG dataJournal of Neuroscience Methods 168:494–499.https://doi.org/10.1016/j.jneumeth.2007.10.012

Brain oscillations and the importance of waveform shapeTrends in Cognitive Sciences 21:137–149.https://doi.org/10.1016/j.tics.2016.12.008

Dynamic predictions: oscillations and synchrony in topdown processingNature Reviews Neuroscience 2:704–716.https://doi.org/10.1038/35094565

Analytical insights on thetagamma coupled neural oscillatorsThe Journal of Mathematical Neuroscience 3:16.https://doi.org/10.1186/21908567316

Neural CrossFrequency coupling: connecting architectures, mechanisms, and functionsTrends in Neurosciences 38:725–740.https://doi.org/10.1016/j.tins.2015.09.001

Position reconstruction from an ensemble of hippocampal place cells: contribution of theta phase codingJournal of Neurophysiology 83:2602–2609.https://doi.org/10.1152/jn.2000.83.5.2602

Crossfrequency coupling in real and virtual brain networksFrontiers in Computational Neuroscience 7:78.https://doi.org/10.3389/fncom.2013.00078

4Hz oscillations synchronize prefrontalamygdala circuits during fear behaviorNature Neuroscience 19:605–612.https://doi.org/10.1038/nn.4251

Sharp edge artifacts and spurious coupling in EEG frequency comodulation measuresJournal of Neuroscience Methods 170:352–357.https://doi.org/10.1016/j.jneumeth.2008.01.020

Assessment of crossfrequency coupling with confidence using generalized linear modelsJournal of Neuroscience Methods 220:64–74.https://doi.org/10.1016/j.jneumeth.2013.08.006

BookCase Studies in Neural Data Analysis: A Guide for the Practicing NeuroscientistThe MIT Press.

Measuring phase synchrony in brain signalsHuman Brain Mapping 8:194–208.https://doi.org/10.1002/(SICI)10970193(1999)8:4<194::AIDHBM4>3.0.CO;2C

An oscillatory hierarchy controlling neuronal excitability and stimulus processing in the auditory cortexJournal of Neurophysiology 94:1904–1911.https://doi.org/10.1152/jn.00263.2005

Phase resetting reduces thetagamma rhythmic interaction to a onedimensional mapJournal of Mathematical Biology 66:1361–1386.https://doi.org/10.1007/s0028501205349

Quantifying phaseamplitude coupling in neuronal network oscillationsProgress in Biophysics and Molecular Biology 105:49–57.https://doi.org/10.1016/j.pbiomolbio.2010.09.007

Testing for nested oscillationJournal of Neuroscience Methods 174:50–61.https://doi.org/10.1016/j.jneumeth.2008.06.035

Inferring spike trains from local field potentialsJournal of Neurophysiology 99:1461–1476.https://doi.org/10.1152/jn.00919.2007

Bifurcation analysis on PhaseAmplitude CrossFrequency coupling in neural networks with dynamic synapsesFrontiers in Computational Neuroscience 11:18.https://doi.org/10.3389/fncom.2017.00018

On highfrequency field oscillations (>100 hz) and the spectral leakage of spiking activityJournal of Neuroscience 33:1535–1539.https://doi.org/10.1523/JNEUROSCI.421712.2013

Testing for nonlinearity in time series: the method of surrogate dataPhysica D: Nonlinear Phenomena 58:77–94.https://doi.org/10.1016/01672789(92)90102S

Measuring phaseamplitude coupling between neuronal oscillations of different frequenciesJournal of Neurophysiology 104:1195–1210.https://doi.org/10.1152/jn.00106.2010

RespirationEntrained brain rhythms are global but often overlookedTrends in Neurosciences 41:186–197.https://doi.org/10.1016/j.tins.2018.01.007

Parametric estimation of crossfrequency couplingJournal of Neuroscience Methods 243:94–102.https://doi.org/10.1016/j.jneumeth.2015.01.032

The brainweb: phase synchronization and largescale integrationNature Reviews Neuroscience 2:229–239.https://doi.org/10.1038/35067550

Gammaphase shifting in awake monkey visual cortexJournal of Neuroscience 30:1250–1257.https://doi.org/10.1523/JNEUROSCI.162309.2010

Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronizationInternational Journal of Psychophysiology 38:301–313.https://doi.org/10.1016/S01678760(00)001720

Inhibitionbased rhythms: experimental and mathematical observations on network dynamicsInternational Journal of Psychophysiology 38:315–336.https://doi.org/10.1016/S01678760(00)001732

Multiple origins of the cortical γ rhythmDevelopmental Neurobiology 71:92–106.https://doi.org/10.1002/dneu.20814

Coherent neuronal ensembles are rapidly recruited when making a lookreach decisionNature Neuroscience 19:327–334.https://doi.org/10.1038/nn.4210

Decoding hidden cognitive states from behavior and physiology using a bayesian approachNeural Computation 31:1751–1788.https://doi.org/10.1162/neco_a_01196
Article and author information
Author details
Funding
National Science Foundation (NSF DMS #1451384)
 Jessica K Nadalin
 Mark A Kramer
National Science Foundation (GRFP)
 Jessica K Nadalin
National Institutes of Health (R21 MH109722)
 Alik S Widge
National Institutes of Health (R01 EB026938)
 Alik S Widge
 Uri T Eden
 Mark A Kramer
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was supported in part by the National Science Foundation Award #1451384, in part by R01 EB026938, in part by R21 MH109722, and in part by the National Science Foundation (NSF) under a Graduate Research Fellowship.
Ethics
Human subjects: All patients were enrolled after informed consent, and consent to publish, was obtained and approval was granted by local Institutional Review Boards at Massachusetts General Hospital and Brigham Women's Hospitals (Partners Human Research Committee), and at Boston University according to National Institutes of Health guidelines (IRB Protocol # 1558X).
Animal experimentation: The animal experimentation received IACUC approval from the University of Minnesota (IACUC Protocol # 180636024A).
Version history
 Received: December 12, 2018
 Accepted: October 6, 2019
 Accepted Manuscript published: October 16, 2019 (version 1)
 Version of Record published: October 30, 2019 (version 2)
Copyright
© 2019, Nadalin et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 4,167
 views

 590
 downloads

 9
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Biochemistry and Chemical Biology
 Neuroscience
In most mammals, conspecific chemosensory communication relies on semiochemical release within complex bodily secretions and subsequent stimulus detection by the vomeronasal organ (VNO). Urine, a rich source of ethologically relevant chemosignals, conveys detailed information about sex, social hierarchy, health, and reproductive state, which becomes accessible to a conspecific via vomeronasal sampling. So far, however, numerous aspects of social chemosignaling along the vomeronasal pathway remain unclear. Moreover, since virtually all research on vomeronasal physiology is based on secretions derived from inbred laboratory mice, it remains uncertain whether such stimuli provide a true representation of potentially more relevant cues found in the wild. Here, we combine a robust lownoise VNO activity assay with comparative molecular profiling of sex and strainspecific mouse urine samples from two inbred laboratory strains as well as from wild mice. With comprehensive molecular portraits of these secretions, VNO activity analysis now enables us to (i) assess whether and, if so, how much sex/strainselective ‘raw’ chemical information in urine is accessible via vomeronasal sampling; (ii) identify which chemicals exhibit sufficient discriminatory power to signal an animal’s sex, strain, or both; (iii) determine the extent to which wild mouse secretions are unique; and (iv) analyze whether vomeronasal response profiles differ between strains. We report both sex and, in particular, strainselective VNO representations of chemical information. Within the urinary ‘secretome’, both volatile compounds and proteins exhibit sufficient discriminative power to provide sex and strainspecific molecular fingerprints. While total protein amount is substantially enriched in male urine, females secrete a larger variety at overall comparatively low concentrations. Surprisingly, the molecular spectrum of wild mouse urine does not dramatically exceed that of inbred strains. Finally, vomeronasal response profiles differ between C57BL/6 and BALB/c animals, with particularly disparate representations of female semiochemicals.

 Neuroscience
Midbrain dopamine neurons impact neural processing in the prefrontal cortex (PFC) through mesocortical projections. However, the signals conveyed by dopamine projections to the PFC remain unclear, particularly at the singleaxon level. Here, we investigated dopaminergic axonal activity in the medial PFC (mPFC) during reward and aversive processing. By optimizing microprismmediated twophoton calcium imaging of dopamine axon terminals, we found diverse activity in dopamine axons responsive to both reward and aversive stimuli. Some axons exhibited a preference for reward, while others favored aversive stimuli, and there was a strong bias for the latter at the population level. Longterm longitudinal imaging revealed that the preference was maintained in reward and aversivepreferring axons throughout classical conditioning in which rewarding and aversive stimuli were paired with preceding auditory cues. However, as mice learned to discriminate reward or aversive cues, a cue activity preference gradually developed only in aversivepreferring axons. We inferred the trialbytrial cue discrimination based on machine learning using anticipatory licking or facial expressions, and found that successful discrimination was accompanied by sharper selectivity for the aversive cue in aversivepreferring axons. Our findings indicate that a group of mesocortical dopamine axons encodes aversiverelated signals, which are modulated by both classical conditioning across days and trialbytrial discrimination within a day.