Unveiling Fast Field Oscillations through Comodulation

Visual Abstract

The term "comodulation" generally denotes the interaction between specific features of two oscillations, such as instantaneous phase, frequency, and amplitude. When oscillations have different frequencies, the interaction is called cross-frequency coupling (CFC). Among other functions, CFC was proposed to underlie communication between brain areas, sensory demodulation, multi-item or sequential working memory, and phase precession (Jensen and Colgin, 2007;Axmacher et al., 2010;Canolty and Knight, 2010;Jirsa and Müller, 2013;Hyafil et al., 2015). In CA1, several results converge in showing that theta phase modulates the amplitude of LG and HG (Colgin et al., 2009;Scheffer-Teixeira et al., 2012;Schomburg et al., 2014;Lasztóczi andKlausberger, 2014, 2016;Lopez-Pigozzi et al., 2016). In addition, theta has been recently shown to modulate the amplitude of an even faster rhythm, which has been called high-frequency oscillations (HFOs;120 -160 Hz;Scheffer-Teixeira et al., 2012;Tort et al., 2013). However, although gamma oscillations have been widely reported, theta-associated HFOs remain relatively unknown. Nevertheless, given the putative importance of neuronal oscillations to physiologic and pathologic brain functioning (Bragin et al., 1999;Wang, 2010;Branka k et al., 2012;Caixeta et al., 2013), it is likely that theta-associated HFOs have functional roles yet to be identified.
Possible reasons for the underreporting of thetaassociated HFOs include (1) lack of proper signal analysis techniques that allow their identification (such as CFC tools); (2) lower importance attributed to the highest part of the local field potential (LFP) spectrum, which is often ignored (many consider oscillations faster than 100 Hz to represent extracellular spikes and not true oscillations; Scheffer-Teixeira et al., 2013); (3) the high prevalence of recordings from the pyramidal cell layer (to capture spiking activity), whereas theta-HFO coupling has been observed in more superficial layers ; and (4) differences in recording techniques, such as electrode impedance, contact surface, and shape (single wires versus tetrodes vs probe contacts) potentially influencing HFO detection.
In the present work, we make a step toward better understanding theta-associated HFOs by outlining how these oscillations can be detected. To that end, we use spectral and CFC analyses to investigate different electrode types, recording locations, and behavioral states. We also contrast HFO activity within theta cycles with LG, HG, and multiunit activity. The results highlight unique properties of HFOs that separate them from gamma oscillations, as well as from spike contamination of LFP signals. This work should help in identifying thetaassociated HFOs in other laboratories worldwide.

Subjects
We used seven male Wistar rats from our breeding colony (2-3 mo old; 300 -400 g). Animals were maintained in a 12:12-h light-dark cycle, and recording sessions were performed during the light phase. All animal manipulations and care were approved by our local institutional ethics committee (CEUA/UFRN, protocol number 060/2011) in accordance with National Institutes of Health guidelines. These same animals were used in a recent publication of ours on cross-frequency phase-phase coupling (Scheffer-Teixeira and .

Surgery
Animals were anesthetized with a mixture of xylazine (10 mg/kg) and ketamine (100 mg/kg); supplementary doses of ketamine were injected if needed. After the skull was exposed and cleaned, target areas (fixation points and recording site) were marked using a stereotaxic apparatus and drilled. Six stainless steel screws were used for cap fixation. Two additional screws in the occipital bone (in contact with the dura mater) served as ground and reference for all electrodes. A thin layer of dental cement was applied around them. The silicon probe was lowered into the target area (CA1 region; AP, -3.6 mm; ML, 2.5 mm; DV, ϳ2 mm, but variable across animals). Ongoing neural activity was monitored during probe placement and impedance measurement. While passing through the CA1 pyramidal layer, successive electrode contacts transiently exhibited low-frequency spikes. After insertion, the probe was secured to the skull with dental cement. The assessment of electrode location was performed by benchmark electrophysiological signatures during recordings sessions (Fig. 2). We used six customdesigned silicon probes manufactured by NeuroNexus (four probes with 4320-m 2 contact areas, one with 703m 2 contact area, and one with 177-m 2 contact area). The probes had 16 channels in a linear vertical design with interelectrode distance of 100 m. Contacts spanned from the parietal cortex above the hippocampus to CA1 layers; because of differences in final probe depth among animals, the span of recorded locations across animals is larger than a single probe span (Fig. 3B). In one animal, we implanted an independent movable microdrive consisting of three tetrodes composed of 12.5-mdiameter platinum-iridium wires and four single 50-mdiameter tungsten wires, which were lowered until the stratum alveus/oriens. After surgery, animals received analgesics and local and systemic antibiotics and were monitored for 1 wk.

Recordings
After 1 wk of recovery, we recorded electrophysiological signals and videorecorded spontaneous behavior of the animals in an open field (1 ϫ 1 m) in ϳ4to 5-h sessions. We used RHA 2116 (Intan Technologies) 16channel recording system. The signals were amplified 200ϫ, bandpass filtered between 1 Hz and 7.5 kHz, and digitized at 25 kHz. For LFP analysis, we filtered the raw data between 1 and 500 Hz and down-sampled to 1000 Hz. For multiunit activity, we filtered the raw data between 800 Hz and 8 kHz and obtained the timestamps of spikes using thresholds established by visual inspection.

Histology
After the recording sessions, animals were killed for histologic validation of electrode positioning. To that end, we cut the brain in coronal sections and used Nissl staining to visualize probe/wire tracking in the parietal cortex and hippocampus.

Behavioral classification
Epochs of active waking (aWK) and REM sleep were identified by inspection of electrophysiological signals and videorecordings. aWK was defined as periods of theta activity and visible movements, and REM sleep was defined by the presence of theta, absence of movements, sleep postures, and preceding slow-wave sleep.

Data analysis
We used custom-written and built-in routines in Matlab (MathWorks). We also used routines from two third-party toolboxes: EEGLAB (Delorme and Makeig, 2004) and CircStat (Berens, 2009).

Power spectral density and 1/f fitting
We computed power spectral densities (PSDs) using the pwelch function (Signal Processing Toolbox; 2-s Hamming windows with 1-s overlap). Time-frequency power decompositions were obtained with the spectrogram function (Signal Processing Toolbox) and used to assist the classification of sleep-wake states. For comparing power of individual fast oscillations (LG, HG, or HFO) between aWk and REM sleep (Fig. 5C), we first fitted a 1/f curve using PSD values around the frequency band of interest and then obtained a normalized peak power value by subtracting the 1/f fit from the actual peak power value (Scheffzük et al., 2013).

Filtering, amplitude, and phase extraction
Filtering was achieved with the eegfilt function of the EEGLAB toolbox. The amplitude envelope and instantaneous phase of filtered LFP signals were obtained as the absolute value and angle of the analytic signal representation, respectively (hilbert function; Signal Processing Toolbox).

Phase-amplitude coupling and comodulation maps
We estimated phase-amplitude coupling (PAC) strength using the modulation index (MI) described in detail elsewhere (Tort et al., 2008(Tort et al., , 2010aScheffer-Teixeira et al., 2012;Caixeta et al., 2013). In brief, the MI measures how much a mean amplitude distribution over phase bins deviates from the uniform distribution. The comodulation map, or "comodulogram," is obtained by expressing the MI computed for multiple frequency pairs (bandpassed LFP signals) by means of a 2D pseudocolor map, in which the x-axis denotes the phase frequency, and the y-axis the amplitude frequency (that is, a warm color in a given comodulogram entry means that the phase of the X frequency modulates the amplitude of the Y frequency). The comodulograms were constructed using 0 -20 Hz (2-Hz steps, 4-Hz bandwidth) and 30 -300 Hz (5-Hz steps, 10-Hz bandwidth) as phase and amplitude frequencies, respectively.
In the 2D pseudocolor maps of Fig. 3C, the x-axis denotes the amplitude-providing frequency (25-200 Hz, 5-Hz steps, 10-Hz bandwidth), the y-axis denotes anatomic depth (zero is pyramidal cell layer of the hippocampus), and color represents the normalized PAC strength between the theta phase and the amplitude of the fast oscillation in x-axis; the normalization is performed for each fast oscillation, obtained by normMI ϭ [MImin(MI)]/[max(MI) -min(MI)], where max(MI) and min(MI) denote the maximum and minimum MI values across depths. In Fig. 8A, the x-axis denotes the phase-providing frequency (0 -16 Hz, 2-Hz steps, 4-Hz bandwidth), the y-axis denotes anatomic depth, and color represents PAC strength between the x-axis frequency and the amplitude of the fast oscillatory activity under study. The continuous time series of theta-HFO MI values shown in Fig. 5A was obtained using 30-s windows with 5-s overlap.
For the comodulograms and phase-amplitude coupling analyses, the amplitude and phase-time series were obtained from the same electrode. Thus, the local theta served as the phase reference. But we note that using either the local LFP or another (fixed) channel to extract the theta phase reference yield similar coupling strength for the fast local oscillation. This is because (1) theta is highly coherent along the parietal-CA1-dentate axis (Lubenov and Siapas, 2009), (2) the estimation of theta phase does not depend on theta amplitude (assuming a reasonable level of signal-to-noise ratio) and is therefore independent of the variability in theta amplitude along the dorsoventral axis, and (3) the modulation index does not take into account the preferred phase of maximal amplitude (i.e., the phase reversal of theta below the pyramidal cell layer does not influence coupling strength). The exception was Fig. 6, which displays the preferred theta phase of maximal amplitude for the different fast oscillations. In this case, for all fast oscillations, the reference theta phase was fixed and taken from stratum oriensalveus (note that the instantaneous theta phase is the same from the pyramidal cell layer up to the parietal cortex; Fig. 2 and Lubenov and Siapas, 2009).

Surrogates
In addition to comparing MI levels between aWK and REM sleep states, we also used surrogate analysis to infer the very existence of PAC (Tort et al., 2010a). To that end, we obtained surrogate MI values using a "split-invertsplice" method: while the phase time series was left intact, we randomly split the amplitude time series and concatenated the two parts inverted. For each comodulogram in Fig. 3A, we computed a distribution of 200 surrogate comodulograms; actual MI values below at least one surrogate MI value were set to zero.

Phase-energy plots
In Fig. 6, we estimated the preferred theta phase (i.e., the phase of maximal amplitude) of each fast oscillation of interest by visual inspection of phase-energy plots (see Fig. 6A for examples). To obtain the phase-energy plot, we used Morlet wavelet spectral decomposition using scaled pseudofrequencies of 20 -200 Hz in 2-Hz steps. We then calculated the mean wavelet amplitude over theta phase bins of 2°. To compensate for the 1/f decay, the mean amplitude of each frequency was normalized by the average across all phase bins. In this analysis, a fixed theta reference at stratum oriens-alveus was used (which provides the same phase reference as in stratum pyramidale but with a better signal-to-noise ratio).

Independent component analysis
To obtain a better estimate for the spatial origin of the fast oscillations recorded with the linear probes, we used independent component analysis (ICA) as described in Schomburg et al. (2014). To that end, we first filtered the LFPs Ͼ30 Hz and then used the KDICA algorithm (Chen, 2006); aWK and REM states were analyzed separately. The resulting independent components (ICs) were used as the amplitude-providing signal for the comodulogram analysis, and a stratum oriens-alveus LFP provided the phase signal. The comodulograms were next inspected and used to classify the ICs as containing LG, HG, or HFO. The absolute value of the voltage loadings of ICs exhibiting isolated gamma or HFO activity were interpreted as an estimate of the spatial source (see Fig. 4 for examples).

Fast oscillatory bursts per theta cycle
For each fast frequency (LG, HG, and HFO), we defined the occurrence of an oscillatory burst as the crossing of the amplitude envelope of a threshold of 4 SD above the mean amplitude. We identified individual theta cycles by the timestamps of consecutive valleys (i.e., a cycle was defined to start/finish at the valley) in bandpass-filtered LFPs from 4 to 20 Hz, a frequency range previously shown to provide a good fit to asymmetric theta cycles (Belluscio et al., 2012). We then quantified the fast oscillatory bursts per theta cycle. In this analysis, we used a channel in stratum oriens-alveus as theta reference for all fast oscillations.

PAC strength controlled for theta amplitude
In Fig. 5E, we compare PAC strength between theta phase and the amplitude of fast oscillations (LG, HG, and HFO) as a function of theta amplitude. This was achieved by first obtaining the phase and amplitude time series for the entire signal. We then z-scored the mean theta amplitude of each theta cycle. Finally, we classified theta cycles by their mean amplitude (in units of z score) and computed PAC strength using only time series periods associated to theta cycles of a fixed theta amplitude bin (bin width ϭ 0.2 z score). We excluded bins containing Ͻ10 s of signal to ensure reliable MI estimation.

Spike-phase coupling
We estimated spike-phase coupling (SPC) strength from the probability distribution of spikes over phase bins (extracted from bandpass-filtered LFPs). We fitted a von Mises distribution and used the concentration parameter (kappa) as a measure of SPC (Berens, 2009). The higher the kappa, the higher the SPC.

Spike-and oscillation-triggered LFP averages and time-energy plots
To construct triggered LFP averages, we used 500-ms windows around the trigger timestamp (ripples and theta peaks in Fig. 2, spike times in Fig. 9). We averaged either raw LFP signals (Fig. 2B left panel and Fig. 9) or filtered signals (Fig. 2B, right, and 2C). The time-energy plots in Fig. 9 were obtained from the spike-triggered averages of the amplitude envelope of filtered LFP signals (20 -250 Hz, 1-Hz steps, 10-Hz bandwidth). In these plots, the amplitude of each frequency was normalized by its maximum value to control for 1/f.

Statistics
We opted to exclude p-value-only based decisions because of recent concerns in the field and as an effort to increase statistical quality (Wasserstein and Lazar, 2016). Therefore, in addition to p-values (calculated using twotailed tests), we based our statistical inferences on 95% confidence intervals (CIs) of the mean, effect size (Cohen's d; calculated as suggested by Lakens [2013]), and more transparent data presentation (Weissgerber et al., 2015). We avoided nested data by analyzing one measure per animal, and unless noted otherwise (e.g., in Fig. 1), the sample size is the number of animals and not electrodes, which could inflate power. Each sample consisted of either the mean or the median of scores (e.g., MI, power, or phase) over all eligible electrodes per animal. Statistical analyses are summarized in Table 1.

Theta-associated HFOs can be detected with different electrode types
We first sought to investigate whether the underreporting of HFOs could be due to differences in electrode types among laboratories. We have previously obtained prominent records of theta-associated HFOs using lowimpedance, single tungsten wires of 50-m diameter (Scheffer-Teixeira et al., 2012). We were interested in knowing whether HFOs can also be detected with electrodes of different geometry, such as planar probe contacts and the thinner wires used in tetrodes. We thus compared LFP recordings obtained from single 50-mdiameter tungsten wires (1963 m 2 contact area), from multisite probes of large (4320 m 2 ) and intermediate (703 m 2 ) contact areas and from tetrodes made of 12.5-mdiameter platinum-iridium wires (123 m 2 contact area).
We could detect prominent theta-HFO coupling in the parietal cortex during REM sleep in all electrode types (Fig. 1A). The highest modulation index was observed for an LFP recorded with the single 50-m-diameter tungsten wire ( Fig. 1A, middle left). We next investigated whether electrode impedance affects HFO detection. We found that electrode impedance (estimated in vivo) has only mild impact on theta-HFO coupling strength (Fig. 1B, C). The scatter plot between CFC strength and impedance shown in Fig. 1B reveals a weak negative correlation (Pearson's r: -0.255, p ϭ 0.022; all electrodes in a same animal along five sessions). The red and blue dots in Fig. 1B highlight changes in impedance and CFC strength in the same electrode between the first and last sessions. Despite the major increase in impedance (2.4-fold, from 802 to 1945 k⍀), the last recording session exhibited theta-HFO coupling as prominent as in the first session. Moreover, although impedance was significantly lower in single 50m-diameter wires than in 12.5-m-diameter tetrode wires (Fig. 1C, left; p ϭ 7.077 ϫ 10 Ϫ14 , t test), these electrodes exhibited comparable levels of theta-HFO coupling (Fig. 1C, right). Only 4 of 20 LFP samples (20%) recorded with the thicker wires had higher theta-HFO coupling than all samples (n ϭ 60) collected with tetrode wires. The median level of theta-HFO coupling was mildly higher for 50-than 12.5-m-diameter wires. Although this difference reached statistical significance (Fig. 1C, right; p ϭ 0.007, t test), the effect size measured using Cohen's d was 0.6 (associated with 76.4% overlap), much smaller than the observed effect size for the difference in impedance, 2.7 (associated with 17.7% overlap). Finally, we computed PSDs and found that both 50-and 12.5-mdiameter wires had comparable levels of theta power; moreover, both electrode types exhibited a clear power bump in the HFO range (120 -160 Hz; Fig. 1D). In all, we conclude that differences in electrode geometry and impedance have only minor influences on HFO detection.

Recording location and behavioral state are the main factors determining HFO detection
In our previous work, we showed that PAC patterns depend on recording depth across the CA1-dentate axis, with theta-HFO coupling being most apparent in more superficial layers than theta-HG coupling (Scheffer-Teixeira et al., 2012). Here, we extend these results by mapping comodulation patterns from multiple, linearly spaced contacts from the parietal cortex to the hippocampus, separating the analyses in two behavioral states: aWK and REM sleep (see Materials and Methods). Fig. 2 shows a representative example of the estimation of contact positions in a multisite linear probe using benchmark electrophysiological signatures, such as sharp-wave reversal, ripple power, theta reversal, and theta power Ylinen et al., 1995). We classified recording depth using as reference the distance from the CA1 pyramidal cell layer (depth ϭ 0; positive and negative depths denote regions above [dorsal] and below [ventral] s. pyramidale, respectively; Fig. 3B). Fig. 3A shows phase-amplitude comodulograms computed at the different recording depths; the left column shows all comodulograms at the same MI scale, and the right column displays each comodulogram on its own MI scale to allow qualitative comparisons of coupling patterns. The comodulograms revealed the same qualitative pattern in each animal and corroborated our previous observations (Scheffer-Teixeira et al., 2012): theta-HFO coupling was detected above the pyramidal cell layer, whereas theta-HG coupling presented a richer picture: it was high in the stratum oriens and pyramidal cell layer, decreased in the stratum radiatum, and achieved maximum strength in the stratum lacunosum-moleculare (Fig. 3A). In both aWK and REM sleep, theta-HFO coupling was most apparent at recordings from the parietal cortex. Interestingly, a careful inspection of Fig. 3A reveals that the transition from theta-HFO to theta-HG coupling occurs at the corpus callosum during aWK (comodulograms at 300 and 400 m); however, during REM sleep, the transition between coupling patterns occurs much deeper, around the CA1 pyramidal cell layer (comodulograms at 100, 0, and -100 m). These results therefore show that recording location across the parietal-hippocampal axis is a crucial factor for the detection of theta-HFO coupling.
During aWK, we observed weak levels of theta-LG coupling in the parietal cortex (co-occurring with theta-HFO coupling) and in stratum radiatum (Fig. 3A). During REM sleep, theta-LG coupling was not as apparent in the . Electrode type and material, from left to right: linear planar silicon probe, tungsten single wire, linear planar silicon probe, platinum-iridium tetrode wire. Each comodulogram was obtained from one recording session in one animal (three different animals were used; the single and tetrode wire recordings were from the same animal and session). B, Left, scatter plot of electrode impedance and theta-HFO modulation index. Data obtained from 16 electrodes positioned in the parietal cortex in the same animal across five sessions. The thick line depicts the linear fit. Right, example comodulograms obtained from the same electrode in sessions 1 and 5 (red and blue circles in the scatter plot; electrode area: 123 m 2 ); red dashed box indicates region used to compute the strength of theta-HFO comodulation. Note a substantial increase in impedance along with a mild decrease in comodulation. C, Electrode impedance (left) and theta-HFO modulation index (right) measured for electrodes of 123 m 2 (n ϭ 12) and 1963 m 2 (n ϭ 4) contact area across five sessions for the same animal as in B. Boxes indicate 25th and 75th percentiles, red trace denotes the median, and dots show individual values. Note that although larger electrodes have much lower impedance, comodulation strength does not differ as much. D, Mean power spectra for the same recordings as in B and C. Thick and dashed lines represent median and 25th-75th percentiles, respectively. Arrows highlight a power peak in the HFO range in both electrode types.
comodulograms, except for a small comodulation "island" in the LG range in stratum radiatum (in -200 m from the pyramidal cell layer; Fig. 3A). Finally, we also observed an unexpected coupling between delta phase and the amplitude of oscillations ϳ100 Hz in some electrode depths during aWK (Fig. 3A 3A, these plots reveal highest theta-HFO coupling strength at the parietal cortex and highest theta-HG coupling near the hippocampal fissure. Theta-LG coupling was weak but apparent at the parietal cortex and stratum radiatum during aWK. Curiously, theta-LG coupling seemed to appear near the pyramidal cell layer during REM sleep. Although at first glance such a finding is very interesting, we suspect it reflects a "contamination" of the LG range by a "tail" (or leakage) of theta-HG coupling, which considerably increases during REM sleep. We also note that a similar phenomenon occurs for the estimation of theta-HFO coupling in the dentate gyrus, which gets heavily contaminated by theta modulation of spike-leaked HFOs (SLHFOs; discussed below). Finally, using the same linear probe data, we also computed phase-amplitude comodulograms for ICs of fast oscillatory activity as introduced in Schomburg et al. (2014). The voltage loadings tended to be highest at the parietal cortex for ICs selectively exhibiting theta-HFO coupling, at stratum radiatum for ICs with selective theta-LG coupling, and at stratum lacunosum-moleculare for ICs exhibiting theta-HG coupling (see Fig. 4 for examples). The left column of Fig. 3A shows that the maximum coupling strength substantially increases from aWK to REM sleep for theta-HG and theta-HFO coupling, but not for theta-LG coupling (note also y-axis scales in Fig. 3C). We proceeded to further investigate quantitative differences between these two theta states. Fig. 5A shows a representative example of variations of theta power (expressed as theta/delta ratio) and theta-HFO coupling strength in time. Note the high similarity of the two timeseries across sleep-awake states and the highest values during REM sleep (see also Fig. 5B). We next compared the power of the different fast oscillations between the two behavioral states. The PSD of aWK periods displayed higher power levels for frequencies Ͼ50 Hz than the PSD of REM epochs; however, the aWk PSD did not exhibit a clear power bump as reported in Fig. 1D (not shown, but see Scheffzük et al., 2011Scheffzük et al., , 2013. To avoid spurious estimates resulting from generalized augmented power (probably due to muscular activity during aWk), we normalized power by subtracting, from the original PSD, the 1/f fit computed using power values around the frequency band of interest (Scheffzük et al., 2013; see Materials and Methods; the boxplot distributions show all electrodes across animals, but for computing CIs and p-values, each sample consisted of the mean over electrodes per animal). Under this framework, neither LG power nor HG power differed between aWK and REM sleep (Fig. 5C Theta power during REM sleep was also significantly higher than during aWK ( Fig. 5D; REM-aWK CI 95 ϭ [0.11; 0.2], p ϭ1.52 ϫ 10 Ϫ4 ). Differences in theta power between aWK and REM could thus explain the higher phaseamplitude coupling observed during REM. In Fig. 5E, however, we show that when controlling for theta amplitude, coupling strength was visibly higher during REM sleep (blue) compared to aWK epochs (green) for both HG and HFO, but not for LG, irrespective of power level. Therefore, our results show that comodulation strength highly depends on behavioral state and that the increase in theta-HG and theta-HFO coupling during REM sleep is not explained by changes in theta power.

Different fast oscillations occupy different theta phases and theta cycles
The results thus far show that the modulation of LG, HG, and HFO amplitude by theta phase is spatially seg-   regated. We next investigated whether these fast oscillations temporally segregate in respect to theta cycle and theta phase (Figs. 6 and 7). During both aWK and REM sleep, we found amplitude peaks for different fast oscillations in different theta phases (Fig. 6A, B). Using theta phase above the pyramidal cell layer as reference (theta peak corresponds to 180°), we found that HG amplitude (recorded at s. lacunosum-moleculare) was maximal near the theta peak, and HFO amplitude (recorded in parietal cortex) was maximal right after the theta peak and HG maximal amplitude. Interestingly, we found that LG activity in parietal cortex (LGp) was maximal at the ascending theta phase, whereas LG in the hippocampus (LGh, recorded at s. radiatum) was maximal at the descending theta phase. We also found evidence for segregation of the theta cycles in which different oscillatory bursts occur (Fig. 7). For each fast oscillation, we obtained the timestamps in which the amplitude envelope was 4 SD above its mean and classified theta cycles according to oscillatory bursts. Possible cases are theta cycles containing no peak of fast oscillation amplitude, an amplitude peak for just one fast oscillation, or amplitude peaks of two or three fast oscillations. We found that about half of all theta cycles did not contain any amplitude peak (52 Ϯ 3% [mean Ϯ SD across animals]). Among all theta cycles, we found that 18 Ϯ 1% contained LG (10% exclusively), 22 Ϯ 2% HG (12% exclusively), and 23 Ϯ 2% HFO (13% exclusively); that is, 73% of theta cycles with a fast oscillatory burst were specific for LG, HG, or HFO (Fig. 7B). A minority of theta cycles contained two or more fast oscillatory bursts of different frequencies. In all, these results show that theta phase modulation of distinct fast oscillations tends to segregate into different theta cycles.

Theta-associated HFOs differ from spike contamination of LFP signals
Previous reports have raised the possibility that fast LFP frequencies (Ͼ100 Hz) stem from spectral leakage of spike waveforms recorded extracellularly (Ray et al., 2008a(Ray et al., , 2008b(Ray et al., , 2008cRay and Maunsell, 2011;Belluscio et al., 2012). Namely, abrupt, nonsinusoidal deflections of the extracellular voltage lead to spurious fast oscillations in filtered signals (Kramer et al., 2008). We previously corroborated that such effects of extracellular spikes do occur and called the spurious LFP activity originated from it "spike-leaked HFOs" (SLHFOs; Scheffer-Teixeira et al., 2013;Tort et al., 2013). But we have also shown that not all LFP activity Ͼ100 Hz is due to spike contamination; in fact, we showed that theta-associated HFOs are genuine oscillations . Here we aimed to reproduce and extend these observations. continued views of the LG range. B, Estimated electrode locations (adapted from Paxinos and Watson, 2004). Depth units are relative to the hippocampal pyramidal layer. C, Left, 2D maps of normalized theta coupling strength as a function of amplitude frequency and anatomic depth (average over four rats). For each amplitude frequency, 0 denotes the minimum modulation index and 1 the maximum. Right, mean theta coupling strength (n ϭ 4 rats) plotted as a function of recording depth separately for each fast oscillation (LG, HG, and HFO) and brain state (aWK and REM). LG, 30 -60 Hz; HG, 60 -100 Hz; HFO, 120 -160 Hz. Spike contamination of the LFP signal in some cases can leak down to 100 Hz (Ray et al., 2008a). Accordingly, in our previous study, we showed that the frequency range of theta-coupled SLHFOs, as assessed by phaseamplitude comodulograms, may also comprise the HFO band . To separate SLHFO activity from genuine HFOs, in this work we filtered LFPs Ͼ150 Hz. Notice therefore that although we cannot avoid contamination of the estimation of theta-HFO coupling by theta-SLHFO coupling (as happens in the dentate gyrus; Fig. 3C), by filtering Ͼ150 Hz we can minimize contamination in the other direction, that is, contamination of the estimation of theta-SLHFO coupling by theta-HFO coupling.
Theta-associated HFOs are best detected above the pyramidal cell layer (Fig. 8A, top left) and appear in the phase-amplitude comodulogram as a circumscribed "island" between ϳ110 and 160 Hz in the amplitude frequency axis (Fig. 8A, top right). On the other hand, SLHFO modulation by theta is most prominent at the pyramidal cell layer (Fig. 8A, bottom left) and in the dentate gyrus (see depths below -500 m in Fig. 3A; see also Fig. 5 in Scheffer-Teixeira et al., 2012), and the phase-amplitude comodulogram displays a wide range of modulated amplitude frequencies (Fig. 8A, bottom right). Importantly, whereas HFO amplitude peaks after the theta peak (s. pyramidale reference), both SLHFO amplitude and CA1 multiunit activity are maximal near the theta trough (  Scheffzük et al., 2011). Each data point corresponds to an electrode (LG, n ϭ 30 electrodes across 7 rats; HG, n ϭ 20; HFO, n ϭ 36). CI 95% for the paired difference between states was computed using the mean value over eligible electrodes for each rat (n ϭ 7 animals; for Cohen's d and p-values, see Table 1). Bottom, modulation index computed for theta phase and the amplitude of each frequency range. Notice larger theta-HG and theta-HFO coupling during REM sleep. D, Power spectral density (left) and mean theta power (right) during aWK and REM. Continuous thick lines represent median, and dashed thin lines 25th-75th percentiles. E, Phase-amplitude coupling strength as a function of theta amplitude for different fast oscillations during aWK and REM sleep. In these analyses, eligible electrodes were selected based on inspection of comodulograms; for LG, we used recordings from s. radiatum, s. pyramidale, and parietal cortex; for HG, from s. lacunosum-moleculare and s. pyramidale; for HFO, from stratum oriens-alveus and parietal cortex. 8B). Interestingly, during REM sleep, the distribution of CA1 spikes over the theta cycle exhibits a second maximum near the theta peak (Poe et al., 2000;Senior et al., 2008), and SLHFO amplitude displays a similar bimodal profile (Fig. 8B, right). Of note, the phase of the second spiking activity peak matches the phase of maximal HFO amplitude within a theta cycle. It is therefore possible that the increase in HFO activity observed during REM sleep (Fig. 5) drives spiking activity in CA1 and leads this second peak in spike probability.
Multiunit activity couples to theta phase in both aWK and REM sleep (Fig. 8A). We next investigated the dependence of spike coupling on LFP frequency and depth across the parietal-hippocampal axis. In Fig. 9A (bottom), we show examples of spike-phase distributions of CA1 multiunit activity for multiple LFP oscillations recorded from the pyramidal cell layer during REM sleep. We observed prominent spike-phase coupling only for theta and SLHFO, but not for LG, HG, or HFO (these results refer to the multiunit spike pool and do not imply lack of coupling for single, well-isolated units). It should be noted that the spike coupling to SLHFO phase is spurious; this effect simply reflects that filtered spike waveforms modulate their own spike timing . When referencing to the phase of LFP oscillations in other depths, we found that CA1 multiunit activity coupled to theta phase in all recorded layers (Fig. 9A, top). We did not find prominent multiunit coupling for LG, HG, or HFO phase in any layer, not even in layers where these oscillations are most prominent. On the other hand, SLHFO phase (spuriously) modulated multiunit activity at the pyramidal cell layer and dentate gyrus, although this effect LGp HFO Parietal max min Norm. Amp. Figure 6. Preferred theta phase of maximal amplitude for different fast oscillations. A, Top, representative phase-energy plots for LFPs recorded in the parietal cortex (left) and stratum radiatum (right) during REM sleep. Bottom, the same in rescaled and zoomed-in view to allow inspection of the LG range. Note different preferred theta phases for LG in parietal cortex (LGp) and hippocampus (LGh). B, Preferred theta phase of maximal amplitude. Linked circles show data from the same animal. Bars depict the circular median over animals. For each fast oscillation and animal, we considered the electrode that displayed the maximum modulated activity at the corresponding frequency band on visual inspection of phase-energy plots (e.g., HFO was taken from the parietal cortex and HG from stratum lacunosum-moleculare). For all fast oscillations, we used the theta phase from an LFP at stratum oriens-alveus.
was also somewhat apparent in other recording locations (Fig. 9A). Finally, we constructed phase-energy plots triggered by spike times during aWK and REM sleep epochs ( Fig.  9B; see Materials and Methods). The bottom traces in Fig.  9B show the spike-triggered averages of raw LFPs, which-consistent with spike-phase distributions-reveal spike timing preferentially occurring near the theta trough. The top plots show spike-triggered averages of the instantaneous amplitude of filtered LFP signals (each frequency was normalized by its maximum value to control for 1/f). Note the prominent SLHFO activity in the aWK plot, that is, a major increase in amplitude for a wide, unbound range of fast frequencies coinciding with the time of spiking (t ϭ 0) near the theta trough (Fig. 9B, left) Figure 7. Different frequency bands use their own theta cycle. A, Examples of raw traces along with filtered LFPs, as labeled. B, Venn diagram of theta cycles classified according to the detection of oscillatory bursts. A burst event was defined to occur when the instantaneous amplitude envelope of the band-filtered signal crossed Ͼ4 SD from the mean. Percentages of theta cycles with oscillatory bursts include intersections. When discounting intersections, burst percentages sum to 48%. LG, yellow; HG, red; HFO, purple. Results are shown only for REM sleep epochs; results for aWK were similar. For each animal and fast oscillation, we analyzed one electrode, selected based on visual inspection of comodulograms. The electrode was located in the parietal cortex for HFO, at s. radiatum for LG, and at s. lacunosum-moleculare for HG. Percentages represent the mean over animals.  . Theta-HFO coupling is strongest above the pyramidal cell layer, which is dominated by theta-SLHFO coupling. Right, mean comodulograms (n ϭ 4 rats) for electrodes located in the parietal cortex (top) and the hippocampal pyramidal cell layer (bottom). Note genuine HFO and spurious SLHFO coupling to theta phase. B, Top, spike probability (black bars; left y-scale) plotted along with the amplitude (right y-scale) of HFO (120 -160 Hz; purple lines, parietal cortex recordings) and SLHFO (Ͼ150 Hz; green lines, pyramidal cell layer recordings) per theta phase during aWK and REM (group data; n ϭ 4 rats). Dashed vertical orange lines indicate 180°and 540°(pyramidal theta peaks). Bottom, same as above, but with values normalized between 0 and 1. Note similar phase distributions between SLHFO amplitude and spike probability, which have a maximum near the trough of the theta wave. phases (Fig. 9B, right), reflecting the bimodal distribution of spiking phases shown above (Fig. 8B). Interestingly, during REM sleep, theta-associated HFOs are also noticeable near the theta peak in this analysis (Fig. 9B, right).

Discussion
In the present work, we characterized the patterns of fast network activity that coexist with theta oscillations in field potential recordings from the hippocampus and parietal cortex. Our results show that phase-amplitude coupling occurs between theta phase and the amplitude of multiple higher-frequency oscillations: LG, HG, and HFO.
These oscillations differ from each other by frequency band, preferred theta phase, theta cycle of occurrence, and anatomic location. We also characterized spurious activity resultant from spike contamination of LFPs (SL-HFO), common in cell-dense regions such as the pyramidal cell layer and dentate gyrus.

Relation to previous findings
Different groups use different nomenclatures for fast field oscillations in the hippocampus. For instance, what we here refer to as LG and HG have also been termed "slow gamma" and "middle gamma," respectively (Colgin  Belluscio et al., 2012;Schomburg et al., 2014). The terms "fast gamma" and "epsilon band" have been used for hippocampal activities Ͼ100 Hz concomitant with theta waves (Belluscio et al., 2012;Schomburg et al., 2014). However, although some groups (Scheffzük et al., 2011;Branka k et al., 2012) have called fast gamma the same HFO activity described in the present work, the "fast gamma/epsilon band" activity reported by others may actually correspond to what we refer to as SLHFO. Namely, in Belluscio et al. (2012) and Schomburg et al. (2014), the authors report fast gamma activity occurring at the trough of the theta cycle, just where the units are most active (see also Lasztóczi and Klausberger, 2017). Moreover, this fast gamma activity exhibited a second local maximum near the theta peak during REM sleep (Belluscio et al., 2012), mirroring the changes in phase distribution displayed by the units from aWK to REM and similarly to what we have reported to occur for SLHFOs in the present work. On the other hand, the fast gamma in Scheffzük et al. (2011Scheffzük et al. ( , 2013 and Branka k et al. (2012) occurred after the peak of the theta wave, and its coupling to theta was most prominent during REM sleep. We prefer not calling spike contaminations of the LFP signal a type of gamma and rather reserve this term for genuine oscillatory activity. In this sense, we are not against calling theta-associated HFOs "fast gamma" as long as this term is uniquely used to refer to this rhythm. Note that this issue illustrates that defining brain oscillations solely based on frequency range is not ideal and may lead to confusion, as previously voiced elsewhere (Tort et al., 2010b;Cole and Voytek, 2017).
Hippocampal gamma oscillations have been the subject of intense research over the last decades Burchell et al., 1998;Csicsvari et al., 2003;Montgomery and Buzsáki, 2007;Colgin et al., 2009;Scheffer-Teixeira et al., 2012;Schomburg et al., 2014;Lasztóczi andKlausberger, 2014, 2016;Scheffer-Teixeira and Tort, 2016). In contrast, theta-associated HFOs have only recently been described (Scheffer-Teixeira et al., 2012;Tort et al., 2013). To try to understand why HFOs remain underreported, we first investigated whether differences in recording techniques could affect HFO detection. Although we found that electrode impedance does influence HFO detection, this influence is rather mild (Fig.  1). In fact, our results demonstrate that HFOs can be detected with different electrode types, such as thin wires used in tetrode recordings, single isolated wires of thicker diameter, and one-sided rectangular probe contacts (Fig.  1). Moreover, HFOs have also been previously detected in surface EEG recordings through steel screws placed at the dura mater over the parietal cortex (Scheffzük et al., 2011;Zhang et al., 2016). On the other hand, our results reveal that recording location and behavioral state are the two main determinants of HFO detection. During aWK, HFOs may be detected in the parietal cortex, where their amplitude is weakly coupled to theta phase (Fig. 3). During REM sleep, however, both HFO power but mainly theta-phase modulation of HFO amplitude become considerably higher, and theta-HFO coupling can also be detected in the hippocampus (Fig. 3A).
In our previous work (Scheffer-Teixeira et al., 2012), we used handmade electrode bundles to simultaneously record from different hippocampal layers. We showed that (1) the amplitudes of HFO and HG couple to theta phase, (2) coupling patterns depend on the hippocampal layer, (3) coupling strength is positively correlated with theta amplitude, and (4) HFO and HG have maximum amplitude in different theta phases. Here we were able to reproduce all these results in an entirely new dataset (different subjects, electrodes, recording system, and laboratory). We used commercial linear probes with a higher number of contacts, which allowed for a more fine-grained analysis of comodulation patterns and a greater span of recording locations than in our first study (Fig. 3). The results show that the recording depth in which comodulation patterns switch between theta-HFO and theta-HG coupling depends on behavioral state, a feature we had not previously recognized: during REM sleep, theta-HFO coupling is apparent down to the pyramidal cell layer, whereas during aWK, theta-HFO coupling is mostly restricted to the parietal cortex. In addition, here we were able to observe coupling between theta and LG in some recorded regions (parietal cortex and stratum radiatum), which was not detected in our previous study. Theta-LG coupling strength, however, was much weaker than theta-HG and theta-HFO coupling. Moreover, as opposed to theta-HG and theta-HFO coupling, theta-LG coupling did not increase during REM sleep.

aWK versus REM sleep
Neuromodulatory content differs between aWK and REM sleep. During REM sleep, the hippocampus exhibits decreased serotoninergic and noradrenergic tonus and a higher cholinergic drive (Kametani and Kawamura, 1990;McCarley et al., 1995;Park et al., 1999;Graves et al., 2001;Lee et al., 2005). Theta activity during REM sleep may relate to higher acetylcholine release in the hippocampus (Kametani and Kawamura, 1990;Monmaur et al., 1997;Lee et al., 2005), which would thus characterize it as type 2 theta (Kramis et al., 1975). Here we found augmented theta amplitude during REM sleep in comparison to aWK (Fig. 5D), a profile similar to that evoked by cholinergic agonists (Markowska et al., 1995;Monmaur et al., 1997). Of note, the higher theta-HG and theta-HFO coupling during REM sleep are not explained by changes in theta power (Fig. 5C-E): during REM sleep theta power increased Ͻ1.5-fold, whereas theta-HG and theta-HFO coupling increased Ͼ10-fold. Interestingly, Newman et al. (2013) reported that scopolamine reduces theta-HG coupling in the entorhinal cortex, whereas theta-LG is unaffected. The dependence of theta-HG coupling on cholinergic transmission may be related to its increase during REM sleep; by the same token, the lack of increase in theta-LG coupling during REM sleep may be related to its insensitivity to scopolamine as observed in Newman et al. (2013). On the other hand, the dependence of theta-HFO coupling on cholinergic modulation is currently unknown. In any event, both acetylcholine and REM sleep have been associated with memory consolidation (Hasselmo, 1999(Hasselmo, , 2006Stickgold and Walker, 2007;Nishida et al., 2009), and it is thus possible that enhanced theta-HFO and theta-HG coupling during REM sleep play functional roles.

Coordination of fast oscillations in different theta cycles and theta phases
Bursts of fast oscillatory activity tend to happen in different theta cycles: most theta cycles containing at least one oscillatory burst were not superimposed (Fig.  7B). Consistent with our findings, Colgin et al. (2009) found that ϳ54% of theta cycles contained at least one gamma burst, and most cycles contained just one type of gamma. In contrast, however, Belluscio et al. (2012) found that 41% of theta cycles contained both low and high gamma. Although the reason for these discrepancies is not clear, they could be due to methodological differences (for instance, our burst threshold was computed using only theta periods; they used the whole recording session) or particularities in behavior or cognitive demands between studies.
The amplitudes of LG, HG, and HFO peaked at different theta phases. HG was maximal near the theta peak and HFO soon afterward (at ϳ180°and ϳ210°respectively; here and below, the phases correspond to theta recorded above s. pyramidale, that is, before the phase reversal). The theta phase of maximal HG amplitude observed here is similar to that previously reported (Tort et al., 2008;Chen et al., 2011;Belluscio et al., 2012;Scheffer-Teixeira et al., 2012;Schomburg et al., 2014;Yamamoto et al., 2014;but see Colgin et al., 2009). As discussed above, the theta phase of maximal HFO amplitude is similar to that found in previous reports describing a same pattern of bounded LFP activity between 110 and 160 Hz (Tort et al., 2008;Sirota et al., 2008;Scheffzük et al., 2011;Branka k et al., 2012;Scheffer-Teixeira et al., 2012;Tort et al., 2013;Zhang et al., 2016) but differs from other studies analyzing recordings from the pyramidal cell layer (Belluscio et al., 2012;Schomburg et al., 2014), probably because of the influence of spike contamination of the LFP signal in the latter.
Interestingly, in the parietal cortex, LG amplitude was maximal at the ascending phase of the theta cycle (ϳ90°), whereas in the hippocampus, LG amplitude peaked at the descending theta phase. Newman et al. (2013) found that LG recorded in the medial entorhinal cortex was maximal at the ascending phase, i.e., similar to our parietal cortex recordings. On the other hand, previous studies have reported different preferred theta phases for LG recorded in the hippocampus (e.g., Colgin et al., 2009;Yamamoto et al., 2014;Lasztóczi and Klausberger, 2017), although most of them seem to find maximal LG at the descending phase, as observed here (Chen et al., 2011;Belluscio et al., 2012;Schomburg et al., 2014;Lasztóczi andKlausberger, 2014, 2016;Lopez-Pigozzi et al., 2016). For example, Chen et al. (2011) found LG at the descending phase but with the exact phase depending on speed, Belluscio et al. (2012) found LG at the descending phase closer to the theta peak, and Schomburg et al. (2014) found LG at the descending phase closer to the theta trough. The reason for such differences is presently un-clear; we suspect it may have to do with traveling wave properties (Lubenov and Siapas, 2009;Patel et al., 2012) of LFP oscillations along with differences in recording location across the antero-posterior or medio-lateral axes of the hippocampus (see Schomburg et al., 2014).
Interestingly, although hippocampal neurons spike at the trough of the theta during aWK, during REM sleep they also tend to spike more at the descending phase of the cycle after the theta peak ( Fig. 8B; Poe et al., 2000;Mizuseki et al., 2011;Belluscio et al., 2012), where we found HFO amplitude to be maximal. Thus, HFOs may potentially be related to spiking activity at this theta phase during REM sleep. Previous work at the single-unit level has shown that some pyramidal cells shift their preferred theta phase of spiking during REM sleep (Poe et al., 2000). Furthermore, Mizuseki et al. (2011) showed that whereas nonshifting pyramidal cells are closer to stratum radiatum, cells that shift the spiking theta phase are closer to stratum oriens, that is, anatomically closer to where thetaassociated HFOs are best detected in the hippocampus.

Conclusions
In summary, here we used comodulation analyses to unveil different types of fast LFP activity in multisite recordings from the parietal cortex to the hippocampus. The results show that theta phase modulates the amplitude of three separated bands: LG, HG, and HFOs. In addition, we showed that the modulation of multiunit activity by theta phase translates into spurious modulation of a wide range of fast LFP activity, termed SLHFO, which overlaps in frequency with genuine oscillations. The bulk of these results corroborate and extend recent work of ours (Scheffer-Teixeira et al., 2012, 2013 and should help placing theta-associated HFOs as a new type of cortical oscillation . As such, HFOs would benefit from a thorough investigation by independent groups, as it has been the case for gamma oscillations. In this sense, we hope our results serve as a guide for identifying HFOs in other laboratories worldwide.