On crossfrequency phasephase coupling between theta and gamma oscillations in the hippocampus
 Cited 39
 Views 3,741
 Annotations
Abstract
Phaseamplitude coupling between theta and multiple gamma subbands is a hallmark of hippocampal activity and believed to take part in information routing. More recently, theta and gamma oscillations were also reported to exhibit phasephase coupling, or n:m phaselocking, suggesting an important mechanism of neuronal coding that has long received theoretical support. However, by analyzing simulated and actual LFPs, here we question the existence of thetagamma phasephase coupling in the rat hippocampus. We show that the quasilinear phase shifts introduced by filtering lead to spurious coupling levels in both white noise and hippocampal LFPs, which highly depend on epoch length, and that significant coupling may be falsely detected when employing improper surrogate methods. We also show that waveform asymmetry and frequency harmonics may generate artifactual n:m phaselocking. Studies investigating phasephase coupling should rely on appropriate statistical controls and be aware of confounding factors; otherwise, they could easily fall into analysis pitfalls.
https://doi.org/10.7554/eLife.20515.001eLife digest
Neuroscientists have long sought to understand how the brain works by analyzing its electrical activity. Placing electrodes on the scalp or lowering them into the brain itself reveals rhythmic waves of activity known as oscillations. These arise when large numbers of neurons fire in synchrony. Recordings reveal that the frequency of these oscillations – the number of cycles of a wave per second, measured in Hertz – can vary between brain regions, and within a single region over time. Moreover, oscillations with different frequencies can coexist and interact with one another.
Within the hippocampus, an area of the brain involved in memory, two types of oscillations dominate: theta waves and gamma waves. Theta waves are relatively slow waves, with a frequency between 5 and 10 Hertz. Gamma waves are faster, with a frequency of up to 100 Hertz. Recent work has suggested that gamma waves and theta waves show a phenomenon called phasephase coupling. Since gamma waves are faster than theta waves, multiple cycles of gamma can occur during a single cycle of theta. Phasephase coupling is the idea that gamma and theta waves align themselves, such that gamma waves always begin at the same relative position within a theta wave. This was thought to help the hippocampus to encode memories.
Using computer simulations and recordings from the rat hippocampus, SchefferTeixeira and Tort have now reexamined the evidence for thetagamma phasephase coupling. The new results suggest that previous reports describing the phenomenon may have relied on inadequate statistical techniques. Using stringent control analyses, SchefferTeixeira and Tort find no evidence for prominent thetagamma phasephase coupling in the hippocampus. Instead, the simulations suggest that what appeared to be statistically significant coupling may in reality be an artifact of the previous analysis.
Phasephase coupling of theta and gamma waves has also been reported in the human hippocampus. The next step therefore is to apply these more robust analysis techniques to data from the human brain. While revisiting previously accepted findings may not always be popular, it will likely be essential if neuroscientists want to accurately understand how new memories are formed.
https://doi.org/10.7554/eLife.20515.002Introduction
Local field potentials (LFPs) exhibit oscillations of different frequencies, which may cooccur and also interact with one another (Jensen and Colgin, 2007; Tort et al., 2010; Hyafil et al., 2015). Crossfrequency phaseamplitude coupling between theta and gamma oscillations has been well described in the hippocampus, whereby the instantaneous amplitude of gamma oscillations depends on the instantaneous phase of theta (SchefferTeixeira et al., 2012; Schomburg et al., 2014). More recently, hippocampal theta and gamma oscillations were also reported to exhibit n:m phasephase coupling, in which multiple gamma cycles are consistently entrained within one cycle of theta (Belluscio et al., 2012; Zheng and Zhang, 2013; Xu et al., 2013, 2015; Zheng et al., 2016). The existence of different types of crossfrequency coupling suggests that the brain may use different coding strategies to transfer multiplexed information.
Coherent oscillations are believed to take part in network communication by allowing opportunity windows for the exchange of information (Varela et al., 2001; Fries, 2005). Standard phase coherence measures the constancy of the phase difference between two oscillations of the same frequency (Lachaux et al., 1999; Hurtado et al., 2004), and has been associated with cognitive processes such as decisionmaking (DeCoteau et al., 2007; Montgomery and Buzsáki, 2007; Nácher et al., 2013). Similarly to coherence, crossfrequency phase–phase coupling, or n:m phaselocking, also relies on assessing the constancy of the difference between two phase time series (Tass et al., 1998). However, in this case the original phase time series are accelerated, so that their instantaneous frequencies can match. Formally, n:m phaselocking occurs when $\mathrm{\Delta}{\phi}_{nm}\left(t\right)=n\ast {\phi}_{B}\left(t\right)m\ast {\phi}_{A}\left(t\right)$ is nonuniform but centered around a preferred value, where $n*{\phi}_{B}\left(m*{\phi}_{A}\right)$ denotes the phase of oscillation B (A) accelerated n (m) times (Tass et al., 1998). For example, the instantaneous phase of theta oscillations at 8 Hz needs to be accelerated five times to match in frequency a 40 Hz gamma. A 1:5 phasephase coupling is then said to occur if theta accelerated five times has a preferred phase lag (i.e., a nonuniform phase difference) in relation to gamma; or, in other words, if five gamma cycles have a consistent phase relationship to one theta cycle.
Crossfrequency phasephase coupling has previously been hypothesized to take part in memory processes (Lisman and Idiart, 1995; Jensen and Lisman, 2005; Lisman, 2005; Schack and Weiss, 2005; Sauseng et al., 2008, 2009; Holz et al., 2010; Fell and Axmacher, 2011). Recent findings suggest that the hippocampus indeed uses such a mechanism (Belluscio et al., 2012; Zheng and Zhang, 2013; Xu et al., 2013, 2015; Zheng et al., 2016). However, by analyzing simulated and actual hippocampal LFPs, in the present work we question the existence of thetagamma phasephase coupling.
Results
Measuring n:m phaselocking
We first certified that we could reliably detect n:m phaselocking when present. To that end, we simulated a system of two Kuramoto oscillators – a ‘theta’ and a ‘gamma’ oscillator – exhibiting variability in instantaneous frequency (see Materials and methods). The mean natural frequency of the theta oscillator was set to 8 Hz, while the mean natural frequency of the gamma oscillator was set to 43 Hz (Figure 1A). When coupled, the mean frequencies aligned to a 1:5 factor by changing to 8.5 Hz and 42.5 Hz, respectively (see Guevara and Glass, 1982; GarcíaAlvarez et al., 2008; Canavier et al., 2009). Figure 1B depicts three versions of accelerated theta phases (m = 3, 5 and 7) along with the instantaneous gamma phase (n = 1) of the coupled oscillators (see Figure 1—figure supplement 1 for the uncoupled case). Also shown are the time series of the difference between gamma and accelerated theta phases ($\mathrm{\Delta}{\phi}_{nm}$). The instantaneous phase difference has a preferred lag only for m = 5; when m = 3 or 7, $\mathrm{\Delta}{\phi}_{nm}$changes over time, precessing forwards (m = 3) or backwards (m = 7) at an average rate of 17 Hz. Consequently, $\mathrm{\Delta}{\phi}_{nm}$ distribution is uniform over 0 and 2π for m = 3 or 7, but highly concentrated for m = 5 (Figure 1C). The concentration (or ‘constancy’) of the phase difference distribution is used as a metric of n:m phaselocking. This metric is defined as the length of the mean resultant vector (R_{n:m}) over unitary vectors whose angle is the instantaneous phase difference ($e}^{i\mathrm{\Delta}{\phi}_{nm}\left(t\right)$), and thereby it varies between 0 and 1. For any pair of phase time series, an R_{n:m}‘curve’ can be calculated by varying m for n = 1 fixed. As shown in Figure 1D, the coupled – but not uncoupled – oscillators exhibited a prominent peak for n:m = 1:5, which shows that R_{n:m} successfully detects n:m phaselocking.
Filteringinduced n:m phaselocking in white noise
We next analyzed whitenoise signals, in which by definition there is no structured activity; in particular, the spectrum is flat and there is no true n:m phaselocking. R_{n:m} values measured from white noise should be regarded as chance levels. We bandpass filtered whitenoise signals to extract the instantaneous phase of theta (θ: 4–12 Hz) and of multiple gamma bands (Figure 2A): slow gamma (γ_{S}: 30–50 Hz), middle gamma (γ_{M}: 50–90 Hz), and fast gamma (γ_{F}: 90–150 Hz). For each frequency pair, we constructed n:m phaselocking curves for epochs of 1 and 10 s, with n = 1 fixed and m varying from 1 to 25 (Figure 2B). In each case, phasephase coupling was high within the ratio of the analyzed frequency ranges: R_{n:m} peaked at m = 4–6 for θ−γ_{S}, at m = 7–11 for θ−γ_{M}, and at m = 12–20 for θ−γ_{F}. Therefore, the existence of a ‘bump’ in the R_{n:m} curve may merely reflect the ratio of the filtered bands and should not be considered as evidence for crossfrequency phasephase coupling: even filtered whitenoise signals exhibit such a pattern.
The bump in the R_{n:m} curve of filtered white noise is explained by the fact that neighboring data points are not independent. In fact, the phase shift between two consecutive data points follows a probability distribution highly concentrated around 2*π*f_{c}*dt, where f_{c} is the filter center frequency and dt the sampling period (Figure 2—figure supplement 1). For instance, for dt = 1 ms (1000 Hz sampling rate), consecutive samples of white noise filtered between 4 and 12 Hz are likely to exhibit phase difference of 0.05 rad (8 Hz center frequency); likewise, signals filtered between 30 and 50 Hz are likely to exhibit phase differences of 0.25 rad (40 Hz center frequency). In turn, the ‘sinusoidality’ imposed by filtering leads to nonzero R_{n:m} values, which peak at the ratio of the center frequencies, akin to the fact that perfect 8 Hz and 40 Hz sine waves have R_{n:m} = 1 at n:m = 1:5. In accordance to this explanation, no R_{n:m} bump occurs when data points of the gamma phase time series are made independent by subsampling with a period longer than a gamma cycle (Figure 2—figure supplement 1), or when extracting phase values from different trials (not shown). As expected, the effect of filteringinduced sinusoidality on R_{n:m} values is stronger for narrower frequency bands (Figure 2—figure supplement 2).
Qualitatively similar results were found for 1 and 10 s epochs; however, R_{n:m} values were considerably lower for the latter (Figure 2B). In fact, for any fixed n:m ratio and frequency pair, R_{n:m} decreased as a function of epoch length (see Figure 2C for θ−γ_{S} and R_{1:5}): the longer the whitenoise epoch the more the phase difference distribution becomes uniform. In other words, as standard phase coherence (Vinck et al., 2010) and phaseamplitude coupling (Tort et al., 2010), phasephase coupling has positive bias for shorter epochs. As a corollary, notice that falsepositive coupling may be detected if control (surrogate) epochs are longer than the original epoch.
Statistical testing of n:m phaselocking
We next investigated the reliability of surrogate methods for detecting n:m phaselocking (Figure 2D). The ‘Original’ R_{n:m} value uses the same time window for extracting theta and gamma phases (Figure 2D, upper panel). A ‘Time Shift’ procedure for creating surrogate epochs has been previously employed (Belluscio et al., 2012; Zheng et al., 2016), in which the time window for gamma phase is randomly shifted between 1 to 200 ms from the time window for theta phase (Figure 2D, upper middle panel). A variant of this procedure is the ‘Random Permutation’, in which the time window for gamma phase is randomly chosen (Figure 2D, lower middle panel). Finally, in the ‘Phase Scramble’ procedure, the timestamps of the gamma phase time series are shuffled (Figure 2D, lower panel); clearly, the latter is the least conservative. For each surrogate procedure, R_{n:m} values were obtained by two approaches: ‘Single Run’ and ‘Pooled’ (Figure 2E). In the first approach, each surrogate run (e.g., a time shift or a random selection of time windows) produces one R_{n:m} value (Figure 2E, top panels). In the second, $\mathrm{\Delta}{\phi}_{nm}$ from several surrogate runs are first pooled, then a single R_{n:m} value is computed from the pooled distribution (Figure 2E, bottom panel). As illustrated in Figure 2E, R_{n:m} computed from a pool of surrogate runs is much smaller than when computed for each individual run. This is due to the dependence of R_{n:m} on the epoch length: pooling instantaneous phase differences across 10 runs of 1 s surrogate epochs is equivalent to analyzing a single surrogate epoch of 10 s. And the longer the analyzed epoch, the more the noise is averaged out and the lower the R_{n:m}. Therefore, pooled surrogate epochs summing up to 10 s of total data have lower R_{n:m} than any individual 1 s surrogate epoch.
No phasephase coupling should be detected in white noise, and therefore Original R_{n:m} values should not differ from properly constructed surrogates. However, as shown in Figure 2F for θ−γ_{S} as an illustrative case (similar results hold for any frequency pair), θ−γ_{S} phasephase coupling in white noise was statistically significantly larger than in phasescrambled surrogates (for either Single Run or Pooled distributions). This was true for surrogate epochs of any length, although the longer the epoch, the lower the actual and the surrogate R_{n:m} values, as expected (compare right and left panels of Figure 2F). Pooled R_{1:5} distributions derived from either timeshifted (Figure 2F) or randomly permutated epochs (not shown) also led to the detection of false positive θ−γ_{S} phasephase coupling. On the other hand, Original R_{n:m} values were not statistically different from chance distributions when these were constructed from Single Run R_{n:m} values for either Time Shift and Random Permutation surrogate procedures (Figure 2F; see also Figure 2—figure supplement 3). We conclude that neither scrambling phases nor pooling individual surrogate epochs should be employed for statistically evaluating n:m phaselocking. Chance distributions should be derived from surrogate epochs of the same length as the original epoch and which preserve phase continuity.
To check if Single Run surrogate distributions are capable of statistically detecting true n:m phaselocking, we next simulated noisy Kuramoto oscillators as in Figure 1, but of mean natural frequencies set to 8 and 40 Hz. Original R_{1:5} values were much greater than the surrogate distribution for coupled – but not uncoupled – oscillators (Figure 3A). This result illustrates that variability in the instantaneous frequency leads to low n:m phaselocking levels for independent oscillators even when their mean frequencies are perfect integer multiples. On the other hand, coupled oscillators have high R_{n:m} because variations of their instantaneous frequencies are mutually dependent. We then proceeded to analyze simulated LFPs from a previously published model network (Kopell et al., 2010). The network has two inhibitory interneurons, called O and I cells, which spike at theta and gamma frequency, respectively (for a motivation of this model, see Tort et al., 2007). Compared to Single Run surrogate distributions, the model LFP exhibited significant n:m phaselocking only when the interneurons were coupled; R_{n:m} levels did not differ from the surrogate distribution for the uncoupled network (Figure 3B). (Note that the R_{n:m} curve also exhibited a peak for both the uncoupled network and Single Run surrogate data, which is due to the low variability in the instantaneous spike frequency of the model cells; without this variability, however, all networks would display perfect n:m phaselocking).
Spurious n:m phaselocking due to nonsinusoidal waveforms
The simulations above show that Single Run surrogates can properly detect n:m phaselocking for oscillators exhibiting variable instantaneous frequency, which is the case of hippocampal theta and gamma oscillations. However, it should be noted that high asymmetry of the theta waveform may also lead to statistically significant R_{n:m} values per se. As illustrated in Figure 4A, a nonsinusoidal oscillation such as a theta sawtooth wave can be decomposed into a sum of sine waves at the fundamental and harmonic frequencies, which have decreasing amplitude (i.e., the higher the harmonic frequency, the lower the amplitude). Importantly, the harmonic frequency components are n:m phaselocked to each other: the first harmonic exhibits a fixed 1:2 phase relationship to the fundamental frequency, the second harmonic a 1:3 relationship, and so on (Figure 4B). Of note, the higher frequency harmonics not only exhibit crossfrequency phasephase coupling to the fundamental theta frequency but also phaseamplitude coupling, since they have higher amplitude at the theta phases where the sharp deflection occurs (Figure 4C left and Figure 4—figure supplement 1; see also Kramer et al., 2008 and Tort et al., 2013).
The gammafiltered component of a theta sawtooth wave of variable peak frequency thus displays spurious gamma oscillations (i.e., theta harmonics) that have a consistent phase relationship to the theta cycle irrespective of variations in cycle length. In randomly permutated data, however, the theta phases associated with spurious gamma differ from cycle to cycle due to the variability in instantaneous theta frequency. As a result, the spurious n:m phasecoupling induced by sharp signal deflections is significantly higher than the Random Permutation/Single Run surrogate distribution (Figure 4C right and Figure 4—figure supplement 2 top row). Interestingly, the significance of this spurious effect is much lower when using the Time Shift procedure (Figure 4—figure supplement 2 bottom row), probably due to the proximity between the original and the timeshifted time series (200 ms maximum distance).
Assessing n:m phaselocking in actual LFPs
We next proceeded to analyze hippocampal CA1 recordings from seven rats, focusing on the periods of prominent theta activity (active waking and REM sleep). We found similar results between white noise and actual LFP data. Namely, R_{n:m} curves peaked at n:m ratios according to the filtered bands, and R_{n:m} values were lower for longer epochs (Figure 5A; compare with Figure 2B). As shown in Figure 5B, Original R_{n:m} values were not statistically different from a proper surrogate distribution (Random Permutation/Single Run) in epochs of up to 100 s (but see Figure 10). Noteworthy, as with whitenoise data (Figure 2F), false positive phasephase coupling would be inferred if an inadequate surrogate method were employed (Time Shift/Pooled) (Figure 5B).
We also found no difference between original and surrogate n:m phaselocking levels when employing the metric described in Sauseng et al. (2009) (Figure 5—figure supplement 1), and when estimating theta phase by interpolating phase values between 4 points of the theta cycle (trough, ascending, peak and descending points) as performed in Belluscio et al. (2012) (Figure 5—figure supplement 2). The latter was somewhat expected since the phasephase coupling results in Belluscio et al. (2012) did not depend on this particular method of phase estimation (see their Figure 6Ce). Moreover, coupling levels did not statistically differ from zero when using the pairwise phase consistency metric described in Vinck et al. (2010) (Figure 5—figure supplement 1).
We further confirmed our results by analyzing data from three additional rats recorded in an independent laboratory (Figure 5—figure supplement 3; see Materials and methods). In addition, we also found similar results in LFPs from other hippocampal layers than s. pyramidale (Figure 5—figure supplement 4), in neocortical LFPs (not shown), in currentsource density (CSD) signals (Figure 5—figure supplement 4), in independent components that isolate activity of specific gamma subbands (Schomburg et al., 2014) (Figure 5—figure supplement 5), and in transient gamma bursts (Figure 5—figure supplement 6).
On diagonal stripes in phasephase plots
Since Original R_{n:m} values were not greater than Single Run surrogate distributions, we concluded that there is lack of convincing evidence for n:m phaselocking in the hippocampal LFPs analyzed here. However, as in previous reports (Belluscio et al., 2012; Zheng et al., 2016), phasephase plots (2D histograms of theta phase vs gamma phase) of actual LFPs displayed diagonal stripes (Figure 6), which seem to suggest phasephase coupling. We next sought to investigate what causes the diagonal stripes in phasephase plots.
In Figure 7 we analyze a representative LFP with prominent theta oscillations at ~7 Hz recorded during REM sleep. Due to the nonsinusoidal shape of theta (Belluscio et al., 2012; Sheremet et al., 2016), the LFP also exhibited spectral peaks at harmonic frequencies (Figure 7A). We constructed phase–phase plots using LFP components narrowly filtered at theta and its harmonics: 14, 21, 28 and 35 Hz. Similarly to the sawtooth wave (Figure 4B), the phasephase plots exhibited diagonal stripes whose number was determined by the harmonic order (i.e., the 1^{st} harmonic exhibited two stripes, the second harmonic three stripes, the third, four stripes and the fourth, five stripes; Figure 7Bi–iv). Interestingly, when the LFP was filtered at a broad gamma band (30–90 Hz), we observed five diagonal stripes, the same number as when narrowly filtering at 35 Hz; moreover, both gamma and 35 Hz filtered signals exhibited the exact same phase lag (Figure 7Biv–v). Therefore, these results indicate that the diagonal stripes in phasephase plots may be influenced by theta harmonics. Under this interpretation, signals filtered at the gamma band would be likely to exhibit as many stripes as expected for the first theta harmonic falling within the filtered band. Consistent with this possibility, we found that the peak frequency of theta relates to the number of stripes (Figure 8).
As in previous studies (Belluscio et al., 2012; Zheng et al., 2016), phasephase plots constructed using data averaged from individual timeshifted epochs exhibited no diagonal stripes (Figure 7Bvi and Figure 8). This is because different time shifts lead to different phase lags; the diagonal stripes of individual surrogate runs that could otherwise be apparent cancel each other out when combining data across multiple runs of different lags (Figure 8—figure supplement 1). Moreover, as in Belluscio et al. (2012), the histogram counts that give rise to the diagonal stripes were deemed statistically significant when compared to the mean and standard deviation over individual counts from timeshifted surrogates (Figure 7—figure supplement 1 and Figure 8).
To gain further insight into what generates the diagonal stripes, we next analyzed whitenoise signals. As shown in Figure 9A, phasephase plots constructed from filtered whitenoise signals also displayed diagonal stripes. Since white noise has no harmonics, these results show that the sinusoidality induced by the filter can by itself lead to diagonal stripes in phasephase plots, in the same way that it leads to a bump in the R_{n:m} curve (Figure 2 and Figure 2—figure supplement 1). Importantly, as in actual LFPs, bin counts in phasephase plots of whitenoise signals were also deemed statistically significant when compared to the distribution of bin counts from timeshifted surrogates (Figure 9A). Since by definition white noise has no n:m phaselocking, we concluded that the statistical analysis of phasephase plots as originally introduced in Belluscio et al. (2012) is too liberal. Nevertheless, we found that phasephase plots of white noise were no longer statistically significant when using the same approach as in Belluscio et al. (2012) but corrected for multiple comparisons (i.e., the number of bins) by the HolmBonferroni method (the FDR correction still led to significant bins; not shown). This result was true for different epoch lengths and also when computing surrogate phasephase plots using the Random Permutation procedure (Figure 9A). Consistently, for all epoch lengths, Original R_{n:m} values fell inside the distribution of Single Run surrogate R_{n:m} values computed using either Time Shift and Random Permutation procedures (Figure 9B).
The observations above suggest that the diagonal stripes in phasephase plots of hippocampal LFPs may actually be caused by filteringinduced sinusoidality, as opposed to being an effect of theta harmonics as we first interpreted. To test this possibility, we next revisited the significance of phasephase plots of actual LFPs. For epochs of up to 100 s, we found similar results as in white noise, namely, bin counts were no longer statistically significant after correcting for multiple comparisons (HolmBonferroni method); this was true when using either the Time Shift or Random Permutation procedures (Figure 10A). Surprisingly, however, when analyzing much longer time series (10 or 20 min of concatenated periods of REM sleep), several bin counts became statistically significant when compared to randomly permutated, but not timeshifted, surrogates (Figure 10A). Moreover, this result reflected in the R_{n:m} curves: the Original R_{n:m} curve fell within the distribution of Time Shift/Single Run surrogate R_{n:m} values for all analyzed lengths, but outside the distribution of Random Permutation/Single Run surrogates for the longer time series (Figure 10B). We believe such a finding relates to what we observed for synthetic sawtooth waves, in which Random Permutation was more sensitive than Time Shift to detect the significance of the artifactual coupling caused by waveform asymmetry (Figure 4—figure supplement 2). In this sense, the n:m phaselocking between fundamental and harmonic frequencies would persist for small time shifts (±200 ms), albeit in different phase relations, while it would not resist the much larger time shifts obtained through random permutations. However, irrespective of this explanation, it should be noted that since the n:m phaselocking metrics cannot separate artifactual from true coupling, the possibility of the latter cannot be discarded. But if this is the case, we consider unlikely that the very low coupling level (~0.03) would have any physiological significance.
We conclude that the diagonal stripes in phasephase plots of both white noise and actual LFPs are mainly caused by a temporary n:m alignment of the phase timeseries secondary to the filteringinduced sinusoidality, and as such they are also apparent in surrogate data (Figure 8—figure supplement 1 and Figure 10—figure supplement 1). However, for actual LFPs there is a second influence, which can only be detected when analyzing very long epoch lengths, and which we believe is due to theta harmonics.
Discussion
Theta and gamma oscillations are hallmarks of hippocampal activity during active exploration and REM sleep (Buzsáki et al., 2003; Csicsvari et al., 2003; Montgomery et al., 2008). Theta and gamma are well known to interact by means of phaseamplitude coupling, in which the instantaneous gamma amplitude waxes and wanes as a function of theta phase (Bragin et al., 1995; SchefferTeixeira et al., 2012; Caixeta et al., 2013). This particular type of crossfrequency coupling has been receiving large attention and related to functional roles (Canolty and Knight, 2010; Hyafil et al., 2015). In addition to phaseamplitude coupling, theta and gamma oscillations can potentially interact in many other ways (Jensen and Colgin, 2007; Hyafil et al., 2015). For example, the power of slow gamma oscillations may be inversely related to theta power (Tort et al., 2008), suggesting amplitudeamplitude coupling. Recently, it has been reported that theta and gamma in hippocampal LFPs would also couple by means of n:m phaselocking (Belluscio et al., 2012; Zheng and Zhang, 2013; Xu et al., 2013, 2015; Zheng et al., 2016). Among other implications, this finding was taken as evidence for network models of working memory (Lisman and Idiart, 1995; Jensen and Lisman, 2005; Lisman, 2005) and for a role of basket cells in generating crossfrequency coupling (Belluscio et al., 2012; Buzsáki and Wang, 2012). However, our results show a lack of convincing evidence for n:m phaselocking in the two hippocampal datasets analyzed here, and further suggest that previous work may have spuriously detected phasephase coupling due to an improper use of surrogate methods, a concern also raised for phaseamplitude coupling (Aru et al., 2015).
Statistical inference of phasephase coupling
When searching for phasephase coupling between theta and gamma, we noticed that our R_{n:m} values differed from those reported in previous studies (Belluscio et al., 2012; Xu et al., 2013, 2015; Zheng et al., 2016). We suspected that this could be due to differences in the duration of the analyzed epochs. We then investigated the dependence of R_{n:m} on epoch length, and found a strong positive bias for shorter epochs. In addition, R_{n:m} values exhibit greater variability across samples as epoch length decreases for both white noise and actual data (e.g., compare in Figure 5B the data dispersion in Original R_{1:5} or R_{1:8} boxplots for different epoch lengths). Since theta and gamma peak frequencies are not constant in these signals, the longer the epoch, the more the theta and gamma peak frequencies are allowed to fluctuate and the more apparent the lack of coupling. On the other hand, $\mathrm{\Delta}{\phi}_{nm}$ distribution becomes less uniform for shorter epochs. The dependence of n:m phasecoupling metrics on epoch length has important implications in designing surrogate epochs for testing the statistical significance of actual R_{n:m} values. Of note, methodological studies on 1:1 phasesynchrony have properly used single surrogate runs of the same length as the original signal (Le Van Quyen et al., 2001; Hurtado et al., 2004). As demonstrated here, spurious detection of phasephase coupling may occur if surrogate epochs are longer than the original epoch. This is the case when one lumps together several surrogate epochs before computing R_{n:m}. When employing proper controls, our results show that R_{n:m} values of real data do not differ from surrogate values in theta epochs of up to 100 s. Moreover, the prominent bump in the R_{n:m} curve disappears when subsampling data at a lower frequency than gamma for both white noise and hippocampal LFPs (see Figure 2—figure supplement 1 and Figure 5—figure supplement 7), which suggests that it is due to the statistical dependence among contiguous data points introduced by the filter (which we referred to as ‘filteringinduced sinusoidality’).
Therefore, even though the n:m phaselocking metric R_{n:m} is theoretically welldefined and varies between 0 and 1, an estimated R_{n:m} value in isolation does not inform if two oscillations exhibit true phasecoupling or not. This can only be inferred after testing the statistical significance of the estimated R_{n:m} value against a proper surrogate distribution (but notice that falsepositive cases may occur due to waveform asymmetry; Figure 4C). While constructing surrogate data renders the metric computationally more expensive, such an issue is not specific for measuring n:m phaselocking but also happens for other metrics commonly used in the analysis of neurophysiological data, such as coherence, spikefield coupling, phaseamplitude coupling, mutual information and directionality measures, among many others (Le Van Quyen et al., 2001; Hurtado et al., 2004; Pereda et al., 2005; Tort et al., 2010).
The recent studies assessing thetagamma phasephase coupling in hippocampal LFPs have not tested the significance of individual R_{n:m} values against chance (Belluscio et al., 2012; Zheng and Zhang, 2013; Xu et al., 2013, 2015; Zheng et al., 2016). Two studies (Belluscio et al., 2012; Zheng et al., 2016) statistically inferred the existence of n:m phaselocking by comparing empirical phasephase plots with those obtained from the average of 1000 timeshifted surrogate runs. Specifically, Belluscio et al. (2012) established a significance threshold for each phasephase bin based on the mean and standard deviation of individual surrogate counts in that bin, and showed that the bin counts leading to diagonal stripes were statistically significant. Here we were able to replicate these results (Figure 7—figure supplement 1 and Figure 8). However, we note that a phasephase bin count is not a metric of n:m phaselocking; it does not inform coupling strength and even coupled oscillators have bins with nonsignificant counts. A bin count would be analogous to a phase difference vector ($e}^{i\mathrm{\Delta}{\phi}_{nm}\left(t\right)$), which is also not a metric of n:m phaselocking per se, but used to compute one. That is, in the same way that the R_{n:m} considers all phase difference vectors, n:m phaselocking can only be inferred when considering all bin counts in a phasephase plot. In this sense, by analyzing the phasephase plot as a whole, it was assumed that the appearance diagonal stripes was due to thetagamma coupling; no such stripes were apparent in phasephase plots constructed from the average over all surrogate runs (see Figure 6A in Belluscio et al., 2012). However, here we showed that single timeshifted surrogate runs do exhibit diagonal stripes (Figure 8—figure supplement 1 and Figure 10—figure supplement 1), that is, similar stripes exist at the level of a Single Run surrogate analysis, in the same way that Single Run surrogates also exhibit a bump in the R_{n:m} curve. Averaging 1000 surrogate phasephase plots destroys the diagonal stripes since different time shifts lead to different phase lags. Moreover, since the average is the sum divided by a scaling factor (the sample size), computing the average phasephase plot is equivalent to computing a single phasephase plot using the pool of all surrogate runs, which is akin to the issue of computing a single R_{n:m} value from a pooled surrogate distribution (Figure 2). Note that even bin counts in phasephase plots of white noise are considered significant under the statistical analysis introduced in Belluscio et al. (2012) (Figure 9A). Nevertheless, this was no longer the case when adapting their original framework to include a HolmBonferroni correction for multiple comparisons (Figure 9A).
Here we showed that the presence of diagonal stripes in phasephase plots is not sufficient to conclude the existence of phasephase coupling. The diagonal stripes are simply a visual manifestation of a maintained phase relationship, and as such they essentially reflect what R_{n:m} measures: that is, the ‘clearer’ the stripes, the higher the R_{n:m}. Therefore, in addition to true coupling, the same confounding factors that influence R_{n:m} also influence phasephase plots, such as filteringinduced sinusoidality and frequency harmonics. Our results suggest that the former is a main factor, because whitenoise signals have no harmonics but nevertheless display stripes in phasephase plots (Figure 9A). In accordance, no stripes are observed in phasephase plots of white noise when subsampling the time series (Figure 10—figure supplement 2; see also Figure 2—figure supplement 1). However, in actual LFPs filtering is not the only influence: (1) for the same filtered gamma band (30–50 Hz), the number of stripes relates to theta frequency (Figure 8); (2) for very long time series (i.e., 10–20 min of concatenated data), the stripes in phasephase plots of actual data – but not of white noise – persist after correcting for multiple comparisons when employing Random Permutation/Single Run surrogates (Figure 10A); (3) a stripedlike pattern remains in phasephase plots of actual LFPs after subsampling the time series (Figure 10—figure supplement 2). Consistently, R_{n:m} values of actual LFPs are greater than those of white noise in 1200 s epochs (~0.03 vs ~0.005, compare the bottom right panels of Figures 9B and 10B). Interestingly, Original R_{n:m} values of actual LFPs are not statistically different from the distribution of Time Shift/Single Run surrogates even for the very long epochs (Figure 10B), which suggests that Random Permutation is more powerful than Time Shift and should therefore be preferred. Though a very weak but true coupling effect cannot be discarded, based on our analysis of sawtooth waves (Figure 4 and Figure 4—figure supplement 2), we believe these results can be explained by theta harmonics, which would remain phaselocked to the fundamental frequency under small time shifts. Sharp signal deflections have been previously recognized to generate artifactual phaseamplitude coupling (Kramer et al., 2008; SchefferTeixeira et al., 2013; Tort et al., 2013; Aru et al., 2015; LozanoSoldevilla et al., 2016). Interestingly, Hyafil (2015) recently suggested that the nonsinusoidality of alpha waves could underlie the 1:2 phaselocking between alpha and beta observed in human EEG (Nikulin and Brismar, 2006; see also Palva et al., 2005). To the best of our knowledge, there is currently no metric capable of automatically distinguishing true crossfrequency coupling from waveforminduced artifacts in collective signals such as LFP, EEG and MEG signals. Ideally, learning how the signal is generated from the activity of different neuronal populations would answer whether true crossfrequency coupling exists or not (Hyafil et al., 2015), but unfortunately this is methodologically challenging.
Lack of evidence vs evidence of nonexistence
One could argue that we did not analyze a proper dataset, or else that prominent phasephase coupling would only occur during certain behavioral states not investigated here. We disagree with these arguments for the following reasons: (1) we could reproduce our results using a second dataset from an independent laboratory (Figure 5—figure supplement 3), and (2) we examined the same behavioral states in which n:m phaselocking was reported to occur (active waking and REM sleep). One could also argue that there exists multiple gammas, and that different gamma types are most prominent in different hippocampal layers (Colgin et al., 2009; SchefferTeixeira et al., 2012; Tort et al., 2013; Schomburg et al., 2014; Lasztóczi and Klausberger, 2014); therefore, prominent thetagamma phasephase coupling could exist in other hippocampal layers not investigated here. We also disagree with this possibility because: (1) we examined the same hippocampal layer in which thetagamma phasephase coupling was reported to occur (Belluscio et al., 2012); moreover, (2) we found similar results in all hippocampal layers (we recorded LFPs using 16channel silicon probes, see Materials and methods) (Figure 5—figure supplement 4) and (3) in parietal and entorhinal cortex recordings (not shown). Furthermore, similar results hold when (4) filtering LFPs within any gamma subband (Figure 5 and Figure 5—figure supplement 1 to 6), (5) analyzing CSD signals (Figure 5—figure supplement 4), or (6) analyzing independent components that maximize activity within particular gamma subbands (Schomburg et al., 2014) (Figure 5—figure supplement 5). Finally, one could argue that gamma oscillations are not continuous but transient, and that assessing phasephase coupling between theta and transient gamma bursts would require a different type of analysis than employed here. Regarding this argument, we once again stress that we used the exact same methodology as originally used to detect thetagamma phasephase coupling (Belluscio et al., 2012). Nevertheless, we also ran analysis only taking into account periods in which gamma amplitude was >2 SD above the mean (‘gamma bursts’) and found no statistically significant phasephase coupling (Figure 5—figure supplement 6).
Following Belluscio et al. (2012), other studies also reported thetagamma phasephase coupling in the rodent hippocampus (Zheng and Zhang, 2013; Xu et al., 2013, 2015; Zheng et al., 2016) and amygdala (Stujenske et al., 2014). In addition, human studies had previously reported thetagamma phasephase coupling in scalp EEG (Sauseng et al., 2008, 2009; Holz et al., 2010). Most of these studies, however, have not tested the statistical significance of coupling levels against chance (Sauseng et al., 2008, 2009; Holz et al., 2010; Zheng and Zhang, 2013; Xu et al., 2013, 2015; Stujenske et al., 2014), while Zheng et al. (2016) based their statistical inferences on the inspection of diagonal stripes in phasephase plots as originally introduced in Belluscio et al. (2012). We further note that epoch length was often not informed in the animal studies. Based on our results, we believe that differences in analyzed epoch length are likely to explain the high variability of R_{n:m} values across different studies, from ~0.4 (Zheng et al., 2016) down to 0.02 (Xu et al., 2013).
Since it is philosophically impossible to prove the absence of an effect, the burden of proof should be placed on demonstrating that a true effect exists. In this sense, and to the best of our knowledge, none of previous research investigating thetagamma phasephase coupling has properly tested R_{n:m} against chance. Many studies have focused on comparing changes in n:m phaselocking levels, but we believe these can be influenced by other variables such as changes in power, which affect the signaltonoise ratio and consequently also the estimation of the phase time series. Interestingly, in their pioneer work, Tass and colleagues used filtered white noise to construct surrogate distributions and did not find significant n:m phaselocking among brain oscillations (Tass et al., 1998, 2003). On the other hand, it is theoretically possible that n:m phaselocking exists but can only be detected by other types of metrics yet to be devised. In any case, our work shows that there is currently no convincing evidence for genuine thetagamma phasephase coupling using the same phaselocking metric (R_{n:m}) as employed in previous studies (Belluscio et al., 2012; Zheng and Zhang, 2013; Xu et al., 2013, 2015; Stujenske et al., 2014; Zheng et al., 2016), at least when examining LFP epochs of up to 100 s of prominent theta activity. For longer epoch lengths, though, we did find that R_{n:m} values of hippocampal LFPs may actually differ from those of randomly permuted, but not timeshifted, surrogates (Figure 10B). While we tend to ascribe such result to the effect of theta harmonics, we note that the possibility of true coupling cannot be discarded. But we are particularly skeptical that the very low levels of coupling strength observed in long LFP epochs would be physiologically meaningful.
Implications for models of neural coding by thetagamma coupling
Lisman and Idiart (1995) proposed an influential model in which theta and gamma oscillations would interact to produce a neural code. The thetagamma coding model has since been improved (Jensen and Lisman, 2005; Lisman, 2005; Lisman and Buzsáki, 2008), but its essence remains the same (Lisman and Jensen, 2013): nested gamma cycles would constitute memory slots, which are parsed at each theta cycle. Accordingly, Lisman and Idiart (1995) hypothesized that working memory capacity (7 ± 2) is determined by the number of gamma cycles per theta cycle.
Both phaseamplitude and phasephase coupling between theta and gamma have been considered experimental evidence for such coding scheme (Lisman and Buzsáki, 2008; Sauseng et al., 2009; Axmacher et al., 2010; Belluscio et al., 2012; Lisman and Jensen, 2013; Hyafil et al., 2015; Rajji et al., 2016). In the case of phaseamplitude coupling, the modulation of gamma amplitude within theta cycles would instruct a reader network when the string of items represented in different gamma cycles starts and terminates. On the other hand, the precise ordering of gamma cycles within theta cycles that is consistent across theta cycles would imply phasephase coupling; indeed, n:m phaselocking is a main feature of computational models of sequence coding by thetagamma coupling (Lisman and Idiart, 1995; Jensen and Lisman, 1996; Jensen et al., 1996). In contrast to these models, however, our results show that the theta phases in which gamma cycles begin/end are not fixed across theta cycles, which is to say that gamma cycles are not precisely timed but rather drift; in other words, gamma is not a clock (Burns et al., 2011).
If thetagamma neural coding exists, our results suggest that the precise location of gamma memory slots within a theta cycle is not required for such a code, and that the ordering of the represented items would be more important than the exact spike timing of the cell assemblies that represent the items (Lisman and Jensen, 2013).
Conclusion
In summary, while absence of evidence is not evidence of absence, our results challenge the hypothesis that thetagamma phasephase coupling exists in the hippocampus. At best, we only found significant R_{n:m} values when examining long LFP epochs (>100 s), but these had very low magnitude (and we particularly attribute their statistical significance to the effects of harmonics). We believe that the evidence in favor of n:m phaselocking in other brain regions and signals could potentially also be explained by simpler effects (e.g., filteringinduced sinusoidality, asymmetrical waveform, and improper statistical tests). While no current technique can differentiate spurious from true phasephase coupling, previous findings should be revisited and, whenever suitable, checked against the confounding factors and the more conservative surrogate procedures outlined here.
Materials and methods
Animals and surgery
Request a detailed protocolAll procedures were approved by our local institutional ethics committee (Comissão de Ética no Uso de Animais  CEUA/UFRN, protocol number 060/2011) and were in accordance with the National Institutes of Health guidelines. We used seven male Wistar rats (2–3 months; 300–400 g) from our breeding colony, kept under 12 hr/12 hr darklight cycle. We recorded from the dorsal hippocampus through either multisite linear probes (n = 6 animals; 4 probes had 16 4320 μm^{2} contacts spaced by 100 μm; 1 probe had 16 703 μm^{2} contacts spaced by 100 μm; 1 probe had 16 177 μm^{2} contacts spaced by 50 μm; all probes from NeuroNexus) or single wires (n = 1 animal; 50 μm diameter) inserted at AP −3.6 mm and ML 2.5 mm. Results shown in the main figures were obtained for LFP recordings from the CA1 pyramidal cell layer, identified by depth coordinate and characteristic electrophysiological benchmarks such as highest ripple power (see Figure 5—figure supplement 4 for an example). Similar results were obtained for recordings from other hippocampal layers (Figure 5—figure supplement 4).
We also analyzed data from three additional rats downloaded from the Collaborative Research in Computational Neuroscience data sharing website (www.crcns.org) (Figure 5—figure supplement 3). These recordings are a generous contribution by György Buzsáki’s laboratory (HC3 dataset, Mizuseki et al., 2013, 2014).
Data collection
Request a detailed protocolRecording sessions were performed in an open field (1 m x 1 m) and lasted 4–5 hr. Raw signals were amplified (200x), filtered between 1 Hz and 7.5 kHz (third order Butterworth filter), and digitized at 25 kHz (RHA2116, IntanTech). The LFP was obtained by further filtering between 1–500 Hz and downsampling to 1000 Hz.
Data analysis
Request a detailed protocolActive waking and REM sleep periods were identified from spectral content (high theta/delta power ratio) and video recordings (movements during active waking; clear sleep posture and preceding slowwave sleep for REM). The results were identical for active waking and REM epochs; throughout this work we only show the latter. The analyzed REM sleep dataset is available at http://dx.doi.org/10.5061/dryad.12t21. MATLAB codes for reproducing our analyses are available at https://github.com/tortlab/phase_phase .
We used builtin and customwritten MATLAB routines. Bandpass filtering was obtained using a least squares finite impulse response (FIR) filter by means of the ‘eegfilt’ function from the EEGLAB Toolbox (Delorme and Makeig, 2004). The filter order was three times the sampling rate divided by the low cutoff frequency. The eegfilt function calls the MATLAB ‘filtfilt’ function, which applies the filter forward and then again backwards to ensure no distortion of phase values. Similar results were obtained when employing other types of filters (Figure 5—figure supplement 8).
The phase time series was estimated through the Hilbert transform. To estimate the instantaneous theta phase of actual data, we filtered the LFP between 4–20 Hz, a bandwidth large enough to capture theta wave asymmetry (Belluscio et al., 2012). Estimating theta phase by the interpolation method described in Belluscio et al. (2012) led to similar results (Figure 5—figure supplement 2).
The CSD signals analyzed in Figure 5—figure supplement 4 were obtained as −A +2B −C, where A, B and C denote LFP signals recorded from adjacent probe sites. In Figure 5—figure supplement 5, the independent components were obtained as described in Schomburg et al. (2014); phaseamplitude comodulograms were computed as described in Tort et al. (2010).
n:m phaselocking
Request a detailed protocolWe measured the consistency of the phase difference between accelerated time series ($\mathrm{\Delta}{\phi}_{nm}\left({t}_{j}\right)=n{\ast \phi}_{\gamma}\left({t}_{j}\right)m\ast {\phi}_{\theta}\left({t}_{j}\right)$). To that end, we created unitary vectors whose angle is the instantaneous phase difference ($e}^{i\mathrm{\Delta}{\phi}_{nm}\left({t}_{j}\right)$), where j indexes the time sample, and then computed the length of the mean vector: ${R}_{n:m}=\Vert \frac{1}{N}{\sum}_{j=1}^{N}{e}^{i\mathrm{\Delta}{\phi}_{nm}\left({t}_{j}\right)}\Vert$, where N is the total number of time samples (epoch length in seconds x sampling frequency in Hz). R_{n:m} equals 1 when $\mathrm{\Delta}{\phi}_{nm}$ is constant for all time samples t_{j}, and 0 when $\mathrm{\Delta}{\phi}_{nm}$ is uniformly distributed. This metric is also commonly referred to as ‘mean resultant length’ or ‘mean radial distance’ (Belluscio et al., 2012; Stujenske et al., 2014; Zheng et al., 2016). Qualitatively similar results were obtained when employing the framework introduced in Sauseng et al. (2009), which computes the mean radial distance using gamma phases in separated theta phase bins, or the pairwise phase consistency metric described in Vinck et al. (2010) (Figure 5—figure supplement 1). Phasephase plots were obtained by first binning theta and gamma phases into 120 bins and next constructing 2D histograms of phase counts, which were smoothed using a Gaussian kernel of σ = 10 bins.
Surrogates
Request a detailed protocolIn all cases, theta phase was kept intact while gamma phase was mocked in three different ways: (1) Time Shift: the gamma phase time series is randomly shifted between 1 and 200 ms; (2) Random Permutation: a contiguous gamma phase time series of the same length as the original is randomly extracted from the same session. (3) Phase Scrambling: the timestamps of the gamma phase time series are randomly shuffled (thus not preserving phase continuity). For each case, R_{n:m} values were computed using either $\mathrm{\Delta}{\phi}_{nm}$ distribution for single surrogate runs (Single Run Distribution) or the pooled distribution of $\mathrm{\Delta}{\phi}_{nm}$ over 100 surrogate runs (Pooled Distribution).
For each animal, behavioral state (active waking or REM sleep) and epoch length, we computed 300 Original R_{n:m} values using different time windows along with 300 mock R_{n:m} values per surrogate method. Therefore, in all figures each boxplot was constructed using the same number of samples (=300 x number of animals). For instance, in Figure 5B we used n = 7 animals x 300 samples per animal = 2100 samples (but see Statistics below). In Figure 2, boxplot distributions for the whitenoise data were constructed using n = 2100.
Simulations
Request a detailed protocolKuramoto oscillators displaying n:m phaselocking were modeled as described in Osipov et al. (2007):
where $\epsilon $ is the coupling strength and ${\omega}_{\theta}$ and ${\omega}_{\gamma}$ are the natural frequencies of theta and gamma, respectively, which followed a Gaussian probability (σ = 5 Hz) at each time step. We used $\epsilon $ = 10, n = 1, m = 5, and dt = 0.001 s. The mean theta and gamma frequencies of each simulation are stated in the main text. For uncoupled oscillators, we set$\epsilon $ = 0.
For implementing the OI cell network (Figure 3B), we simulated the model previously described in Kopell et al. (2010). We used the same parameters as in Figure 3A of Kopell et al. (2010), with white noise (σ = 0.001) added to the I cell drive to create variations in spike frequency. NEURON (https://www.neuron.yale.edu/) codes for the model are available at ModelDB (https://senselab.med.yale.edu/).
The sawtooth wave in Figure 4C was simulated using dt = 0.001 s. Its instantaneous frequency followed a Gaussian distribution with mean = 8 Hz and σ = 5 Hz; white noise (σ = 0.1) was added to the signal.
In Figures 3 and 4C, boxplot distributions for simulated data were constructed using n = 300.
Statistics
For white noise data (Figure 2F), given the large sample size (n = 2100) and independence among samples, we used oneway ANOVA with Bonferroni posthoc test. For statistical analysis of real data (Figure 5B), we avoided nested design and inflation of power and used the mean R_{n:m} value per animal. In this case, due to the reduced sample size (n = 7) and lack of evidence of normal distribution (ShapiroWilk normality test), we used the Friedman’s test and Nemenyi posthoc test. In Figures 3 and 4C, we tested if R_{n:m} values of simulated data were greater than the distribution of surrogate values using onetailed ttests.
References

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

3
Crossfrequency phasephase coupling between θ and γ oscillations in the hippocampusJournal of Neuroscience 32:423–435.https://doi.org/10.1523/JNEUROSCI.412211.2012

4
Gamma (40100 Hz) oscillation in the hippocampus of the behaving ratJournal of Neuroscience 15:47–60.

5
Is gammaband activity in the local field potential of V1 cortex a "clock" or filtered noise?Journal of Neuroscience 31:9658–9664.https://doi.org/10.1523/JNEUROSCI.066011.2011
 6

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

8
Ketamine alters oscillatory coupling in the hippocampusScientific Reports 3:2348.https://doi.org/10.1038/srep02348
 9

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

14
EEGLAB: an open source toolbox for analysis of singletrial EEG dynamics including independent component analysisJournal of Neuroscience Methods 134:9–21.https://doi.org/10.1016/j.jneumeth.2003.10.009

15
The role of phase synchronization in memory processesNature Reviews Neuroscience 12:105–118.https://doi.org/10.1038/nrn2979

16
A mechanism for cognitive dynamics: neuronal communication through neuronal coherenceTrends in Cognitive Sciences 9:474–480.https://doi.org/10.1016/j.tics.2005.08.011
 17
 18
 19

20
Statistical method for detection of phaselocking episodes in neural oscillationsJournal of Neurophysiology 91:1883–1898.https://doi.org/10.1152/jn.00853.2003

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

22
Misidentifications of specific forms of crossfrequency coupling: three warningsFrontiers in Neuroscience 9:370.https://doi.org/10.3389/fnins.2015.00370
 23
 24

25
Hippocampal sequenceencoding driven by a cortical multiitem working memory bufferTrends in Neurosciences 28:67–72.https://doi.org/10.1016/j.tins.2004.12.001

26
Crossfrequency coupling between neuronal oscillationsTrends in Cognitive Sciences 11:267–269.https://doi.org/10.1016/j.tics.2007.05.003

27
Gamma and theta rhythms in biophysical models of hippocampal circuitsIn: V Cutsuridis, B Graham, S Cobb, I Vida, editors. Hippocampal Microcircuits: A Computational Modeler S Resource Book, 5. New York: SpringerVerlag. pp. 423–457.https://doi.org/10.1007/9781441909961_15

28
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

29
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
 30

31
Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchronyJournal of Neuroscience Methods 111:83–98.https://doi.org/10.1016/S01650270(01)003727
 32
 33

34
A neural coding scheme formed by the combined function of gamma and theta oscillationsSchizophrenia Bulletin 34:974–980.https://doi.org/10.1093/schbul/sbn060
 35

36
Neuronal oscillations with nonsinusoidal morphology produce spurious phasetoamplitude coupling and directionalityFrontiers in Computational Neuroscience 10:87.https://doi.org/10.3389/fncom.2016.00087
 37
 38
 39

40
Theta and gamma coordination of hippocampal networks during waking and rapid eye movement sleepJournal of Neuroscience 28:6731–6741.https://doi.org/10.1523/JNEUROSCI.122708.2008
 41
 42

43
Synchronization in Oscillatory NetworksSpringer Science & Business Media.https://doi.org/10.1007/9783540712695

44
Phase synchrony among neuronal oscillations in the human cortexJournal of Neuroscience 25:3962–3972.https://doi.org/10.1523/JNEUROSCI.425004.2005

45
Nonlinear multivariate analysis of neurophysiological signalsProgress in Neurobiology 77:1–37.https://doi.org/10.1016/j.pneurobio.2005.10.003

46
Ordering information in working memory and modulation of gamma by theta oscillations in humansCerebral Cortex, 10.1093/cercor/bhv326, 26759480.
 47

48
Brain oscillatory substrates of visual shortterm memory capacityCurrent Biology 19:1846–1852.https://doi.org/10.1016/j.cub.2009.08.062
 49
 50

51
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
 52

53
Movement enhances the nonlinearity of hippocampal thetaJournal of Neuroscience 36:4218–4230.https://doi.org/10.1523/JNEUROSCI.356415.2016
 54

55
Detection of N:m phase locking from noisy data: application to magnetoencephalographyPhysical Review Letters 81:3291–3294.https://doi.org/10.1103/physrevlett.81.3291
 56
 57
 58

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

61
The brainweb: phase synchronization and largescale integrationNature reviews. Neuroscience 2:229–239.https://doi.org/10.1038/35067550
 62
 63
 64
 65
 66
Decision letter

Frances K SkinnerReviewing Editor; University Health Network, Canada
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
[Editors’ note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration. The first decision letter after peer review is shown below.]
Thank you for submitting your work entitled "Lack of evidence for crossfrequency phasephase coupling between theta and gamma oscillations in the hippocampus" for consideration by eLife. Your article has been favorably evaluated by Timothy Behrens (Senior Editor) and four reviewers, one of whom is a member of our Board of Reviewing Editors.
Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews, we regret to inform you that your work will not be considered further for publication in eLife at this time. Below we provide detailed feedback for your consideration if you choose to resubmit, and we encourage this possibility.
While all of the reviewers felt that this was an important and interesting paper, they also thought that it was presented in a confusing manner, was limited, and could lead to confusion rather than clarification for the field. Given the contrarian nature of the manuscript as presented, this would not be helpful. The reviewers' overall comments are provided here:
1) The paper presentation should be changed to be less confrontational and completely clear on what is being said and done. A wellcrafted manuscript with clarity and sufficient explanations is lacking in the present submission.
A) Specifically, a more detailed review of the literature, a more careful presentation and discussion of the simulation results, and a more careful comparison with the procedures utilized in [Belluscio et al. 2012] is needed.
B) We note that phasecoupling is well defined theoretically from the tools of dynamical systems, but what is not clear is a selective measure of crossfrequency phase coupling. Defining an improved phasephase coupling measure that detects true coupling and ignores artifactual coupling might be possible. Such a measure – even if it's not the optimal measure – might serve as a "patch" to existing approaches.
C) Finding the optimal phasephase coupling measure (and assessing its statistical properties) would be a challenge. That is, a consistent and operational way to define how crossfrequency phase coupling can be measured (such that white noise or asymmetrical oscillation would not qualify). Clearly stating the lack of a clear measure as an open problem formally would be of value. In essence, the manuscript would need to make these aspects clear (i.e., present limitations of their work and others, possible fixes/challenges associated with having phasephase coupling measures), and in this way can identify a methodological issue that's driving incorrect conclusions in the literature, and so be of service to the field.
2) Given that a clear measure of crossfrequency phasecoupling does not exist, a change of title should be considered.
3) A repository (of data and algorithms) should complement the work, and would be beneficial for the community moving forth as analyses, clear understanding and interpretation of artificial and biological data would be available.
[Editors’ note: what now follows is the decision letter after the authors submitted for further consideration.]
Thank you for submitting your article "On crossfrequency phasephase coupling between theta and gamma oscillations in the hippocampus" for consideration by eLife. Your article has been favorably evaluated by Timothy Behrens (Senior Editor) and four reviewers, one of whom is a member of our Board of Reviewing Editors.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
In this resubmitted paper on phasephase coupling, all of the reviewers thought that several concerns had been addressed and suggestions implemented (such as the title change). The reviewers agreed that this is an important article that highlights a methodological issue impacting our understanding of the brain's activity and function. The authors are commended for sharing their analysis methods (on GitHub) and data.
The authors address three points in their paper: a statistical issue (using appropriates surrogate to detect n:m coupling), a conceptual issue (whether n:m coupling can be dissociated from asymmetrical oscillations), and an experimental one (whether n:m coupling is present in hippocampus). The first and last points are well made. However, this is not the case for the second, and other aspects remain unclear as presented. The authors are encouraged to revise and edit their paper to clarify all of the following issues.
1) A clear distinction between true phase versus estimated phase needs to be made as it is confusing as presented. Specifically:
A) As the authors' state, the theoretical quantity of n:m phase locking is welldefined. However, the estimated quantity of n:m phase locking as deduced from noisy neural data is fraught with difficulty. The authors show here that simply applying the theoretical n:m phase locking equation to data produces estimates that – without appropriate statistical tests – lead to spurious results. This is an important conclusion. Future work might seek to improve or correct (somehow) the estimation of the theoretical quantity, or propose a new theoretical quantity that is more easily estimated.
B) The point in the Discussion (subsection “Statistical inference of phasephase coupling”, last paragraph) that true n:m coupling can be disentangled from asymmetrical oscillation by visual inspection of the LFP or EEG trace, precisely because they will look very much alike does not seem appropriate. This distinction may be hard to make at the macroscopic level (LFP/EEG) and could rather be searched for in how this signal is generated from the activity of distinct neural subpopulations. If the same gamma oscillations that are shown to phaselock to theta oscillations can be dissociated neurally from theta oscillations (and indeed the authors did a tremendous job in previous publications at showing that hippocampal theta and gamma rely on distinct interconnected networks), then the evidence is rather of favor of genuine n:m locking. In the converse case, if no population selectively engaged in the faster oscillation can be identified, I would say that the asymmetrical oscillation is the most plausible explanation.
2) More detail of why the stripes in phasephase plot can be produced by asymmetry of the theta waves needs to be provided.
It is not clear how the authors can conclude for the dataset at hand that stripes in phasephase plot are due to the asymmetry of the theta waves and not to genuine gamma activity (Results, eighth paragraph). The authors seem to somehow tie this possible caveat to the stripes in phasephase plots. However, stripes are just a visual inspection for maintained relationship between the two phases, just as (properly controlled) Rnm should measure.
3) The discussion about the caveat of asymmetrical oscillations is unsatisfying for the following reasons:
A) This point is made in the middle of the analysis of experimental data. Since this has nothing to do with the data at hand, in which no evidence for phase coupling was found, this obscures the message. A separate section should be devoted to the caveat of the asymmetrical oscillations, possibly before the analysis of experimental data.
B) The point should be clarified, by not only focussing on the presence of stripes for asymmetrical oscillations, but by explaining conceptually why this dissociation emerges. Possibly it could be illustrated by showing how an asymmetrical wave resembles the sum of multiples sinusoids at the fundamental and harmonics frequencies n:m locked to each other (and possibly with also some phaseamplitude coupling, see Kramer et al.).
4) Any differences between the measure of the authors and that of Belluscio's needs to be made clear. As specified by one of the reviewers:
It is not quite clear whether the methods used are exactly the same. In Belluscio et al., I think the counts were performed first for each surrogate (corresponding to a specific timeshift), and then the distribution of counts over the distinct timeshift was computed for each bin (i.e. Single Run analysis, which should be the correct method). By contrast, it seems that here – not quite clear in the text – surrogates for distinct time shifts were merged before counts were done in Figure 5 and Figure 4—figure supplement 1 (i.e. Pooled analysis, which of course averages the count and creates many false positives). [Note that I use Single Run vs. Pooled terminology to refer to whether the metrics was computed for each surrogate independently (Single Run) or for merged data, the metrics applying yo either Rnm or counts.] I hope that I am being wrong here. Otherwise that would seriously limit the conclusions of the author about invalidating Belluscio's findings.
Note that there is still a statistical problem with the Single Run analysis of phasephase counts, which relates to the absence of a correction for the number of tested bins (either for a given theta phase or across all theta phases). In any case, the pvalue in Belluscio et al. reaches very low value (some p<10^10), which would probably resist proper corrections. If this is the case, that would provide statistically sound evidence for n:m coupling (asymmetrical theta being still a possible confound). Do the authors observe such low values in their own dataset? Could they remove statistical significance in the stripes in white noise using an appropriate correction? Or does this deceiving statistical significance come from yet another explanation?
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "On crossfrequency phasephase coupling between theta and gamma oscillations in the hippocampus" for further consideration at eLife. Your revised article has been favorably evaluated by Timothy Behrens (Senior editor), a Reviewing editor, and three reviewers.
The manuscript has been improved but there are remaining issues that need to be addressed. All of the reviewers appreciated the authors' revisions and thought that this work could serve as a useful cautionary tale for neuroscientists as well as a useful starting point for continuing work. Even though it was felt that the results were not completely satisfying, it was also felt that it would be of great benefit in preventing researchers from wandering down the wrong analysis path. It was further noted that if this paper were presented as a statistics paper it could be improved more but then it would probably be less read and absorbed by the neuroscience community for some time.
An overall suggestion is that the authors should be more cautious in their interpretation and clearer about there being room for further progress in editing their manuscript.
Specific issues to address are:
1) The Abstract would benefit from significant revisions considering the various changes already done, and ones to be done. Sentences that should be targeted are: "filtered white noise has similar n:m phaselocking levels as actual data", and "the diagonal stripes in thetagamma phasephase histograms of actual data can be explained by theta harmonics".
2) Validity of Belluscio's analysis – overall conclusion:
A) In the rebuttal, the authors consider that the analysis by Belluscio is intrinsically flawed rather than being simply not well controlled statistically. I tend to defend the alternative option (clearly here I do not refer to the analysis of Mean PhasePhase Plot, which the authors have convincingly showed that is flawed, but to the tests on phase counts; I also want to reassure the authors that I really do not have any personal motivation to go one way or the other). My intuition is that the presence of one single significant count in the phasephase plot, if appropriately controlled for multiple comparisons, would provide a valid statistical measure. In essence it would not simply detect the presence of stripes (which surrogate run also feature) but measure whether their amplitude is larger than those of surrogates. This could be shown in a quite straightforward way by looking whether significant points in phasephase plots obtained from white noise persist when a correction for multiple comparisons is applied. The authors seem to agree that they would not. (As for the type of correction, the Bonferroni method that the authors refer to looks too conservative as the bin counts clearly are nonindependent; less stringent correction such as FDR or HolmBonferroni may be preferred). A count is not a metric for n:m coupling for sure, but it can inform of the particular concentration of the high frequency phase at one specific phase of the high frequency phase. More importantly, what matters here is not a metric of the coupling strength, but a reliable statistical test that selectively detects coupling, and my intuition is that counts in the phasephase plot may provide one. If a simple correction can be applied to give a sound statistical test, this could allow Belluscio and colleagues to look back at their data and see whether the significant counts indeed resist correction for multiple comparisons.
B) About the very low pvalues in Belluscio's data, the authors suggest that "it is possible that, by analyzing a longer epoch length, the influence of theta harmonics becomes more apparent and would lead to lower pvalues, while the effect of the filteringinduced sinusoidality is washed out for both actual and surrogate epochs (of the same length)."
Well the manuscript previously demonstrates that n:m coupling measures as well as phasephase plots cannot tear apart asymmetrical waves from true crossfrequency coupling. Thus, if indeed pvalues remain lower than threshold when controlled appropriately, it could equally be due to asymmetrical theta or true thetagamma phase coupling, but the latter could not be dismissed.
3) The reason for stripes in the phasephase plot of hippocampal data:
A) In the manuscript authors state that "the diagonal stripes in phasephase plots are due to theta harmonics and not to genuine gamma activity." If it were true harmonics, that would imply a consistent phase relationship between theta and gamma over sustained periods (not just over small periods as for white noise), i.e. would test positive for n:m coupling. So in my opinion here the effect is rather due to the temporary n:m alignment of phases, just as for white noise. In other words, stripes emerge also when analyzing white noise despite there being no harmonics in the signal.
B) In the rebuttal, the authors defend that "Believing in two different and genuine gamma activities – one coupled at 5 cycles per theta, the other at 4 cycles per theta – would go against the parsimonious principle." (this is when theta frequency evolves). However, as established by theoretical studies of coupled oscillators, n:m coupling is not a fixed property of the network but emerges as a combination of oscillator intrinsic dynamics and the characteristics of the connectivity pattern. See for example Arnold tongue: if we assume there is just one fixed gamma, a fluctuating theta and fixed connectivity, it is perfectly normal that n:m coupling will shift from 1:5 to 1:4 if the lower oscillator accelerates and the ratio of frequency goes from around 5 to around 4.
4) Asymmetrical theta vs. phasecoupled thetagamma in general
I am coming back to the authors' comment on this in the rebuttal, although this is no longer present in the manuscript. This is a comment for the authors benefit and speaks to the overall comment mentioned above. Sharp edges in LFP/EEG are not by itself an indication that we are measuring a single asymmetrical oscillation as the superposition of a slow and a fast n:m coupled oscillation can also give rise to sharp edges (Figure 4C; note that this can even be obtained with just two oscillations). Thus, sharp edges is no more a selective feature of asymmetrical oscillations than n:m coherence is selective of phasecoupled oscillators. In other words, visual inspection cannot tell us more than statistics.
5) It looks like "Random Perm" is more powerful test than "Time Shift", as mentioned in the first paragraph of the subsection “Lack of evidence for n:m phaselocking in actual LFPs”. If this is the case then it could be stated as a conclusion of the work that "Random perm" should be preferred over "Time Shift" (and of course of the "Scrambling" procedure), as it is less likely to miss existing effects (lower False Rejections rate).
6) Conclusion: the very new paper by LozanoSoldevilla and colleagues (Frontiers Comp Neuro) provides another example of spurious crossfrequency coupling measures due to asymmetrical oscillations, could be worth referencing.
https://doi.org/10.7554/eLife.20515.036Author response
[Editors’ note: the author responses to the first round of peer review follow.]
While all of the reviewers felt that this was an important and interesting paper, they also thought that it was presented in a confusing manner, was limited, and could lead to confusion rather than clarification for the field. Given the contrarian nature of the manuscript as presented, this would not be helpful. The reviewers' overall comments are provided here:
1) The paper presentation should be changed to be less confrontational and completely clear on what is being said and done. A wellcrafted manuscript with clarity and sufficient explanations is lacking in the present submission.
We have made several modifications in our manuscript to let it less confrontational, including a change in the title. We have also edited our Results section to clarify points raised by reviewer 3.
A) Specifically, a more detailed review of the literature, a more careful presentation and discussion of the simulation results, and a more careful comparison with the procedures utilized in [Belluscio et al. 2012] is needed.
We have included new references mentioned by the reviewers, and performed new analyses and simulations to address important comments by reviewer 3; this led to the elaboration of new figures and text. To further conform to the reviewer’s comments, we have deleted former Figure 1 and related text, and restructured the way we present and discuss our results. Finally, we now provide p value maps to allow straight comparison with the procedure performed in Belluscio et al. (2012).
B) We note that phasecoupling is well defined theoretically from the tools of dynamical systems, but what is not clear is a selective measure of crossfrequency phase coupling. Defining an improved phasephase coupling measure that detects true coupling and ignores artifactual coupling might be possible. Such a measure – even if it's not the optimal measure – might serve as a "patch" to existing approaches.
C) Finding the optimal phasephase coupling measure (and assessing its statistical properties) would be a challenge. That is, a consistent and operational way to define how crossfrequency phase coupling can be measured (such that white noise or asymmetrical oscillation would not qualify). Clearly stating the lack of a clear measure as an open problem formally would be of value. In essence, the manuscript would need to make these aspects clear (i.e., present limitations of their work and others, possible fixes/challenges associated with having phasephase coupling measures), and in this way can identify a methodological issue that's driving incorrect conclusions in the literature, and so be of service to the field.
We will answer to points B and C together:
The definition of crossfrequency phase coupling follows the exact same rationale as the definition of phasecoupling, and the n:m phasecoupling metric is the natural expansion of the 1:1 phase coupling metric. We believe both metrics are well defined. The current concept of n:m phaselocking reflected in the metric’s definition (i.e., the concentration of the phase difference distribution) is used not only by experimentalists but also by theoretical groups. Please note that we did not create the n:m phaselocking metric in this work but studied the exact same metric in current widespread use in the literature, which led influential conclusions about how the brain might work. In our paper, we identify the lack of proper control/surrogate analyses as the methodological issue driving incorrect conclusions, which should prompt other groups to revisittheir findings. We believe these results are of enough importance to many neuroscientists working with brain rhythms and their interactions (e.g., Schack et al., 2002; Schack et al., 2005; Schack and Weiss, 2005; Palva et al., 2005; Holz et al., 2010; Sauseng et al., 2008; Sauseng et al., 2009; Zheng et al., 2013; Xu et al., 2013; Stujenske et al., 2014; Xu et al., 2015; Aru et al., 2015; Hyafill et al., 2015; Chaieb et al., 2015; Zheng et al., 2016, among many others).
The points above seem to implicitly assume that n:m phasephase coupling would exist in the brain but the current metric is not able to detect it. While we have difficulties in imagining how this would exactly happen (since the metric and the definition of n:m phase coupling are intrinsically related), we do not disagree with such possibility and have edited our text to discuss it. In addition, we performed new simulations to show that cases of significant coupling may reflect not only true coupling but also artifactual coupling. We tried but were unable to devise a way to separate these cases; in the revised version, we have expanded our Results and Discussion sections to comment on these limitations.
But please note that our results point to lack of statistically significant coupling in actual LFPs upon proper surrogate procedures, and therefore such issues do not interfere with the main conclusion of our paper.
2) Given that a clear measure of crossfrequency phasecoupling does not exist, a change of title should be considered.
We have changed our title to “On crossfrequency phasephase coupling between theta and gamma oscillations in the hippocampus”. But please note that above and in our answers below we defend that there exists a current, welldefined measure of n:m phaselocking, and that the lack of evidence for n:m phasecoupling reported in our work may not be due to a “fault” of this metric but simply to lack of existence of such effect.
3) A repository (of data and algorithms) should complement the work, and would be beneficial for the community moving forth as analyses, clear understanding and interpretation of artificial and biological data would be available.
We have made our codes available at https://github.com/tortlab/phase_phase. In addition, we also made the analyzed data available at Dryad (http://datadryad.org/). If the manuscript is accepted for publication, the dataset will have the doi: 10.5061/dryad.12t21.
[Editors' note: the author responses to the rereview follow.]
In this resubmitted paper on phasephase coupling, all of the reviewers thought that several concerns had been addressed and suggestions implemented (such as the title change). The reviewers agreed that this is an important article that highlights a methodological issue impacting our understanding of the brain's activity and function. The authors are commended for sharing their analysis methods (on GitHub) and data.
The authors address three points in their paper: a statistical issue (using appropriates surrogate to detect n:m coupling), a conceptual issue (whether n:m coupling can be dissociated from asymmetrical oscillations), and an experimental one (whether n:m coupling is present in hippocampus). The first and last points are well made. However, this is not the case for the second, and other aspects remain unclear as presented. The authors are encouraged to revise and edit their paper to clarify all of the following issues.
1) A clear distinction between true phase versus estimated phase needs to be made as it is confusing as presented. Specifically:
A) As the authors' state, the theoretical quantity of n:m phase locking is welldefined. However, the estimated quantity of n:m phase locking as deduced from noisy neural data is fraught with difficulty. The authors show here that simply applying the theoretical n:m phase locking equation to data produces estimates that – without appropriate statistical tests – lead to spurious results. This is an important conclusion. Future work might seek to improve or correct (somehow) the estimation of the theoretical quantity, or propose a new theoretical quantity that is more easily estimated.
We understand the reviewer’s point: for any given signal the n:m phaselocking metric will provide an estimated value irrespective of the existence of true phasephase coupling or not; the former can only be inferred when comparing the estimated value with a chance/surrogate distribution. But we note that this issue is not specific for n:m phaselocking estimates, as it also occurs for other metrics commonly used in Neuroscience. For instance, estimates of mutual information, coherence, spikefield coupling, phaseamplitude coupling, directionally measures – to mention a few – all suffer from short sample bias. In the rerevised manuscript, we now devote a paragraph to discuss such issue and highlight the difference between estimated and true phasephase coupling (please see subsection “Statistical inference of phasephase coupling”, second paragraph).
B) The point in the Discussion (subsection “Statistical inference of phasephase coupling”, last paragraph) that true n:m coupling can be disentangled from asymmetrical oscillation by visual inspection of the LFP or EEG trace, precisely because they will look very much alike does not seem appropriate. This distinction may be hard to make at the macroscopic level (LFP/EEG) and could rather be searched for in how this signal is generated from the activity of distinct neural subpopulations. If the same gamma oscillations that are shown to phaselock to theta oscillations can be dissociated neurally from theta oscillations (and indeed the authors did a tremendous job in previous publications at showing that hippocampal theta and gamma rely on distinct interconnected networks), then the evidence is rather of favor of genuine n:m locking. In the converse case, if no population selectively engaged in the faster oscillation can be identified, I would say that the asymmetrical oscillation is the most plausible explanation.
The passage of our Discussion mentioned by the reviewer reads, former version):
“To the best of our knowledge, there is currently no metric capable of automatically distinguishing true crossfrequency coupling from waveforminduced artifacts in collective signals such as LFP, EEG and MEG signals. Until one is developed, we believe the best practice is to judge statistically significant coupling levels on a casebycase basis by visual inspection of raw signals.”
We meant that if theta and gamma oscillations genuinely n:m phasecouple, one should be able to observe this effect in the unfiltered LFP (i.e., one should see the actual gamma cycles coupled to theta upon visual inspection); on the other hand, if n:m coupling is spurious and due to an asymmetric theta wave, one would observe sharp deflections in the unfiltered LFP instead of genuine gamma cycles. We have previously discussed such “manual” solution in Kramer et al. (2008) for the cases of spurious phaseamplitude coupling. On page 356 of that paper, we wrote:
“To ameliorate the potential problem of spurious frequency comodulation, we recommend careful inspection of the unfiltered data to distinguish between sharp edges and true high frequency oscillations.” And later: “First, we recommend visual inspection of the unfiltered data at times corresponding to increases in the amplitude envelope of the high frequency activity. At these times, do the unfiltered data exhibit high frequency oscillations or sharp edges?”
But we acknowledge that a problem arises when the signal has both gamma oscillations and sharp deflections, in which case it is very difficult – if not impossible – to realize by visual inspection which effect is at play, and therefore we agree with the reviewer that the sentence may not be appropriate.
We further agree with the reviewer that understanding how the signal is generated by the activity of distinct neural subpopulations would be ideal to sort apart true from artifactual coupling. But this constitutes a huge problem in Neuroscience by itself: after all, what gives rise to LFP oscillations? After decades of research, we still debate about what generates theta and gamma at the circuitry and neuronal levels; we currently have models but not a complete understanding of oscillatory genesis; the available experimental data are not yet fully conclusive. Computer simulations – as the ones mentioned by the reviewer – are useful tools to generate insights into potential mechanisms, but they cannot be considered ultimate proofs. Thus, the gold standard solution referred to by the reviewer (of searching how the signal is generated by the activity of neural subpopulations) is unfortunately not a practical one but methodologically very challenging.
To address this issue, in the rerevised version we removed the sentence the reviewer did not find appropriate, and instead we now spell out the ideal suggestion put forward by the reviewer (subsection “Statistical inference of phasephase coupling”, last paragraph).
2) More detail of why the stripes in phasephase plot can be produced by asymmetry of the theta waves needs to be provided.
It is not clear how the authors can conclude for the dataset at hand that stripes in phasephase plot are due to the asymmetry of the theta waves and not to genuine gamma activity (Results, eighth paragraph). The authors seem to somehow tie this possible caveat to the stripes in phasephase plots. However, stripes are just a visual inspection for maintained relationship between the two phases, just as (properly controlled) Rnm should measure.
We agree with the reviewer that stripes in phasephase plots are only a visual manifestation of a maintained phase relationship, and in essence the stripes reflect what R_{n:m} measures. The stripes may be due to (1) theta harmonics (i.e., theta asymmetry), due to (2) genuinely coupled gamma activity, or even due to (3) the sinusoidality imposed by the filter. Similarly, all these 3 effects can lead to a bump in the R_{n:m} curve (see Figure 3 for genuine coupling, Figure 4 for theta asymmetry, and Figures 2 and Figure 2—figure supplement 1 for the filteringimposed sinusoidality in white noise).
Therefore, the reviewer is right in that it is not possible to differentiate genuinely coupled gamma from theta harmonics (or filtering effects) solely based on visual inspection of phasephase plots, and the same is also true for the visual inspection of a bump in the R_{n:m} curve. But please note that we did not mean to imply that the presence of stripes per se was a proof of theta harmonics and absence of genuine coupling. Rather, our conclusion of absence of genuine coupling stemmed from a conjunction of observations, along with the Occam's principle of parsimony:
1) As correctly pointed out by the reviewer, the R_{n:m} metric and the stripes are essentially related: the more prominent the stripes are in phasephase plots, the higher the R_{n:m}. Therefore, were the stripes due to genuinely coupled gamma activity, they should be associated with R_{n:m} values different from chance, which was not the case (Figure 5 and its figure supplements).
2) As shown in Figure 7, the number of stripes perfectly relates to the order of the first theta harmonics falling within the filtered bandwidth. That is, instead of assuming the existence of multiple types of gamma activities coupled to theta, we deemed more parsimonious to attribute the different number of stripes shown in Figure 7 to different orders of the theta harmonics.
3) Please notice in Figure 8 that the filter bandwidth for gamma is the same (3050 Hz). Yet, the phasephase plot of the LFP in 8A has 5 stripes while the phasephase plot of the LFP in 8B has 4 stripes. Instead of assuming that theta would couple to a genuine 1:5 gamma (i.e., 5 gamma cycles per theta cycle) and to a genuine 1:4 gamma (4 gamma cycles per theta cycle), we believe the stripes are straightforwardly explained by the fact that theta frequency in 8A example is 7.1 Hz while in 8B theta frequency is 8.4 Hz; therefore, it takes a factor of 5 and 4, respectively, to make these theta frequencies larger than the low cutoff frequency of the gamma filter (i.e., 7.1 Hz x 5 = 35.5 Hz > 30 Hz and 8.4 x 4 = 33.6 Hz > 30 Hz), which would explain why we see 5 and 4 stripes, respectively. Believing in two different and genuine gamma activities – one coupled at 5 cycles per theta, the other at 4 cycles per theta – would go against the parsimonious principle.
But while our observations argue against genuine gamma coupling, the distinction between the contribution of theta harmonics and of filteringinduced sinusoidality is trickier, and actually both effects may be at play. In brief, though, the effect of filteringinduced sinusoidality is short lived and would not account for the stripes seen in very long epochs such as the 20min long signal analyzed in Figure 7.
In the rerevised manuscript, we elaborated new text to clarify these points (subsection “Statistical inference of phasephase coupling”, fourth paragraph).
3) The discussion about the caveat of asymmetrical oscillations is unsatisfying for the following reasons:
A) This point is made in the middle of the analysis of experimental data. Since this has nothing to do with the data at hand, in which no evidence for phase coupling was found, this obscures the message. A separate section should be devoted to the caveat of the asymmetrical oscillations, possibly before the analysis of experimental data.
We thank the reviewer for the nice suggestion, which we have followed in the rerevised version. Namely, we now divide our Results section into subsections (in the former version we had a single text), and devote one subsection – before we present the experimental results – to show the caveat of asymmetrical/nonsinusoidal oscillations. This has led to major changes throughout the manuscript, such as the elaboration of a new figure and associated text (see below), and reordering of previous figures and text. For convenience, we submitted as “related manuscript file” a PDF version of our manuscript with all changes highlighted.
B) The point should be clarified, by not only focussing on the presence of stripes for asymmetrical oscillations, but by explaining conceptually why this dissociation emerges. Possibly it could be illustrated by showing how an asymmetrical wave resembles the sum of multiples sinusoids at the fundamental and harmonics frequencies n:m locked to each other (and possibly with also some phaseamplitude coupling, see Kramer et al.).
To address this point, we have elaborated a new figure in which we illustrate that the Fourier decomposition of a nonsinusoidal wave is characterized by the sum of multiple sinusoids at harmonic frequencies that are n:m phaselocked to each other; please see new Figure 4 and new subsection of the Results section (c.f. our answer above). We mention that an asymmetric wave shape may also lead to spurious phaseamplitude coupling in the first paragraph of the subsection “Spurious n:m phaselocking due to nonsinusoidal waveforms” and in the last paragraph of the subsection “Statistical inference of phasephase coupling”, and explicitly show this effect in new Figure 4—figure supplement 1.
4) Any differences between the measure of the authors and that of Belluscio's needs to be made clear. As specified by one of the reviewers:
It is not quite clear whether the methods used are exactly the same. In Belluscio et al., I think the counts were performed first for each surrogate (corresponding to a specific timeshift), and then the distribution of counts over the distinct timeshift was computed for each bin (i.e. Single Run analysis, which should be the correct method). By contrast, it seems that here – not quite clear in the text – surrogates for distinct time shifts were merged before counts were done in Figure 5 and Figure 4—figure supplement 1 (i.e. Pooled analysis, which of course averages the count and creates many false positives). [Note that I use Single Run vs. Pooled terminology to refer to whether the metrics was computed for each surrogate independently (Single Run) or for merged data, the metrics applying yo either Rnm or counts.] I hope that I am being wrong here. Otherwise that would seriously limit the conclusions of the author about invalidating Belluscio's findings.
The statistical analysis of the phasephase plots was done exactly as described in Belluscio et al. and correctly summarized by the reviewer. Namely, for each phasephase bin, each single surrogate run contributed with one count (please note that we did not merge single surrogates before counting). Then, for each phasephase bin we constructed a distribution of counts, and, based on the mean and standard deviation of this distribution, we inferred the statistical significance of the actual count for that bin in the phasephase plot.
In the rerevised version, we have edited the legend and panel labels of Figure 7, Figure 7—figure supplement 1, Figure 8, Figure 8—figure supplement 1 (current numeration) and subsection “Lack of evidence for n:m phaselocking in actual LFPs”, second paragraph and subsection “Statistical inference of phasephase coupling”, third paragraph to let it clearer that the mean and SD are computed over individual bin counts, and not over pooled counts.
Importantly, please note that we do not consider the framework employed in Belluscio et al. as a Single Run analysis. We reserve the term “Single Run” to denote an analysis in which each single surrogate run generates a single measure of n:m phaselocking. In this sense, we don’t see a single histogram count as a metric of n:m phaselocking, in the same way that we don’t see a single vector of phase difference (${e}^{i\ufffd{\text{\phi}}_{nm}\left(t\right)}$) as a metric of n:m phaselocking, simply because neither a single histogram count nor a single phase difference vector measures the level of n:m phaselocking between the time series. Rather, for each single surrogate run, a metric of n:m phaselocking is obtained by computing the length of the mean vector over all phase difference vectors (R_{n:m}). Therefore, sticking to this analogy, the counts in phasephase plots would correspond to the phase difference vectors used to compute the R_{n:m}, but not to the R_{n:m} metric itself. Again, a count is not an n:m phaselocking metric. In the case of the phasephase plot, a Single Run surrogate would be an entire phasephase plot constructed from a single surrogate (i.e., timeshifted) time series, because it is at this level that we infer the existence of n:m phaselocking – and not at the level of a single bin count. In other words, Belluscio et al. inferred the existence of n:m phaselocking by the visual inspection of stripes. And, as shown in Figure 8—figure supplement 1, phasephase plots from single surrogate runs do exhibit stripes similar to the original phasephase plot. This relates to the fact that single surrogate R_{n:m} values are similar to original R_{n:m} values, that is, there is no difference between actual and surrogate data at the level of an n:m phaselocking metric, which justifies our conclusions.
In the rerevised version, we substantially edited a previous Discussion paragraph to clarify these points (please see subsection “Statistical inference of phasephase coupling”, third paragraph).
Supplementary remarks: In the right panel of Figure 6A, Belluscio et al. show the mean phasephase plot over all surrogate runs. In our paper we note that plotting the mean number of counts across all surrogates is qualitatively identical to plotting the pool of counts over all surrogate runs, up to a scaling factor, as shown in Author response image 1.
Clearly, the scaling factor is the number of surrogates (n=1000). The equivalence above motivated us to state that the procedure in Belluscio et al. is akin to what we call as a pooled surrogate analysis: there is no visual inspection of phasephase plots for single surrogate runs but only the inspection of their average, which is equivalent to the visual inspection of their pool (of note, this improper visual analysis was also performed in Zheng et al. 2016). At any event, we stress that there is no proper Single Run surrogate analysis performed at the level of an n:m phaselocking metric such as the R_{n:m} in Belluscio et al.
Note that there is still a statistical problem with the Single Run analysis of phasephase counts, which relates to the absence of a correction for the number of tested bins (either for a given theta phase or across all theta phases). In any case, the pvalue in Belluscio et al. reaches very low value (some p<10^10), which would probably resist proper corrections. If this is the case, that would provide statistically sound evidence for n:m coupling (asymmetrical theta being still a possible confound). Do the authors observe such low values in their own dataset? Could they remove statistical significance in the stripes in white noise using an appropriate correction? Or does this deceiving statistical significance come from yet another explanation?
Indeed, in most phasephase bins the counts would not survive a Bonferroni correction. (Since we used 120x120 phasephase bins, a correction would bring the significance threshold down to 10^{6}; we don’t know how many phasephase bins Belluscio et al. employed and whether their counts would survive a Bonferroni correction). But please notice that this is actually one more argument against the existence of phasephase coupling. Forgetting for the moment the issues with this framework mentioned above, were the described effect robust, we should have been able to replicate it.
We suspect the pvalues in Belluscio et al. are lower than ours due to differences in epoch length; they stated they analyzed 11 sessions in three rats, but the total analyzed length is not informed. Intuitively, we believe it is possible that, by analyzing a longer epoch length, the influence of theta harmonics becomes more apparent and would lead to lower pvalues, while the effect of the filteringinduced sinusoidality is washed out for both actual and surrogate epochs (of the same length).
As for this statement: “If this is the case, that would provide statistically sound evidence for n:m coupling (asymmetrical theta being still a possible confound).”
As reviewers ourselves, we appreciate the voluntary and timeconsuming work of the reviewer, and thank him/her for the relevant suggestions that helped us improving our paper. But while we respect the reviewer’s opinion, we feel compelled to disagree with this particular statement. Based on arguments provided above, we do not think that a proper statistical procedure was performed in Belluscio et al. After all, what is the n:m phaselocking metric being measured in Belluscio et al. framework? A single phasephase bin count is not a metric, it does not inform how coupled the time series are, and even coupled oscillators have nonsignificant bin counts. Please note that Belluscio et al. state “The diagonal stripes in the phase–phase plots indicate a statistically reliable relationship between the respective phases of theta and gamma oscillations during both RUN and REM”.Thus, the “metric” in Belluscio et al. is in essence the visual inspection of stripes. But single run surrogates also exhibit stripes (Figure 8—figure supplement 1), in the same way that single run surrogates also exhibit a bump in the R_{n:m} curve. The stripes are trivially averaged out when one computes the mean over timeshifted phasephase plots, which relates to a pooled analysis (c.f. Author response image 1); it is not too surprising that actual counts of some phasephase bins are deemed statistically significant under this framework. Please note that a single surrogate phasephase plot – which has stripes – would also be deemed statistically significant if one were to apply this very same framework. Therefore, we do not consider this framework a proper one. Instead, Belluscio et al. should have inferred the statistical significance of coupling at the R_{n:m} level, which is an actual metric of n:m phaselocking.
In the rerevised version, we have substantially edited a previous Discussion paragraph to explain why an histogram count would not constitute an n:m phaselocking metric, and why we do not consider the analysis framework in Belluscio et al. a proper one. Please see subsection “Statistical inference of phasephase coupling”, third paragraph.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
[…] An overall suggestion is that the authors should be more cautious in their interpretation and clearer about there being room for further progress in editing their manuscript.
Specific issues to address are:
1) The Abstract would benefit from significant revisions considering the various changes already done, and ones to be done. Sentences that should be targeted are: "filtered white noise has similar n:m phaselocking levels as actual data", and "the diagonal stripes in thetagamma phasephase histograms of actual datacan be explained by theta harmonics".
We thank the reviewer for calling attention to this issue. In the revised version we have edited the Abstract to reflect changes done in the previous and current version. The two mentioned sentences no longer appear in the revised version.
2) Validity of Belluscio's analysis – overall conclusion:
A) In the rebuttal, the authors consider that the analysis by Belluscio is intrinsically flawed rather than being simply not well controlled statistically. I tend to defend the alternative option (clearly here I do not refer to the analysis of Mean PhasePhase Plot, which the authors have convincingly showed that is flawed, but to the tests on phase counts; I also want to reassure the authors that I really do not have any personal motivation to go one way or the other). My intuition is that the presence of one single significant count in the phasephase plot, if appropriately controlled for multiple comparisons, would provide a valid statistical measure. In essence it would not simply detect the presence of stripes (which surrogate run also feature) but measure whether their amplitude is larger than those of surrogates. This could be shown in a quite straightforward way by looking whether significant points in phasephase plots obtained from white noise persist when a correction for multiple comparisons is applied. The authors seem to agree that they would not. (As for the type of correction, the Bonferroni method that the authors refer to looks too conservative as the bin counts clearly are nonindependent; less stringent correction such as FDR or HolmBonferroni may be preferred). A count is not a metric for n:m coupling for sure, but it can inform of the particular concentration of the high frequency phase at one specific phase of the high frequency phase. More importantly, what matters here is not a metric of the coupling strength, but a reliable statistical test that selectively detects coupling, and my intuition is that counts in the phasephase plot may provide one. If a simple correction can be applied to give a sound statistical test, this could allow Belluscio and colleagues to look back at their data and see whether the significant counts indeed resist correction for multiple comparisons.
Indeed, we tended not to consider Belluscio et al. analysis a proper one; despite a bin count not being a metric of n:m coupling strength (c.f. our previous answers and acknowledged by the reviewer), we predicted difficulties in interpreting significant bin counts. For instance, a phasephase plot of 120 x 120 bins (i.e., 3^{o} bins) has 14400 bin counts. If a single bin count out of these 14400 bins is deemed statistically significant (after an appropriate correction), should we consider the time series n:m phasecoupled? We particularly did not consider that having only one bin deemed statistically significant would suffice. But then, what would be the minimal number of significant bins for assuming n:m phasecoupling? And would the appearance of striped patterns be necessary? Or would a few significant bins scattered at random locations suffice? Again, especially considering that there is a welldefined metric of n:m phasecoupling (R_{n:m}), we had particular difficulty in considering the bin count analysis a proper one. That said, though, we do respect the opposite opinion manifested by the reviewer, and have followed the suggestion of subjecting the pvalues of phasephase bin counts to a correction for multiple comparisons. These analyses show that the statistically significant bins in phasephase plots of white noise under the original framework of Belluscio et al. do not persist after correcting for multiple comparisons using HolmBonferroni (we considered the FDR correction too liberal since it led to significant bins in white noise). The results were true when building surrogate distribution using both Time Shift and Random Permutation procedures. This is now shown in new Figure 9A and stated/discussed in the fourth paragraph of the subsection “On diagonal stripes in phasephase plots” and in the third paragraph of the subsection “Statistical inference of phasephase coupling”. Please note that we created a new Results subsection to discuss the origin of stripes in phasephase plots (“On diagonal stripes in phasephase plots”), which led to editions and rearrangement of previous paragraphs and figures (in specific, former Figures 5B, 6 and 7 are now new Figures 6, 7 and 8, respectively).
In actual LFP data, we also found no significant bins after correcting for multiple comparisons when employing the same ± 200 ms timeshift surrogate procedure as in Belluscio et al. Interestingly, however, for very long epochs (510 minutes of concatenated data) we did find significant bins even after HolmBonferroni correction when employing randomly permutated – but not timeshifted – surrogates. We believe this is because small timeshifts as performed in Belluscio et al. ( ± 200 ms) would preserve phaselocking between the fundamental and harmonic frequencies (albeit in different phase relations) while the larger time shifts obtained through random permutations would not. These results are now shown in new Figure 10 stated/discussed in the last two paragraphs of the subsection “On diagonal stripes in phasephase plots” and in the last paragraph of the subsection “Statistical inference of phasephase coupling”. Please note that in the revised manuscript we leave room for this result being due to a weak but true coupling effect (subsection “On diagonal stripes in phasephase plots”, fifth paragraph, subsection “Statistical inference of phasephase coupling”, last paragraph, and subsection “Lack of evidence vs evidence of nonexistence”, last paragraph). Please see also our answers below.
B) About the very low pvalues in Belluscio's data, the authors suggest that "it is possible that, by analyzing a longer epoch length, the influence of theta harmonics becomes more apparent and would lead to lower pvalues, while the effect of the filteringinduced sinusoidality is washed out for both actual and surrogate epochs (of the same length)."
Well the manuscript previously demonstrates that n:m coupling measures as well as phasephase plots cannot tear apart asymmetrical waves from true crossfrequency coupling. Thus, if indeed pvalues remain lower than threshold when controlled appropriately, it could equally be due to asymmetrical theta or true thetagamma phase coupling, but the latter could not be dismissed.
As stated above and shown in new Figure 10, pvalues in phasephase plots of actual data are no longer considered statistically significant after appropriate correction for multiple comparisons when employing the surrogate procedure introduced in Belluscio et al. ( ± 200 ms time shifts), which in principle would settle this issue. However, we note that bin counts are deemed statistically significant even after correction when employing randomly permutated surrogates, especially for very long epochs (Figure 10A). We agree with the reviewer that n:m coupling measures cannot separate true coupling from asymmetrical waves; while we indeed tend to ascribe these results to harmonic influences, in the revised manuscript we now state our interpretation more cautiously and leave open the possibility of true thetagamma phase coupling (subsection “On diagonal stripes in phasephase plots”, fifth paragraph, and the Conclusion).
3) The reason for stripes in the phasephase plot of hippocampal data:
A) In the manuscript authors state that "the diagonal stripes in phasephase plots are due to theta harmonics and not to genuine gamma activity." If it were true harmonics, that would imply a consistent phase relationship between theta and gamma over sustained periods (not just over small periods as for white noise), i.e. would test positive for n:m coupling. So in my opinion here the effect is rather due to the temporary n:m alignment of phases, just as for white noise. In other words, stripes emerge also when analyzing white noise despite there being no harmonics in the signal.
We thank the reviewer for the insightful comment. Indeed, after performing new analyses on white noise and actual LFPs (which led to new Figure 9, Figure 10, Figure 10—figure supplement 1 and Figure 10—figure supplement 2), we came to agree with the reviewer that the temporary n:m phase alignment due to the sinusoidality induced by filtering is a main factor in determining the diagonal stripes in phasephase plots. These conclusions are backed by data shown in our new figures and explicitly stated in subsection “On diagonal stripes in phasephase plots”, last three paragraphs and subsection “Statistical inference of phasephase coupling”, last paragraph. However, our new analyses now reveal that in actual LFP signals there is a second influence, which is not accounted for solely by filtering. This effect is only statistically detected when analyzing very long epoch lengths (>100 s) and when contrasting original results (either phasephase plots or Rnm curves) against randomly permutated data (please see new Figure 10). Based on our analysis of sawtooth wave (in which random permutation was more sensitive to detect coupling than timeshifted surrogates, Figure 4—figure supplement 2), and given that the effect of filteringinduced sinusoidality is controlled for in surrogate data, we particularly believe that the weak (Rnm=0.03) but significant coupling detected in very long time series is due to theta harmonics (discussed in subsection “Statistical inference of phasephase coupling”, last paragraph). But in the revised manuscript we leave open the possibility of this low level of coupling as denoting true coupling (subsection “On diagonal stripes in phasephase plots”, fifth paragraph, subsection “Statistical inference of phasephase coupling”, last paragraph, subsection “Lack of evidence vs. evidence of nonexistence”, last paragraph, and the Conclusion).
B) In the rebuttal, the authors defend that "Believing in two different and genuine gamma activities – one coupled at 5 cycles per theta, the other at 4 cycles per theta – would go against the parsimonious principle." (this is when theta frequency evolves). However, as established by theoretical studies of coupled oscillators, n:m coupling is not a fixed property of the network but emerges as a combination of oscillator intrinsic dynamics and the characteristics of the connectivity pattern. See for example Arnold tongue: if we assume there is just one fixed gamma, a fluctuating theta and fixed connectivity, it is perfectly normal that n:m coupling will shift from 1:5 to 1:4 if the lower oscillator accelerates and the ratio of frequency goes from around 5 to around 4.
We thank the reviewer for calling our attention to these facts, point taken. We just note that the “parsimony argument” properly criticized by the reviewer was not present in the manuscript text but only stated in the rebuttal letter.
4) Asymmetrical theta vs. phasecoupled thetagamma in general
I am coming back to the authors' comment on this in the rebuttal, although this is no longer present in the manuscript. This is a comment for the authors benefit and speaks to the overall comment mentioned above. Sharp edges in LFP/EEG are not by itself an indication that we are measuring a single asymmetrical oscillation as the superposition of a slow and a fast n:m coupled oscillation can also give rise to sharp edges (Figure 4C; note that this can even be obtained with just two oscillations). Thus, sharp edges is no more a selective feature of asymmetrical oscillations than n:m coherence is selective of phasecoupled oscillators. In other words, visual inspection cannot tell us more than statistics.
We appreciate the thoughtful comments by the reviewer. (As stated by the reviewer, this issue is no longer present in our manuscript.)
5) It looks like "Random Perm" is more powerful test than "Time Shift", as mentioned in the first paragraph of the subsection “Lack of evidence for n:m phaselocking in actual LFPs”. If this is the case then it could be stated as a conclusion of the work that "Random perm" should be preferred over "Time Shift" (and of course of the "Scrambling" procedure), as it is less likely to miss existing effects (lower False Rejections rate).
We agree (and actually our new analyses of very long LFP epochs support such conclusion). This conclusion is now explicitly stated in the last paragraph of the subsection “Statistical inference of phasephase coupling”.
6) Conclusion: the very new paper by LozanoSoldevilla and colleagues (Frontiers Comp Neuro) provides another example of spurious crossfrequency coupling measures due to asymmetrical oscillations, could be worth referencing.
We now cite LozanoSoldevilla et al. et al. at the passage mentioned by the reviewer (currently in the last paragraph of the subsection “Statistical inference of phasephase coupling”).
https://doi.org/10.7554/eLife.20515.037Article and author information
Author details
Funding
Conselho Nacional de Desenvolvimento Científico e Tecnológico
 Robson SchefferTeixeira
 Adriano BL Tort
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior
 Robson SchefferTeixeira
 Adriano BL Tort
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors are indebted to the reviewers for many constructive comments and helpful suggestions. The authors are grateful to Jurij Brankačk and Andreas Draguhn for donation of NeuroNexus probes. Supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). The authors declare no competing financial interests.
Ethics
Animal experimentation: All procedures were approved by the institutional ethics committee of Federal University of Rio Grande do Norte (Comissão de Ética no Uso de Animais  CEUA/UFRN, protocol number 060/2011) and were in accordance with the National Institutes of Health guidelines.
Reviewing Editor
 Frances K Skinner, University Health Network, Canada
Publication history
 Received: August 10, 2016
 Accepted: December 5, 2016
 Accepted Manuscript published: December 7, 2016 (version 1)
 Version of Record published: December 29, 2016 (version 2)
Copyright
© 2016, SchefferTeixeira 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

 3,741
 Page views

 697
 Downloads

 39
 Citations
Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Neuroscience
 Physics of Living Systems

 Neuroscience