Distinct directional couplings between slow and fast gamma power to the phase of theta oscillations in the rat hippocampus

It is well‐established that theta (~4–10 Hz) and gamma (~25–100 Hz) oscillations interact in the rat hippocampus. This cross‐frequency coupling might facilitate neuronal coordination both within and between brain areas. However, it remains unclear whether the phase of theta oscillations controls the power of slow and fast gamma activity or vice versa. We here applied spectral Granger causality, phase slope index and a newly developed cross‐frequency directionality (CFD) measure to investigate directional interactions between local field potentials recorded within and across hippocampal subregions of CA1 and CA3 of freely exploring rats. Given the well‐known CA3 to CA1 anatomical connection, we hypothesized that interregional directional interactions were constrained by anatomical connection, and within‐frequency and cross‐frequency directional interactions were always from CA3 to CA1. As expected, we found that CA3 drove CA1 in the theta band, and theta phase‐to‐gamma power coupling was prominent both within and between CA3 and CA1 regions. The CFD measure further demonstrated that distinct directional couplings with respect to theta phase was different between slow and fast gamma activity. Importantly, CA3 slow gamma power phase‐adjusted CA1 theta oscillations, suggesting that slow gamma activity in CA3 entrains theta oscillations in CA1. In contrast, CA3 theta phase controls CA1 fast gamma activity, indicating that communication at CA1 fast gamma is coordinated by CA3 theta phase. Overall, these findings demonstrate dynamic directional interactions between theta and slow/fast gamma oscillations in the hippocampal network, suggesting that anatomical connections constrain the directional interactions.

In the entorhinal-hippocampal network, the gamma-band activity is often divided into slow gamma (~25-55 Hz) and fast gamma (~60-100 Hz) oscillations (Colgin et al., 2009). Slow and fast gamma oscillations emerge in different phases of the theta rhythms and thus provide a mechanism to temporally segregate potentially interfering information (Hasselmo, Bodelon, & Wyble, 2002a). It has demonstrated that slow gamma activity was associated with prospective coding while fast gamma was related to retrospective coding (Zheng, Bieri, Hsiao, & Colgin, 2016). Moreover, fast and slow gamma oscillations have recently been shown to be associated with different running speeds in rats. Fast gamma activity driven by medial entorhinal cortex is predominant during high running speeds, whereas slow gamma arriving from CA3 is prevalent at low running speeds (Ahmed & Mehta, 2012;Zheng, Bieri, Trettel, & Colgin, 2015). Taken together, slow and fast gamma oscillations are suggested to correspond to distinct functional states in the entorhinal-hippocampal network (Colgin, 2015a).
While there is strong evidence that gamma power is coupled to the theta phase, it remains unclear if the theta phase causally modulates the gamma power or whether burst of gamma power phase adjusts the theta oscillations. Local field potentials (LFPs) recordings from hippocampal structures in behaving rats provide an excellent opportunity to investigate the directional interactions between theta oscillations in relation to slow and fast gamma activity. To this end, we analysed LFPs recorded from CA1 and CA3 in rats exploring an open field environment. We used spectral Granger causality (GC) and phase slope index (PSI) to investigate the dominant within-frequency information flow between CA3 and CA1 (Dhamala, Rangarajan, & Ding, 2008;Nolte et al., 2008). While PSI is mainly a measure estimating the dominant direction of interaction, GC allows for estimating bidirectional interaction (Nolte, Ziehe, Krämer, Popescu, & Müller, 2010). To reconcile PSI and GC, we therefore also estimated interregional GC differences and compared it to PSI. More critically, a novel measure of cross-frequency directionality (CFD), which can evaluate the directional coupling between the phase of slow oscillations and the power of fast oscillations in a robust manner, was applied to these data (Jiang, Bahramisharif, Gerven, & Jensen, 2015). The CFD measure was validated by means of simulation studies in Jiang et al., 2015, which could be associated with delays. We speculate that these delays are a consequence of neural transmission, albeit this still needs to be demonstrated by further electrophysiological recordings.
Along the tri-synaptic pathway within the hippocampus, it is well-established that CA3 projects to CA1 but not the reverse, because pyramidal cells of CA3 provide a major input to CA1 through Schaffer collaterals (Amaral & Witter, 1989;Andersen, Bland, Myhrer, & Schwartzkroin, 1979). As hypothesized based on anatomical connection constrains, we predicted that theta and gamma within-frequency and cross-frequency directional interactions were always from CA3 to CA1.

| Animals
Six male Long Evans rats weighing between 350 and 500 g were used in this study. Recording methods are similar to those described previously (Bieri, Bobbitt, & Colgin, 2014;Zheng et al., 2015). In brief, electrode drives containing independently moveable tetrodes ("hyperdrives") (Gothard, Skaggs, Moore, & McNaughton, 1996) were implanted in CA1 and CA3. The rats' light-dark cycle was inverted (lights off from 8:00 to 20:00 and lights on from 20:00 to 8:00) to maintain a behavioural test facility (Beeler, Prendergast, & Zhuang, 2006), and the behavioural sessions were executed during the dark cycle. Behavioural training and data collection started at least one week after recovery from the surgery. The rats were food-deprived to the level of about 90% of freefeeding weight during data collection. All experiments were approved by the University of Texas at Austin and conducted according to the protocol of the United States National Institutes of Health Guide for the Care and Use of Laboratory Animals, in accordance with the Society for Neuroscience's Policies on the Use of Animals in Neuroscience Research.

| Tetrode placement
Over the course of a few weeks after surgery, tetrodes were slowly lowered towards either the CA1 or CA3 stratum pyramidale. For each hyperdrive, one tetrode was placed at the corpus callosum or higher and used as a reference, which typically is used for such recordings as it is considered relatively silent (Bieri et al., 2014;Zheng et al., 2016Zheng et al., , 2015. To make sure that reference tetrodes were placed in a silent location, they were recorded continuously against the ground. After the experiment finished, all recording locations were histologically verified. The final recording sites were located in or close to the CA1 and CA3 stratum pyramidale ( Figure 1). CA1 tetrodes were selected as the tetrodes that were located approximately in the middle of the proximodistal axis of CA1. CA3 tetrodes were selected as those tetrodes located as close as possible to the middle of the proximodistal axis of CA3 (i.e. as close as possible to the middle of CA3b).

| Data collection
Data collection began when cells were recorded approximately at the proper depth for the region of interest with amplitudes exceeding ~4-5 times the noise levels. EEG characteristics (e.g. polarity and amplitude of sharp waves in hippocampus, theta modulation) additionally helped establish recording locations. We used a Neuralynx data acquisition system (Neuralynx) to record the data. A unity gain, multichannel headstage (HS-54; Neuralynx) was connected to the recording drive. Continuous LFP recordings were sampled at 2000 Hz and digitally 0.1-500 Hz bandpass filtered. By using a breakout board (MDR-50 breakout board; Neuralynx), we duplicated the reference signal and the reference signal was recorded against the ground continuously. LFPs were then obtained by differentiating against the referenced signal.

| Behaviour
Rats restarted behavioural training at least one week after recovering from surgery. Six rats were trained to run in a 60 cm × 60 cm open field enclosure in three 10-min sessions each day. Small pieces of cookies were randomly scattered throughout the enclosure to motivate rats to run. To make sure rats were familiar with the environment, data acquisition was conducted after two days of familiarization training. Following each recording session, rats had about 10 min of rest in a towel-lined, elevated flowerpot. In each recording session, the LFP recording from the whole session was divided into 1s epochs and averaged running speed was calculated for each epoch. To make sure the rats were active exploring, only the epochs with highest 33.3% (speed rank > 67.7%) running speed were used for later data analysis.

| Spectral analysis
The analyses were done using the FieldTrip toolbox (Oostenveld, Fries, Maris, & Schoffelen, 2011) and in-house MATLAB scripts (MATLAB and Statistics Toolbox Release 2014b, The MathWorks, Inc., Natick). Spectral coherence, GC and PSI estimations were computed by a fast Fourier transform (FFT) using a multitaper approach (7 Slepian tapers) (Mitra & Pesaran, 1999). The 1 s epoch lengths resulted in 1 Hz spectral resolution and a ± 2 Hz spectral smoothing by multitapering. For GC analysis, non-parametric spectral matrix factorization was applied to the cross-spectral density. Non-parametric GC analysis is superior to the parametric approach since it does not require the autoregressive model order to be estimated (Dhamala et al., 2008). Additionally, PSI was applied to assess the within-frequency directionality. PSI is a robust method to estimate the direction of information influx by computing the slope of phase differences in a pre-specified frequency range (Nolte et al., 2008). We used 2 Hz bandwidth to calculate the phase slope.

| Cross-frequency coupling (CFC) and directionality (CFD) analysis
Let x s denote the raw signal at segment s. Let y v,s denote the power envelope of segment s where v is the frequency of the fast oscillation. We define X s and Y v,s as the FFTs of x s and y v,s , respectively. Let v,s = X s (Y v,s ) * be the cross-spectrum between X s and Y v,s where "*" denotes the complex conjugate. Each of these complex-valued vectors is centred at frequencies f ∈ 0,Δf ,2Δf , … , n FFT 2 Δf where Δf = Fs∕n FFT is the frequency resolution with sampling frequency Fs and length of the Fourier transform n FFT . We use notation X s (f ) ,Y v,s (f ) and v,s (f ) to denote elements of these vectors centred at f. The measure of cross-frequency coherence is based on the coherence between the power envelope of a high-frequency signal and low-frequency of the raw signal centred at f: where S is the number of data segments.
(1) CFD is computed to evaluate the directionality of interactions between neuronal oscillations, which is based on the PSI between the phase of slower oscillations and the power envelope of faster oscillations (Jiang et al., 2015). PSI is a robust method to quantify directionality because it allows one to infer whether one signal is leading or lagging a second signal by considering the slope of phase differences in a pre-specified frequency range (Nolte et al., 2008). The assumption is that a constant lag in the time domain can be translated into phase differences, which will change linearly with frequency in the considered range. Let the complex coherency be defined as.
It should be noted that Equation (2) and Equation (1) are related because Equation (1) is the modulus squared of Equation (2). The CFD between signal x and the power envelope of the signal y v at frequency tile (v,f ) is defined as: where β is the bandwidth used to calculate the phase slope and Im denotes the imaginary part.
For the CFC and CFD calculations, the high-frequency power envelope was extracted using a sliding time window approach. This was implemented by applying a discrete Fourier transform to successive segments of the data after multiplying with a Hanning taper (5 cycles long with respect to the frequency of interest). High-frequency power envelope extraction was done from 20 to 100 Hz in steps of 2 Hz. To calculate CFD, the bandwidth β for estimating the phase slope index was set to 2 Hz at central phase frequency from 2 to 20 Hz in 1 Hz steps. It was worth to mention that we used zero padding to increase the frequency resolution to 0.5 Hz, so the number of frequency bins to compute the phase slope index was 4. CFC and CFD values are normalized by dividing the absolute maximum value of all intraregional and interregional interactions per animal, thus all CFC and CFD values are in the range (0, 1] and [−1, 1], respectively.

| Statistic testing
To assess the significance of within-frequency GC differences, we applied a non-parametric cluster permutation approach to account for multiple comparison correction (Maris & Oostenveld, 2007). First, for every frequency bin, we computed a two-sided two-sample t test between CA3 to CA1 and CA1 to CA3 Granger influences, resulting in t statistic map. Cluster candidates were determined by t values that exceeded the 95% percentile of the t statistic (p < .05). To form clusters, at least two neighbouring t map candidates were required and cluster scores were computed as the summed t values within the cluster. Next, we circularly shifted a random number of time points in CA1 while keeping the original CA3, and recomputed the GC, permuted t values and cluster scores. This procedure was done 1,000 times, resulting in 1,000 maximum cluster scores in the cluster level reference distribution. By comparing to the cluster permutation distribution, observed cluster scores higher than the 97.5th percentile or lower than the 2.5th percentile were considered statistically significant at p < .05. For PSI statistic assessment, the procedure was similar, but we used one-sample t test as the test statistic instead because PSI is symmetric. To evaluate the CFC and CFD statistics, we used R Package "ARTool" for nonparametric two-way repeated measure ANOVAs (https :// cran.r-proje ct.org/web/packa ges/ARToo l/index.html) due to the small number of animals (n = 6). The ARTool relies on a preprocessing step that "aligns" data before applying averaged ranks, after which point common ANOVA procedures can be used (Wobbrock, Findlater, Gergle, & Higgins, 2011).

| Within-frequency interaction in the hippocampal network
First, we quantified the power spectra in the CA1 and CA3. During the active exploration, there was a prominent peak in the 6-10 Hz theta band in the power spectra in both CA1 and CA3 (Figure 2a). The interregional CA1 and CA3 synchronization was then computed by the coherence metric. The coherence spectrum generally revealed two distinct regimes: one in the theta band (~4-10 Hz) and a less prominent broad peak in the slow gamma band (~30-50 Hz) (Figure 2b). Next, we investigated frequency-specific directional influences by computing the spectral GC and PSI, respectively (Dhamala et al., 2008;Nolte et al., 2008). The GC measure from CA3 to CA1 revealed a distinct peak in the theta (~8 Hz) (Figure 2c). A cluster permutation approach was applied to statically access the difference in GC values while controlling for multiple comparisons over frequency bins. This was done by keeping the original CA3 while circularly shifting a random number of time points in CA1 1,000 times and recalculating the GC difference t-maps to obtain a reference distribution. This analysis confirmed that the GC influence from CA3 to CA1 was significantly stronger than the influence from CA1 to CA3 in the theta band ( Figure  2c). Likewise, PSI showed that theta mediated information flow was from CA3 to CA1 (Figure 2d), similar to our GC findings. In short, this within-frequency directionality analysis demonstrates that CA3 activity drives CA1 activity in the theta band more than the reversed direction.

| General cross-frequency interactions in the hippocampal network
Next, we quantified the cross-frequency couplings in the hippocampal network. This was first done by calculating the grand average of the interactions both within and between CA1 and CA3 regions (CA1 phase to CA1 power; CA3 phase to CA3 power; CA1 phase to CA3 power; CA3 phase to CA1 power). CFC provides a measure for determining how much the power envelope of faster oscillations correlates with the phase of slower oscillations. The grand average CFC map showed that the phase of theta oscillations (4-10 Hz) was strongly coupled to gamma power (~30-100 Hz) (Figure 3a). We next estimated the CFD. The CFD quantifies whether the phase of slower oscillations drives the power of faster oscillations (positive CFD) or conversely whether the power of faster oscillations drives the phase of slower oscillations (negative CFD). Visual inspection of the grand averaged CFD map revealed two distinct gamma patterns: negative CFD in the range of theta (4-10 Hz) phase to 30-50 Hz gamma power and positive CFD in the range of theta (4-10 Hz) phase to 60-90 Hz gamma power ( Figure  3b). Thus, we defined these distinct two gamma frequency bands as slow gamma (30-50 Hz) and fast gamma (60-90 Hz), respectively.

| Hippocampal subregional crossfrequency interactions
We then investigated hippocampal subregional cross-frequency interactions by quantifying CA1 and CA3 intraregional and interregional cross-frequency interactions. Since slow and fast gamma power was directionally coupled to theta phase differentially in the hippocampal network, we evaluated theta phase (4-10 Hz) to respectively slow gamma (30-50 Hz) and fast gamma power (60-90 Hz) for CFC and CFD effects separately. The measures we calculated by averaging the corresponding values in the defined phase and power frequency ranges identified in Figure 3. To assess the statistical significance, two-way repeated measure ANOVA analyses were performed for intraregional and interregional interactions. We first quantified the CA1 and CA3 intraregional cross-frequency interactions, using a two-way repeated F I G U R E 2 Power spectrum, coherence, GC and PSI in CA1 and CA3 regions. (a) CA1 (red line) and CA3 (blue line) power spectra. (b) Coherence between CA1 and CA3 (blue line). (c) GC spectra between CA3 and CA1 in both directions (CA1 to CA3, blue line; CA3 to CA1, red line). Frequency ranges with significant differences marked by black dashed line (p < .05). (d) PSI between CA3 and CA1. Frequency ranges with significant differences marked by black dashed line (p < .05). Note that PSI is symmetric since PSI from CA3 to CA1 is the negative number of PSI from CA1to CA3. The positive PSI suggests that information flow is from CA3 to CA1 and vice versa. The shaded area represents standard deviations [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 3 Grand averaged CFC and CFD in the hippocampal
network obtained by averaging CA1 phase to CA1 power, CA3 phase to CA3 power, CA1 phase to CA3 power, and CA3 phase to CA1 power interactions. measure ANOVA analyses with factors region (CA1 vs. CA3) and gamma range (slow gamma vs. fast gamma power). We only found significant main effect of region (F 1,20 = 16.08, p < .005) (Figure 4a), indicating CA3 intraregional CFC between theta phase and gamma power (mean = 0.24, SD = 0.13) was significantly higher than CA1 intraregional CFC (mean = 0.09, SD = 0.06). With respect to the CA1 and CA3 intraregional CFD, there were no significant main or interaction effects (Figure 4b).
With regard to interregional directional phase-topower coupling, the CA1 and CA3 interregional ANOVA in CFD revealed the significant main effect of direction (F 1,20 = 8.71, p = .008) (Figure 5b). Specifically, the negative CFD of CA1 phase to CA3 slow gamma power coupling (mean = −0.071, SD = 0.070), implicating CA3 slow gamma power was driving CA1 theta phase. Conversely, the F I G U R E 4 CA1 and CA3 intraregional CFC and CFD. (a) CA1 and CA3 intraregional CFC. Left panel: CA1 and CA3 intraregional CFC phase-power comodugrams; Right panel: Mean theta phase (4-10 Hz) to slow gamma power (30-50 Hz) or fast gamma power (60-90 Hz) CA1 and CA3 intraregional CFC, which were obtained by averaging the CFC values in the defined phase and power frequency ranges. The error bars represent standard deviations across six rats. (b) Similar to A but for CA1 and CA3 intraregional CFD. *** p < .005 [Colour figure can be viewed at wileyonlinelibrary.com] positive CA3 phase to CA1 fast gamma power (mean = 0.26, SD = 0.13), implicating CA3 theta phase controlling CA1 fast gamma power. Additionally, there appeared to be two separate CFD patterns in the 60-90 Hz fast gamma range: positive 4-8 theta to fast gamma CFD and negative 8-10 Hz theta to fast gamma CFD. However, after closer inspection, this was driven by one rat and thus not general. The statistics revealed that the main effect of gamma range was also significant (F 1,20 = 5.35, p = .03) with positive CFD for fast gamma (mean = 0.24, SD = 0.14) and negative CFD of slow gamma (mean = −0.07, SD = 0.05). There was no significant interaction between direction and gamma range factors in the CA1 and CA3 interregional CFD. Of note, when further examining each CA1 and CA3 interregional CFD individually, we found significant negative CA1 theta phase to CA3 slow gamma power CFD (Wilcoxon signed-rank test, p = .03) and positive CA3 phase to CA1 fast gamma power CFD (Wilcoxon signed-rank test, p = .03), suggesting they were the main driving factors underlying the significant main effects of direction and gamma. Overall, these indicate distinct directional couplings between theta phase and slow/fast gamma power in the rat hippocampus CA1-CA3 circuit. In particular, the CA1-CA3 circuit cross-frequency dynamics are mainly reflected by CA3 slow gamma driving CA1 theta phase and CA1 theta phase controlling CA3 fast gamma.

| DISCUSSION
We here have provided novel insights into the neuronal dynamics supporting both within-frequency and cross-frequency directed information flow between hippocampal subregions. This was done by analysing LFP recordings from hippocampal subregions CA1 and CA3 in exploring rats. We found that CA3 drove CA1 in the theta band. We then confirmed the presence of prominent coupling between theta phase and gamma power within and between CA3 and CA1 regions (Belluscio et al., 2012;Bragin et al., 1995;Colgin et al., 2009). Importantly, we demonstrated distinct directional functional couplings between slow and fast gamma power to the phase of theta oscillations in the rat hippocampus CA1-CA3 network: CA3 slow gamma activity is in control of CA1 theta oscillation while CA3 theta phase controls CA1 fast gamma activity.

| Distinct directional couplings in slow and fast gamma power to theta phase
What might be the purpose of the coupling between CA3 slow gamma power to the CA1 theta phase from a functional point of view? Theta-modulated slow gamma has been proposed to facilitate memory retrieval (Colgin, 2015b). Memory retrieval is thought to be supported by CA3 due to its extensive recurrent collaterals (Brun et al., 2002;Steffenach, Sloviter, Moser, & Moser, 2002;Treves & Rolls, 1991). Phase synchronization between CA3 and CA1 in the slow gamma band might facilitate the transfer of retrieved memory representations from CA3 to CA1 (Carr, Karlsson, & Frank, 2012;Colgin et al., 2009). One possibility is that bursts of slow gamma activity in CA3 may phase-reset theta activity in CA1 to ensure that memory representations reflected by gammaband synchronization are effectively transmitted from CA3 to CA1 within discrete theta cycles. This is consistent with the notion that related information is packaged together within individual theta cycles (Colgin, 2013).
What is the functional role of CA3 theta oscillation entraining CA1 fast gamma activity? Theta-modulated fast gamma has been proposed to facilitate memory encoding (Colgin, 2015b). One explanation is that CA3 theta phase coordinates multi-item memory information represented by CA1 fast gamma activity. This theta-gamma code is used to format memory encoding, whereby different information is represented in different fast gamma subcycles of a theta cycle.

| What makes an oscillation distinct and how should we define their ranges?
While gamma oscillations in the hippocampus are often reported and discussed as if they are well defined, the exact gamma sub-band ranges are often opaque and inconsistent across studies (Belluscio et al., 2012;Bieri et al., 2014;Fernandez-Ruiz et al., 2017;Schomburg et al., 2014;Zheng et al., 2015). This is mainly due to methodological differences, agreed-upon phenomenon. For example, Belluscio et al. defined the ranges based on distinct gamma bands associated with different phases of theta waves (Belluscio et al., 2012). Schomburg et al. (2014) identified distinct gamma sub-bands with phase-amplitude coupling comodugrams. Fernandez-Ruiz et al. (2017) determined different gamma via current source density analysis as well as independent component analysis decomposition of the multi-electrode LFP. Here, we defined gamma sub-bands with phase-amplitude directional coupling comodugrams. Although it is difficult to determine the best practice of defining gamma sub-band ranges, future studies should try to reconcile as much as possible.

| Concerns on artifactual coupling
We characterized the relationship between theta oscillations (4-10 Hz) and gamma oscillations . In particular, we demonstrated CFC between theta and gamma oscillations by quantifying theta phase-to-gamma power coupling. A concern when interpreting CFC results pertains to whether fast oscillations are associated with distinct neuronal activity in the gamma band or if the coupling is explained by the non-sinusoidal shape of theta oscillations (Aru et al., 2015). In the latter case, coupling would be artifactual. The point has been made that non-sinusoidal wave shapes in the theta band can create spurious phase-amplitude coupling (Kramer, Tort, & Kopell, 2008;Lozano-Soldevilla, Huurne, & Oostenveld, 2016). To check this potential confound, we examine the sharpness of theta activity in relation to CFC (Cole et al., 2017). Sharpness defines the asymmetry ratio between the sharpness of the oscillatory signal peaks compared with the troughs as follows: If the signal is perfectly symmetric and sinusoidal, the sharpness ratio is equal to 1. The sharpness ratios of CA1and CA3 are bigger than 1 (CA1: 1.09 ± 0.08; CA3: 1.47 ± 0.19;), indicating they are non-sinusoidal and raising the question whether the slow/fast gamma power were independent oscillators or the by-product of the non-sinusoidal properties. To test between these possibilities, we conducted complementary time-frequency representations (TFRs) of induced power locked to theta oscillation peaks analysis. If slow/fast gamma power is phase locked to the theta rhythm as the CFC suggested, specific theta phase segments (i.e. peaks or thoughts) should be associated with power modulations. As illustrated in Figure 6, we observed (1) Slow gamma appeared around the theta peak while fast gamma exhibited around the theta trough In CA1; (2) Slow gamma was phase locked to the theta peak in CA3. Therefore, slow and fast gamma power indeed modulate within the theta cycle. Moreover, we checked the bicoherence in CA1 and CA3 (Figure 7). Bicoherence has been suggested to assess the level of phase-power and not of phase-phase coupling as commonly accepted (Hyafil, Giraud, Fontolan, & Gutkin, 2015;Shahbazi Avarvand et al., 2018). In the bicoherence comodugrams, true phase-power coupling is characterized by strong bicoherence outside the diagonal regions, while spurious phase-power coupling due to sharp peaks and harmonics is reflected by bicoherence in the diagonal regions (Kovach, Oya, & Kawasaki, 2018). Therefore, the theta-theta bicoherence in both CA3 and CA1 might be related to the harmonic of theta, while the theta-slow/fast gamma bicoherence might suggest the true phase-power coupling. Lastly, both the slow and fast gamma oscillations were associated with fast neuronal oscillations as seen in the spike-field recordings (Colgin et al., 2009). Taken together, these reduce concerns of the CFC in the rat hippocampus being created by non-sinusoidal theta oscillations (see (Jensen, Spaak, & Park, 2016) for discussion).

| Relation to other studies
At the single neuron level, Csicsvari et al. demonstrated CA3 pyramidal neurons discharging CA3 and CA1 interneurons at latencies, in which CA3 pyramidal neurons fired significantly earlier than CA1 interneurons (Csicsvari et al., 2003). Here, we found CA3 slow gamma led CA1 theta activity and CA3 theta activity led CA3 fast gamma. Taken F I G U R E 6 Grand average of time-frequency representations (TFRs) of induced power locked to theta oscillation peaks, which were identified after applying a 4-10 Hz bandpass filter to the data. (a) Top panel: Grand averaged TFRs time-locked to theta peaks (t = 0 s) in CA1. TFRs were calculated for each 0.2 s time window (−0.1 to 0.1s) around the theta peak and then averaged. The colour bar represents the relative power change normalized to the whole 0.2s time window. Bottom panel: Grand mean representation of theta peak-triggered trace over the 0.2 s time window around the theta peak in CA1. (b) Similar to A but for CA3 [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 7 Grand average of CA1 and CA3 bicoherence.
Prominent theta-theta bicoherence in both CA3 and CA1 indicates the second harmonic of theta. Additionally, there are relatively strong theta-slow/fast gamma bicoherence in CA3 and CA1 [Colour figure can be viewed at wileyonlinelibrary.com] together, these might suggest the firing order of CA1 and CA3 neurons might be reflected at the general neuron population level in terms of theta activity, slow gamma activity and fast gamma activity between CA3 and CA1 region. Noting that, slow gamma driving theta activity was also found in the hippocampal circuit previously. For example, both vitro and vivo LFP recordings in rats showed CFC between CA3 gamma (25-40 Hz) and subicular theta (Jackson et al., 2014). To assess if the theta-gamma coupling was associated with a particular delay, they shifted the theta phase across different lags while keeping the gamma amplitude timing constant. Directed CA3 slow gamma to subicular theta interaction was found to be the dominant directional interaction when rats were exposed to a novel open field environment. Moreover, Vaidya and Johnston (2013) demonstrated gamma-to-theta power conversion in the dendrites of CA1 pyramidal neurons. Gamma frequency synaptic bursts could generate theta-frequency components important for oscillatory synchrony. In our study, this gamma-totheta power conversion might occur at slow gamma band between CA3 and CA1 since we found CA3 slow gamma power driving CA1 theta activity. Lastly, Nandi, Swiatek, Kocsis, & Ding (2019) investigated the hippocampus -prefrontal cortex and dentate gyrus (DG) -the Ammon's horn (CA1) interregional directional interactions. The ground truth was provided by the known anatomical connections predicting hippocampus → PFC and DG → CA1. They found that (1) hippocampal high-gamma amplitude was significantly coupled to PFC theta phase, but not vice versa; (2) DG high-gamma amplitude was significantly coupled to CA1 theta phase, but not vice versa. Similarly, we found that the theta and gamma within-frequency and cross-frequency directional interactions were always from CA3 to CA1, suggesting that anatomical connections were constraining the directional connectivity in the hippocampal CA3-CA1 network.
In conclusion, our analysis reveals complex directional interactions between theta and slow/fast gamma oscillations in the hippocampal network. In particular, CA3 slow gamma activity entrains the onset of CA1 theta cycles while CA3 theta oscillation controls CA1 fast gamma activity. These findings provide novel insight into how information flow is controlled in the hippocampus. In future studies, it would be of great interest to study these directional interactions in other behavioural states such as during anesthetized, sleep and memory tasks.