Gamma oscillations coordinate different theta rhythms in the hippocampus

Theta-gamma cross-frequency coupling (CFC) is thought to route information flow in the brain. How this idea copes with the co-existence of multiple theta rhythm generators is not well understood. We have analysed multiple theta and gamma activities in the hippocampus to unveil the dynamic synchronization of theta oscillations across hippocampal layers, and its differential coupling to layer-specific gamma frequency bands. We found that theta-gamma CFC is stronger between oscillations originated in the same hippocampal layer. Interestingly, strong CFC was linked to theta phase locking across layers in a behaviourally related manner, being higher during memory retrieval and encoding. Systematic analysis of cross-frequency directionality indicated that the amplitude of gamma oscillations sets the phase of theta in all layer-specific theta-gamma pairs. These results suggest, contrary to an extended assumption, that layer- and band-specific gamma-oscillations coordinate theta rhythms. This mechanism may explain how anatomically distributed computations, organized in theta waves, can be bound together.


Introduction
Brain oscillations of different frequencies are thought to reflect a multi-scale organization in which information can be bound or segregated in oscillatory cycles [1][2][3] .
Interactions between different oscillations, known as cross-frequency coupling (CFC) 4-6 , have been measured in multiple brain regions during perception [7][8][9] , attention 10,11 and memory formation 7,[12][13][14][15] . It is believed that these interactions play a role in the coordination of local computations and large-scale network communication [1][2][3][16][17][18] . In the hippocampus, theta and gamma oscillations are the most prominent rhythms recorded in freely moving animals 19,20 , and it has been proposed that information transmission between the CA1 region and its afferent regions in CA3 and the entorhinal cortex (EC) is organized in separated gamma frequency channels that are synchronized by the phase of the slower CA1 theta rhythm 18 . The output activity from CA1 has been shown to organize in yet another separated gamma band specifically overlapping the pyramidal cell layer 21,22 and the coupling of the three layerspecific gamma oscillations to the underlying theta rhythm was shown to depend on the behavioural state [21][22][23][24] . However, theta oscillations originating in different anatomical layers are also known to coexist in the hippocampus 19,20,[25][26][27][28][29][30][31] , and therefore theta-gamma interactions need to be understood in the context of multiple rhythm generators.
In addition to the classical medial septum/diagonal band of Broca input imposing a global rhythmicity to the hippocampus and EC, important rhythm generators are located in EC layers II (EC2) and III (EC3), whose activity reach the DG and hippocampus proper though the perforant and temporoammonic pathways respectively, and from CA3 activity reaching CA1 stratum radiatum through the Schafer collaterals 19 . Importantly, although theta oscillation in the hippocampus are most commonly studied as a unique coherent oscillation across hippocampal layers, exhibiting a characteristic amplitude/phase vs. depth variation 19 , the frequency and phase of the CA3 theta rhythm generator was shown to change relatively independently from the EC theta inputs 32 . How these multiple theta rhythm generators and layer-specific gamma oscillations interact in the hippocampus is not well understood. One intriguing possibility is that theta and gamma activities of different laminar origin may represent independent communication channels with possibility to coordinate distributed processes.
Here we investigated the function of layer-specific synchronization of oscillatory activity in the hippocampus of rats freely exploring known and novel environments and resolving a T-maze. Using high density electrophysiological recordings aided by source separation techniques we characterized the dynamic properties of three different theta and three different gamma dipoles in the hippocampus with origins in the CA3 Schaffer collateral layer, the EC3 projection to the stratum laconusum-moleculare and EC2 projection to the mid-molecular layer of the DG, respectively, and found strong support for the existence of relatively independent theta-gamma frameworks. We further characterized theta-gamma interactions between the different layers and established an association with the synchronization state in the hippocampal network. Finally, we investigated the functional role of the characterized theta-gamma coordinated time frames for contextual learning.

Layer-specific theta and gamma oscillations
We performed electrophysiological recordings using linear array electrodes across the dorsal hippocampus in five rats (see Supplementary Methods, Figure 1 The power spectra of these signals exhibited a clear peak at theta frequency (6)(7)(8) and prominent broadband gamma activity (Figure 1e). The dominant theta current sinks and sources calculated from the three layer-specific LFPs (Figure 1f) showed phase differences consistent with the firing properties of principal neurons in their respective upstream afferent layers 38 . Entorhinal principal cells in EC2 and EC3 have been shown to fire in anti-phase, relative to the theta oscillation and, accordingly, theta current sinks ( Figure 1f) and large amplitude gamma oscillations (Figure 1g) in PP-IC and lm-IC were shifted 180º. CA3 and   EC3 neurons have been previously shown to fire phase locked to discrete gamma band   oscillations in the downstream Sch-IC and lm-IC, respectively 21,39 , with gamma oscillations segregated in the phases of the slower theta wave recorded in CA1 18,21,22 . In good agreement, layer-specific gamma oscillations were segregated in the theta cycle, with lm-IC close to the theta peak and followed by Sch-ICin the descending phase to trough of the cycle and PP-IC in the transition from the trough to the ascending phase of the theta cycle ( Figure 1g).
Overall, these results support the use of multichannel recordings and source separation tools to investigate interactions between theta and gamma current generators in multiple layers of the hippocampal formation. (a) Depth profiles of the electrophysiological signals recorded in the dorsal hippocampus (32 recordings, sites spaced every 100 μm) evoked by an electric pulse stimulating the perforant pathway (left panel) or during resting activity (right panel). Black traces represent the LFPs and color maps the corresponding current source density (CSD). Evoked activity was used to consistently localize the electrodes during implantation. Or, stratum oriens; pyr, pyramidal layer; rad, stratum radiatum; lm, stratum lacunosum-moleculare; gc, granule cell layer; hil, hilus. (b and c) Voltage-and CSD-loadings of the three layer-specific LFPs extracted with the ICA, with maximum loadings overlapping the corresponding afferent layers in the str. radiatum (Sch-IC), lacunosum-moleculare (lm-IC) and the molecular layer of the dentate gyrus (PP-IC). (d) Position of the recording electrode in one representative animal (arrows). The histological section is immunostained with GFAP antibodies. (e) Power spectrums of the three ICs show a clear peak at the theta frequency and a broadband gamma activity. (f) Representative low-pass filtered CSD traces of each IC. (g) Broadband gamma activity of each IC along the theta cycle recorded in the CA1 pyramidal layer.
(h) Histogram of the delay difference between the theta rhythms of the ICs and the LFP from CA1 pyramidal cell layer. The trough and the peak of the cycle have a value of 0 and π radians respectively. Distribution have their maximum at the preferred delay between theta oscillations (marked with dots), and their width reflect the wide variability in delays across cycles.

Different theta frameworks coexist in the dorsal hippocampus
The frequency and phase of the CA3 theta rhythm generator has been shown to change relatively independently from the EC theta inputs 32 . However, theta oscillations in the hippocampus are commonly viewed as a unique coherent oscillation. Taking    (a) Power spectrum of the IC-LFPs during high (blue, SI > 0.8) and low (red, SI < 0.4) theta synchronization epochs. A strong increase and right-shift of the theta peak can be seen during theta synchronization. (b) Theta power correlates with SI in lm-IC and PP-IC (black lines represent statistically significant linear correlations; R=0.42/0.32, p<0.01/0.05, respectively), (c) Gamma power is not correlated with SI. (d) Theta frequency significantly correlates with SI in all IC-LFPs (R=0.39/0.80/0.86, p<0.01/0.0001/ 0.0001, respectively). Horizontal dashed lines mark the same frequency in all panels. Here and in the next figures, box-and-whisker plots, represent the complete data range (whiskers), the first and third quartiles (box) and the median (red line). Points are drawn as outliers (red crosses) if they are larger than q3 + 1.5(q3 -q1) or smaller than q1 -1.5(q3 -q1), where q1 and q3 are the 25th and 75th percentiles, respectively.

Theta-gamma CFC reflects layer-specific interactions
The existence of different theta frameworks opens the possibility to multiple thetagamma interactions (Figure 4a). We therefore computed the theta-gamma CFC within and between hippocampal layers. We first measured the amplitude of gamma oscillations in the three IC-LFPs referenced to the phase of theta recorded in the CA1 pyramidal layer, as is usually done 14,15,18,[21][22][23][24]40 . The results (Figure 4b) confirmed previous findings showing coupling between CA1 theta and a slow gamma band of CA3 origin (Sch-IC; maximal modulation at 37.5 ± 5 Hz, CA1S) 18,21,22 and a medium gamma band of EC3 origin (lm-IC; 82.5 ± 4 Hz medium gamma, CA1M) 18,21,22 . They also revealed an additional theta-nested fast gamma band (130 ± 10 Hz,) in the mid-molecular layer of the DG (DGF) overlapping the terminal field of EC2 inputs, compatible with the previously found theta-gamma CFC in the DG 4 . DGF and CA1M in PP-IC and lm-IC, respectively, were locked close to 180º phase ( Figure 4d) as shown before 38 . CA1M activity was closely followed by CA1S in the theta cycle (Figure 4d) 21 . The key new finding in our analysis was the existence of thetagamma coupling dominant within each layer, this is, the CFC between the theta and gamma oscillations recorded in the same IC-LFP was stronger than any between-layer combination (Figure 4c and 4d and Figure S3). Theta-gamma CFC reflects layer-specific interactions and further supports the existence of independent theta-gamma synchronization frameworks in the hippocampus 32 , rather than a single theta framework that implies a unique carrier theta wave to which the gamma activity is multiplexed in segregated theta-gamma channels. (a) Representative theta-and gamma-filtered traces of lm-IC showing how the gamma envelope is phase-locked at the trough of the theta oscillation recorded in the same IC-LFP, but not with the oscillation recorded in the CA1 pyramidal layer (pyr CA1). (b) CFC between gamma amplitude recorded in the different IC-LFPs and the theta phase of pyr CA1. (c) Same as (b) but using the theta phase recorded in the corresponding IC-LFPs. Same color scales, representing the modulation index (MI), have been used for the same IC-LFPs in panels (b) and (c). (d) Strength of the CFC across all gamma and theta oscillations recorded in the three IC-LFPs and pyr CA1. MI is color-coded. The location and width of the rectangles indicate the theta phase at which gamma amplitude is coupled. Theta waveforms (black traces) are extracted as the average of all theta cycles in the corresponding signals. The highest CFC strength was always found between theta and gamma oscillations of the same anatomical layer (See Figure S3 for statistics).

Synchronization between theta frameworks is associated to local theta-gamma CFC
Knowing that theta-gamma CFC is layer-specific ( Figure 4) and that different layers can oscillate relatively independently ( Figure 2 and 3) 32 , we explored theta-gamma CFC accounting for the different synchronization states. This analysis unveiled a striking correlation between the CFC and theta synchronization ( Figure 5a). Strong theta-gamma modulation was associated to high SI values, while weak or nearly absent CFC was found in periods of low SI (Figure 5b). This result indicated that within-layer CFC was associated to the synchronization between layers and allowed us to hypothesize that CFC could represent a mechanism to synchronize theta frameworks.
Theta and gamma oscillations reflect the extracellularly added synaptic and active dendritic currents of two processes occurring at different timescales, the second being tightly paced by inhibitory neurons 37 . We asked which of these processes was driving the CFC between the two frequencies. We computed the cross-frequency directionality index 41 (CFD) (Supplementary Methods), based on the phase-slope index, which computes the phase difference between two signals as an indication of directed interactions between both frequencies. An increase of the phase difference between the theta phase and the gamma amplitude with frequency gives rise to a positive slope of the phase spectrum (i.e. a positive CFD value) when the phase of the slow oscillation sets the amplitude of the fast, and negative otherwise. As shown in Figure 5c for the averaged data, and Figure S4   (c) Within-layer CFD analysis reveals maximum negative values (amplitude-phase coupling) for those pairs of theta-gamma oscillations with the highest CFC (encircled area). These results suggested a driving role of gamma activity over the theta phase, which is maintained during epochs of high and low SI.

Behavioural modulation of theta synchronization and CFC
The strength of CFC in the hippocampus has been shown to correlate with learning 15 and we found here that it correlates with theta coherence across hippocampal layers ( Figure   5a). We hypothesized that the CFC represents a mechanism to synchronize theta oscillations and coordinate information transmission across hippocampal layers. This hypothesis allows us to predict that CFC and theta synchronization would predominate during learning, maximally when the animal updates an existing memory with novel information, a condition likely requiring coordination between incoming sensory inputs, the retrieval of the stored contextual information and updating the memory with new information 42-44 . Therefore, in our final set of experiments we tested this prediction using two behavioural paradigms. In the first one, novel information was added on top of a stored contextual memory. In the second, a hippocampus-dependent delayed spatial alternation task was used in which the animal needed to remember the arm visited in the previous trial and update the memory with the choice made in the current trial 45-47 . In the first task, after habituation to an open field (8 min session 1 per day during 8-10 days), we introduced a novel tactile stimulus in the floor of the otherwise unchanged field (novelty session, see Supplementary Methods). We computed and compared theta synchrony and CFC between the novelty session and the habituation session the day before. When the animal entered the arena, the theta SI was high and comparable in both conditions during the first two minutes of exploration (Figure 6a, t1). As the animal explored the context, synchronization remained high during novelty, but rapidly decayed in the known environment (Figure 6a, t2).
Consistent with the notion of information transmission to update an existing memory, by the end of the exploration time both conditions decreased to the same level of theta synchronization (Figure 6a, t3), when the introduced tactile stimulus lost its novelty. The simultaneously computed CFC MI correlated with theta synchronization during the complete session in both conditions, as shown in Figure 6a and 6b. Importantly, CFC strength was higher during the novelty sessions, when the theta synchronization was also higher, and decreased towards the end of the session in parallel with SI ( Figure 6b).
In the second experiment, rats learned in an 8-shaped T-maze to alternate between the left or right arms on successive trials for water reward (Fig. 6d). In this task, the central arm is associated with memory recall, decision making and encoding of the current decision 14,46-48 . We first computed and compared theta-gamma CFC between the side and central arms and found significantly increased modulation in the central arm for the three IC-LFPs (Fig.   6e). This result is consistent with previously shown increased CFC in the same task in the CA1 radiatum and lacunosum-moleculare IC-LFPs 21 , although in that case the authors used a common theta reference recorded in CA1 pyramidal layer. Theta synchronization across hippocampal layers in the T-maze showed comparable average values as during open field explorations ( Fig. 6f) and, again, differentiated between side and central arms (Fig. 6g).
Importantly, the SI was higher in the central arm (Fig. 6g) consistent with the higher CFC MI found in the same behavioural epochs (Fig. 6e). Overall these results support the above prediction and the idea that CFC synchronizes theta frameworks in the hippocampus facilitating information transmission to update a contextual memory. After the habituation period (control), the animals were exposed to a different floor (sand paper) located inside the familiar open field, providing a new tactile stimulus (novelty). An example of the exploration trajectories is shown (bottom). (b) Time evolution of the SI (mean ± s.e.m. across all subjects) during exploration: before (blue) and after (red) the introduction of the novel tactile stimulus. Both conditions have a maximum SI value at the beginning of the task (t1), corresponding to the initial exploration, followed by a decay in control but not in novelty (t2, inset * p<0.05, t-test). Both conditions decrease to the same SI level by the end of the exploration time (t3). (c) CFC ratio computed as the ratio between the MI in the novelty condition with respect to the control one (* p<0.05, t-test). (d) Example of running trajectories during the T-maze task. (e) Ratio between the MI at the center of the maze and that at the sides (* p<0.05, t-test). (f) Average value of SI in the complete T-maze task (mean ± s.e.m. across all subjects). (g) SI ratio between CENTER and SIDE for each possible combination of IC-LFPs (* p<0.05, t-test).

Discussion
Overall, our results provide functional evidence supporting independent theta oscillations in the hippocampus whose coordination can be seen as a mechanism to channel information between hippocampal layers and binds distributed computations. Less synchronized theta states may secure relatively independent processing in local circuits of the hippocampal formation. Theta coordination (phase locking) is achieved by increasing theta frequency in all layers, most notably in those receiving EC inputs, and strongly correlates with the strength of theta-gamma CFC, so that coherent theta oscillation across hippocampal layers is associated with stronger CFC. Furthermore, directionality analysis suggests that band-and layer-specific gamma activity contributes to the synchronization of theta oscillations across layers. We thus hypothesize that the CFC represents a mechanism operated by gamma activity to coordinate theta oscillations in separated regions. In a network with multiple connected nodes, theta-phase locking between specific nodes will further contribute to the directionality of the information flow, habilitating targets between which communication is permitted in defined time windows. We have provided evidence supporting this hypothesis by showing that CFC and the coordination between the theta generators recorded in the hippocampus increase in the mnemonic process.
Extensive previous research has demonstrated the existence of multiple theta rhythms and current generators in the hippocampus and EC (reviewed in 19 . While septal activity is required for theta rhythmicity, and lesions targeting the medial septum eliminate theta oscillation in both structures, intrinsic hippocampal activity from CA3 and extrinsic EC inputs do also contribute to the recorded theta oscillations 19 . Surgical removal of the EC unveils a theta oscillation that depends on the integrity of CA3 and is highly coherent across hippocampal layers 4 . In the presence of an intact EC, however, the coherence between theta signals in the stratum radiatum and lacunosum-moleculare, receiving the inputs from CA3 and EC3 respectively, is reduced 32 . The relative independency of theta oscillations in these layers is supported by our findings, showing that frequency and phase can be regulated independently in these theta generators (Sch-IC and lm-IC, respectively, Figure 2 and 3).
Furthermore, the key new finding of our experiments is that theta coherence between the generators is not fixed; it rather changes dynamically and is regulated by behaviour.
Accordingly, periods of perfect phase locking (SI = 1) between Sch-IC and lm-IC were also frequent ( Figure 2), and enhanced during particular behaviours ( Figure 6). Brain oscillations can be seen as rhythmic changes in neuronal excitability that can define sequential information packages in the framework of assemblies 16,49 . Therefore, the dynamic variation in theta synchrony between hippocampal layers found in this study, likely reflects the existence of multiple theta-coordinated time frames with phase differences between oscillations having a large impact on the timing of principal cells firing in the respective layers. Synchronization of theta frameworks will, in turn, coordinate, though not necessarily synchronize 50 , firing sequences in consecutive hippocampal stations.
Interactions between the phase of the theta oscillation and the amplitude (power) of the gamma activity have been extensively documented and proposed as an effective mechanism to integrate activity across different spatial and temporal scales [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] . In the hippocampus and despite the evidence supporting the existence of multiple theta generators (see above), theta-gamma CFC has been almost exclusively investigated as the statistical dependence between the gamma activity in different frequency bands and a single theta reference, most commonly recorded in the pyramidal layer of the CA1 region. Our analysis demonstrates that phase-amplitude CFC between theta and gamma oscillations in the hippocampus varies for different theta generators, and is stronger when both oscillations are recorded in the same anatomical layer (Figure 4 and Figure S3). This observation, together with previous and important evidence demonstrating that firing of principal cells in CA3 and EC3 is phase-locked to downstream gamma oscillations recorded in the CA1 stratum radiatum (CA1S) and lacusosum-moleculare (CA1M), respectively 18,21,22 , suggest that CFC is a local phenomenon driven by upstream afferences organizing downstream gamma oscillations. Our CFD analysis further supports this interpretation, since it shows a predominant amplitude-to-phase coupling for the same layer-specific gamma frequency bands. This does not mean that 4-11 Hz oscillations in the hippocampus are generated by gamma activity, on the contrary, our data suggest that gamma activities, reflecting the interplay of inhibitory-excitatory networks 51-56 , impose phase shifts on the on-going theta oscillations in their corresponding layers. Therefore, we hypothesize that local gammagenerating circuits, driven by afferences from their respective upstream layers, might not be activated at a particular theta phase, as commonly interpreted, but rather, might be actually setting the phase of the local theta oscillation. This interpretation would also explain phasephase coupling between CA1 theta and CA1S and CA1M 57 , as the consequence of theta phase driven by layer-specific gamma activity entrained by upstream inputs in CA3 and EC3, respectively.
The proposed new scenario provides a mechanism to coordinate distributed computations organized in theta waves, by synchronizing theta oscillations in connected regions through theta-gamma CFC. The highly significant positive correlation between CFC strength and theta synchronization found in our study ( Figure 5) supports this view. In this way, our results link theta-gamma CFC 6 and coherence-based communication 49 , the first being the mechanism to align higher excitability windows to open communication channels between defined network nodes. A number of different studies investigating learning and memory processes in the entorhinal-hippocampal network have indeed, but separately, reported increases in theta-gamma CFC or synchronization in different frequency bands associated with memory performance [1][2][3]16 . The association of both phenomena in our study has been possible by the separation of pathway-specific LFPs, identifying the CFC as a local phenomenon and separating theta activities, that otherwise would have render a mixed readout of theta-gamma interactions, highly dependent on electrode implantation coordinates and interindividual variability 37 . Here we have further demonstrated that the strong association between CFC and synchronization between the three theta oscillations predominantly occur when the animal updates a memory with novel information (Figure 6).
The behavioural test used were selected to contrast behavioural epochs in which the animal explores a known environment, recalling its representation from memory, and simultaneously needs to encode new information (novel stimuli, last arm visited) to update the memory. In this condition, several processing streams likely coexist and need to coordinate in the entorhinal-hippocampal network 42-44 . It is in these conditions that we found an increased association between CFC strength and theta synchronization across layers, suggesting its role in coordinating distributed computation during memory formation.
We have suggested a causal link between CFC strength and theta synchronization through the gamma-driven adjustment of theta phase and frequency (Figure 3 and 5). While dissecting the precise circuit mechanisms supporting the transfer function from gamma activity to theta phase is out of the scope of the present work, several possibilities exist.
Computational works have demonstrated that theta-gamma CFC emerges from the interactions between functionally distinct interneuron populations, as the basket and orienslacunosum-moleculare (OLM) cells, interconnected in a network of principal cells receiving an external theta rhythm generator, such as the septal input 53-56 . The on-going theta oscillation is thus modulated by the inhibitory network, whose activity is known to be reflected in the gamma activity in in vivo extracellular recordings 40,51,52 , and is associated to layer-specific inputs 18,21,22 . Subsets of interneurons can phase-lock to different hippocampal rhythms 58,59 and, interestingly, recent findings showed in the CA1 region that some interneurons can specifically phase-lock to CA1S and others to CA1M, supporting the idea that different classes of interneurons drive slow and medium gamma oscillations 22,24,60 . Thus, an appealing mechanism for gamma-modulation of theta phase would be the control of different interneuron classes by layer specific inputs, which would entrain specific gamma networks modulating principal cell excitability and firing in response to on-going theta inputs, advancing or delaying theta phase. Spiking resonance in principal cells may contribute to this mechanism too, since optogenetic activation of basket interneurons (parvalbumin expressing cells) in the hippocampus and neocortex has been shown to pace pyramidal cell firing in the theta range, by virtue of postinhibitory rebound of Ih activity 51 . In that experiment, theta-band firing of excitatory neurons required rhythmic activation of basket cells, as white noise activation effectively modulated their activity but did not entrained pyramidal theta-band firing 51 , suggesting that feed-forward activation of interneurons from upstream layers or an external rhythmic input (i.e. cholinergic or GABAergic inputs form the septum), are required for resonance amplification. Thus, intrinsic cellular properties and network mechanism may interact to support gamma-dependent coordination of theta phases across hippocampal layers.
Interactions between slow and fast brain oscillations have been measured in multiple brain regions during perception, attention, learning and memory formation 1-3 . Despite its ubiquitous presence in fundamental cognitive processes, its function is largely unknown. Our results provide a mechanism for binding distributed computations packed on theta waves and routing the information flow based on theta-gamma cross-frequency coupling. Important questions remain to be answered. How theta synchronization in the hippocampus relates to hippocampal-neocortical interactions 61,62 known to be favoured at theta and beta frequencies 63,64 and modulated by synaptic plasticity in the hippocampus 65,66 ? The conditions triggering the coordination between theta-gamma frameworks are not well understood, nor are the precise cognitive processes they subserve, but given that theta-gamma uncoupling seems to represent an early electrophysiological signature of hippocampal network dysfunction in Alzheimer's disease [67][68][69][70]

Animals and Surgery
Five male Long-Evans rats, with a weigh of 250-300 g. were trained in different behavioral tasks, with a multichannel electrode recording the electrophysiological activity in the hippocampus. All of them were implanted with a multisite silicon probe (Neuronexus Technologies, Michigan, USA) connected in turn to a jumper consisting of two corresponding connectors joined by 5 cm of flexible cable. An Ag/AgCl wire (World Precision Instruments, Florida, USA) electrode was placed in contact with the skin on the sides of the surgery area, and used as ground. We adjusted the final position of both electrodes using as a reference the typical evoked potentials at the dentate gyrus 1 , so that a maximal population spike in the dentate gyrus was recorded.
After the surgery, the rats were left for at least 10 days until they recovered completely.
During the first 72 hours, they were injected subcutaneously with analgesic twice per day (Buprenorphine, dose 2-5 μg/kg, RB Pharmaceutical Ltd., Berkshire, UK). During 1 week, they had as well antibiotic dissolved in the water (Enrofloxacin, dose 10 mg/kg, Syva, León, Spain). The behavioral tasks were not started until the animals showed no signs of discomfort with the manipulation of the implants.

Data acquisition and preprocessed
During each behavioral task, local field potentials (LFPs) from all sites were acquired simultaneously for each subject at 20 kHz. For LFP analysis, data were downsampled to 2500Hz and filtered with a high-pass filter (>0.5 Hz), and with a Notch filter at 50 Hz and 100 Hz to eliminate the net noise and its first harmonic. To increase the similarities between data sets of different subjects, they were normalized, imposing zero mean and variance one for each complete time series. Four of the previous subjects carried out a memory task in a modified T-maze as has been described previously 2 . It consisted in several tracks in 8-like shape (132/102/80 cm long/wide/high with track wide 8 cm, Figure 6d). The starting point was located at the beginning of the center rail (Figure 6d, start) and the rat was forced to run across that arm In this work, the electrophysiological recordings acquired in the T-maze were split into two groups of epochs as a function of the location during the task: at the center arm, before the Tjunction, or at one side after the reward (Figure 6d, center and side, respectively).

Current source density analysis of LFPs
The first approach to achieve the information of the sources contributing to the LFPs was the use of current source density 3 (CSD) analysis. This method reduces the volume conductions effects through the recording sites computing the second spatial derivate (Laplacian). In this way, it measures the transmembrane currents, providing a spatiotemporal distribution of the sinks and sources (inward and outward currents, respectively).
For probes with common separation between channels, we used the one-dimensional approach, which calculates the CSD from the voltage distribution of the adjacent sites: where ( ) is the LFP recorded at the m-th site, ℎ is the distance between channels and is the conductivity of the extracellular space.
The CSD does not discriminate contributions from different pathways. To assess that question, an independent component analysis (ICA) should be performed.

Independent Component Analysis of LFPs
To obtain the characteristic activity in each of the regions of interest from the LFP, we performed an ICA approach. This type of methodology aims to solve the problem of separating N statistically independent sources that have been mixed in N output channels.
The approach performs a blind separation of patterns, because the different distributions of the sources are unknown. Moreover, it assumes spatial immobility of the sources or, in other terms, a fixed location of the axon terminals. The contribution of their synaptic currents to the LFP conforms the different ICs to unravel.
Its potential to find sources associated with known hippocampal pathways from LFPs records has been previously demonstrated and employed [4][5][6][7]  As the number of generators with significant variance and spatial distributions associated to different current sources 10 is relatively low (between 4 and 7 using electrodes with 32 channels), the ICA algorithm can be optimized by a dimension reduction using the principal component analysis. This method allows an improvement in the computational cost and reduces the number of noisy generators due to the presence of artifacts in the time-series 6 .
Those components with explained variance <1% (respect to the original LFP variance) were considered as noise and rejected. The mathematical validation and practical limitations of ICA algorithms have been extensively investigated through realistic modelling 11,12 .
Thanks to the laminar structure of the hippocampus, with a separation of the neural terminations into different stratums, the ICA allows the isolation of determined postsynaptic activity from other currents. Thereby, it has been extracted two ICs in CA1, located in str.
radiatum (Sch-IC) and str. lacunosum-moleculare (lm-IC), respectively, and another in the dentate gyrus, in the molecular layer (pp-IC). In the case of Sch-IC, the measured activity corresponds to the Schaffer collateral input from CA3 to the pyramidal cells in CA1, while the information in lm-IC is the input from entorhinal cortex layer III to the previous cells in CA1. Regarding the pp-IC, it contains the activity of the perforant pathway, the main input to the hippocampus from the layer II of the entorhinal cortex.

Synchronization
For each theta cycle, a synchronization index (SI) was computed, which measured the variance of phase delay between IC-LFPs in one cycle with respect to the previous one. To do that and considering our IC-LFPS as ( ) (where k=1,2,3), we have bandpass-filtered the signals at theta frequency (6-10 Hz) and extracted their phase using the Hilbert transform.
We denominated ( ) as the time delay from the c trough in ( ) to the closes trough of ( ). Therefore, the SI from i to j can be computed by follows: The maximum time delay ( ) was settled so that at least 99% of the total cycles in every ∆ ( ) have a lower value ( =21 ms in our case).
To assess the synchronization between several oscillations at the same time, we considered the total delay difference as the summation of the delay differences between pairs: To avoid desynchronized states due to a lack of theta oscillations, only those cycles with theta two times higher than delta power where taken into account.

Cross-Frequency Coupling and directionality
To assess the cross-frequency coupling (CFC) between the phase and the amplitude, we have used the modulation index (MI), proposed by Tort and colleagues 13 , which is based on entropy measurements. Given two signals, ( ) and ( ), filtered at two frequencies, one slow and the other faster, we can extract the phase ( ( )) and the amplitude ( ( )), respectively, by the Hilbert transform. Next, each whole cycle in ( ) is divided in bins of the same size, calling < > ( ) the mean amplitude value at the phase bin . From them, we can calculate the entropy measure , defined by: where = 18, and is given by The value of MI is obtained normalizing by the maximum entropy ( ), given by the uniform distribution = 1⁄ (i.e. = log ): A value of MI near to 0 indicates lack of phase-to-amplitude modulation, while larger MI values show a high phase-to-amplitude modulation. The statistical significance has been assessed following the steps proposed by Canolty and colleagues 14 , by a surrogate analysis (n=100 surrogates) in which each surrogate is built by random shifts between the phase and the amplitude of both signals. The values of all series are approximated to a gaussian distribution, which mean value is considered as a significance threshold.
To cover all values represented in each comodulogram of MIs, we have varied both the values of the slow and fast frequencies. Each point of the x-axis has been filtered with a bandwidth of 2 Hz, and a step of 0.5 Hz, while for the y-axis, the bandwidth was fixed at 18 Hz, at least two times the frequency where the maximum theta value is expected 15 .
Although the MI is a measurement of the degree of interaction between the phase and the amplitude, the directionality of this coupling can be assessed. On the one hand, the theta phase would modulate the amplitude at gamma frequencies while, on the other hand, the gamma activity could be leading the phase. To identify who is the driver, we have used the cross-frequency directionality (CFD) index 16 . It is based on the phase-slope index 17  is the complex coherency, is the theta frequency under study, is the number of segments in which the signal has been divided and β is the bandwidth for which the phase slope is measured, and it has been fixed at 2 Hz, 4 times the resolution (∆ = 0.5 Hz).
To provide statistical significance, a new surrogate test (n=100 surrogates) has been developed following the same steps than in the MI analysis. As the result of the CFD can be positive or negative, two thresholds have been established, as the average of the module of all surrogate results, and its inverse. To emphasize the directionality in the region of higher CFC, the MI comodulogram has been redefined as a mask, with values from 0 to 1 (minimum and maximum value MI). Applying this mask to the CFD comodulogram, areas without phase-amplitude coupling will be attenuated, while the main cluster remains constant. Figure S1: Phase difference distribution, computed as the time delay between IC-LFPs theta troughs in one cycle minus the delay in the previous one. In all cases, it presents a Gaussian distribution centred in zero. Figure S2: High synchronization epochs preferentially last three consecutive theta cycles. Each pair of bars compares the probability (mean ± s.e.m) of, given a theta cycle with SI>0.8, find at least a certain number of consecutive cycles with strong synchronization (SI>0.8), and the probability expected by chance. The minimum number of cycles with a significant difference between measured and expected measurements is three theta cycles (p<0.001, t-test), and goes on until six cycles (p<0.01/0.05/0.05, respectively, t-test). Figure S3: Maximum theta-gamma CFC corresponds to oscillations recorded in the same IC-LFP, higher than any between-layer combination. Left/middle/right panel represents the MI (mean ± s.e.m.) between slow/medium/fast gamma recorded from Sch-IC/lm-IC/PP-IC, and the theta phase of all IC-LFPs and the LFP from pyrCA1. Significant stronger MI values were found in all cases when both theta and gamma have the same origin (*/*** p<0.05/0.001, one-way ANOVA, followed by Bonferroni correction). Moreover, there is a trend when lm-IC and PP-IC gamma modulation, with higher values when the coupling is measured using a theta oscillation from a common origin in entorhinal cortex (PP-IC and lm-IC, respectively), compared to the theta rhythm of Sch-IC, generated in CA3. Figure S4: CFD for individual animals measured along all the recording time. Significant negative values (gamma amplitude drives theta phase) were found in 4 of 5 subjects in Sch-IC and lm-IC for those pair of frequencies with maximum CFC (37.5 ± 5 Hz in Sch-IC and 82.5 ± 4 Hz in lm-IC), while negative CFDs were significant in PP-IC for all subjects.